Feature/classical conditionals (#152)

# Description

Adding support so that `simulate` can: 
- run circuits with classical logic (#151)
- run circuits with a single qubit (#156)

Also included some bugfixes of MPSxMPO (I had missed calling `_flush()`
in the newer functions).

# Related issues

#151 
#156 

# Checklist

- [x] I have run the tests on a device with GPUs.
- [x] I have performed a self-review of my code.
- [x] I have commented hard-to-understand parts of my code.
- [x] I have made corresponding changes to the public API documentation.
- [x] I have added tests that prove my fix is effective or that my
feature works.
- [x] I have updated the changelog with any user-facing changes.

docs/changelog.rst

@@ -1,6 +1,20 @@
 Changelog
 ~~~~~~~~~
 
+Unreleased
+----------
+
+* API breaking changes
+    * Removed ``use_kahypar`` option from ``Config``. It can still be set via the ``simulate`` option ``compilation_params``.
+
+* New feature: ``simulate`` now accepts pytket circuits with ``Measure``, ``Reset``, ``Conditional``, ``ClassicalExpBox`` and more classical operations. You can now retrieve classical bit values using ``get_bits``.
+* When calling ``simulate``, the gates on the circuit are no longer sorted by default. Use ``compilation_params["sort_gates"] = True`` to recover this behaviour, which is now deprecated.
+* ``StructuredState`` now supports simulation of single qubit circuits.
+* Some bugfixes on MPSxMPO relating to measurement and relabelling qubits. The bug was caused due to these functions not guaranteeing the MPO was applied before their action.
+* Documentation fixes:
+    * ``apply_qubit_relabelling`` now appears in the documentation.
+    * ``add_qubit`` removed from documentation of MPSxMPO, since it is not supported.
+
 0.7.1 (July 2024)
 -----------------
 

docs/modules/structured_state.rst

@@ -26,6 +26,7 @@ Classes
     .. automethod:: apply_gate
     .. automethod:: apply_unitary
     .. automethod:: apply_scalar
+    .. automethod:: apply_qubit_relabelling
     .. automethod:: vdot
     .. automethod:: sample
     .. automethod:: measure
@@ -52,7 +53,6 @@ Classes
 .. autoclass:: pytket.extensions.cutensornet.structured_state.MPSxMPO()
 
     .. automethod:: __init__
-    .. automethod:: add_qubit
 
 
 Miscellaneous

examples/python/mps_tutorial.py

@@ -1,7 +1,7 @@
 import numpy as np
 from time import time
 import matplotlib.pyplot as plt
-from pytket import Circuit
+from pytket import Circuit, OpType
 from pytket.circuit.display import render_circuit_jupyter
 
 from pytket.extensions.cutensornet.structured_state import (
@@ -145,6 +145,43 @@
 print("Is the inner product correct?")
 print(np.isclose(np.vdot(my_state, other_state), inner_product))
 
+# ### Mid-circuit measurements and classical control
+# Mid-circuit measurements and classical control is supported (only in `MPSxGate` as of v0.8.0). For instance, we can implement the teleportation protocol on a pytket circuit and simulate it:
+
+circ = Circuit()
+alice = circ.add_q_register("alice", 2)
+alice_bits = circ.add_c_register("alice_bits", 2)
+bob = circ.add_q_register("bob", 1)
+# Initialise Alice's first qubit in some arbitrary state
+circ.Rx(0.42, alice[0])
+orig_state = circ.get_statevector()
+# Create a Bell pair shared between Alice and Bob
+circ.H(alice[1]).CX(alice[1], bob[0])
+# Apply a Bell measurement on Alice's qubits
+circ.CX(alice[0], alice[1]).H(alice[0])
+circ.Measure(alice[0], alice_bits[0])
+circ.Measure(alice[1], alice_bits[1])
+# Apply conditional corrections on Bob's qubits
+circ.add_gate(OpType.X, [bob[0]], condition_bits=[alice_bits[1]], condition_value=1)
+circ.add_gate(OpType.Z, [bob[0]], condition_bits=[alice_bits[0]], condition_value=1)
+# Reset Alice's qubits
+circ.add_gate(OpType.Reset, [alice[0]])
+circ.add_gate(OpType.Reset, [alice[1]])
+# Display the circuit
+render_circuit_jupyter(circ)
+
+# We can now simulate the circuit and check that the qubit has been successfully teleported.
+
+print(
+    f"Initial state:\n {np.round(orig_state[0],2)}|00>|0> + {np.round(orig_state[4],2)}|10>|0>"
+)
+with CuTensorNetHandle() as libhandle:
+    state = simulate(libhandle, circ, SimulationAlgorithm.MPSxGate, Config())
+    print(
+        f"Teleported state:\n {np.round(state.get_amplitude(0),2)}|00>|0> + {np.round(state.get_amplitude(1),2)}|00>|1>"
+    )
+    print(f"Measurement outcomes:\n {state.get_bits()}")
+
 # ### Two-qubit gates acting on non-adjacent qubits
 # Standard MPS algorithms only support simulation of two-qubit gates acting on neighbour qubits. In our implementation, however, two-qubit gates between arbitrary qubits may be applied, as shown below.
 

pytket/extensions/cutensornet/structured_state/classical.py

@@ -0,0 +1,135 @@
+# Copyright 2019-2024 Quantinuum
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+##
+#     http://www.apache.org/licenses/LICENSE-2.0
+##
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from typing import Union, Any
+
+from pytket.circuit import (
+    Op,
+    OpType,
+    Bit,
+    BitRegister,
+    SetBitsOp,
+    CopyBitsOp,
+    RangePredicateOp,
+    ClassicalExpBox,
+    LogicExp,
+    BitWiseOp,
+    RegWiseOp,
+)
+
+
+ExtendedLogicExp = Union[LogicExp, Bit, BitRegister, int]
+
+
+def apply_classical_command(
+    op: Op, bits: list[Bit], args: list[Any], bits_dict: dict[Bit, bool]
+) -> None:
+    """Evaluate classical commands and update the `bits_dict` accordingly."""
+    if isinstance(op, SetBitsOp):
+        for b, v in zip(bits, op.values):
+            bits_dict[b] = v
+
+    elif isinstance(op, CopyBitsOp):
+        output_bits = bits
+        input_bits = args[: len(output_bits)]
+        for i, o in zip(input_bits, output_bits):
+            assert isinstance(i, Bit)
+            bits_dict[o] = bits_dict[i]
+
+    elif isinstance(op, RangePredicateOp):
+        assert len(bits) == 1
+        res_bit = bits[0]
+        input_bits = args[:-1]
+        # The input_bits encode a "value" int in little-endian
+        val = from_little_endian([bits_dict[b] for b in input_bits])  # type: ignore
+        # Check that the value is in the range
+        bits_dict[res_bit] = val >= op.lower and val <= op.upper
+
+    elif isinstance(op, ClassicalExpBox):
+        the_exp = op.get_exp()
+        result = evaluate_logic_exp(the_exp, bits_dict)
+
+        # The result is an int in little-endian encoding. We update the
+        # output register accordingly.
+        for b in bits:
+            bits_dict[b] = (result % 2) == 1
+            result = result >> 1
+        assert result == 0  # All bits consumed
+
+    elif op.type == OpType.Barrier:
+        pass
+
+    else:
+        raise NotImplementedError(f"Commands of type {op.type} are not supported.")
+
+
+def evaluate_logic_exp(exp: ExtendedLogicExp, bits_dict: dict[Bit, bool]) -> int:
+    """Recursive evaluation of a LogicExp."""
+
+    if isinstance(exp, int):
+        return exp
+    elif isinstance(exp, Bit):
+        return 1 if bits_dict[exp] else 0
+    elif isinstance(exp, BitRegister):
+        return from_little_endian([bits_dict[b] for b in exp])
+    else:
+
+        arg_values = [evaluate_logic_exp(arg, bits_dict) for arg in exp.args]
+
+        if exp.op in [BitWiseOp.AND, RegWiseOp.AND]:
+            return arg_values[0] & arg_values[1]
+        elif exp.op in [BitWiseOp.OR, RegWiseOp.OR]:
+            return arg_values[0] | arg_values[1]
+        elif exp.op in [BitWiseOp.XOR, RegWiseOp.XOR]:
+            return arg_values[0] ^ arg_values[1]
+        elif exp.op in [BitWiseOp.EQ, RegWiseOp.EQ]:
+            return int(arg_values[0] == arg_values[1])
+        elif exp.op in [BitWiseOp.NEQ, RegWiseOp.NEQ]:
+            return int(arg_values[0] != arg_values[1])
+        elif exp.op == BitWiseOp.NOT:
+            return 1 - arg_values[0]
+        elif exp.op == BitWiseOp.ZERO:
+            return 0
+        elif exp.op == BitWiseOp.ONE:
+            return 1
+        # elif exp.op == RegWiseOp.ADD:
+        #     return arg_values[0] + arg_values[1]
+        # elif exp.op == RegWiseOp.SUB:
+        #     return arg_values[0] - arg_values[1]
+        # elif exp.op == RegWiseOp.MUL:
+        #     return arg_values[0] * arg_values[1]
+        # elif exp.op == RegWiseOp.POW:
+        #     return int(arg_values[0] ** arg_values[1])
+        # elif exp.op == RegWiseOp.LSH:
+        #     return arg_values[0] << arg_values[1]
+        elif exp.op == RegWiseOp.RSH:
+            return arg_values[0] >> arg_values[1]
+        # elif exp.op == RegWiseOp.NEG:
+        #     return -arg_values[0]
+        else:
+            # TODO: Currently not supporting RegWiseOp's DIV, EQ, NEQ, LT, GT, LEQ,
+            # GEQ and NOT, since these do not return int, so I am unsure what the
+            # semantic is meant to be.
+            # TODO: Similarly, it is not clear what to do with overflow of ADD, etc.
+            # so I have decided to not support them for now.
+            raise NotImplementedError(
+                f"Evaluation of {exp.op} not supported in ClassicalExpBox ",
+                "by pytket-cutensornet.",
+            )
+
+
+def from_little_endian(bitstring: list[bool]) -> int:
+    """Obtain the integer from the little-endian encoded bitstring (i.e. bitstring
+    [False, True] is interpreted as the integer 2)."""
+    return sum(1 << i for i, b in enumerate(bitstring) if b)

pytket/extensions/cutensornet/structured_state/general.py

@@ -19,7 +19,14 @@
 
 import numpy as np  # type: ignore
 
-from pytket.circuit import Command, Qubit
+from pytket.circuit import (
+    Command,
+    Op,
+    OpType,
+    Qubit,
+    Bit,
+    Conditional,
+)
 from pytket.pauli import QubitPauliString
 
 try:
@@ -28,6 +35,7 @@
     warnings.warn("local settings failed to import cupy", ImportWarning)
 
 from pytket.extensions.cutensornet import CuTensorNetHandle
+from .classical import apply_classical_command, from_little_endian
 
 # An alias for the CuPy type used for tensors
 try:
@@ -47,7 +55,6 @@ def __init__(
         float_precision: Type[Any] = np.float64,
         value_of_zero: float = 1e-16,
         leaf_size: int = 8,
-        use_kahypar: bool = False,
         k: int = 4,
         optim_delta: float = 1e-5,
         loglevel: int = logging.WARNING,
@@ -84,9 +91,6 @@ def __init__(
                 ``np.float64`` precision (default) and ``1e-7`` for ``np.float32``.
             leaf_size: For ``TTN`` simulation only. Sets the maximum number of
                 qubits in a leaf node when using ``TTN``. Default is 8.
-            use_kahypar: Use KaHyPar for graph partitioning (used in ``TTN``) if this
-                is True. Otherwise, use NetworkX (worse, but easy to setup). Defaults
-                to False.
             k: For ``MPSxMPO`` simulation only. Sets the maximum number of layers
                 the MPO is allowed to have before being contracted. Increasing this
                 might increase fidelity, but it will also increase resource requirements
@@ -153,7 +157,6 @@ def __init__(
             raise ValueError("Maximum allowed leaf_size is 65.")
 
         self.leaf_size = leaf_size
-        self.use_kahypar = use_kahypar
         self.k = k
         self.optim_delta = 1e-5
         self.loglevel = loglevel
@@ -176,29 +179,85 @@ def copy(self) -> Config:
 class StructuredState(ABC):
     """Class representing a Tensor Network state."""
 
-    @abstractmethod
-    def is_valid(self) -> bool:
-        """Verify that the tensor network state is valid.
-
-        Returns:
-            False if a violation was detected or True otherwise.
-        """
-        raise NotImplementedError(f"Method not implemented in {type(self).__name__}.")
+    _lib: CuTensorNetHandle
+    _cfg: Config
+    _logger: logging.Logger
+    _bits_dict: dict[Bit, bool]  # Tracks the state of the classical variables
 
-    @abstractmethod
     def apply_gate(self, gate: Command) -> StructuredState:
-        """Applies the gate to the StructuredState.
+        """Apply the command to the `StructuredState`.
+
+        Note:
+            Only one-qubit gates and two-qubit gates are supported.
 
         Args:
-            gate: The gate to be applied.
+            gate: The command to be applied.
 
         Returns:
             ``self``, to allow for method chaining.
 
         Raises:
             RuntimeError: If the ``CuTensorNetHandle`` is out of scope.
             ValueError: If the command introduced is not a unitary gate.
-            ValueError: If gate acts on more than 2 qubits.
+            ValueError: If the command acts on more than 2 qubits.
+        """
+        self._logger.debug(f"Applying {gate}.")
+        self._apply_command(gate.op, gate.qubits, gate.bits, gate.args)
+        return self
+
+    def _apply_command(
+        self, op: Op, qubits: list[Qubit], bits: list[Bit], args: list[Any]
+    ) -> None:
+        """The implementation of `apply_gate`, acting on the unwrapped Command info."""
+        if op.type == OpType.Measure:
+            q = qubits[0]
+            b = bits[0]
+            self._bits_dict[b] = self.measure({q}, destructive=False)[q] != 0
+
+        elif op.type == OpType.Reset:
+            assert len(qubits)
+            q = qubits[0]
+            # Measure and correct if outcome is |1>
+            outcome_1 = self.measure({q}, destructive=False)[q] != 0
+            if outcome_1:
+                self._apply_command(Op.create(OpType.X), [q], [], [q])
+
+        elif op.is_gate():  # Either a unitary gate or a not supported "gate"
+            try:
+                unitary = op.get_unitary()
+            except:
+                raise ValueError(f"The command {op.type} introduced is not supported.")
+
+            # Load the gate's unitary to the GPU memory
+            unitary = unitary.astype(dtype=self._cfg._complex_t, copy=False)
+            unitary = cp.asarray(unitary, dtype=self._cfg._complex_t)
+
+            if len(qubits) not in [1, 2]:
+                raise ValueError(
+                    "Gates must act on only 1 or 2 qubits! "
+                    + f"This is not satisfied by {op.type}."
+                )
+
+            self.apply_unitary(unitary, qubits)
+
+        elif isinstance(op, Conditional):
+            input_bits = args[: op.width]
+            tgt_value = op.value
+            # The input_bits encode a "value" int in little-endian
+            var_value = from_little_endian([self._bits_dict[b] for b in input_bits])  # type: ignore
+            # If the condition is apply the command in the body
+            if var_value == tgt_value:
+                self._apply_command(op.op, qubits, bits, args)
+
+        else:  # A purely classical operation
+            apply_classical_command(op, bits, args, self._bits_dict)
+
+    @abstractmethod
+    def is_valid(self) -> bool:
+        """Verify that the tensor network state is valid.
+
+        Returns:
+            False if a violation was detected or True otherwise.
         """
         raise NotImplementedError(f"Method not implemented in {type(self).__name__}.")
 
@@ -388,6 +447,10 @@ def get_amplitude(self, state: int) -> complex:
         """
         raise NotImplementedError(f"Method not implemented in {type(self).__name__}.")
 
+    def get_bits(self) -> dict[Bit, bool]:
+        """Returns the dictionary of bits and their values."""
+        return self._bits_dict.copy()
+
     @abstractmethod
     def get_qubits(self) -> set[Qubit]:
         """Returns the set of qubits that ``self`` is defined on."""