RGD optimization method for quantum state tomography

.gitignore

@@ -3,6 +3,7 @@
 __pycache__/
 *.py[cod]
 *$py.class
+src/qibo/tomography_RGD/testData/
 
 # C extensions
 *.so

src/qibo/tomography_RGD/BasicTools.py

@@ -0,0 +1,76 @@
+#
+#   some tools
+#   1. generate all labels
+#   2. plot the results
+#
+
+import matplotlib.pyplot as plt
+
+
+# --------------------------------- #
+#   to generate all symbols         #
+# --------------------------------- #
+def Generate_All_labels(Nk, symbols=["I", "X", "Y", "Z"]):
+    """generate all possible labels
+
+    Args:
+        Nk (int): number of qubits
+        symbols (list, optional): the possible choice of each qubit site. Defaults to ['I', 'X', 'Y', 'Z'].
+
+    Returns:
+        list: list of all possible labels
+    """
+
+    symList = symbols
+
+    for i in range(1, Nk):
+        print("   the {}-th qubit".format(i))
+
+        sym_Generated = []
+        for symNow in symList:
+            sym_Generated = sym_Generated + ["".join([symNow, s]) for s in symbols]
+            # print(sym_Generated)
+        symList = sym_Generated
+    print("  totol number of labels {}".format(len(symList)))
+
+    return symList
+
+
+def Plt_Err_Time(worker):
+    """plot the Error w.r.t. optimization run time
+
+    Args:
+        worker (class): the optimization class instance
+    """
+
+    Target_Err_Xk = worker.Target_Err_Xk
+    step_Time = worker.step_Time
+
+    RunT = []
+    Ttot = 0
+    for ti in range(-1, len(step_Time) - 1):
+        Each_Step = step_Time[ti]
+        Ttot += Each_Step
+        RunT.append(Ttot)
+
+    mk_list = ["+", "^", "o", "x", ">", "<", 2, 3]
+    ln_list = ["-", "-.", "--", "--", ":", "-"]
+
+    fig, axNow = plt.subplots(1, 1, figsize=(8, 6))
+
+    info = "{} qubits with sampling {} labels".format(worker.n, worker.num_labels)
+    axNow.plot(
+        RunT,
+        Target_Err_Xk,
+        marker=mk_list[0],
+        linestyle=ln_list[0],
+        label="{}".format(info),
+    )
+
+    axNow.set_xlabel("  Run time (sec)", fontsize=14)
+    axNow.set_ylabel(r"$\left\Vert X_k -\rho \right\Vert_F$", fontsize=14)
+
+    axNow.set_title("Error w.r.t. run time", y=1.0, fontsize=14)
+
+    plt.legend(loc="upper left")
+    plt.show()

src/qibo/tomography_RGD/Run_Tomo.py

