Merge pull request #170 from Zhaoyilunnn/master

Finalize algorithms module

quafu/algorithms/ansatz.py

@@ -17,6 +17,7 @@
 
 import numpy as np
 from quafu.circuits.quantum_circuit import QuantumCircuit
+from quafu.elements import Parameter
 from quafu.synthesis.evolution import ProductFormula
 
 from .hamiltonian import Hamiltonian
@@ -52,8 +53,10 @@ def __init__(self, hamiltonian: Hamiltonian, num_qubits: int, num_layers: int =
         self._evol = ProductFormula()
 
         # Initialize parameters
-        self._beta = np.zeros(num_layers)
-        self._gamma = np.zeros(num_layers)
+        self._beta = np.array([Parameter(f"beta_{i}", 0.0) for i in range(num_layers)])
+        self._gamma = np.array(
+            [Parameter(f"gamma_{i}", 0.0) for i in range(num_layers)]
+        )
 
         # Build circuit structure
         super().__init__(num_qubits)
@@ -122,7 +125,10 @@ def __init__(self, num_qubits: int, layer: int):
             layer: Number of layers.
         """
         self._layer = layer
-        self._theta = np.zeros((layer + 1, num_qubits))
+        self._theta = np.array(
+            [Parameter(f"theta_{i}", 0.0) for i in range((layer + 1) * num_qubits)]
+        )
+        self._theta = np.reshape(self._theta, (layer + 1, num_qubits))
         super().__init__(num_qubits)
 
     def _build(self):
@@ -153,25 +159,34 @@ def __init__(
         self._transformer = InterfaceProvider.get(interface)
         self._layers = layers
 
-        # FIXME(zhaoyilun): don't use this default value
-        self._weights = np.empty((1, 1))
+        self._weights = None
 
         self._backend = backend
         super().__init__(num_qubits)
 
-    def __call__(self, features):
+    def __call__(self, inputs):
         """Compute outputs of QNN given input features"""
         from .estimator import Estimator
 
         estimator = Estimator(self, backend=self._backend)
-        return self._transformer.execute(self, features, estimator=estimator)
+        return self._transformer.execute(self, inputs, estimator=estimator)
 
     def _build(self):
         """Essentially initialize weights using transformer"""
         self.add_gates(self._layers)
 
-        self._weights = self._transformer.init_weights((1, self.num_parameters))
+        self._weights = self._transformer.init_weights((1, self.num_tunable_parameters))
 
     @property
     def weights(self):
         return self._weights
+
+    @property
+    def num_tunable_parameters(self):
+        num_tunable_params = 0
+        for g in self.gates:
+            paras = g.paras
+            for p in paras:
+                if hasattr(p, "tunable") and p.tunable:
+                    num_tunable_params += 1
+        return num_tunable_params

quafu/algorithms/estimator.py

@@ -13,15 +13,16 @@
 # limitations under the License.
 """Pre-build wrapper to calculate expectation value"""
 from typing import List, Optional
+
 from ..circuits.quantum_circuit import QuantumCircuit
+from ..simulators import simulate
 from ..tasks.tasks import Task
 from .hamiltonian import Hamiltonian
-from ..simulators import simulate
 
 
 def execute_circuit(circ: QuantumCircuit, observables: Hamiltonian):
     """Execute circuit on quafu simulator"""
-    sim_res = simulate(circ, hamiltonian= observables)
+    sim_res = simulate(circ, hamiltonian=observables)
     expectations = sim_res["pauli_expects"]
     return sum(expectations)
 
@@ -44,6 +45,7 @@ def __init__(
             task_options: options to config a task instance
         """
         self._circ = circ
+        self._circ.get_parameter_grads()  # parameter shift currently requires calling this for initialization
         self._backend = backend
         self._task = None
         if backend != "sim":
@@ -85,7 +87,7 @@ def run(self, observables: Hamiltonian, params: List[float]):
             Expectation value
         """
         if params is not None:
-            self._circ.update_params(params)
+            self._circ._update_params(params)
 
         if self._backend == "sim":
             return self._run_simulation(observables)

quafu/algorithms/gradients/__init__.py

@@ -12,5 +12,6 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+from .gradiant import grad_adjoint, grad_finit_diff, grad_para_shift
 from .param_shift import ParamShift
 from .vjp import compute_vjp, jacobian, run_circ

quafu/algorithms/gradient.py → quafu/algorithms/gradients/gradiant.py

@@ -1,45 +1,51 @@
-from ..circuits.quantum_circuit import QuantumCircuit
-from ..simulators.simulator import SVSimulator
-from ..elements import Parameter, ParameterExpression
 import numpy as np
-from ..exceptions import CircuitError
-from ..elements.matrices import XMatrix, YMatrix, ZMatrix
-from ..elements import QuantumGate, ControlledGate
+from quafu.circuits.quantum_circuit import QuantumCircuit
+from quafu.elements import ControlledGate, Parameter, ParameterExpression, QuantumGate
+from quafu.elements.matrices import XMatrix, YMatrix, ZMatrix
+from quafu.exceptions import CircuitError
+from quafu.simulators.simulator import SVSimulator
+
 
 def assemble_grads(para_grads, gate_grads):
     grads = []
     for var in para_grads:
         grad_p = para_grads[var]
-        fullgrad = 0.
+        fullgrad = 0.0
         for pos_g in grad_p:
-            pos, gp = pos_g 
+            pos, gp = pos_g
             gg = gate_grads[pos[0]][pos[1]]
             fullgrad += gg * gp
         grads.append(fullgrad)
-    
+
     return grads
 
-def grad_para_shift(qc:QuantumCircuit, hamiltonian, backend=SVSimulator()):
+
+def grad_para_shift(qc: QuantumCircuit, hamiltonian, backend=SVSimulator()):
     """
     Parameter shift gradients. Each gate must have one parameter
     """
     para_grads = qc._calc_parameter_grads()
-    gate_grads= [[] for _ in qc.gates]
+    gate_grads = [[] for _ in qc.gates]
 
     for i, op in enumerate(qc.gates):
         if len(op.paras) > 0:
-            if isinstance(op.paras[0], Parameter) or isinstance(op.paras[0],ParameterExpression):
+            if isinstance(op.paras[0], Parameter) or isinstance(
+                op.paras[0], ParameterExpression
+            ):
                 if op.name not in ["RX", "RY", "RZ"]:
-                    raise CircuitError("It seems the circuit can not apply parameter-shift rule to calculate gradient.You may need compile the circuit first")
-                op.paras[0] = op.paras[0] + np.pi/2
+                    raise CircuitError(
+                        "It seems the circuit can not apply parameter-shift rule to calculate gradient. You may need compile the circuit first"
+                    )
+                op.paras[0] = op.paras[0] + np.pi / 2
                 res1 = sum(backend.run(qc, hamiltonian=hamiltonian)["pauli_expects"])
                 op.paras[0] = op.paras[0] - np.pi
                 res2 = sum(backend.run(qc, hamiltonian=hamiltonian)["pauli_expects"])
                 op.paras[0]._undo(2)
-                gate_grads[i].append((res1 - res2) / 2.)
-        
+                gate_grads[i].append((res1 - res2) / 2.0)
+
     return assemble_grads(para_grads, gate_grads)
 
+
 def grad_finit_diff(qc, hamiltonian, backend=SVSimulator()):
     variables = qc.variables
     grads = []
@@ -50,7 +56,7 @@ def grad_finit_diff(qc, hamiltonian, backend=SVSimulator()):
         res2 = sum(backend.run(qc, hamiltonian=hamiltonian)["pauli_expects"])
         v.value += 1e-10
         grads.append((res1 - res2) / (2 * 1e-10))
-    
+
     return grads
 
 
@@ -60,37 +66,37 @@ def grad_gate(op):
     """
     if isinstance(op, ControlledGate):
         if op._targ_name == "RX":
-            circ = QuantumCircuit(max(op.pos)+1)
+            circ = QuantumCircuit(max(op.pos) + 1)
             deriv_mat = -0.5j * XMatrix @ op._get_targ_matrix()
             circ << QuantumGate("dRX", op.targs, [], deriv_mat)
             cdim = 1 << (len(op.ctrls))
             proj_mat = np.zeros((cdim, cdim))
-            proj_mat[cdim-1, cdim-1] = 1.
+            proj_mat[cdim - 1, cdim - 1] = 1.0
             circ << QuantumGate("projCtrl", op.ctrls, [], proj_mat)
             return circ.wrap()
-        
+
         elif op._targ_name == "RY":
-            circ = QuantumCircuit(max(op.pos)+1)
+            circ = QuantumCircuit(max(op.pos) + 1)
             deriv_mat = -0.5j * YMatrix @ op._get_targ_matrix()
             circ << QuantumGate("dRY", op.targs, [], deriv_mat)
             cdim = 1 << (len(op.ctrls))
             proj_mat = np.zeros((cdim, cdim))
-            proj_mat[cdim-1, cdim-1] = 1.
+            proj_mat[cdim - 1, cdim - 1] = 1.0
             circ << QuantumGate("projCtrl", op.ctrls, [], proj_mat)
             return circ.wrap()
-        
+
         elif op._targ_name == "RZ":
-            circ = QuantumCircuit(max(op.pos)+1)
+            circ = QuantumCircuit(max(op.pos) + 1)
             deriv_mat = -0.5j * ZMatrix @ op._get_targ_matrix()
             circ << QuantumGate("dRZ", op.targs, [], deriv_mat)
             cdim = 1 << (len(op.ctrls))
             proj_mat = np.zeros((cdim, cdim))
-            proj_mat[cdim-1, cdim-1] = 1.
+            proj_mat[cdim - 1, cdim - 1] = 1.0
             circ << QuantumGate("projCtrl", op.ctrls, [], proj_mat)
             return circ.wrap()
         else:
             raise NotImplementedError
-        
+
     else:
         if op.name == "RX":
             deriv_mat = -0.5j * XMatrix @ op.matrix
@@ -103,30 +109,37 @@ def grad_gate(op):
             return QuantumGate("dRZ", op.pos, [], deriv_mat)
         else:
             raise NotImplementedError
-
+
+
 def grad_adjoint(qc, hamiltonian, psi_in=np.array([], dtype=complex)):
     """
     Reverse mode gradient: arXiv:2009.02823
     """
     para_grads = qc._calc_parameter_grads()
     backend = SVSimulator()
-    lam = backend.run(qc, psi = psi_in)["statevector"]
+    lam = backend.run(qc, psi=psi_in)["statevector"]
     phi = np.copy(lam)
     lam = backend._apply_hamil(hamiltonian, lam)
     begin = 0
     end = len(qc.gates)
-    gate_grads= [[] for _ in range(end)]
+    gate_grads = [[] for _ in range(end)]
     for i, op in enumerate(qc.gates):
-        if len(op.paras) > 0 and (isinstance(op.paras[0], Parameter) or isinstance(op.paras[0],ParameterExpression)):
+        if len(op.paras) > 0 and (
+            isinstance(op.paras[0], Parameter)
+            or isinstance(op.paras[0], ParameterExpression)
+        ):
             begin = i
             break
-    
+
     for i in range(begin, end)[::-1]:
         op = qc.gates[i]
         phi = backend._apply_op(op.dagger(), phi)
-        if len(op.paras) > 0 and (isinstance(op.paras[0], Parameter) or isinstance(op.paras[0],ParameterExpression)):
-                mu = np.copy(phi)
-                mu = backend._apply_op(grad_gate(op), mu)
-                gate_grads[i].append(np.real(2. * np.inner(lam.conj(), mu)))
+        if len(op.paras) > 0 and (
+            isinstance(op.paras[0], Parameter)
+            or isinstance(op.paras[0], ParameterExpression)
+        ):
+            mu = np.copy(phi)
+            mu = backend._apply_op(grad_gate(op), mu)
+            gate_grads[i].append(np.real(2.0 * np.inner(lam.conj(), mu)))
         lam = backend._apply_op(op.dagger(), lam)
-    return assemble_grads(para_grads, gate_grads)
+    return assemble_grads(para_grads, gate_grads)

quafu/algorithms/gradients/param_shift.py

@@ -19,6 +19,7 @@
 
 from ..estimator import Estimator
 from ..hamiltonian import Hamiltonian
+from .gradiant import grad_para_shift
 
 
 class ParamShift:
@@ -34,7 +35,7 @@ def __call__(self, obs: Hamiltonian, params: List[float]):
             estimator (Estimator): estimator to calculate expectation values
             params (List[float]): params to optimize
         """
-        return self.grad(obs, params)
+        return self.new_grad(obs, params)
 
     def _gen_param_shift_vals(self, params):
         """Given a param list with n values, replicate to 2*n param list"""
@@ -45,6 +46,7 @@ def _gen_param_shift_vals(self, params):
         minus_params = params - offsets * np.pi / 2
         return plus_params.tolist() + minus_params.tolist()
 
+    # TODO: delete after 0.4.1
     def grad(self, obs: Hamiltonian, params: List[float]):
         """grad.
 
@@ -61,3 +63,12 @@ def grad(self, obs: Hamiltonian, params: List[float]):
         num_shift_params = len(res)
         grads = (res[: num_shift_params // 2] - res[num_shift_params // 2 :]) / 2
         return grads
+
+    def new_grad(self, obs: Hamiltonian, params: List[float]):
+        """Calculate the gradients of given the circuit based on the parameter shift rule
+        Args:
+            obs (Hamiltonian): observables for measurement.
+            params (List[float]): parameters to apply to the circuit.
+        """
+        self._est._circ._update_params(params)
+        return grad_para_shift(self._est._circ, obs)

quafu/algorithms/interface/torch.py

@@ -17,7 +17,9 @@
 
 import numpy as np
 import torch
+from quafu.algorithms.ansatz import QuantumNeuralNetwork
 from quafu.algorithms.estimator import Estimator
+from torch import nn
 
 from quafu import QuantumCircuit
 
@@ -56,14 +58,7 @@ def execute(
             "estimator": estimator,
         }
 
-        if method == "external":
-            return ExecuteCircuits.apply(parameters, kwargs)
-        if method == "internal":
-            from ..ansatz import QuantumNeuralNetwork
-
-            assert isinstance(circ, QuantumNeuralNetwork)
-            return ExecuteCircuits.apply(circ.weights, kwargs)
-        raise NotImplementedError(f"Unsupported execution method: {method}")
+        return ExecuteCircuits.apply(parameters, kwargs)
 
 
 class ExecuteCircuits(torch.autograd.Function):
@@ -91,3 +86,42 @@ def backward(ctx, grad_out):
         vjp = compute_vjp(jac, grad_out.numpy())
         vjp = torch.from_numpy(vjp)
         return vjp, None
+
+
+class ModuleWrapper(nn.Module):
+    """
+    A wrapper class to transform quafu circuit to a torch module
+    """
+
+    def __init__(self, qnn: QuantumNeuralNetwork):
+        """
+        Initialization of quafu torch module
+
+        Args:
+            circ (QuantumCircuit): the original parameterized quantum circuit
+        """
+        super().__init__()
+        self._qnn = qnn
+        if qnn.weights is not None:
+            self.weights = nn.parameter.Parameter(qnn.weights)
+        else:
+            self.weights = None
+
+    def forward(self, inputs: torch.Tensor):
+        """
+        Args:
+            inputs (torch.Tensor): raw input data or output from previous
+                classical/quantum layers.
+        """
+        # if weights are not empty, it will be combined with inputs to form
+        # the complete parameter vector and feed to the quantum circuit
+        bsz, _ = inputs.shape  # FIXME: currently we assume 2-D inputs
+
+        # use the last dimension since it is currently initialized as (1, D)
+        if self.weights is not None:
+            weight_dim = self.weights.size(-1)
+            weights_expanded = self.weights.expand(bsz, weight_dim)
+            inputs_to_circ = torch.cat((inputs, weights_expanded), dim=1)
+        else:
+            inputs_to_circ = inputs
+        return self._qnn(inputs_to_circ)