feat:add parameter

quafu/__init__.py

@@ -1,4 +1,5 @@
 from .circuits.quantum_circuit import QuantumCircuit
+from .algorithms.hamiltonian import Hamiltonian
 from .circuits.quantum_register import QuantumRegister, Qubit
 from .results.results import ExecResult, SimuResult
 from .simulators.simulator import simulate
@@ -9,6 +10,7 @@
     "QuantumCircuit",
     "QuantumRegister",
     "Qubit",
+    "Hamiltonian",
     "ExecResult",
     "Task",
     "User",

quafu/algorithms/__init__.py

@@ -5,4 +5,5 @@
 from .hamiltonian import Hamiltonian
 from .templates.amplitude import AmplitudeEmbedding
 from .templates.angle import AngleEmbedding
+from .templates.amplitude import AmplitudeEmbedding
 from .templates.basic_entangle import BasicEntangleLayers

quafu/algorithms/ansatz.py

@@ -21,7 +21,6 @@
 
 from .hamiltonian import Hamiltonian
 from .interface_provider import InterfaceProvider
-from .templates import AngleEmbedding
 
 
 class Ansatz(QuantumCircuit, ABC):
@@ -44,10 +43,11 @@ def _build(self):
 class QAOAAnsatz(Ansatz):
     """QAOA Ansatz"""
 
-    def __init__(self, hamiltonian: Hamiltonian, num_layers: int = 1):
+    def __init__(self, hamiltonian: Hamiltonian, num_qubits: int, num_layers: int = 1):
         """Instantiate a QAOAAnsatz"""
-        self._pauli_list = hamiltonian.pauli_list
-        self._coeffs = hamiltonian.coeffs
+        # self._pauli_list = hamiltonian.pauli_list
+        # self._coeffs = hamiltonian.coeffs
+        self._h = hamiltonian
         self._num_layers = num_layers
         self._evol = ProductFormula()
 
@@ -56,7 +56,6 @@ def __init__(self, hamiltonian: Hamiltonian, num_layers: int = 1):
         self._gamma = np.zeros(num_layers)
 
         # Build circuit structure
-        num_qubits = len(self._pauli_list[0])
         super().__init__(num_qubits)
 
     @property
@@ -68,7 +67,7 @@ def parameters(self):
         """Return complete parameters of the circuit"""
         paras_list = []
         for layer in range(self._num_layers):
-            for _ in range(len(self._pauli_list)):
+            for _ in range(len(self._h.paulis)):
                 paras_list.append(self._gamma[layer])
             for _ in range(self.num):
                 paras_list.append(self._beta[layer])
@@ -86,8 +85,8 @@ def _build(self):
 
         for layer in range(self._num_layers):
             # Add H_D layer
-            for pauli_str in self._pauli_list:
-                gate_list = self._evol.evol(pauli_str, self._gamma[layer])
+            for pauli in self._h.paulis:
+                gate_list = self._evol.evol(pauli, self._gamma[layer])
                 for g in gate_list:
                     self.add_ins(g)
 
@@ -156,10 +155,6 @@ def __init__(self, num_qubits: int, layers: List[Any], interface="torch"):
         self._weights = np.empty((1, 1))
         super().__init__(num_qubits)
 
-    def __call__(self, features):
-        """Compute outputs of QNN given input features"""
-        return self._transformer.execute(self, features)
-
     def _build(self):
         """Essentially initialize weights using transformer"""
         for layer in self._layers:

quafu/algorithms/estimator.py

@@ -13,24 +13,17 @@
 # limitations under the License.
 """Pre-build wrapper to calculate expectation value"""
 from typing import List, Optional
-
-import numpy as np
-from quafu.algorithms.hamiltonian import Hamiltonian
-from quafu.simulators.simulator import simulate
-from quafu.tasks.tasks import Task
-
-from quafu import QuantumCircuit
+from ..circuits.quantum_circuit import QuantumCircuit
+from ..tasks.tasks import Task
+from .hamiltonian import Hamiltonian
+from ..simulators.simulator import simulate
 
 
 def execute_circuit(circ: QuantumCircuit, observables: Hamiltonian):
     """Execute circuit on quafu simulator"""
-    sim_state = simulate(circ, output="state_vector").get_statevector()
-
-    expectation = np.matmul(
-        np.matmul(sim_state.conj().T, observables.get_matrix()), sim_state
-    ).real
-
-    return expectation
+    sim_res = simulate(circ, output="state_vector")
+    expectations = sim_res.expect_paulis(observables)
+    return sum(expectations)
 
 
 class Estimator:
