diff --git a/.gitignore b/.gitignore index 1027fc4ff3..ec085d5b59 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,7 @@ __pycache__/ *.py[cod] *$py.class +src/qibo/tomography_RGD/testData/ # C extensions *.so diff --git a/src/qibo/tomography_RGD/BasicTools.py b/src/qibo/tomography_RGD/BasicTools.py new file mode 100755 index 0000000000..e76a5b2bdd --- /dev/null +++ b/src/qibo/tomography_RGD/BasicTools.py @@ -0,0 +1,57 @@ +# +# some tools +# 1. generate all labels +# 2. plot the results +# + +# --------------------------------- # +# to generate all symbols # +# --------------------------------- # +from itertools import product + +import matplotlib.pyplot as plt + + +def Generate_All_labels(nqubits, symbols=["I", "X", "Y", "Z"]): + all_labels = ["".join(elem) for elem in product(symbols, repeat=nqubits)] + return all_labels + + +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() diff --git a/src/qibo/tomography_RGD/RGD.py b/src/qibo/tomography_RGD/RGD.py new file mode 100644 index 0000000000..2bbeed1fc0 --- /dev/null +++ b/src/qibo/tomography_RGD/RGD.py @@ -0,0 +1,845 @@ +import time + +import numpy as np +from numpy import linalg as LA +from scipy.sparse.linalg import svds + + +class BasicParameterInfo: + def __init__(self, params_dict): + + # + # read in params_dict + # + Nr = params_dict.get("Nr", 1) + trace = params_dict.get("trace", 1.0) + target_state = params_dict.get("target_state", None) + target_DM = params_dict.get("target_DM") + + label_list = params_dict.get("labels") + measurement_list = params_dict.get("measurement_list") + # projector_list = params_dict.get("projector_list") + # symProj_list = params_dict.get("symProj_list") + + num_iterations = params_dict.get("num_iterations") + convergence_check_period = params_dict.get("convergence_check_period", 10) + + relative_error_tolerance = params_dict.get("relative_error_tolerance", 0.001) + seed = params_dict.get("seed", 0) + + # + # basic system information + # + n = len(label_list[0]) + d = 2**n + + self.n = n # number of qubits + self.num_elements = d + self.Nr = Nr # the rank of the target_density_matrix + + self.target_DM = target_DM # target state --> density matrix + + self.num_labels = len(label_list) + self.measurement_list = measurement_list + # self.projector_list = projector_list + # self.symProj_list = symProj_list + + # + # numerical book keeping + # + + self.process_idx = 0 + + self.num_iterations = num_iterations + + self.convergence_check_period = ( + convergence_check_period # how often to check convergence + ) + self.relative_error_tolerance = ( + relative_error_tolerance # user-decided relative error + ) + self.seed = seed + + self.Err_relative_st = [] + + self.Target_Err_Xk = [] # Frobenius norm difference from matrix Xk + self.Err_relative_Xk = [] + + self.iteration = 0 + self.converged = False + self.convergence_iteration = 0 + + +class Pauli_Op: + def __init__(self, symProj_list): + + self.symProj_list = symProj_list + + def Amea(self, Xmatrix): + + yM = [np.trace(Proj @ Xmatrix).real for Proj in self.symProj_list] + return np.array(yM) + + def A_dagger(self, yMea): + r"""implements A_dagger(yMea) = \sum_j yMea_j * (Pauli operator)_j + + self.symProj_list (list): list of Pauli operators in the form of SymbolicHamiltonian + + Args: + yMea (ndarray): the input vector as the coefficient array for the Pauli matrices + + Returns: + Adag (SymbolicHamiltonian): calculated matrix A_dagger(yMea) + """ + # from functools import reduce + # xx = [yP*Proj for yP, Proj in zip(yMea, self.symProj_list)] + # yy = reduce(lambda x, y: x+y, xx) + Adag = sum(map(lambda a, b: a * b, yMea, self.symProj_list)) + + return Adag + + def A_dagger_ym_vec(self, yMea, vec): + """to obtain A^+(yMea) @ vec + + self.symProj_list (list): list of Pauli operators in the form of SymbolicHamiltonian + + Args: + yMea (ndarray): the input vector for the A^+ operator + symProj_list (list): list of Pauli operators in the form of SymbolicHamiltonian + vec (ndarray): a vector which is supposed to be the right singular vector + + Returns: + Ad_vec (ndarray): the A^+(y) @ vec + + """ + + # Adag = A_dagger(yMea, self.symProj_list) + # Ad_vec = Adag.matrix @ vec + + # from functools import reduce + # qq = [yP*Proj for yP, Proj in zip(yMea, self.symProj_list)] + # yy = [yP.matrix @ vec for yP in qq] + # Ad_vec = reduce(lambda x, y: x+y, yy) + + xx = [Proj.matrix @ vec for Proj in self.symProj_list] + Ad_vec = sum(map(lambda a, b: a * b, yMea, xx)) + + return Ad_vec + + +class BasicWorkerRGD(BasicParameterInfo, Pauli_Op): + """ + Basic worker class. Implements RGD, which performs + accelerated gradient descent on the factored space U + (where UU^* = rho = density matrix) + """ + + def __init__( + self, + target_DM, + labels, + measurement_list, + symProj_list, + Nr=1, + convergence_check_period=1, + num_iterations=150, + seed=0, + relative_error_tolerance=0.001, + ): + # + # read in params_dict + # + + # BasicParameterInfo.__init__(self, params_dict) + + # + # basic system information + # + n = len(labels[0]) + d = 2**n + + self.n = n # number of qubits + self.num_elements = d + + self.Nr = Nr # the rank of the target_density_matrix + + self.target_DM = target_DM # target state --> density matrix + + self.num_labels = len(labels) + self.measurement_list = measurement_list + + Pauli_Op.__init__(self, symProj_list) + + # + # numerical book keeping + # + + self.num_iterations = num_iterations + + self.convergence_check_period = ( + convergence_check_period # how often to check convergence + ) + self.relative_error_tolerance = ( + relative_error_tolerance # user-decided relative error + ) + + self.seed = seed + + self.Err_relative_st = [] + + self.Target_Err_Xk = [] # Frobenius norm difference from matrix Xk + self.Err_relative_Xk = [] + + self.iteration = 0 + self.converged = False + self.convergence_iteration = 0 + + def initialize_RndSt(self): + """ + Random initialization of the state density matrix. + """ + + print("\n ****** initial choice same as the MiFGD ***** \n") + print("self.seed = {}".format(self.seed)) + self.InitChoice = "random initial (same as MiFGD)" + + seed = self.seed + + stateIni = [] + for xx in range(self.Nr): + # real_state = np.random.RandomState(seed).randn(self.num_elements) + # imag_state = np.random.RandomState(seed).randn(self.num_elements) + # seed += 1 + + real_state = np.random.RandomState().randn(self.num_elements) + imag_state = np.random.RandomState().randn(self.num_elements) + + state = real_state + 1.0j * imag_state + + state_norm = np.linalg.norm(state) + state = state / state_norm + + stateIni.append(state) + + stateIni = np.array(stateIni).T + + self.state = stateIni + self.uk = stateIni + self.vk = stateIni + self.ukh = stateIni.T.conj() + self.vkh = np.transpose(stateIni.conj()) + + self.s0_choice = 0 + if self.s0_choice == 0: + # ss = np.random.random(self.Nr) + # ss = ss / np.sum(ss) + ss = np.ones(self.Nr) + + elif self.s0_choice == 1: + ss = np.ones(self.Nr) / self.Nr + + self.sDiag = np.diag(ss) + + X0 = stateIni @ self.sDiag @ self.vkh + + self.u0 = stateIni + self.v0 = stateIni + self.s0 = np.diag(ss) + self.X0 = X0 + + self.Xk = X0 + + self.check_Hermitian_x_y(self.uk, self.vk) + + def initialize_yAd(self): + """initialize the initial density matrix for the algorithm""" + + self.InitChoice = "paper" + + # ------------------------------------------------------------ # + + if self.Ch_svd >= 0: # need to compute X0 + + print(" start running self.A_dagger") + tt1 = time.time() + # y = self.yIni ## this y = coef * measurement_list -> order follows self.label_list + + Xk = self.A_dagger(self.yIni).matrix * self.coef + + tt2 = time.time() + + print(" self.A_dagger done --> time = {}".format(tt2 - tt1)) + print(" Initial SVD --> choice = {}".format(self.Ch_svd)) + + uk, sDiag, vkh = Hr_Hermitian( + Xk, self.Nr, self.Ch_svd + ) # Ch_svd =0: full SVD + + elif self.Ch_svd < 0: # don't need to compute X0 + + print( + " ********* using randomized-SVD to construct X0 = uk @ sDiag @ vkh **************************" + ) + + tt2a = time.time() + uk, sDiag, vkh = self.rSVD(k=self.Nr, s=2, p=15) + tt2b = time.time() + print(" self.rSVD --> time = {}".format(tt2b - tt2a)) + + # Compare_2_SVD(u1, sDiag1, v1h, uk, sDiag, vkh) + + self.check_Hermitian_x_y(uk, vkh.T.conj()) + + if self.EigV_positive == 1: # U == V + vk = uk + vkh = uk.T.conj() + + else: # U != V + vk = vkh.T.conj() + + ukh = uk.T.conj() + Xk = uk @ sDiag @ vkh # only keep rank r singular vectors + + self.u0 = uk + self.v0 = vk + self.s0 = sDiag + self.X0 = Xk + + self.Xk = Xk + self.uk = uk + self.vk = vk + + self.ukh = ukh + self.vkh = vkh + self.sDiag = sDiag + + def rSVD(self, k, s=3, p=12, matrixM_Hermitian=1): + """randomized SVD + + ref: https://arxiv.org/pdf/1810.06860.pdf --> A ~ QQ* A (B = Q* A ) + with modification: if Remain_Hermitian --> A ~ QQ* A QQ* (B = Q* A Q) + + Args: + k (int): = Nr # rank parameter + s (int, optional): specify oversampling. Defaults to 3. + p (int, optional): specify power parameter. Defaults to 12. + matrixM_Hermitian (int, optional): specify the matrix is Hermitian or not. Defaults to 1. + + Returns: + ndarray: (u) left singular vectors + ndarray: np.diag(s) the diagonal matrix with diagonal elements being s + ndarray: (vh) complex conjugate of the right singular vectors + """ + + qDs = min(k + s, self.num_elements) # k + s small dimension for Q + + Omega = np.random.RandomState().randn(self.num_elements, qDs) + + Mw = self.A_dagger_ym_vec(self.yIni, Omega) + + Q, R = LA.qr(Mw) + + for ii in range(p): + tt1 = time.time() + + ATQ = self.A_dagger_ym_vec(self.yIni, Q) + + G, R = LA.qr(ATQ) + tt2 = time.time() + + AG = self.A_dagger_ym_vec(self.yIni, G) + + Q, R = LA.qr(AG) + + tt3 = time.time() + print( + " *** {}-th (A A*) done: Time --> ATQ: {}, AG: {} ***".format( + ii, tt2 - tt1, tt3 - tt2 + ) + ) + + if matrixM_Hermitian == 0: + # B = Q.T.conj() @ M + # B = self.A_dagger_vec(Q).T.conj() + B = self.A_dagger_ym_vec(self.yIni, Q).T.conj() * self.coef + + uB, s, vh = LA.svd(B) + u = Q @ uB[:, :k] # k = Nr + s = s[:k] # k = Nr + vh = vh[:k, :] # k = Nr + + elif matrixM_Hermitian == 1: + # B = Q.T.conj() @ M @ Q + # B = Q.T.conj() @ self.A_dagger_vec(Q) + B = Q.T.conj() @ self.A_dagger_ym_vec(self.yIni, Q) * self.coef + + uB, s, vh = LA.svd(B, hermitian=True) + u = Q @ uB[:, :k] # k = Nr + s = s[:k] # k = Nr + + vh = vh[:k, :] @ Q.T.conj() # k = Nr + + return u, np.diag(s), vh + + def check_Hermitian_x_y(self, x, y, Hermitian_criteria=1e-13): + """to check wheter x and y are the same, + where x and y are supposed to be the left and right singular vectors + if x = y, then the original SVDed matrix is Hermitian + + Args: + x (ndarray): the 1st vector to compare, supposed to be the left singular vector + y (ndarray): the 2nd vector to compare, supposed to be the right singular vector + Hermitian_criteria (float, optional): the numerical error minimum criteria. Defaults to 1e-13. + """ + + # if LA.norm(x - y) < Hermitian_criteria: + if np.allclose(x, y): + + self.Remain_Hermitian = 1 # True for remaining Hermitian + self.EigV_positive = 1 + + # elif LA.norm(Col_flip_sign(x) - Col_flip_sign(y)) < Hermitian_criteria: + elif np.allclose(Col_flip_sign(x), Col_flip_sign(y)): + self.Remain_Hermitian = 1 + self.EigV_positive = 0 + + else: + self.EigV_positive = None + self.Remain_Hermitian = 0 # False for remaining Hermitian + self.st_Non_Hermitian.append( + self.iteration + ) # start num for non-Hermitian + + self.EigV_pm_List.append(self.EigV_positive) + self.Hermitian_List.append(self.Remain_Hermitian) + + # print( + # " Hermitian ? Remain_Hermiain = {}, EigV_postive = {}".format( + # self.Remain_Hermitian, self.EigV_positive + # ) + # ) + + def Get_measured_y(self): + """to get the coefficients of all sampled Pauli operators for the A^+ operator""" + + coef = np.sqrt(self.num_elements / self.num_labels) + self.yIni = np.array(self.measurement_list) * coef + self.coef = coef + + def computeRGD(self, InitX, Ch_svd=0): + """Basic workflow of gradient descent iteration. + + 1. initializes state dnesity matrix. + 2. makes a step (defined differently for each "Worker" class below) until + convergence criterion is satisfied. + + Args: + Ch_svd (int, optional): _description_. Defaults to 0. + """ + + self.Ch_svd = Ch_svd # choice for initial SVD + + self.EigV_pm_List = ( + [] + ) # Eig-value is positive or not (=1: positive, = 0, non-positive, = None, no EigV) + self.Hermitian_List = [] + self.st_Non_Hermitian = [] + self.step_Time = {} + + tt1 = time.time() + + self.iteration = -1 + + self.Get_measured_y() + print(" self.coef = {}".format(self.coef)) + print(" InitX = {}".format(InitX)) + + if InitX == 0: + self.initialize_RndSt() # initial by random + elif InitX == 1: + self.initialize_yAd() # initial by A^\dagger(y) + + tt2 = time.time() + self.step_Time[-1] = tt2 - tt1 + + print( + " -------- X0 --> (-1)-th iteration [ SVD choice = {} ] \n".format( + Ch_svd + ) + ) + print(" X0 step (RGD) --> time = {}\n".format(tt2 - tt1)) + + X0_target_Err = LA.norm(self.Xk - self.target_DM, "fro") + self.Target_Err_Xk.append(X0_target_Err) # for X0 + self.InitErr = X0_target_Err + + print( + " X0 --> Tr(X0) = {}\n".format( + np.trace(self.X0) + ) + ) + print(" X0 --> Target Err = {}\n".format(X0_target_Err)) + print(" RGD max num_iterations = {}\n".format(self.num_iterations)) + + self.uGv_list = [] + self.Rec_Alpha = [] + self.Alpha_sc = [] + self.sDiag_list = [] + + for self.iteration in range(self.num_iterations): + if not self.converged: + self.stepRGD() + else: + break + if self.convergence_iteration == 0: + self.convergence_iteration = self.iteration + + def stepRGD(self): + """each iteration step of the RGD algorithm""" + + self.Xk_old = self.Xk + self.uk_old = self.uk + + print(" -------- {}-th iteration \n".format(self.iteration)) + tt1 = time.time() + + # ----------------------------------------- # + # calc of Axk = mapA(Xk) # + # ----------------------------------------- # + + Axk = self.Amea(self.Xk) * self.coef ## 1st quickest + self.zm = ( + self.yIni - Axk + ) ## [yIni = coef*measurement_list] order follows self.label_list + + # ----------------------------------------------------- # + # calculate (n1, n2) by direct calculation of Gk # + # ----------------------------------------------------- # + + print( + " {}-th step -> self.Remain_Hermitian = {}, self.EigV_positive = {}".format( + self.iteration, self.Remain_Hermitian, self.EigV_positive + ) + ) + + uGv = self.calc_PtG_2_uG_Gv() + + Alpha = self.calc_n1_n2() + + self.Hr_tangentWk(Alpha, uGv) + + if self.iteration % 100 == 0: + print(" **** {}-th iteratation\n".format(self.iteration)) + + tt7 = time.time() + self.convergence_check() + tt8 = time.time() + + print(" convergence check --> time = {}\n".format(tt8 - tt7)) + print(" stepRGD --> time = {}\n".format(tt8 - tt1)) + + self.step_Time[self.iteration] = tt8 - tt1 + + def convergence_check(self): + """ + Check whether convergence criterion is satisfied by comparing + the relative error of the current estimate and the target density matrices. + Utilized in "step" function below. + """ + if ( + # self.process_idx == 0 + # and + self.iteration % self.convergence_check_period + == 0 + ): + # compute relative error + + XkErrRatio = LA.norm(self.Xk - self.Xk_old) / LA.norm(self.Xk_old) + + xx = LA.norm(Col_flip_sign(self.uk) - Col_flip_sign(self.uk_old)) + + ukErrRatio = xx / LA.norm(self.uk_old) + + print(" min(uk +/- uk_old) = {}\n".format(xx)) + + # self.Rec_st_Err.append( LA.norm(self.uk - self.uk_old) ) + self.Err_relative_Xk.append(XkErrRatio) + + self.Err_relative_st.append(ukErrRatio) + + if ( + min(ukErrRatio, XkErrRatio) <= self.relative_error_tolerance + ): # using state or Xk to check convergence + self.converged = True + self.convergence_iteration = self.iteration + print( + " ********* XkErrRatio = {} < StpTol ******".format(XkErrRatio) + ) + + if self.target_DM is not None: + Fro_diff = LA.norm(self.Xk - self.target_DM, "fro") + + if np.isnan(Fro_diff): + print( + " ********* {}-th Fro_diff = NaN = {}".format( + self.iteration, Fro_diff + ) + ) + # break + return + + self.Target_Err_Xk.append(Fro_diff) + + print( + " relative_errU = {}, relative_errX = {}".format( + ukErrRatio, XkErrRatio + ) + ) + print(" Target_Error = {}\n".format(Fro_diff)) + + if Fro_diff > 1e2: + raise ValueError(" Target_Err_Xk too big !!! Stop and check \n") + + def calc_PtG_2_uG_Gv(self): + """to project Gk onto the tangnet space and + calculate the related objects + + Returns: + ndarray: (uGv) uk@ Gk @ vk + """ + + if self.EigV_positive == 1: # Hermitian & EigV > 0 --> U = V + + Gu = self.A_dagger_ym_vec(self.zm, self.uk) + + self.Gu = self.coef * Gu + uGv = self.ukh @ self.Gu + + self.PtGk = ( + self.uk @ self.Gu.T.conj() + + self.Gu @ self.ukh + - self.uk @ uGv @ self.ukh + ) + + else: # U != V + + Gu = self.A_dagger_ym_vec(self.zm, self.uk) + Gv = self.A_dagger_ym_vec(self.zm, self.vk) + + print( + " np.allclose( Col_flip_sign(Gu), Col_flip_sign(Gv) ) = {}".format( + np.allclose(Col_flip_sign(Gu), Col_flip_sign(Gv)) + ) + ) + + self.Gv = self.coef * Gv + self.Gu = self.coef * Gu + uGv = self.ukh @ self.Gv + + self.PtGk = ( + self.uk @ self.Gu.T.conj() + + self.Gv @ self.vkh + - self.uk @ uGv @ self.vkh + ) + + return uGv + + def calc_n1_n2(self): + # ----------------------------- # + # calc AptGk = A(PtGk) # + # ----------------------------- # + + AptGk = self.coef * self.Amea(self.PtGk) + + n2 = LA.norm(AptGk) + n1 = LA.norm(self.PtGk) + + # --------------------------------- # + # calculate (n1, n2) # + # --------------------------------- # + + if n1 == 0.0: + print("Projected Gradient PtGk norm = n1 = 0 = {}".format(n1)) + print(" --> should achieve the minimum | or not successful ??") + return + + Alpha = (n1 / n2) ** 2 + + return Alpha + + def Hr_tangentWk(self, Alpha, uGv): + """Hard thresholding Hr(.) of the intermediate matrix Jk + according to the 4th step in the RGD algorithm in the paper, + and update the SVD of the new predicted self.Xk = Hr(Jk) + + Args: + Alpha (float): the step size in the 2nd step of the RGD algorithm + uGv (ndarray): uk @ Gk @ vk + + Returns: + ndarray: (uk) updated left singular vectors of self.Xk + ndarray: (sDiag) updated singular values of self.Xk + ndarray: (vk) updated right singular vectors of self.Xk + """ + + Nr = self.Nr + + # ------------------------------------------------- # + # re-scaling Alpha according to uGv or sDiag # + # ------------------------------------------------- # + max_scaling_time = 1 + for ii_sc in range(max_scaling_time): + # print(" ************ {}-th scaling for the Alpha".format(ii_sc)) + + D2s = Alpha * uGv + self.sDiag + + if self.EigV_positive == 1: # Hermitian & all EigVal > 0 + Y2 = Alpha * (self.Gu - self.uk @ uGv) + + q2, r2 = LA.qr(Y2) + # print("k={}, (q2.shape, r2.shape) = ({}, {})".format(k, q2.shape, r2.shape)) + U2r = np.concatenate((self.uk[:, :Nr], q2[:, :Nr]), axis=1) + # print("U2r shape = {}".format(U2r.shape)) + + D2s = np.concatenate((D2s, r2[:Nr, :Nr].T.conj()), axis=1) + # print(D2s.shape) + + else: # not Hermitian OR Eig-val < 0, i.e. U != V + Y1 = Alpha * (self.Gu - self.vk @ uGv.T.conj()) + Y2 = Alpha * (self.Gv - self.uk @ uGv) + + q1, r1 = LA.qr(Y1) + # print("k={}, (q1.shape, r1.shape) = ({}, {})".format(k, q1.shape, r1.shape)) + q2, r2 = LA.qr(Y2) + # print("k={}, (q2.shape, r2.shape) = ({}, {})".format(k, q2.shape, r2.shape)) + U2r = np.concatenate((self.uk[:, :Nr], q2[:, :Nr]), axis=1) + # print("U2r shape = {}".format(U2r.shape)) + V2r = np.concatenate((self.vk[:, :Nr], q1[:, :Nr]), axis=1) + # print("V2r shape = {}".format(V2r.shape)) + + D2s = np.concatenate((D2s, r1[:Nr, :Nr].T.conj()), axis=1) + # print(D2s.shape) + + D2s = np.concatenate( + (D2s, np.concatenate((r2[:Nr, :Nr], np.zeros((Nr, Nr))), axis=1)), + axis=0, + ) + # print(D2s.shape) + + print( + " Hr_tangentWk --> Remain_Hermitian = {}, EigV_positve = {}".format( + self.Remain_Hermitian, self.EigV_positive + ) + ) + print( + " D2s - D2s.T.conj() = {}".format(LA.norm(D2s - D2s.T.conj())) + ) + print( + " np.allclose(D2s, D2s.T.conj()) = {}\n".format( + np.allclose(D2s, D2s.T.conj()) + ) + ) + + if np.allclose(D2s, D2s.T.conj()) == True: + Mu2r, s2r, Mv2rh = LA.svd(D2s, full_matrices=True, hermitian=True) + + else: # not Hermitian + Mu2r, s2r, Mv2rh = LA.svd(D2s, full_matrices=True) + + Mu2r_Nr = Mu2r[:, :Nr] + Mv2r_Nr = Mv2rh[:Nr, :].T.conj() + + # print(" s2r = {}".format(s2r)) + # print(" max(s2r) = {}".format(max(s2r))) + if max(s2r) < 1: + break + else: + Alpha = Alpha * 0.5 + + # --------------------------------------------------------- # + # END of re-scaling Alpha according to uGv or sDiag # + # --------------------------------------------------------- # + + self.sDiag = np.diag(s2r[:Nr]) + + if self.EigV_positive == 1: # original U = V + + uk = U2r.dot(Mu2r_Nr) + vk = U2r.dot(Mv2r_Nr) # since may Mu2r_Nr != Mv2r_Nr + + else: # original U != V + + uk = U2r.dot(Mu2r_Nr) + vk = V2r.dot(Mv2r_Nr) + + # ------------------------------------- # + # update the Xk, uk, vk # + # ------------------------------------- # + + self.check_Hermitian_x_y(uk, vk) + + self.uk = uk + self.ukh = np.transpose(np.conj(uk)) + + if self.EigV_positive == 1: # U == V + self.vkh = np.transpose(np.conj(uk)) + self.vk = uk + + else: + self.vkh = np.transpose(np.conj(vk)) + self.vk = vk + + self.Xk = self.uk.dot(self.sDiag) @ self.vkh + + +def Col_flip_sign(xCol): + """to change the sign of each column + such that the largest element of each column is positive. + + The purpose of this function is for comparison, such as that will be + used in check_Hermitian_x_y + + Args: + xCol (ndarray): the input column vector + + Returns: + ndarray: the update xCol with sign changed + """ + + if xCol.shape[0] < xCol.shape[1]: + print(" *** ERROR: This is a row, not a column ***") + return + + ID_max_abs = [np.argmax(np.abs(xCol[:, ii])) for ii in range(xCol.shape[1])] + sign = np.sign([xCol[ID_max_abs[ii], ii] for ii in range(xCol.shape[1])]) + + xCol = np.multiply(sign, xCol) + + return xCol + + +def Hr_Hermitian(X, r, Choice=0): + """hard theresholding the matrix X to rank r approximation + + Args: + X (ndarray): the input matrix + r (int): the target rank for the final approximated matrix + Choice (int, optional): the method of doing SVD. Defaults to 0. + + Returns: + ndarray: (u) the final left singular vector + ndarray: np.diag(s)[:r, :r] = the matrix with diagonal elements the singular values + ndarray: (vh) the complex conjugate of the right singular vector + """ + + print(" Choice for Hr_Hermitian = {}\n".format(Choice)) + + if Choice == 0: + u, s, vh = LA.svd(X, hermitian=True) + + u = np.array(u[:, :r]) # only keep the first r columns + vh = np.array(vh[:r, :]) # only keep the first r rows + elif Choice == 1: + u, s, vh = svds(X, k=r) + + return u, np.diag(s)[:r, :r], vh diff --git a/src/qibo/tomography_RGD/main_Tomography_RGD.py b/src/qibo/tomography_RGD/main_Tomography_RGD.py new file mode 100644 index 0000000000..c85e4b7bb7 --- /dev/null +++ b/src/qibo/tomography_RGD/main_Tomography_RGD.py @@ -0,0 +1,97 @@ +# import projectors + +from functools import reduce + +import numpy as np +import RGD +from BasicTools import Plt_Err_Time +from qibochem.measurement import expectation + +from qibo.hamiltonians import SymbolicHamiltonian +from qibo.models.encodings import ghz_state +from qibo.symbols import I, X, Y, Z + +symbol_map = { + "X": X, + "Y": Y, + "Z": Z, + "I": I, +} + + +def generate_random_label(n, symbols=["I", "X", "Y", "Z"]): + num_symbols = len(symbols) + label = "".join([symbols[i] for i in np.random.randint(0, num_symbols, size=n)]) + return label + + +if __name__ == "__main__": + + ############################################################ + ### Example of creating and running an experiment + ############################################################ + num_labels = 5 # number of Pauli labels + n = 3 # number of qubits + num_shots = 200 # number of shot measurements + + # generate circuit & shot measurements + # + Nr = 1 # rank of the target density matrix + + stateGHZ = ghz_state(n) # generate GHZ state circuit + target_state_GHZ = stateGHZ.execute().state() # get the state vector of the circuit + target_density_matrix = np.outer(target_state_GHZ, target_state_GHZ.conj()) + + # labels = [generate_random_label(n) for i in range(num_labels)] + + # labels = ["YXY", "IXX", "ZYI", "XXX", "YZZ", "ZYY", "IXX", "XYY", "XZI"] + # labels = ['YZY', 'YYZ', 'XIY', 'IZY', 'YYY','XZI','IXZ','IIY','XXY','YZZ'] + # labels = ['XXYY', 'YXZY', 'XZIZ', 'ZYYZ', 'YXYX', 'ZXXZ', 'ZIYZ', 'YXYY', 'YXYY', 'ZYYY'] + # labels = ['ZIYZ', 'YXZX', 'ZIXX', 'YXYX', 'ZYXX'] + labels = ["IIZ", "ZIZ", "IZI", "ZII", "IIZ", "IZZ"] + + Proj_list = [] + coef_exact = [] + coef_shots = [] + for label in labels: + # print(f"label = {label}") + + qubitPauli = [symbol_map[Ps](i) for i, Ps in enumerate(label)] + symbolPauli = SymbolicHamiltonian(reduce(lambda x, y: x * y, qubitPauli)) + + _circuit = stateGHZ.copy() + + coef_Pauli_exact = expectation(_circuit, symbolPauli) + if label == "III": + coef_Pauli_shots = 1.0 + else: + coef_Pauli_shots = symbolPauli.expectation_from_circuit( + _circuit, nshots=num_shots + ) + + Proj_list.append(symbolPauli) + coef_exact.append(coef_Pauli_exact) + coef_shots.append(coef_Pauli_shots.real) + + # + # system parameters as the input to the RGD worker + # + params_dict = { + "Nr": Nr, + "target_DM": target_density_matrix, + "labels": labels, + "measurement_list": coef_shots, # measurement_list, + "symProj_list": Proj_list, # symProj_list, + # "num_iterations": 150, + # "convergence_check_period": 1, + } + + Ch_svd = -1 # choice for initial SVD (0: LA.svd, 1: svds, -1: rSVD) + InitX_RGD = 1 # method of choosing initial X0 + + worker = RGD.BasicWorkerRGD( + target_density_matrix, labels, coef_shots, Proj_list, Nr + ) + worker.computeRGD(InitX_RGD, Ch_svd) + + Plt_Err_Time(worker) # plot the error and time evolution of the RGD algorithm diff --git a/src/qibo/tomography_RGD/main_sample.ipynb b/src/qibo/tomography_RGD/main_sample.ipynb new file mode 100644 index 0000000000..daa1db4e5e --- /dev/null +++ b/src/qibo/tomography_RGD/main_sample.ipynb @@ -0,0 +1,359 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "6fccc5ce", + "metadata": {}, + "source": [ + "# The sample codes for running RGD tomography\n", + "\n", + "## define the functions" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "6deb5a68", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import RGD\n", + "from BasicTools import Plt_Err_Time\n", + "\n", + "from functools import reduce\n", + "\n", + "from qibo.models.encodings import ghz_state\n", + "from qibo.hamiltonians import SymbolicHamiltonian\n", + "from qibochem.measurement import expectation\n", + "\n", + "\n", + "from qibo.symbols import I, X, Y, Z\n", + "\n", + "symbol_map = {\n", + " \"X\": X,\n", + " \"Y\": Y,\n", + " \"Z\": Z,\n", + " \"I\": I,\n", + "}\n", + "\n", + "def generate_random_label(n, symbols=[\"I\", \"X\", \"Y\", \"Z\"]):\n", + " num_symbols = len(symbols)\n", + " label = \"\".join([symbols[i] for i in np.random.randint(0, num_symbols, size=n)])\n", + " return label" + ] + }, + { + "cell_type": "markdown", + "id": "36a048f7", + "metadata": {}, + "source": [ + "## choose the Pauli string (operator) for doing tomography" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7b2befca", + "metadata": {}, + "outputs": [], + "source": [ + "num_labels = 5 # number of Pauli labels\n", + "n = 3 # number of qubits\n", + "# labels = [generate_random_label(n) for i in range(num_labels)]\n", + "\n", + "labels = [\"YXY\", \"IXX\", \"ZYI\", \"XXX\", \"YZZ\", \"ZYY\", \"IXX\", \"XYY\", \"XZI\"]\n", + "# labels = ['YZY', 'YYZ', 'XIY', 'IZY', 'YYY','XZI','IXZ','IIY','XXY','YZZ']\n", + "# labels = ['XXYY', 'YXZY', 'XZIZ', 'ZYYZ', 'YXYX', 'ZXXZ', 'ZIYZ', 'YXYY', 'YXYY', 'ZYYY']\n", + "# labels = ['ZIYZ', 'YXZX', 'ZIXX', 'YXYX', 'ZYXX']\n", + "# labels = ['IZZ', 'IIZ','ZIZ'] # NotImplementedError: Observable is not a Z Pauli string.\n" + ] + }, + { + "cell_type": "markdown", + "id": "1a350d88", + "metadata": {}, + "source": [ + "## generate circuit & shot measurements" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "b68f34e9", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "[Qibo 0.2.17|INFO|2025-07-28 23:08:52]: Using qibojit (numba) backend on /CPU:0\n" + ] + } + ], + "source": [ + "Nr = 1 # rank of the target density matrix\n", + "\n", + "stateGHZ = ghz_state(n) # generate GHZ state circuit\n", + "target_state_GHZ = stateGHZ.execute().state() # get the state vector of the circuit\n", + "target_density_matrix = np.outer(target_state_GHZ, target_state_GHZ.conj())\n", + "\n", + "num_shots = 200 # number of shot measurements" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "ace35edd", + "metadata": {}, + "outputs": [], + "source": [ + "Proj_list = []\n", + "coef_exact = []\n", + "coef_shots = []\n", + "for label in labels:\n", + " # print(f\"label = {label}\")\n", + "\n", + " qubitPauli = [symbol_map[Ps](i) for i, Ps in enumerate(label)]\n", + " symbolPauli = SymbolicHamiltonian(reduce(lambda x, y: x * y, qubitPauli))\n", + "\n", + " _circuit = stateGHZ.copy()\n", + "\n", + " coef_Pauli_exact = expectation(_circuit, symbolPauli)\n", + " if label == \"III\":\n", + " coef_Pauli_shots = 1.0\n", + " else:\n", + " coef_Pauli_shots = symbolPauli.expectation_from_circuit(_circuit, nshots=num_shots)\n", + "\n", + " Proj_list.append(symbolPauli)\n", + " coef_exact.append(coef_Pauli_exact)\n", + " coef_shots.append(coef_Pauli_shots.real)\n" + ] + }, + { + "cell_type": "markdown", + "id": "4fad7418", + "metadata": {}, + "source": [ + "## system parameters as the input to the RGD worker" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "cd20c28a", + "metadata": {}, + "outputs": [], + "source": [ + "params_dict = {\n", + " \"Nr\": Nr,\n", + " \"target_DM\": target_density_matrix,\n", + " \"labels\": labels,\n", + " \"measurement_list\": coef_shots, #measurement_list,\n", + " \"symProj_list\": Proj_list, #symProj_list,\n", + " # \"num_iterations\": 150,\n", + " \"convergence_check_period\": 1,\n", + "}\n", + "\n", + "Ch_svd = -1 # choice for initial SVD (0: LA.svd, 1: svds, -1: rSVD)\n", + "InitX_RGD = 1 # method of choosing initial X0" + ] + }, + { + "cell_type": "markdown", + "id": "8ee4e9af", + "metadata": {}, + "source": [ + "## running the RGD once collecting the input information" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "259cea35", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "[Qibo 0.2.17|WARNING|2025-07-28 23:09:58]: Calculating the dense form of a symbolic Hamiltonian. This operation is memory inefficient.\n", + "[Qibo 0.2.17|WARNING|2025-07-28 23:09:58]: Calculating the dense form of a symbolic Hamiltonian. This operation is memory inefficient.\n", + "[Qibo 0.2.17|WARNING|2025-07-28 23:09:58]: Calculating the dense form of a symbolic Hamiltonian. This operation is memory inefficient.\n", + "[Qibo 0.2.17|WARNING|2025-07-28 23:09:58]: Calculating the dense form of a symbolic Hamiltonian. This operation is memory inefficient.\n", + "[Qibo 0.2.17|WARNING|2025-07-28 23:09:58]: Calculating the dense form of a symbolic Hamiltonian. This operation is memory inefficient.\n", + "[Qibo 0.2.17|WARNING|2025-07-28 23:09:58]: Calculating the dense form of a symbolic Hamiltonian. This operation is memory inefficient.\n", + "[Qibo 0.2.17|WARNING|2025-07-28 23:09:58]: Calculating the dense form of a symbolic Hamiltonian. This operation is memory inefficient.\n", + "[Qibo 0.2.17|WARNING|2025-07-28 23:09:58]: Calculating the dense form of a symbolic Hamiltonian. This operation is memory inefficient.\n", + "[Qibo 0.2.17|WARNING|2025-07-28 23:09:58]: Calculating the dense form of a symbolic Hamiltonian. This operation is memory inefficient.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "num_iterations = 150\n", + " self.coef = 0.9428090415820634\n", + " InitX = 1\n", + " ********* using randomized-SVD to construct X0 = uk @ sDiag @ vkh **************************\n", + " *** 0-th (A A*) done: Time --> ATQ: 6.198883056640625e-05, AG: 5.698204040527344e-05 ***\n", + " *** 1-th (A A*) done: Time --> ATQ: 4.9114227294921875e-05, AG: 4.506111145019531e-05 ***\n", + " *** 2-th (A A*) done: Time --> ATQ: 4.315376281738281e-05, AG: 4.1961669921875e-05 ***\n", + " *** 3-th (A A*) done: Time --> ATQ: 5.316734313964844e-05, AG: 4.363059997558594e-05 ***\n", + " *** 4-th (A A*) done: Time --> ATQ: 4.291534423828125e-05, AG: 4.00543212890625e-05 ***\n", + " *** 5-th (A A*) done: Time --> ATQ: 4.100799560546875e-05, AG: 4.00543212890625e-05 ***\n", + " *** 6-th (A A*) done: Time --> ATQ: 3.886222839355469e-05, AG: 3.910064697265625e-05 ***\n", + " *** 7-th (A A*) done: Time --> ATQ: 4.100799560546875e-05, AG: 4.00543212890625e-05 ***\n", + " *** 8-th (A A*) done: Time --> ATQ: 3.981590270996094e-05, AG: 4.100799560546875e-05 ***\n", + " *** 9-th (A A*) done: Time --> ATQ: 3.981590270996094e-05, AG: 3.910064697265625e-05 ***\n", + " *** 10-th (A A*) done: Time --> ATQ: 3.910064697265625e-05, AG: 3.886222839355469e-05 ***\n", + " *** 11-th (A A*) done: Time --> ATQ: 3.933906555175781e-05, AG: 3.886222839355469e-05 ***\n", + " *** 12-th (A A*) done: Time --> ATQ: 4.00543212890625e-05, AG: 3.790855407714844e-05 ***\n", + " *** 13-th (A A*) done: Time --> ATQ: 3.981590270996094e-05, AG: 3.933906555175781e-05 ***\n", + " *** 14-th (A A*) done: Time --> ATQ: 3.910064697265625e-05, AG: 3.886222839355469e-05 ***\n", + " self.rSVD --> time = 0.00697016716003418\n", + " -------- X0 --> (-1)-th iteration [ SVD choice = -1 ] \n", + "\n", + " X0 step (RGD) --> time = 0.007112026214599609\n", + "\n", + " X0 --> Tr(X0) = (2.6806883554733316+4.593228591279375e-17j)\n", + "\n", + " X0 --> Target Err = 1.6920764698162736\n", + "\n", + " RGD max num_iterations = 150\n", + "\n", + " -------- 0-th iteration \n", + "\n", + " 0-th step -> self.Remain_Hermitian = 1, self.EigV_positive = 1\n", + " Hr_tangentWk --> Remain_Hermitian = 1, EigV_positve = 1\n", + " D2s - D2s.T.conj() = 0.0\n", + " np.allclose(D2s, D2s.T.conj()) = True\n", + "\n", + " **** 0-th iteratation\n", + "\n", + " min(uk +/- uk_old) = 0.09935648824662649\n", + "\n", + " relative_errU = 0.0993564882466265, relative_errX = 0.6216376939638819\n", + " Target_Error = 0.11009543694651341\n", + "\n", + " convergence check --> time = 9.870529174804688e-05\n", + "\n", + " stepRGD --> time = 0.0010688304901123047\n", + "\n", + " -------- 1-th iteration \n", + "\n", + " 1-th step -> self.Remain_Hermitian = 1, self.EigV_positive = 1\n", + " Hr_tangentWk --> Remain_Hermitian = 1, EigV_positve = 1\n", + " D2s - D2s.T.conj() = 9.811794782561298e-19\n", + " np.allclose(D2s, D2s.T.conj()) = True\n", + "\n", + " min(uk +/- uk_old) = 0.06974724727423012\n", + "\n", + " relative_errU = 0.06974724727423015, relative_errX = 0.09960414871315858\n", + " Target_Error = 0.09774914657106261\n", + "\n", + " convergence check --> time = 5.817413330078125e-05\n", + "\n", + " stepRGD --> time = 0.0008833408355712891\n", + "\n", + " -------- 2-th iteration \n", + "\n", + " 2-th step -> self.Remain_Hermitian = 1, self.EigV_positive = 1\n", + " Hr_tangentWk --> Remain_Hermitian = 1, EigV_positve = 1\n", + " D2s - D2s.T.conj() = 6.706169350291565e-19\n", + " np.allclose(D2s, D2s.T.conj()) = True\n", + "\n", + " min(uk +/- uk_old) = 0.0013605475291729539\n", + "\n", + " relative_errU = 0.0013605475291729545, relative_errX = 0.006773174987836168\n", + " Target_Error = 0.09638833925949807\n", + "\n", + " convergence check --> time = 5.1975250244140625e-05\n", + "\n", + " stepRGD --> time = 0.0008130073547363281\n", + "\n", + " -------- 3-th iteration \n", + "\n", + " 3-th step -> self.Remain_Hermitian = 1, self.EigV_positive = 1\n", + " Hr_tangentWk --> Remain_Hermitian = 1, EigV_positve = 1\n", + " D2s - D2s.T.conj() = 0.0\n", + " np.allclose(D2s, D2s.T.conj()) = True\n", + "\n", + " min(uk +/- uk_old) = 0.00076725236824781\n", + "\n", + " ********* XkErrRatio = 0.001120030733734077 < StpTol ******\n", + " relative_errU = 0.0007672523682478105, relative_errX = 0.001120030733734077\n", + " Target_Error = 0.09638335176146144\n", + "\n", + " convergence check --> time = 5.2928924560546875e-05\n", + "\n", + " stepRGD --> time = 0.0007750988006591797\n", + "\n" + ] + } + ], + "source": [ + "worker = RGD.BasicWorkerRGD(params_dict)\n", + "worker.computeRGD(InitX_RGD, Ch_svd)" + ] + }, + { + "cell_type": "markdown", + "id": "aacd934b", + "metadata": {}, + "source": [ + "## plot the error and time evolution of the RGD algorithm" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "b8de7e36", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAArsAAAIpCAYAAACizfDKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABxZElEQVR4nO3dd3xUVd4G8OfOJJM+kwLpjRIgtBAChAhICwTECDZaKLI2FMuKNbtS91VQwbKKuotSlA4KoiCI9BogEEAIPZCQSksnbea8f4TMMqSQMsmdTJ7vZ8dl7tzym5MhPDk551xJCCFARERERGSGFHIXQERERERUXxh2iYiIiMhsMewSERERkdli2CUiIiIis8WwS0RERERmi2GXiIiIiMwWwy4RERERmS2GXSIiIiIyWwy7RERERGS2GHaJiMjs7Nq1C5IkYebMmXKXQkQyY9glogZ15coVSJJU5cPf31/uMhstSZLQr18/uctoEE3pvRJR7VnIXQARNU2tWrXCuHHjKnzN0dGxYYshs9OjRw/Ex8ejWbNmcpdCRDJj2CUiWbRu3Zq/YqZ6Y2tri3bt2sldBhGZAA5jICKTV/br6uTkZEyYMAHu7u5QKBTYtWuXwdjMAwcOYPDgwXB0dIQkSfrj8/LyMGPGDLRr1w7W1tZwdnbGsGHDsH///nLXmjlzJiRJwq5du7BkyRJ07doVtra2Vf66XKfTwcXFBR07djTYfuvWLSgUCkiShD///NPgtWeeeQaSJOHq1auVnrdsn8uXL2P+/Plo3749rKys8Mwzz5Tbt6wdAGD37t0Gw0KWLFlS6TUexN/fH/7+/sjMzMQrr7wCHx8fWFhY6M9Z9npF+vXrZ/B1AAzbd8WKFejSpQtsbGzg4eGB119/HXfu3HlgTdV5r5WN2S2rNysrCy+99BI8PDxgZ2eHhx9+GMeOHQMApKSkYNy4cXB1dYWNjQ0GDx6MCxcuVFhLQkICnnvuOfj6+sLKygoeHh545plnqvy6ElHDYs8uETUKN2/eRFhYGJydnTF69GgUFBRArVYjOzsbAHDgwAF8+OGH6N+/P1544QUkJiYCAAoKCjBgwAAcPnwYXbt2xd///nekp6dj9erV2Lp1K1auXImnn3663PU++eQT7Ny5E8OHD8fgwYOhVCorrU2hUKBv375Yv349MjIy4OrqCqA0iAkhAAA7d+5EeHi4/pidO3eiRYsW8PPze+B7f/XVV3Ho0CEMGzYMkZGR+vPfy9/fHzNmzMCsWbPg5+dnEIi7dOnywGtUpbCwEAMGDEBubi4ee+wxWFhYwM3NrU7n/Oqrr7BlyxYMHz4cAwYMwJYtW/Dvf/8bN27cwPLly6s8tq7vtaioCIMGDUJBQQFGjRqF9PR0rFmzBuHh4Thw4AAiIiLg4eGBcePG4eLFi/j1118xbNgwxMfHG3wOYmJiEBERgby8PDz66KMICAjAlStXsHz5cvz+++84ePAgWrZsWdsmIiJjEUREDSghIUEAEK1atRIzZsyo8PH7778bHANAABCTJk0SJSUlBq/t3LlT//qiRYvKXW/WrFkCgIiKihI6nU6//dixY0KlUglHR0eRnZ2t3z5jxgwBQNjZ2YmTJ09W+339+9//FgDE6tWr9dteffVVYWdnJ3r27CnCwsL02y9duiQAiL/97W9VnnPixIkCgPD29hZXr16tVh0ARN++fatd94P4+fkJACIiIkLk5+dX+Lqfn1+Fx/bt21fc/89MWftqNBpx9uxZ/fb8/HzRpk0boVAoRHJycrVqq+q9ln0uZsyYUeH7efrpp0VxcbF++0cffSQACEdHR/HGG28YfFZeeuklAUD89NNP+m1FRUXC399fODg4iGPHjhlcY+/evUKpVIpHH320Wu+DiOoXhzEQkSwuXbqEWbNmVfjYsmVLuf1VKhU+/vjjSntYu3btikmTJpXbvnTpUlhaWmLu3LkGv1IPDg7GxIkTkZmZiQ0bNpQ77oUXXkCnTp2q/X769+8PANixY4d+286dO9G7d28MHjwYR44cQW5urn47gGqvJPD222/D19e32rXUh48//hg2NjZGO9/rr7+Otm3b6p/b2NhgzJgx0Ol0iI2NNdp1KjNv3jxYWPzvl5tjxowBAJSUlOD//u//DD4rZa+dOHFCv+23337DlStX8PbbbyM4ONjg3L1798bw4cOxefNm/W8eiEg+HMZARLKIiIioMNRWpkWLFlXOrO/evXu5bdnZ2bh8+TICAwPh7e1d7vX+/ftj4cKFiIuLw/jx4w1e69GjR7VrA4AOHTqgefPm+iB7/fp1nD59GuPHj0ePHj0we/Zs7N27F0OHDtXvUxaQH6SmtRibtbV1jYJ/dYSEhJTbVvY1yszMNOq17ufk5FTuhwcPDw8AQEBAAGxtbSt8LSUlRb/t0KFDAIBz585VONEyLS0NOp0O58+fR7du3YxZPhHVEMMuETUKDxojWtHrZb1qlR1bFmIq6n2r6ZjUskl0a9euRUpKCvbv3w8hBAYMGIBOnTrB2toaO3fuxNChQ7Fr1y60bt26wgBekbqOj60rV1fXchPN6kqtVpfbVtbTqtVqjXqtmly7qteKi4v1227dugUADxxfnJeXV+s6icg4OIyBiBqFB4Wtil4vCy7p6ekVHpOWlmawX02uV5GyntqdO3di165d0Gg0CA4OhpWVFcLCwrBz505cuHABycnJ1e7VrW0txlTV9RUKBUpKSip8LSsrq75Kkl3ZZ+bXX3+FEKLSR9++fWWulIgYdonIbKnVarRs2RIXL15EcnJyudd37doFoO6rFZS5d9zuzp070bdvX/0Y4wEDBuD48eNYv349gOqP160phUJR7z2j93JyckJGRka5wJuXl1fpcl3G0tDv9V6hoaEAgIMHD8pyfSKqPoZdIjJrEydORHFxMaKjo/XLgAHAyZMnsWTJEmg0GowYMcIo12rXrh3c3d3x66+/Ij4+HgMGDNC/1r9/f2i1WsybN0//vMyNGzdw9uxZ3Lhxo9rXunTpEs6ePWvwq3UAcHZ2xrVr1yo97uzZszh79my1r/Mg3bt3R3FxscGv84UQiI6Orvdf4T/ovdan4cOHw9fXF59++in27NlT7vXi4mLs27dPhsqI6H4cs0tEsrh48WKVd1B77733YG1tXefrvPPOO9i0aRN+/PFHxMfHY+DAgcjIyMDq1atRUlKChQsXwsHBoc7XKdO/f3+sXLlS/+cyPXr0gJ2dHa5fv462bdvqxwsDpWvOzpo1CzNmzKj2XeUGDhyIq1evIiEhweCmDgMGDMCaNWswYsQIBAcHQ6lU4rHHHkPnzp0BAIGBgQBgEPzr4pVXXsHixYvx3HPPYdu2bWjevDn27t2LzMxMBAUFGaxgYGwPeq/1ycrKCuvWrcPQoUPRt29f/djsshuF7N27Fy4uLkb9wYKIaodhl4hkUbb0WGX+/ve/GyXsWltbY8eOHfjoo4+wevVqfPbZZ7C1tUXfvn3xj3/8A717967zNe5VFnabNWtmsIKBpaUlevXqhT/++KPehjAAwBdffAGgdCjFr7/+Cp1OB29v73oLgB07dsSWLVsQHR2NdevWwd7eHo888gjmzZuHkSNH1ss1yzT0e71f9+7dceLECXzyySfYvHkz9u/fDysrK3h5eWHEiBH6JcuISF6SMNaP90REREREJoZjdomIiIjIbDHsEhEREZHZYtglIiIiIrPFsEtEREREZothl4iIiIjMFsMuEREREZktrrN7H51Oh5SUFDg4OMh+P3oiIiIiKk8IgZycHHh6ekKhqLrvlmH3PikpKfDx8ZG7DCIiIiJ6gKSkJHh7e1e5D8PufcpuG5qUlAS1Wi1zNURERER0v+zsbPj4+FTrdu8Mu/cpG7qgVqsZdomIiIhMWHWGnHKCGhERERGZLYZdIiIiIjJbDLtEREREZLY4ZrcWhBAoKSmBVquVuxQiMjFKpRIWFhZcupCIyEQw7NZQUVERUlNTkZ+fL3cpRGSibG1t4eHhAZVKJXcpRERNHsNuDeh0OiQkJECpVMLT0xMqlYq9N0SkJ4RAUVERrl+/joSEBAQEBDxwsXMiIqpfDLs1UFRUBJ1OBx8fH9ja2spdDhGZIBsbG1haWuLq1asoKiqCtbW13CURETVp7HKoBfbUEFFV+D2CiMh08DsyEREREZkthl0iIiIiMlsMu1TvJEnChg0bKn39ypUrkCQJcXFxDVbT/fz9/fH5559Xuc/MmTPRpUuXBqmnvu3atQuSJCEzMxMAsGTJEjg6OspaU2Xur7U6+vXrh7///e91uq4ptwkREVUfw66MMrIL8Nm288jILqj3a33zzTfo3Lkz1Go11Go1wsLC8Pvvv9f7davDx8cHqamp6NixI4DahZu6OnLkCF544QX98wcFdHMzatQonD9/vkGutWDBAgQGBsLGxgZt27bFDz/80CDXJSKipomrMcgoI6cQX2y/gEHt3eCqrt8Z297e3pg7dy4CAgIghMDSpUsxfPhwHD9+HB06dKjXaz+IUqmEu7u7rDU0b95c1uvLzcbGBjY2NvV+nW+++QbR0dFYuHAhunfvjsOHD+P555+Hk5MTIiMj6/36RETU9LBnt46EEMgvKqnVo6C49A5sBcXaWh0vhKh2nZGRkXjkkUcQEBCANm3a4IMPPoC9vT0OHTpU6TFarRZTp06Fo6MjXFxc8M4772DixIkYMWKEfp+Kfv3fpUsXzJw502Bbamoqhg4dChsbG7Rs2RLr1q3Tv3bvMIYrV66gf//+AAAnJydIkoRnnnkGALBu3Tp06tQJNjY2cHFxQXh4OPLy8iqsvVu3bpg3b57++YgRI2BpaYnc3FwAwLVr1yBJEi5evFjuffj7+wMAHn/8cUiSpH9e5scff4S/vz80Gg1Gjx6NnJycStvw6tWriIyMhJOTE+zs7NChQwds3rxZ377PPvssWrRooe/l/OKLLwyOf+aZZzBixAh8+OGHcHNzg6OjI2bPno2SkhK8/fbbcHZ2hre3NxYvXlyuPVetWoWHHnoI1tbW6NixI3bv3l1pnff/yr5syEZV7zUnJwdRUVGws7ODh4cHPvvsswcOH/jxxx/x4osvYtSoUWjZsiVGjx6NF154AR999FGlx9zv5s2bGDNmDLy8vGBra4tOnTph5cqV5fYrKSnBK6+8Ao1Gg2bNmmHatGkGf2cKCwvx1ltvwcvLC3Z2dggNDcWuXbsqve6JEyfQv39/ODg4QK1WIyQkBEePHq123UREJA/27NbRnWIt2k/fWqdzPPXtwVodd2Z2BGxVNf8SarVarF27Fnl5eQgLC6t0v/nz52PJkiVYtGgRAgMDMX/+fKxfvx4DBgyo8TWnTZuGuXPn4osvvsCPP/6I0aNH49SpUwgMDDTYz8fHBz/99BOefPJJnDt3Dmq1GjY2NkhNTcWYMWPw8ccf4/HHH0dOTg727t1baeDv27cvdu3ahbfeegtCCOzduxeOjo7Yt28fhgwZgt27d8PLywutW7cud+yRI0fg6uqKxYsXY8iQIVAqlfrXLl26hA0bNuC3337D7du3MXLkSMydOxcffPBBhXVMmTIFRUVF2LNnD+zs7HDmzBnY29sDKL1Jibe3N9auXQsXFxccOHAAL7zwAjw8PDBy5Ej9OXbs2AFvb2/s2bMH+/fvx7PPPosDBw7g4YcfRkxMDFavXo0XX3wRgwYNgre3t/64t99+G59//jnat2+PTz/9FJGRkUhISICLi0u1vmYPeq9Tp07F/v37sXHjRri5uWH69Ok4duxYleOaCwsLy607a2Njg8OHD6O4uBiWlpYPrKugoAAhISF49913oVarsWnTJowfPx6tWrVCjx499PstXboUzz77LA4fPoyjR4/ihRdegK+vL55//nkAwCuvvIIzZ85g1apV8PT0xPr16zFkyBCcOnUKAQEB5a4bFRWF4OBgfPPNN1AqlYiLi6tWvUREJC/27DYhp06dgr29PaysrDB58mSsX78e7du3r3T/zz//HNHR0XjiiScQGBiIb7/9FhqNplbXfvrpp/Hcc8+hTZs2+Ne//oVu3brhyy+/LLefUqmEs7MzAMDV1RXu7u7QaDRITU1FSUkJnnjiCfj7+6NTp054+eWX9cHxfv369cO+ffug1Wpx8uRJqFQqREVF6Xvudu3ahb59+1Z4bNmQBkdHR7i7uxsMcdDpdFiyZAk6duyIPn36YPz48di+fXul7zsxMRG9evVCp06d0LJlSzz66KN4+OGHAQCWlpaYNWsWunXrhhYtWiAqKgqTJk3CmjVrDM7h7OyMf//732jbti3+9re/oW3btsjPz8c//vEPBAQEIDo6GiqVCvv27TM47pVXXsGTTz6JwMBAfPPNN9BoNPj+++8rrfV+Vb3XnJwcLF26FPPmzcPAgQPRsWNHLF68GFqttspzRkRE4LvvvkNsbCyEEDh69Ci+++47FBcX48aNG9Wqy8vLC2+99Ra6dOmCli1b4tVXX8WQIUPKtZuPjw8+++wztG3bFlFRUXj11Vfx2WefASj9uixevBhr165Fnz590KpVK7z11lvo3bu3QS/5vRITExEeHo527dohICAATz/9NIKCgqpVMxERyYc9u3VkY6nEmdkR1d7/ek4hrucUAgDOpGZj+i+n8dbgNmhmbwVvJ2v4ONuhuYNVta9dE23btkVcXByysrKwbt06TJw4Ebt3764w8GZlZSE1NRWhoaH6bRYWFujWrVuNhk+Uub8HOSwsrEarLwQFBWHgwIHo1KkTIiIiMHjwYDz11FNwcnKqcP8+ffogJycHx48fx4EDB9C3b1/069cPc+fOBQDs3r0bb7/9do3fh7+/PxwcHPTPPTw8kJGRUen+r732Gl566SX88ccfCA8Px5NPPonOnTvrX1+wYAEWLVqExMRE3LlzB0VFReV6Rjt06GBwkwI3Nzf9ZD6g9AcEFxeXcnXc2+ZlX7v4+HijvNfLly+juLjYoCdVo9Ggbdu2VZ5z2rRpSEtLQ8+ePSGEgJubGyZOnIiPP/642jdi0Gq1+PDDD7FmzRokJyejqKgIhYWF5e5q2LNnT4PbeYeFhWH+/PnQarU4deoUtFot2rRpY3BMYWFhpT3fU6dOxXPPPYcff/wR4eHhePrpp9GqVatq1UxERPJhz24dSZIEW5VFtR9+Lnbo5u+MIB9HtPdQAwC8HG3Q3lON9p4auDpYQQJgqVQ88Fz3/kNeHSqVCq1bt0ZISAjmzJmDoKCgcmNEa0qhUJQLv8XFxXU6Z0WUSiW2bduG33//He3bt8eXX36Jtm3bIiEhocL9HR0dERQUhF27dmH37t3o168fHn74YRw/fhznz5/HhQsXKu3Zrcr9v7aWJAk6na7S/Z977jlcvnwZ48ePx6lTpwx6tFetWoW33noLzz77LP744w/ExcVh0qRJKCoqeuA1a1pHbdTHNWxsbLBo0SLk5+fjypUrSExM1Ifq6k4S/OSTT/DFF1/g3Xffxc6dOxEXF4eIiIhy7VaV3NxcKJVKxMbGIi4uTv+Ij4+v9O/EzJkzcfr0aQwbNgw7duxA+/btsX79+mpfk4iI5MGwK5NbeUVIupVvsO3a7Tu4kJGLCxm5uJVX/X+4a0un06GwsLDC1zQaDTw8PBATE6PfVlJSgtjYWIP9mjdvjtTUVP3z7OzsCgPo/RPhDh06VG68bhmVSgUA5X4lLkkSevXqhVmzZuH48eNQqVRVho2+ffti586d2LNnD/r16wdnZ2cEBgbigw8+gIeHR7levXtZWlo+8Ffy1eXj44PJkyfj559/xptvvomFCxcCAPbv34+HHnoIL7/8MoKDg9G6dWtcunTJKNcEDNu87GtXWZvXVMuWLWFpaYkjR47ot2VlZVV7+TJLS0t4e3tDqVRi1apVePTRR6vds7t//34MHz4c48aNQ1BQEFq2bFnhde/97AKl7REQEAClUong4GBotVpkZGSgdevWBo+qVgZp06YN3njjDfzxxx944oknKh3yQEREpoPDGGTibKdCVz8nvNyvlX7ZMUcbS/0QBgulcX8OiY6OxtChQ+Hr64ucnBysWLECu3btwtatlU+ue/311/XLlbVr1w6ffvppubVvBwwYgCVLliAyMhKOjo6YPn26wYSuMmvXrkW3bt3Qu3dvLF++HIcPH650/Kifnx8kScJvv/2GRx55BDY2Njh9+jS2b9+OwYMHw9XVFTExMbh+/XqV4a1fv3748ssv0bx5c7Rr106/7auvvsLTTz9dZXv5+/tj+/bt6NWrF6ysrCodLvEgf//73zF06FC0adMGt2/fxs6dO/U1BwQE4IcffsDWrVvRokUL/Pjjjzhy5AhatGhRq2vdb8GCBQgICEBgYCA+++wz3L59G3/729+Mcm4HBwdMnDhRvyKEq6srZsyYAYVCUeVvHM6fP4/Dhw8jNDQUt2/fxqeffoq//voLS5curfa1AwICsG7dOhw4cABOTk749NNPkZ6eXm44TmJiIqZOnYoXX3wRx44dw5dffon58+cDKA2tUVFRmDBhAubPn4/g4GBcv34d27dvR+fOnTFs2DCDc925cwdvv/02nnrqKbRo0QLXrl3DkSNH8OSTT9ag1YiISA7s2ZWJpVIBPxc7vDOkHVo2twMA5BdpYW2phI3KApZGDrsZGRmYMGEC2rZti4EDB+LIkSPYunUrBg0aVOkxb775JsaPH4+JEyciLCwMDg4OePzxxw32iY6ORt++ffHoo49i2LBhGDFiRIXjGGfNmoVVq1ahc+fO+OGHH7By5cpKJ8d5eXlh1qxZeO+99+Dm5oZXXnkFarUae/bswSOPPII2bdrg/fffx/z58zF06NBK6+/Tpw90Op3BcIV+/fpBq9WiX79+VbbX/PnzsW3bNvj4+CA4OLjKfaui1WoxZcoUBAYGYsiQIWjTpg2+/vprAMCLL76IJ554AqNGjUJoaChu3ryJl19+udbXut/cuXMxd+5cBAUFYd++fdi4cSOaNWtmtPN/+umnCAsLw6OPPorw8HD06tULgYGB5VZbuJdWq8X8+fMRFBSEQYMGoaCgAAcOHCi3vFtV3n//fXTt2hURERHo168f3N3dDZbDKzNhwgTcuXMHPXr0wJQpU/D6668b3Dhk8eLFmDBhAt588020bdsWI0aMwJEjR+Dr61vuXEqlEjdv3sSECRPQpk0bjBw5EkOHDsWsWbOqXTcREclDErWZbWTGsrOzodFokJWVBbVabfBaQUEBEhIS0KJFiyr/Qa+pvMISXLpeuv5ry2Z2sLc23eWMnnnmGWRmZjapu4s1NleuXEGLFi1w/PjxBr29cV5eHry8vDB//nw8++yzDXZdU1Rf3yuIiKhUVXntfhzGYAJUFgrYWCpxp1iLm3lFJh12icocP34cZ8+eRY8ePZCVlYXZs2cDAIYPHy5zZURERP/DsGsCLJUKeDvZ4kJGDrLvlKBYqzP6MAai+jBv3jycO3cOKpUKISEh2Lt3r1GHShAREdUVw66JsFEpYauyQH5RCW7nFeknrZmaJUuWyF0CPYC/v3+t1kKuqeDg4HKrcxAREZkadh+aEBe70iW3buUVNUhYISIiIjJ3DLu1UF9BVGNjCaVCQpFWh5yCknq5BhHVP/6wSkRkOhh2a6DsjlL5+fkP2LN2FAoJzralvbs3G+CmEkRUP8q+R9x/FzoiImp4HLNbA0qlEo6OjsjIyAAA2Nra1viWvQ9ip9Qho6QI2blFyLYGVBblb9BARKZJCIH8/HxkZGTA0dGxwhusEBFRw2LYraGyW4mWBd76kJ1TiIISHe7csoDGhj1DRI2No6NjlbcdJiKihsOwW0OSJMHDwwOurq4oLi6ul2tcO5+Bmb+egaOtJVY+HwaVBUebEDUWlpaW7NElIjIhDLu1pFQq6+0ftH4dvFH863mcTi/A7kuZiAzyrJfrEBEREZk7dhmaIEulAqN7+AIAlh26KnM1RERERI0Xw66JGt3dBwoJiEm4hQvpOXKXQ0RERNQoMeyaKE9HGwwMdAMALI9JlLkaIiIiosaJYdeEjevpBwD46dg15BfxJhNERERENcWwa8L6tG4GX2db5BSU4LcTqXKXQ0RERNToMOyaMIVCwtjQuxPVYjhRjYiIiKimTDrs7tmzB5GRkfD09IQkSdiwYcMDjyksLMQ///lP+Pn5wcrKCv7+/li0aFH9F1tPng7xhkqpwMlrWTh5LVPucoiIiIgaFZMOu3l5eQgKCsKCBQuqfczIkSOxfft2fP/99zh37hxWrlyJtm3b1mOV9cvF3gpDO5XeiWn5IU5UIyIiIqoJk76pxNChQzF06NBq779lyxbs3r0bly9fhrOzMwDA39+/nqprOON6+uGXuBT8ciIZ/xgWyFsIExEREVWTSffs1tTGjRvRrVs3fPzxx/Dy8kKbNm3w1ltv4c6dO5UeU1hYiOzsbIOHqenm54Q2bvYoKNbh52PX5C6HiIiIqNEwq7B7+fJl7Nu3D3/99RfWr1+Pzz//HOvWrcPLL79c6TFz5syBRqPRP3x8fBqw4uqRJEm/DNnymEQIIWSuiIiIiKhxMKuwq9PpIEkSli9fjh49euCRRx7Bp59+iqVLl1bauxsdHY2srCz9IykpqYGrrp7Hg71gq1LiYkYuYhJuyV0OERERUaNgVmHXw8MDXl5e0Gg0+m2BgYEQQuDatYp//W9lZQW1Wm3wMEUO1pYY3sUTAO+oRkRERFRdZhV2e/XqhZSUFOTm5uq3nT9/HgqFAt7e3jJWZhxRoaVDGbb8lYrrOYUyV0NERERk+kw67Obm5iIuLg5xcXEAgISEBMTFxSExsbRnMzo6GhMmTNDvP3bsWLi4uGDSpEk4c+YM9uzZg7fffht/+9vfYGNjI8dbMKqOXhoE+TiiWCuw5qhpDrcgIiIiMiUmHXaPHj2K4OBgBAcHAwCmTp2K4OBgTJ8+HQCQmpqqD74AYG9vj23btiEzMxPdunVDVFQUIiMj8e9//1uW+uvDuLt3VFt5OBFaHSeqEREREVVFEpzabyA7OxsajQZZWVkmOX63oFiLHh/8ieyCEix+pjv6t3OVuyQiIiKiBlWTvGbSPbtUnrWlEk+FlC6PtuzQVZmrISIiIjJtDLuNUFTP0qEMO85l4NrtfJmrISIiIjJdDLuNUKvm9niolQuEAFYd5kQ1IiIiosow7DZSZcuQrTqShGKtTuZqiIiIiEwTw24jNbiDG5o7WOFGbiH+OJ0udzlEREREJolht5GyVCowqhsnqhERERFVhWG3ERsT6guFBBy8fBMXM3IffAARERFRE8Ow24h5OdpgwN11dlfEJD5gbyIiIqKmh2G3kSubqLYuNgkFxVqZqyEiIiIyLQy7jdzDbZrD28kG2QUl+PVEitzlEBEREZkUht1GTqmQMDa09CYTyziUgYiIiMgAw64ZGNnNB5ZKCSeSMvFXcpbc5RARERGZDIZdM9DM3gpDOnoAAJbHcBkyIiIiojIMu2Yi6u5Qhg3HU5BdUCxzNURERESmgWHXTIS2cEaAqz3uFGux4Xiy3OUQERERmQSGXTMhSZK+d3fZoasQQshcEREREZH8GHbNyONdvWFjqcT59FwcvXpb7nKIiIiIZMewa0Y0NpZ4LMgTQGnvLhEREVFTx7BrZqJ6lg5l+P1UGm7mFspcDREREZG8GHbNTGdvR3T21qBIq8Pa2Gtyl0NEREQkK4ZdMzQu1A8AsCImETodJ6oRERFR08Wwa4YeDfKAg7UFEm/lY8+F63KXQ0RERCQbhl0zZKuywJNdvQEAy2MSZa6GiIiISD4Mu2Zq3N2Jatvj05GSeUfmaoiIiIjkwbBrplq7OiC0hTN0Alh1JEnucoiIiIhkwbBrxsb1LJ2otupwIoq1OpmrISIiImp4DLtmLKKDO5rZq5CRU4g/z6TLXQ4RERFRg2PYNWMqCwVGdvMBwIlqRERE1DQx7Jq5MT18IUnAvos3kHAjT+5yiIiIiBoUw66Z83G2Rb82zQEAK2KuylwNERERUcNi2G0CyiaqrY29hoJirczVEBERETUcht0moF9bV3g52iAzvxibTqbKXQ4RERFRg2HYbQKUCgljepRNVONQBiIiImo6GHabiJHdfWChkHAsMRNnUrLlLoeIiIioQTDsNhGuDtaI6OAOAFjG3l0iIiJqIhh2m5Conr4AgF+OJyO3sETmaoiIiIjqH8NuExLW0gUtm9shr0iL9ceT5S6HiIiIqN4x7DYhkiQhKrR0GbLlh65CCCFzRURERET1i2G3iXmqqzesLBQ4m5aDY4m35S6HiIiIqF4x7DYxGltLRAZ5AgCWHUqUuRoiIiKi+sWw2wSV3VFt06lU3MorkrkaIiIiovrDsNsEBXlr0NFLjaISHdbFJsldDhEREVG9Memwu2fPHkRGRsLT0xOSJGHDhg3VPnb//v2wsLBAly5d6q2+xspgolpMInQ6TlQjIiIi82TSYTcvLw9BQUFYsGBBjY7LzMzEhAkTMHDgwHqqrPEb3sUTDlYWuHozH/sv3ZC7HCIiIqJ6YSF3AVUZOnQohg4dWuPjJk+ejLFjx0KpVNaoN7gpsVVZ4ImuXlh68CqWHbqKPgHN5S6JiIiIyOhMume3NhYvXozLly9jxowZ1dq/sLAQ2dnZBo+mIuruRLU/4zOQllUgczVERERExmdWYffChQt47733sGzZMlhYVK/Tes6cOdBoNPqHj49PPVdpOtq4OaCHvzO0OoFVR7gMGREREZkfswm7Wq0WY8eOxaxZs9CmTZtqHxcdHY2srCz9Iympaa1OENXTFwCw6nASSrQ6mashIiIiMi6THrNbEzk5OTh69CiOHz+OV155BQCg0+kghICFhQX++OMPDBgwoNxxVlZWsLKyauhyTcaQju5wsVMhLbsA289mIKKDu9wlERERERmN2fTsqtVqnDp1CnFxcfrH5MmT0bZtW8TFxSE0NFTuEk2SlYUST3crHbqx7NBVmashIiIiMi6T7tnNzc3FxYsX9c8TEhIQFxcHZ2dn+Pr6Ijo6GsnJyfjhhx+gUCjQsWNHg+NdXV1hbW1dbjsZGtvDF//Zcwl7L9zA1Zt58HOxk7skIiIiIqMw6Z7do0ePIjg4GMHBwQCAqVOnIjg4GNOnTwcApKamIjGRE6vqytfFFg/fXXpsRQzbk4iIiMyHJITg7bPukZ2dDY1Gg6ysLKjVarnLaTB/nE7DCz/GwsnWEgejB8LaUil3SUREREQVqkleM+meXWo4A9q5wkNjjdv5xdjyV5rc5RAREREZBcMuAQAslAqM6VG6DBknqhEREZG5YNglvVHdfaBUSDh69TbOpjWdO8kRERGR+WLYJT03tTUGt3cDACw/xIlqRERE1Pgx7JKBcT39AADrjycjr7BE5mqIiIiI6oZhlwyEtXRBi2Z2yC0swS9xKXKXQ0RERFQnDLtkQKGQEBX6v4lqXJmOiIiIGjOGXSrnya7eUFkocCY1G8eTMuUuh4iIiKjWGHapHCc7FR7t7AGAE9WIiIiocWPYpQqVTVT77WQKMvOLZK6GiIiIqHYYdqlCwT6OCPRQo7BEh3Wx1+Quh4iIiKhWGHapQpIkYVzP0olqy2MSOVGNiIiIGiWGXarU8C5esLeyQMKNPBy4dFPucoiIiIhqjGGXKmVvZYERwZ4AgOUxV2WuhoiIiKjmGHapSlGhpRPV/jidjozsApmrISIiIqoZhl2qUqCHGiF+TijRCaw6kiR3OUREREQ1wrBLD1Q2UW3l4USUaHUyV0NERERUfQy79EBDO3rAydYSqVkF2HnuutzlEBEREVUbwy49kLWlEk938wHAiWpERETUuDDsUrWM7VE6lGH3+etIupUvczVERERE1cOwS9Xi38wOfQKaQYjSm0wQERERNQYMu1RtZcuQrT2ahMISrczVEBERET0Ywy5VW3igK9zUVriZV4Qtf6XJXQ4RERHRAzHsUrVZKBUY3b107O7yQxzKQERERKaPYZdqZEwPXygVEg5fuYXz6Tlyl0NERERUJYZdqhF3jTXCA10BAMsPcRkyIiIiMm0Mu1RjZRPVfj6WjPyiEpmrISIiIqocwy7VWO/WzeDnYoucwhJsjEuRuxwiIiKiSjHsUo0pFJL+JhPLYq5CCCFzRUREREQVY9ilWnm6mw9UFgr8lZyNk9ey5C6HiIiIqEIMu1QrznYqDOvkAQBYxolqREREZKIYdqnWokJLhzL8ejIFWfnFMldDREREVB7DLtVaiJ8T2rk7oKBYh5+OXZO7HCIiIqJyGHap1iRJQlTP0mXIlnOiGhEREZkghl2qkxFdPGGrUuLS9TwcunxL7nKIiIiIDDDsUp04WFtiRLAXgNJlyIiIiIhMCcMu1VnZRLWtf6UhI6dA5mqIiIiI/odhl+qsg6cGwb6OKNEJrD3KiWpERERkOhh2ySjGhZZOVFsRkwitjhPViIiIyDQw7JJRDOvsAY2NJZIz72D3+Qy5yyEiIiICwLBLRmJtqcTTId4AgGWHEmWuhoiIiKiUSYfdPXv2IDIyEp6enpAkCRs2bKhy/59//hmDBg1C8+bNoVarERYWhq1btzZMsYSxdyeq7TyXgaRb+TJXQ0RERGTiYTcvLw9BQUFYsGBBtfbfs2cPBg0ahM2bNyM2Nhb9+/dHZGQkjh8/Xs+VEgC0bG6PXq1dIASw6gh7d4mIiEh+kmgkt72SJAnr16/HiBEjanRchw4dMGrUKEyfPr1a+2dnZ0Oj0SArKwtqtboWlTZtv59KxUvLj6GZvQoH3hsIlYVJ/zxFREREjVBN8ppZJxGdToecnBw4OztXuk9hYSGys7MNHlR74e3d4OpghRu5Rdh6Ok3ucoiIiKiJM+uwO2/ePOTm5mLkyJGV7jNnzhxoNBr9w8fHpwErND+WSgVGdy9tw+W8oxoRERHJzGzD7ooVKzBr1iysWbMGrq6ule4XHR2NrKws/SMpKakBqzRPo3v4QiEBhy7fwsWMHLnLISIioibMLMPuqlWr8Nxzz2HNmjUIDw+vcl8rKyuo1WqDB9WNp6MNBrRzAwAsj+FENSIiIpKP2YXdlStXYtKkSVi5ciWGDRsmdzlN1riepcuQ/RR7DXeKtDJXQ0RERE2VSYfd3NxcxMXFIS4uDgCQkJCAuLg4JCaW9hZGR0djwoQJ+v1XrFiBCRMmYP78+QgNDUVaWhrS0tKQlZUlR/lN2sMBzeHjbIPsghL8eiJF7nKIiIioiTLpsHv06FEEBwcjODgYADB16lQEBwfrlxFLTU3VB18A+O9//4uSkhJMmTIFHh4e+sfrr78uS/1NmUIhYWwPPwCcqEZERETyaTTr7DYUrrNrPDdzC9FzznYUawV+faU3Onlr5C6JiIiIzADX2SWT4GJvhaEdPQCwd5eIiIjkwbBL9Wpcz9KhDL/EpSC7oFjmaoiIiKipYdiletXd3wlt3Oxxp1iLn2OvyV0OERERNTEMu1SvJElCVGjZRLVEcIg4ERERNSSGXap3j3f1go2lEhcycnE44Zbc5RAREVETwrBL9U5tbYnhXTwBAMt4RzUiIiJqQAy71CDKJqpt+SsVN3ILZa6GiIiImgqGXWoQHb00CPJxRLFWYM3RJLnLISIioiaCYZcaTFSoLwBgRUwidDpOVCMiIqL6x7BLDSaysyfU1ha4dvsOdl+4Lnc5RERE1AQw7FKDsVEp8WSINwBg+SHeUY2IiIjqH8MuNaiyNXd3nM1AcuYdmashIiIic8ewSw2qtas9wlq6QCeAVYe5DBkRERHVL4ZdanBRPUsnqq06koRirU7maoiIiMicMexSgxvc3h3N7K1wPacQ286ky10OERERmTGGXWpwKgsFRnf3AQAs40Q1IiIiqkcMuySL0T18IEnAgUs3cel6rtzlEBERkZli2CVZeDvZYkBbVwClN5kgIiIiqg8MuySbsolq62KvoaBYK3M1REREZI4Ydkk2fdu4wsvRBll3ivHbyVS5yyEiIiIzxLBLslEqJIwNLe3d5UQ1IiIiqg8MuySrkd18YKmUEJeUib+Ss+Quh4iIiMwMwy7JqrmDFSI6uAMAlnOiGhERERkZwy7JLirUDwDwS1wycgqKZa6GiIiIzAnDLsmuZ0tntHa1R36RFhuOJ8tdDhEREZkRhl2SnSRJiNJPVEuEEELmioiIiMhcMOySSXiiqzesLRU4l56Do1dvy10OERERmQmGXTIJGhtLPBbkCQBYzmXIiIiIyEgYdslkjOtZOlFt86k03MwtlLkaIiIiMgcMu2QyOns7opOXBkVaHdbFXpO7HCIiIjIDDLtkUsb1LJ2otuJwInQ6TlQjIiKiumHYJZMSGeQJB2sLXL2Zj70Xb8hdDhERETVyDLtkUmxVFniyqzcATlQjIiKiumPYJZNTtubun/HpSM26I3M1RERE1JgZLex6e3sjMzPTWKejJizAzQE9WjhDJ4CVh5PkLoeIiIgaMaOF3ZSUFBQVFemfR0VF4datW8Y6PTUxZcuQrTqciGKtTuZqiIiIqLGqt2EMGzduZE8v1dqQDu5wsVMhI6cQ2+PT5S6HiIiIGimO2SWTpLJQYGR3HwDA8phEmashIiKixsqoYfe3337D5cuXjXlKasLG9vCFJAF7L9xAwo08ucshIiKiRshoYTc4OBgvv/wyAgIC4OTkhDt37uCjjz7C8uXLcfr0aeh0HHdJNePjbIu+bZoDAFbEcBkyIiIiqjkLY50oNjYWJSUlOH36NI4dO4bjx4/j2LFjWLFiBfLz82FlZYVOnTohJibGWJekJmBcqB92nbuOtbHX8ObgtrC2VMpdEhERETUiRgu7AGBhYYGgoCAEBQVh0qRJAAAhBM6fP4/Y2FjExcUZ83LUBPRv5wpPjTVSsgqw+VQqnrh7wwkiIiKi6qj3CWqSJKFt27YYO3YsPv744xodu2fPHkRGRsLT0xOSJGHDhg0PPGbXrl3o2rUrrKys0Lp1ayxZsqR2hZNJUCokjOlRepMJTlQjIiKimjLp1Rjy8vIQFBSEBQsWVGv/hIQEDBs2DP3790dcXBz+/ve/47nnnsPWrVvruVKqT6O6+8BCISH26m3Ep2bLXQ4RERE1IjUaxjB79uw6X1CSJEybNq1a+w4dOhRDhw6t9rm//fZbtGjRAvPnzwcABAYGYt++ffjss88QERFRq3pJfq5qawzu4IbNp9Kw7NBVfPB4J7lLIiIiokaiRmF35syZdb5gTcJuTR08eBDh4eEG2yIiIvD3v/+90mMKCwtRWFiof56dzZ5DUzQu1A+bT6Vhw/FkRD8SCHsrow43JyIiIjNVo2EMOp2uzg+tVltf7wVpaWlwc3Mz2Obm5obs7GzcuXOnwmPmzJkDjUajf/j4+NRbfVR7Ya1c0LKZHfKKtNhwPFnucoiIiKiRqFH32IABA+p8QUmSsH379jqfx1iio6MxdepU/fPs7GwGXhMkSRLGhvri/zbFY9mhq4gK9YUkSXKXRURERCauxj27Qog6Perz5hLu7u5IT0832Jaeng61Wg0bG5sKj7GysoJarTZ4kGl6KsQbVhYKnE3LwbHETLnLISIiokagRj27u3btqqcyjCMsLAybN2822LZt2zaEhYXJVBEZk6OtCpFBnlgXew3LD11FiJ+T3CURERGRiTPppcdyc3MRFxenvxlFQkIC4uLikJhYut5qdHQ0JkyYoN9/8uTJuHz5Mt555x2cPXsWX3/9NdasWYM33nhDjvKpHkSFlq65+9upVNzOK5K5GiIiIjJ1dQq7R48exYQJExAcHIxu3bph1KhRWLhwIbKysoxS3NGjRxEcHIzg4GAAwNSpUxEcHIzp06cDAFJTU/XBFwBatGiBTZs2Ydu2bQgKCsL8+fPx3XffcdkxM9LFxxEdPNUoKtFhXew1ucshIiIiEycJIURtDty4cSOeeuoplJSUlHvN3t4es2fPrnLJL1OVnZ0NjUaDrKwsjt81UStiEvGP9afg72KLHW/2g0LBiWpERERNSU3yWq17dqdPnw4hBD777DOkpKQgMzMTf/31F+bOnQtXV1e8+eabmDx5cm1PT1Sp4V08YW9lgSs383Hg0k25yyEiIiITVuuwGx8fj9GjR+P111+Hu7s71Go12rdvj3feeQenTp3ChAkTsHDhQqxYscKY9RLBzsoCT3T1AgAsO3RV5mqIiIjIlNU67FpYWMDX17fC12xsbPD999+jXbt2+OKLL2pdHFFlokL9AADb4tORllUgczVERERkqmoddtu0aYODBw9WfmKFAsOGDcNff/1V20sQVaqtuwO6+ztBqxNYfSRJ7nKIiIjIRNU67P7zn//Erl278P3331e6z507d+Dg4FDbSxBVaVzP0t7dlYcTUaKtv5uVEBERUeNV67D71FNPYeLEiXjhhRcwduxYHD582OD1uLg4rFixAk8++WSdiySqyJCO7nC2UyEtuwA7zmbIXQ4RERGZoDqts/v999/j7bffxk8//YSwsDA0b94cXbp0QZs2bRASEoJevXph/vz5xqqVyICVhRJPd/MGACyLSXzA3kRERNQU1Xqd3XudP38eP/74IzZt2oSTJ09Cpyv9lbIkSXBxcUGXLl0MHu3bt69z4fWF6+w2Lldv5qHvJ7sAALvf7gc/Fzt5CyIiIqJ6V5O8ZpSwe687d+4gLi4Ox44d0z/OnDmD4uLi0gtKErRarTEvaVQMu43PhEWHsef8dbzYtyWihwbKXQ4RERHVs5rkNQtjX9zGxgZhYWEICwvTbysqKsKpU6cQGxuLuLg4Y1+Smrhxob7Yc/461h69hqmD2sDKQil3SURERGQijB52K6JSqRASEoKQkJCGuBw1MQPaucJDY43UrAJs+SsNw7t4yV0SERERmYg6TVAjMgUWSgVGdy+9wQnvqEZERET3YtglszCquw+UCglHrtzG2bRsucshIiIiE8GwS2bBXWONQYFuAIAVXIaMiIiI7mLYJbNRdke1n48lI6+wROZqiIiIyBQYLexeunQJe/bsMdbpiGrsoVYu8HexRW5hCX6JS5G7HCIiIjIBRgu7n3/+Ofr372+s0xHVmEIhISq0tHd3ecxVGHkJaSIiImqEOIyBzMpTId5QWShwOiUbcUmZcpdDREREMmPYJbPiZKfCo508AADLOVGNiIioyWPYJbMTdXei2q8nUpCZXyRzNURERCQnhl0yO119HdHO3QGFJTqsi70mdzlEREQkI4ZdMjuSJOmXIVsRk8iJakRERE0Ywy6ZpRHBXrBTKXH5Rh4OXropdzlEREQkE4ZdMkv2VhYYEewFAFgWc1XmaoiIiEguDLtktsrW3P3jdDoysgtkroaIiIjkwLBLZqu9pxpdfR1RohNYfSRJ7nKIiIhIBgy7ZNbKJqqtPJwIrY4T1YiIiJoao4VdlUoFW1tbY52OyCge6eQBR1tLpGQVYOfZDLnLISIiogZmtLA7f/585OTkGOt0REZhbanE0yHeADhRjYiIqCniMAYye2PvTlTbff46km7ly1wNERERNSSGXTJ7LZrZoU9AMwgBrDicKHc5RERE1IAYdqlJiAr1BQCsOZKEohKdzNUQERFRQ2HYpSZhYKAb3NRWuJlXhC2n0+Quh4iIiBqI0cNufj7HRJLpsVQqMKp7ae/uskOcqEZERNRUGD3sTpgwwdinJDKKMT18oFRIOJxwCxfSuXIIERFRU2D0sJuamorZs2eX237nzh2MGTPG2JcjqjYPjQ0GtnMFACyP4UQ1IiKipsDoYXfdunX473//iw0bNui3JScno3fv3rh6lb8+JnlF3b2j2k+x15BfVCJzNURERFTfjB52PTw8sGbNGjz33HM4ffo0Dh48iJCQEHTu3Bm7du0y9uWIaqRP62bwc7FFTmEJfj2RInc5REREVM+MEnb/9re/4d///jd2796NzMxMPPTQQ5gzZw4GDx6MQYMG4d1338XixYuhUqmMcTmiWlMoJIztUTZRjUMZiIiIzJ1Rwq61tTVWr16NyMhIuLi4wM/PD7/++isKCgowbtw4PProo8a4DJFRPBXiDZVSgVPJWTh5LVPucoiIiKgeGSXsfv3119i/fz+ys7Nx/vx5fPbZZwgODkafPn2wdetWtG3bFvb29ggNDTXG5YjqxMXeCo90cgfAZciIiIjMnYWxT9iqVSu0atUKTzzxhH5bVlYWTpw4gZMnTxr7ckS1EtXTDxviUrDxRAr++Uh7aGwt5S6JiIiI6kGNenZ///33Wl1Eo9Hg4YcfxiuvvFKr4xcsWAB/f39YW1sjNDQUhw8frnL/zz//HG3btoWNjQ18fHzwxhtvoKCgoFbXJvPUzc8Jbd0cUFCsw8/Hr8ldDhEREdWTGoXd4cOH47vvvquvWiq0evVqTJ06FTNmzMCxY8cQFBSEiIgIZGRkVLj/ihUr8N5772HGjBmIj4/H999/j9WrV+Mf//hHg9ZNpk2SJIzrWTpRbXlMIoQQMldERERE9aFGYdfLywsvvvgi3n///Wofc/78+RoXda9PP/0Uzz//PCZNmoT27dvj22+/ha2tLRYtWlTh/gcOHECvXr0wduxY+Pv7Y/DgwRgzZswDe4Op6RkR7AVblRIXM3Jx6PItucshIiKielCjsBsTE4Pg4GDMmTMHEyZMQElJ5Yvynzp1CqNHj0aHDh1qXVxRURFiY2MRHh7+v4IVCoSHh+PgwYMVHvPQQw8hNjZWH24vX76MzZs345FHHqlw/8LCQmRnZxs8qGlwsLbE8C5eAIDlMZyoRkREZI5qFHZdXV2xZ88ePPLII1i2bBmGDBlSLhweOXIEw4cPR5cuXbBmzRp07dq11sXduHEDWq0Wbm5uBtvd3NyQlpZW4TFjx47F7Nmz0bt3b1haWqJVq1bo169fpcMY5syZA41Go3/4+PjUul5qfKJCS4cybD2dhus5hTJXQ0RERMZW46XHbG1t8csvv2Dy5MnYsWMHevfujaSkJOzZswcRERHo2bMnfv31V/Tq1QtbtmxBTExMfdRdqV27duHDDz/E119/jWPHjuHnn3/Gpk2b8K9//avC/aOjo5GVlaV/JCUlNWi9JK+OXhp08XFEsVZgzVF+7YmIiMxNrZYeUygU+Prrr9GiRQu8++67aNeuHQoKCiCEwMCBAzFt2jQ8/PDDdS6uWbNmUCqVSE9PN9ienp4Od3f3Co+ZNm0axo8fj+eeew4A0KlTJ+Tl5eGFF17AP//5TygUhvneysoKVlZWda6VGq9xPf0Ql5SJFTGJmNy3FZQKSe6SiIiIyEhqfVOJjRs3Yt26dQCAO3fuAAA++OADbNu2zShBFwBUKhVCQkKwfft2/TadToft27cjLCyswmPy8/PLBVqlUgkAnHFPFXq0swc0NpZIzryD3ecrXuWDiIiIGqcah93Vq1cjKCgIjz/+OGJjYzFy5Ej88MMP0Gg0mDVrFpYtW2bUAqdOnYqFCxdi6dKliI+Px0svvYS8vDxMmjQJADBhwgRER0fr94+MjMQ333yDVatWISEhAdu2bcO0adMQGRmpD71E97K2VOKpEG8AwPJDiTJXQ0RERMZUo2EM7dq1w4ULF6BUKjF+/Hj84x//QJs2bQAAXbt2xdChQzFx4kQkJSUZBNC6GDVqFK5fv47p06cjLS0NXbp0wZYtW/ST1hITEw16ct9//31IkoT3338fycnJaN68OSIjI/HBBx8YpR4yT2NDffH9vgTsOJeBa7fz4e1kK3dJREREZASSqMHv9q2srPDMM8/gvffeQ4sWLcq9npKSgkceeQSnTp3CCy+8gK+//hqS1LjGP2ZnZ0Oj0SArKwtqtVrucqgBjV14CAcu3cQr/VvjrYi2cpdDRERElahJXqvRMIbLly/jP//5T4VBFwA8PT2xb98+DBw4EP/5z38wfPhw/XheIlM3rqcfAGDVkSQUlehkroaIiIiMocZ3UHsQe3t7bN68GePHj8dvv/2Gfv361bY2ogY1qL0bmjtY4UZuIf44U/E6zkRERNS41Ho1hqpYWFhg6dKliI6OxtGjR+vjEkRGZ6lUYHT30puKcKIaERGReaiXsFvmgw8+wLfffluflyAyqtE9fKGQgIOXb+JiRq7c5RAREVEd1Wg1htmzZ9fqIvceJ0kSpk2bVqvzENU3L0cbDGjnij/jM7A85ipmRHaQuyQiIiKqgxqtxnD/zRpqdUFJglarrfN56gtXY6Cd5zIwafERqK0tEPOPcNiouD4zERGRKam31Rh0Ol2dH6YcdIkAoG9Ac3g72SC7oAS/nkyRuxwiIiKqgxoNYxgwYECdLyhJksHtf4lMjUIhYWyoLz7ecg7LYxIxspuP3CURERFRLdW4Z1cIUaeHTsf1S8n0jezmA0ulhBNJmfgrOUvucoiIiKiWatSzu2vXrnoqg8i0NLO3wpCOHvj1RAqWHbqKuU92lrskIiIiqoV6XXqMqDEbF+oLAPglLgXZBcUyV0NERES1wbBLVIkeLZwR4GqPO8VarD+WLHc5REREVAsMu0SVkCQJUXd7d5cduooarNJHREREJoJhl6gKT4R4w8ZSiQsZuThy5bbc5RAREVENMewSVUFtbYnHgjwBlPbuEhERUePCsEv0AON6+gEAfv8rFTdyC2WuhoiIiGqCYZfoATp5axDkrUGxVmDt0Wtyl0NEREQ1wLBLVA1RoaW9uysOX4VOx4lqREREjQXDLlE1RAZ5Qm1tgaRbd7DnwnW5yyEiIqJqYtglqgYblRJPhngDAJYdSpS5GiIiIqouhl2iaipbc3fH2XSkZN6RuRoiIiKqDoZdompq7eqAni2doRPAqsPs3SUiImoMGHaJaqBsotqqI0ko1upkroaIiIgehGGXqAYiOrijmb0VMnIK8eeZdLnLISIiogdg2CWqAZWFAqO6352oFsM7qhEREZk6hl2iGhrd3ReSBOy/eBOXr+fKXQ4RERFVgWGXqIZ8nG3Rv60rAGBFDCeqERERmTKGXaJaKFuGbG3sNRQUa2WuhoiIiCrDsEtUC/3ausLL0QZZd4qx6WSq3OUQERFRJRh2iWpBqZAw9m7vLieqERERmS6GXaJaerqbNywUEo4nZuJ0Spbc5RAREVEFGHaJasnVwRoRHd0BAMs5UY2IiMgkMewS1cG4u3dU23A8GTkFxTJXQ0RERPdj2CWqg54tndGquR3yi7TYcDxZ7nKIiIjoPgy7RHUgSRKi7vbuLo9JhBBC5oqIiIjoXgy7RHX0ZFdvWFsqcDYtB7FXb8tdDhEREd2DYZeojjS2lojs7AmAE9WIiIhMDcMukRGM61k6lGHTyVTcyiuSuRoiIiIqw7BLZASdvTXo6KVGkVaHtUeT5C6HiIiI7mLYJTICSZL0y5CtOJwInY4T1YiIiEwBwy6RkTzWxRMOVha4ejMf+y7ekLscIiIiQiMJuwsWLIC/vz+sra0RGhqKw4cPV7l/ZmYmpkyZAg8PD1hZWaFNmzbYvHlzA1VLTZWtygJPdPUCACyPuSpzNURERAQ0grC7evVqTJ06FTNmzMCxY8cQFBSEiIgIZGRkVLh/UVERBg0ahCtXrmDdunU4d+4cFi5cCC8vrwaunJqiqLsT1f6Mz0BaVoHM1RAREZHJh91PP/0Uzz//PCZNmoT27dvj22+/ha2tLRYtWlTh/osWLcKtW7ewYcMG9OrVC/7+/ujbty+CgoIauHJqitq4OaCHvzO0OoGVh7kMGRERkdxMOuwWFRUhNjYW4eHh+m0KhQLh4eE4ePBghcds3LgRYWFhmDJlCtzc3NCxY0d8+OGH0Gq1Fe5fWFiI7OxsgwdRXUT19AUArDqSiBKtTuZqiIiImjaTDrs3btyAVquFm5ubwXY3NzekpaVVeMzly5exbt06aLVabN68GdOmTcP8+fPxf//3fxXuP2fOHGg0Gv3Dx8fH6O+DmpYhHd3hYqdCenYh/oyveLgNERERNQyTDru1odPp4Orqiv/+978ICQnBqFGj8M9//hPffvtthftHR0cjKytL/0hK4hqpVDdWFko83a30hyZOVCMiIpKXSYfdZs2aQalUIj093WB7eno63N3dKzzGw8MDbdq0gVKp1G8LDAxEWloaiorK39nKysoKarXa4EFUV2N7+EKSgL0XbuDKjTy5yyEiImqyTDrsqlQqhISEYPv27fptOp0O27dvR1hYWIXH9OrVCxcvXoRO97+xkufPn4eHhwdUKlW910wEAL4utng4oDmA0ptMEBERkTxMOuwCwNSpU7Fw4UIsXboU8fHxeOmll5CXl4dJkyYBACZMmIDo6Gj9/i+99BJu3bqF119/HefPn8emTZvw4YcfYsqUKXK9BWqixt1dhmzt0SQUFFc8QZKIiIjql4XcBTzIqFGjcP36dUyfPh1paWno0qULtmzZop+0lpiYCIXif5ndx8cHW7duxRtvvIHOnTvDy8sLr7/+Ot5991253gI1UQPaucJTY42UrAL8/lcqHg/2lrskIiKiJkcSQgi5izAl2dnZ0Gg0yMrK4vhdqrN/b7+AT7edR4ifE3566SG5yyEiIjILNclrJj+MgagxG93dBxYKCbFXbyM+lWs4ExERNTSGXaJ65Kq2xuAOpUNuuAwZERFRw2PYJapnUaGlE9XWH0tGbmGJzNUQERE1LQy7RPXsoVYuaNnMDnlFWvwSlyx3OURERE0Kwy5RPZMkCWNDfQEAyw4lgnNCiYiIGg7DLlEDeCrEG1YWCsSnZuN4Uqbc5RARETUZDLtEDcDRVoVHO3sCAJYd4kQ1IiKihsKwS9RAonqWDmX47WQqbucVyVwNERFR08CwS9RAgn0c0d5DjaISHX46dk3ucoiIiJoEhl2iBiJJEsb1LF2GbHlMInQ6TlQjIiKqbwy7RA1oeBdP2FtZIOFGHg5evil3OURERGaPYZeoAdlZWeDxYC8AnKhGRETUEBh2iRpY2US1P86kIz27QOZqiIiIzBvDLlEDa+euRjc/J2h1AquPJMldDhERkVlj2CWSQdlEtZWHE1Gi1clcDRERkfli2CWSwZCO7nCytURqVgF2nM2QuxwiIiKzxbBLJANrSyVGdvMBULoMGREREdUPhl0imYzpUTpRbc+F60i8mS9zNUREROaJYZdIJv7N7NAnoBmEAFYcZu8uERFRfWDYJZJR2US1NUeTUFiilbkaIiIi88OwSySjge1c4a62xq28Imz5K03ucoiIiMwOwy6RjCyUCozucXei2iEOZSAiIjI2hl0imY3u7gulQsLhK7dwLi1H7nKIiIjMCsMukczcNdYID3QFAKyIuSpzNUREROaFYZfIBJRNVPv5WDLyCktkroaIiMh8MOwSmYBerZrBz8UWOYUl2HgiRe5yiIiIzAbDLpEJUCgkRIWW3mRi2aGrEELIXBEREZF5YNglMhFPhfhAZaHA6ZRsnLiWJXc5REREZoFhl8hEONupMKyTB4DS3l0iIiKqO4ZdIhMyrmfpUIZfT6QgK79Y5mqIiIgaP4ZdIhPS1dcJ7dwdUFiiw7pj1+Quh4iIqNFj2CUyIZIkIeruMmTLYzhRjYiIqK4YdolMzOPBXrBTKXH5eh4OXr4pdzlERESNGsMukYmxt7LA8GAvAMDyQ4kyV0NERNS4MewSmaBxoaVDGbaeTkNGToHM1RARETVeDLtEJqi9pxpdfR1RohNYcyRJ7nKIiIgaLYZdIhMVdbd3d+XhJGh1nKhGRERUGwy7RCZqWGcPONpaIjnzDnady5C7HCIiokaJYZfIRFlbKvFUV28AvKMaERFRbTHsEpmwsjV3d52/jqRb+TJXQ0RE1Pgw7BKZsBbN7NC7dTMIAaw8zGXIiIiIaophl8jERYX6AgDWHE1CUYlO5mqIiIgal0YRdhcsWAB/f39YW1sjNDQUhw8frtZxq1atgiRJGDFiRP0WSFSPwtu7wdXBCjdyi7D1dJrc5RARETUqJh92V69ejalTp2LGjBk4duwYgoKCEBERgYyMqmenX7lyBW+99Rb69OnTQJUS1Q9LpQKje5T27nKiGhERUc2YfNj99NNP8fzzz2PSpElo3749vv32W9ja2mLRokWVHqPVahEVFYVZs2ahZcuWDVgtUf0Y3d0HCgmISbiFixk5cpdDRETUaJh02C0qKkJsbCzCw8P12xQKBcLDw3Hw4MFKj5s9ezZcXV3x7LPPPvAahYWFyM7ONngQmRpPRxsMDHQDACw7xIlqRERE1WXSYffGjRvQarVwc3Mz2O7m5oa0tIrHLu7btw/ff/89Fi5cWK1rzJkzBxqNRv/w8fGpc91E9aFsotpPx64hv6hE5mqIiIgaB5MOuzWVk5OD8ePHY+HChWjWrFm1jomOjkZWVpb+kZSUVM9VEtXOwwHN4etsi5yCEvx2IlXucoiIiBoFC7kLqEqzZs2gVCqRnp5usD09PR3u7u7l9r906RKuXLmCyMhI/TadrnSpJgsLC5w7dw6tWrUyOMbKygpWVlb1UD2RcSkUEsaG+mLu72exLOYqRnbnbyGIiIgexKR7dlUqFUJCQrB9+3b9Np1Oh+3btyMsLKzc/u3atcOpU6cQFxenfzz22GPo378/4uLiOESBGr2nQ7yhUipw8loWTl7LlLscIiIik2fSPbsAMHXqVEycOBHdunVDjx498PnnnyMvLw+TJk0CAEyYMAFeXl6YM2cOrK2t0bFjR4PjHR0dAaDcdqLGyMXeCkM7ueOXuBQsP5SIzk85yl0SERGRSTP5sDtq1Chcv34d06dPR1paGrp06YItW7boJ60lJiZCoTDpDmoio4oK9cMvcSn45UQy/jEsEBobS7lLIiIiMlmSEELIXYQpyc7OhkajQVZWFtRqtdzlEJUjhEDE53twPj0XMyPb45leLeQuiYiIqEHVJK+xS5SokZEkCeN6+gEAlsUkgj+vEhERVY5hl6gRGhHsBRtLJS5m5CIm4Zbc5RAREZkshl2iRkhtbYkRwZ4AgOUxvKMaERFRZRh2iRqpqNDSoQxb/krF9ZxCmashIiIyTQy7RI1URy8NgnwcUawVWHOUd/4jIiKqCMMuUSM2LtQXALDycCK0Ok5UIyIiuh/DLlEj9mhnT6itLXDt9h3sOX9d7nKIiIhMDsMuUSNmo1LiqZDS22Avj7kqczVERESmh2GXqJGL6lk6lGHH2QwkZ96RuRoiIiLTwrBL1Mi1am6PsJYu0AlgJZchIyIiMsCwS2QGyu6otupIEoq1OpmrISIiMh0Mu0RmYHAHNzR3sMKN3EL8cTpd7nKIiIhMBsMukRmwVCowqhsnqhEREd2PYZfITIwJ9YVCAg5cuolL13PlLoeIiMgkMOwSmQkvRxv0b+sKAFh+iBPViIiIAIZdIrNSNlFtXWwSCoq1MldDREQkP4ZdIjPycJvm8HayQXZBCX49kSJ3OURERLJj2CUyI0qFhDE9Sm8ysYxr7hIRETHsEpmbkd18YKmUcCIpE38lZ8ldDhERkawYdonMTHMHK0R0cAfAZciIiIgYdonMUNlEtV/iUpBdUCxzNURERPJh2CUyQ6EtnNHa1R75RVpsOJ4sdzlERESyYdglMkOSJCEq9O5EtUNXIYSQuSIiIiJ5MOwSmaknunrDxlKJ8+m5OHr1ttzlEBERyYJhl8hMaWws8ViQJ4DS3l0iIqKmiGGXyIxF9SwdyvD7qTTczC2UuRoiIqKGx7BLZMY6ezuis7cGRVod1sZek7scIiKiBsewS2TmyiaqrYhJhE7HiWpERNS0MOwSmbnIIE84WFsg8VY+9l68IXc5REREDYphl8jM2aos8GRXbwCcqEZERE0Pwy5RE1A2lGF7fDpSMu/IXA0REVHDYdglagIC3BwQ2sIZOgGsOpIkdzlEREQNhmGXqImI6ukHAFh1OBHFWp3M1RARETUMhl2iJmJIB3c0s1chI6cQ2+PT5S6HiIioQTDsEjURKgsFRnbzAQAsO5QoczVEREQNg2GXqAkZ08MXkgTsu3gDCTfy5C6HiIio3jHsEjUhPs626NemOQBgRQyXISMiIvPHsEvUxIy7O1Ftbew1FBRrZa6GiIiofjHsEjUx/dq6wsvRBpn5xdh0MlXucoiIiOoVwy5RE6NUSBjTo3Si2nIOZSAiIjPHsEvUBI3s7gMLhYRjiZk4k5ItdzlERET1hmGXqAlydbBGRAd3AOzdJSIi89Yowu6CBQvg7+8Pa2trhIaG4vDhw5Xuu3DhQvTp0wdOTk5wcnJCeHh4lfsTNVVRPX0BABuOJyO3sETmaoiIiOqHyYfd1atXY+rUqZgxYwaOHTuGoKAgREREICMjo8L9d+3ahTFjxmDnzp04ePAgfHx8MHjwYCQnJzdw5USmLaylC1o2t0NekRbrj/PvBxERmSdJCCHkLqIqoaGh6N69O7766isAgE6ng4+PD1599VW89957Dzxeq9XCyckJX331FSZMmFDu9cLCQhQWFuqfZ2dnw8fHB1lZWVCr1cZ7I0Qm6Pt9CfjXb2fQzt0Bv7/eB5IkyV0SERHRA2VnZ0Oj0VQrr5l0z25RURFiY2MRHh6u36ZQKBAeHo6DBw9W6xz5+fkoLi6Gs7Nzha/PmTMHGo1G//Dx8TFK7USNwVNdvWFlocDZtBwcS7wtdzlERERGZ9Jh98aNG9BqtXBzczPY7ubmhrS0tGqd491334Wnp6dBYL5XdHQ0srKy9I+kpKQ6103UWGhsLREZ5AkAWHYoUeZqiIiIjM+kw25dzZ07F6tWrcL69ethbW1d4T5WVlZQq9UGD6KmpOyOaptOpeJWXpHM1RARERmXSYfdZs2aQalUIj093WB7eno63N3dqzx23rx5mDt3Lv744w907ty5PsskatSCvDXo4KlGUYkO62L5mw0iIjIvJh12VSoVQkJCsH37dv02nU6H7du3IywsrNLjPv74Y/zrX//Cli1b0K1bt4YolajRkiRJ37u7IiYROp1Jz1klIiKqEZMOuwAwdepULFy4EEuXLkV8fDxeeukl5OXlYdKkSQCACRMmIDo6Wr//Rx99hGnTpmHRokXw9/dHWloa0tLSkJubK9dbIDJ5jwV5wsHKAldu5mP/pRtyl0NERGQ0Jh92R40ahXnz5mH69Ono0qUL4uLisGXLFv2ktcTERKSmpur3/+abb1BUVISnnnoKHh4e+se8efPkegtEJs/OygKPd/UCACw7xDuqERGR+TD5dXYbWk3WbSMyJ+fSchDx+R4oFRL2vzsA7pqKJ3USERHJzWzW2SWihtPW3QHd/Z2g1QmsOsJlyIiIyDww7BKRXtlEtVWHk1Ci1clcDRERUd0x7BKR3pCO7nC2UyEtuwDbz2bIXQ4REVGdMewSkZ6VhRJPd/MGwIlqRERkHhh2ichAVA8/SBKw98INXL2ZJ3c5REREdcKwS0QGfF1s8XBAcwClN5kgIiJqzBh2iaicqFBfAMCao0koKNbKXA0REVHtMewSUTkD2rnCQ2ON2/nF2PJXmtzlEBER1RrDLhGVY6FUYHT30t5dTlQjIqLGjGGXiCo0uocPlAoJR6/extm0bLnLISIiqhWGXSKqkJvaGoPbuwEAlh/iRDUiImqcGHaJqFJRoaV3VFt/PBl5hSUyV0NERFRzDLtEVKmHWrmgRTM75BaW4MdDV/HZtvPIyC6QuywiIqJqY9glokopFJJ+GbLVR5LwxfYLyMgplLkqIiKi6mPYJaIqPdnVGyoLBRJu8G5qRETU+FjIXQARma6M7AJk5BSiV2sX7Dx7HQAwduEh2KiUsFBIsLZQwsZKCUulApZKBVRKBSyVUulzi/ue3/dnlcV9z5UKWFpUvO/9r6vu3efudSwUEpQKCZIkydxqBJR+dpbHJCIq1Beuamu5yyGiJoxhl4gqtTwmEV9sv2CwLbugBNkFpjlZTZJQPnRXEqwtDF6/P1j/L7CXPbe4/zV9oL/v2LvnslBUdN7/BXYLMw/mGTmF+GL7BQxq78aw24Txh56mw5S/1gy7RFSpqFBfDLq7/Nihyzfxf5vi8eqA1vB1sUWJVsDeygL2VhYo0upQrNWhRCv0fy4u0aH43ufau89L/vfcYH+tQLFWd8/rd5/rzydQoit7vfS1Ep0wqFcIoKikdJ/GQHVP77SF4m5wtqggdFcQrC2U0j093KUh2uC5Urob8ivqKb+7r8V9z/XbJFgqDP+sUJhvMKf6wx96mg5T/loz7BJRpVzV1uW+aUV0cEdHL41MFRnS6QSKdXfD792QXHRPUL43NBeXlH/t3uBcfP9rdwO2wfN7rlOsu+fPWh2K7j1PyX3P7x57vyKtDkVaoPQ/ps1CIT0wRAMCQpT2rt8pLn1P7/50AmprFQDA2lIBG5USACBBwt3/lT6XpHv+XLq9rOdb0v+n9DjJ4Li75yr78z07V7SfpD8PDHrWparOf18dUjXPj3uOK1dnBeev8P1UUNv/9r3//Pe05X11VNpe950f9xxX7mtSRR2VnT/pVj4A4MDFG7h2Ox9kvhJvln59M/OLZK6kPIZdImq0FAoJVgolrCwAWMldTdWEECjRCX0vtWGPtw5F9wXrEoPgXUmwLvnf85KKQvd95733XCW68j3tZT3xFfWMl+gESnRa3Cmu2fs+nZJjpBakxuzD38/KXQI1kF3nrsPRtvQHXFcHK5Po5WXYJaJqcXWwwusDA+DqYOKp0kRJkqTvBYVK7mqqJoSAVif04bekkmEl9wf269lFuJlXiBKdDldv5uOnY8kY3sUTno42gBCwt7aAg7UlhCi9htBfD/f8WdxTByDu9hYDpfv878//2152XFX7iXtOanhdYVDDvcfB4HyV71eujirOD/2x99RbxX6V1lFBu/3vPd53vgfUYfi+DM9//3UrrPe+tki6nY+UzMrX4/bUWMPLyabS16nxSL59BylZhl/r7/Yl4Lt9CQCA1wcG4I1BbeQozQDDLhFVi6va2iS+aVH9kyQJFkoJFkrABspaneOv5Cz8dCwZz/dpaTLDXqhhlK3iApR+Dt77+RTmPtFJ/zkwld4+qrvqfK1NAcMuERERGU1FY/07emn4Q48Zaixfa95UgoiIjI7DXojIVLBnl4iIjI7DXgjgDz1NiSl/rSVx76hzQnZ2NjQaDbKysqBWq+Uuh4iIiIjuU5O8xmEMRERERGS2GHaJiIiIyGwx7BIRERGR2WLYJSIiIiKzxbBLRERERGaLYZeIiIiIzBbDLhERERGZLYZdIiIiIjJbDLtEREREZLYYdomIiIjIbDHsEhEREZHZYtglIiIiIrPFsEtEREREZothl4iIiIjMloXcBZgaIQQAIDs7W+ZKiIiIiKgiZTmtLLdVhWH3Pjk5OQAAHx8fmSshIiIioqrk5ORAo9FUuY8kqhOJmxCdToeUlBQ4ODhAkqQq983OzoaPjw+SkpKgVqsbqELzwLarG7Zf3bD9ao9tVzdsv9pj29WNubWfEAI5OTnw9PSEQlH1qFz27N5HoVDA29u7Rseo1Wqz+ODIgW1XN2y/umH71R7brm7YfrXHtqsbc2q/B/XoluEENSIiIiIyWwy7RERERGS2GHbrwMrKCjNmzICVlZXcpTQ6bLu6YfvVDduv9th2dcP2qz22Xd005fbjBDUiIiIiMlvs2SUiIiIis8WwS0RERERmi2GXiIiIiMwWwy4RERERma0mFXYXLFgAf39/WFtbIzQ0FIcPH65y/7Vr16Jdu3awtrZGp06dsHnzZoPXhRCYPn06PDw8YGNjg/DwcFy4cEH/+q5duyBJUoWPI0eO6Pc7efIk+vTpA2tra/j4+ODjjz827hs3AlNsuytXrlT4+qFDh4zfAHXU0O0HAOfPn8fw4cPRrFkzqNVq9O7dGzt37jTYJzExEcOGDYOtrS1cXV3x9ttvo6SkxDhv2khMte0q+uytWrXKOG/aiORov2PHjmHQoEFwdHSEi4sLXnjhBeTm5hrs0xg+e4Dptl9j+PwZu+1+/vlnDB48GC4uLpAkCXFxceXOUVBQgClTpsDFxQX29vZ48sknkZ6ebrBPU/3sGav9GsNnrxzRRKxatUqoVCqxaNEicfr0afH8888LR0dHkZ6eXuH++/fvF0qlUnz88cfizJkz4v333xeWlpbi1KlT+n3mzp0rNBqN2LBhgzhx4oR47LHHRIsWLcSdO3eEEEIUFhaK1NRUg8dzzz0nWrRoIXQ6nRBCiKysLOHm5iaioqLEX3/9JVauXClsbGzEf/7zn/pvlGoy1bZLSEgQAMSff/5psF9RUVH9N0oNyNF+QggREBAgHnnkEXHixAlx/vx58fLLLwtbW1uRmpoqhBCipKREdOzYUYSHh4vjx4+LzZs3i2bNmono6Oj6bZAaMNW2E0IIAGLx4sUGn717z2EK5Gi/5ORk4eTkJCZPnizOnj0rDh8+LB566CHx5JNP6s/RGD57Qphu+wlh+p+/+mi7H374QcyaNUssXLhQABDHjx8vd57JkycLHx8fsX37dnH06FHRs2dP8dBDD+lfb8qfPWO0nxCm/9mrSJMJuz169BBTpkzRP9dqtcLT01PMmTOnwv1Hjhwphg0bZrAtNDRUvPjii0IIIXQ6nXB3dxeffPKJ/vXMzExhZWUlVq5cWeE5i4qKRPPmzcXs2bP1277++mvh5OQkCgsL9dveffdd0bZt25q/yXpiqm1XFnYr+gtrSuRov+vXrwsAYs+ePfp9srOzBQCxbds2IYQQmzdvFgqFQqSlpen3+eabb4RarTb4PMrJVNtOiNJv+OvXr6/ze6xPcrTff/7zH+Hq6iq0Wq1+n5MnTwoA4sKFC0KIxvHZE8J0208I0//8Gbvt7lXZ9/7MzExhaWkp1q5dq98WHx8vAIiDBw8KIZruZ+9edWk/IUz/s1eRJjGMoaioCLGxsQgPD9dvUygUCA8Px8GDBys85uDBgwb7A0BERIR+/4SEBKSlpRnso9FoEBoaWuk5N27ciJs3b2LSpEkG13n44YehUqkMrnPu3Dncvn275m/WyEy57co89thjcHV1Re/evbFx48Yav8f6JFf7ubi4oG3btvjhhx+Ql5eHkpIS/Oc//4GrqytCQkL01+nUqRPc3NwMrpOdnY3Tp08bpwHqwJTbrsyUKVPQrFkz9OjRA4sWLYIwoWXL5Wq/wsJCqFQqKBT/++fFxsYGALBv3z79dUz5sweYdvuVMdXPX320XXXExsaiuLjY4Dzt2rWDr6+v/jxN9bNXHdVpvzKm+tmrTJMIuzdu3IBWqzX4cAOAm5sb0tLSKjwmLS2tyv3L/r8m5/z+++8REREBb2/vB17n3mvIyZTbzt7eHvPnz8fatWuxadMm9O7dGyNGjDCpwCtX+0mShD///BPHjx+Hg4MDrK2t8emnn2LLli1wcnKq8jr3XkNOptx2ADB79mysWbMG27Ztw5NPPomXX34ZX375Zd3etBHJ1X4DBgxAWloaPvnkExQVFeH27dt47733AACpqalVXufea8jNlNsPMO3PX320XXWkpaVBpVLB0dGx0vM01c9edVSn/QDT/uxVxkLuApqKa9euYevWrVizZo3cpTQ6lbVds2bNMHXqVP3z7t27IyUlBZ988gkee+yxhi7TpAghMGXKFLi6umLv3r2wsbHBd999h8jISBw5cgQeHh5yl2iyqtt206ZN0x8THByMvLw8fPLJJ3jttdfkKt0kdOjQAUuXLsXUqVMRHR0NpVKJ1157DW5ubga9lVSx6rYfP38kl8b42WsS33maNWsGpVJZbkZheno63N3dKzzG3d29yv3L/r+651y8eDFcXFzKhbDKrnPvNeRkym1XkdDQUFy8ePGB+zUUudpvx44d+O2337Bq1Sr06tULXbt2xddffw0bGxssXbq0yuvcew05mXLbVSQ0NBTXrl1DYWFhzd5oPZHz7+7YsWORlpaG5ORk3Lx5EzNnzsT169fRsmXLKq9z7zXkZsrtVxFT+vzVR9tVh7u7O4qKipCZmVnpeZrqZ686qtN+FTGlz15lmkTYValUCAkJwfbt2/XbdDodtm/fjrCwsAqPCQsLM9gfALZt26bfv0WLFnB3dzfYJzs7GzExMeXOKYTA4sWLMWHCBFhaWpa7zp49e1BcXGxwnbZt2xr8ylQuptx2FYmLizOpXku52i8/Px8AyvWkKRQK6HQ6/XVOnTqFjIwMg+uo1Wq0b9++tm/ZaEy57SoSFxcHJycnWFlZ1eBd1h+5/+4Cpb/+tLe3x+rVq2FtbY1Bgwbpr2PKnz3AtNuvIqb0+auPtquOkJAQWFpaGpzn3LlzSExM1J+nqX72qqM67VcRU/rsVUrGyXENatWqVcLKykosWbJEnDlzRrzwwgvC0dFRPyNz/Pjx4r333tPvv3//fmFhYSHmzZsn4uPjxYwZMypcQsbR0VH88ssv4uTJk2L48OHlljASQog///xTABDx8fHl6srMzBRubm5i/Pjx4q+//hKrVq0Stra2Jrf0mCm23ZIlS8SKFStEfHy8iI+PFx988IFQKBRi0aJF9dQStSNH+12/fl24uLiIJ554QsTFxYlz586Jt956S1haWoq4uDghxP+W4Bk8eLCIi4sTW7ZsEc2bNzepJXhMte02btwoFi5cKE6dOiUuXLggvv76a2FrayumT5/egK3zYHL93f3yyy9FbGysOHfunPjqq6+EjY2N+OKLL/SvN4bPnhCm236N4fNXH2138+ZNcfz4cbFp0yYBQKxatUocP37cYEnAyZMnC19fX7Fjxw5x9OhRERYWJsLCwvSvN+XPnjHarzF89irSZMKuEKXfQHx9fYVKpRI9evQQhw4d0r/Wt29fMXHiRIP916xZI9q0aSNUKpXo0KGD2LRpk8HrOp1OTJs2Tbi5uQkrKysxcOBAce7cuXLXHTNmTLl16u514sQJ0bt3b2FlZSW8vLzE3Llz6/ZG64Eptt2SJUtEYGCgsLW1FWq1WvTo0cNgyRRTIkf7HTlyRAwePFg4OzsLBwcH0bNnT7F582aDfa5cuSKGDh0qbGxsRLNmzcSbb74piouLjfvm68gU2+73338XXbp0Efb29sLOzk4EBQWJb7/91mC5KFMhR/uNHz9eODs7C5VKJTp37ix++OGHcnU1hs+eEKbZfo3l82fstlu8eLEAUO4xY8YM/T537twRL7/8snBychK2trbi8ccfNwhzQjTdz54x2q+xfPbuJwlh4utFEBERERHVUpMYs0tERERETRPDLhERERGZLYZdIiIiIjJbDLtEREREZLYYdomIiIjIbDHsEhEREZHZYtglIiIiIrPFsEtEREREZothl4ioEXrmmWcgSRKuXLkidynVNmLECAQGBkKr1cpdioFx48bBz88PBQUFcpdCRPWAYZeImoyygHjvw8HBASEhIfj4449RWFgod4l6S5YsgSRJWLJkidylGMXu3bvxyy+/YMaMGVAqlXKXY2D69OlITk7G559/LncpRFQPGHaJqMl59tlnMWPGDEyfPh2jR4/G1atX8e6772L48OFyl1Ztc+bMQXx8PLy8vOQupVqmTZsGPz8/jBw5Uu5SymnTpg2GDx+OuXPnIi8vT+5yiMjIGHaJqMl57rnnMHPmTMyaNQsLFy7EuXPn4Onpia1bt2Lnzp1yl1ctHh4eaNeuHSwtLeUu5YFOnz6NvXv3Yty4cVAoTPOfnXHjxiErKwurVq2SuxQiMjLT/K5DRNSAXFxcMGLECABAbGysfvuuXbsgSRJmzpxZ7pgrV65AkiQ888wzBtv9/f3h7++P3NxcvP766/D09ISVlRU6d+6MdevWVaueZ555BpMmTQIATJo0yWDYxb373D9m9956Dxw4gP79+8PBwQHNmzfHyy+/jDt37gAANm3ahLCwMNjZ2cHNzQ3vvPMOSkpKKqzll19+wcCBA+Hk5ARra2t07NgR8+bNq9G428WLFwMAnn766XKvZWVlYfr06Wjfvj3s7e2hVqvRunVrTJw4EVevXjXYVwiBRYsWoVevXlCr1bC1tUW3bt2waNGiCq8rhMDixYvRp08fODo6wtbWFgEBAXjxxReRmJhosO+wYcNga2trNsNGiOh/LOQugIjIlFhY1P3bYnFxMQYPHozbt2/jySefRH5+PlatWoWRI0diy5YtGDx4cJXHjxgxApmZmfjll18wfPhwdOnSpUbXj4mJwUcffYSIiAi8+OKL2LlzJ7755htkZ2cjMjISzzzzDIYPH46wsDBs2rQJn3zyCezt7TF9+nSD80RHR2Pu3Lnw8vLCE088AY1Gg7179+Ltt99GTEwM1q5dW616tm/fDjs7O3Ts2NFguxACERERiImJQa9evTBkyBAoFApcvXoVGzduxPjx4+Hn56ffNyoqCitXrkRAQADGjh0LlUqFbdu24dlnn8WZM2cwb948/bl1Oh1GjRqFdevWwcvLC2PGjIFarcaVK1ewZs0aDB06FL6+vvr9VSoVQkJCcPDgQeTl5cHOzq5GbU5EJkwQETUREydOFADEwYMHDbbfuHFDeHp6CgDi8OHD+u07d+4UAMSMGTPKnSshIUEAEBMnTjTY7ufnJwCI4cOHi8LCQv32P//8UwAQERER1ap18eLFAoBYvHhxle8lISGhXL0AxIYNG/Tbi4qKROfOnYUkSaJZs2YG7zE7O1u4uroKZ2dnUVRUpN/+xx9/6OvNzc3Vb9fpdGLy5MkCgFi3bt0D30dOTo5QKBSiV69e5V47efKkACBGjBhR7rWCggKRk5Ojf/7f//5XABCTJk0yqLOwsFBERkYKAOLo0aP67V9++aUAIAYOHCjy8/MNzp2fny9u3rxZ7ppvvPGGACB27NjxwPdFRI0HhzEQUZPz3XffYebMmZgxYwaef/55tGvXDikpKXjttdfQvXt3o1zjs88+g0ql0j8fOHAg/Pz8cOTIEaOcvyr9+/c3mGxnaWmJp556CkIIREZGGrxHBwcHPProo7h16xauXbum3/7VV18BAP773/8a9HJKkoS5c+dCkiSsXLnygbWkpKRAp9PBzc2t0n1sbGzKbbOysoK9vb1BPXZ2dliwYIHBOGWVSoUPPvgAAAzq+frrr6FUKvHNN9+UO7+NjQ2cnZ3LXbOsxnvbgYgaPw5jIKIm5/vvvy+37c033zT4NXhdODo6okWLFuW2e3t74+DBg0a5RlUqGvbg4eHxwNdSUlL0dR86dAh2dnaVjoe1sbHB2bNnH1jLzZs3AZS2yf0CAwPRuXNnrFy5EteuXcOIESPQr18/dOnSxWAiW35+Pk6dOgVPT0989NFH5c5TXFwMAPp6cnNzER8fj9atWyMgIOCBNZYpC8A3btyo9jFEZPoYdomoyTl48CB69uyJoqIinDhxAi+//DLmz5+PwMBAPPvss3U+v0ajqXC7hYUFdDpdnc//IGq1usJrP+i1stAIALdu3UJJSQlmzZpV6XWqs0xXWa9qRTdssLCwwI4dOzBz5kz89NNPePPNNwEAzZs3xyuvvIJ//vOfUCqVuH37NoQQSE5OrlY9WVlZAFDjZdnKJvDZ2trW6DgiMm0cxkBETZZKpUL37t2xefNmODk54bXXXkNycrL+9bLexYpWKigLVOZKrVbDxcUFQohKHwkJCQ88T/PmzQGUhueKuLi44Msvv0RycjLOnDmDr776Cs7OzpgxYwY+/vhjfS0AEBISUmU9ZcvGlf2wce/XsjrKaiyrmYjMA8MuETV5zZs3x4wZM5Cfn2/Qc+jk5ASg4tB0/Pjxeq2p7C5jct1aNzQ0FDdv3sSFCxfqdB5PT0+4uLjg3LlzVe4nSRICAwMxZcoUbNu2DQCwceNGAKXjigMDAxEfH4/MzMwHXtPe3h7t27dHQkJCjeovq7FTp07VPoaITB/DLhERgBdffBGenp5YvHixvseybdu2cHBwwMaNGw16JtPT0/F///d/9VpP2fjRpKSker1OZV577TUAwN/+9jf9uNt7paWlIT4+/oHnkSQJffr0QUJCAq5fv27w2pUrVwzWCS6Tnp4OALC2tjaoJz8/H88//3yFwycSEhIMzjVlyhRotVqD9YXLFBQUVNjTHBMTAw8PjxqN8yUi08cxu0REKA1W7733Hl577TXMnj0bixcvhkqlwquvvooPP/wQXbt2xfDhw5GTk4Nff/0Vffv2xaVLl+qtnrCwMNjY2ODzzz/H7du39b9af//99+vtmvcaMmQIpk2bhn/9619o3bo1hgwZAj8/P9y8eRMXL17E3r178X//938IDAx84Lkef/xxbNiwAdu2bcPYsWP12+Pi4vDEE0+gR48eaN++Pdzd3ZGcnIwNGzZAoVDgjTfe0O/74osv4tChQ1i6dCn279+P8PBweHp6Ij09HWfPnkVMTAxWrFgBf39/AMBLL72E3bt3Y82aNQgICMBjjz0GtVqNxMREbN26Fd9//73+RiIAcOnSJSQkJOCll14yWhsSkYlo+NXOiIjkUdk6u2UKCgqEl5eXUCqV4ty5c0IIIbRarZg5c6bw8fERKpVKtGnTRnzxxRfi8uXLla6z6+fnV+H5+/btK2rybXfTpk2ie/fuwsbGRr9+7v3vpaJ1ditaF7iqdXtnzJghAIidO3eWe23btm0iMjJSNG/eXFhaWgp3d3cRFhYm/vWvf4nExMRqvY87d+4IZ2dnMXToUIPtSUlJ4r333hM9e/YUrq6uQqVSCV9fX/HEE09U+jVavXq1CA8PF05OTsLS0lJ4eXmJfv36ifnz54vr168b7KvT6cR3330nevbsKezs7IStra0ICAgQkydPLlf7zJkzBQARFxdXrfdERI2HJIQQcgVtIiJqGqZNm4a5c+fi4sWL+ruimYqSkhIEBASgRYsW2LFjh9zlEJGRccwuERHVu3feeQfOzs76G0CYkqVLl+Lq1atGW2eZiEwLwy4REdU7BwcH/Pjjj/D395dthYnKSJKEhQsXomvXrnKXQkT1gMMYiIiIiMhssWeXiIiIiMwWwy4RERERmS2GXSIiIiIyWwy7RERERGS2GHaJiIiIyGwx7BIRERGR2WLYJSIiIiKzxbBLRERERGaLYZeIiIiIzNb/AyeNK9DTgxRIAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "Plt_Err_Time(worker) # plot the error and time evolution of the RGD algorithm\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14023014", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "qiskit10", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/qibo/tomography_RGD/methods_ParmBasic.py b/src/qibo/tomography_RGD/methods_ParmBasic.py new file mode 100755 index 0000000000..b8b73641b0 --- /dev/null +++ b/src/qibo/tomography_RGD/methods_ParmBasic.py @@ -0,0 +1,76 @@ +# +# get basic parameter information for tomography +# + + +class BasicParameterInfo: + def __init__(self, params_dict): + + # + # read in params_dict + # + Nr = params_dict.get("Nr", 1) + trace = params_dict.get("trace", 1.0) + target_state = params_dict.get("target_state", None) + target_DM = params_dict.get("target_DM") + + label_list = params_dict.get("labels") + measurement_list = params_dict.get("measurement_list") + projector_list = params_dict.get("projector_list") + + num_iterations = params_dict.get("num_iterations") + convergence_check_period = params_dict.get("convergence_check_period", 10) + + relative_error_tolerance = params_dict.get("relative_error_tolerance", 0.001) + seed = params_dict.get("seed", 0) + + # + # basic system information + # + n = len(label_list[0]) + d = 2**n + + self.n = n # number of qubits + self.num_elements = d + self.Nr = Nr # the rank of the target_density_matrix + + self.trace = trace + self.target_state = target_state + self.target_DM = target_DM # target state --> density matrix + + self.num_labels = len(label_list) + self.measurement_list = measurement_list + self.projector_list = projector_list + + # + # numerical book keeping + # + + self.process_idx = 0 + + self.num_iterations = num_iterations + + self.convergence_check_period = ( + convergence_check_period # how often to check convergence + ) + self.relative_error_tolerance = ( + relative_error_tolerance # user-decided relative error + ) + self.seed = seed + + self.Err_relative_st = [] + + self.Target_Err_st = [] # Frobenius norm difference from State + self.Target_Err_Xk = [] # Frobenius norm difference from matrix Xk + + self.target_error_list = [] + self.target_relative_error_list = [] + + self.fidelity_Xk_list = [] # Fidelity between (Xk, target) + self.Err_relative_Xk = [] + + self.iteration = 0 + self.converged = False + self.convergence_iteration = 0 + + self.fidelity_list = []