Merge pull request #137 from lss0208/master

feat:add iqcs features

CMakeLists.txt

@@ -6,8 +6,9 @@ project(qfvm LANGUAGES CXX C)
 
 if(NOT CMAKE_BUILD_TYPE)
     message(STATUS "Build type not set - defaulting to Release")
-    set(CMAKE_BUILD_TYPE Release)
+    set(CMAKE_BUILD_TYPE Debug)
 endif()
+set(CMAKE_BUILD_TYPE Debug)
 
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CUDA_ARCHITECTURES 70;75;80;90)
@@ -40,6 +41,7 @@ else()
 	endif()
 endif()
 
+list ( APPEND PRJ_COMPILE_OPTIONS -g)
 #install eigen
 set(EIGEN3_ROOT  ${CMAKE_SOURCE_DIR}/thirdparty/Eigen3)
 set(EIGEN3_INCLUDE_DIR ${EIGEN3_ROOT}/include/eigen3)

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/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/algorithms/templates/amplitude.py

@@ -191,7 +191,7 @@ def _apply_uniform_rotation_dagger(gate, alpha, control_wires, target_wire):
 
     if gray_code_rank == 0:
         if np.all(theta[..., 0] != 0.0):
-            gate_list.append(gate(pos=target_wire, paras=theta[0]))
+            gate_list.append(gate(target_wire, theta[0]))
         return gate_list
 
     code = gray_code(gray_code_rank)
@@ -204,7 +204,7 @@ def _apply_uniform_rotation_dagger(gate, alpha, control_wires, target_wire):
 
     for i, control_index in enumerate(control_indices):
         if np.all(theta[..., i] != 0.0):
-            gate_list.append(gate(pos=target_wire, paras=theta[i]))
+            gate_list.append(gate(target_wire, theta[i]))
         gate_list.append(qeg.CXGate(control_wires[control_index], target_wire))
     return gate_list
 

quafu/algorithms/templates/angle.py

@@ -43,7 +43,7 @@ def __init__(self, features, num_qubits: int, rotation="X"):
     def _build(self):
         gate_list = []
         for i in range(self.num_qubits):
-            gate = self.op(pos=i, paras=self.features[i])
+            gate = self.op(i, self.features[i])
             gate_list.append(gate)
         return gate_list
 

quafu/algorithms/templates/basic_entangle.py

@@ -59,7 +59,7 @@ def _build(self):
         gate_list = []
         for layer in range(repeat):
             for i in range(self.num_qubits):
-                gate = self.op(pos=i, paras=self.weights[layer][i])
+                gate = self.op(i, self.weights[layer][i])
                 gate_list.append(gate)
 
             # if num_qubits equals two, it just need to apply CNOT one time