refactored execute. The observable is now passed as lists of operators, coefficients and qubits. Implemented functions for handling conversion from Pauli string to quimb operator and coefficients. Fixed: error in parameters handing in expectation. Fixed GATE_MAP to match the common observables

src/qibotn/backends/quimb.py

@@ -1,21 +1,17 @@
-import re
+from collections import Counter
-import warnings
-from collections import Counter, defaultdict
 
 import numpy as np
-import jax
+
-import jax.numpy as jnp
+from qibo.models import Circuit
-import quimb as qu
-import quimb.tensor as qtn
 from qibo.backends import NumpyBackend
 from qibo.config import raise_error
-from qibo.result import QuantumState
-
 from qibotn.backends.abstract import QibotnBackend
 from qibotn.result import TensorNetworkResult
-
 from qibo.gates.abstract import ParametrizedGate
 
+import quimb as qu
+import quimb.tensor as qtn
+
 
 GATE_MAP = {
         "h": "H",
@@ -23,20 +19,13 @@ GATE_MAP = {
         "y": "Y",
         "z": "Z",
         "s": "S",
-
-        "sdg": "SDG",
         "t": "T",
-        "tdg": "TDG",
-        "sx": "SX",
-        "sxdg": "SXDG",
 
         "rx": "RX",
         "ry": "RY",
         "rz": "RZ",
 
-        "u1": "U1",
+        "u3": "U3", # TODO: check
-        "u2": "U2",
-        "u3": "U3",
 
         "cx": "CX",
         "cnot": "CNOT",
@@ -45,19 +34,18 @@ GATE_MAP = {
 
         "iswap": "ISWAP",
         "swap": "SWAP",
+
         "ccx": "CCX",
-        "toffoli": "CCX",
+        "ccy": "CCY",
         "ccz": "CCZ",
+
+        "toffoli": "TOFFOLI",
         "cswap": "CSWAP",
-        "fredkin": "CSWAP",
+        "fredkin": "FREDKIN",
-        "crx": "CRX",
+
-        "cry": "CRY",
+        "fsim": "fsim",
-        "crz": "CRZ",
+
-        "fsim": "FSIM",
+        "m": "measure"
-        "rxx": "RXX",
-        "ryy": "RYY",
-        "rzz": "RZZ",
-        "m": None,  # measurement, skip
     }
 
 class QuimbBackend(QibotnBackend, NumpyBackend):
@@ -73,7 +61,7 @@ class QuimbBackend(QibotnBackend, NumpyBackend):
 
     def configure_tn_simulation(
         self,
-        ansatz: str = "any",
+        ansatz: str = None,
         max_bond_dimension: int = 10,
         n_most_frequent_states: int = 100,
     ):
@@ -82,9 +70,9 @@ class QuimbBackend(QibotnBackend, NumpyBackend):
 
         Args:
             ansatz : str, optional
-                The tensor network ansatz to use. Currently, only "MPS" or "any" is supported. In the second case
+                The tensor network ansatz to use. Default is `None` and, in this case, a
-                the generic Circuit Quimb class is used.
+                generic Circuit Quimb class is used.            
-            max_bond_dimension : int, optional
+                max_bond_dimension : int, optional
                 The maximum bond dimension for the MPS ansatz. Default is 10.
 
         Notes:
@@ -95,20 +83,32 @@ class QuimbBackend(QibotnBackend, NumpyBackend):
         self.max_bond_dimension = max_bond_dimension
         self.n_most_frequent_states = n_most_frequent_states
 
-    def setup_backend_specifics(self, qimb_backend="numpy", optimizer="auto-hq"):
+    def setup_backend_specifics(self, qimb_backend="numpy", contractions_optimizer="auto-hq"):
         """Setup backend specifics.
         Args:
             qimb_backend: str
                 The backend to use for the quimb tensor network simulation.
-            optimizer: str, optional
+            contractions_optimizer: str, optional
-                The optimizer to use for the quimb tensor network simulation.
+                The contractions_optimizer to use for the quimb tensor network simulation.
         """
+        if qimb_backend == "jax":
+            import jax.numpy as jnp
+            self.np = jnp
+        elif qimb_backend == "numpy":
+            import numpy as np
+            self.np = np
+        elif qimb_backend == "torch":
+            import torch
+            self.np = torch
+        else:
+            raise_error(ValueError, f"Unsupported quimb backend: {qimb_backend}")
+
         self.backend = qimb_backend
-        self.optimizer = optimizer
+        self.contractions_optimizer = contractions_optimizer
 
     def execute_circuit(
         self,
-        circuit,
+        circuit: Circuit,
         initial_state=None,
         nshots=None,
         return_array=False,
@@ -180,7 +180,7 @@ class QuimbBackend(QibotnBackend, NumpyBackend):
             measured_probabilities = None
 
         statevector = (
-            circ_quimb.to_dense(backend=self.backend, optimize=self.optimizer)
+            circ_quimb.to_dense(backend=self.backend, optimize=self.contractions_optimizer)
             if return_array
             else None
         )
@@ -192,193 +192,53 @@ class QuimbBackend(QibotnBackend, NumpyBackend):
             prob_type="default",
             statevector=statevector,
         )
-  
+
-    def expectation(self, circuit, observable):
+    def expectation(self, circuit: Circuit, operators_list: list[str], sites_list: list[str], coeffs_list: list[str]):
         """
-        Compute the expectation value of a Qibo-friendly ``observable`` on the Tensor Network constructed from a Qibo ``circuit``.
+        Compute the expectation value of a symbolic Hamiltonian on a quantum circuit using tensor network contraction.
-        This method takes a Qibo-style symbolic Hamiltonian (e.g., `X(0)*Z(1) + 2.0*Y(2)*Z(0)`)
+        This method takes a Qibo circuit, converts it to a Quimb tensor network circuit, and evaluates the expectation value
-        as the observable, converts it into a Quimb observable and computes its expectation
+        of a Hamiltonian specified by three lists of strings: operators, sites, and coefficients.
-        value using the provided circuit.
+        The expectation value is computed by summing the contributions from each term in the Hamiltonian, where each term's
-
+        expectation is calculated using Quimb's `local_expectation` function.
-        Args:
+        
-            circuit: A Qibo quantum circuit object on which the expectation value
+        Parameters
-                is computed.
+        ----------
-            observable: The observable whose expectation value we want to compute.
+        circuit : qibo.models.Circuit
-                This must be provided in the symbolic Hamiltonian form supported by Qibo
+            The quantum circuit to evaluate, provided as a Qibo circuit object.
-                (e.g., `X(0)*Y(1)` or `Z(0)*Z(1) + 1.5*Y(2)`).
+        operators_list : list of str
-
+            List of operator strings representing the symbolic Hamiltonian terms.
-        Returns:
+        sites_list : list of str
-            float: The expectation value (real part).
+            List of strings, each specifying the qubits (sites) the corresponding operator acts on.
+        coeffs_list : list of str
+            List of strings representing the coefficients for each Hamiltonian term.
+        Returns
+        -------
+        float
+            The real part of the expectation value of the Hamiltonian on the given circuit state.
         """
-
+        quimb_circuit = self._qibo_circuit_to_quimb(circuit, quimb_circuit_type=qtn.Circuit)
-        '''Convert Qibo observables to Quimb'''
+            
-        operators_list, sites_list, coeffs_list = self._qiboobs_to_quimbobs(observable)
-
-        '''Convert Qibo circuit to Quimb circuit'''
-        parameters = circuit.get_parameters()
-        quimb_circuit = self._qibo_circuit_to_quimb(
-            circuit, quimb_circuit_type=qtn.Circuit, to_backend=jnp.array, convert_eager=True
-        )
-        quimb_parameters = {
-            key: jnp.asarray(parameters[i]) for i, key in enumerate(quimb_circuit.get_params().keys())
-        }
-        quimb_circuit.set_params(quimb_parameters)
-
-        '''Compute expectation value'''
         expectation_value = 0.0
-        for ops, sites, coeffs in zip(operators_list, sites_list, coeffs_list):
+        for opstr, sitesstr, coeffstr in zip(operators_list, sites_list, coeffs_list):
+            
+            ops = self._string_to_quimb_operator(opstr)
+            coeff = self._parse_coefficient(coeffstr)
+            sites = tuple(int(q) for q in sitesstr)
+            
             exp_values = quimb_circuit.local_expectation(
                 ops,
                 where=sites,
                 backend=self.backend,
-                optimize=self.optimizer
+                optimize=self.contractions_optimizer
             )
-            expectation_value = expectation_value + coeffs * exp_values
+
+            expectation_value = expectation_value + coeff * exp_values
         
-        return jnp.real(expectation_value)
+        return np.real(expectation_value)
-
-    def expectation_old(self, circuit, observable):
-        """Compute the expectation value of a Qibo-friendly ``observable`` on the Tensor Network constructed from a Qibo ``circuit``.
-
-        This method takes a Qibo-style symbolic Hamiltonian (e.g., `X(0)*Z(1) + 2.0*Y(2)*Z(0)`)
-        as the observable, converts it into a Quimb observable and computes its expectation
-        value using the provided circuit. In case of multiple terms on the same group of qubits, they can be computed in a single contraction.
-        A grouping procedure is applied to optimize the number of contractions performed.
-
-        Args:
-            circuit: A Qibo quantum circuit object on which the expectation value
-                is computed.
-            observable: The observable whose expectation value we want to compute.
-                This must be provided in the symbolic Hamiltonian form supported by Qibo
-                (e.g., `X(0)*Y(1)` or `Z(0)*Z(1) + 1.5*Y(2)`).
-
-        Returns:
-            float: The expectation value (real part).
-        """
-
-        # Map the Qibo observable to Quimb operators and group local operators on the same sites
-        # for computing them in a single contraction. This does not work with CircuitMPS for some now
-        # for Quimb 1.11.1
-        operators_list, sites_list, coeffs_list = self._qiboobs_to_quimbobs(observable)
-        sites_list_grouped, operators_list_grouped, coeffs_list_grouped = (
-            self._group_by_tuples(sites_list, operators_list, coeffs_list)
-        )
-
-        if self.ansatz == "MPS":
-            if len(sites_list) - len(sites_list_grouped) > 10:
-                warnings.warn(
-                    "More than 10 local operators on the same sites are not being grouped as this is not compatible with CircuitMPS. Expected value computation can be more efficient without an MPS ansatz."
-                )
-            circ_ansatz = qtn.circuit.CircuitMPS
-            circ = circ_ansatz.from_openqasm2_str(circuit.to_qasm())
-            expectation_value = 0.0
-            for ops, sites, coeffs in zip(operators_list, sites_list, coeffs_list):
-                exp_values = circ.local_expectation(
-                    ops, where=sites, backend=self.backend, optimize=self.optimizer
-                )
-                expectation_value += np.dot(coeffs, exp_values)
-            return np.real(expectation_value)
-
-        else:
-            circ_ansatz = qtn.circuit.Circuit
-            circ = circ_ansatz.from_openqasm2_str(circuit.to_qasm())
-
-            expectation_value = 0.0
-            for ops, sites, coeffs in zip(
-                operators_list_grouped, sites_list_grouped, coeffs_list_grouped
-            ):
-                exp_values = circ.local_expectation(
-                    ops, where=sites, backend=self.backend, optimize=self.optimizer
-                )
-                expectation_value += np.dot(coeffs, exp_values)
-            return np.real(expectation_value)
-
-    def _qiboobs_to_quimbobs(self, hamiltonian):
-        """
-        Convert a Qibo SymbolicHamiltonian into a Quimb-compatible decomposition.
-
-        Returns three lists:
-        - operators_list: Quimb operators (tensor products of Pauli matrices).
-        - sites_list: tuples of qubit indices the operators act on.
-        - coeffs_list: coefficients for each term.
-        """
-
-        factor_pattern = re.compile(r"([^\d]+)(\d+)")
-
-        operators_list = []
-        sites_list = []
-        coeffs_list = []
-
-        for term in hamiltonian.terms:
-            coeff = term.coefficient
-            term_ops = []
-            term_sites = []
-
-            for factor in term.factors:
-                match = factor_pattern.match(str(factor))
-                if not match:
-                    raise ValueError(
-                        f"Factor '{str(factor)}' does not match the expected format."
-                    )
-
-                operator_name = match.group(1)
-                qubit_index = int(match.group(2))
-
-                # Build the single-qubit operator
-                if operator_name not in {"X", "Y", "Z", "I"}:
-                    raise ValueError(f"Unsupported operator {operator_name}")
-                op = qu.pauli(operator_name)
-
-                term_ops.append(op)
-                term_sites.append(qubit_index)
-
-            # Build the tensor product if more than one factor
-            if term_ops:
-                full_op = term_ops[0]
-                for op in term_ops[1:]:
-                    full_op = full_op & op
-            else:
-                # Identity term (just coefficient)
-                full_op = qu.eye(2)
-
-            operators_list.append(full_op)
-            sites_list.append(tuple(term_sites))
-            coeffs_list.append(coeff)
-
-        return operators_list, sites_list, coeffs_list
-
-    def _group_by_tuples(self, A, B, C):
-        """
-        Groups the elements of B and C by the unique tuples in A.
-
-        Parameters:
-            A (list of tuples): key tuples (can contain duplicates)
-            B (list): values aligned with A
-            C (list): values aligned with A
-
-        Returns:
-            (A_new, B_new, C_new):
-                A_new: list of unique tuples
-                B_new: list of lists of grouped values from B
-                C_new: list of lists of grouped values from C
-        """
-
-        grouped_B = defaultdict(list)
-        grouped_C = defaultdict(list)
-
-        for a, b, c in zip(A, B, C):
-            grouped_B[a].append(b)
-            grouped_C[a].append(c)
-
-        A_new = list(grouped_B.keys())
-        B_new = list(grouped_B.values())
-        C_new = list(grouped_C.values())
-
-        return A_new, B_new, C_new
 
     def _qibo_circuit_to_quimb(self, qibo_circ, quimb_circuit_type=qtn.Circuit, **circuit_kwargs):
         """
-        Convert a Qibo Circuit to a Quimb Circuit.
+        Convert a Qibo Circuit to a Quimb Circuit. Measurement gates are ignored. If are given gates not supported by Quimb, an error is raised.
 
         Parameters
         ----------
@@ -400,29 +260,71 @@ class QuimbBackend(QibotnBackend, NumpyBackend):
         for gate in qibo_circ.queue:
             gname = getattr(gate, "name", None)
             qname = GATE_MAP.get(gname, None)
+            if qname == "measure":
+                continue
             if qname is None:
-                continue  # skip measurements and unknown gates
+                raise_error(ValueError, f"Gate {gname} not supported in Quimb backend.")
 
             params = getattr(gate, "parameters", ())
             qubits = getattr(gate, "qubits", ())
             
-            # Check if the gate is parametrized
             is_parametrized = (
                 isinstance(gate, ParametrizedGate)
                 and getattr(gate, "trainable", True)
-            )
+                )
-            # import pdb; pdb.set_trace()
             if is_parametrized:
                 circ.apply_gate(
                     qname,
                     *params,
                     *qubits,
-                    parametrized= is_parametrized
+                    parametrized=is_parametrized
                 )
             else:
                 circ.apply_gate(
                     qname,
                     *params,
                     *qubits,
-                                )
+                )
         return circ
+
+    def _parse_coefficient(self, s):
+        """Parse a coefficient from string to float, int, or complex.
+        Args:
+            s: str
+                The string representation of the coefficient.
+        Returns:
+            The coefficient as float, int, or complex.
+        """ 
+        try:
+            return float(s)
+        except ValueError:
+            pass
+        if s == "j":
+            return 1j
+        elif s == "-j":
+            return -1j
+        else:
+            try:
+                return complex(s)
+            except ValueError:
+                raise ValueError(f"Cannot parse coefficient: {s}")
+
+    def _string_to_quimb_operator(self, op_str):
+        """
+        Convert a Pauli string (e.g. 'xzy') to a Quimb operator using '&' chaining.
+
+        Parameters
+        ----------
+        op_str : str
+            A string like 'xzy', where each character is one of 'x', 'y', 'z', 'i'.
+
+        Returns
+        -------
+        qu_op : quimb.Qarray
+            The corresponding Quimb operator.
+        """
+        op_str = op_str.lower()
+        op = qu.pauli(op_str[0])
+        for c in op_str[1:]:
+            op = op & qu.pauli(c)
+        return op