@@ -91,11 +84,6 @@ def run(self, observables: Hamiltonian, params: List[float]):
         Returns:
             Expectation value
         """
-        if observables.num_qubits != self._circ.num:
-            raise ValueError(
-                "The number of qubits in the observables does not match the circuit"
-            )
-
         if params is not None:
             self._circ.update_params(params)
 

quafu/algorithms/gradients/vjp.py

@@ -24,9 +24,8 @@
 
 def _generate_expval_z(num_qubits: int):
     obs_list = []
-    base_pauli = "I" * num_qubits
     for i in range(num_qubits):
-        pauli = base_pauli[:i] + "Z" + base_pauli[i + 1 :]
+        pauli = "Z" + str(i)
         obs_list.append(Hamiltonian.from_pauli_list([(pauli, 1)]))
     return obs_list
 

quafu/algorithms/hamiltonian.py

@@ -11,97 +11,189 @@
 # 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.
-"""Quafu Hamiltonian class"""
 
-from __future__ import annotations
+from typing import Iterable
+import scipy.sparse as sp
+import numpy as np
 
-from collections.abc import Iterable
+from quafu.exceptions.quafu_error import QuafuError
 
-import numpy as np
-from quafu.elements.matrices import IdMatrix, XMatrix, YMatrix, ZMatrix
-from quafu.exceptions import QuafuError
+IMat = sp.coo_matrix(np.array([[1.0, 0.0], [0.0, 1.0]], dtype=complex))
 
-PAULI_MAT = {"I": IdMatrix, "X": XMatrix, "Y": YMatrix, "Z": ZMatrix}
+XMat = sp.coo_matrix(np.array([[0.0, 1.0], [1.0, 0.0]], dtype=complex))
 
+YMat = sp.coo_matrix(np.array([[0.0, -1.0j], [1.0j, 0.0]], dtype=complex))
 
-class Hamiltonian:
-    """TODO"""
+ZMat = sp.coo_matrix(np.array([[1.0, 0.0], [0.0, -1.0]], dtype=complex))
 
-    def __init__(self, pauli_str_list: list[str], coeffs: np.ndarray) -> None:
-        """
-        Args:
-            pauli_str_list: List of Pauli strs, e.g., ['IIIZZ', "IIZIZ", ...]
-            coeffs: List of efficients
-        """
-        self._pauli_str_list = pauli_str_list
-        self._coeffs = coeffs
+PauliMats = {"X": XMat, "Y": YMat, "Z": ZMat, "I": IMat}
+
+
+class PauliOp:
+    def __init__(self, paulis: str, coeff: complex = 1.0):
+        paulist = paulis.split(" ")
+        self.paulistr = ""
+        self.pos = []
+        for p in paulist:
+            assert p[0] in "XYZ"
+            self.paulistr += p[0]
+            self.pos.append(int(p[1:]))
+        self.coeff = coeff
+
+    def __repr__(self):
+        repstr = ""
+        if self.coeff != 1.0:
+            repstr = str(self.coeff) + "*"
+        for i in range(len(self.pos)):
+            repstr += self.paulistr[i]
+            repstr += str(self.pos[i])
+            repstr += "*"
+        return repstr[:-1]
+
+    def __str__(self):
+        return self.__repr__()
+
+    def __mul__(self, obj):
+        pass
+
+    def __rmul__(self, obj):
+        pass
+
+    def commutator(self, obj):
+        pass
 
-    @property
-    def num_qubits(self):
-        """Get the number of qubits of the system"""
-        return len(self.pauli_list[0])
+    def get_matrix(self, qnum, big_endian=False):
+        if qnum - 1 < max(self.pos):
+            raise ValueError("The support of the paulis exceed the total qubit number")
 
-    @property
-    def pauli_list(self):
-        """Get pauli string list"""
-        return self._pauli_str_list
+        pos = np.array(self.pos)
+        if not big_endian:
+            pos = qnum - 1 - pos
+        inds = np.argsort(pos)
+        iq = 0
+        ip = 0
+        mat = 1.0
+        while iq < qnum:
+            if ip < len(pos):
+                if iq == pos[inds[ip]]:
+                    opstr = self.paulistr[inds[ip]]
+                    mat = sp.kron(mat, PauliMats[opstr])
+                    iq += 1
+                    ip += 1
+                else:
+                    mat = sp.kron(mat, PauliMats["I"])
+                    iq += 1
+            else:
+                mat = sp.kron(mat, PauliMats["I"])
+                iq += 1
 
-    @property
-    def coeffs(self):
-        """Get coefficients of each pauli string"""
-        return self._coeffs
+        return self.coeff * mat
+
+
+class Hamiltonian:
+    def __init__(self, paulis: list[PauliOp]):
+        self.paulis = paulis
 
     @staticmethod
-    def from_pauli_list(pauli_list: Iterable[tuple[str, complex]]) -> Hamiltonian:
-        """
+    def from_pauli_list(pauli_list: Iterable[tuple[str, complex]]):
+        """Generate Hamiltonian from Pauli string list
+
         Args:
-            pauli: The supported format of pauli list is [(<pauli-str>, <coefficient>)],
-                e.g., [('IIIZZ', 1), ("IIZIZ", 1), ...)], 0th qubit is farthest right
+            pauli_list: e.g., [("Z0 Z1", 1), ("Z0 Z2", 1), ...]
         """
