{ "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 choosing 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", "Φ=0.682353\n", "\n", "\n", "Φ=1.839246\n", "\n", "\n", "Φ=5.325111\n", "\n", "\n", "Φ=1.540778\n", "\n", "\n", "Φ=6.134872\n", "\n", "\n", "Φ=5.133685\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.077922\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=3.115285\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.773545\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.310632\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.118361\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.116664\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.298164\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.383843\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.603772\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.871176\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.91304\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.896901\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.873668\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.198433\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.864696\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.929913\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.252653\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=3.269561\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.223128\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.226287\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.467155\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.570584\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.129614\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=3.067041\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=1.169704\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=3.33072\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.346625\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=0.247317\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=5.549437\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "Rx\n", "\n", "\n", "Φ=4.853279\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=pcvl.InterferometerShape.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,0,0,1,2> 0.078508
|0,0,0,0,0,3> 0.053524
|2,1,0,0,0,0> 0.051962
|0,0,0,0,2,1> 0.042042
|3,0,0,0,0,0> 0.036251
|0,0,1,1,0,1> 0.033647
|1,0,1,0,1,0> 0.032665
|1,1,0,0,0,1> 0.032312
|1,2,0,0,0,0> 0.028953
|2,0,0,0,1,0> 0.028039
|0,1,2,0,0,0> 0.027845
|0,1,1,1,0,0> 0.026747
|0,0,1,0,2,0> 0.02671
|0,0,0,0,3,0> 0.026532
|1,0,1,0,0,1> 0.026327
|1,0,0,1,1,0> 0.024339
|0,0,2,1,0,0> 0.023886
|1,0,0,0,0,2> 0.023375
|1,1,0,1,0,0> 0.021942
|1,0,0,1,0,1> 0.021592
" ], "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 choosing 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,0,0,1,2> 0.078508
|0,0,0,0,0,3> 0.053524
|2,1,0,0,0,0> 0.051962
|0,0,0,0,2,1> 0.042042
|3,0,0,0,0,0> 0.036251
|0,0,1,1,0,1> 0.033647
|1,0,1,0,1,0> 0.032665
|1,1,0,0,0,1> 0.032312
|1,2,0,0,0,0> 0.028953
|2,0,0,0,1,0> 0.028039
|0,1,2,0,0,0> 0.027845
|0,1,1,1,0,0> 0.026747
|0,0,1,0,2,0> 0.02671
|0,0,0,0,3,0> 0.026532
|1,0,1,0,0,1> 0.026327
|1,0,0,1,1,0> 0.024339
|0,0,2,1,0,0> 0.023886
|1,0,0,0,0,2> 0.023375
|1,1,0,1,0,0> 0.021942
|1,0,0,1,0,1> 0.021592
" ], "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": "iVBORw0KGgoAAAANSUhEUgAAAl4AAAHHCAYAAABuoFaQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABijklEQVR4nO3deVwU9f8H8NfuAruAXIJcgqCAkhcoKqECHhSaR1qmlopaWimZRlZSv7RTM4/8mmeW5p1laaVpKpp3ahoepYgHigeXyqly7H5+f+BurtwKO7D7ej4e+3jI7Gdm3zvrzrx25jOfkQkhBIiIiIioxsmlLoCIiIjIVDB4ERERERkIgxcRERGRgTB4ERERERkIgxcRERGRgTB4ERERERkIgxcRERGRgTB4ERERERkIgxcRERGRgRht8JLJZPjggw+qdZnffvstZDIZkpKSqnW51W3GjBlo0qQJFAoFAgMDqzz/H3/8AZlMhvXr11d/cXXIBx98AJlMhoyMDKlLISPk7e2NESNGSF1GndKlSxe0bNmy2paXlJQEmUyGmTNnVthWuz2434OfoXbb+ccff1RbjVWhfT/ffvutblppddeULl26oEuXLrq/Db0vGTFiBLy9vQ3yWo+iRoOXNqiU9fjzzz9r8uUf2tSpU7Fx40apy3go27Ztw9tvv41OnTph2bJlmDp1aplt16xZgzlz5hiuuPt06dKl3P8b2kffvn0hk8nwf//3f2UuKzExETKZDDExMQD+29BoH1ZWVmjUqBH69OmDZcuWIT8/31BvswQp17mhJCUlYeTIkfDx8YFKpYKrqyvCwsIwZcoUvXaV3YneuHEDb731Fpo1awaVSoX69esjMjISmzZtKrV9eno6xo8fD39/f1haWsLZ2RkdOnTAO++8g9zc3Gp5j4ak/X88atSoUp9/7733dG3u/5EwYsQIve+Bra0tAgICMGvWrBLfgX379qFnz55o2LAhVCqV7vuyZs2aGn1vpqAufuevXbuGDz74APHx8VKXUkJtrq2yzAzxIh999BEaN25cYrqvr68hXr7Kpk6digEDBqBfv35604cNG4bBgwdDqVRKU1gl7Ny5E3K5HN988w0sLCzKbbtmzRqcOnUKEyZMMExx93nvvff0diRHjhzB3Llz8e677+Kxxx7TTW/dujUSExOxdu1afPLJJ6UuS7tzGDp0qN70hQsXol69esjPz8fVq1fx+++/48UXX8ScOXOwadMmeHp61sA7K5+U69wQzp07h/bt28PS0hIvvvgivL29cf36dRw7dgzTp0/Hhx9+WKXlJSQkoHv37khPT8fIkSPRrl07ZGZmYvXq1ejTpw8mTpyIGTNm6NrfvHkT7dq1Q3Z2Nl588UX4+/vjxo0bOHHiBBYuXIgxY8agXr161f22a5xKpcKPP/6IBQsWlPher127FiqVCnfv3i0xn1KpxNdffw0AyMzMxI8//oiJEyfiyJEj+O677wAAP/zwAwYNGoTAwECMHz8eDg4OuHjxIvbs2YMlS5bghRdeqPk3WAf83//9HyZNmlRum7CwMNy5c0fvM5L6O1+Zuh907do1fPjhh/D29q7SWZNt27ZVsbqqK6+2JUuWQKPR1HgNj8ogwatnz55o166dIV6qRikUCigUCqnLKFdaWhosLS0rDF1Se+KJJ/T+VqlUmDt3Lp544gm9Q9UAMGTIELz//vv4888/8fjjj5dY1tq1a+Hv74+2bdvqTR8wYACcnJx0f0+ePBmrV69GVFQUnnvuuVp7xLUu++KLL5Cbm4v4+Hh4eXnpPZeWllalZRUWFmLAgAG4desW9uzZg+DgYN1zb7zxBoYMGYKZM2eiXbt2GDRoEADgm2++weXLl7F//3507NhRb3nZ2dm1/ntRlh49euCXX37Bli1b8PTTT+umHzhwABcvXsSzzz6LH3/8scR8ZmZmej9Ixo4di+DgYKxbtw6zZ8+Gu7s7PvjgAzRv3hx//vlnifVT1c+sqvLy8mBtbV2jr1FdzMzMYGZW/i5TLpdDpVIZqKLKqUzdj+r27duwsrKS/Ptlbm4u6etXluR9vAoLC1G/fn2MHDmyxHPZ2dlQqVSYOHGiblpaWhpeeukluLi4QKVSISAgAMuXL6/wdco69/vg+W+ZTIa8vDwsX75cd4heew6/rD5eCxYsQIsWLaBUKuHu7o7o6GhkZmbqtdGeVvn333/RtWtXWFlZoWHDhvj8888rrB0AioqK8PHHH8PHxwdKpRLe3t5499139U4ZyGQyLFu2DHl5ebra7z/X/2A9mzdvxqVLl3RtH1w/Go0Gn376KTw8PKBSqdC9e3ecO3euxLIOHTqEHj16wM7ODlZWVggPD8f+/fsr9b4qY8iQIQBQ6mmPo0ePIiEhQdemMssaNWoUDh06hO3bt1dqnoyMDAwcOBC2trZwdHTE+PHjSz26sGrVKgQFBcHS0hL169fH4MGDkZycrHu+rHUuhICTk5PuVClQvO7t7e2hUCj0/i9Nnz4dZmZmeqfMzpw5gwEDBqB+/fpQqVRo164dfvnllxL1ZWZmYsKECfD09IRSqYSvry+mT5+u9wvx/j4vX331le7/W/v27XHkyJEK19X58+fh4eFRInQBgLOzc4Xz3+/HH3/EqVOnMGnSJL3QBRT/CFq8eDHs7e31+nKeP38eCoWi1IBua2tb4U7x0qVLGDt2LJo1awZLS0s4OjriueeeK/Gd124L9u/fj5iYGDRo0ADW1tbo378/0tPT9doKIfDJJ5/Aw8MDVlZW6Nq1K/75558qrYuGDRsiLCysxHdg9erVaNWqVaX7Pcnlct0PG+17On/+PNq3b1/qTrMyn5m3tzd69+6Nbdu2ITAwECqVCs2bN8dPP/2k1067znbv3o2xY8fC2dkZHh4euucrsx3VOnr0KDp27AhLS0s0btwYixYt0nu+oKAAkydPRlBQEOzs7GBtbY3Q0FDs2rWrzPfxxRdfwMvLC5aWlggPD8epU6f0nq9MX6kH+3iV9Z3Pzc2FtbU1xo8fX2IZV65cgUKhwLRp08p9rczMTIwYMQJ2dnawt7fH8OHDS11fpdW9fft2dO7cGfb29qhXrx6aNWuGd999V/ce2rdvDwAYOXJkiX2Jdl929OhRhIWFwcrKSjfvg328tNRqNd599124urrC2toaffv21ds2AmX3ebx/mRXVVtp+Pi8vD2+++aZuu9esWTPMnDkTQgi9djKZDK+99ho2btyIli1bQqlUokWLFti6dateu5ycHEyYMAHe3t5QKpVwdnbGE088gWPHjpWovSwGCV5ZWVnIyMjQe9y4cQNAcULt378/Nm7ciIKCAr35Nm7ciPz8fAwePBgAcOfOHXTp0gUrV67EkCFDMGPGDNjZ2WHEiBH43//+Vy21rly5EkqlEqGhoVi5ciVWrlyJV155pcz2H3zwAaKjo+Hu7o5Zs2bh2WefxeLFi/Hkk0+isLBQr+2tW7fQo0cPXT8Lf39/vPPOO9iyZUuFdY0aNQqTJ09G27Zt8cUXXyA8PBzTpk3TrRtt7aGhoVAqlbraw8LCSl3ee++9h8DAQDg5OenaPtgP4bPPPsOGDRswceJExMbG4s8//ywRcHbu3ImwsDBkZ2djypQpmDp1KjIzM9GtWzccPny4wvdVGY0bN0bHjh3x/fffQ61W6z2n3RFV5XTIsGHDAFT+sPjAgQNx9+5dTJs2DU899RTmzp2Ll19+Wa/Np59+iqioKPj5+WH27NmYMGEC4uLiEBYWptsYlrXOZTIZOnXqhD179uiWd+LECWRlZQGAXojdu3cv2rRpoztd9s8//+Dxxx/H6dOnMWnSJMyaNQvW1tbo168fNmzYoJvv9u3bCA8Px6pVqxAVFYW5c+eiU6dOiI2N1Qt8WmvWrMGMGTPwyiuv4JNPPkFSUhKeeeaZEv+nH+Tl5YXk5GTs3LmzUuu2PL/++isAICoqqtTn7ezs8PTTT+PMmTO6HwReXl5Qq9VYuXLlQ73mkSNHcODAAQwePBhz587Fq6++iri4OHTp0gW3b98u0X7cuHE4fvw4pkyZgjFjxuDXX3/Fa6+9ptdm8uTJeP/99xEQEKC78OXJJ59EXl5elWp74YUX8Ouvv+pCd1FREX744Ycqnwo8f/48AMDR0RFA8TqLi4vDlStXqrSc+yUmJmLQoEHo2bMnpk2bBjMzMzz33HOl/rgZO3Ys/v33X0yePFl3Cqyq29GnnnoKQUFB+Pzzz+Hh4YExY8Zg6dKlujbZ2dn4+uuv0aVLF0yfPh0ffPAB0tPTERkZWWrfoBUrVmDu3LmIjo5GbGwsTp06hW7duiE1NfWh1wlQ9ne+Xr166N+/P9atW1dim7Z27VoIIcr9MSmEwNNPP42VK1di6NCh+OSTT3DlyhUMHz68wpr++ecf9O7dG/n5+fjoo48wa9Ys9O3bV7edeeyxx/DRRx8BAF5++eVS9yU3btxAz549ERgYiDlz5qBr167lvuann36KzZs345133sHrr7+O7du3IyIiAnfu3Kmw3vtVprb7CSHQt29ffPHFF+jRowdmz56NZs2a4a233ip1u7dv3z6MHTsWgwcPxueff467d+/i2Wef1eUVAHj11VexcOFCPPvss1iwYAEmTpwIS0tLnD59uvJvRNSgZcuWCQClPpRKpa7d77//LgCIX3/9VW/+p556SjRp0kT395w5cwQAsWrVKt20goICERISIurVqyeys7N10wGIKVOm6P4ePny48PLyKlHjlClTxIOrwdraWgwfPrzM93Px4kUhhBBpaWnCwsJCPPnkk0KtVuvazZs3TwAQS5cu1U0LDw8XAMSKFSt00/Lz84Wrq6t49tlnS7zW/eLj4wUAMWrUKL3pEydOFADEzp079d6ntbV1ucvT6tWrV6nrZNeuXQKAeOyxx0R+fr5u+v/+9z8BQJw8eVIIIYRGoxF+fn4iMjJSaDQaXbvbt2+Lxo0biyeeeKJSdQghxA8//CAAiF27dpX6/Pz58wUA8fvvv+umqdVq0bBhQxESEqLXVvuZpqenl7qsW7duCQCif//+5dakXU7fvn31po8dO1YAEMePHxdCCJGUlCQUCoX49NNP9dqdPHlSmJmZ6U0va53PmDFDKBQK3f/huXPnCi8vL9GhQwfxzjvv6N6vvb29eOONN3Tzde/eXbRq1UrcvXtXN02j0YiOHTsKPz8/3bSPP/5YWFtbi7Nnz+q97qRJk4RCoRCXL18WQghx8eJFAUA4OjqKmzdv6tr9/PPPpX5HH3Tq1ClhaWkpAIjAwEAxfvx4sXHjRpGXl1eibXh4uGjRokWZywoMDBR2dnblvt7s2bMFAPHLL78IIYRISUkRDRo0EACEv7+/ePXVV8WaNWtEZmZmucvRun37dolpBw8eLPHd1W4LIiIi9P7vv/HGG0KhUOheT7uN6NWrl167d999VwAodTvzIAAiOjpa3Lx5U1hYWIiVK1cKIYTYvHmzkMlkIikpqdT/89ptQXp6ukhPTxfnzp0TU6dOFTKZTLRu3VrX7ptvvhEAhIWFhejatat4//33xd69e/W2aeXx8vISAMSPP/6om5aVlSXc3NxEmzZtSqyzzp07i6KiIt30h9mOzpo1SzctPz9fBAYGCmdnZ1FQUCCEEKKoqEhv2yVE8ffexcVFvPjii7pp2v/vlpaW4sqVK7rphw4dEgD0vmul7Su8vLz0PkPttvP+7VhZ33ntfm/Lli1601u3bi3Cw8NLtL/fxo0bBQDx+eef66YVFRWJ0NBQAUAsW7aszLq/+OKLcrePQghx5MiREsvR0n4GixYtKvW5+2vXro+GDRvq7Z+///57AUD873//0017cF2WtczyantwP69dT5988oleuwEDBgiZTCbOnTunm6b9Dtw/7fjx4wKA+PLLL3XT7OzsRHR0dInXrgqDHPGaP38+tm/frve4/yhPt27d4OTkhHXr1umm3bp1C9u3b9f13QCA3377Da6urnj++ed108zNzfH6668jNzcXu3fvNsTb0dmxYwcKCgowYcIEyOX/rcrRo0fD1tYWmzdv1mtfr149vf4WFhYW6NChAy5cuFDu6/z2228AUCKhv/nmmwBQ4nWqy8iRI/VOP4SGhgKArt74+HgkJibihRdewI0bN3RHM/Py8tC9e3fs2bOn2jo6Dho0CObm5nqnWnbv3o2rV69W+jSjlvZoUU5OTqXaR0dH6/09btw4AP99Lj/99BM0Gg0GDhyod1TX1dUVfn5+5Z7e0AoNDYVarcaBAwcAFB/ZCg0NRWhoKPbu3QsAOHXqFDIzM3Wfw82bN7Fz504MHDgQOTk5ekeTIyMjkZiYiKtXrwIo7kAdGhoKBwcHvRojIiKgVqv1jrYBxevbwcFBrz4AFf5fbdGiBeLj4zF06FAkJSXhf//7H/r16wcXFxcsWbKkwvVwv5ycHNjY2JTbRvt8dnY2AMDFxQXHjx/Hq6++ilu3bmHRokV44YUX4OzsjI8//rjE6YUHWVpa6v5dWFiIGzduwNfXF/b29qWeSnj55Zf1TuNoP8dLly4B+G8bMW7cOL12D9PR2sHBAT169MDatWsBFB+V7NixY6mndbXy8vLQoEEDNGjQAL6+vnj33XcREhKidzT0xRdfxNatW9GlSxfs27cPH3/8MUJDQ+Hn56f7/1gRd3d39O/fX/e3ra0toqKi8PfffyMlJUWv7ejRo/X6ylZ1O2pmZqZ3FsLCwgKvvPIK0tLScPToUQDFp6K12y6NRoObN2+iqKgI7dq1K/Vz7NevHxo2bKj7u0OHDggODtZ9x2tCREQE3N3dsXr1at20U6dO4cSJEyUuFHrQb7/9BjMzM4wZM0Y3TaFQ6LZN5bG3twcA/Pzzzw+9fVYqlaV2DypLVFSU3nd5wIABcHNzq9H1CxSvJ4VCgddff11v+ptvvgkhRImzTREREfDx8dH93bp1a9ja2upt9+zt7XHo0CFcu3btoesySPDq0KEDIiIi9B73H5o0MzPDs88+i59//lnXZ+mnn35CYWGhXvC6dOkS/Pz89L6cAHRXwWk3doaifb1mzZrpTbewsECTJk1K1OPh4VHiXLuDgwNu3bpV4evI5fISV4G6urrC3t6+xt53o0aN9P7W7oi19SYmJgIAhg8frtu4ax9ff/018vPzdafLHpWjoyMiIyOxYcMGXf+qNWvWwMzMDAMHDqzSsrSnairaqWv5+fnp/e3j4wO5XK7rI5OYmAghBPz8/Eqsh9OnT1eqg3Lbtm1hZWWlC1na4BUWFoa//voLd+/e1T3XuXNnAMVXEAoh8P7775d4Xe3QDdrXTkxMxNatW0u0i4iI0GunVdFnX56mTZti5cqVyMjIwIkTJzB16lSYmZnh5Zdfxo4dOyqcX8vGxqbCcKx9/v7P0s3NDQsXLsT169eRkJCAuXPnokGDBpg8eTK++eabcpd3584dTJ48WdcfxMnJCQ0aNEBmZmap/5crWk/a7+aD/4caNGigF2wr64UXXsD27dtx+fJlbNy4scLTjCqVSvdjd8+ePUhOTsb+/fvRpEkTvXaRkZH4/fffkZmZiT179iA6OhqXLl1C7969K/X/19fXt8S2rWnTpgBQon/cg1e4V3U76u7uXqJDfmmvtXz5crRu3RoqlQqOjo5o0KABNm/eXOrn+ODno11mTY7ZKJfLMWTIEGzcuFF3Gnv16tVQqVR47rnnyp330qVLcHNzK3GF7oPrsDSDBg1Cp06dMGrUKLi4uGDw4MH4/vvvqxTCGjZsWKWO9A+uX5lMBl9f3xofE/PSpUtwd3cvsa0vKzM8+H0GSu6jP//8c5w6dQqenp7o0KEDPvjggwp/kD7IIFc1VsbgwYOxePFibNmyBf369cP3338Pf39/BAQEVMvyy+oU+eD59ZpU1hWRFf0K1zLUIHhaFdWr/aLOmDGjzEuOq/PS/aFDh2LTpk3YtGkT+vbtix9//BFPPvkkGjRoUKXlaDvNPuxwJg9+DhqNBjKZDFu2bCl1nVVmHZibmyM4OBh79uzBuXPnkJKSgtDQULi4uKCwsBCHDh3C3r174e/vr3u/2vU/ceJEREZGlrpc7XvUaDR44okn8Pbbb5faTrvj0nrU/6vaZbRq1QqtWrVCSEgIunbtitWrV+vCXkUee+wxxMfH4/Lly6VuEIHivnAA0Lx58xLPyWQyNG3aFE2bNkWvXr3g5+eH1atXlzkeFlB8NHPZsmWYMGECQkJCYGdnB5lMhsGDB5e6Y6qO9VQVffv2hVKpxPDhw5Gfn1/hjw6FQlHp9Q0AVlZWuiOtTk5O+PDDD7Fly5ZK9R2qrPuPKtaUVatWYcSIEejXrx/eeustODs76zqsa/u41QZRUVGYMWMGNm7ciOeffx5r1qxB7969YWdnV2OvaWlpiT179mDXrl3YvHkztm7dinXr1qFbt27Ytm1bpa7cr4nPsLx9tKFGE6jM93ngwIEIDQ3Fhg0bsG3bNsyYMQPTp0/HTz/9hJ49e1bqdWpN8AoLC4ObmxvWrVuHzp07Y+fOnXjvvff02nh5eeHEiRPQaDR6R73OnDmje74sDg4OpV7xUdrRosoGHO3rJSQk6P2CLCgowMWLF6u0wavodTQaDRITE/XGuEpNTUVmZma577s8jxrktIdkbW1tq+29lqdv376wsbHBmjVrYG5ujlu3blX5NCMAXcfrssLKgxITE/V+pZ87dw4ajUZ39YyPjw+EEGjcuHGJAPOg8tZ5aGgopk+fjh07dsDJyQn+/v6QyWRo0aIF9u7di71796J379669tr/c+bm5hWufx8fH+Tm5hrkcyqNdjiZ69evV3qe3r17Y+3atVixYkWpA+hmZ2fj559/hr+/f4UhukmTJnBwcKjw9devX4/hw4dj1qxZuml3794t8+q6imi/m4mJiXrbiPT09EodPXyQpaUl+vXrh1WrVqFnz556w6VUt6p8Ztqjr/f//z579iwAVDiSeFW3o9euXSsxDMWDr7V+/Xo0adIEP/30k15NDw7iq6U9en+/s2fPVsso6OV951u2bIk2bdpg9erV8PDwwOXLl/Hll19WuEztBRG5ubl6P+wSEhIqVZNcLkf37t3RvXt3zJ49G1OnTsV7772HXbt2ISIiotp/5D+4foUQOHfuHFq3bq2bVt4++v7/F1WpzcvLCzt27CjRbaEymaE8bm5uGDt2LMaOHYu0tDS0bdsWn376aaWDl+TDSWjJ5XIMGDAAv/76K1auXImioiK904wA8NRTTyElJUWvL1hRURG+/PJL1KtXD+Hh4WUu38fHB1lZWbpfyEDxBuX+vg5a1tbWldrQRkREwMLCAnPnztVLxN988w2ysrLQq1evCpdRGU899RQAlLjqcPbs2QDw0K9jbW39SKcCg4KC4OPjg5kzZ5Y6IviDl9U/KktLS/Tv3x+//fYbFi5cCGtra70xjSpjzZo1+PrrrxESEoLu3btXap758+fr/a3dMGq/ZM888wwUCgU+/PDDEkc6hBB6V8SUt85DQ0ORn5+POXPmoHPnzroNjPYK22vXrun6WgHFl/p36dIFixcvLnXneP/6HzhwIA4ePIjff/+9RLvMzEwUFRWVuw4qa+/evaVe+ajty1GZUyFaAwYMQPPmzfHZZ5/hr7/+0ntOo9FgzJgxuHXrlt7O9NChQ6VeLXj48GHcuHGjwtdXKBQlPsMvv/zyoY+MR0REwNzcHF9++aXech9lJPOJEydiypQpeP/99x96GfeLi4srdXpVPrNr167pbUuzs7OxYsUKBAYGwtXVtdx5q7odLSoqwuLFi3V/FxQUYPHixWjQoAGCgoIA/Hfk4v7lHTp0CAcPHiy1ho0bN+r6QwLF/18OHTpU6R1peSrazg4bNgzbtm3DnDlz4OjoWKnXfOqpp1BUVISFCxfqpqnV6kqFtps3b5aYpj1joe3qow21D/uD40ErVqzQ6zawfv16XL9+Xe+9+vj44M8//9Qb3WDTpk0lhp2oSm1PPfUU1Go15s2bpzf9iy++gEwmq/Lnq1arS3yWzs7OcHd3r9IdUQxyxGvLli26hHm/jh076iXZQYMG4csvv8SUKVPQqlUrvaM7QHFH1sWLF2PEiBE4evQovL29sX79euzfvx9z5swpt8/O4MGD8c4776B///54/fXXcfv2bSxcuBBNmzYt0dkyKCgIO3bs0A0w2Lhx4xLjCAHF/TRiY2Px4YcfokePHujbty8SEhKwYMECtG/fvsIOkpUVEBCA4cOH46uvvkJmZibCw8Nx+PBhLF++HP369avwUt6yBAUFYd26dYiJiUH79u1Rr1499OnTp9Lzy+VyfP311+jZsydatGiBkSNHomHDhrh69Sp27doFW1tb3ZAA1WXo0KFYsWIFfv/9dwwZMqTcwRfXr1+PevXqoaCgQDdy/f79+xEQEIAffvih0q958eJF9O3bFz169MDBgwexatUqvPDCC7rT4D4+Pvjkk08QGxuLpKQk9OvXDzY2Nrh48SI2bNiAl19+WTcWXXnrPCQkBGZmZkhISNAbriIsLEy3gb0/eAHFobBz585o1aoVRo8ejSZNmiA1NRUHDx7ElStXcPz4cQDAW2+9hV9++QW9e/fGiBEjEBQUhLy8PJw8eRLr169HUlJStRw9mT59Oo4ePYpnnnlG92v22LFjWLFiBerXr1+iU3l6enqpdyRo3LgxhgwZgvXr16N79+7o3Lmz3sj1a9aswbFjx/Dmm2+WGFJl9erV6N+/P4KCgmBhYYHTp09j6dKlUKlUuvGGytK7d2+sXLkSdnZ2aN68OQ4ePIgdO3bohl6oqgYNGmDixImYNm0aevfujaeeegp///03tmzZ8tDrOyAgoNq6YADA008/jcaNG6NPnz7w8fFBXl4eduzYgV9//RXt27ev1DahadOmeOmll3DkyBG4uLhg6dKlSE1NxbJlyyqct6rbUXd3d0yfPh1JSUlo2rQp1q1bh/j4eHz11Ve6ATR79+6Nn376Cf3790evXr1w8eJFLFq0CM2bNy/1R6Kvry86d+6MMWPG6H78ODo6lnlqvioq2s6+8MILePvtt7FhwwaMGTOmUoOA9unTB506dcKkSZOQlJSkGzetMj+kP/roI+zZswe9evWCl5cX0tLSsGDBAnh4eOj6j/r4+MDe3h6LFi2CjY0NrK2tERwcXOodaCqjfv36uu9wamoq5syZA19fX4wePVrXZtSoUVi/fj169OiBgQMH4vz581i1apVeZ/eq1tanTx907doV7733HpKSkhAQEIBt27bh559/xoQJE0osuyI5OTnw8PDAgAEDEBAQgHr16mHHjh04cuSI3lHyCj3SNZEVKG84CZRyOahGoxGenp6lXv6plZqaKkaOHCmcnJyEhYWFaNWqVamXleKB4SSEEGLbtm2iZcuWwsLCQjRr1kysWrWq1EuEz5w5I8LCwnSXxWsvcX1wOAmtefPmCX9/f2Fubi5cXFzEmDFjxK1bt/TalHXpfFnDXDyosLBQfPjhh6Jx48bC3NxceHp6itjYWL1hBLTLq+xwErm5ueKFF14Q9vb2AoCuDu0lwD/88INee+2l1w+u77///ls888wzwtHRUSiVSuHl5SUGDhwo4uLiKlWHEBUPJ6FVVFQk3NzcBADx22+/ldpG+5lqHyqVSnh4eIjevXuLpUuXllhnZdEu599//xUDBgwQNjY2wsHBQbz22mvizp07Jdr/+OOPonPnzsLa2lpYW1sLf39/ER0dLRISEnRtylrnWu3btxcAxKFDh3TTrly5IgAIT0/PUus8f/68iIqKEq6ursLc3Fw0bNhQ9O7dW6xfv16vXU5OjoiNjRW+vr7CwsJCODk5iY4dO4qZM2fqLsPXfsYzZswo8TqlfacetH//fhEdHS1atmwp7OzshLm5uWjUqJEYMWKEOH/+vF5b7WXppT26d++ua5eWliZiYmKEr6+vUCqVwt7eXkREROiGkLjfiRMnxFtvvSXatm0r6tevL8zMzISbm5t47rnnxLFjx8qtXYjiIQe025d69eqJyMhIcebMmRKXumu3BUeOHNGbv7ThBNRqtfjwww+Fm5ubsLS0FF26dBGnTp0q8/L5B+HecBLlKW84iYqsXbtWDB48WPj4+AhLS0uhUqlE8+bNxXvvvac3BEBZvLy8RK9evcTvv/8uWrduLZRKpfD39y+x/ShrnWlVZTv6119/iZCQEKFSqYSXl5eYN2+eXjuNRiOmTp0qvLy8hFKpFG3atBGbNm0qsb29///7rFmzhKenp1AqlSI0NFQ3XIzWww4nUdF3XojioZMAiAMHDpS6bkpz48YNMWzYMGFrayvs7OzEsGHDxN9//13hcBJxcXHi6aefFu7u7sLCwkK4u7uL559/vsRQMz///LNo3ry5MDMz01tmecPAlDWcxNq1a0VsbKxwdnYWlpaWolevXuLSpUsl5p81a5Zo2LChUCqVolOnTuKvv/4qsczyaittf5qTkyPeeOMN4e7uLszNzYWfn5+YMWOG3vAuQpT9Pbv/M87PzxdvvfWWCAgIEDY2NsLa2loEBASIBQsWlLo+yiK794JERERV5u3tjZYtW5Z503KqWP/+/XHy5MlS7wxCxqfW9PEiIiIyNdevX8fmzZt1d9Qg41drrmokIiIyFRcvXsT+/fvx9ddfw9zcvNxb05Fx4REvIiIiA9u9ezeGDRuGixcvYvny5RVe/UnGg328iIiIiAyER7yIiIiIDITBi4iIiMhATK5zvUajwbVr12BjY2Pwex8SERHRwxFCICcnB+7u7nq3DaxrTC54Xbt2DZ6enlKXQURERA8hOTkZHh4eUpfx0EwueGlvK5ScnAxbW1uJqyEiIqLKyM7OhqenZ7m3B6wLTC54aU8v2traMngRERHVMXW9m1DdPUlKREREVMcweBEREREZCIMXERERkYEweBEREREZCIMXERERkYEweBEREREZCIMXERERkYEweBEREREZCIMXERERkYEweBEREREZCIMXERERkYEweBEREREZiMndJLum5BepkZ6TL3UZZKQcrCxgreTXlYioruOWvJr8cy0bzyw4IHUZZKRU5nK81LkxXgn3ga3KXOpyiIjoITF4VRMZAKUZz9xS9RMA7hZqMH/Xeaw5dBnjuvlhyOONoDRTSF0aERFVkUwIIaQuwpCys7NhZ2eHrKws2NraSl0OUYWEENj+byqmbz2D8+l5AADP+paY+GQz9GntDrlcJnGFREQ1z1j23wxeRHVEkVqDH45ewRfbzyLtXn/CFu62mNTTH6F+DSSujoioZhnL/pvBi6iOuV1QhGX7k7Doj/PIyS8CAIT6OeGdHv5o2dBO4uqIiGqGsey/GbyI6qibeQWYt/McVv6ZhEJ18df46UB3THyyGTzrW0lcHRFR9TKW/TeDF1Edl3zzNmZuS8DP8dcAAOYKGYY+7oVx3fxQ39pC4uqIiKqHsey/GbyIjMSpq1mYvvUM9iZmAABslGZ4tYsPXuzUGJYWvAKSiOo2Y9l/M3gRGZm9ien4bMsZ/HMtGwDgbKPEG080xXNBHjBTcMgTIqqbjGX/zeBFZIQ0GoFfT1zDjN8TcOXWHQCATwNrvN3DH082d4FMxiEoiKhuMZb9N4MXkRHLL1Jj9Z+X8eXORNy6XQgAaOflgEk9/dHOu77E1RERVZ6x7L8ZvIhMQPbdQizefR7f7LuIu4UaAMATzV3wTo9m8HW2kbg6IqKKGcv+m8GLyISkZt/FnB1nse5IMjQCkMuAQe09MSGiKVxsVVKXR0RUJmPZfzN4EZmgc2k5+HxrArb9mwqAN+EmotrPWPbfDF5EJuyvpJuYtuUMjl66BQBwsDLHa938MJQ34SaiWsZY9t8MXkQmjjfhJqK6wFj23wxeRASg+Cbc649ewRc7ziI1mzfhJqLaxVj23wxeRKTnToEaS/df5E24iahWMZb9N4MXEZWKN+EmotrEWPbfDF5EVK7km7cxa1sCNvIm3EQkIWPZfzN4EVGl8CbcRCQlY9l/M3gRUZXwJtxEJAVj2X8zeBFRlfEm3ERkaMay/2bwIqKHVtpNuIO8HBDLm3ATUTUzlv03gxcRPbLsu4X4avcFfL3vAm/CTUQ1wlj23wxeRFRteBNuIqopxrL/lrwn7Pz58+Ht7Q2VSoXg4GAcPny43PaZmZmIjo6Gm5sblEolmjZtit9++81A1RJReVxsVZj2TGtseyMMTzZ3gUYAaw8nI3zGLsz4/Qyy7xZKXSIRkaQkDV7r1q1DTEwMpkyZgmPHjiEgIACRkZFIS0srtX1BQQGeeOIJJCUlYf369UhISMCSJUvQsGFDA1dOROXxdbbBV1HtsP7VEAR5OeBuoQbzd53HE7N3Iy3nrtTlERFJRtJTjcHBwWjfvj3mzZsHANBoNPD09MS4ceMwadKkEu0XLVqEGTNm4MyZMzA3N3+o1zSWQ5VEdYX2Jtwf/PIPrmXdxZQ+zTGyU2OpyyKiOsZY9t+SHfEqKCjA0aNHERER8V8xcjkiIiJw8ODBUuf55ZdfEBISgujoaLi4uKBly5aYOnUq1Gp1ma+Tn5+P7OxsvQcRGY5MJsOTLVwxopM3AGD32XRpCyIikpBkwSsjIwNqtRouLi56011cXJCSklLqPBcuXMD69euhVqvx22+/4f3338esWbPwySeflPk606ZNg52dne7h6elZre+DiConrGkDAMCfF27gbmHZP5aIiIyZ5J3rq0Kj0cDZ2RlfffUVgoKCMGjQILz33ntYtGhRmfPExsYiKytL90hOTjZgxUSk1czFBi62Stwt1OBI0k2pyyEikoRkwcvJyQkKhQKpqal601NTU+Hq6lrqPG5ubmjatCkUiv/uC/fYY48hJSUFBQUFpc6jVCpha2ur9yAiw5PJZAjzKz7qtYenG4nIREkWvCwsLBAUFIS4uDjdNI1Gg7i4OISEhJQ6T6dOnXDu3DloNBrdtLNnz8LNzQ0WFhY1XjMRPRrt6Ub28yIiUyXpqcaYmBgsWbIEy5cvx+nTpzFmzBjk5eVh5MiRAICoqCjExsbq2o8ZMwY3b97E+PHjcfbsWWzevBlTp05FdHS0VG+BiKqgs68T5DLgbGourmfdkbocIiKDM5PyxQcNGoT09HRMnjwZKSkpCAwMxNatW3Ud7i9fvgy5/L9s6Onpid9//x1vvPEGWrdujYYNG2L8+PF45513pHoLRFQFDtYWaO1hj/jkTOw9m4GB7XmxCxGZFt4yiIgM6ovtZ/G/uET0auWG+UPaSl0OEdURxrL/rlNXNRJR3aft57XvXAaK1JoKWhMRGRcGLyIyqAAPO9hZmiPrTiGOX8mSuhwiIoNi8CIigzJTyNHZ1wkAh5UgItPD4EVEBhfWtDh4cVgJIjI1DF5EZHDafl4nrmQi83bpgx8TERkjBi8iMjg3O0s0dakHjSjuZE9EZCoYvIhIEuHaUewTeLqRiEwHgxcRSUJ7unFPYjpMbDhBIjJhDF5EJIn23vWhMpcjNTsfCak5UpdDRGQQDF5EJAmVuQKPN3EEwGEliMh0MHgRkWTC/O7182LwIiITweBFRJIJb1YcvI5cvIXbBUUSV0NEVPMYvIhIMk2crNHQ3hIFag0OXbgpdTlERDWOwYuIJCOTyXRHvXi6kYhMAYMXEUlK28+LHeyJyBQweBGRpDr6OsJMLsOFjDwk37wtdTlERDWKwYuIJGWrMkfbRg4AeLqRiIwfgxcRSS6sqRMABi8iMn4MXkQkufCmzgCAg+dvoKBII3E1REQ1h8GLiCTXwt0WjtYWyM0vwrHLt6Quh4ioxjB4EZHk5HIZQv2KTzfy6kYiMmYMXkRUK4Q1vTesRCKDFxEZLwYvIqoVQu+N53XqajbSc/IlroaIqGYweBFRrdDARokW7rYAgH3neNSLiIwTgxcR1Rra0427Exi8iMg4MXgRUa0Rfi947U3MgEYjJK6GiKj6MXgRUa3RtpEDrC0UuJFXgH+uZUtdDhFRtWPwIqJaw8JMjo6+94aV4NWNRGSEGLyIqFbR9fPieF5EZIQYvIioVgm/N6zEsUu3kHO3UOJqiIiqF4MXEdUqjRyt0NjJGkUagQPnb0hdDhFRtWLwIqJaJ+ze7YN4upGIjA2DFxHVOuHN7t0+6Gw6hOCwEkRkPBi8iKjWebyJIywUcly5dQcXMvKkLoeIqNoweBFRrWNlYYb2jR0AFB/1IiIyFgxeRFQrhfn9d7qRiMhYMHgRUa2k7ed18MIN3C1US1wNEVH1YPAiolqpmYsNXGyVuFuowV9Jt6Quh4ioWjB4EVGtJJPJEOqnHcU+TeJqiIiqB4MXEdVa4U21/bwyJK6EiKh61IrgNX/+fHh7e0OlUiE4OBiHDx8us+23334LmUym91CpVAaslogMpbOvE2QyICE1B9ez7khdDhHRI5M8eK1btw4xMTGYMmUKjh07hoCAAERGRiItrexTC7a2trh+/brucenSJQNWTESG4mBtgQAPewDAXh71IiIjIHnwmj17NkaPHo2RI0eiefPmWLRoEaysrLB06dIy55HJZHB1ddU9XFxcDFgxERlS2L3TjbsTOawEEdV9kgavgoICHD16FBEREbppcrkcEREROHjwYJnz5ebmwsvLC56ennj66afxzz//lNk2Pz8f2dnZeg8iqju0/bz2JWZAreHtg4iobpM0eGVkZECtVpc4YuXi4oKUlJRS52nWrBmWLl2Kn3/+GatWrYJGo0HHjh1x5cqVUttPmzYNdnZ2uoenp2e1vw8iqjkBHnawVZkh604hjl/JlLocIqJHIvmpxqoKCQlBVFQUAgMDER4ejp9++gkNGjTA4sWLS20fGxuLrKws3SM5OdnAFRPRozBTyNHZzwkAsDuBpxuJqG6TNHg5OTlBoVAgNTVVb3pqaipcXV0rtQxzc3O0adMG586dK/V5pVIJW1tbvQcR1S26YSXYz4uI6jhJg5eFhQWCgoIQFxenm6bRaBAXF4eQkJBKLUOtVuPkyZNwc3OrqTKJSGLaDvbHkzORebtA4mqIiB6e5KcaY2JisGTJEixfvhynT5/GmDFjkJeXh5EjRwIAoqKiEBsbq2v/0UcfYdu2bbhw4QKOHTuGoUOH4tKlSxg1apRUb4GIapibnSWautSDRgD7znFYCSKqu8ykLmDQoEFIT0/H5MmTkZKSgsDAQGzdulXX4f7y5cuQy//Lh7du3cLo0aORkpICBwcHBAUF4cCBA2jevLlUb4GIDCDMrwHOpuZid0I6erd2l7ocIqKHIhNCmNT12dnZ2bCzs0NWVhb7exHVIXsT0zHsm8NwsVXiz9jukMlkUpdERAZkLPtvyU81EhFVRnvv+lCZy5GanY+zqblSl0NE9FAYvIioTlCZK/B4E0cAwO6zZd9SjIioNmPwIqI6I8zv3rASvG8jEdVRDF5EVGdoh5U4fPEmbhcUSVwNEVHVMXgRUZ3h08AaDe0tUaDW4NCFm1KXQ0RUZQxeRFRnyGQy3VGv3Wc5ij0R1T0MXkRUp+huH8TgRUR1EIMXEdUpHX0doZDLcCEjD8k3b0tdDhFRlTB4EVGdYqsyR1AjBwA83UhEdQ+DFxHVOWFNnQDwdCMR1T0MXkRU52g72B84fwOFao3E1RARVR6DFxHVOS3d7VDf2gK5+UU4dumW1OUQEVUagxcR1TlyuQyhfsWnG9nPi4jqEgYvIqqTdMNKJDJ4EVHdweBFRHVS6L37Np66mo2M3HyJqyEiqhwGLyKqkxrYKNHC3RYAsJdHvYiojmDwIqI6K0w3in2GxJUQEVUOgxcR1Vlhfv/dPkijERJXQ0RUMQYvIqqzgrwcYG2hwI28Avx7PVvqcoiIKsTgRUR1loWZHCE+HFaCiOoOBi8iqtPCmxWfbmTwIqK6gMGLiOq08Hv9vI5duoWcu4USV0NEVD4GLyKq0xo5WqGxkzWKNAIHzt+QuhwionIxeBFRnRd27/ZBe3i6kYhqOQYvIqrztON57T6bDiE4rAQR1V4MXkRU5z3exBEWCjmu3LqDixl5UpdDRFQmBi8iqvOslWZo5+0AgFc3ElHtxuBFREYhvOl/o9gTEdVWDF5EZBS0/bwOXriBu4VqiashIiodgxcRGQV/Vxs42yhxt1CDv5JuSV0OEVGpGLyIyCjIZDLdUa89iTzdSES1E4MXERkN3bASCQxeRFQ7MXgRkdEI9XWCTAYkpOYgJeuu1OUQEZXA4EVERsPB2gKtPewB8OpGIqqdGLyIyKhoh5XYzX5eRFQLMXgRkVEJb1p838Z9iRlQa3j7ICKqXRi8iMioBHjYw1Zlhqw7hTh+JVPqcoiI9DB4EZFRMVPI0dmv+KgX+3kRUW3D4EVERifM714/LwYvIqplGLyIyOhox/M6npyJzNsFEldDRPSfWhG85s+fD29vb6hUKgQHB+Pw4cOVmu+7776DTCZDv379arZAIqpT3O0t4edcDxoB7DuXIXU5REQ6kgevdevWISYmBlOmTMGxY8cQEBCAyMhIpKWllTtfUlISJk6ciNDQUANVSkR1iXZYCfbzIqLaRPLgNXv2bIwePRojR45E8+bNsWjRIlhZWWHp0qVlzqNWqzFkyBB8+OGHaNKkiQGrJaK6Qnf7oLPpEILDShBR7SBp8CooKMDRo0cRERGhmyaXyxEREYGDBw+WOd9HH30EZ2dnvPTSSxW+Rn5+PrKzs/UeRGT8OjSuD5W5HKnZ+Tibmit1OUREACQOXhkZGVCr1XBxcdGb7uLigpSUlFLn2bdvH7755hssWbKkUq8xbdo02NnZ6R6enp6PXDcR1X4qcwWCGzsC4OlGIqo9JD/VWBU5OTkYNmwYlixZAicnp0rNExsbi6ysLN0jOTm5hqskotri/tONRES1gZmUL+7k5ASFQoHU1FS96ampqXB1dS3R/vz580hKSkKfPn100zQaDQDAzMwMCQkJ8PHx0ZtHqVRCqVTWQPVEVNuFN22AjwEcvngTtwuKYGUh6SaPiEjaI14WFhYICgpCXFycbppGo0FcXBxCQkJKtPf398fJkycRHx+ve/Tt2xddu3ZFfHw8TyMSkR6fBtZoaG+JArUGhy7clLocIiJpj3gBQExMDIYPH4527dqhQ4cOmDNnDvLy8jBy5EgAQFRUFBo2bIhp06ZBpVKhZcuWevPb29sDQInpREQymQxhTRtg7eHL2H02HV39naUuiYhMnOTBa9CgQUhPT8fkyZORkpKCwMBAbN26Vdfh/vLly5DL61RXNCKqRcKbOmHt4cvsYE9EtYJMmNgAN9nZ2bCzs0NWVhZsbW2lLoeIalj23UK0+Wg71BqBvW93hWd9K6lLIqKHYCz7bx5KIiKjZqsyR9tG9gCAPYk86kVE0mLwIiKjF+Z3b1iJBAYvIpIWgxcRGb3wZsXB68D5GyhUaySuhohMGYMXERm9lu52qG9tgdz8Ihy7dEvqcojIhDF4EZHRk8tlCPUrvtsF+3kRkZQYvIjIJOj6eXFYCSKSEIMXEZmE0KbFR7xOXc1GRm6+xNUQkali8CIik+Bso0Jzt+Kxf/bydCMRSYTBi4hMRljT4tONe85mSFwJEZkqBi8iMhnh94LX3sR0aDQmddMOIqolGLyIyGQEeTnA2kKBjNwC/Hs9W+pyiMgEMXgRkcmwMJMjxKe4kz2vbiQiKTB4EZFJCW/K4EVE0mHwIiKTEt7UGQBw7NIt5NwtlLgaIjI1DF5EZFIaOVrB29EKRRqBA+dvSF0OEZkYBi8iMjn/DSvB041EZFgMXkRkcrTDSuw+mw4hOKwEERkOgxcRmZzHmzjCXCHDlVt3cDEjT+pyiMiEVCl4aTQaLF26FL1790bLli3RqlUr9O3bFytWrOCvRiKqM6yVZmjvXR8ATzcSkWFVOngJIdC3b1+MGjUKV69eRatWrdCiRQtcunQJI0aMQP/+/WuyTiKiahV23+lGIiJDMatsw2+//RZ79uxBXFwcunbtqvfczp070a9fP6xYsQJRUVHVXiQRUXULb9oAn205gz8v3MTdQjVU5gqpSyIiE1DpI15r167Fu+++WyJ0AUC3bt0wadIkrF69ulqLIyKqKf6uNnC2UeJOoRp/Jd2SuhwiMhGVDl4nTpxAjx49yny+Z8+eOH78eLUURURU02Qy2X/DSiTydCMRGUalg9fNmzfh4uJS5vMuLi64dYu/Gomo7uB4XkRkaJUOXmq1GmZmZXcJUygUKCoqqpaiiIgMIdTXCTIZcCYlBylZd6Uuh4hMQKU71wshMGLECCiVylKfz8/Pr7aiiIgMwcHaAq097HE8ORN7EtMxsJ2n1CURkZGrdPCKioqCTCarsA0RUV0S7ueE48mZ2H2WwYuIal6VhpMgIjI24c0aYO7Oc9iXmAG1RkAhL/8HJhHRo6h0H68BAwZg69atHKGeiIxKgIc9bFRmyLpTiONXMqUuh4iMXKWD161bt9CrVy80atQIkydPxoULF2qyLiIigzBTyBHq5wSAVzcSUc2rdPCKi4vDhQsX8NJLL2HVqlXw8/NDt27dsGbNGnasJ6I6LcyPw0oQkWFU6SbZXl5e+OCDD3DhwgVs374d7u7uGD16NNzc3BAdHY2jR4/WVJ1ERDVGO55XfHImsm4XSlwNERmzKgWv+3Xr1g2rVq1CSkoKpk2bhu+++w7BwcHVWRsRkUG421vCz7keNALYdy5D6nKIyIg9dPACgIsXL2LmzJmYOnUqsrKyEBERUV11EREZlPao1+6zaRJXQkTGrMrB6+7du1i1ahW6desGPz8/rFixAi+99BIuXryIrVu31kSNREQ1Llx3+6AMXr1NRDWm0uN4HT58GEuXLsW6detw9+5d9O/fH1u3bkX37t0rHFiViKi269C4PpRmcqRk38XZ1Fw0c7WRuiQiMkKVDl6PP/44AgIC8PHHH2PIkCFwcHCoybqIiAxKZa7A400csftsOvacTWfwIqIaUelTjb1798b+/fvx2muvMXQRkVHS9vPak8hhJYioZlQ6eG3evBm5ubk1WQsRkaTCmxYPpHro4k3cKVBLXA0RGaNKBy92NiUiY+fToB4a2luioEiDPy/ekLocIjJCVbqqsaY60c+fPx/e3t5QqVQIDg7G4cOHy2z7008/oV27drC3t4e1tTUCAwOxcuXKGqmLiEyLTCZD2L2jXrsTeLqRiKpfpTvXA0DTpk0rDF83b96sUgHr1q1DTEwMFi1ahODgYMyZMweRkZFISEiAs7Nzifb169fHe++9B39/f1hYWGDTpk0YOXIknJ2dERkZWaXXJiJ6UHjTBlh7OJn9vIioRshEJc8hyuVyzJkzB3Z2duW2Gz58eJUKCA4ORvv27TFv3jwAgEajgaenJ8aNG4dJkyZVahlt27ZFr1698PHHH1fYNjs7G3Z2dsjKyoKtrW2VaiUi45d9txBtPtoOtUZg79td4VnfSuqSiAjGs/+u0hGvwYMHl3oU6mEVFBTg6NGjiI2N1U2Ty+WIiIjAwYMHK5xfCIGdO3ciISEB06dPr7a6iMh02arM0baRPY4k3cKexHQMCfaSuiQiMiKV7uNVE/27MjIyoFar4eLiojfdxcUFKSkpZc6XlZWFevXqwcLCAr169cKXX36JJ554otS2+fn5yM7O1nsQEZUnzO/e7YPYz4uIqlmdvKrRxsYG8fHxOHLkCD799FPExMTgjz/+KLXttGnTYGdnp3t4enoatlgiqnO043kdOH8DhWqNxNUQkTGpdPDSaDTVepoRAJycnKBQKJCamqo3PTU1Fa6urmXOJ5fL4evri8DAQLz55psYMGAApk2bVmrb2NhYZGVl6R7JycnV+h6IyPi0amiH+tYWyM0vwt+XM6Uuh4iMSJVvkl2dLCwsEBQUhLi4ON00jUaDuLg4hISEVHo5Go0G+fn5pT6nVCpha2ur9yAiKo9cLkNn33vDSpxNk7gaIjImkgYvAIiJicGSJUuwfPlynD59GmPGjEFeXh5GjhwJAIiKitLrfD9t2jRs374dFy5cwOnTpzFr1iysXLkSQ4cOleotEJERCtfePuhshsSVEJExqdJVjTVh0KBBSE9Px+TJk5GSkoLAwEBs3bpV1+H+8uXLkMv/y4d5eXkYO3Ysrly5AktLS/j7+2PVqlUYNGiQVG+BiIxQ6L2BVE9ezUJGbj6c6iklroiIjEGlx/EyFsYyDggR1byn/rcX/17PxpxBgejXpqHU5RCZNGPZf0t+qpGIqLbSXt24+yyHlSCi6sHgRURUBu19G/cmpkOjMamTA0RUQxi8iIjK0M6rPqwsFMjILcC/1zn4MhE9OgYvIqIyWJjJ0dHHEQBPNxJR9WDwIiIqx3/DSjB4EdGjY/AiIiqHtoP90Uu3kHO3UOJqiKiuY/AiIiqHl6M1vB2tUKQROHj+htTlEFEdx+BFRFQBDitBRNWFwYuIqAJhfvf6eSWmw8TGnCaiasbgRURUgRAfR5grZEi+eQdJN25LXQ4R1WEMXkREFbBWmqGdV30AwO6ENImrIaK6jMGLiKgSwptpTzdmSFwJEdVlDF5ERJWg7ed18PwN5BepJa6GiOoqBi8iokp4zM0GDWyUuFOoxl9Jt6Quh4jqKAYvIqJKkMlkuqNeHFaCiB4WgxcRUSWFNXUCwNsHEdHDY/AiIqqkUL8GkMmAMyk5SM2+K3U5RFQHMXgREVVSfWsLtPawBwBsOnFd2mKIqE5i8CIiqoLngjwAAN/svYCCIo3E1RBRXcPgRURUBQOCPNDARolrWXexMf6q1OUQUR3D4EVEVAUqcwVGhzYGACz64zzUGt67kYgqj8GLiKiKXgj2gp2lOS5k5GHrqRSpyyGiOoTBi4ioiuopzTCiozcAYP6ucxCCR72IqHIYvIiIHsKIjt6wslDg3+vZ+IPjehFRJTF4ERE9BAdrCwwJbgQAWLDrnMTVEFFdweBFRPSQRoU2gYVCjiNJt3D44k2pyyGiOoDBi4joIbnYqjCgXfG4XvN51IuIKoHBi4joEbwa5gO5rPjG2aeuZkldDhHVcgxeRESPoJGjFfoGuAMAFvzBo15EVD4GLyKiRzSmiy8AYMupFJxLy5W4GiKqzRi8iIgeUTNXGzzR3AVCAIt2n5e6HCKqxRi8iIiqwdguPgCAjX9fxZVbtyWuhohqKwYvIqJq0KaRAzr5OqJII7BkzwWpyyGiWorBi4iomkTf6+v13ZFkpOfkS1wNEdVGDF5ERNUkxMcRgZ72yC/SYOn+i1KXQ0S1EIMXEVE1kclkiO5afNRr5cFLyLpTKHFFRFTbMHgREVWj7v7OaOZig9z8Iqw8mCR1OURUyzB4ERFVI7lchjH3rnBcuj8JdwrUEldERLUJgxcRUTXr3doNnvUtcTOvAN8duSx1OURUizB4ERFVMzOFHK+GFx/1+mrPBRQUaSSuiIhqi1oRvObPnw9vb2+oVCoEBwfj8OHDZbZdsmQJQkND4eDgAAcHB0RERJTbnohICs+29YCzjRLXs+5i499XpS6HiGoJyYPXunXrEBMTgylTpuDYsWMICAhAZGQk0tLSSm3/xx9/4Pnnn8euXbtw8OBBeHp64sknn8TVq9ywEVHtoTJXYHRoEwDAwt3nodYIiSsiotpAJoSQdGsQHByM9u3bY968eQAAjUYDT09PjBs3DpMmTapwfrVaDQcHB8ybNw9RUVEVts/OzoadnR2ysrJga2v7yPUTEZUlL78IHT/biaw7hZj3Qhv0bu0udUlEdZax7L8lPeJVUFCAo0ePIiIiQjdNLpcjIiICBw8erNQybt++jcLCQtSvX7+myiQieijWSjOM7OQNAJi/6zwk/p1LRLWApMErIyMDarUaLi4uetNdXFyQkpJSqWW88847cHd31wtv98vPz0d2drbeg4jIUEZ09IaVhQKnr2fjj4R0qcshIolJ3sfrUXz22Wf47rvvsGHDBqhUqlLbTJs2DXZ2drqHp6engaskIlNmb2WBoY97AQDm7TrHo15EJk7S4OXk5ASFQoHU1FS96ampqXB1dS133pkzZ+Kzzz7Dtm3b0Lp16zLbxcbGIisrS/dITk6ultqJiCprVOfGsFDIcfTSLRy+eFPqcohIQpIGLwsLCwQFBSEuLk43TaPRIC4uDiEhIWXO9/nnn+Pjjz/G1q1b0a5du3JfQ6lUwtbWVu9BRGRIzrYqPNfOAwAw/4/zEldDRFKS/FRjTEwMlixZguXLl+P06dMYM2YM8vLyMHLkSABAVFQUYmNjde2nT5+O999/H0uXLoW3tzdSUlKQkpKC3Nxcqd4CEVGFXgnzgUIuw56z6Th5JUvqcohIIpIHr0GDBmHmzJmYPHkyAgMDER8fj61bt+o63F++fBnXr1/XtV+4cCEKCgowYMAAuLm56R4zZ86U6i0QEVWokaMV+gYUDyex4I9zEldDRFKRfBwvQzOWcUCIqO45m5qDJ7/YA5kM2P5GGHydbaQuiajOMJb9t+RHvIiITEVTFxs82dwFQgAL/7ggdTlEJAEGLyIiAxrb1RcAsDH+KpJv3pa4GiIyNAYvIiIDCvS0R2dfJ6g1Akv28qgXkalh8CIiMrCxXX0AAN8dSUZazl2JqyEiQ2LwIiIysJAmjmjTyB4FRRos3ZckdTlEZEAMXkREBiaTyRDdpbiv16o/LyHrdqHEFRGRoTB4ERFJoJu/M/xdbZCbX4QVB5OkLoeIDITBi4hIAnK5DGO6FPf1Wrr/Im4XFElcEREZAoMXEZFEerVyg5ejFW7dLsR3h5OlLoeIDIDBi4hIImYKOV4NLz7q9dWeCygo0khcERHVNAYvIiIJPdO2IVxslUjJvosNf1+RuhwiqmEMXkREElKaKTA6tAkAYOEf56HWmNTtc4lMDoMXEZHEnu/QCPZW5ki6cRu/nbwudTlEVIMYvIiIJGatNMPIjo0BAPN3nYMQPOpFZKwYvIiIaoHhHb1gbaHAmZQc7EpIk7ocIqohDF5ERLWAvZUFhj7uBQCYt5NHvYiMFYMXEVEt8VLnxrAwk+PY5UwcunhT6nKIqAYweBER1RLOtioMbOcBoLivFxEZHwYvIqJa5JUwHyjkMuxNzMCJK5lSl0NE1YzBi4ioFvGsb4WnA9wBAAt2nZe4GiKqbgxeRES1jPbm2Vv/SUFiao7E1RBRdWLwIiKqZfxcbBDZwgUAsHA3j3oRGRMGLyKiWmhsF18AwM/x15B887bE1RBRdWHwIiKqhQI87RHq5wS1RuCrPRekLoeIqgmDFxFRLaU96rXur2Sk5dyVuBoiqg4MXkREtdTjTeqjbSN7FBRp8M2+i1KXQ0TVgMGLiKiWkslkiO5afNRr1cFLyLpdKHFFRPSoGLyIiGqxbv7O8He1QV6BGssPJkldDhE9IgYvIqJaTCaTYey9o17L9l/E7YIiiSsiokfB4EVEVMv1auUGb0cr3LpdiLWHk6Uuh4geAYMXEVEtp5DL8Gp48Wj2S/ZcQH6RWuKKiOhhMXgREdUB/ds2hKutCinZd7Hh2FWpyyGih8TgRURUByjNFBgd1gRA8W2EitQaiSsioofB4EVEVEc838ETDlbmuHTjNn47lSJ1OUT0EBi8iIjqCCsLM4zs1BgAsGDXOQghJK6IiKqKwYuIqA4ZHuINawsFzqTkYOeZNKnLIaIqYvAiIqpD7KzMMTTECwAwj0e9iOocBi8iojrmpc6NYWEmx9+XM/HnhZtSl0NEVcDgRURUxzjbqDConScAYMEf5ySuhoiqgsGLiKgOejmsCRRyGfYmZuB4cqbU5RBRJUkevObPnw9vb2+oVCoEBwfj8OHDZbb9559/8Oyzz8Lb2xsymQxz5swxXKFERLWIZ30rPB3oDoBHvYjqEkmD17p16xATE4MpU6bg2LFjCAgIQGRkJNLSSr9S5/bt22jSpAk+++wzuLq6GrhaIqLaZWwXH8hkwO//pCIxNUfqcoioEiQNXrNnz8bo0aMxcuRING/eHIsWLYKVlRWWLl1aavv27dtjxowZGDx4MJRKpYGrJSKqXXydbRDZvPhH6MI/zktcDRFVhmTBq6CgAEePHkVERMR/xcjliIiIwMGDB6Uqi4ioThnbtfjm2T8fv4bkm7clroaIKiJZ8MrIyIBarYaLi4vedBcXF6SkVN+tMPLz85Gdna33ICIyFq097BHq5wS1RmDxHh71IqrtJO9cX9OmTZsGOzs73cPT01PqkoiIqlV0V18AwPd/XUFa9l2JqyGi8kgWvJycnKBQKJCamqo3PTU1tVo7zsfGxiIrK0v3SE5OrrZlExHVBsGN6yPIywEFRRp8s++i1OUQUTkkC14WFhYICgpCXFycbppGo0FcXBxCQkKq7XWUSiVsbW31HkRExkQmkyH6Xl+vVX9eQubtAokrIqKySHqqMSYmBkuWLMHy5ctx+vRpjBkzBnl5eRg5ciQAICoqCrGxsbr2BQUFiI+PR3x8PAoKCnD16lXEx8fj3DmOYUNEpq1rM2f4u9ogr0CN5QcuSV0OEZVB0uA1aNAgzJw5E5MnT0ZgYCDi4+OxdetWXYf7y5cv4/r167r2165dQ5s2bdCmTRtcv34dM2fORJs2bTBq1Cip3gIRUa1QfNSruK/XsgMXkZdfJHFFRFQamTCxW9tnZ2fDzs4OWVlZPO1IREZFrRHoPusPJN24jf/r9RhGhTaRuiSiamMs+2+jv6qRiMhUKOQyjOlS3Ndryd4LyC9SS1wRET2IwYuIyIj0b+MBNzsVUrPz8dOxq1KXQ0QPYPAiIjIiFmZyjL53inHR7vMoUmskroiI7sfgRURkZAZ38ER9awtcunEbm09er3gGIjIYBi8iIiNjZWGGFzt5AwAW7DoPjcakrqEiqtUYvIiIjNCwEG/UU5ohITUHO8+kSV0OEd3D4EVEZITsLM0x9HEvAMC8XedgYiMHEdVaDF5EREbqpc6NoTSTIz45Ewcv3JC6HCICgxcRkdFqYKPEoPaeAIr7ehGR9Bi8iIiM2MthTWAml2HfuQzEJ2dKXQ6RyWPwIiIyYh4OVng6sCEAYMGucxJXQ0QMXkRERm5MlyaQyYBt/6bibGqO1OUQmTQGLyIiI+frbIMeLVwBAAv/YF8vIikxeBERmYCxXXwBAL8cv4bLN25LXA2R6WLwIiIyAa087BDWtAHUGoHFe3jUi0gqDF5ERCYiuosPAOCHv64gLfuuxNUQmSYGLyIiE9GhcX2083JAgVqDr/ddlLocIpPE4EVEZCJkMhmiuxb39Vr15yVk3i6QuCIi08PgRURkQro0a4DH3Gxxu0CNbw8kSV0Okclh8CIiMiHFR72K+3ot25+E3PwiiSsiMi0MXkREJqZnSzc0drJG1p1C9PlyH74/koyCIo3UZRGZBAYvIiITo5DL8Em/lnCwMsfFjDy8/eMJdJ35B1YcTMLdQrXU5REZNZkQQkhdhCFlZ2fDzs4OWVlZsLW1lbocIiLJ5OUXYc2hy/hq7wWk5+QDABrYKDE6tDGGBHvBWmkmcYVE/zGW/TeDFxGRibtbqMYPfyVj0e4LuJp5BwBgb2WOkR0bY0RHb9hZmUtcIZHx7L8ZvIiICABQqNZgw99XsfCP87iYkQcAqKc0w7AQL7zUuTGc6iklrpBMmbHsvxm8iIhIj1oj8NvJ65i/6xzOpOQAAFTmcjzfoRFeDmsCNztLiSskU2Qs+28GLyIiKpVGIxB3Jg3zdp3D8eRMAIC5QoYBQZ4YE+6DRo5W0hZIJsVY9t8MXkREVC4hBPady8C8nedw6OJNAMVXRvYNcMfYLj7wc7GRuEIyBcay/2bwIiKiSjuSdBPzdp7D7rPpAACZDOjRwhXRXX3RsqGdxNWRMTOW/TeDFxERVdnJK1mYv+sctv6TopvWtVkDvNbNF0Fe9SWsjIyVsey/GbyIiOihnU3NwYJd5/DL8WvQ3NubPN6kPsZ180NHH0fIZDJpCySjYSz7bwYvIiJ6ZEkZeVi0+zx+PHYFheri3Uqgpz3GdfNFN39nBjB6ZMay/2bwIiKianMt8w6+2nMBaw9fRv69+z8+5maL6K4+6NnSDQo5Axg9HGPZfzN4ERFRtUvPycfX+y5g1cFLyCsovv9jkwbWGNvFF08HusNcwVsFU9UYy/6bwYuIiGpM5u0CfHsgCcv2JyHrTiEAwMPBEq+G+2BAkAdU5gqJK6S6wlj23wxeRERU43Lzi7Dqz0v4eu8FZOQWAACcbZR4OawJXghuBCsL3pCbymcs+28GLyIiMpi7hWp8d/gyFu+5gOtZdwEA9a0t8GInbwwL8YadJW/ITaUzlv03gxcRERlcQZEGPx27goW7z+PSjdsAABulGYZ39MaLnRujvrWFxBVSbWMs+28GLyIikkyRWoPNJ69j3s5zSEzLBQBYmivwQnDxDbldbFUSV0i1hbHsvxm8iIhIchqNwLZ/UzF/1zmcvJoFALBQyPFcOw+8Gu4Dz/q8IbepM5b9N4MXERHVGkII7D6bjvm7zuFI0i0AxTfk7hfYEGO7+sCnQT2JKySpGMv+u1YMpDJ//nx4e3tDpVIhODgYhw8fLrf9Dz/8AH9/f6hUKrRq1Qq//fabgSolIqKaJJPJ0KWZM354tSPWvfw4Qv2coNYI/HjsCiJm70b06mP491q21GUSPTTJg9e6desQExODKVOm4NixYwgICEBkZCTS0tJKbX/gwAE8//zzeOmll/D333+jX79+6NevH06dOmXgyomIqCYFN3HEypeCsTG6EyIec4EQwOaT1/HU3L146dsjOHb5ltQlElWZ5Kcag4OD0b59e8ybNw8AoNFo4OnpiXHjxmHSpEkl2g8aNAh5eXnYtGmTbtrjjz+OwMBALFq0qMLXM5ZDlUREpuZMSjbm7zqPTSeuQbvn6uTriDHhvvB2Yh+w2sTCTA5nm+q9MMJY9t+SjlhXUFCAo0ePIjY2VjdNLpcjIiICBw8eLHWegwcPIiYmRm9aZGQkNm7cWGr7/Px85Ofn6/7OzuYhaiKiusjf1RZfPt8Gb0T4YeEf57Hh76vYf+4G9p+7IXVp9IC2jezx09hOUpdRK0kavDIyMqBWq+Hi4qI33cXFBWfOnCl1npSUlFLbp6SklNp+2rRp+PDDD6unYCIiklyTBvUw47kAjI/ww+LdF/Bz/FXdDbmpduC9OMtm9PdoiI2N1TtClp2dDU9PTwkrIiKi6uDhYIWP+7XEx/1aSl0KUaVJGrycnJygUCiQmpqqNz01NRWurq6lzuPq6lql9kqlEkqlsnoKJiIiInoEkh4LtLCwQFBQEOLi4nTTNBoN4uLiEBISUuo8ISEheu0BYPv27WW2JyIiIqotJD/VGBMTg+HDh6Ndu3bo0KED5syZg7y8PIwcORIAEBUVhYYNG2LatGkAgPHjxyM8PByzZs1Cr1698N133+Gvv/7CV199JeXbICIiIqqQ5MFr0KBBSE9Px+TJk5GSkoLAwEBs3bpV14H+8uXLkMv/OzDXsWNHrFmzBv/3f/+Hd999F35+fti4cSNatuQ5fiIiIqrdJB/Hy9CMZRwQIiIiU2Is+29e70lERERkIAxeRERERAbC4EVERERkIAxeRERERAbC4EVERERkIAxeRERERAbC4EVERERkIAxeRERERAbC4EVERERkIJLfMsjQtAP1Z2dnS1wJERERVZZ2v13Xb7hjcsErJycHAODp6SlxJURERFRVOTk5sLOzk7qMh2Zy92rUaDS4du0abGxsIJPJqnXZ2dnZ8PT0RHJycp2+j5Sx4OdRu/DzqF34edQ+/EzKJ4RATk4O3N3dIZfX3Z5SJnfESy6Xw8PDo0Zfw9bWll+aWoSfR+3Cz6N24edR+/AzKVtdPtKlVXcjIxEREVEdw+BFREREZCAMXtVIqVRiypQpUCqVUpdC4OdR2/DzqF34edQ+/ExMg8l1riciIiKSCo94ERERERkIgxcRERGRgTB4ERERERkIgxcRERGRgTB4VZP58+fD29sbKpUKwcHBOHz4sNQlmaxp06ahffv2sLGxgbOzM/r164eEhASpy6J7PvvsM8hkMkyYMEHqUkzW1atXMXToUDg6OsLS0hKtWrXCX3/9JXVZJkmtVuP9999H48aNYWlpCR8fH3z88cd1/n6EVDYGr2qwbt06xMTEYMqUKTh27BgCAgIQGRmJtLQ0qUszSbt370Z0dDT+/PNPbN++HYWFhXjyySeRl5cndWkm78iRI1i8eDFat24tdSkm69atW+jUqRPMzc2xZcsW/Pvvv5g1axYcHBykLs0kTZ8+HQsXLsS8efNw+vRpTJ8+HZ9//jm+/PJLqUujGsLhJKpBcHAw2rdvj3nz5gEovh+kp6cnxo0bh0mTJklcHaWnp8PZ2Rm7d+9GWFiY1OWYrNzcXLRt2xYLFizAJ598gsDAQMyZM0fqskzOpEmTsH//fuzdu1fqUghA79694eLigm+++UY37dlnn4WlpSVWrVolYWVUU3jE6xEVFBTg6NGjiIiI0E2Ty+WIiIjAwYMHJayMtLKysgAA9evXl7gS0xYdHY1evXrpfVfI8H755Re0a9cOzz33HJydndGmTRssWbJE6rJMVseOHREXF4ezZ88CAI4fP459+/ahZ8+eEldGNcXkbpJd3TIyMqBWq+Hi4qI33cXFBWfOnJGoKtLSaDSYMGECOnXqhJYtW0pdjsn67rvvcOzYMRw5ckTqUkzehQsXsHDhQsTExODdd9/FkSNH8Prrr8PCwgLDhw+XujyTM2nSJGRnZ8Pf3x8KhQJqtRqffvophgwZInVpVEMYvMioRUdH49SpU9i3b5/UpZis5ORkjB8/Htu3b4dKpZK6HJOn0WjQrl07TJ06FQDQpk0bnDp1CosWLWLwksD333+P1atXY82aNWjRogXi4+MxYcIEuLu78/MwUgxej8jJyQkKhQKpqal601NTU+Hq6ipRVQQAr732GjZt2oQ9e/bAw8ND6nJM1tGjR5GWloa2bdvqpqnVauzZswfz5s1Dfn4+FAqFhBWaFjc3NzRv3lxv2mOPPYYff/xRoopM21tvvYVJkyZh8ODBAIBWrVrh0qVLmDZtGoOXkWIfr0dkYWGBoKAgxMXF6aZpNBrExcUhJCREwspMlxACr732GjZs2ICdO3eicePGUpdk0rp3746TJ08iPj5e92jXrh2GDBmC+Ph4hi4D69SpU4nhVc6ePQsvLy+JKjJtt2/fhlyuvytWKBTQaDQSVUQ1jUe8qkFMTAyGDx+Odu3aoUOHDpgzZw7y8vIwcuRIqUszSdHR0VizZg1+/vln2NjYICUlBQBgZ2cHS0tLiaszPTY2NiX611lbW8PR0ZH97iTwxhtvoGPHjpg6dSoGDhyIw4cP46uvvsJXX30ldWkmqU+fPvj000/RqFEjtGjRAn///Tdmz56NF198UerSqIZwOIlqMm/ePMyYMQMpKSkIDAzE3LlzERwcLHVZJkkmk5U6fdmyZRgxYoRhi6FSdenShcNJSGjTpk2IjY1FYmIiGjdujJiYGIwePVrqskxSTk4O3n//fWzYsAFpaWlwd3fH888/j8mTJ8PCwkLq8qgGMHgRERERGQj7eBEREREZCIMXERERkYEweBEREREZCIMXERERkYEweBEREREZCIMXERERkYEweBEREREZCIMXERmct7d3lQdPHTFiBPr166f7u0uXLpgwYUK11lUTZDIZNm7cKHUZRFRLMHgRmaARI0ZAJpPpHo6OjujRowdOnDghdWmV9tNPP+Hjjz+WuowKXb9+HT179pS6DCKqJRi8iExUjx49cP36dVy/fh1xcXEwMzND7969pS6r0urXrw8bGxupy6iQq6srlEql1GUQUS3B4EVkopRKJVxdXeHq6orAwEBMmjQJycnJSE9P17U5efIkunXrBktLSzg6OuLll19Gbm6u7nnt6b+ZM2fCzc0Njo6OiI6ORmFhoa5NWloa+vTpA0tLSzRu3BirV6+usDa1Wo2YmBjY29vD0dERb7/9Nh68u9mDpxq9vb3xySefICoqCvXq1YOXlxd++eUXpKen4+mnn0a9evXQunVr/PXXX3rL2bdvH0JDQ2FpaQlPT0+8/vrryMvL01vu1KlT8eKLL8LGxgaNGjXSu6F0QUEBXnvtNbi5uUGlUsHLywvTpk3TPf/gqcbqWKdEVHcxeBERcnNzsWrVKvj6+sLR0REAkJeXh8jISDg4OODIkSP44YcfsGPHDrz22mt68+7atQvnz5/Hrl27sHz5cnz77bf49ttvdc+PGDECycnJ2LVrF9avX48FCxYgLS2t3HpmzZqFb7/9FkuXLsW+fftw8+ZNbNiwocL38cUXX6BTp074+++/0atXLwwbNgxRUVEYOnQojh07Bh8fH0RFRelC3Pnz59GjRw88++yzOHHiBNatW4d9+/aVeI+zZs1Cu3bt8Pfff2Ps2LEYM2YMEhISAABz587FL7/8gu+//x4JCQlYvXo1vL29S62vutYpEdVhgohMzvDhw4VCoRDW1tbC2tpaABBubm7i6NGjujZfffWVcHBwELm5ubppmzdvFnK5XKSkpOiW4+XlJYqKinRtnnvuOTFo0CAhhBAJCQkCgDh8+LDu+dOnTwsA4osvviizPjc3N/H555/r/i4sLBQeHh7i6aef1k0LDw8X48eP1/3t5eUlhg4dqvv7+vXrAoB4//33ddMOHjwoAIjr168LIYR46aWXxMsvv6z32nv37hVyuVzcuXOn1OVqNBrh7OwsFi5cKIQQYty4caJbt25Co9GU+l4AiA0bNgghqmedElHdxiNeRCaqa9euiI+PR3x8PA4fPozIyEj07NkTly5dAgCcPn0aAQEBsLa21s3TqVMnaDQa3dEeAGjRogUUCoXubzc3N90RrdOnT8PMzAxBQUG65/39/WFvb19mXVlZWbh+/TqCg4N108zMzNCuXbsK31Pr1q11/3ZxcQEAtGrVqsQ0bX3Hjx/Ht99+i3r16ukekZGR0Gg0uHjxYqnLlclkcHV11S1jxIgRiI+PR7NmzfD6669j27ZtZdZXHeuUiOo2M6kLICJpWFtbw9fXV/f3119/DTs7OyxZsgSffPJJpZdjbm6u97dMJoNGo6m2Oqvi/lpkMlmZ07T15ebm4pVXXsHrr79eYlmNGjUqdbna5WiX0bZtW1y8eBFbtmzBjh07MHDgQERERGD9+vXV8j4efD0iqtt4xIuIABTv3OVyOe7cuQMAeOyxx3D8+HG9jub79++HXC5Hs2bNKrVMf39/FBUV4ejRo7ppCQkJyMzMLHMeOzs7uLm54dChQ7ppDy6jurRt2xb//vsvfH19SzwsLCwqvRxbW1sMGjQIS5Yswbp16/Djjz/i5s2bJdpVxzolorqNwYvIROXn5yMlJQUpKSk4ffo0xo0bh9zcXPTp0wcAMGTIEKhUKgwfPhynTp3Crl27MG7cOAwbNkx3yq4izZo1Q48ePfDKK6/g0KFDOHr0KEaNGgVLS8ty5xs/fjw+++wzbNy4EWfOnMHYsWPLDWsP65133sGBAwfw2muvIT4+HomJifj5559LdHYvz+zZs7F27VqcOXMGZ8+exQ8//ABXV9dST6dWxzolorqNwYvIRG3duhVubm5wc3NDcHCw7iq7Ll26AACsrKzw+++/4+bNm2jfvj0GDBiA7t27Y968eVV6nWXLlsHd3R3h4eF45pln8PLLL8PZ2bnced58800MGzYMw4cPR0hICGxsbNC/f/+Hfatlat26NXbv3o2zZ88iNDQUbdq0weTJk+Hu7l7pZdjY2ODzzz9Hu3bt0L59eyQlJeG3336DXF5y81pd65SI6i6ZEA8MjkNERERENYJHvIiIiIgMhMGLiIiIyEAYvIiIiIgMhMGLiIiIyEAYvIiIiIgMhMGLiIiIyEAYvIiIiIgMhMGLiIiIyEAYvIiIiIgMhMGLiIiIyEAYvIiIiIgMhMGLiIiIyED+H2m2TPPYS6DzAAAAAElFTkSuQmCC", "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, until 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 }