{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# MPS techniques for Boson Sampling\n", "\n", "In this notebook, we explain how to use the MPS (Matrix Product State) backend to simulate a linear circuit. MPS simulation is based on a type of tensor network simulation, which gives an approximation of the output states [[1]](#References) [[2]](#References). It does the computation on each component of the circuits one-by-one, and not on the whole unitary. The states are represented by tensors, which are then updated at each component. These tensors can be seen as a big set of matrices, and the approximation is done by chosing the dimension of these matrices, called the *bond dimension*. For this example, we simulate a simple boson sampling problem, with 6 modes and 3 photons [[3]](#References)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import perceval as pcvl\n", "import matplotlib.pyplot as plt" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Definition of the circuit\n", "\n", "Just like in the Boson Sampling notebook, we generate a Haar-random unitary and its decomposition in a circuit :" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "6 modes triangular Boson Sampler :\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Φ=4.209463\n", "\n", "\n", "Φ=1.323728\n", "\n", "\n", "Φ=1.854704\n", "\n", "\n", "Φ=5.076331\n", "\n", "\n", "Φ=1.600842\n", "\n", "\n", "Φ=4.730231\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.733616\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.322967\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.287612\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.955033\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=6.096243\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.478647\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.824055\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.84601\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.370167\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=2.522208\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.464975\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.272642\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.315043\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.406417\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.707762\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=2.772472\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.322713\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.10431\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.931496\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=3.83593\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.799656\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.282961\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.927045\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.760345\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.402043\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.512277\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.878572\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.821695\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.891583\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.355083\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "0\n", "1\n", "2\n", "3\n", "4\n", "5\n", "0\n", "1\n", "2\n", "3\n", "4\n", "5\n", "" ], "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = 6\n", "unitary = pcvl.Matrix.random_unitary(m)\n", "mzi = (pcvl.BS() // (0, pcvl.PS(phi=pcvl.Parameter(\"φ_a\")))\n", " // pcvl.BS() // (1, pcvl.PS(phi=pcvl.Parameter(\"φ_b\"))))\n", "linear_circuit = pcvl.Circuit.decomposition(unitary, mzi,\n", " phase_shifter_fn=pcvl.PS,\n", " shape=\"triangle\")\n", "\n", "print(m, \" modes triangular Boson Sampler :\")\n", "pcvl.pdisplay(linear_circuit, compact = True)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## MPS simulation\n", "\n", "Let us now define the MPS simulation, using the MPS backend in Perceval." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "mps = pcvl.MPSBackend()\n", "mps.set_circuit(linear_circuit)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We now choose the size of the matrices (the bond dimension) for our simulation." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "chi = 8\n", "mps.set_cutoff(chi)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "And finally, we get the output probability distribution from a given input state with 3 photons." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
state probability
|0,0,1,0,2,0> 0.065989
|1,0,1,0,0,1> 0.058044
|2,1,0,0,0,0> 0.057078
|0,1,1,0,1,0> 0.050927
|0,0,1,1,1,0> 0.046433
|1,1,1,0,0,0> 0.040445
|1,0,0,0,2,0> 0.040187
|0,0,1,0,0,2> 0.036448
|2,0,0,1,0,0> 0.036259
|0,1,1,1,0,0> 0.033582
|1,0,2,0,0,0> 0.032595
|1,1,0,0,1,0> 0.031704
|1,0,1,0,1,0> 0.026459
|3,0,0,0,0,0> 0.025905
|0,1,0,2,0,0> 0.022392
|0,0,2,0,0,1> 0.022339
|0,0,1,1,0,1> 0.021238
|0,1,1,0,0,1> 0.020432
|0,0,0,2,1,0> 0.019627
|0,1,0,1,1,0> 0.019612
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n = 3\n", "input_state = pcvl.BasicState([1]*n + [0]*(m-n))\n", "mps.set_input_state(input_state)\n", "probs = mps.prob_distribution()\n", "pcvl.pdisplay(probs, max_v=20)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Certification of the MPS method :\n", "\n", "As we make an approximation by chosing the bond dimension, we have to check when does this approximation becomes good enough. Unfortunately, there is no formula giving the minimal size for a given approximation error. What we can do though is to compute the *Total Variance Distance* (TVD) between an ideal simulation of Boson Sampling, and an approximated one. To compute the ideal one, we can for instance use the *SLOS* backend on Perceval :" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
state probability
|0,0,1,0,2,0> 0.065989
|1,0,1,0,0,1> 0.058044
|2,1,0,0,0,0> 0.057078
|0,1,1,0,1,0> 0.050927
|0,0,1,1,1,0> 0.046433
|1,1,1,0,0,0> 0.040445
|1,0,0,0,2,0> 0.040187
|0,0,1,0,0,2> 0.036448
|2,0,0,1,0,0> 0.036259
|0,1,1,1,0,0> 0.033582
|1,0,2,0,0,0> 0.032595
|1,1,0,0,1,0> 0.031704
|1,0,1,0,1,0> 0.026459
|3,0,0,0,0,0> 0.025905
|0,1,0,2,0,0> 0.022392
|0,0,2,0,0,1> 0.022339
|0,0,1,1,0,1> 0.021238
|0,1,1,0,0,1> 0.020432
|0,0,0,2,1,0> 0.019627
|0,1,0,1,1,0> 0.019612
" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "slos = pcvl.SLOSBackend()\n", "slos.set_circuit(pcvl.Unitary(unitary))\n", "slos.set_input_state(input_state)\n", "probs_slos = slos.prob_distribution()\n", "pcvl.pdisplay(probs_slos, max_v=20)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We also have to define the TVD function." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def tvd(probs1, probs2):\n", " tvd = 0\n", " for state, prob in probs1.items():\n", " tvd += abs(prob - probs2[state])\n", " return tvd\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now we compute the TVD between the two simulations for different bond dimensions, going from 1 to 10." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAHHCAYAAABuoFaQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAA9hAAAPYQGoP6dpAABo0klEQVR4nO3dd1wT9/8H8FcSIEG2IFMExIETFJSiolhpcVerVetAadUW0Wqp/VW0rg6tddSvdda9a+uuqyp1S8WiOFrEBeJiOVgqI7nfH5bUyFbIAXk9H497POTyubtXLubyzt3nPpEIgiCAiIiIiCqcVOwARERERLqChRcRERGRlrDwIiIiItISFl5EREREWsLCi4iIiEhLWHgRERERaQkLLyIiIiItYeFFREREpCUsvIiIiIi0pNoWXhKJBNOmTSvXda5ZswYSiQTx8fHlut7yNnv2bNStWxcymQweHh5lXv7o0aOQSCTYunVr+YerQqZNmwaJRILU1FSxo1A15OzsjGHDhokdo0rx8/ND06ZNy2198fHxkEgkmDNnTolt848HL3r5Ncw/dh49erTcMpZF/vNZs2aNel5huSuKn58f/Pz81H9r+7Nk2LBhcHZ21sq2XkeFFl75hUpR059//lmRm39lM2bMwM6dO8WO8UoOHjyI//u//0Pbtm2xevVqzJgxo8i2mzZtwvz587UX7gV+fn7F/t/In3r27AmJRIIvv/yyyHVdu3YNEokEoaGhAP470ORPNWrUQJ06ddCjRw+sXr0a2dnZ2nqaBYi5z7UlPj4eQUFBcHV1hUKhgK2tLdq3b4+pU6dqtCvth+iDBw/w+eefo2HDhlAoFKhZsyYCAgKwZ8+eQtunpKRg7NixcHNzg6GhIaytrdG6dWt88cUXyMzMLJfnqE35/4+HDx9e6OOTJk1St3nxS8KwYcM03gempqZwd3fH3LlzC7wHTp48iS5dusDBwQEKhUL9ftm0aVOFPjddUBXf8/fu3cO0adMQHR0tdpQCKnO20tLTxka++uoruLi4FJhfr149bWy+zGbMmIG+ffuiV69eGvOHDBmCAQMGQC6XixOsFP744w9IpVKsXLkSBgYGxbbdtGkTLl++jHHjxmkn3AsmTZqk8UFy9uxZLFiwABMnTkSjRo3U85s3b45r165h8+bN+OabbwpdV/6Hw+DBgzXmL1myBMbGxsjOzsbdu3fx+++/44MPPsD8+fOxZ88eODo6VsAzK56Y+1wbrl+/jlatWsHQ0BAffPABnJ2dcf/+fZw7dw6zZs3C9OnTy7S+2NhYdOrUCSkpKQgKCoKXlxceP36MjRs3okePHhg/fjxmz56tbv/w4UN4eXkhPT0dH3zwAdzc3PDgwQNcvHgRS5YsQXBwMIyNjcv7aVc4hUKBbdu2YfHixQXe15s3b4ZCocCzZ88KLCeXy7FixQoAwOPHj7Ft2zaMHz8eZ8+exc8//wwA+PXXX9G/f394eHhg7NixsLCwQFxcHI4fP47ly5dj4MCBFf8Eq4Avv/wSEyZMKLZN+/bt8fTpU43XSOz3fGlyv+zevXuYPn06nJ2dy3TV5ODBg2VMV3bFZVu+fDlUKlWFZ3hdWim8unTpAi8vL21sqkLJZDLIZDKxYxQrOTkZhoaGJRZdYnvrrbc0/lYoFFiwYAHeeustjVPVADBo0CBMnjwZf/75J954440C69q8eTPc3NzQsmVLjfl9+/aFlZWV+u8pU6Zg48aNCAwMxHvvvVdpz7hWZT/88AMyMzMRHR0NJycnjceSk5PLtK7c3Fz07dsXjx49wvHjx+Ht7a1+7NNPP8WgQYMwZ84ceHl5oX///gCAlStXIiEhAadOnUKbNm001peenl7p3xdF6dy5M3bv3o39+/fjnXfeUc8/ffo04uLi0KdPH2zbtq3Acnp6ehpfSEaNGgVvb29s2bIF8+bNg729PaZNm4bGjRvjzz//LLB/yvqalVVWVhaMjIwqdBvlRU9PD3p6xX9kSqVSKBQKLSUqndLkfl1PnjxBjRo1RH9/6evri7r90hK9j1dubi5q1qyJoKCgAo+lp6dDoVBg/Pjx6nnJycn48MMPYWNjA4VCAXd3d6xdu7bE7RR17ffl698SiQRZWVlYu3at+hR9/jX8ovp4LV68GE2aNIFcLoe9vT1CQkLw+PFjjTb5l1X++ecfdOzYETVq1ICDgwO+//77ErMDQF5eHr7++mu4urpCLpfD2dkZEydO1LhkIJFIsHr1amRlZamzv3it/+U8e/fuxa1bt9RtX94/KpUK3377LWrXrg2FQoFOnTrh+vXrBdZ15swZdO7cGWZmZqhRowY6dOiAU6dOlep5lcagQYMAoNDLHlFRUYiNjVW3Kc26hg8fjjNnzuDQoUOlWiY1NRX9+vWDqakpLC0tMXbs2ELPLmzYsAGenp4wNDREzZo1MWDAANy+fVv9eFH7XBAEWFlZqS+VAs/3vbm5OWQymcb/pVmzZkFPT0/jktmVK1fQt29f1KxZEwqFAl5eXti9e3eBfI8fP8a4cePg6OgIuVyOevXqYdasWRrfEF/s8/LTTz+p/7+1atUKZ8+eLXFf3bhxA7Vr1y5QdAGAtbV1icu/aNu2bbh8+TImTJigUXQBz78ELVu2DObm5hp9OW/cuAGZTFZogW5qalrih+KtW7cwatQoNGzYEIaGhrC0tMR7771X4D2ffyw4deoUQkNDUatWLRgZGaF3795ISUnRaCsIAr755hvUrl0bNWrUQMeOHfH333+XaV84ODigffv2Bd4DGzduRLNmzUrd70kqlaq/2OQ/pxs3bqBVq1aFfmiW5jVzdnZG9+7dcfDgQXh4eEChUKBx48bYvn27Rrv8fXbs2DGMGjUK1tbWqF27tvrx0hxH80VFRaFNmzYwNDSEi4sLli5dqvF4Tk4OpkyZAk9PT5iZmcHIyAi+vr44cuRIkc/jhx9+gJOTEwwNDdGhQwdcvnxZ4/HS9JV6uY9XUe/5zMxMGBkZYezYsQXWcefOHchkMsycObPYbT1+/BjDhg2DmZkZzM3NMXTo0EL3V2G5Dx06hHbt2sHc3BzGxsZo2LAhJk6cqH4OrVq1AgAEBQUV+CzJ/yyLiopC+/btUaNGDfWyL/fxyqdUKjFx4kTY2trCyMgIPXv21Dg2AkX3eXxxnSVlK+xzPisrC5999pn6uNewYUPMmTMHgiBotJNIJBg9ejR27tyJpk2bQi6Xo0mTJjhw4IBGu4yMDIwbNw7Ozs6Qy+WwtrbGW2+9hXPnzhXIXhStFF5paWlITU3VmB48eADgeYXau3dv7Ny5Ezk5ORrL7dy5E9nZ2RgwYAAA4OnTp/Dz88P69esxaNAgzJ49G2ZmZhg2bBj+97//lUvW9evXQy6Xw9fXF+vXr8f69evx0UcfFdl+2rRpCAkJgb29PebOnYs+ffpg2bJlePvtt5Gbm6vR9tGjR+jcubO6n4Wbmxu++OIL7N+/v8Rcw4cPx5QpU9CyZUv88MMP6NChA2bOnKneN/nZfX19IZfL1dnbt29f6PomTZoEDw8PWFlZqdu+3A/hu+++w44dOzB+/HiEhYXhzz//LFDg/PHHH2jfvj3S09MxdepUzJgxA48fP8abb76JyMjIEp9Xabi4uKBNmzb45ZdfoFQqNR7L/yAqy+WQIUOGACj9afF+/frh2bNnmDlzJrp27YoFCxZg5MiRGm2+/fZbBAYGon79+pg3bx7GjRuH8PBwtG/fXn0wLGqfSyQStG3bFsePH1ev7+LFi0hLSwMAjSL2xIkTaNGihfpy2d9//4033ngDMTExmDBhAubOnQsjIyP06tULO3bsUC/35MkTdOjQARs2bEBgYCAWLFiAtm3bIiwsTKPgy7dp0ybMnj0bH330Eb755hvEx8fj3XffLfB/+mVOTk64ffs2/vjjj1Lt2+L89ttvAIDAwMBCHzczM8M777yDK1euqL8QODk5QalUYv369a+0zbNnz+L06dMYMGAAFixYgI8//hjh4eHw8/PDkydPCrQfM2YMLly4gKlTpyI4OBi//fYbRo8erdFmypQpmDx5Mtzd3dU3vrz99tvIysoqU7aBAwfit99+UxfdeXl5+PXXX8t8KfDGjRsAAEtLSwDP91l4eDju3LlTpvW86Nq1a+jfvz+6dOmCmTNnQk9PD++9916hX25GjRqFf/75B1OmTFFfAivrcbRr167w9PTE999/j9q1ayM4OBirVq1St0lPT8eKFSvg5+eHWbNmYdq0aUhJSUFAQEChfYPWrVuHBQsWICQkBGFhYbh8+TLefPNNJCUlvfI+AYp+zxsbG6N3797YsmVLgWPa5s2bIQhCsV8mBUHAO++8g/Xr12Pw4MH45ptvcOfOHQwdOrTETH///Te6d++O7OxsfPXVV5g7dy569uypPs40atQIX331FQBg5MiRhX6WPHjwAF26dIGHhwfmz5+Pjh07FrvNb7/9Fnv37sUXX3yBTz75BIcOHYK/vz+ePn1aYt4XlSbbiwRBQM+ePfHDDz+gc+fOmDdvHho2bIjPP/+80OPeyZMnMWrUKAwYMADff/89nj17hj59+qjrFQD4+OOPsWTJEvTp0weLFy/G+PHjYWhoiJiYmNI/EaECrV69WgBQ6CSXy9Xtfv/9dwGA8Ntvv2ks37VrV6Fu3brqv+fPny8AEDZs2KCel5OTI/j4+AjGxsZCenq6ej4AYerUqeq/hw4dKjg5ORXIOHXqVOHl3WBkZCQMHTq0yOcTFxcnCIIgJCcnCwYGBsLbb78tKJVKdbuFCxcKAIRVq1ap53Xo0EEAIKxbt049Lzs7W7C1tRX69OlTYFsvio6OFgAIw4cP15g/fvx4AYDwxx9/aDxPIyOjYteXr1u3boXukyNHjggAhEaNGgnZ2dnq+f/73/8EAMKlS5cEQRAElUol1K9fXwgICBBUKpW63ZMnTwQXFxfhrbfeKlUOQRCEX3/9VQAgHDlypNDHFy1aJAAQfv/9d/U8pVIpODg4CD4+Phpt81/TlJSUQtf16NEjAYDQu3fvYjPlr6dnz54a80eNGiUAEC5cuCAIgiDEx8cLMplM+PbbbzXaXbp0SdDT09OYX9Q+nz17tiCTydT/hxcsWCA4OTkJrVu3Fr744gv18zU3Nxc+/fRT9XKdOnUSmjVrJjx79kw9T6VSCW3atBHq16+vnvf1118LRkZGwtWrVzW2O2HCBEEmkwkJCQmCIAhCXFycAECwtLQUHj58qG63a9euQt+jL7t8+bJgaGgoABA8PDyEsWPHCjt37hSysrIKtO3QoYPQpEmTItfl4eEhmJmZFbu9efPmCQCE3bt3C4IgCImJiUKtWrUEAIKbm5vw8ccfC5s2bRIeP35c7HryPXnypMC8iIiIAu/d/GOBv7+/xv/9Tz/9VJDJZOrt5R8junXrptFu4sSJAoBCjzMvAyCEhIQIDx8+FAwMDIT169cLgiAIe/fuFSQSiRAfH1/o//n8Y0FKSoqQkpIiXL9+XZgxY4YgkUiE5s2bq9utXLlSACAYGBgIHTt2FCZPniycOHFC45hWHCcnJwGAsG3bNvW8tLQ0wc7OTmjRokWBfdauXTshLy9PPf9VjqNz585Vz8vOzhY8PDwEa2trIScnRxAEQcjLy9M4dgnC8/e9jY2N8MEHH6jn5f9/NzQ0FO7cuaOef+bMGQGAxnutsM8KJycnjdcw/9j54nGsqPd8/ufe/v37NeY3b95c6NChQ4H2L9q5c6cAQPj+++/V8/Ly8gRfX18BgLB69eoic//www/FHh8FQRDOnj1bYD358l+DpUuXFvrYi9nz94eDg4PG5/Mvv/wiABD+97//qee9vC+LWmdx2V7+nM/fT998841Gu759+woSiUS4fv26el7+e+DFeRcuXBAACD/++KN6npmZmRASElJg22WhlTNeixYtwqFDhzSmF8/yvPnmm7CyssKWLVvU8x49eoRDhw6p+24AwL59+2Bra4v3339fPU9fXx+ffPIJMjMzcezYMW08HbXDhw8jJycH48aNg1T6364cMWIETE1NsXfvXo32xsbGGv0tDAwM0Lp1a9y8ebPY7ezbtw8AClTon332GQAU2E55CQoK0rj84OvrCwDqvNHR0bh27RoGDhyIBw8eqM9mZmVloVOnTjh+/Hi5dXTs378/9PX1NS61HDt2DHfv3i31ZcZ8+WeLMjIyStU+JCRE4+8xY8YA+O912b59O1QqFfr166dxVtfW1hb169cv9vJGPl9fXyiVSpw+fRrA8zNbvr6+8PX1xYkTJwAAly9fxuPHj9Wvw8OHD/HHH3+gX79+yMjI0DibHBAQgGvXruHu3bsAnneg9vX1hYWFhUZGf39/KJVKjbNtwPP9bWFhoZEPQIn/V5s0aYLo6GgMHjwY8fHx+N///odevXrBxsYGy5cvL3E/vCgjIwMmJibFtsl/PD09HQBgY2ODCxcu4OOPP8ajR4+wdOlSDBw4ENbW1vj6668LXF54maGhofrfubm5ePDgAerVqwdzc/NCLyWMHDlS4zJO/ut469YtAP8dI8aMGaPR7lU6WltYWKBz587YvHkzgOdnJdu0aVPoZd18WVlZqFWrFmrVqoV69eph4sSJ8PHx0Tgb+sEHH+DAgQPw8/PDyZMn8fXXX8PX1xf169dX/38sib29PXr37q3+29TUFIGBgTh//jwSExM12o4YMUKjr2xZj6N6enoaVyEMDAzw0UcfITk5GVFRUQCeX4rOP3apVCo8fPgQeXl58PLyKvR17NWrFxwcHNR/t27dGt7e3ur3eEXw9/eHvb09Nm7cqJ53+fJlXLx4scCNQi/bt28f9PT0EBwcrJ4nk8nUx6bimJubAwB27dr1ysdnuVxeaPegogQGBmq8l/v27Qs7O7sK3b/A8/0kk8nwySefaMz/7LPPIAhCgatN/v7+cHV1Vf/dvHlzmJqaahz3zM3NcebMGdy7d++Vc2ml8GrdujX8/f01phdPTerp6aFPnz7YtWuXus/S9u3bkZubq1F43bp1C/Xr19d4cwJQ3wWXf7DTlvztNWzYUGO+gYEB6tatWyBP7dq1C1xrt7CwwKNHj0rcjlQqLXAXqK2tLczNzSvsedepU0fj7/wP4vy8165dAwAMHTpUfXDPn1asWIHs7Gz15bLXZWlpiYCAAOzYsUPdv2rTpk3Q09NDv379yrSu/Es1JX2o56tfv77G366urpBKpeo+MteuXYMgCKhfv36B/RATE1OqDsotW7ZEjRo11EVWfuHVvn17/PXXX3j27Jn6sXbt2gF4fgehIAiYPHlyge3mD92Qv+1r167hwIEDBdr5+/trtMtX0mtfnAYNGmD9+vVITU3FxYsXMWPGDOjp6WHkyJE4fPhwicvnMzExKbE4zn/8xdfSzs4OS5Yswf379xEbG4sFCxagVq1amDJlClauXFns+p4+fYopU6ao+4NYWVmhVq1aePz4caH/l0vaT/nvzZf/D9WqVUujsC2tgQMH4tChQ0hISMDOnTtLvMyoUCjUX3aPHz+O27dv49SpU6hbt65Gu4CAAPz+++94/Pgxjh8/jpCQENy6dQvdu3cv1f/fevXqFTi2NWjQAAAK9I97+Q73sh5H7e3tC3TIL2xba9euRfPmzaFQKGBpaYlatWph7969hb6OL78++eusyDEbpVIpBg0ahJ07d6ovY2/cuBEKhQLvvfdescveunULdnZ2Be7QfXkfFqZ///5o27Ythg8fDhsbGwwYMAC//PJLmYowBweHMnWkf3n/SiQS1KtXr8LHxLx16xbs7e0LHOuLqhlefj8DBT+jv//+e1y+fBmOjo5o3bo1pk2bVuIX0pdp5a7G0hgwYACWLVuG/fv3o1evXvjll1/g5uYGd3f3cll/UZ0iX76+XpGKuiOypG/h+bQ1CF6+kvLmv1Fnz55d5C3H5Xnr/uDBg7Fnzx7s2bMHPXv2xLZt2/D222+jVq1aZVpPfqfZVx3O5OXXQaVSQSKRYP/+/YXus9LsA319fXh7e+P48eO4fv06EhMT4evrCxsbG+Tm5uLMmTM4ceIE3Nzc1M83f/+PHz8eAQEBha43/zmqVCq89dZb+L//+79C2+V/cOV73f+r+eto1qwZmjVrBh8fH3Ts2BEbN25UF3sladSoEaKjo5GQkFDoARF43hcOABo3blzgMYlEggYNGqBBgwbo1q0b6tevj40bNxY5Hhbw/Gzm6tWrMW7cOPj4+MDMzAwSiQQDBgwo9IOpPPZTWfTs2RNyuRxDhw5FdnZ2iV86ZDJZqfc3ANSoUUN9ptXKygrTp0/H/v37S9V3qLRePKtYUTZs2IBhw4ahV69e+Pzzz2Ftba3usJ7fx60yCAwMxOzZs7Fz5068//772LRpE7p37w4zM7MK26ahoSGOHz+OI0eOYO/evThw4AC2bNmCN998EwcPHizVnfsV8RoW9xmtrdEESvN+7tevH3x9fbFjxw4cPHgQs2fPxqxZs7B9+3Z06dKlVNupNIVX+/btYWdnhy1btqBdu3b4448/MGnSJI02Tk5OuHjxIlQqlcZZrytXrqgfL4qFhUWhd3wUdraotAVO/vZiY2M1vkHm5OQgLi6uTAe8krajUqlw7do1jTGukpKS8Pjx42Kfd3Fet5DLPyVrampabs+1OD179oSJiQk2bdoEfX19PHr0qMyXGQGoO14XVay87Nq1axrf0q9fvw6VSqW+e8bV1RWCIMDFxaVAAfOy4va5r68vZs2ahcOHD8PKygpubm6QSCRo0qQJTpw4gRMnTqB79+7q9vn/5/T19Uvc/66ursjMzNTK61SY/OFk7t+/X+plunfvjs2bN2PdunWFDqCbnp6OXbt2wc3NrcQium7durCwsChx+1u3bsXQoUMxd+5c9bxnz54VeXddSfLfm9euXdM4RqSkpJTq7OHLDA0N0atXL2zYsAFdunTRGC6lvJXlNcs/+/ri/++rV68CQIkjiZf1OHrv3r0Cw1C8vK2tW7eibt262L59u0amlwfxzZd/9v5FV69eLZdR0It7zzdt2hQtWrTAxo0bUbt2bSQkJODHH38scZ35N0RkZmZqfLGLjY0tVSapVIpOnTqhU6dOmDdvHmbMmIFJkybhyJEj8Pf3L/cv+S/vX0EQcP36dTRv3lw9r7jP6Bf/X5Qlm5OTEw4fPlyg20Jpaobi2NnZYdSoURg1ahSSk5PRsmVLfPvtt6UuvEQfTiKfVCpF37598dtvv2H9+vXIy8vTuMwIAF27dkViYqJGX7C8vDz8+OOPMDY2RocOHYpcv6urK9LS0tTfkIHnB5QX+zrkMzIyKtWB1t/fHwYGBliwYIFGRbxy5UqkpaWhW7duJa6jNLp27QoABe46nDdvHgC88naMjIxe61Kgp6cnXF1dMWfOnEJHBH/5tvrXZWhoiN69e2Pfvn1YsmQJjIyMNMY0Ko1NmzZhxYoV8PHxQadOnUq1zKJFizT+zj8w5r/J3n33XchkMkyfPr3AmQ5BEDTuiClun/v6+iI7Oxvz589Hu3bt1AeY/Dts7927p+5rBTy/1d/Pzw/Lli0r9MPxxf3fr18/RERE4Pfffy/Q7vHjx8jLyyt2H5TWiRMnCr3zMb8vR2kuheTr27cvGjdujO+++w5//fWXxmMqlQrBwcF49OiRxofpmTNnCr1bMDIyEg8ePChx+zKZrMBr+OOPP77ymXF/f3/o6+vjxx9/1Fjv64xkPn78eEydOhWTJ09+5XW8KDw8vND5ZXnN7t27p3EsTU9Px7p16+Dh4QFbW9tily3rcTQvLw/Lli1T/52Tk4Nly5ahVq1a8PT0BPDfmYsX13fmzBlEREQUmmHnzp3q/pDA8/8vZ86cKfUHaXFKOs4OGTIEBw8exPz582FpaVmqbXbt2hV5eXlYsmSJep5SqSxV0fbw4cMC8/KvWOR39ckval/1C8fL1q1bp9FtYOvWrbh//77Gc3V1dcWff/6pMbrBnj17Cgw7UZZsXbt2hVKpxMKFCzXm//DDD5BIJGV+fZVKZYHX0traGvb29mX6RRStnPHav3+/usJ8UZs2bTQq2f79++PHH3/E1KlT0axZM42zO8DzjqzLli3DsGHDEBUVBWdnZ2zduhWnTp3C/Pnzi+2zM2DAAHzxxRfo3bs3PvnkEzx58gRLlixBgwYNCnS29PT0xOHDh9UDDLq4uBQYRwh43k8jLCwM06dPR+fOndGzZ0/ExsZi8eLFaNWqVYkdJEvL3d0dQ4cOxU8//YTHjx+jQ4cOiIyMxNq1a9GrV68Sb+UtiqenJ7Zs2YLQ0FC0atUKxsbG6NGjR6mXl0qlWLFiBbp06YImTZogKCgIDg4OuHv3Lo4cOQJTU1P1kADlZfDgwVi3bh1+//13DBo0qNjBF7du3QpjY2Pk5OSoR64/deoU3N3d8euvv5Z6m3FxcejZsyc6d+6MiIgIbNiwAQMHDlRfBnd1dcU333yDsLAwxMfHo1evXjAxMUFcXBx27NiBkSNHqseiK26f+/j4QE9PD7GxsRrDVbRv3159gH2x8AKeF4Xt2rVDs2bNMGLECNStWxdJSUmIiIjAnTt3cOHCBQDA559/jt27d6N79+4YNmwYPD09kZWVhUuXLmHr1q2Ij48vl7Mns2bNQlRUFN599131t9lz585h3bp1qFmzZoFO5SkpKYX+IoGLiwsGDRqErVu3olOnTmjXrp3GyPWbNm3CuXPn8NlnnxUYUmXjxo3o3bs3PD09YWBggJiYGKxatQoKhUI93lBRunfvjvXr18PMzAyNGzdGREQEDh8+rB56oaxq1aqF8ePHY+bMmejevTu6du2K8+fPY//+/a+8v93d3cutCwYAvPPOO3BxcUGPHj3g6uqKrKwsHD58GL/99htatWpVqmNCgwYN8OGHH+Ls2bOwsbHBqlWrkJSUhNWrV5e4bFmPo/b29pg1axbi4+PRoEEDbNmyBdHR0fjpp5/UA2h2794d27dvR+/evdGtWzfExcVh6dKlaNy4caFfEuvVq4d27dohODhY/eXH0tKyyEvzZVHScXbgwIH4v//7P+zYsQPBwcGlGgS0R48eaNu2LSZMmID4+Hj1uGml+SL91Vdf4fjx4+jWrRucnJyQnJyMxYsXo3bt2ur+o66urjA3N8fSpUthYmICIyMjeHt7F/oLNKVRs2ZN9Xs4KSkJ8+fPR7169TBixAh1m+HDh2Pr1q3o3Lkz+vXrhxs3bmDDhg0and3Lmq1Hjx7o2LEjJk2ahPj4eLi7u+PgwYPYtWsXxo0bV2DdJcnIyEDt2rXRt29fuLu7w9jYGIcPH8bZs2c1zpKX6LXuiSxBccNJoJDbQVUqleDo6Fjo7Z/5kpKShKCgIMHKykowMDAQmjVrVuhtpXhpOAlBEISDBw8KTZs2FQwMDISGDRsKGzZsKPQW4StXrgjt27dX3xaff4vry8NJ5Fu4cKHg5uYm6OvrCzY2NkJwcLDw6NEjjTZF3Tpf1DAXL8vNzRWmT58uuLi4CPr6+oKjo6MQFhamMYxA/vpKO5xEZmamMHDgQMHc3FwAoM6Rfwvwr7/+qtE+/9brl/f3+fPnhXfffVewtLQU5HK54OTkJPTr108IDw8vVQ5BKHk4iXx5eXmCnZ2dAEDYt29foW3yX9P8SaFQCLVr1xa6d+8urFq1qsA+K0r+ev755x+hb9++gomJiWBhYSGMHj1aePr0aYH227ZtE9q1aycYGRkJRkZGgpubmxASEiLExsaq2xS1z/O1atVKACCcOXNGPe/OnTsCAMHR0bHQnDdu3BACAwMFW1tbQV9fX3BwcBC6d+8ubN26VaNdRkaGEBYWJtSrV08wMDAQrKyshDZt2ghz5sxR34af/xrPnj27wHYKe0+97NSpU0JISIjQtGlTwczMTNDX1xfq1KkjDBs2TLhx44ZG2/zb0gubOnXqpG6XnJwshIaGCvXq1RPkcrlgbm4u+Pv7q4eQeNHFixeFzz//XGjZsqVQs2ZNQU9PT7CzsxPee+894dy5c8VmF4TnQw7kH1+MjY2FgIAA4cqVKwVudc8/Fpw9e1Zj+cKGE1AqlcL06dMFOzs7wdDQUPDz8xMuX75c5O3zL8O/w0kUp7jhJEqyefNmYcCAAYKrq6tgaGgoKBQKoXHjxsKkSZM0hgAoipOTk9CtWzfh999/F5o3by7I5XLBzc2twPGjqH2WryzH0b/++kvw8fERFAqF4OTkJCxcuFCjnUqlEmbMmCE4OTkJcrlcaNGihbBnz54Cx9sX/7/PnTtXcHR0FORyueDr66seLibfqw4nUdJ7XhCeD50EQDh9+nSh+6YwDx48EIYMGSKYmpoKZmZmwpAhQ4Tz58+XOJxEeHi48M477wj29vaCgYGBYG9vL7z//vsFhprZtWuX0LhxY0FPT09jncUNA1PUcBKbN28WwsLCBGtra8HQ0FDo1q2bcOvWrQLLz507V3BwcBDkcrnQtm1b4a+//iqwzuKyFfZ5mpGRIXz66aeCvb29oK+vL9SvX1+YPXu2xvAuglD0++zF1zg7O1v4/PPPBXd3d8HExEQwMjIS3N3dhcWLFxe6P4oi+XeDREREZebs7IymTZsW+aPlVLLevXvj0qVLhf4yCFU/laaPFxERka65f/8+9u7dq/5FDar+Ks1djURERLoiLi4Op06dwooVK6Cvr1/sT9NR9cIzXkRERFp27NgxDBkyBHFxcVi7dm2Jd39S9cE+XkRERERawjNeRERERFrCwouIiIhIS3Suc71KpcK9e/dgYmKi9d8+JCIiolcjCAIyMjJgb2+v8bOBVU6ZRv2qAAsXLlQPcte6dWuNgSML88MPPwgNGjRQD4o5bty4QgezLMrt27eLHdSVEydOnDhx4lR5p9u3b79u6SEqUc945f+MwtKlS+Ht7Y358+cjICAAsbGxsLa2LtB+06ZNmDBhAlatWoU2bdrg6tWrGDZsGCQSifp3C0uS/7NCt2/fhqmpabk+HyIiIqoY6enpcHR0LPbnAasCUe9q9Pb2RqtWrdQ/YKlSqeDo6IgxY8ZgwoQJBdqPHj0aMTExGj/q+tlnn+HMmTM4efJkqbaZnp4OMzMzpKWlsfAiIiKqIqrL57doF0lzcnIQFRUFf3///8JIpfD39y/yF+TbtGmDqKgoREZGAgBu3ryJffv2oWvXrkVuJzs7G+np6RoTERERkRhEu9SYmpoKpVIJGxsbjfk2Nja4cuVKocsMHDgQqampaNeuHQRBQF5eHj7++GNMnDixyO3MnDkT06dPL9fsRERERK+iSt0WcPToUcyYMQOLFy/GuXPnsH37duzduxdff/11kcuEhYUhLS1NPd2+fVuLiYmIiIj+I9oZLysrK8hkMiQlJWnMT0pKKvKnEyZPnowhQ4Zg+PDhAIBmzZohKysLI0eOxKRJkwq9vVQul0Mul5f/EyAiIiIqI9HOeBkYGMDT01Ojo7xKpUJ4eDh8fHwKXebJkycFiiuZTAYAEPEeASIiIqJSEXU4idDQUAwdOhReXl5o3bo15s+fj6ysLAQFBQEAAgMD4eDggJkzZwIAevTogXnz5qFFixbw9vbG9evXMXnyZPTo0UNdgBERERFVVqIWXv3790dKSgqmTJmCxMREeHh44MCBA+oO9wkJCRpnuL788ktIJBJ8+eWXuHv3LmrVqoUePXrg22+/FespEBEREZWaqON4iaG6jANCRESkS6rL53eVuquRiIiIqCpj4UVERESkJSy8iIiIiLSEhRcRERGRloh6V2N1kp2nREpGttgxqJqyNJLD0IBDphARVXUsvMrJ3/fS8e7i02LHoGrKRK6H7/s2R5dmdmJHISKi18DCq5xIAMj1eOWWyp8gABnZeQjeeA6j/Fzx2dsNIZNKxI5FRESvgON4EVVyeUoVvtt/BStOxgEA2jeohQUDPGBew0DkZERE2lNdPr95ioaoktOTSfFl98b43wAPKPSlOH41BT0WnsQ/99LFjkZERGXEwouoinjHwwHbg9vCsaYhbj98ineXnMKu6LtixyIiojJg4UVUhTS2N8Vvo9vBt74VnuWqMPbnaHyz5x/kKVViRyMiolJg4UVUxZjXMMCaoNYI9nMFAKw4GYfAVZF4kMnhTIiIKjsWXkRVkEwqwRed3bB4UEvUMJDh9I0H6LnwFC7dSRM7GhERFYOFF1EV1rWZHXaGtIWLlRHuPn6KPktPY2vUHbFjERFREVh4EVVxDWxMsDOkLd50s0ZOngrjf72AqbsuI5f9voiIKh0WXkTVgJmhPlYEeuGTTvUBAGsjbmHQ8jNIzngmcjIiInoRCy+iakIqlSD0rQZYHugFY7keIuMfosePJ3Eu4ZHY0YiI6F8svIiqmbca22DX6LZwrWWEpPRsDFj2JzZHJogdi4iIwMKLqFpyrWWMXaPboXMTW+QoVQjbfglh2y8hO08pdjQiIp3GwouomjKW62HJ4Jb4PKAhJBJgc2QCBvz0JxLT2O+LiEgsLLyIqjGJRIKQjvWwelgrmCr0cD7hMbr/eBJn4x+KHY2ISCex8CLSAX4NrfHbmHZwszVBamY23v/pT6yLiIcgCGJHIyLSKSy8iHSEk6URto9qgx7u9shTCZiy62+M//UinuWy3xcRkbaw8CLSITUM9LBggAcmdW0EqQTYdu4O3lsagbuPn4odjYhIJ7DwItIxEokEI9rXxfoPvWFRQx+X7qahx48ncfpGqtjRiIiqPRZeRDqqbT0r/DamHZrYm+JhVg6GrIzEihM32e+LiKgCsfAi0mG1LWpgW3AbvNvSAUqVgG/2xmDsz9F4msN+X0REFYGFF5GOU+jLMPc9d0zr0Rh6Ugl2X7iH3otPIeHBE7GjERFVOyy8iAgSiQTD2rpg43BvWBkb4EpiBnosPInjV1PEjkZEVK2w8CIiNe+6lvhtTDu4O5oj7Wkuhq6OxOKj19nvi4ionLDwIiINdmaG+OWjNzCglSMEAfj+QCxCNp1DZnae2NGIiKo8Fl5EVIBcT4bv+jTHjN7NoC+TYN+lRPRedApxqVliRyMiqtIqReG1aNEiODs7Q6FQwNvbG5GRkUW29fPzg0QiKTB169ZNi4mJdMNA7zr4eaQPrE3kuJaciZ4LTyI8JknsWEREVZbohdeWLVsQGhqKqVOn4ty5c3B3d0dAQACSk5MLbb99+3bcv39fPV2+fBkymQzvvfeelpMT6QZPJwvsGdMOXk4WyHiWhw/X/oX5h69CpWK/LyKishK98Jo3bx5GjBiBoKAgNG7cGEuXLkWNGjWwatWqQtvXrFkTtra26unQoUOoUaMGCy+iCmRtqsCmEW8g0McJADD/8DWMXP8X0p/lipyMiKhqEbXwysnJQVRUFPz9/dXzpFIp/P39ERERUap1rFy5EgMGDICRkVFFxSQiAAZ6Unz1TlPM7tscBnpSHI5JRq+Fp3AtKUPsaEREVYaohVdqaiqUSiVsbGw05tvY2CAxMbHE5SMjI3H58mUMHz68yDbZ2dlIT0/XmIjo1b3n5YitH/vA3kyBm6lZ6LXoFA5cvi92LCKiKkH0S42vY+XKlWjWrBlat25dZJuZM2fCzMxMPTk6OmoxIVH11Ly2OX4b0w5v1K2JrBwlPt5wDt8fuAIl+30RERVL1MLLysoKMpkMSUmad0klJSXB1ta22GWzsrLw888/48MPPyy2XVhYGNLS0tTT7du3Xzs3EQGWxnJs+NAbw9u5AAAWH72BoDVn8fhJjsjJiIgqL1ELLwMDA3h6eiI8PFw9T6VSITw8HD4+PsUu++uvvyI7OxuDBw8utp1cLoepqanGRETlQ08mxZfdG+N/Azyg0Jfi+NUU9Fx4CjH3eUmfiKgwol9qDA0NxfLly7F27VrExMQgODgYWVlZCAoKAgAEBgYiLCyswHIrV65Er169YGlpqe3IRPSSdzwcsD24LRxrGiLh4RO8u/g0dl+4J3YsIqJKR0/sAP3790dKSgqmTJmCxMREeHh44MCBA+oO9wkJCZBKNevD2NhYnDx5EgcPHhQjMhEVorG9KX4b3Q5jNp/HiWup+GTzeVy68xhfdHaDnkz073hERJWCRNCxX79NT0+HmZkZ0tLSeNmRqAIoVQLmHIzFkqM3AABtXC2xcGBL1DQyEDkZEVVl1eXzm19DiahcyaQSfNHZDYsHtUQNAxlO33iAHj+exOW7aWJHIyISHQsvIqoQXZvZYWdIW7hYGeHu46fos+Q0tkXdETsWEZGoWHgRUYVpYGOCnSFt8aabNbLzVPjs1ws4cqXw32ElItIFLLyIqEKZGepjRaAX3m9dBwAw7be/8SxXKXIqIiJxsPAiogonlUowqVsjWJvIcevBE6w8GSd2JCIiUbDwIiKtMJbrYVK3RgCAH/+4hruPn4qciIhI+1h4EZHW9HS3R2vnmniWq8KMvTFixyEi0joWXkSkNRKJBNN6NoFUAuy9dB+nrqeKHYmISKtYeBGRVjW2N8WQN5wAAFN3/41cpUrkRERE2sPCi4i0LvSthqhpZIDryZlYezpe7DhERFrDwouItM6shj6+6NwQADD/8DUkpz8TORERkXaw8CIiUbzn6Qj32mbIzM7Dd/uviB2HiEgrWHgRkSikUgm+eqcpJBJg+/m7+Cv+odiRiIgqHAsvIhKNu6M5+ns5AgCm7PobSpUgciIioorFwouIRPV5QEOYKvTwz/10bIpMEDsOEVGFYuFFRKKyNJZjfMDzjvZzfo/Fw6wckRMREVUcFl5EJLqBreugkZ0p0p7mYvbvsWLHISKqMCy8iEh0ejIpvnqnCQDg57MJuHjnsbiBiIgqCAsvIqoUWjnXRO8WDhCE5x3tVexoT0TVEAsvIqo0wrq4wchAhujbj7H13B2x4xARlTsWXkRUaVibKjDOvwEAYNb+K0h7mityIiKi8sXCi4gqlaFtnOFaywgPsnLww6GrYschIipXLLyIqFIx0JNiWs/nHe3XRcQj5n66yImIiMoPCy8iqnR869dCl6a2UAnA1N1/QxDY0Z6IqgcWXkRUKU3q1ggKfSki4x5i94V7YschIioXLLyIqFKqbVEDIX71AAAz9sUgKztP5ERERK+PhRcRVVoj2tdFnZo1kJSejR//uC52HCKi18bCi4gqLYW+DFN7NAYArDx5EzdSMkVORET0elh4EVGl1qmRDd50s0auUsA0drQnoiqOhRcRVXpTujeGgUyKE9dScfCfJLHjEBG9MhZeRFTpOVsZYWT7ugCAr377B89ylSInIiJ6NSy8iKhKGNXRFfZmCtx9/BRLjt4QOw4R0Sth4UVEVUINAz182f15R/slx24g4cETkRMREZWd6IXXokWL4OzsDIVCAW9vb0RGRhbb/vHjxwgJCYGdnR3kcjkaNGiAffv2aSktEYmpS1NbtK1niZw8Fb7e+4/YcYiIykzUwmvLli0IDQ3F1KlTce7cObi7uyMgIADJycmFts/JycFbb72F+Ph4bN26FbGxsVi+fDkcHBy0nJyIxCCRSDCtRxPoSSU49E8SjsQWfqwgIqqsJIKI92Z7e3ujVatWWLhwIQBApVLB0dERY8aMwYQJEwq0X7p0KWbPno0rV65AX1//lbaZnp4OMzMzpKWlwdTU9LXyE5E4vtnzD1acjIOLlREOjPOFXE8mdiQiqmDV5fNbtDNeOTk5iIqKgr+//39hpFL4+/sjIiKi0GV2794NHx8fhISEwMbGBk2bNsWMGTOgVBZ9h1N2djbS09M1JiKq2sb614eVsRxxqVlYeTJO7DhERKUmWuGVmpoKpVIJGxsbjfk2NjZITEwsdJmbN29i69atUCqV2LdvHyZPnoy5c+fim2++KXI7M2fOhJmZmXpydHQs1+dBRNpnotDHxK5uAIAfw6/jftpTkRMREZWO6J3ry0KlUsHa2ho//fQTPD090b9/f0yaNAlLly4tcpmwsDCkpaWpp9u3b2sxMRFVlN4tHODlZIGnuUrM2HdF7DhERKUiWuFlZWUFmUyGpCTNUaiTkpJga2tb6DJ2dnZo0KABZLL/+nM0atQIiYmJyMnJKXQZuVwOU1NTjYmIqj6JRILp7zSBVAL8duEeIm48EDsSEVGJRCu8DAwM4OnpifDwcPU8lUqF8PBw+Pj4FLpM27Ztcf36dahUKvW8q1evws7ODgYGBhWemYgqlyb2Zhjk7QQAmLb7b+QqVSUsQUQkLlEvNYaGhmL58uVYu3YtYmJiEBwcjKysLAQFBQEAAgMDERYWpm4fHByMhw8fYuzYsbh69Sr27t2LGTNmICQkRKynQEQi++ztBrCooY/YpAysj7gldhwiomLpibnx/v37IyUlBVOmTEFiYiI8PDxw4MABdYf7hIQESKX/1YaOjo74/fff8emnn6J58+ZwcHDA2LFj8cUXX4j1FIhIZOY1DPB/nd0Qtv0Sfjh0FT3c7VHLRC52LCKiQok6jpcYqss4IET0H6VKQO/Fp3DxThr6etbGnPfcxY5EROWsunx+V6m7GomICiOTSjC9ZxMAwNaoO4i69UjkREREhWPhRUTVQos6FujnVRsAMHX3ZShVOnUyn4iqCBZeRFRt/F9nN5go9HD5bjp+PpsgdhwiogJYeBFRtWFlLMdnbzUAAMz+PRaPsgof34+ISCwsvIioWhn8hhPcbE3w+Eku5hyMFTsOEZEGFl5EVK3oyaSY9m9H+02RCbh8N03kRERE/2HhRUTVzht1LdHT3R6CAEzZdRkqdrQnokqChRcRVUsTuzZCDQMZziU8xo7zd8WOQ0QEgIUXEVVTtmYKfNKpPgBg5v4rSH+WK3IiIiIWXkRUjX3Q1gV1rYyQmpmN/x2+JnYcIiIWXkRUfRno/dfRfs3peFxNyhA5ERHpOhZeRFSttW9QCwFNbKBUCZi662/o2M/TElElw8KLiKq9L7s1hlxPioibD7D30n2x4xCRDmPhRUTVnmPNGhjlVw8A8O3eGGRl54mciIh0FQsvItIJH3WoC8eahrif9gyLjlwXOw4R6SgWXkSkExT6Mkzp/ryj/fITN3EzJVPkRESki1h4EZHO8G9kDb+GtZCrFDD9t3/Y0Z6ItI6FFxHpDIlEgqk9msBAJsWxqyk4HJMsdiQi0jEsvIhIp7hYGWG4rwsA4Ks9f+NZrlLkRESkS1h4EZHOCelYD7amCtx++BTLjt0UOw4R6RAWXkSkc4zkepjUrREAYPHR67j98InIiYhIV7DwIiKd1L25Hd6oWxPZeSp8uzdG7DhEpCNYeBGRTpJIJJjesylkUgkO/J2I41dTxI5ERDqAhRcR6ayGtiYY6uMMAJj229/IyVOJG4iIqj0WXkSk08a9VR9Wxga4mZKF1afixI5DRNUcCy8i0mmmCn1M6PK8o/2C8GtITHsmciIiqs5YeBGRznu3hQNa1jFHVo4SM/ezoz0RVRwWXkSk86RSCb56pykkEmBX9D2cuflA7EhEVE2x8CIiAtDUwQwDW9cBAEzd/TfylOxoT0Tlj4UXEdG/xr/dEOY19HElMQMb/rwldhwiqoZYeBER/cvCyACfBzQEAMw9dBWpmdkiJyKi6oaFFxHRCwa0qoOmDqbIeJaH7w9cETsOEVUzlaLwWrRoEZydnaFQKODt7Y3IyMgi265ZswYSiURjUigUWkxLRNWZTPp8RHsA+OWvOzif8EjkRERUnYheeG3ZsgWhoaGYOnUqzp07B3d3dwQEBCA5ObnIZUxNTXH//n31dOsW+2IQUfnxdLJAn5a1ATzvaK9SCSInIqLqQvTCa968eRgxYgSCgoLQuHFjLF26FDVq1MCqVauKXEYikcDW1lY92djYaDExEemCL7o0hIlcDxfvpOGXv26LHYeIqglRC6+cnBxERUXB399fPU8qlcLf3x8RERFFLpeZmQknJyc4OjrinXfewd9//62NuESkQ6xNFBj3VgMAwKwDV/D4SY7IiYioOhC18EpNTYVSqSxwxsrGxgaJiYmFLtOwYUOsWrUKu3btwoYNG6BSqdCmTRvcuXOn0PbZ2dlIT0/XmIiISiPQxwkNbIzx6Eku5h26KnYcIqoGRL/UWFY+Pj4IDAyEh4cHOnTogO3bt6NWrVpYtmxZoe1nzpwJMzMz9eTo6KjlxERUVenLpJjWswkAYMOft/D3vTSRExFRVSdq4WVlZQWZTIakpCSN+UlJSbC1tS3VOvT19dGiRQtcv3690MfDwsKQlpamnm7fZl8NIiq9Nq5W6N7cDioBmLrrbwgCO9oT0asTtfAyMDCAp6cnwsPD1fNUKhXCw8Ph4+NTqnUolUpcunQJdnZ2hT4ul8thamqqMRERlcWkbo1gqC/DX7ceYWf0XbHjEFEVJvqlxtDQUCxfvhxr165FTEwMgoODkZWVhaCgIABAYGAgwsLC1O2/+uorHDx4EDdv3sS5c+cwePBg3Lp1C8OHDxfrKRBRNWdnZogxneoBAGbsu4KMZ7kiJyKiqkpP7AD9+/dHSkoKpkyZgsTERHh4eODAgQPqDvcJCQmQSv+rDx89eoQRI0YgMTERFhYW8PT0xOnTp9G4cWOxngIR6YAP27ng17/uIC41CwvCr2FSNx5ziKjsJIKOdVhIT0+HmZkZ0tLSeNmRiMrkaGwyhq0+Cz2pBAfGtUc9a2OxIxHpjOry+S36pUYioqrCr6E1/BtZI08lYPGRwm/oISIqDgsvIqIy+KRTfQDArgv3cOfRE5HTEFFVw8KLiKgMmtc2R7t6VlCqBCw/flPsOERUxbDwIiIqo2A/VwDAz2dvIzUzW+Q0RFSVsPAiIiqjNq6WcK9thuw8Fdaejhc7DhFVISy8iIjKSCKRqM96rT0dz3G9iKjUWHgREb2Ctxvbom4tI6Q/y8PmyASx4xBRFcHCi4joFUilEnzc4flZrxUn4pCdpxQ5ERFVBSy8iIheUS8PB9iaKpCckY3t5/gbjkRUMhZeRESvyEBPiuG+LgCAZcduQKnSqR8CIaJXwMKLiOg1vN+6Dsxr6CP+wRMcuJwodhwiquRYeBERvQYjuR6G+jgDAJYcuw4d+/lbIiojFl5ERK9paBtnGOrLcPluOk5cSxU7DhFVYiy8iIheU00jAwxo7QgAWHL0hshpiKgyY+FFRFQORvjWhZ5UgoibD3A+4ZHYcYiokmLhRURUDuzNDdGrhQMAnvUioqKx8CIiKicfd6gLiQQ4+E8SridniB2HiCohFl5EROWknrUJ3m5sAwBYeuymyGmIqDJi4UVEVI7yf0Zo5/m7uPf4qchpiKiyYeFFRFSOWtSxgE9dS+SpBCw/wbNeRKSJhRcRUTkL9nt+1uvnyNt4mJUjchoiqkxYeBERlTPf+lZo6mCKp7lKrD0dL3YcIqpEWHgREZUziUSC4A71AABrTscjKztP5EREVFmw8CIiqgCdm9rCxcoIaU9zsTkyQew4RFRJsPAiIqoAMqkEH7WvCwBYcSIOOXkqkRMRUWXAwouIqIL0bukAaxM5EtOfYef5u2LHIaJKgIUXEVEFkevJMNzXBQCw9PgNKFWCyImISGwsvIiIKtBAbyeYKvRwMyULh/5JFDsOEYmMhRcRUQUyluthaBtnAM9/PFsQeNaLSJex8CIiqmDD2jhDoS/FhTtpOH3jgdhxiEhELLyIiCqYpbEc/b0cATw/60VEuouFFxGRFoxoXxcyqQQnr6fi4p3HYschIpGUqfBSqVRYtWoVunfvjqZNm6JZs2bo2bMn1q1bx34LRETFqG1RA++42wPgWS8iXVbqwksQBPTs2RPDhw/H3bt30axZMzRp0gS3bt3CsGHD0Lt371cOsWjRIjg7O0OhUMDb2xuRkZGlWu7nn3+GRCJBr169XnnbRETa8vG/P5594O9E3EjJFDkNEYmh1IXXmjVrcPz4cYSHh+P8+fPYvHkzfv75Z1y4cAGHDx/GH3/8gXXr1pU5wJYtWxAaGoqpU6fi3LlzcHd3R0BAAJKTk4tdLj4+HuPHj4evr2+Zt0lEJIYGNibwb2QDQQB+OnZT7DhEJIJSF16bN2/GxIkT0bFjxwKPvfnmm5gwYQI2btxY5gDz5s3DiBEjEBQUhMaNG2Pp0qWoUaMGVq1aVeQySqUSgwYNwvTp01G3bt0yb5OISCzB/5712n7+DhLTnomchoi0rdSF18WLF9G5c+ciH+/SpQsuXLhQpo3n5OQgKioK/v7+/wWSSuHv74+IiIgil/vqq69gbW2NDz/8sMRtZGdnIz09XWMiIhKLp5MFWrvURK5SwIoTPOtFpGtKXXg9fPgQNjY2RT5uY2ODR48elWnjqampUCqVBdZrY2ODxMTCR3g+efIkVq5cieXLl5dqGzNnzoSZmZl6cnR0LFNGIqLyNurfs16bIhPw+EmOyGmISJtKXXgplUro6ekV+bhMJkNeXl65hCpKRkYGhgwZguXLl8PKyqpUy4SFhSEtLU093b59u0IzEhGVpEODWmhkZ4onOUqsi7gldhwi0qKiK6mXCIKAYcOGQS6XF/p4dnZ2mTduZWUFmUyGpKQkjflJSUmwtbUt0P7GjRuIj49Hjx491PNUKhUAQE9PD7GxsXB1ddVYRi6XF5mZiEgMEokEwX6u+GTzeaw+FYfhvi6oYVDqwzERVWGlPuMVGBgIa2trjct2L07W1tYIDAws08YNDAzg6emJ8PBw9TyVSoXw8HD4+PgUaO/m5oZLly4hOjpaPfXs2RMdO3ZEdHQ0LyMSUZXRtakt6tSsgUdPcrHlLM/EE+mKUn/FWrNmTYUECA0NxdChQ+Hl5YXWrVtj/vz5yMrKQlBQEIDnBZ+DgwNmzpwJhUKBpk2baixvbm4OAAXmExFVZnoyKT7qUBeTdlzG8uM3MfgNJ+jL+GMiRNVdqd/lffv2xYEDB8p9hPr+/ftjzpw5mDJlCjw8PBAdHY0DBw6oO9wnJCTg/v375bpNIqLKoE/L2rAyluNe2jPsir4ndhwi0gKJUMpKqlOnTjh69Cjs7e0RFBSEYcOGVckxtNLT02FmZoa0tDSYmpqKHYeIdNySozcw68AV1LM2xsFx7SGVSsSORFQpVZfP71Kf8QoPD8fNmzfx4YcfYsOGDahfvz7efPNNbNq06ZU61hMRETD4jTowUejhenImDscklbwAEVVpZepQ4OTkhGnTpuHmzZs4dOgQ7O3tMWLECNjZ2SEkJARRUVEVlZOIqFoyUehjyBtOAIDFR2+Ue3cOIqpcSn2psSgZGRnYtGkTJk6ciLS0tAofy+t1VZdTlURUfaRkZKPdrD+QnafC5hFvwMfVUuxIRJVOdfn8fq1baOLi4jBnzhzMmDEDaWlpGj/9Q0REpVPLRI5+Xs+Hw1ly7IbIaYioIpW58Hr27Bk2bNiAN998E/Xr18e6devw4YcfIi4uDgcOHKiIjERE1d7I9nUhk0pw/GoKLt9NEzsOEVWQUhdekZGR+Pjjj2FnZ4cRI0bA1tYWBw4cwM2bNzFlyhQOXkpE9Boca9ZA9+Z2AHjWi6g6K/UAqm+88Qbc3d3x9ddfY9CgQbCwsKjIXEREOufjDq7YFX0P+y/dR1xqFlysjMSORETlrNRnvLp3745Tp05h9OjRLLqIiCpAIztTvOlmDZUA/HT8pthxiKgClLrw2rt3LzIzMysyCxGRzgv2cwUAbIu6g6T0ZyKnIaLyVurCi2PLEBFVvFbONeHlZIEcpQqrTsaJHYeIylmZ7mqUSPhTFkREFW1Ux+dnvTb8eQtpT3JFTkNE5anUnesBoEGDBiUWXw8fPnytQEREuq5jQ2s0tDFBbFIGNpy5hZCO9cSORETlpEyF1/Tp02FmZlZRWYiICM+vLgT7uWLclmisOhmHD9q6wNBAJnYsIioHZSq8BgwYAGtr64rKQkRE/+re3A5zDsbizqOn+DXqNgJ9nMWORETloNR9vNi/i4hIe/RkUnzUvi4AYNmxm8hVqkRORETlgXc1EhFVUu95OcLSyAB3Hz/Fnov3xI5DROWg1IWXSqXiZUYiIi1S6MvwQTsXAMCSozegUvELMFFVV+YfySYiIu0Z/IYTjOV6uJqUiSOxyWLHIaLXxMKLiKgSMzPUx6A36gB4ftaLiKo2Fl5ERJXch21dYKAnxV+3HiEyjmMlElVlLLyIiCo5a1MF+nrWBgAsOXpd5DRE9DpYeBERVQEjfetCKgGOxKYg5n662HGI6BWx8CIiqgKcrYzQtZkdAPb1IqrKWHgREVURH3d4/uPZey7eQ8KDJyKnIaJXwcKLiKiKaOpghg4NakElAD+d4FkvoqqIhRcRURUS7Pf8rNcvf91BcsYzkdMQUVmx8CIiqkK8XWqiRR1z5OSpsPpUvNhxiKiMWHgREVUhEokEo/zqAQA2RNxC+rNckRMRUVmw8CIiqmI6uVmjvrUxMrLzsPHPBLHjEFEZsPAiIqpipFKJ+g7HlSfj8CxXKXIiIiotFl5ERFVQTw97OJgbIjUzG1uj7ogdh4hKiYUXEVEVpC+TYoSvCwDgp+M3kadUiZyIiEqDhRcRURXVv1Ud1DQyQMLDJ9h76b7YcYioFCpF4bVo0SI4OztDoVDA29sbkZGRRbbdvn07vLy8YG5uDiMjI3h4eGD9+vVaTEtEVDkYGsgwrI0zgOc/IyQIgriBiKhEohdeW7ZsQWhoKKZOnYpz587B3d0dAQEBSE5OLrR9zZo1MWnSJERERODixYsICgpCUFAQfv/9dy0nJyIS31AfZxgZyHAlMQNHr6aIHYeISiARRP6K5O3tjVatWmHhwoUAAJVKBUdHR4wZMwYTJkwo1TpatmyJbt264euvvy6xbXp6OszMzJCWlgZTU9PXyk5EVBl8u/cfLD8Rh9YuNfHLRz5ixyGqENXl81vUM145OTmIioqCv7+/ep5UKoW/vz8iIiJKXF4QBISHhyM2Nhbt27cvtE12djbS09M1JiKi6uTDdnWhL5MgMu4hom49FDsOERVD1MIrNTUVSqUSNjY2GvNtbGyQmJhY5HJpaWkwNjaGgYEBunXrhh9//BFvvfVWoW1nzpwJMzMz9eTo6Fiuz4GISGy2Zgr0aVkbwPO+XkRUeYnex+tVmJiYIDo6GmfPnsW3336L0NBQHD16tNC2YWFhSEtLU0+3b9/WblgiIi0Y2b4uJBLgcEwyYhMzxI5DREXQE3PjVlZWkMlkSEpK0piflJQEW1vbIpeTSqWoV+/5b5V5eHggJiYGM2fOhJ+fX4G2crkccrm8XHMTEVU2dWsZo0tTW+y7lIilx27gh/4eYkciokKIesbLwMAAnp6eCA8PV89TqVQIDw+Hj0/pO4iqVCpkZ2dXREQioiojuMPzL6S7L9zD7YdPRE5DRIUR/VJjaGgoli9fjrVr1yImJgbBwcHIyspCUFAQACAwMBBhYWHq9jNnzsShQ4dw8+ZNxMTEYO7cuVi/fj0GDx4s1lMgIqoUmtU2g299KyhVAlacuCl2HCIqhKiXGgGgf//+SElJwZQpU5CYmAgPDw8cOHBA3eE+ISEBUul/9WFWVhZGjRqFO3fuwNDQEG5ubtiwYQP69+8v1lMgIqo0gju44sS1VPx89jbGdKoPK2N2tSCqTEQfx0vbqss4IEREhREEAb0WncKFO2kY3bEexgc0FDsSUbmoLp/fol9qJCKi8iORSBDs97yv17qIeGQ8yxU5ERG9iIUXEVE183ZjG9StZYT0Z3nYHJkgdhwiegELLyKiakYqleDjDq4AgBUn4pCdpxQ5ERHlY+FFRFQN9fJwgJ2ZAskZ2dh+7q7YcYjoXyy8iIiqIQM9KYb71gUALDt2A0qVTt1HRVRpsfAiIqqmBrRyhHkNfcQ/eIL9l++LHYeIwMKLiKjaMpLrYaiPM4DnP56tY6MHEVVKLLyIiKqxYW2cYagvw9/30nHiWqrYcYh0HgsvIqJqzMLIAO+3rgPg+VkvIhIXCy8iompuuK8L9GUSRNx8gPMJj8SOQ6TTWHgREVVz9uaG6OXhAIBnvYjExsKLiEgHfNShLiQS4OA/SbienCF2HCKdxcKLiEgH1LM2wduNbQAAS47eFDkNke5i4UVEpCPyfzx7V/Rd3H38VOQ0RLqJhRcRkY7wcDRHG1dL5KkErDjBs15EYmDhRUSkQ4L9nv949s+Rt/EwK0fkNES6h4UXEZEOaVfPCk0dTPE0V4k1p+PFjkOkc1h4ERHpEIlEglH/9vVaezoeWdl5Iici0i0svIiIdExAE1u4WBkh7WkuNp1JEDsOkU5h4UVEpGNkUgk+7lAXALD8xE1k5ylFTkSkO1h4ERHpoN4tasPOTIHkjGxsi7ordhwincHCi4hIBxnoSTHC9/lZr6XHbiBPqRI5EZFuYOFFRKSjBrR2RE0jAyQ8fIK9l+6LHYdIJ7DwIiLSUTUM9PBBW2cAwOIjN6BSCeIGItIBLLyIiHTYEB9nGMv1EJuUgfAryWLHIar2WHgREekwM0N9DPFxAgAsPHIdgsCzXkQViYUXEZGO+6CtC+R6Uly4/RgRNx6IHYeoWmPhRUSk42qZyDGglSMAYNHR6yKnIareWHgRERFGtK8LPakEp64/wPmER2LHIaq2WHgRERFqW9RArxYOAIDFR2+InIao+mLhRUREAICPO7hCIgEO/ZOE2MQMseMQVUssvIiICABQz9oYXZraAgCWsK8XUYWoFIXXokWL4OzsDIVCAW9vb0RGRhbZdvny5fD19YWFhQUsLCzg7+9fbHsiIiq9UX71AAC7L9xDwoMnIqchqn5EL7y2bNmC0NBQTJ06FefOnYO7uzsCAgKQnFz4QH5Hjx7F+++/jyNHjiAiIgKOjo54++23cfcuf+SViOh1NXUwQ4cGtaASgKXH2deLqLxJBJFHy/P29karVq2wcOFCAIBKpYKjoyPGjBmDCRMmlLi8UqmEhYUFFi5ciMDAwBLbp6enw8zMDGlpaTA1NX3t/ERE1U1k3EP0WxYBA5kUJ77oCBtThdiRiKrN57eoZ7xycnIQFRUFf39/9TypVAp/f39ERESUah1PnjxBbm4uatasWVExiYh0SmuXmmjlbIEcpQorTtwUOw5RtSJq4ZWamgqlUgkbGxuN+TY2NkhMTCzVOr744gvY29trFG8vys7ORnp6usZERETFG9XxeV+vjWcS8CgrR+Q0RNWH6H28Xsd3332Hn3/+GTt27IBCUfip8JkzZ8LMzEw9OTo6ajklEVHV49egFhrbmeJJjhJrTseLHYeo2hC18LKysoJMJkNSUpLG/KSkJNja2ha77Jw5c/Ddd9/h4MGDaN68eZHtwsLCkJaWpp5u375dLtmJiKoziUSCkH/Peq05HY/M7DyRExFVD6IWXgYGBvD09ER4eLh6nkqlQnh4OHx8fIpc7vvvv8fXX3+NAwcOwMvLq9htyOVymJqaakxERFSyzk1tUdfKCGlPc7HpzC2x4xBVC6JfagwNDcXy5cuxdu1axMTEIDg4GFlZWQgKCgIABAYGIiwsTN1+1qxZmDx5MlatWgVnZ2ckJiYiMTERmZmZYj0FIqJqSSaV4GM/VwDA8hNxeJarFDkRUdUneuHVv39/zJkzB1OmTIGHhweio6Nx4MABdYf7hIQE3L9/X91+yZIlyMnJQd++fWFnZ6ee5syZI9ZTICKqtnp5OMDeTIGUjGxsO3dH7DhEVZ7o43hpW3UZB4SISFvWnIrDtN/+gWNNQxz5zA96MtG/s5MOqi6f33z3EBFRsfq3qgNLIwPcfvgUey7eL3kBIioSCy8iIiqWoYEMH7RzAQAsPnodKpVOXSghKlcsvIiIqERDfJxgItfD1aRMHI5JKnkBIioUCy8iIiqRqUIfgW2cAACLjt6AjnUPJio3LLyIiKhUgtq6QKEvxYXbj3H6xgOx4xBVSSy8iIioVKyM5RjQqg4AYNGR6yKnIaqaWHgREVGpjWxfF3pSCU7feIBzCY/EjkNU5bDwIiKiUrM3N8S7LR0AAIuP3BA5DVHVw8KLiIjK5OMOrpBIgMMxSbiSmC52HKIqhYUXERGVSd1axujazA4AsOQoz3oRlQULLyIiKrNR//549m8X7uHWgyyR0xBVHSy8iIiozJrYm6Fjw1pQCcDSYzfFjkNUZbDwIiKiVxLSsR4AYFvUHSSmPRM5DVHVwMKLiIheiZdzTbR2qYkcpQorTvCsF1FpsPAiIqJXln/Wa+OZBDzKyhE5DVHlx8KLiIheWfv6VmjqYIqnuUqsPh0vdhyiSo+FFxERvTKJRIIQv+dnvdacikNmdp7IiYgqNxZeRET0WgKa2MK1lhHSn+Vh45+3xI5DVKmx8CIiotcilUoQ/O9Zr+Un4vAsVylyIqLKi4UXERG9tnc87OFgbojUzGz8GnVH7DhElRYLLyIiem36Mik+6lAXALDs2A3kKlUiJyKqnFh4ERFRuejn5QgrYwPcefQUv124J3YcokqJhRcREZULhb4MH7Z7ftZr8dEbUKkEkRMRVT4svIiIqNwMfqMOTBR6uJ6ciYP/JIkdh6jSYeFFRETlxkShj2FtnAEAi49ehyDwrBfRi1h4ERFRuQpq6wJDfRku3knDyeupYschqlRYeBERUbmqaWSA91vXAQAsOnJd5DRElQsLLyIiKncj2rtAXybBnzcfIurWQ7HjEFUaLLyIiKjc2ZkZok/L2gCAxUduiJyGqPJg4UVERBXiow6ukEqA8CvJ+OdeuthxiCoFFl5ERFQhXKyM0K25PQBgyTGe9SICWHgREVEFGuXnCgDYe/Ee4lKzRE5DJD4WXkREVGEa2Zmik5s1VMLz33Ak0nWiF16LFi2Cs7MzFAoFvL29ERkZWWTbv//+G3369IGzszMkEgnmz5+vvaBERPRKRnWsBwDYdu4O7qc9FTkNkbhELby2bNmC0NBQTJ06FefOnYO7uzsCAgKQnJxcaPsnT56gbt26+O6772Bra6vltERE9Co8nSzwRt2ayFUKWH48Tuw4RKIStfCaN28eRowYgaCgIDRu3BhLly5FjRo1sGrVqkLbt2rVCrNnz8aAAQMgl8u1nJaIiF5VyL9nvTZHJuBBZrbIaYjEI1rhlZOTg6ioKPj7+/8XRiqFv78/IiIiym072dnZSE9P15iIiEi72tWzQvPaZniaq8Sa0/FixyESjWiFV2pqKpRKJWxsbDTm29jYIDExsdy2M3PmTJiZmaknR0fHcls3ERGVjkQiwSi/52e91pyOR8azXJETEYlD9M71FS0sLAxpaWnq6fbt22JHIiLSSW83tkE9a2NkPMvDhj8TxI5DJArRCi8rKyvIZDIkJSVpzE9KSirXjvNyuRympqYaExERaZ9UKkFwh+fjeq08eRPPcpUiJyLSPtEKLwMDA3h6eiI8PFw9T6VSITw8HD4+PmLFIiKiCtTTwx4O5oZIzczBL3/xCgTpHlEvNYaGhmL58uVYu3YtYmJiEBwcjKysLAQFBQEAAgMDERYWpm6fk5OD6OhoREdHIycnB3fv3kV0dDSuX78u1lMgIqIy0JdJ8XGHugCAZcduIlepEjkRkXbpibnx/v37IyUlBVOmTEFiYiI8PDxw4MABdYf7hIQESKX/1Yb37t1DixYt1H/PmTMHc+bMQYcOHXD06FFtxyciolfwnpcj/hd+HXcfP8Wu6Hvo61lb7EhEWiMRBEEQO4Q2paenw8zMDGlpaezvRUQkkqXHbuC7/VdQt5YRDn3aATKpROxIVMlVl8/van9XIxERVT6DvOvAVKGHmylZOPh3+Q0hRFTZsfAiIiKtM1HoY1gbZwDAoqPXoWMXX0iHsfAiIiJRDGvrAkN9GS7fTcfxa6lixyHSChZeREQkippGBhjoXQcAsOgI704n3cDCi4iIRDPCty70ZRJExj3E2fiHYschqnAsvIiISDS2Zgr1cBKLedaLdAALLyIiEtVH7V0hlQBHYlPw9700seMQVSgWXkREJCpnKyN0b24PAFh89IbIaYgqFgsvIiISXbDf8x/P3nfpPm6mZIqchqjisPAiIiLRNbIzhX8jawjC81HtiaorFl5ERFQpjOpYDwCw/dxd3H38VOQ0RBWDhRcREVUKLetYwKeuJfJUApYfvyl2HKIKwcKLiIgqjZB/z3r9fDYBqZnZIqchKn8svIiIqNJoW88S7rXN8CxXhdWn4sSOQ1TuWHgREVGlIZFI1H291p2+hfRnuSInIipfLLyIiKhSeauRDepbGyMjOw/rI26JHYeoXLHwIiKiSkUqlWBUx+fjeq06GYenOUqRExGVHxZeRERU6fRobo/aFoZ4kJWDLWcTxI5DVG5YeBERUaWjJ5Pi4w7Pz3r9dPwmcvJUIiciKh8svIiIqFLq61kbtUzkuJf2DDuj74odh6hcsPAiIqJKSaEvwwhfFwDA0qM3oFQJIicien0svIiIqNIa6O0EM0N93EzNwoHLiWLHIXptLLyIiKjSMpbrYVgbZwDAoiPXIQg860VVGwsvIiKq1Ia1cUYNAxn+uZ+Oo1dTxI5D9FpYeBERUaVmYWSAQd51AACLj1wXOQ3R62HhRUREld5w37owkElxNv4RIuMeih2H6JWx8CIiokrPxlSBvl61ATzv60VUVbHwIiKiKuHj9q6QSoBjV1Nw+W6a2HGIXgkLLyIiqhLqWNZAT3d7AMDiozzrRVUTCy8iIqoygv3qAQD2X07E9eRMkdMQlR0LLyIiqjIa2prgrcY2EARg6bEbYschKjMWXkREVKWM8nv+49k7z9/FnUdPRE5DVDYsvIiIqEppUccCbetZIk8lYPnxm2LHISqTSlF4LVq0CM7OzlAoFPD29kZkZGSx7X/99Ve4ublBoVCgWbNm2Ldvn5aSEhFRZRDyb1+vn8/eRkpGtshpiEpP9MJry5YtCA0NxdSpU3Hu3Dm4u7sjICAAycnJhbY/ffo03n//fXz44Yc4f/48evXqhV69euHy5ctaTk5ERGLxcbWEh6M5svNUWHUqTuw4RKUmEUT+xVFvb2+0atUKCxcuBACoVCo4OjpizJgxmDBhQoH2/fv3R1ZWFvbs2aOe98Ybb8DDwwNLly4tcXvp6ekwMzNDWloaTE1Ny++JEBGRVh36Jwkj1v0FY7kedoa0gUJfJnYk+peBnhTWJopyXWd1+fzWE3PjOTk5iIqKQlhYmHqeVCqFv78/IiIiCl0mIiICoaGhGvMCAgKwc+fOQttnZ2cjO/u/09Dp6emvH5yIiETXyc0aDW1MEJuUAf95x8WOQy9oWccc20e1FTtGpSRq4ZWamgqlUgkbGxuN+TY2Nrhy5UqhyyQmJhbaPjExsdD2M2fOxPTp08snMBERVRpSqQRhXd0wbks0nuYoxY5DL9CXid6TqdIStfDShrCwMI0zZOnp6XB0dBQxERERlRe/htaInvK22DGISk3UwsvKygoymQxJSUka85OSkmBra1voMra2tmVqL5fLIZfLyycwERER0WsQ9VyggYEBPD09ER4erp6nUqkQHh4OHx+fQpfx8fHRaA8Ahw4dKrI9ERERUWUh+qXG0NBQDB06FF5eXmjdujXmz5+PrKwsBAUFAQACAwPh4OCAmTNnAgDGjh2LDh06YO7cuejWrRt+/vln/PXXX/jpp5/EfBpEREREJRK98Orfvz9SUlIwZcoUJCYmwsPDAwcOHFB3oE9ISIBU+t+JuTZt2mDTpk348ssvMXHiRNSvXx87d+5E06ZNxXoKRERERKUi+jhe2lZdxgEhIiLSJdXl85v3exIRERFpCQsvIiIiIi1h4UVERESkJSy8iIiIiLSEhRcRERGRlrDwIiIiItISFl5EREREWsLCi4iIiEhLWHgRERERaYnoPxmkbfkD9aenp4uchIiIiEor/3O7qv/gjs4VXhkZGQAAR0dHkZMQERFRWWVkZMDMzEzsGK9M536rUaVS4d69ezAxMYFEIinXdaenp8PR0RG3b9+u0r8jVV3w9ahc+HpULnw9Kh++JsUTBAEZGRmwt7eHVFp1e0rp3BkvqVSK2rVrV+g2TE1N+aapRPh6VC58PSoXvh6VD1+TolXlM135qm7JSERERFTFsPAiIiIi0hIWXuVILpdj6tSpkMvlYkch8PWobPh6VC58PSofvia6Qec61xMRERGJhWe8iIiIiLSEhRcRERGRlrDwIiIiItISFl5EREREWsLCq5wsWrQIzs7OUCgU8Pb2RmRkpNiRdNbMmTPRqlUrmJiYwNraGr169UJsbKzYsehf3333HSQSCcaNGyd2FJ119+5dDB48GJaWljA0NESzZs3w119/iR1LJymVSkyePBkuLi4wNDSEq6srvv766yr/e4RUNBZe5WDLli0IDQ3F1KlTce7cObi7uyMgIADJycliR9NJx44dQ0hICP78808cOnQIubm5ePvtt5GVlSV2NJ139uxZLFu2DM2bNxc7is569OgR2rZtC319fezfvx///PMP5s6dCwsLC7Gj6aRZs2ZhyZIlWLhwIWJiYjBr1ix8//33+PHHH8WORhWEw0mUA29vb7Rq1QoLFy4E8Pz3IB0dHTFmzBhMmDBB5HSUkpICa2trHDt2DO3btxc7js7KzMxEy5YtsXjxYnzzzTfw8PDA/PnzxY6lcyZMmIBTp07hxIkTYkchAN27d4eNjQ1WrlypntenTx8YGhpiw4YNIiajisIzXq8pJycHUVFR8Pf3V8+TSqXw9/dHRESEiMkoX1paGgCgZs2aIifRbSEhIejWrZvGe4W0b/fu3fDy8sJ7770Ha2trtGjRAsuXLxc7ls5q06YNwsPDcfXqVQDAhQsXcPLkSXTp0kXkZFRRdO5HsstbamoqlEolbGxsNObb2NjgypUrIqWifCqVCuPGjUPbtm3RtGlTsePorJ9//hnnzp3D2bNnxY6i827evIklS5YgNDQUEydOxNmzZ/HJJ5/AwMAAQ4cOFTuezpkwYQLS09Ph5uYGmUwGpVKJb7/9FoMGDRI7GlUQFl5UrYWEhODy5cs4efKk2FF01u3btzF27FgcOnQICoVC7Dg6T6VSwcvLCzNmzAAAtGjRApcvX8bSpUtZeIngl19+wcaNG7Fp0yY0adIE0dHRGDduHOzt7fl6VFMsvF6TlZUVZDIZkpKSNOYnJSXB1tZWpFQEAKNHj8aePXtw/Phx1K5dW+w4OisqKgrJyclo2bKlep5SqcTx48excOFCZGdnQyaTiZhQt9jZ2aFx48Ya8xo1aoRt27aJlEi3ff7555gwYQIGDBgAAGjWrBlu3bqFmTNnsvCqptjH6zUZGBjA09MT4eHh6nkqlQrh4eHw8fERMZnuEgQBo0ePxo4dO/DHH3/AxcVF7Eg6rVOnTrh06RKio6PVk5eXFwYNGoTo6GgWXVrWtm3bAsOrXL16FU5OTiIl0m1PnjyBVKr5USyTyaBSqURKRBWNZ7zKQWhoKIYOHQovLy+0bt0a8+fPR1ZWFoKCgsSOppNCQkKwadMm7Nq1CyYmJkhMTAQAmJmZwdDQUOR0usfExKRA/zojIyNYWlqy350IPv30U7Rp0wYzZsxAv379EBkZiZ9++gk//fST2NF0Uo8ePfDtt9+iTp06aNKkCc6fP4958+bhgw8+EDsaVRAOJ1FOFi5ciNmzZyMxMREeHh5YsGABvL29xY6lkyQSSaHzV69ejWHDhmk3DBXKz8+Pw0mIaM+ePQgLC8O1a9fg4uKC0NBQjBgxQuxYOikjIwOTJ0/Gjh07kJycDHt7e7z//vuYMmUKDAwMxI5HFYCFFxEREZGWsI8XERERkZaw8CIiIiLSEhZeRERERFrCwouIiIhIS1h4EREREWkJCy8iIiIiLWHhRURERKQlLLyISOucnZ3LPHjqsGHD0KtXL/Xffn5+GDduXLnmqggSiQQ7d+4UOwYRVRIsvIh00LBhwyCRSNSTpaUlOnfujIsXL4odrdS2b9+Or7/+WuwYJbp//z66dOkidgwiqiRYeBHpqM6dO+P+/fu4f/8+wsPDoaenh+7du4sdq9Rq1qwJExMTsWOUyNbWFnK5XOwYRFRJsPAi0lFyuRy2trawtbWFh4cHJkyYgNu3byMlJUXd5tKlS3jzzTdhaGgIS0tLjBw5EpmZmerH8y//zZkzB3Z2drC0tERISAhyc3PVbZKTk9GjRw8YGhrCxcUFGzduLDGbUqlEaGgozM3NYWlpif/7v//Dy79u9vKlRmdnZ3zzzTcIDAyEsbExnJycsHv3bqSkpOCdd96BsbExmjdvjr/++ktjPSdPnoSvry8MDQ3h6OiITz75BFlZWRrrnTFjBj744AOYmJigTp06Gj8onZOTg9GjR8POzg4KhQJOTk6YOXOm+vGXLzWWxz4loqqLhRcRITMzExs2bEC9evVgaWkJAMjKykJAQAAsLCxw9uxZ/Prrrzh8+DBGjx6tseyRI0dw48YNHDlyBGvXrsWaNWuwZs0a9ePDhg3D7du3ceTIEWzduhWLFy9GcnJysXnmzp2LNWvWYNWqVTh58iQePnyIHTt2lPg8fvjhB7Rt2xbnz59Ht27dMGTIEAQGBmLw4ME4d+4cXF1dERgYqC7ibty4gc6dO6NPnz64ePEitmzZgpMnTxZ4jnPnzoWXlxfOnz+PUaNGITg4GLGxsQCABQsWYPfu3fjll18QGxuLjRs3wtnZudB85bVPiagKE4hI5wwdOlSQyWSCkZGRYGRkJAAQ7OzshKioKHWbn376SbCwsBAyMzPV8/bu3StIpVIhMTFRvR4nJychLy9P3ea9994T+vfvLwiCIMTGxgoAhMjISPXjMTExAgDhhx9+KDKfnZ2d8P3336v/zs3NFWrXri2888476nkdOnQQxo4dq/7byclJGDx4sPrv+/fvCwCEyZMnq+dFREQIAIT79+8LgiAIH374oTBy5EiNbZ84cUKQSqXC06dPC12vSqUSrK2thSVLlgiCIAhjxowR3nzzTUGlUhX6XAAIO3bsEAShfPYpEVVtPONFpKM6duyI6OhoREdHIzIyEgEBAejSpQtu3boFAIiJiYG7uzuMjIzUy7Rt2xYqlUp9tgcAmjRpAplMpv7bzs5OfUYrJiYGenp68PT0VD/u5uYGc3PzInOlpaXh/v378Pb2Vs/T09ODl5dXic+pefPm6n/b2NgAAJo1a1ZgXn6+CxcuYM2aNTA2NlZPAQEBUKlUiIuLK3S9EokEtra26nUMGzYM0dHRaNiwIT755BMcPHiwyHzlsU+JqGrTEzsAEYnDyMgI9erVU/+9YsUKmJmZYfny5fjmm29KvR59fX2NvyUSCVQqVbnlLIsXs0gkkiLn5efLzMzERx99hE8++aTAuurUqVPoevPXk7+Oli1bIi4uDvv378fhw4fRr18/+Pv7Y+vWreXyPF7eHhFVbTzjRUQAnn+4S6VSPH36FADQqFEjXLhwQaOj+alTpyCVStGwYcNSrdPNzQ15eXmIiopSz4uNjcXjx4+LXMbMzAx2dnY4c+aMet7L6ygvLVu2xD///IN69eoVmAwMDEq9HlNTU/Tv3x/Lly/Hli1bsG3bNjx8+LBAu/LYp0RUtbHwItJR2dnZSExMRGJiImJiYjBmzBhkZmaiR48eAIBBgwZBoVBg6NChuHz5Mo4cOYIxY8ZgyJAh6kt2JWnYsCE6d+6Mjz76CGfOnEFUVBSGDx8OQ0PDYpcbO3YsvvvuO+zcuRNXrlzBqFGjii3WXtUXX3yB06dPY/To0YiOjsa1a9ewa9euAp3dizNv3jxs3rwZV65cwdWrV/Hrr7/C1ta20Mup5bFPiahqY+FFpKMOHDgAOzs72NnZwdvbW32XnZ+fHwCgRo0a+P333/Hw4UO0atUKffv2RadOnbBw4cIybWf16tWwt7dHhw4d8O6772LkyJGwtrYudpnPPvsMQ4YMwdChQ+Hj4wMTExP07t37VZ9qkZo3b45jx47h6tWr8PX1RYsWLTBlyhTY29uXeh0mJib4/vvv4eXlhVatWiE+Ph779u2DVFrw8Fpe+5SIqi6JILw0OA4RERERVQie8SIiIiLSEhZeRERERFrCwouIiIhIS1h4EREREWkJCy8iIiIiLWHhRURERKQlLLyIiIiItISFFxEREZGWsPAiIiIi0hIWXkRERERawsKLiIiISEtYeBERERFpyf8DHjpBZkZpV+0AAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "TVD = []\n", "for chi in range(1,11):\n", " mps.set_cutoff(chi)\n", " mps.set_input_state(input_state)\n", " probs_mps = mps.prob_distribution()\n", " TVD.append(tvd(probs_mps, probs_slos))\n", "\n", "plt.plot(TVD)\n", "plt.xlabel('Bond dimension')\n", "plt.ylabel('TVD')\n", "plt.title('Evolution of the TVD between SLOS and MPS probability distributions');" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the TVD decreases as the size of the matrices increases, untill reaching 0 for a bond dimension of 7." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "\n", "> [1] Ulrich Schollwöck. The density-matrix renormalization group in the age of matrix product states. [Annals of Physics](https://doi.org/10.1016/j.aop.2010.09.012), 326(1):96–192, jan 2011.\n", "\n", "> [2] Changhun Oh, Kyungjoo Noh, Bill Fefferman, and Liang Jiang. Classical simulation of lossy boson sampling using matrix product operators. [Physical Review A](https://doi.org/10.1103/PhysRevA.104.022407), 104(2), aug 2021.\n", "\n", "> [3] Hui Wang, et al. Boson Sampling with 20 Input Photons and a 60-Mode Interferometer in a $10^{14}$-Dimensional Hilbert Space. [Physical Review Letters](https://link.aps.org/doi/10.1103/PhysRevLett.123.250503), 123(25):250503, December 2019. Publisher: American Physical Society." ] } ], "metadata": { "language_info": { "name": "python" } }, "nbformat": 4, "nbformat_minor": 2 }