From ed18b83fb7669a722d813bb8c6f0f817cdb21b5d Mon Sep 17 00:00:00 2001 From: Michal Adamkiewicz Date: Fri, 7 Jun 2019 19:57:00 -0700 Subject: [PATCH 01/14] Added compressed sensing tomography --- forest/benchmarking/tomography.py | 40 +++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index dd1506a6..43859546 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -163,6 +163,45 @@ def linear_inv_state_estimate(results: List[ExperimentResult], rho = pinv(measurement_matrix) @ expectations return unvec(rho) +def compressed_sensing_state_estimate(results: List[ExperimentResult], + qubits: List[int]) -> np.ndarray: + """ + Estimate a quantum state using compressed sensing + + [FLAMMIA] Quantum Tomography via Compressed Sensing: Error Bounds, Sample Complexity, and Efficient Estimators + Steven T. Flammia et. al. + (2012). + https://arxiv.org/pdf/1205.2300.pdf + + :param results: A tomographically complete list of results. + :param qubits: All qubits that were tomographized. This specifies the order in + which qubits will be kron'ed together. + :return: A point estimate of the quantum state rho. + """ + pauli_num = len(results) + for i in range(pauli_num): + #Convert the Pauli term into a tensor + r = results[i] + p_tensor = lifted_pauli(r.setting.out_operator, qubits) + e = r.expectation + #A[i] = e * scale_factor + pauli_list.append(p_tensor) + expectation_list.append(e) + + s = cp.Variable((d,d),complex = True) + obj = cp.Minimize(cp.norm(s, 'nuc')) + constraints = [cp.trace(s) == 1] + + for i in range(len(pauli_list)): + trace_bool = cp.trace(cp.matmul(pauli_list[i], s)) - expectation_list[i] == 0 + constraints.append(trace_bool) + + # Form and solve problem. + prob = cp.Problem(obj, constraints) + prob.solve() + rho = s.value + return rho + def iterative_mle_state_estimate(results: List[ExperimentResult], qubits: List[int], epsilon=.1, entropy_penalty=0.0, beta=0.0, tol=1e-9, maxiter=10_000) \ @@ -611,3 +650,4 @@ def _grad_cost(A, n, estimate, eps=1e-6): p = np.clip(p, a_min=eps, a_max=None) eta = n / p return unvec(-A.conj().T @ eta) + From f6aae0463d5b7c172262a573b61824592a6ae3c7 Mon Sep 17 00:00:00 2001 From: Michal Adamkiewicz Date: Fri, 7 Jun 2019 19:57:24 -0700 Subject: [PATCH 02/14] added fucntion to be used as debugger --- forest/benchmarking/tomography.py | 44 +++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index 43859546..76ed939e 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -651,3 +651,47 @@ def _grad_cost(A, n, estimate, eps=1e-6): eta = n / p return unvec(-A.conj().T @ eta) +def tomographize(program: Program, qubits: List[int], num_shots=1000, t_type='compressed_sensing', pauli_num=None): + """ + Runs tomography on the state generated by program and estimates the state. Can be used as a debugger + + :param program: which program to run the tomography on + :param quibits: whihc qubits to run the tomography debugger on + :param num_shots: the number of times to run each tomography experiment to get the expected value + :param t_type: which tomography type to use. Possible values: "linear_inv", "mle", "compressed_sensing" + :param pauli_num: the number of pauli matrices to use in the tomography + :return: the density matrix as an ndarray + """ + + # if no pauli_num is specified use the maximum + if pauli_num==None: + pauli_num=len(Qubits) + + #Generate experiments + qubit_experiments = generate_state_tomography_experiment(program=program, qubits=qubits) + + exp_list = [] + #Experiment holds all 2^n possible pauli matrices for the given number of qubits + #Now take pauli_num random pauli matrices as per the paper's advice + if (pauli_num > len(qubit_experiments)): + print("Cannot sample more Pauli matrices thatn d^2!") + return None + exp_list = random.sample(list(qubit_experiments), pauli_num) + input_exp = PyQuilTomographyExperiment(settings=exp_list, program=program) + + #Group experiments if possible to minimize QPU runs + #experiment = group_experiments(experiment) + + #NOTE: Change qvm depending on whether we are simulating qvm + qc = get_qc('%dq-qvm' % len(qubits)) + qc.compiler.quil_to_native_quil = basic_compile + + results = list(measure_observables(qc=qc, tomo_experiment=input_exp, n_shots=num_shots)) + + if t_type == 'compressed_sensing': + return compressed_sensing_state_estimate(results=results, qubits=qubits) + elif t_type == 'mle': + return itterative_mle_state_estimate(results=results, qubits=qubits) + elif t_type == "linear_inv": + return linear_inv_state_estimate(results=results, qubits=qubits) + From d8e65d2efe47fc6cd4ac913d2d498138b3e60a0a Mon Sep 17 00:00:00 2001 From: Michal Adamkiewicz Date: Fri, 7 Jun 2019 20:12:44 -0700 Subject: [PATCH 03/14] added examples files showing how to use debugger --- examples/tomography_debugger.ipynb | 223 +++++++++++++++++++++++++++++ 1 file changed, 223 insertions(+) create mode 100644 examples/tomography_debugger.ipynb diff --git a/examples/tomography_debugger.ipynb b/examples/tomography_debugger.ipynb new file mode 100644 index 00000000..5540ddeb --- /dev/null +++ b/examples/tomography_debugger.ipynb @@ -0,0 +1,223 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Tomography Debugger\n", + "State tomography involves measuring a quantum state repeatedly in the bases given by `itertools.product(['X', 'Y', 'Z], repeat=n_qubits)`. From these measurements, we can reconstruct a density matrix $\\rho$ using a varaiety of methods described in forest.benchmarking.tomography under the heading \"state tomography\". This is all done automaticly in using the forest.benchmarking.tomography.tomographize function allowing it to be use effectivly as a quantum debugger." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import time\n", + "\n", + "from pyquil import Program\n", + "from pyquil.gates import *\n", + "from forest.benchmarking.tomography import tomographize" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Construct a state with a `Program`\n", + "We'll construct a two-qubit graph state by Hadamarding all qubits and then applying a controlled-Z operation across edges of our graph. In the two-qubit case, there's only one edge. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "qubits = [0, 1]\n", + "\n", + "program = Program()\n", + "for qubit in qubits:\n", + " program += H(qubit)\n", + "program += CZ(qubits[0], qubits[1])\n", + "program += RY(-np.pi/2, qubits[0])\n", + "program += X(qubits[1])\n", + "program += CNOT(qubits[1], qubits[0])\n", + "\n", + "print(program)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run the tomography debugger and print output" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "start_linear = time.time()\n", + "m = 10\n", + "rho_linear = tomographize(program, qubits, pauli_num=10, t_type=\"linear_inv\")\n", + "end_linear = time.time() - start_linear\n", + "\n", + "print(\"Linear tomography took %gs\" % end_linear)\n", + "print(\"Recovered density matrix:\\n\")\n", + "print(rho_linear)\n", + "\n", + "start_compressed = time.time()\n", + "rho_compressed = tomographize(program, qubits, pauli_num=10, t_type=\"compressed_sensing\")\n", + "end_compressed = time.time() - start_compressed\n", + "print(\"Compressed tomography took %gs\" % end_compressed)\n", + "print(\"Recovered density matrix:\\n\")\n", + "print(rho_compressed)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Compare results to true output obtained using wavefunction simulator" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pyquil.api import WavefunctionSimulator\n", + "wf_sim = WavefunctionSimulator()\n", + "wf = wf_sim.wavefunction(program)\n", + "psi = wf.amplitudes\n", + "\n", + "rho_true = np.outer(psi, psi.T.conj())\n", + "print(np.around(rho_true, decimals=3))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualize using Hinton plots" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "from forest.benchmarking.plotting import hinton\n", + "fig, (ax1, ax2, ax3) = plt.subplots(1, 3)\n", + "hinton(rho_true, ax=ax1)\n", + "hinton(rho_linear, ax=ax2)\n", + "hinton(rho_compressed, ax=ax3)\n", + "ax1.set_title('Analytical Linear')\n", + "ax2.set_title('Estimated Linear')\n", + "ax3.set_title('Estimated Compressed')\n", + "fig.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Calculate matrix norm between true and estimated rho" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Linear norm:\")\n", + "print(np.linalg.norm(rho_linear - rho_true))\n", + "print(\"Compressed norm:\")\n", + "print(np.linalg.norm(rho_compressed - rho_true))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot graph of results for various measurement values" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "max_pauli_num = 4 ** len(qubits)\n", + "num_trials = 5\n", + "\n", + "linear_norms = []\n", + "compressed_norms = []\n", + "\n", + "print(\"Analyzing performance of linear vs. compressed on program:\")\n", + "print(program)\n", + "\n", + "for i in range(1, max_pauli_num):\n", + " print(\"Running iteration %d/%d\" % (i, max_pauli_num))\n", + " linear_norm_mean = 0.0\n", + " compressed_norm_mean = 0.0\n", + " for j in range(num_trials):\n", + " rho_linear = tomographize(program, qubits, pauli_num=i, t_type=\"linear_inv\")\n", + " rho_compressed = tomographize(program, qubits, pauli_num=i, t_type=\"compressed_sensing\")\n", + " linear_norm_mean += np.linalg.norm(rho_linear - rho_true)\n", + " compressed_norm_mean += np.linalg.norm(rho_compressed - rho_true)\n", + " \n", + " linear_norm_mean /= num_trials\n", + " compressed_norm_mean /= num_trials\n", + " \n", + " linear_norms.append(linear_norm_mean)\n", + " compressed_norms.append(compressed_norm_mean)\n", + "\n", + "plt.plot(linear_norms, label='linear')\n", + "plt.plot(compressed_norms, label='compressed')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} From a56627c5b35182001d3ca736cc13a7040fb8c435 Mon Sep 17 00:00:00 2001 From: Michal Adamkiewicz Date: Fri, 7 Jun 2019 23:32:36 -0700 Subject: [PATCH 04/14] added lasso and made qc be passable --- forest/benchmarking/tomography.py | 64 +++++++++++++++++++++++++++++-- 1 file changed, 60 insertions(+), 4 deletions(-) diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index 76ed939e..65d0554c 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -3,6 +3,7 @@ from operator import mul from typing import Callable, Tuple, List, Sequence import warnings +import random import numpy as np from scipy.linalg import logm, pinv, eigh @@ -179,6 +180,12 @@ def compressed_sensing_state_estimate(results: List[ExperimentResult], :return: A point estimate of the quantum state rho. """ pauli_num = len(results) + qubit_num = len(qubits) + d = 2 ** qubit_num + pauli_list = [] + expectation_list = [] + y = np.zeros((m,1)) + for i in range(pauli_num): #Convert the Pauli term into a tensor r = results[i] @@ -202,6 +209,53 @@ def compressed_sensing_state_estimate(results: List[ExperimentResult], rho = s.value return rho +def lasso_state_estimate(results: List[ExperimentResult], + qubits: List[int]) -> np.ndarray: + """ + Estimate a quantum state using compressed sensing + + [FLAMMIA] Quantum Tomography via Compressed Sensing: Error Bounds, Sample Complexity, and Efficient Estimators + Steven T. Flammia et. al. + (2012). + https://arxiv.org/pdf/1205.2300.pdf + + :param results: A tomographically complete list of results. + :param qubits: All qubits that were tomographized. This specifies the order in + which qubits will be kron'ed together. + :return: A point estimate of the quantum state rho. + """ + pauli_num = len(results) + qubit_num = len(qubits) + d = 2 ** qubit_num + pauli_list = [] + expectation_list = [] + y = np.zeros((m,1)) + + mu = 4 * pauli_num / np.sqrt(1000 * pauli_num) + + for i in range(m): + #Convert the Pauli term into a tensor + r = results[i] + p_tensor = lifted_pauli(r.setting.out_operator, qubits) + e = r.expectation + y[i] = e + pauli_list.append(p_tensor) + expectation_list.append(e) + + x = cp.Variable((d,d),complex = True) + A = cp.vstack([cp.trace(cp.matmul(pauli_list[i], x)) * np.sqrt(d / m) for i in range(m)]) + + #Minimize trace norm + obj = cp.Minimize(0.5 * cp.norm((A - y), 2) + mu * cp.norm(x, 'nuc')) + constraints = [cp.trace(x) == 1] + + # Form and solve problem. + prob = cp.Problem(obj, constraints) + prob.solve() + rho = x.value + + return rho + def iterative_mle_state_estimate(results: List[ExperimentResult], qubits: List[int], epsilon=.1, entropy_penalty=0.0, beta=0.0, tol=1e-9, maxiter=10_000) \ @@ -651,14 +705,15 @@ def _grad_cost(A, n, estimate, eps=1e-6): eta = n / p return unvec(-A.conj().T @ eta) -def tomographize(program: Program, qubits: List[int], num_shots=1000, t_type='compressed_sensing', pauli_num=None): +def tomographize(qc, program: Program, qubits: List[int], num_shots=1000, t_type='lasso', pauli_num=None): """ Runs tomography on the state generated by program and estimates the state. Can be used as a debugger + :param qc: the quantum computer to run the debugger on :param program: which program to run the tomography on :param quibits: whihc qubits to run the tomography debugger on :param num_shots: the number of times to run each tomography experiment to get the expected value - :param t_type: which tomography type to use. Possible values: "linear_inv", "mle", "compressed_sensing" + :param t_type: which tomography type to use. Possible values: "linear_inv", "mle", "compressed_sensing", "lasso" :param pauli_num: the number of pauli matrices to use in the tomography :return: the density matrix as an ndarray """ @@ -683,7 +738,6 @@ def tomographize(program: Program, qubits: List[int], num_shots=1000, t_type='co #experiment = group_experiments(experiment) #NOTE: Change qvm depending on whether we are simulating qvm - qc = get_qc('%dq-qvm' % len(qubits)) qc.compiler.quil_to_native_quil = basic_compile results = list(measure_observables(qc=qc, tomo_experiment=input_exp, n_shots=num_shots)) @@ -691,7 +745,9 @@ def tomographize(program: Program, qubits: List[int], num_shots=1000, t_type='co if t_type == 'compressed_sensing': return compressed_sensing_state_estimate(results=results, qubits=qubits) elif t_type == 'mle': - return itterative_mle_state_estimate(results=results, qubits=qubits) + return iterative_mle_state_estimate(results=results, qubits=qubits) elif t_type == "linear_inv": return linear_inv_state_estimate(results=results, qubits=qubits) + elif t_type == "lasso": + return lasso_state_estimate(results=results, qubits=qubits) From 03568b5b60621a5e68648ff091b8ea6f9080612f Mon Sep 17 00:00:00 2001 From: Michal Adamkiewicz Date: Fri, 7 Jun 2019 23:35:18 -0700 Subject: [PATCH 05/14] adde lasso to example and proper qc --- examples/tomography_debugger.ipynb | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/examples/tomography_debugger.ipynb b/examples/tomography_debugger.ipynb index 5540ddeb..19a6320a 100644 --- a/examples/tomography_debugger.ipynb +++ b/examples/tomography_debugger.ipynb @@ -17,7 +17,7 @@ "import numpy as np\n", "import time\n", "\n", - "from pyquil import Program\n", + "from pyquil import Program, get_qc\n", "from pyquil.gates import *\n", "from forest.benchmarking.tomography import tomographize" ] @@ -62,9 +62,12 @@ "metadata": {}, "outputs": [], "source": [ + "qc = get_qc('%dq-qvm' % len(qubits))\n", + "\n", + "\n", "start_linear = time.time()\n", "m = 10\n", - "rho_linear = tomographize(program, qubits, pauli_num=10, t_type=\"linear_inv\")\n", + "rho_linear = tomographize(qc, program, qubits, pauli_num=10, t_type=\"linear_inv\")\n", "end_linear = time.time() - start_linear\n", "\n", "print(\"Linear tomography took %gs\" % end_linear)\n", @@ -72,11 +75,18 @@ "print(rho_linear)\n", "\n", "start_compressed = time.time()\n", - "rho_compressed = tomographize(program, qubits, pauli_num=10, t_type=\"compressed_sensing\")\n", + "rho_compressed = tomographize(qc, program, qubits, pauli_num=10, t_type=\"compressed_sensing\")\n", "end_compressed = time.time() - start_compressed\n", "print(\"Compressed tomography took %gs\" % end_compressed)\n", "print(\"Recovered density matrix:\\n\")\n", - "print(rho_compressed)" + "print(rho_compressed)\n", + "\n", + "start_lasso = time.time()\n", + "rho_lasso = tomographize(qc, program, qubits, pauli_num=10, t_type=\"lasso\")\n", + "end_lasso = time.time() - start_lasso\n", + "print(\"Compressed tomography took %gs\" % end_lasso)\n", + "print(\"Recovered density matrix:\\n\")\n", + "print(rho_lasso)" ] }, { @@ -142,7 +152,9 @@ "print(\"Linear norm:\")\n", "print(np.linalg.norm(rho_linear - rho_true))\n", "print(\"Compressed norm:\")\n", - "print(np.linalg.norm(rho_compressed - rho_true))" + "print(np.linalg.norm(rho_compressed - rho_true))\n", + "print(\"Lasso norm:\")\n", + "print(np.linalg.norm(rho_lasso - rho_true))" ] }, { @@ -174,8 +186,8 @@ " linear_norm_mean = 0.0\n", " compressed_norm_mean = 0.0\n", " for j in range(num_trials):\n", - " rho_linear = tomographize(program, qubits, pauli_num=i, t_type=\"linear_inv\")\n", - " rho_compressed = tomographize(program, qubits, pauli_num=i, t_type=\"compressed_sensing\")\n", + " rho_linear = tomographize(qc, program, qubits, pauli_num=i, t_type=\"linear_inv\")\n", + " rho_compressed = tomographize(qc, program, qubits, pauli_num=i, t_type=\"compressed_sensing\")\n", " linear_norm_mean += np.linalg.norm(rho_linear - rho_true)\n", " compressed_norm_mean += np.linalg.norm(rho_compressed - rho_true)\n", " \n", From 5c43b9fc60d407d82d349b2374033cb8dfa1e841 Mon Sep 17 00:00:00 2001 From: Michal Adamkiewicz Date: Fri, 7 Jun 2019 23:36:26 -0700 Subject: [PATCH 06/14] added group expirents --- forest/benchmarking/tomography.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index 65d0554c..e4b8bee2 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -735,7 +735,7 @@ def tomographize(qc, program: Program, qubits: List[int], num_shots=1000, t_type input_exp = PyQuilTomographyExperiment(settings=exp_list, program=program) #Group experiments if possible to minimize QPU runs - #experiment = group_experiments(experiment) + input_exp = group_experiments(input_exp) #NOTE: Change qvm depending on whether we are simulating qvm qc.compiler.quil_to_native_quil = basic_compile From 2072af726aa42bcd0079247cb14ea9238352f057 Mon Sep 17 00:00:00 2001 From: Michal Adamkiewicz Date: Fri, 7 Jun 2019 23:41:57 -0700 Subject: [PATCH 07/14] fixed imports --- forest/benchmarking/tomography.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index e4b8bee2..78160dcc 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -11,10 +11,12 @@ import forest.benchmarking.distance_measures as dm from forest.benchmarking.superoperator_tools import vec, unvec, proj_choi_to_physical from forest.benchmarking.utils import n_qubit_pauli_basis +from forest.benchmarking.compilation import basic_compile from pyquil import Program from pyquil.operator_estimation import ExperimentSetting, \ TomographyExperiment as PyQuilTomographyExperiment, ExperimentResult, SIC0, SIC1, SIC2, SIC3, \ plusX, minusX, plusY, minusY, plusZ, minusZ, TensorProductState, zeros_state +from pyquil.operator_estimation import group_experiments from pyquil.paulis import sI, sX, sY, sZ, PauliSum, PauliTerm, is_identity from pyquil.unitary_tools import lifted_pauli as pauli2matrix, lifted_state_operator as state2matrix From af18643c2a638cdb1c461175f9473376b8732c4a Mon Sep 17 00:00:00 2001 From: Michal Adamkiewicz Date: Fri, 7 Jun 2019 23:47:44 -0700 Subject: [PATCH 08/14] added docs --- docs/tomography.rst | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/docs/tomography.rst b/docs/tomography.rst index 4f46e50a..999c5af3 100644 --- a/docs/tomography.rst +++ b/docs/tomography.rst @@ -77,6 +77,14 @@ Finally, we analyze our data with one of the analysis routines:: [-0.1869-0.0077j -0.1794-0.0188j -0.169 -0.0169j 0.2202-0.j ]] Purity = (0.6889520199999999+4.597017211338539e-17j) +Debugger +~~~~~~~~~ + +The above steps can be automated to create a basic debugger that can be used to +peek into the state of a program running on a qc. This can be done using the +tomographize function:: + + rho = tomographize(qc, program, qubits, pauli_num=10, t_type="compressed_sensing") API Reference ------------- @@ -96,4 +104,5 @@ API Reference iterative_mle_state project_density_matrix estimate_variance + tomographize From f05afe14a71af9ea5436cd504413eac628100711 Mon Sep 17 00:00:00 2001 From: Michal Adamkiewicz Date: Sat, 8 Jun 2019 00:06:41 -0700 Subject: [PATCH 09/14] added missing imports --- forest/benchmarking/tomography.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index 78160dcc..0d997637 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -7,6 +7,7 @@ import numpy as np from scipy.linalg import logm, pinv, eigh +import cvxpy as cp import forest.benchmarking.distance_measures as dm from forest.benchmarking.superoperator_tools import vec, unvec, proj_choi_to_physical @@ -16,7 +17,7 @@ from pyquil.operator_estimation import ExperimentSetting, \ TomographyExperiment as PyQuilTomographyExperiment, ExperimentResult, SIC0, SIC1, SIC2, SIC3, \ plusX, minusX, plusY, minusY, plusZ, minusZ, TensorProductState, zeros_state -from pyquil.operator_estimation import group_experiments +from pyquil.operator_estimation import group_experiments, measure_observables from pyquil.paulis import sI, sX, sY, sZ, PauliSum, PauliTerm, is_identity from pyquil.unitary_tools import lifted_pauli as pauli2matrix, lifted_state_operator as state2matrix @@ -191,7 +192,7 @@ def compressed_sensing_state_estimate(results: List[ExperimentResult], for i in range(pauli_num): #Convert the Pauli term into a tensor r = results[i] - p_tensor = lifted_pauli(r.setting.out_operator, qubits) + p_tensor = pauli2matrix(r.setting.out_operator, qubits) e = r.expectation #A[i] = e * scale_factor pauli_list.append(p_tensor) @@ -238,7 +239,7 @@ def lasso_state_estimate(results: List[ExperimentResult], for i in range(m): #Convert the Pauli term into a tensor r = results[i] - p_tensor = lifted_pauli(r.setting.out_operator, qubits) + p_tensor = pauli2matrix(r.setting.out_operator, qubits) e = r.expectation y[i] = e pauli_list.append(p_tensor) From 998c5db8efd718ebdcf46aa5530bee513a3aa45d Mon Sep 17 00:00:00 2001 From: Joshua Combes Date: Wed, 26 Jun 2019 08:22:23 +1000 Subject: [PATCH 10/14] from Kyle's review Qubits -> qubits --- forest/benchmarking/tomography.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index 88c74f76..2356bec4 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -723,7 +723,7 @@ def tomographize(qc, program: Program, qubits: List[int], num_shots=1000, t_type # if no pauli_num is specified use the maximum if pauli_num==None: - pauli_num=len(Qubits) + pauli_num=len(qubits) #Generate experiments qubit_experiments = generate_state_tomography_experiment(program=program, qubits=qubits) From 0e41c6eba819934e55f0ef8a31d8850678570c00 Mon Sep 17 00:00:00 2001 From: Joshua Combes Date: Wed, 26 Jun 2019 09:51:55 +1000 Subject: [PATCH 11/14] edits from the code review and other get the code working again due to changes in forest benchmarking --- examples/tomography_debugger.ipynb | 133 ++++++++++++++++++++++------- forest/benchmarking/tomography.py | 51 ++++++----- 2 files changed, 133 insertions(+), 51 deletions(-) diff --git a/examples/tomography_debugger.ipynb b/examples/tomography_debugger.ipynb index e2853951..4de5bd5b 100644 --- a/examples/tomography_debugger.ipynb +++ b/examples/tomography_debugger.ipynb @@ -79,30 +79,39 @@ "name": "stdout", "output_type": "stream", "text": [ - "Linear tomography took 1.5077s\n", + "Linear tomography took 0.460453s\n", "Recovered density matrix:\n", "\n", - "[[ 0.249-8.57061622e-18j 0.25 -7.50000000e-03j 0.003+1.30000000e-02j\n", - " 0.012-5.00000000e-03j]\n", - " [ 0.25 +7.50000000e-03j 0.251-1.42843604e-17j 0.005-5.00000000e-03j\n", - " -0.017-3.00000000e-03j]\n", - " [ 0.003-1.30000000e-02j 0.005+5.00000000e-03j 0.249-2.85687207e-18j\n", - " -0.25 +7.50000000e-03j]\n", - " [ 0.012+5.00000000e-03j -0.017+3.00000000e-03j -0.25 -7.50000000e-03j\n", - " 0.251+2.85687207e-18j]]\n" - ] - }, - { - "ename": "NameError", - "evalue": "name 'm' is not defined", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", - "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[1;32m 12\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 13\u001b[0m \u001b[0mstart_compressed\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 14\u001b[0;31m \u001b[0mrho_compressed\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtomographize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mqc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprogram\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mqubits\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpauli_num\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_type\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"compressed_sensing\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 15\u001b[0m \u001b[0mend_compressed\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mstart_compressed\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 16\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Compressed tomography took %gs\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mend_compressed\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/forest-benchmarking/forest/benchmarking/tomography.py\u001b[0m in \u001b[0;36mtomographize\u001b[0;34m(qc, program, qubits, num_shots, t_type, pauli_num)\u001b[0m\n\u001b[1;32m 747\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 748\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mt_type\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'compressed_sensing'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 749\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mcompressed_sensing_state_estimate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresults\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mresults\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mqubits\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mqubits\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 750\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0mt_type\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'mle'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 751\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0miterative_mle_state_estimate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresults\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mresults\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mqubits\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mqubits\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;32m~/forest-benchmarking/forest/benchmarking/tomography.py\u001b[0m in \u001b[0;36mcompressed_sensing_state_estimate\u001b[0;34m(results, qubits)\u001b[0m\n\u001b[1;32m 182\u001b[0m \u001b[0mpauli_list\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 183\u001b[0m \u001b[0mexpectation_list\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 184\u001b[0;31m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mm\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 185\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 186\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpauli_num\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", - "\u001b[0;31mNameError\u001b[0m: name 'm' is not defined" + "[[ 3.50000000e-03-1.36512468e-18j 2.93058701e-19+1.50000000e-03j\n", + " 4.50000000e-03-3.50000000e-03j 5.50000000e-03-5.00000000e-04j]\n", + " [-2.74342100e-17-1.50000000e-03j -3.50000000e-03+2.14325778e-18j\n", + " -1.45000000e-02+5.00000000e-04j 4.50000000e-03+1.55000000e-02j]\n", + " [ 4.50000000e-03+3.50000000e-03j -1.45000000e-02-5.00000000e-04j\n", + " 4.96500000e-01+1.13981602e-18j 6.83320553e-19-5.00000000e-04j]\n", + " [ 5.50000000e-03+5.00000000e-04j 4.50000000e-03-1.55000000e-02j\n", + " -6.83320553e-19+5.00000000e-04j 5.03500000e-01+5.30658834e-19j]]\n", + "Compressed tomography took 0.499431s\n", + "Recovered density matrix:\n", + "\n", + "[[-2.53127791e-08-3.50638546e-17j -8.70056892e-03-1.64999653e-02j\n", + " 1.21925031e-08-4.21747408e-04j 8.39998555e-03+2.40000386e-02j]\n", + " [-8.70056892e-03+1.64999653e-02j 4.71000042e-01+6.24351317e-18j\n", + " 1.10522052e-02+2.00001772e-03j -5.00000056e-01-1.94216460e-02j]\n", + " [ 1.21917260e-08+4.21747408e-04j 1.10522052e-02-2.00001772e-03j\n", + " 3.82506459e-08+3.77619537e-16j -1.16478627e-02-5.50002269e-03j]\n", + " [ 8.39998555e-03-2.40000386e-02j -5.00000056e-01+1.94216460e-02j\n", + " -1.16478627e-02+5.50002269e-03j 5.29000052e-01+7.61888035e-17j]]\n", + "Compressed tomography took 0.483477s\n", + "Recovered density matrix:\n", + "\n", + "[[ 3.44521584e-01+3.17512724e-17j 3.22732238e-03+2.82322119e-03j\n", + " -3.43649473e-01-1.89699327e-02j -1.18441881e-02+1.23272052e-02j]\n", + " [ 3.22732238e-03-2.82322119e-03j 1.48123284e-01-4.03748157e-17j\n", + " 3.35119267e-04+8.78427542e-03j -1.54915266e-01+4.86985721e-03j]\n", + " [-3.43649473e-01+1.89699327e-02j 3.35119267e-04-8.78427542e-03j\n", + " 3.44227157e-01+4.59209697e-17j 7.46896847e-03-6.41707610e-03j]\n", + " [-1.18441881e-02-1.23272052e-02j -1.54915266e-01-4.86985721e-03j\n", + " 7.46896847e-03+6.41707610e-03j 1.63110900e-01-4.77300410e-17j]]\n" ] } ], @@ -143,9 +152,20 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[ 0. +0.j -0. -0.j 0. +0.j 0. +0.j]\n", + " [-0. +0.j 0.5+0.j -0. +0.j -0.5+0.j]\n", + " [ 0. +0.j -0. -0.j 0. +0.j 0. +0.j]\n", + " [ 0. +0.j -0.5-0.j 0. +0.j 0.5+0.j]]\n" + ] + } + ], "source": [ "from pyquil.api import WavefunctionSimulator\n", "wf_sim = WavefunctionSimulator()\n", @@ -165,7 +185,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -190,9 +210,22 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Linear norm:\n", + "1.0050482575478652\n", + "Compressed norm:\n", + "0.07078083250036442\n", + "Lasso norm:\n", + "0.9749834380713842\n" + ] + } + ], "source": [ "print(\"Linear norm:\")\n", "print(np.linalg.norm(rho_linear - rho_true))\n", @@ -211,9 +244,51 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Analyzing performance of linear vs. compressed on program:\n", + "H 0\n", + "H 1\n", + "CZ 0 1\n", + "RY(-pi/2) 0\n", + "X 1\n", + "CNOT 1 0\n", + "\n", + "Running iteration 1/16\n", + "Running iteration 2/16\n", + "Running iteration 3/16\n", + "Running iteration 4/16\n", + "Running iteration 5/16\n", + "Running iteration 6/16\n", + "Running iteration 7/16\n", + "Running iteration 8/16\n", + "Running iteration 9/16\n", + "Running iteration 10/16\n", + "Running iteration 11/16\n", + "Running iteration 12/16\n", + "Running iteration 13/16\n", + "Running iteration 14/16\n", + "Running iteration 15/16\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd4VFX+x/H3Sa8kIY2ShITeCRCQGlBEQcCKAioroui6lrWsrruu3d111XVXV38qq4IVRGQVFSyIIlKULhB6SQg1JCSkkjLn98cZIEBCBjKTOzP5vp4nT6bc3PudED5z5txzz1Faa4QQQngXH6sLEEII4XwS7kII4YUk3IUQwgtJuAshhBeScBdCCC8k4S6EEF5Iwl0IIbyQhLsQQnghCXchhPBCflYdOCYmRicnJ1t1eCGE8EirVq06rLWOrWs7y8I9OTmZlStXWnV4IYTwSEqpTEe2k24ZIYTwQhLuQgjhhSTchRDCC1nW516TiooKsrOzKSsrs7qURikoKIiEhAT8/f2tLkUIUU9uFe7Z2dmEh4eTnJyMUsrqchoVrTW5ublkZ2eTkpJidTlCiHpyq26ZsrIyoqOjJdgtoJQiOjpaPjUJ4SXcKtwBCXYLye9eCO/hduEuhBDeKqfwGG8s2sGyHbkuP5aE+2nCwsIA2LdvH2PHjrW4GiGEp6ussrFw80Fuf28l/f/+HX+fv5lFW3Ncftw6T6gqpd4GRgOHtNZda3j+BuCP9rtFwB1a63VOrdICLVq0YPbs2S49RmVlJX5+bnVOWwjhJJm5xcxauYfZq7I5ePQY0aEBTB6UwnVpCbSNC3f58R1JlunAK8C7tTy/CxiitT6ilBoJTAUucE551tm9ezejR49mw4YNTJ8+nblz51JSUsKOHTu46qqreO655wD45ptvePzxxzl27Bht2rRh2rRphIWF8dRTT/H5559TWlrKgAEDeOONN1BKMXToUAYMGMCSJUu4/PLLeeCBByx+pUIIZymrqOKrDQf4aMUelu3MxUfBkPaxPHl5Ihd1jCfAr+E6S+oMd631j0qp5LM8v7Ta3eVAQv3Lgic/30jGvqPO2NUJnVs04fExXc7rZ9euXcuaNWsIDAykQ4cO3H333QQHB/PMM8+wYMECQkND+cc//sGLL77IY489xl133cVjjz0GwMSJE/niiy8YM2YMAPn5+SxatMhpr0sIYa0Newv4aMUePlu7l6NllSQ1DeEPl7Tnmt4JNI8ItqQmZ/cJ3ALMd/I+3cKwYcOIiIgAoHPnzmRmZpKfn09GRgYDBw4EoLy8nP79+wPw/fff89xzz1FSUkJeXh5dunQ5Ee7jxo2z5kUIIZymoKSCz9btZeYve8jYf5QAPx9Gdm3GuD6J9EuJxsfH2tFnTgt3pdSFmHAfdJZtbgNuA0hKSjrr/s63he0qgYGBJ277+vpSWVmJ1prhw4czY8aMU7YtKyvjd7/7HStXriQxMZEnnnjilPHjoaGhDVa3EMJ5bDbN8p25fLRyD/M3HKC80kaXFk146oouXNGjJREh7nN1t1PCXSnVHXgTGKm1rnWMj9Z6KqZPnrS0NO2MY1upX79+3HnnnWzfvp22bdtSUlJCdnY2cXFxAMTExFBUVMTs2bNl5I0QTnassoqN+46yJiufPXklhAf5ERHsT0SwP5EhAUSG+BNpvx8R4k+gn+95H2t/QSmzV2bz8apssvJKaBLkx/g+iVyXlkjXlhFOfFXOU+9wV0olAXOAiVrrrfUvyXPExsYyffp0JkyYwLFjxwB45plnaN++PVOmTKFbt24kJyfTp08fiysVwrNprck+UsqaPfmsyTrCmqx8MvYdpbzKBkB4oB/F5ZXYztJkDPb3JTLEv9obgD+RwQFEhJx23347ItifjftMX/qirTnYNPRvHc39w9szomszgvzP/82iISitz96AVkrNAIYCMcBB4HHAH0Br/bpS6k3gGuD4BPKVWuu0ug6clpamT1+sY9OmTXTq1OkcX4JwJvk3EO6g+Fglv2YXsGaPCfI1WfkcLjINqCB/H7onRNIzKZKeiVH0TIokvkkQNpum8FglR0sryC+pIL+0nPySCgpKzVd+ycn7+aUVFNi3KSitoKzCVmstzZoEMbZ3AtemJdAq2vouVaXUKkcy1pHRMhPqeP5W4NZzqE0IIU6w2TQ7Dxex2h7ia/fks+XA0ROt8NYxoaS3j6FnUhQ9EyPp0Cwcf98zhxT6+KgTrfLEpudWQ1lFlf0N4NQ3grgmgQxuF4uvxSdHz4dcQSOEaFD5JeX27hXTxbJ2Tz6FZZUAhAf5kZoYyfCL2tEzKZLUhEiiQgNcXlOQvy9B/r7ENwly+bEaioS7EKJB/Lg1h2e+zGDrwSIAfBR0aNaEMT1a0DMxkp5JUbSOCbV8CKG3kHAXQrhUWUUVz87fzPSlu2kTG8pDIzrQMzGK7gkRhAZKBLmK/GaFEC6zcV8B985cy7ZDRUwakMzDIzu6/SgTbyHhLoRwuiqb5s3FO3nhmy1EhQTwzuS+DGkfa3VZjYpM+dsI/fDDD4wePdrqMoSX2ptfyvX/Xc7f529mWMd4vr43XYLdAtJydxMy/a/wBp+t3ctfPt2AzaZ5bmx3ru2dICt8WURa7jV499136d69Oz169GDixIlkZmYybNgwunfvzrBhw8jKygJg0qRJ3HHHHVx44YW0bt2aRYsWMXnyZDp16sSkSZNO7C8sLIwHHniAXr16MWzYMHJyzET9Q4cO5c9//jNDhgzhpZdeIicnh2uuuYY+ffrQp08flixZAsCiRYtITU0lNTWVnj17UlhYyP79+0lPTyc1NZWuXbuyePFiwExB3L9/f3r16sW1115LUZEZmfDVV1/RsWNHBg0axJw5cxrwtymcTWvNocIylm4/zDtLdzNr5R6Kj1VaWlNBSQX3zFjD72eupX18OPN/n851aYkS7BZy36bi/IfhwHrn7rNZNxj57Fk32bhxI3/9619ZsmQJMTEx5OXlcdNNN/Gb3/yGm266ibfffpt77rmHTz/9FIAjR46wcOFC5s6dy5gxY1iyZAlvvvkmffr0Ye3ataSmplJcXEyvXr345z//yVNPPcWTTz7JK6+8Apw6/e/111/Pfffdx6BBg8jKyuLSSy9l06ZNvPDCC7z66qsMHDiQoqIigoKCmDp1KpdeeimPPPIIVVVVlJSUcPjw4RqnIH7ooYeYMmUKCxcupG3btjIrpYcwIX6MbQeL2HaokK0Hi9h+qJBth4rIL6k4ZdunP8/g6l4tubFfK9rFu34hiOqW7jjMH2at41DhMR4Y3p47hrbBr4aLjETDct9wt8jChQsZO3YsMTExADRt2pRly5adaO1OnDiRhx566MT2Y8aMQSlFt27diI+Pp1u3bgB06dKF3bt3k5qaio+Pz4lAvfHGG7n66qtP/Hz1oF2wYAEZGRkn7h89epTCwkIGDhzI/fffzw033MDVV19NQkICffr0YfLkyVRUVHDllVeSmprKokWLapyCePPmzaSkpNCuXbsTNUydOtUVvz5xHrTWHDhaxraDRWw9WMj2Q0VsO1TEtoOFHC072SKPDPGnfVw4l3VrTru4MNrFhdMuPozsIyW8tyyTGb/s4Z1lmfRNacqN/Voxokszly4Ocayyihe/2crUxTtJiQ7lkzsG0CMx0mXHE+fGfcO9jha2q2it6/woWf3541MB+/j4nDItsI+PD5WVNX9Urv7z1af/tdlsLFu2jODgUyf3f/jhhxk1ahTz5s2jX79+LFiwgPT0dH788Ue+/PJLJk6cyIMPPkhUVFSNUxCvXbtWPh67Aa01+wrK2Haw8ERrfNuhIrYfLKKwWrdKdGgAbePCuDy1Be3jw2lrD/KYsIAa/x3jmwTRu1VTHh19jI9XZfPBz5ncM2MNMWEBjOuTyIS+SSREhTj1tWw5UMi9H61l0/6j3HBBEo+M6kRIgPvGSWMk/xqnGTZsGFdddRX33Xcf0dHR5OXlMWDAAGbOnMnEiRP54IMPGDSo1inra2Sz2Zg9ezbjx4/nww8/rPXnL7nkEl555RUefPBBgBPdOjt27KBbt25069aNZcuWsXnzZoKDg2nZsiVTpkyhuLiY1atX88gjj9Q4BXHHjh3ZtWsXO3bsoE2bNmeEv3CdnMJjLNx8kG8zDrF8Zy5F1UI8JiyQdnFhXN2rJW3jw+2t8TCiwwLPssfaRYcF8tshbbhtcGsWbcvhg+WZvPbDDl77YQcXdojjxn6tSG9fv3lSbDbN9KW7efarzYQH+vHWTWkM6xR/3vsTriPhfpouXbrwyCOPMGTIEHx9fenZsycvv/wykydP5vnnnyc2NpZp06ad0z5DQ0PZuHEjvXv3JiIigo8++qjG7V5++WXuvPNOunfvTmVlJenp6bz++uv8+9//5vvvv8fX15fOnTszcuRIZs6cyfPPP4+/vz9hYWG8++67Z52CeOrUqYwaNYqYmBgGDRrEhg0b6v27EmfSWrMjp4hvMw7xbcYB1uzJR2toGRnMFakt6NyiielOiQtz2ZwpPj6KCzvEcWGHOPbmlzLj5yxmrtjDd9NXkNg0mOv7tuK6tIRzfhM5UFDGg7PXsXjbYYZ1jOPZa7oTG35+b0TC9eqc8tdVGtOUv2FhYSdGrbg7b/03cKXKKhurs/L5NuMACzYdYtfhYgC6tYxgeOd4Lu4UT6fm4ZZ2jZVX2vgm4wDvL89k+c48Anx9GNmtGTf2a0Vaq6g6a5u3fj9/mrOe8kobj47uzIS+MhLGKk6b8lcIcabiY5Us3pbDtxmHWLj5IEdKKgjw9aF/m2gmD0rh4k5xli2MXJMAPx9Gd2/B6O4t2HawkA9+zuKTVdl8tnYfHZuFc0O/VlzVsyVhp831UlhWwRNzM/hkdTY9EiL417hUWseGWfQqxLmQlrs4hfwb1O7g0TIWbDrIgoyDLNmRS3mljYhgfy7qGMfwzvGkt489IxzdWUl5JXPX7uO95Zls3HeU0ABfruxphlN2at6EFbvzuO+jtezLL+WuC9ty97B2Nc6jLhqWx7bcHRmtIlzDqjd6d6W1ZsvBQhZkHOTbjIOsyy4AIKlpCBP7teLiTvH0SY7y2DHdIQF+jO+bxLg+iazdk8/7y7OYvSqbD37OonPzJmw+cJSWUcF8/Nv+9G51jqtfCMu5VbgHBQWRm5tLdHS0BHwD01qTm5tLUJD3LFZwriqqbBwoKGPX4WK+33KIBZsOsievFIDUxEgevLQDwzvH0y4uzKv+PpVSZpWjpCgeHd2J2auy+XTtXsb1SeSRUZ096tOIOMmtumUqKirIzs6mrKzMkpoau6CgIBISEvD397e6FJcoLKtgb34p+/JL2XuklL35Zea2/bGDR8tOLO0W4OfDoLYxDO8cz7COccR50Qo9wrN5ZLeMv78/KSkpVpchPFCVTZNTeIy9+SUnQ/tI6Ynw3ptfemIpt+P8fRXNI4JpGRnMgDYxtIwKpmVkEAlRIfRMipSLcoRHk79eD7Qzp4gnPs/Az0fRKjqE5OhQkmNCSY4OoWVksMf2AZ+LPXklTF+6mw17C9hXUMr+/DIqbad+Co0I9qdFZDAJUcFckNKUFpHBtIwKNo9FBhMTFihLugmvJeHuYVZl5nHrOyvRQPOIYJbvzKWkvOrE834+ioSoYHvYh5rwt99OiAr2+NEOG/cV8MainXy5fj8+CnokRNIrKYoW3U0LvKU9wJtHBBEe5J3dS0I4QsLdg3y14QC/n7mG5hFBTL+5L8kxoWhtuiN255aw+3Axu3OLycwtYXduMSt25VFcLfh97cHfKtq08k2LP4RW0aEkRoW4dJKp+tBas2xHLq//uJMft+YQFujHLYNSmDwwhWYR0hcuRE0k3D3EtCW7eOqLDFITI3nzN2knLh1XShHXJIi4JkH0TTl1uJrWmsNF5WTmFp8R/msyj5wyWZWPgpZRwbSNDWNYp3hGdG1GzHnOceIsVTbN1xsP8MaiHazLLiAmLJAHL+3Ajf1aEREsrXIhzqbO0TJKqbeB0cAhrXXXGp5XwEvAZUAJMElrvbquA9c0WkacyWbT/G3eJt78aRfDO8fz8vieBAfUf4FhrTV5xeXszi0x4X/YvAGs31vArsPF+Cjom9KUUd2ac2mXZg06WqSsooo5q/fy38U72XW4mOToEKakt+aaXgmyuLJo9BwdLeNIuKcDRcC7tYT7ZcDdmHC/AHhJa31BXQeWcK9bWUUVD8xax5fr93NT/1Y8NqZLvWb0c8TxC3fmrT/AvPX72X6oCKWgT6umXNatGSO6NndZV8jRsgreX57JtCW7ySk8RveECH47pA2Xdmnm8tcthKdwWrjbd5YMfFFLuL8B/KC1nmG/vwUYqrXef7Z9SrifXX5JOVPeXcmK3Ud45LJO3Do4xZILZ7YdLOTL9fuZv/4AWw4WAtC7VRSXdWvOyK7NaBFZ//lTDh4t4+2fdvHBz1kUHatkcLsY7hjShv5t5GI2IU7XkOPcWwJ7qt3Ptj921nAXtduTV8Kkab+wJ6+U/0zoyZgeLSyrpV18OPfGh3Pvxe3ZfqiI+ev3M2/DAZ7+IoOnv8igZ1Ikl3VtzoiuzUhsem4LQuzIKWLqop38b81eKm02RnVvwe3prenaMsJFr0aIxsMZ4V5T06rGjwNKqduA2wCSkpKccGjvsz67gJunr6C8sor3bunLBa2jrS7phLZxYdw9rB13D2vHrsPFzFu/n/kb9vPXeZv467xN9EiIYGS35lzWtTlJ0bUH/eqsI7yxaAffZBwkwNeHcX0SmTK49Vl/RghxbqRbxo18v/kQd364mqiQAN6Z3Ie2cQ270PH5yswtZv6GA8xfv//E5FpdWzZhZNfmjOrW/MSQzR+25PDaoh38siuPiGB/ftO/FTcNSLZ8VI4QnqQh+9xHAXdx8oTqy1rrvnXtU8L9VDN+yeIvn26gU/Nw3r6pj8fOZbInr4SvNhxg3ob9rMnKB6BT8yZordl8oJAWEUHcMrg14/skEioTUglxzpw5WmYGMBSIAQ4CjwP+AFrr1+1DIV8BRmCGQt6sta4ztSXcDa01L367lf8s3M6Q9rH83w29vCb09uaX8pW9RX+s0sakAclcntrC46+SFcJKTm25u4KEu1n67OFPfmXOmr2MS0vkmau6SvAJIc7KI2eFbEyOllVwx/urWLI9l/uHt+fui9rKsD8hhNNIuFtgf0EpN09bwfZDRbxwbQ/G9k6wuiQhhJeRcG9gmw8cZdLbKyg6Vsm0m/swuF2s1SUJIbyQhHsDWrL9ML99bxUhgb7Mur0/nVs0sbokIYSXknBvIHNWZ/PHT36ldUwY027u45TL9oUQojYS7i6mteb/ftjB819voX/raF6f2FumqxVCuJyEez2UV9rILy2noKSCIyUV5JeUk19SQX5puf1+BZm5xSzdkcuVqS14bmwPt10QQwjhXRpFuGutsWmz+EOVTVOl9cnbNo1NayptmrKKKhPO9pA+UlJOQan5nm8P65O3y09Z5eh0/r6KiOAAokL8ZaijEKLBeVy4/7DlEE99kYHteEhXHQ9rqLLZ7GENlTYbNhsngvx8+SiIDAkgMsSfyGB/mjUJokOzcKJCAogM9icy1HyPOr5NiD+RIQGEBvhKmAshLONx4d4k2J9OzZvgqxR+PgofH4WvUvj62r/7nPalzDZ+tTx2/Of9fBSB/j5E2IM6KiSAiBB/wgP98JGFIoQQHsbjwr1XUhS9ro+yugwhhHBrcnZPCCG8kIS7EEJ4IQl3IYTwQhLuQgjhhSTchRDCC0m4CyGEF5JwF0IILyThLoQQXkjCXQghvJCEuxBCeCHPC/eCvbDqHTi8DfT5TwjmcjYb7P8Vfvkv7FtrdTVCiEbG4+aWYdci+Pweczs0FpL6Q6sB5nuzbuDja01dVZVwYB3sXgKZSyBrGZQVmOciEuHOXyAgxJrahBCNjueFe48J0DINspZC5lLIXAab5prnAsIhsa8J+1YDoEUv8A9yTR2Vx2DfGtj9k6ljz89QXmSea9oGOl8BrQaCXyB8PAmWvAQX/sk1tQghxGkcCnel1AjgJcAXeFNr/expzycB7wCR9m0e1lrPc3Ktxw8Gse3NV+9J5rGCbBPyxwN/4dPmcd8AaNnb3rofaII/6DwXpa4ohewVZv+7fzK3K8vMc7GdoMd4+5vKQAhvdurPbvoclvwbUq+HqFbnd3whhDgHStfRb62U8gW2AsOBbGAFMEFrnVFtm6nAGq31a0qpzsA8rXXy2fablpamV65cWc/ya1GcC3uWmyDOWmb6vHUVKB+I73qyG6fVAAiLq3kfxwpNazxzqelq2bsKbBWAMt0/rQZC8kCzn9CYs9dTkA2v9IG2F8O495z+coUQjYdSapXWOq2u7RxpufcFtmutd9p3PBO4Asioto0GjjeJI4B951auk4VGQ8dR5gvgWJFpaWctM2G96h34+XXzXHTbk0EfFGla/7uXwP519jcEX2iRCv3uMIGe1A+CI8+tnogEGHw/LHwGdv4ArYc68cUKIcSZHGm5jwVGaK1vtd+fCFygtb6r2jbNgW+AKCAUuFhrveps+3Vpy70uleWwf+3Jln31k5/Hu3KOd7EkXgCBYfU/ZkUZvNoX/EPgt4vB17/++xRCNDrObLnXtMbc6e8IE4DpWut/KqX6A+8ppbpqrW2nFXUbcBtAUlKSA4d2Eb8A0/+e2Be41wxbzNkEpfnQshf4Bzv/mP5BMOLvMPN6WPEW9Put848hhBB2joxzzwYSq91P4Mxul1uAWQBa62VAEHBGR7TWeqrWOk1rnRYbG3t+FbuCjw/EdzF96K4I9uM6XAZtLoLv/wbFh113HCFEo+dIuK8A2imlUpRSAcB4YO5p22QBwwCUUp0w4Z7jzEK9glIw4lmoKIbvnrK6GiGEF6sz3LXWlcBdwNfAJmCW1nqjUuoppdTl9s0eAKYopdYBM4BJuq7O/MYqtgNc8FtY/a4ZJy+EEC5Q5wlVV7H0hKrVygrgP72haWuY/LVp0QshhAMcPaHqeXPLeIOgCLj4CTOOfv3HVlcjhPBCEu5W6XG9mR7hm0fNBVNCCOFEEu5W8fGBy56HogOw+J9WVyOE8DIS7lZKSIPUG2DZq5C7w+pqhBBeRMLdasMeB99A+PrPVlcihPAiEu5WC4+HoX+ErV/B1m+srkYI4SUk3N1B39shuh189bCZ90YIIepJwt0d+AWYK1fzdsDPr1ldjRDCC0i4u4t2F0P7kbDoOSg8YHU1QggPJ+HuTi79K1SVw4InrK5ECOHhJNzdSXQb6H8XrJsBe36xuhohhAeTcHc3gx+A8OYw70Ezz7wQQpwHCXd3ExgGw582K0Wtfd/qaoQQHkrC3R11GwuJ/WDBk2Z1KCGEOEcS7u5IKbjsOSjJhUX/sLoaIYQHknB3V817QO9J8PMbcGiz1dUIITyMhLs7u+hR0wc//yGQha2EEOdAwt2dhUbDhX+BXYtg8xdWVyOE8CAS7u4ubTLEdTazRlaUuu44Bdlm6uG595hlAIUQHs3P6gJEHXz9YOQ/4J0xsPQ/MOQh5+376D7I+Aw2/s8s+XdcQTZcP8scWwjhkeR/rydISYfOV8LiF6HHBIhMPP99Hd0Pm+aaQM9aZh6L7woX/QU6XwWZP8Hnv4ev/2RWihJCeCQJd09xydOw9Wv49lG4dvq5/WzhwZOBnrkU0BDXxfTnd7kSYtqd3DamLRzeBstegZj20HeKM1+FEKKBSLh7isgkGHQf/PA3SLsFUgafffuiQ/ZA/xR2/wRoiO0IQ/9kAj22Q+0/O/wpyN0O8/9o5rtpc5FTX4oQwvWUtmiIXVpaml65cqUlx/ZYFaXwSl8IDIfbfzyzT7z48MkW+u6fQNtM67vL1SbQ4zo5fqxjhfDWpab//dZvz/5mIIRoMEqpVVrrtLq2c2i0jFJqhFJqi1Jqu1Lq4Vq2uU4plaGU2qiU+vBcCxYO8A820wIf2girppnHinNh5TR49wp4oR18cZ85UTr4D3DHMrjzF7jwT+cW7GDeQK6faRYS+fA6cxwhhMeos+WulPIFtgLDgWxgBTBBa51RbZt2wCzgIq31EaVUnNb60Nn2Ky3386S1CfL966BFT9j1I+gqaNra3kK/CuK7mCkMnGHPCpg+ChLSYOKnJuyFEJZxZsu9L7Bda71Ta10OzASuOG2bKcCrWusjAHUFu6gHpczQyMoyOLIbBv4ebl8Md6+GYY9Cs67OC3aAxD5wxauQucR8KpArZYXwCI6cUG0J7Kl2Pxu44LRt2gMopZYAvsATWuuvnFKhOFNcJ3hoJ/iHODfIa9P9WsjdZiYxi2kHg+51/TGFEPXiSLjXlB6nN9/8gHbAUCABWKyU6qq1PmW+WqXUbcBtAElJSedcrKgmILRhjzfkYTi81SwBGN0WOo1u2OMLIc6JI90y2UD1q2YSgH01bPOZ1rpCa70L2IIJ+1NoradqrdO01mmxsbHnW7Owgo8PXPkatOwFc6aYPn8hhNtyJNxXAO2UUilKqQBgPDD3tG0+BS4EUErFYLppdjqzUOEG/INh/IcQHAUzJkDhAasrEkLUos5w11pXAncBXwObgFla641KqaeUUpfbN/sayFVKZQDfAw9qrWXsnDcKbwYTZpoVomZMcO1kZkKI8yYXMYnzs/lLmHkDdL4Cxk4z3TZCCJdz6kVMQpyh4ygY/iRkfAo//N3qaoQQp5G5ZcT5G3AP5GyFH58z0xx0v9bqioQQdtJyF+dPKRj9L2g1ED67E/b8YnVFQgg7CXdRP34BMO59aNICZl4P+VmuP2bpEVjxllnAZM0Hrj+eEB5IumVE/YU0NSs3vXkxfDgObvnGTDzmTFUVsH0BrJsBW+ZDVTkoXygvhp43OPdYQngBabkL54htD9e9AzlbYPYtYKuq/z61hn1rYf7D8M+OMGM87F5i5rO/bREMfgD2rZE1X4WogbTchfO0uRAuew6+fAC+eRRG/O389nN0P6yfBWtnQM4m8A2ADiPNEoNtLwZff7PdsUJzMjdzGXQY4bzXIYQXkHAXztXnVrNM3/JXTWu+9yTHfq68BDZ/Ybpddv5gFhpJ6AujXoSuV5urYk+X0Ad8A820xxLuQpxCwl043yV/NcsfS3dsAAAUSklEQVT0ffkARKVA6yE1b2ezQdZS00LP+BTKiyAiySw00mO8WeLvbPyDIOkC2P2j81+DEB5Owl04n68fjH0b3roEZk2EWxeahbePy91hWujrPoKCLAgIg85XQuoESBpwble7JqfD989ASZ45sSuEACTchasERcD1H8F/LzLL9N34CexYCOtmQvYvoHyg9VAY9pi52jUg5PyOk5JuZjPa/RN0vrzOzYVoLCTchetEJZtZJN8ZAy+nmsdiO8Hwp6DbtWZsfH217AX+oabfXcJdiBMk3IVrJfWDa98xfetdx0LzHs5dPcrXH1r1h92LnbdPIbyAjHMXrtfxMrjkGWiR6pplAZMHQ85mKDzo/H0L4aEk3IXnS0k336X1LsQJEu7C8zXvAYERpt9dCAFIuAtv4OMLrQZIy12IaiTchXdISYe8nVCQbXUlQrgFCXfhHVIGm++7pPUuBEi4C28R1wWCm0rXjBB2Eu7CO/j4QPIgc1LVokXfhXAnEu7Ce6SkQ8EeOLLb6kqEsJyEu/Aex8e7y5BIISTchReJaQ9h8dLvLgQOhrtSaoRSaotSartS6uGzbDdWKaWVUmnOK1EIByllpiKQfnch6g53pZQv8CowEugMTFBKda5hu3DgHuBnZxcphMNS0qHooFkNSohGzJGWe19gu9Z6p9a6HJgJXFHDdk8DzwFlTqxPiHNzYrz7ImvrEMJijoR7S2BPtfvZ9sdOUEr1BBK11l84sTYhzl1UCkQkSr+7aPQcCfea5mg90aGplPIB/gU8UOeOlLpNKbVSKbUyJyfH8SqFcNSJfvfFZo1WIRopR8I9G0isdj8B2FftfjjQFfhBKbUb6AfMremkqtZ6qtY6TWudFhsbe/5VC3E2KelQmgeHMqyuRAjLOBLuK4B2SqkUpVQAMB6Ye/xJrXWB1jpGa52stU4GlgOXa61XuqRiIepyot9dxruLxqvOcNdaVwJ3AV8Dm4BZWuuNSqmnlFKyaKVwPxEJpu9d+t1FI+bQGqpa63nAvNMee6yWbYfWvywh6iklHTZ+CrYqM9+7EI2MXKEqvFNKOhwrgP3rrK5ECEtIuAvvlGzvd5euGdFISbgL7xQeDzEd5KSqaLQk3IX3SkmHzGVQVWF1JUI0OAl34b1SBkNFMexdbXUlQjQ4CXfhvU70u0vXjGh8JNyF9wppCvHdpN9dNEoS7sK7paTDnl+gQiYrFY2LhLvwbimDobIMsldYXYkQDUrCXXi3VgNA+ch4d9HoSLgL7xYUAc1Tpd9dNDoS7sL7pQyG7JVQXmx1JUI0GAl34f1S0sFWAVnLra5EiAYj4S68X2I/8PGTfnfRqEi4C+8XGAYt06TfXTQqEu6icUgZDPvWQtlRqysRokFIuIvGISUddBVkLbO6EiEahIS7aBwS+oJvoHTNiEZDwl00Dv5BkNhXwl00GhLuovFISYcD66Ekz+pKhHA5CXfReKSkAxoyl1hdiRAuJ+EuGo8WvcA/RLpmvM2B9fDd01BZbnUlbsXP6gKEaDB+AZDUH3bJxUxeo/gwfDgOju4FWyUMf9LqityGtNxF45IyGHI2QdEhqysR9WWrgjlTTMC3uwSWvAQ7vre6KrfhULgrpUYopbYopbYrpR6u4fn7lVIZSqlflVLfKaVaOb9UIZwgOd18l6kIPN+Pz8OOhXDZc3DtOxDTHv53uwl7UXe4K6V8gVeBkUBnYIJSqvNpm60B0rTW3YHZwHPOLlQIp2jeAwKbeHa/e2k+bFsA3/8N3r0CPrsTbDarq2pY2xfAD89CjwnQ6yYICIGxb5nfzad3gNZWV2g5R/rc+wLbtdY7AZRSM4ErgIzjG2itq38WWg7c6MwihXAaXz+zgIen9LvbbJC7Dfb8bJYLzF4BOZvNc8oHmraGnT9AkwS48E+WltpgCrLhkykQ1xlGvQhKmcebdYNLnoH5D8LPr0O/O6yt02KOhHtLYE+1+9nABWfZ/hZgfn2KEsKlUtJh61dQsBciWlpdzamOFZq557NX2MP8FygrMM8FR0FCH+g21lxx27IXBISZlvuiZ6FFT+gwwtr6Xa2yHGbdBFUVcN27psVeXd8ppqvm28fMm3jzHtbU6QYcCXdVw2M1fuZRSt0IpAFDann+NuA2gKSkJAdLFMLJkgeb77sXQ4/x1tWhNeTtPBnie36BQxmgbYCC2I7Q+UpzZW3iBRDd9mQrtbpR/4SDG2DObXDb9xDdpsFfSoP59lHYu9Lex972zOeVgitehdcGwOxb4PZFEBDa8HW6AUfCPRtIrHY/Adh3+kZKqYuBR4AhWutjNe1Iaz0VmAqQlpYmnWLCGvFdTSt4148NG+7lJbBvtQnx44FekmueC2wCCWnQcTQk9jFTFAdHOrZf/2AY9z68MQQ+uhFu+dZMc+xtNsyxd7f8DrpcWft2odFw9VRzPmL+H+GKVxquRjfiSLivANoppVKAvcB44PrqGyilegJvACO01jLGTLg3Hx9IHtSw/e6rpsOXfzArQoFphbcfYbpZEi+A2A7g43v++49MMicU378G5t4NY9+uuZXvqXK2mteV0BcudmAse+shMPh+WPxPaHMhdL3G9TW6mTrDXWtdqZS6C/ga8AXe1lpvVEo9BazUWs8FngfCgI+V+YPK0lpf7sK6haiflCGw6XM4shuikl17rK3fwBf3me6g/neaQA9p6vzjtLkIhj0GC56Alr1hwF3OP4YVyoth1m/ALxCunW4uRnPE0D+ZT2ef32s+CUU1rhHaDl2hqrWeB8w77bHHqt2+2Ml1CeFax/vdd/3o2nDfvw4+nmS6gsZ/6PrukoH3wt7V5oRi8+72+XQ8mNbmjTFnM0ycc24nwH394Zo34fXB8MmtcPN8M1qqkZArVEXjFNsBQuNc2zVTsNdcGh8cBdfPaph+cKXgyv8z3T4fTzLDBj3Zqmnw60emFd7monP/+ahkGP0vc35j0bNOL8+dSbiLxkkpMxXBrh9dc8FL2VH48Do4VgQ3zIImzZ1/jNoEhsP4D8ywwY8mQkVZwx3bmfatMSdE2wyD9AfPfz/dxkLqjfDjC55zfYMTSLiLxislHYoOQO525+63qsK0mg9tguvegfguzt2/I2LawVWvm9E58+sRjFYpPWL62UPj4Or/mpPg9THyH2aI6JwpUJzrnBrdnIS7aLxO9Lsvct4+tYYvH4Ad35nugLbDnLfvc9VpNAz+A6x+14zW8RQ2G/zvt3B0v3lzDI2u/z4Dw8wIopJcmHtXo5ieQMJdNF5NW5vL9p35UX3Jv2H1OzDofuh9k/P2e74u/LPp1pj3oLny1RMs+be5gvjSv5mx/87SvIcZRrllHqx403n7dVMS7qLxOt7vvnuxcybe2jDHDEPscjVc9Gj99+cMPr5mxEh4c9P/7u5THe9aDAufNr/DvlOcv/9+d0Db4fD1I3Bwo/P370Yk3EXjljzYfFTP2VS//WQtN10Jif3gytfq30fsTCFNzRWspXnw8c3mnIA7KjwAsyebkT6Xv+yai7CUMv8+QRHmWOUlzj+Gm3Cjv0AhLJBSbbz7+crdATMmmDHY4z8E/yDn1OZMzbvDmJch8yf49nGrqzlTVaV54ykvMhOCBYa77lhhsXD1G2bs/Nd/dt1xLCbhLhq3yCQzFvp8+91L8uCDa83tG2Y75+Sfq/QYB31vh+WvwvrZVldzqu+ehKylMOYliOvk+uO1uQgG/t6Mo8+Y6/rjWUDCXYiUdNj9k1m27VxUlMHM682FQhNmeMZsjJf+1awj+9ldcGCD1dUYm76ApS9D2i3Q/bqGO+6FfzGLps+9C/L31L29h5FwFyI5HY4VwIFfHf8Zmw0++x1kLYOrXoOkfq6rz5l8/c10uUERZgbJ0iPW1pO3Ez79nZmLfsTfG/bYfgFmsjVblZkuuaqyYY/vYhLuQpxPv/v3z8CGT2DY454342B4PIx7z3zimHObdUv0VZSaC5WUMm84foENX0PT1mY1p6ylsPiFhj++C0m4CxHezCyu7Gi/++p3zVSyvW6CQfe5tjZXSewLI5+Fbd/Aon9YU8P8h+DAejP3upUzNvYYB93Hm99D5lLr6nAyCXchwPS7Zy6te5jgjoVmCtk2F5kVkDx5zvS0WyD1BjOh1pYGXhlzzQfmTXLwA9D+0oY9dk1GvWBOrH8yxfquKieRcBcCzHj3imIzWVVtDm4063fGdjTdCL7+DVefKyhl3qCa9zDdM7k7Gua4BzbAl/eb3/lQNxmKGBgO17xl5hqae49XTE8g4S4EnDq/e02O7ocPrjPrcd4wC4KaNFxtrnR8iT4fP5h5g5nF0pXKCmDWRAiKNHO9uNP86i17mXMom+Z61lw8tZBwFwLM+PT4rjWH+7EimDHOfFy//iOISGj4+lzp+BJ9h7e4blKtY0Wmf/3T38GRTLh2GoTFOf849dX/LtPl9tXDZlbP+rLZzLUQh7dD1s+m+2vN+w0yz48bvW0KYbGUdFj5NlQeOzlyw1YFn9xigmnCR6YLwxudsUTf3ee+j5I8yNtlhjcesX8/fr+42pw2lzwDrQY4rXSn8vGBK1+H1wbA7Ftgynfm0w2YN72yAjNdRUmemc6hJPe0r7xTv5fmga5hNNKAu507KVoNJNyFOC55MCz/P8heYRbQ1tq04LZ+BZe9AO0vsbpC16q+RF+z7maR6eq0NvO/nB7cx++XFZy6fXgLaJpifm9RKWbYYWwHa+a3Pxfh8WYu/A/GwhtDQPmcDGpbLWPhffwhJNrM4xMSDXEd7fftX8FNT30+NNblL0PCXYjjWg0w/5F3LTbhvvw1+GWq+ajuihkK3c3xJfr+uwVm3wxDHob8TLOIeN5O872i2kRbyhciE01odx1rvje1h3hU8skWrydqNxxG/MP0v4c0hZALzh7WgeFuN3JKaYvOCqelpemVKz1kfmnReEwdCn7B0P93ZorcTmPMyBh3muXR1Q5vg/8OM1ft+gaawD7e8m6acvJ+ZJLnjxjyQEqpVVrrOvt0pOUuRHXJg02L/ZMppu/56qmNK9jBLNF3769QXmzmgW9sr99LyL+aENWlDAFbhRnJMWGmZ3ct1EdwpJnCWILdY0nLXYjqUgbDBXdAn1vNvN9CeCiH3paVUiOUUluUUtuVUg/X8HygUuoj+/M/K6WSnV2oEA3CL9DMuRLT1upKhKiXOsNdKeULvAqMBDoDE5RSnU/b7BbgiNa6LfAvwKKZiIQQQoBjLfe+wHat9U6tdTkwE7jitG2uAN6x354NDFPKzcYFCSFEI+JIuLcEqi9Tkm1/rMZttNaVQAHgxuuNCSGEd3Mk3GtqgZ8+ON6RbVBK3aaUWqmUWpmTk+NIfUIIIc6DI+GeDSRWu58A7KttG6WUHxAB5J2+I631VK11mtY6LTZWRiIIIYSrOBLuK4B2SqkUpVQAMB44fbnwucBN9ttjgYXaqktfhRBC1D3OXWtdqZS6C/ga8AXe1lpvVEo9BazUWs8F3gLeU0ptx7TYx7uyaCGEEGfn0EVMWut5wLzTHnus2u0y4FrnliaEEOJ8WTZxmFIqB8g8zx+PAQ47sRxX86R6PalW8Kx6PalW8Kx6PalWqF+9rbTWdZ60tCzc60MptdKRWdHchSfV60m1gmfV60m1gmfV60m1QsPUK7MCCSGEF5JwF0IIL+Sp4T7V6gLOkSfV60m1gmfV60m1gmfV60m1QgPU65F97kIIIc7OU1vuQgghzsLjwr2uueXdhVIqUSn1vVJqk1Jqo1Lq91bX5AillK9Sao1S6gurazkbpVSkUmq2Umqz/Xfc3+qazkYpdZ/972CDUmqGUirI6pqqU0q9rZQ6pJTaUO2xpkqpb5VS2+zfo6ys8bhaan3e/rfwq1Lqf0qpSCtrrK6meqs99wellFZKxTj7uB4V7g7OLe8uKoEHtNadgH7AnW5ca3W/BzZZXYQDXgK+0lp3BHrgxjUrpVoC9wBpWuuumCu93e0q7unAiNMeexj4TmvdDvjOft8dTOfMWr8FumqtuwNbgT81dFFnMZ0z60UplQgMB7JccVCPCnccm1veLWit92utV9tvF2LC5/Spkt2KUioBGAW8aXUtZ6OUagKkY6a9QGtdrrXOt7aqOvkBwfaJ9UI4c/I9S2mtf+TMyf6qr9PwDnBlgxZVi5pq1Vp/Y59uHGA5ZoJDt1DL7xbMwkYPUcMMus7gaeHuyNzybse+7GBP4GdrK6nTvzF/bDarC6lDayAHmGbvQnpTKRVqdVG10VrvBV7AtND2AwVa62+srcoh8Vrr/WAaK0CcxfU4ajIw3+oizkYpdTmwV2u9zlXH8LRwd2jeeHeilAoDPgHu1Voftbqe2iilRgOHtNarrK7FAX5AL+A1rXVPoBj36TI4g72v+gogBWgBhCqlbrS2Ku+klHoE0yX6gdW11EYpFQI8AjxW17b14Wnh7sjc8m5DKeWPCfYPtNZzrK6nDgOBy5VSuzHdXRcppd63tqRaZQPZWuvjn4RmY8LeXV0M7NJa52itK4A5wACLa3LEQaVUcwD790MW13NWSqmbgNHADW4+5XgbzBv9Ovv/twRgtVKqmTMP4mnh7sjc8m7BvobsW8AmrfWLVtdTF631n7TWCVrrZMzvdaHW2i1bl1rrA8AepVQH+0PDgAwLS6pLFtBPKRVi/7sYhhufAK6m+joNNwGfWVjLWSmlRgB/BC7XWpdYXc/ZaK3Xa63jtNbJ9v9v2UAv+9+103hUuNtPmByfW34TMEtrvdHaqmo1EJiIaQGvtX9dZnVRXuRu4AOl1K9AKvA3i+uplf0TxmxgNbAe8//Ora6oVErNAJYBHZRS2UqpW4BngeFKqW2YUR3PWlnjcbXU+goQDnxr/7/2uqVFVlNLva4/rnt/ehFCCHE+PKrlLoQQwjES7kII4YUk3IUQwgtJuAshhBeScBdCCC8k4S6EEF5Iwl0IIbyQhLsQQnih/wd/+O7+nXnyiQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "import matplotlib.pyplot as plt\n", "\n", diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index 2356bec4..16b811e6 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -164,7 +164,9 @@ def linear_inv_state_estimate(results: List[ExperimentResult], def compressed_sensing_state_estimate(results: List[ExperimentResult], qubits: List[int]) -> np.ndarray: """ - Estimate a quantum state using compressed sensing + Estimate a quantum state using compressed sensing. + + Specifically by constrained trace minimization a.k.a. the matrix Dantzig selector. [FLAMMIA] Quantum Tomography via Compressed Sensing: Error Bounds, Sample Complexity, and Efficient Estimators Steven T. Flammia et. al. @@ -176,23 +178,24 @@ def compressed_sensing_state_estimate(results: List[ExperimentResult], which qubits will be kron'ed together. :return: A point estimate of the quantum state rho. """ - pauli_num = len(results) - qubit_num = len(qubits) - d = 2 ** qubit_num + num_pauli = len(results) + num_qubit = len(qubits) + d = 2 ** num_qubit pauli_list = [] expectation_list = [] - y = np.zeros((m,1)) + y = np.zeros((num_pauli, 1)) - for i in range(pauli_num): + for i in range(num_pauli): #Convert the Pauli term into a tensor r = results[i] - p_tensor = pauli2matrix(r.setting.out_operator, qubits) + p_tensor = pauli2matrix(r.setting.observable, qubits) e = r.expectation #A[i] = e * scale_factor pauli_list.append(p_tensor) expectation_list.append(e) - s = cp.Variable((d,d),complex = True) + # The objective and constraint in terms of Eqn (3) and (34) in [FLAMMIA + s = cp.Variable((d, d), complex = True) obj = cp.Minimize(cp.norm(s, 'nuc')) constraints = [cp.trace(s) == 1] @@ -221,28 +224,32 @@ def lasso_state_estimate(results: List[ExperimentResult], which qubits will be kron'ed together. :return: A point estimate of the quantum state rho. """ - pauli_num = len(results) - qubit_num = len(qubits) - d = 2 ** qubit_num + num_pauli = len(results) + num_qubit = len(qubits) + d = 2 ** num_qubit pauli_list = [] expectation_list = [] - y = np.zeros((m,1)) - - mu = 4 * pauli_num / np.sqrt(1000 * pauli_num) - - for i in range(m): + y = np.zeros((num_pauli, 1)) + + + for i in range(num_pauli): #Convert the Pauli term into a tensor r = results[i] - p_tensor = pauli2matrix(r.setting.out_operator, qubits) + p_tensor = pauli2matrix(r.setting.observable, qubits) e = r.expectation y[i] = e pauli_list.append(p_tensor) expectation_list.append(e) - - x = cp.Variable((d,d),complex = True) - A = cp.vstack([cp.trace(cp.matmul(pauli_list[i], x)) * np.sqrt(d / m) for i in range(m)]) - - #Minimize trace norm + + x = cp.Variable((d, d), complex=True) + A = cp.vstack([cp.trace(cp.matmul(pauli_list[i], x)) * np.sqrt(d / num_pauli) + for i in range(num_pauli)]) + + # look at section V. A of [FLAMMIA] for more information related to mu + mu = 4 * num_pauli / np.sqrt(1000 * num_pauli) + + # The equation below is Eqn. (4) and Eqn. (35) from [FLAMMIA] + # Minimize trace norm obj = cp.Minimize(0.5 * cp.norm((A - y), 2) + mu * cp.norm(x, 'nuc')) constraints = [cp.trace(x) == 1] From 30e86856c3ad5f33aa567385cf62d7051672cebe Mon Sep 17 00:00:00 2001 From: Joshua Combes Date: Wed, 26 Jun 2019 12:43:25 +1000 Subject: [PATCH 12/14] list comprehension --- forest/benchmarking/tomography.py | 29 +++++++++++------------------ 1 file changed, 11 insertions(+), 18 deletions(-) diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index 16b811e6..a01b42a0 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -224,29 +224,22 @@ def lasso_state_estimate(results: List[ExperimentResult], which qubits will be kron'ed together. :return: A point estimate of the quantum state rho. """ - num_pauli = len(results) - num_qubit = len(qubits) - d = 2 ** num_qubit - pauli_list = [] - expectation_list = [] - y = np.zeros((num_pauli, 1)) - + n_pauli = len(results) + n_qubits = len(qubits) + d = 2 ** n_qubits - for i in range(num_pauli): - #Convert the Pauli term into a tensor - r = results[i] - p_tensor = pauli2matrix(r.setting.observable, qubits) - e = r.expectation - y[i] = e - pauli_list.append(p_tensor) - expectation_list.append(e) + # Convert the Pauli term into a matrix + pauli_list = [pauli2matrix(r.setting.observable, qubits) for r in results] + # shape of y is (num_pauli,1) + y = np.array([[r.expectation] for r in results]) x = cp.Variable((d, d), complex=True) - A = cp.vstack([cp.trace(cp.matmul(pauli_list[i], x)) * np.sqrt(d / num_pauli) - for i in range(num_pauli)]) + A = cp.vstack([cp.trace(cp.matmul(pauli_list[i], x)) * np.sqrt(d / n_pauli) + for i in range(n_pauli)]) # look at section V. A of [FLAMMIA] for more information related to mu - mu = 4 * num_pauli / np.sqrt(1000 * num_pauli) + num_experiments = (results[0].total_counts) * n_pauli + mu = 4 * n_pauli / np.sqrt(num_experiments) # The equation below is Eqn. (4) and Eqn. (35) from [FLAMMIA] # Minimize trace norm From 27d8456c9ea38f3043d0189f0daae0e368aef915 Mon Sep 17 00:00:00 2001 From: Joshua Combes Date: Wed, 26 Jun 2019 13:45:59 +1000 Subject: [PATCH 13/14] more code tidying --- forest/benchmarking/tomography.py | 67 +++++++++++++++---------------- 1 file changed, 32 insertions(+), 35 deletions(-) diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index a01b42a0..1bf66391 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -168,45 +168,39 @@ def compressed_sensing_state_estimate(results: List[ExperimentResult], Specifically by constrained trace minimization a.k.a. the matrix Dantzig selector. - [FLAMMIA] Quantum Tomography via Compressed Sensing: Error Bounds, Sample Complexity, and Efficient Estimators - Steven T. Flammia et. al. - (2012). - https://arxiv.org/pdf/1205.2300.pdf + [QTvCS] Quantum Tomography via Compressed Sensing: Error Bounds, Sample Complexity, and ... + Flammia et. al. + New J. Phys. 14, 095022 (2012) + https://dx.doi.org/10.1088/1367-2630/14/9/095022 + https://arxiv.org/pdf/1205.2300.pdf :param results: A tomographically complete list of results. :param qubits: All qubits that were tomographized. This specifies the order in which qubits will be kron'ed together. :return: A point estimate of the quantum state rho. """ - num_pauli = len(results) - num_qubit = len(qubits) - d = 2 ** num_qubit - pauli_list = [] - expectation_list = [] - y = np.zeros((num_pauli, 1)) - - for i in range(num_pauli): - #Convert the Pauli term into a tensor - r = results[i] - p_tensor = pauli2matrix(r.setting.observable, qubits) - e = r.expectation - #A[i] = e * scale_factor - pauli_list.append(p_tensor) - expectation_list.append(e) - - # The objective and constraint in terms of Eqn (3) and (34) in [FLAMMIA - s = cp.Variable((d, d), complex = True) - obj = cp.Minimize(cp.norm(s, 'nuc')) - constraints = [cp.trace(s) == 1] - - for i in range(len(pauli_list)): - trace_bool = cp.trace(cp.matmul(pauli_list[i], s)) - expectation_list[i] == 0 - constraints.append(trace_bool) + n_pauli = len(results) + n_qubits = len(qubits) + d = 2 ** n_qubits + + + # Convert the Pauli term into a matrix + pauli_list = [pauli2matrix(r.setting.observable, qubits) for r in results] + # a list of expectations + y_list = [r.expectation for r in results] + + # The objective and constraint in terms of Eqn (3) and (34) in [QTvCS] + x = cp.Variable((d, d), complex = True) + obj = cp.Minimize(cp.norm(x, 'nuc')) + # A[i] = y[i] * scale_factor + constraints = [cp.trace(cp.matmul(pauli_list[i], x)) - y_list[i] == 0 for i in range(n_pauli)] + constraints.insert(0, cp.trace(x) == 1) # Form and solve problem. prob = cp.Problem(obj, constraints) prob.solve() - rho = s.value + rho = x.value + return rho def lasso_state_estimate(results: List[ExperimentResult], @@ -214,10 +208,11 @@ def lasso_state_estimate(results: List[ExperimentResult], """ Estimate a quantum state using compressed sensing - [FLAMMIA] Quantum Tomography via Compressed Sensing: Error Bounds, Sample Complexity, and Efficient Estimators - Steven T. Flammia et. al. - (2012). - https://arxiv.org/pdf/1205.2300.pdf + [QTvCS] Quantum Tomography via Compressed Sensing: Error Bounds, Sample Complexity, and ... + Flammia et. al. + New J. Phys. 14, 095022 (2012) + https://dx.doi.org/10.1088/1367-2630/14/9/095022 + https://arxiv.org/pdf/1205.2300.pdf :param results: A tomographically complete list of results. :param qubits: All qubits that were tomographized. This specifies the order in @@ -230,18 +225,20 @@ def lasso_state_estimate(results: List[ExperimentResult], # Convert the Pauli term into a matrix pauli_list = [pauli2matrix(r.setting.observable, qubits) for r in results] + # shape of y is (num_pauli,1) y = np.array([[r.expectation] for r in results]) x = cp.Variable((d, d), complex=True) + # Eqn 1 of [QTvCS] A = cp.vstack([cp.trace(cp.matmul(pauli_list[i], x)) * np.sqrt(d / n_pauli) for i in range(n_pauli)]) - # look at section V. A of [FLAMMIA] for more information related to mu + # look at section V. A of [QTvCS] for more information related to mu num_experiments = (results[0].total_counts) * n_pauli mu = 4 * n_pauli / np.sqrt(num_experiments) - # The equation below is Eqn. (4) and Eqn. (35) from [FLAMMIA] + # The equation below is Eqn. (4) and Eqn. (35) from [QTvCS] # Minimize trace norm obj = cp.Minimize(0.5 * cp.norm((A - y), 2) + mu * cp.norm(x, 'nuc')) constraints = [cp.trace(x) == 1] From 1101793b4a960619cea215ea3f30dfe1a6b5aec8 Mon Sep 17 00:00:00 2001 From: Joshua Combes Date: Wed, 26 Jun 2019 14:20:10 +1000 Subject: [PATCH 14/14] doc strings update --- forest/benchmarking/tomography.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index 1bf66391..6ff7fea7 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -164,9 +164,9 @@ def linear_inv_state_estimate(results: List[ExperimentResult], def compressed_sensing_state_estimate(results: List[ExperimentResult], qubits: List[int]) -> np.ndarray: """ - Estimate a quantum state using compressed sensing. + Estimate a quantum state using compressed sensing via the matrix Dantzig selector. - Specifically by constrained trace minimization a.k.a. the matrix Dantzig selector. + The matrix Dantzig selector is constrained trace minimization. See [QTvCS] for more information. [QTvCS] Quantum Tomography via Compressed Sensing: Error Bounds, Sample Complexity, and ... Flammia et. al. @@ -183,7 +183,6 @@ def compressed_sensing_state_estimate(results: List[ExperimentResult], n_qubits = len(qubits) d = 2 ** n_qubits - # Convert the Pauli term into a matrix pauli_list = [pauli2matrix(r.setting.observable, qubits) for r in results] # a list of expectations @@ -206,7 +205,10 @@ def compressed_sensing_state_estimate(results: List[ExperimentResult], def lasso_state_estimate(results: List[ExperimentResult], qubits: List[int]) -> np.ndarray: """ - Estimate a quantum state using compressed sensing + Estimate a quantum state using compressed sensing via a matrix Lasso. + + For more information see https://en.wikipedia.org/wiki/Lasso_(statistics) and the reference + [QTvCS]. [QTvCS] Quantum Tomography via Compressed Sensing: Error Bounds, Sample Complexity, and ... Flammia et. al.