diff --git a/docs/tomography.rst b/docs/tomography.rst index c682f4cb..1a2b5c23 100644 --- a/docs/tomography.rst +++ b/docs/tomography.rst @@ -75,6 +75,14 @@ Finally, we analyze our data with one of the analysis routines:: [ 0.23 +0.027j 0.175-0.j 0.277-0.j -0.173+0.004j] [-0.203+0.01j -0.168+0.019j -0.173-0.004j 0.229-0.j ]] +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") State Tomography @@ -88,6 +96,7 @@ State Tomography linear_inv_state_estimate iterative_mle_state_estimate estimate_variance + tomographize Process Tomography diff --git a/examples/tomography_debugger.ipynb b/examples/tomography_debugger.ipynb new file mode 100644 index 00000000..4de5bd5b --- /dev/null +++ b/examples/tomography_debugger.ipynb @@ -0,0 +1,355 @@ +{ + "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": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import time\n", + "\n", + "from pyquil import Program, get_qc\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": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "H 0\n", + "H 1\n", + "CZ 0 1\n", + "RY(-pi/2) 0\n", + "X 1\n", + "CNOT 1 0\n", + "\n" + ] + } + ], + "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": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Linear tomography took 0.460453s\n", + "Recovered density matrix:\n", + "\n", + "[[ 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" + ] + } + ], + "source": [ + "qc = get_qc('%dq-qvm' % len(qubits))\n", + "\n", + "\n", + "start_linear = time.time()\n", + "m = 10\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", + "print(\"Recovered density matrix:\\n\")\n", + "print(rho_linear)\n", + "\n", + "start_compressed = time.time()\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)\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)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Compare results to true output obtained using wavefunction simulator" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "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", + "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": 5, + "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": 6, + "metadata": {}, + "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", + "print(\"Compressed norm:\")\n", + "print(np.linalg.norm(rho_compressed - rho_true))\n", + "print(\"Lasso norm:\")\n", + "print(np.linalg.norm(rho_lasso - rho_true))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot graph of results for various measurement values" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "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", + "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(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", + " 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 +} diff --git a/forest/benchmarking/tomography.py b/forest/benchmarking/tomography.py index 0728ecef..6ff7fea7 100644 --- a/forest/benchmarking/tomography.py +++ b/forest/benchmarking/tomography.py @@ -3,20 +3,23 @@ from operator import mul from typing import Callable, Tuple, List, Sequence, Iterator import warnings +import random import numpy as np +import cvxpy as cp from scipy.linalg import logm, pinv from pyquil import Program from pyquil.unitary_tools import lifted_pauli as pauli2matrix, lifted_state_operator as state2matrix +from forest.benchmarking.compilation import basic_compile import forest.benchmarking.distance_measures as dm from forest.benchmarking.utils import all_traceless_pauli_terms from forest.benchmarking.operator_tools import vec, unvec, proj_choi_to_physical from forest.benchmarking.operator_tools.project_state_matrix import project_state_matrix_to_physical from forest.benchmarking.observable_estimation import ExperimentSetting, ObservablesExperiment, \ ExperimentResult, SIC0, SIC1, SIC2, SIC3, plusX, minusX, plusY, minusY, plusZ, minusZ, \ - TensorProductState, zeros_state + TensorProductState, zeros_state, estimate_observables, group_settings MAXITER = "maxiter" OPTIMAL = "optimal" @@ -49,7 +52,7 @@ def generate_state_tomography_experiment(program: Program, qubits: List[int]): To collect data, try:: from forest.benchmarking.operator_estimation import estimate_observables - results = list(measure_observables(qc=qc, experiment, num_shots=100_000)) + results = list(estimate_observables(qc=qc, experiment, num_shots=100_000)) :param program: The program to prepare a state to tomographize :param qubits: The qubits to tomographize @@ -158,6 +161,97 @@ def linear_inv_state_estimate(results: List[ExperimentResult], dim = 2**len(qubits) return unvec(rho) + np.eye(dim)/dim +def compressed_sensing_state_estimate(results: List[ExperimentResult], + qubits: List[int]) -> np.ndarray: + """ + Estimate a quantum state using compressed sensing via 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. + 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. + """ + 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 = x.value + + return rho + +def lasso_state_estimate(results: List[ExperimentResult], + qubits: List[int]) -> np.ndarray: + """ + 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. + 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. + """ + 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] + + # 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 [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 [QTvCS] + # 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) \ @@ -612,3 +706,50 @@ 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) + +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", "lasso" + :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 = ObservablesExperiment(settings=exp_list, program=program) + + #Group experiments if possible to minimize QPU runs + input_exp = group_settings(input_exp) + + #NOTE: Change qvm depending on whether we are simulating qvm + qc.compiler.quil_to_native_quil = basic_compile + + results = list(estimate_observables(qc=qc, obs_expt=input_exp, num_shots=num_shots)) + + if t_type == 'compressed_sensing': + return compressed_sensing_state_estimate(results=results, qubits=qubits) + elif t_type == 'mle': + 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) +