-
         pauli_list = list(pauli_list)
 
         size = len(pauli_list)
         if size == 0:
             raise QuafuError("Pauli list cannot be empty.")
 
-        coeffs = np.zeros(size, dtype=complex)
+        pauli_op_list = [PauliOp(p[0], coeff=p[1]) for p in pauli_list]
+
+        return Hamiltonian(pauli_op_list)
 
-        pauli_str_list = []
-        for i, (pauli_str, coef) in enumerate(pauli_list):
-            pauli_str_list.append(pauli_str)
-            coeffs[i] = coef
+    def __repr__(self):
+        return "+".join([str(pauli) for pauli in self.paulis])
 
-        return Hamiltonian(pauli_str_list, coeffs)
+    def __str__(self):
+        return self.__repr__()
+
+    def get_matrix(self, qnum, big_endian=False):
+        mat = 0.0
+        for pauli in self.paulis:
+            mat += pauli.get_matrix(qnum, big_endian)
+
+        return mat
 
     # TODO(zhaoyilun): delete this in the future
     def to_legacy_quafu_pauli_list(self):
         """Transform to legacy quafu pauli list format,
         this is a temperal function and should be deleted later"""
         res = []
-        for pauli_str in self._pauli_str_list:
-            for i, pauli in enumerate(pauli_str[::-1]):
-                if pauli in ["X", "Y", "Z"]:
-                    res.append([pauli, [i]])
+        for pauli_str in self.paulis:
+            for i, pos in enumerate(pauli_str.pos):
+                res.append([pauli_str.paulistr[i], [pos]])
         return res
 
-    def _get_pauli_mat(self, pauli_str: str):
-        """Calculate the matrix of a pauli string"""
-        mat = None
-        for pauli in pauli_str[::-1]:
-            mat = PAULI_MAT[pauli] if mat is None else np.kron(PAULI_MAT[pauli], mat)
-        return mat
-
-    def matrix_generator(self):
-        """Generating matrix for each Pauli str"""
-        for i, pauli_str in enumerate(self._pauli_str_list):
-            yield self._coeffs[i] * self._get_pauli_mat(pauli_str)
-
-    def get_matrix(self):
-        """Generate matrix of Hamiltonian"""
 
-        dim = 2**self.num_qubits
-        matrix = np.zeros((dim, dim), dtype=complex)
-        for mat in self.matrix_generator():
-            matrix += mat
-        return matrix
+def intersec(a, b):
+    inter = []
+    aind = []
+    bind = []
+    for i in range(len(a)):
+        for j in range(len(b)):
+            if a[i] == b[j]:
+                inter.append(a[i])
+                aind.append(i)
+                bind.append(j)
+
+    return inter, aind, bind
+
+
+def diff(a, b):
+    diff = []
+    aind = []
+    for i in range(len(a)):
+        if a[i] not in b:
+            diff.append(a[i])
+            aind.append(i)
+
+    return diff, aind
+
+
+def merge_paulis(obslist):
+    measure_basis = []
+    targ_basis = []
+    for obs in obslist:
+        if len(measure_basis) == 0:
+            measure_basis.append(obs)
+            targ_basis.append(len(measure_basis) - 1)
+        else:
+            added = 0
+            for mi in range(len(measure_basis)):
+                measure_base = measure_basis[mi]
+                interset, intobsi, intbasei = intersec(obs.pos, measure_base.pos)
+                diffset, diffobsi = diff(obs.pos, measure_base.pos)
+                if not len(interset) == 0:
+                    if all(
+                        np.array(list(obs.paulistr))[intobsi]
+                        == np.array(list(measure_base.paulistr))[intbasei]
+                    ):
+                        measure_base.paulistr += "".join(
+                            np.array(list(obs.paulistr))[diffobsi]
+                        )
+                        measure_base.pos.extend(diffset)
+                        targ_basis.append(mi)
+                        added = 1
+                        break
+                else:
+                    measure_base.paulistr += obs.paulistr
+                    measure_base.pos.extend(obs.pos)
+                    targ_basis.append(mi)
+                    added = 1
+                    break
+
+            if not added:
+                measure_basis.append(obs)
+                targ_basis.append(len(measure_basis) - 1)
+
+    return measure_basis, targ_basis