test all

quafu/__init__.py

@@ -2,7 +2,7 @@
 from .algorithms.hamiltonian import Hamiltonian
 from .circuits.quantum_register import QuantumRegister, Qubit
 from .results.results import ExecResult, SimuResult
-from .simulators.simulator import simulate
+from .simulators import simulate
 from .tasks.tasks import Task
 from .users.userapi import User
 

quafu/algorithms/__init__.py

@@ -1,9 +1,9 @@
 """Algorithm module"""
 
-# from .ansatz import AlterLayeredAnsatz, QAOAAnsatz, QuantumNeuralNetwork
-# from .estimator import Estimator
-# 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
+from .ansatz import AlterLayeredAnsatz, QAOAAnsatz, QuantumNeuralNetwork
+from .estimator import Estimator
+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/estimator.py

@@ -16,13 +16,13 @@
 from ..circuits.quantum_circuit import QuantumCircuit
 from ..tasks.tasks import Task
 from .hamiltonian import Hamiltonian
-from ..simulators.simulator import simulate
+from ..simulators import simulate
 
 
 def execute_circuit(circ: QuantumCircuit, observables: Hamiltonian):
     """Execute circuit on quafu simulator"""
-    sim_res = simulate(circ, output="state_vector")
-    expectations = sim_res.expect_paulis(observables)
+    sim_res = simulate(circ, hamiltonian= observables)
+    expectations = sim_res["pauli_expects"]
     return sum(expectations)
 
 
@@ -67,7 +67,7 @@ def _run_real_machine(self, observables: Hamiltonian):
 
     def _run_simulation(self, observables: Hamiltonian):
         """Run using quafu simulator"""
-        # sim_state = simulate(self._circ, output="state_vector").get_statevector()
+        # sim_state = simulate(self._circ).get_statevector()
         # expectation = np.matmul(
         #     np.matmul(sim_state.conj().T, observables.get_matrix()), sim_state
         # ).real

quafu/algorithms/gradient.py

@@ -0,0 +1,132 @@
+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
+
+def assemble_grads(para_grads, gate_grads):
+    grads = []
+    for var in para_grads:
+        grad_p = para_grads[var]
+        fullgrad = 0.
+        for pos_g in grad_p:
+            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()):
+    """
+    Parameter shift gradients. Each gate must have one parameter
+    """
+    para_grads = qc._calc_parameter_grads()
+    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 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
+                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.)
+
+    return assemble_grads(para_grads, gate_grads)
+
+def grad_finit_diff(qc, hamiltonian, backend=SVSimulator()):
+    variables = qc.variables
+    grads = []
+    for v in variables:
+        v.value += 1e-10
+        res1 = sum(backend.run(qc, hamiltonian=hamiltonian)["pauli_expects"])
+        v.value -= 2 * 1e-10
+        res2 = sum(backend.run(qc, hamiltonian=hamiltonian)["pauli_expects"])
+        v.value += 1e-10
+        grads.append((res1 - res2) / (2 * 1e-10))
+
+    return grads
+
+
+def grad_gate(op):
+    """
+    TODO:support more gates
+    """
+    if isinstance(op, ControlledGate):
+        if op._targ_name == "RX":
+            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.
+            circ << QuantumGate("projCtrl", op.ctrls, [], proj_mat)
+            return circ.wrap()
+
+        elif op._targ_name == "RY":
+            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.
+            circ << QuantumGate("projCtrl", op.ctrls, [], proj_mat)
+            return circ.wrap()
+
+        elif op._targ_name == "RZ":
+            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.
+            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
+            return QuantumGate("dRX", op.pos, [], deriv_mat)
+        elif op.name == "RY":
+            deriv_mat = -0.5j * YMatrix @ op.matrix
+            return QuantumGate("dRY", op.pos, [], deriv_mat)
+        elif op.name == "RZ":
+            deriv_mat = -0.5j * ZMatrix @ op.matrix
+            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"]
+    phi = np.copy(lam)
+    lam = backend._apply_hamil(hamiltonian, lam)
+    begin = 0
+    end = len(qc.gates)
+    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)):
+            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)))
+        lam = backend._apply_op(op.dagger(), lam)
+    return assemble_grads(para_grads, gate_grads)

quafu/algorithms/optimizer.py