@@ -0,0 +1,196 @@
+#
+#   prepare the state, do the measurements
+#   and directly do the tomography
+#
+#
+
+import measurements
+import methodsMiFGD_core
+import methodsRGD_core
+import numpy as np
+import projectors
+import qutip as qu
+from BasicTools import Generate_All_labels, Plt_Err_Time
+
+# from states import GHZState, HadamardState, RandomState
+from qibo_states import GHZState, HadamardState, RandomState
+
+from qibo import quantum_info
+
+if __name__ == "__main__":
+
+    ############################################################
+    ### Example of creating and running an experiment
+    ############################################################
+
+    # n = 3;    labels = projectors.generate_random_label_list(50, n)
+    n = 4
+    labels = projectors.generate_random_label_list(120, n)
+
+    # labels = ['YXY', 'IXX', 'ZYI', 'XXX', 'YZZ']
+    # labels = ['YZYX', 'ZZIX', 'XXIZ', 'XZIY', 'YXYI', 'ZYYX', 'YXXX', 'IIYY', 'ZIXZ', 'IXXI', 'YZXI', 'ZZYI', 'YZXY', 'XYZI', 'XZXI', 'XZYX', 'YIXI', 'IZYY', 'ZIZX', 'YXXY']
+    # labels = ['IIIX', 'IYIY', 'YYXI', 'ZZYY', 'ZYIX', 'XIII', 'XXZI', 'YXZI', 'IZXX', 'YYIZ', 'XXIY', 'XXZY', 'ZZIY', 'YIYX', 'YYZZ', 'YZXZ', 'YZYZ', 'ZXYY', 'IXIZ', 'XZII']
+    # labels = Generate_All_labels(n)
+
+    num_labels = len(labels)
+
+    circuit_Choice = 1
+    if circuit_Choice == 1:  #  generate from circuit
+        Nr = 1
+
+        # state   = GHZState(n)
+        # state   = HadamardState(n)
+        state = RandomState(n)
+
+        target_density_matrix = state.get_state_matrix()
+        target_state = state.get_state_vector()
+        # print(state.get_state_vector())
+
+        #
+        # DO the shot measurement
+        #
+
+        state.create_circuit()
+        data_dict_list = state.execute_measurement_circuits(labels)
+        # print(data_dict_list)
+
+        #
+        # shot measurement results -->  coefficient for each Pauli operator
+        #
+
+        data_dict = {}
+        for ii in range(num_labels):
+            label = data_dict_list[ii]["label"]
+            data_dict[label] = data_dict_list[ii]["count_dict"]
+
+        measurement_list = measurements.MeasurementStore.calc_measurement_list(
+            labels, data_dict
+        )
+
+    elif circuit_Choice == 2:  # directly generate density matrix via qutip
+        Nr = 1
+        rho = qu.rand_dm_ginibre(2**n, dims=[[2] * n, [2] * n], rank=Nr)
+        target_density_matrix = rho.full()
+
+    elif circuit_Choice == 3:  # directly generate density matrix via qibo
+        Nr = 3
+        target_density_matrix = random_density_matrix(2**n, Nr)
+
+    # ---------------------------------------------------------------- #
+    #  construct Pauli matrix Projectors                               #
+    # ---------------------------------------------------------------- #
+
+    projector_store_path = "./testData/qiskit"
+    # projector_store_path = './testData/qibo'
+    projector_store = projectors.ProjectorStore(labels)
+
+    ## store method (1)
+    # projector_store.populate(projector_store_path)
+    # projector_dict = projectors.ProjectorStore.load(projector_store_path, labels)
+
+    ## store method (2)
+    num_cpus, saveP_bulk, Partition_Pj = projector_store.mpPool_map(
+        projector_store_path
+    )
+    projector_dict = projectors.ProjectorStore.load_PoolMap(
+        projector_store_path, labels
+    )
+
+    # ---------------------------------------------------------------- #
+    #  calculate exact coefficient for each Pauli operator             #
+    # ---------------------------------------------------------------- #
+
+    projector_list = [projector_dict[label] for label in labels]
+    yProj_Exact = methodsRGD_core.Amea(
+        projector_list, target_density_matrix, num_labels, 1
+    )  # argv[-1] = coef
+
+    if circuit_Choice == 1:  #  generated from circuit
+        #
+        #   comparison: manual check the result of shot measurements
+        #
+        ml = np.array(measurement_list, dtype=float)
+        ind = np.where(np.abs(yProj_Exact) > 0.5)
+
+        print(yProj_Exact[ind])
+        print(ml[ind])
+    elif circuit_Choice == 2 or circuit_Choice == 3:  # directly generate density matrix
+        measurement_list = yProj_Exact
+        target_state = None
+
+    #
+    #   system parameters
+    #
+    params_dict = {
+        "Nr": Nr,
+        "target_DM": target_density_matrix,
+        "labels": labels,
+        "measurement_list": measurement_list,
+        "projector_list": projector_list,
+        "num_iterations": 150,
+        "convergence_check_period": 1,
+    }
+    # params_dict['target_state'] = target_state
+
+    # ----------------------------------------------------------------- #
+    #   do the tomography optimization                                  #
+    # ----------------------------------------------------------------- #
+
+    #
+    #   MiFGD numerical parameters
+    #
+
+    # Call_MiFGD = 1      #  = 1: call the MiFGD optimization to calculate
+    # muList = [4.5e-5]
+
+    InitX_MiFGD = 1  # 0: random start,  1: MiFGD specified init
+    mu = 4.5e-5
+    eta = 0.01
+    Option = 2
+
+    # Num_mu = 1      #  number of mu for running MiFGD
+    # pm_MiFGD = [Call_MiFGD, InitX_MiFGD, muList, Num_mu]
+    # Call_MiFGD, InitX_MiFGD, muList, Num_mu = pm_MiFGD
+
+    params_MiFGD = {"mu": mu, "eta": eta, "Option": Option}
+    # params_dict = {**params_dict, **params_MiFGD}
+
+    # Rpm_MiFGD = [InitX_MiFGD, mu, eta, Option]
+    # Frec_MiFGD, wc, RunTime = Run_MiFGD(params_dict, Rpm_MiFGD)
+
+    worker2 = methodsMiFGD_core.BasicWorker(params_dict, params_MiFGD)
+    worker2.compute(InitX_MiFGD)
+
+    #
+    #   RGD numerical parameters
+    #
+
+    print("\n +++++++++++++++++     do the RGD tomography     +++++++++++++++++\n")
+
+    # Call_RGD = 1        #  = 1: call the RGD   optimization to calculate
+
+    Md_tr = 0  #   Method if including trace = 1
+    Md_alp = 0  #   method for scaling alpha
+    Md_sig = 0  #   method for scaling singular value
+    Ch_svd = (
+        -1
+    )  #   choice for initial SVD  (0: LA.svd, 1: svds;  2: power_Largest_EigV, -1: rSVD)
+    InitX_RGD = 1
+    #   method of choosing initial X0
+
+    # Rpm    = [InitX_RGD, Md_tr, Md_alp, Md_sig, Ch_svd]
+    # pm_RGD = [Call_RGD, Rpm]
+    # Call_RGD, Rpm   = pm_RGD
+
+    # InitX_RGD, Md_tr, Md_alp, Md_sig, Ch_svd = Rpm
+    # InitX_RGD = 1
+    # Rpm = InitX_RGD, Md_tr, Md_alp, Md_sig, Ch_svd
+
+    # exec(open('RGD_optRun.py').read())
+    # Frec_RGD, wc, RunTime = Run_RGD(params_dict, Rpm)
+
+    worker = methodsRGD_core.BasicWorkerRGD(params_dict)
+    # worker.computeRGD(InitX_RGD, Ch_svd, Md_tr, Md_alp, Md_sig)
+    worker.computeRGD(InitX_RGD, Ch_svd)
+
+    Plt_Err_Time(worker)