@@ -0,0 +1,72 @@
+import numpy as np
+
+def adam(func, x, gradsf, args=(), call_back=None, eta=0.01,beta1=0.9, beta2=0.999, epsilon=1e-8, maxiter=1000, tol=1e-8, verbose=False):
+    n_para = len(x)
+    mt = np.zeros(n_para)
+    vt = np.zeros(n_para)
+
+    traj = []
+    if verbose:
+        print((" "*5).join(["step".ljust(10), "loss".ljust(10), "grad_norm".ljust(10)]))
+
+    grads_norm = 0.
+    for j in range(maxiter):
+        traj.append(func(x, *args))
+        if verbose:
+            print((" "*5).join([("%d" %j).ljust(10), ("%.5f" %traj[j]).ljust(10), ("%.5f" %grads_norm).ljust(10)]))
+        if j > 0 and abs(traj[j]-traj[j-1]) <= tol:
+            return x, traj[-1], traj
+
+        grads = np.array(gradsf(x, *args))
+        grads_norm = np.linalg.norm(grads)
+        mt = beta1 * mt + (1.-beta1) * grads
+        vt = beta1 * vt + (1.-beta2) * grads**2
+        mtt = mt / (1-beta1**(j+2))
+        vtt = vt / (1-beta2**(j+2))
+        x = x - eta * mtt / (np.sqrt(vtt) + epsilon)
+        if call_back:
+            call_back(x, *args)
+    traj.append(func(x, *args))
+    return x, traj[-1], traj
+
+
+def spsa_grad(func, x, k, args=(), spsa_iter=10, c=0.1, gamma=0.101):
+
+    dim = len(x)
+    ck = c/(k)**gamma
+    gx = 0.0
+    for i in range(spsa_iter):
+        Delta = 2*np.round(np.random.rand(dim))-1
+        x1 = x + ck*Delta
+        x2 = x - ck*Delta
+        y1 = func(x1, *args) 
+        y2 = func(x2, *args)
+        gx += (y1 - y2) / (2*ck*Delta)
+    gx = gx / spsa_iter
+    return gx
+
+def spsa(func, x, args=(), call_back=None, spsa_iter=10, max_iter=1000, a=0.1, c=0.1, A=100, alpha=0.602, gamma=0.101, tol=1e-8, verbose=False):
+    '''SPSA minimize
+       c: at a level of standard deviation of func
+       A: <=10% of max_iter '''
+    traj = [func(x, *args)]
+    if verbose:
+            print((" "*5).join(["step".ljust(10), "loss".ljust(10), "grad_norm".ljust(10)]))
+
+    grads_norm = 0.
+    for k in range(max_iter):
+        if verbose:
+            print((" "*5).join([("%d" %k).ljust(10), ("%.5f" %traj[k]).ljust(10), ("%.5f" %grads_norm).ljust(10)]))
+        if k > 0 and abs(traj[k]-traj[k-1]) <= tol:
+            return x, traj[-1], traj
+
+        ak = a/(k+1+A)**alpha
+        grads = spsa_grad(func, x, k+1, args, spsa_iter, c=c, gamma=gamma)
+        grads_norm = np.linalg.norm(grads)
+        x = x - ak * grads
+        if call_back:
+            call_back(x, *args)
+        traj.append(func(x, *args))
+
+    return x, traj[-1], traj
+

quafu/circuits/quantum_circuit.py

@@ -22,14 +22,15 @@
 from quafu.elements.instruction import Instruction
 from quafu.elements import Measure, Reset
 from quafu.elements.parameters import ParameterType, ParameterExpression, Parameter
-from quafu.elements.pulses import QuantumPulse
-from quafu.elements.oracle import OracleGate, ControlledOracle
+from ..elements.quantum_gate import CircuitWrapper, ControlledCircuitWrapper
 from ..elements import (
     Barrier,
     ControlledGate,
     Delay,
     QuantumGate,
     XYResonance,
+    KrausChannel, OracleGate, ControlledOracle, QuantumPulse,
+    CircuitWrapper, ControlledCircuitWrapper
 )
 from ..exceptions import CircuitError
 from .classical_register import ClassicalRegister
@@ -41,14 +42,15 @@ class QuantumCircuit:
     Representation of quantum circuit.
     """
 
-    def __init__(self, qnum: int, cnum: Optional[int] = None, *args, **kwargs):
+    def __init__(self, qnum: int, cnum: Optional[int] = None, name="", *args, **kwargs):
         """
         Initialize a QuantumCircuit object
 
         Args:
             qnum (int): Total qubit number used
             cnum (int): Classical bit number, equals to qubit number in default
         """
+        self.name = name
         self.qregs = [QuantumRegister(qnum)] if qnum > 0 else []
         cnum = self.num if cnum is None else cnum
         self.cregs = [ClassicalRegister(cnum)] if cnum > 0 else []
@@ -118,7 +120,7 @@ def __lshift__(self, operation: Instruction):
         if max_pos >= self.num:
             raise CircuitError("Operation act on qubit that not allocated")
         self.add_ins(operation)
-        if isinstance(operation, OracleGate):
+        if isinstance(operation, CircuitWrapper):
             self._has_wrap = True
         return self
 
@@ -149,7 +151,7 @@ def add_ins(self, ins: Instruction):
         """
         Add instruction to circuit, with NO checking yet.
         """
-        if isinstance(ins, (QuantumGate, Delay, Barrier, XYResonance)):
+        if isinstance(ins, (QuantumGate, Delay, Barrier, XYResonance, KrausChannel)):
             # TODO: Delay, Barrier added by add_gate for backward compatibility.
             #       Figure out better handling in the future.
             self.add_gate(ins)
@@ -164,8 +166,11 @@ def add_noise(self, channel:str, channel_args, qubits=[], gates=[]):
                 raise  ValueError("Invalid gate name")
 
         newinstructions = []
+        newgates = []
         for op in self.instructions:
             newinstructions.append(op)
+            if isinstance(op, (QuantumGate, Delay, Barrier, XYResonance)):
+                newgates.append(op)
             if isinstance(op, QuantumGate):
                 add_q = False
                 add_g = False
@@ -185,10 +190,22 @@ def add_noise(self, channel:str, channel_args, qubits=[], gates=[]):
                 if add_q and add_g:
                     for q in op.pos:
                         if q in qubits:
-                            newinstructions.append(Instruction.gate_classes[channel](q, *channel_args))
+                            newinstructions.append(Instruction.ins_classes[channel](q, *channel_args))
+                            newgates.append(Instruction.ins_classes[channel](q, *channel_args))
         self.instructions = newinstructions
+        self._gates = newgates
         return self
 
+    @property
+    def noised(self):
+        if self._has_wrap:
+            self.unwrap()
+
+        for op in self.instructions:
+            if isinstance(op, KrausChannel):
+                return True
+        return False
+
     def get_parameter_grads(self):
         if self._has_wrap:
             print("warning: The circuit has wrapped gates, it will unwarp automaticllay")
@@ -555,7 +572,7 @@ def add_controls(self, ctrlnum, ctrls:List[int]=[], targs:List[int]=[], inplace=
         return qc
 
     def join(self, 
-             qc : ["QuantumCircuit", OracleGate, ControlledOracle], 
+             qc : ["QuantumCircuit",CircuitWrapper, ControlledCircuitWrapper], 
              qbits:List[int]=[], 
              inplace=True)->"QuantumCircuit":
 
@@ -624,17 +641,29 @@ def dagger(self, inplace=False):
                     nq << op
             nq.measures = self.measures
             return nq
+
+    def wrap(self, qbits=[]):
+        #TODO:use OracleGate
+        name = self.name if self.name else "Oracle"
+        return CircuitWrapper(name, self, qbits)
 
     def unwrap(self):
         instructions = []
+        gates = []
         for op in self.instructions:
-            if isinstance(op, OracleGate) or isinstance(op, ControlledOracle):
+            if isinstance(op, CircuitWrapper):
                 circ = op.circuit.unwrap()
                 for op_ in circ.instructions:
                     instructions.append(op_)
+                    if isinstance(op_, (QuantumGate, Delay, Barrier, XYResonance, KrausChannel)):
+                        gates.append(op_)
             else:
                 instructions.append(op)
+                if isinstance(op, (QuantumGate, Delay, Barrier, XYResonance, KrausChannel)):
+                    gates.append(op)
+
         self.instructions = instructions
+        self._gates = gates
         self._has_wrap = False
         return self
 

quafu/elements/__init__.py

@@ -1,6 +1,9 @@
 from .classical_element import Cif
 from .instruction import Barrier, Instruction, Measure, Reset
 from .pulses import Delay, QuantumPulse, XYResonance
-from .quantum_gate import ControlledGate, QuantumGate
+from .quantum_gate import ControlledGate, QuantumGate, CircuitWrapper, ControlledCircuitWrapper
 from .unitary import UnitaryDecomposer
 from .utils import extract_float, reorder_matrix
+from .oracle import OracleGate, ControlledOracle
+from .parameters import Parameter, ParameterExpression, ParameterType
+from .noise import KrausChannel, UnitaryChannel