From 56781871dd518f836fd543954834e82f18c22208 Mon Sep 17 00:00:00 2001 From: Gabriel Matos Date: Wed, 10 Jan 2024 15:43:11 +0000 Subject: [PATCH] Implement experimental general `get_params_to_statevector` * Supports general operations beyond parameterized gate application * Allows for specification of parameters using dictionaries * Adds a classical register for storage of intermediate computations (e.g. mid-circuit-measurements) * Allows the user to define new gates identified by a string In addition: * Introduce unstable experimental API in experimental.py * Add `./examples/experimental` folder for notebooks using unstable API * Add `experimental/noise_channel_monte_carlo.ipynb` with example of simulating noise using Monte-Carlo/quantum trajectories * Add test_experimental.py for testing unstable API * Add `ParameterizedGateFunction`, `UnparameterizedGateFunction` types --- docs/examples.md | 4 + docs/experimental.rst | 8 + docs/experimental/statetensor.rst | 13 + .../statetensor/get_default_gates.rst | 4 + .../statetensor/get_default_operations.rst | 4 + docs/experimental/statetensor/get_params.rst | 4 + .../get_params_to_statetensor_func.rst | 4 + docs/experimental/statetensor/parse_op.rst | 4 + .../statetensor/wrap_parameterised_tensor.rst | 4 + docs/index.rst | 1 + docs/statetensor/apply_gate.rst | 2 +- .../noise_channel_monte_carlo.ipynb | 628 ++++++++++++++++++ qujax/__init__.py | 2 + qujax/densitytensor.py | 4 +- qujax/experimental/statetensor.py | 392 +++++++++++ qujax/statetensor.py | 11 +- qujax/typing.py | 8 +- tests/test_experimental.py | 96 +++ 18 files changed, 1184 insertions(+), 9 deletions(-) create mode 100644 docs/experimental.rst create mode 100644 docs/experimental/statetensor.rst create mode 100644 docs/experimental/statetensor/get_default_gates.rst create mode 100644 docs/experimental/statetensor/get_default_operations.rst create mode 100644 docs/experimental/statetensor/get_params.rst create mode 100644 docs/experimental/statetensor/get_params_to_statetensor_func.rst create mode 100644 docs/experimental/statetensor/parse_op.rst create mode 100644 docs/experimental/statetensor/wrap_parameterised_tensor.rst create mode 100644 examples/experimental/noise_channel_monte_carlo.ipynb create mode 100644 qujax/experimental/statetensor.py create mode 100644 tests/test_experimental.py diff --git a/docs/examples.md b/docs/examples.md index ea7ac98..677a8fe 100644 --- a/docs/examples.md +++ b/docs/examples.md @@ -12,6 +12,10 @@ Below are some use-case notebooks. These both illustrate the flexibility of quja - [classification.ipynb](https://github.com/CQCL/qujax/blob/develop/examples/classification.ipynb) - train a quantum circuit for binary classification using data re-uploading. - [generative_modelling.ipynb](https://github.com/CQCL/qujax/blob/develop/examples/generative_modelling.ipynb) - uses a parameterised quantum circuit as a generative model for a real life dataset. Trains via stochastic gradient Langevin dynamics on the maximum mean discrepancy between statetensor and dataset. +Experimental (i.e. uses an unstable API which might change in future versions): + +- [noise_channel_monte_carlo.ipynb](https://github.com/CQCL/qujax/blob/develop/examples/experimental/noise_channel_monte_carlo.ipynb) - statevector simulation of circuit noise using the Monte-Carlo/quantum trajectories approach. + The [pytket](https://github.com/CQCL/pytket) repository also contains `tk_to_qujax` implementations for some of the above at [pytket-qujax_classification.ipynb](https://github.com/CQCL/pytket/blob/main/examples/pytket-qujax-classification.ipynb), [pytket-qujax_heisenberg_vqe.ipynb](https://github.com/CQCL/pytket/blob/main/examples/pytket-qujax_heisenberg_vqe.ipynb) and [pytket-qujax_qaoa.ipynb](https://github.com/CQCL/pytket/blob/main/examples/pytket-qujax_qaoa.ipynb). \ No newline at end of file diff --git a/docs/experimental.rst b/docs/experimental.rst new file mode 100644 index 0000000..54e4209 --- /dev/null +++ b/docs/experimental.rst @@ -0,0 +1,8 @@ +Experimental +======================= + +.. toctree:: + :titlesonly: + + Pure state simulation + diff --git a/docs/experimental/statetensor.rst b/docs/experimental/statetensor.rst new file mode 100644 index 0000000..904a550 --- /dev/null +++ b/docs/experimental/statetensor.rst @@ -0,0 +1,13 @@ +Pure state simulation +======================= + +.. toctree:: + :titlesonly: + + statetensor/get_default_gates + statetensor/get_default_operations + statetensor/get_params_to_statetensor_func + statetensor/get_params + statetensor/parse_op + statetensor/wrap_parameterised_tensor + diff --git a/docs/experimental/statetensor/get_default_gates.rst b/docs/experimental/statetensor/get_default_gates.rst new file mode 100644 index 0000000..4a67da7 --- /dev/null +++ b/docs/experimental/statetensor/get_default_gates.rst @@ -0,0 +1,4 @@ +get_default_gates +============================================== + +.. autofunction:: qujax.experimental.statetensor.get_default_gates diff --git a/docs/experimental/statetensor/get_default_operations.rst b/docs/experimental/statetensor/get_default_operations.rst new file mode 100644 index 0000000..9abad9e --- /dev/null +++ b/docs/experimental/statetensor/get_default_operations.rst @@ -0,0 +1,4 @@ +get_default_operations +============================================== + +.. autofunction:: qujax.experimental.statetensor.get_default_operations diff --git a/docs/experimental/statetensor/get_params.rst b/docs/experimental/statetensor/get_params.rst new file mode 100644 index 0000000..ce5b878 --- /dev/null +++ b/docs/experimental/statetensor/get_params.rst @@ -0,0 +1,4 @@ +get_params +============================================== + +.. autofunction:: qujax.experimental.statetensor.get_params diff --git a/docs/experimental/statetensor/get_params_to_statetensor_func.rst b/docs/experimental/statetensor/get_params_to_statetensor_func.rst new file mode 100644 index 0000000..62408f0 --- /dev/null +++ b/docs/experimental/statetensor/get_params_to_statetensor_func.rst @@ -0,0 +1,4 @@ +get_params_to_statetensor_func +===================================== + +.. autofunction:: qujax.experimental.statetensor.get_params_to_statetensor_func diff --git a/docs/experimental/statetensor/parse_op.rst b/docs/experimental/statetensor/parse_op.rst new file mode 100644 index 0000000..c9c9474 --- /dev/null +++ b/docs/experimental/statetensor/parse_op.rst @@ -0,0 +1,4 @@ +parse_op +============================================== + +.. autofunction:: qujax.experimental.statetensor.parse_op diff --git a/docs/experimental/statetensor/wrap_parameterised_tensor.rst b/docs/experimental/statetensor/wrap_parameterised_tensor.rst new file mode 100644 index 0000000..9383e53 --- /dev/null +++ b/docs/experimental/statetensor/wrap_parameterised_tensor.rst @@ -0,0 +1,4 @@ +wrap_parameterised_tensor +============================================== + +.. autofunction:: qujax.experimental.statetensor.wrap_parameterised_tensor diff --git a/docs/index.rst b/docs/index.rst index e3d50bb..9a1c982 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -59,6 +59,7 @@ Contents Pure state simulation Mixed state simulation Utility functions + Experimental .. toctree:: :caption: Links: diff --git a/docs/statetensor/apply_gate.rst b/docs/statetensor/apply_gate.rst index a7486b0..e100783 100644 --- a/docs/statetensor/apply_gate.rst +++ b/docs/statetensor/apply_gate.rst @@ -1,4 +1,4 @@ apply_gate ============================================== -.. autofunction:: qujax.apply_gate +.. autofunction:: qujax.experimental.statetensor.apply_gate diff --git a/examples/experimental/noise_channel_monte_carlo.ipynb b/examples/experimental/noise_channel_monte_carlo.ipynb new file mode 100644 index 0000000..a55f54a --- /dev/null +++ b/examples/experimental/noise_channel_monte_carlo.ipynb @@ -0,0 +1,628 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Noise channel w/ quantum trajectories (Monte Carlo)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Imports" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n" + ] + } + ], + "source": [ + "from copy import copy\n", + "from typing import Union\n", + "\n", + "import numpy as np\n", + "\n", + "from jax import numpy as jnp\n", + "from jax import vmap, jit, value_and_grad\n", + "from jax.random import PRNGKey, choice\n", + "\n", + "import qujax\n", + "from qujax import (\n", + " get_params_to_densitytensor_func,\n", + " get_statetensor_to_expectation_func,\n", + " get_densitytensor_to_expectation_func,\n", + " all_zeros_statetensor,\n", + ")\n", + "from qujax.experimental.statetensor import get_params_to_statetensor_func, get_default_gates\n", + "\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Motivation\n", + "\n", + "Here, we will aim to simulate a quantum circuit suffering from a noise channel using only statevector simulation, as opposed to using density matrices as is done in [this notebook](https://github.com/CQCL/qujax/blob/develop/examples/noise_channel.ipynb). \n", + "\n", + "A reason to do this is that density matrices require quadratically more resources to store than statevectors. Indeed, for $N$ qubits, a generic statevector requires $2^N$ complex numbers to specify, while for density matrices it requires $2^{2N}$ complex numbers. This makes it so that for a fixed amount of memory, statevector simulation can be performed for twice the number of qubits, as the following plot illustrates:\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGwCAYAAABFFQqPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACJlElEQVR4nOzdd1xV9f/A8ddlbxAHQ5YLFJmKOXOLYu6FVubIysrUHJWVOSpHue1n2VAb35Q0MTN37tRKBRwoLhQHgjiYsu49vz+uXiQcoMBlvJ+PBw+9n3Puue9zD3DffM77vI9KURQFIYQQQohKyEDfAQghhBBC6IskQkIIIYSotCQREkIIIUSlJYmQEEIIISotSYSEEEIIUWlJIiSEEEKISksSISGEEEJUWkb6DqCs02g0XL16FWtra1Qqlb7DEUIIIUQhKIpCamoqzs7OGBg8fN5HEqHHuHr1Kq6urvoOQwghhBBP4NKlS7i4uDx0uSRCj2FtbQ1o30gbGxs9RyOEEEKIwkhJScHV1VX3Of4wkgg9xr3TYTY2NpIICSGEEOXM48papFhaCCGEEJWWJEJCCCGEqLQkERJCCCFEpSU1QsVAo9GQnZ2t7zCEqJSMjY0xNDTUdxhCiHJKEqGnlJ2dTWxsLBqNRt+hCFFp2dnZ4ejoKL2+hBBFJonQU1AUhfj4eAwNDXF1dX1kwyYhRPFTFIWMjAwSExMBcHJy0nNEQojyRhKhp5Cbm0tGRgbOzs5YWFjoOxwhKiVzc3MAEhMTqVGjhpwmE0IUiUxhPAW1Wg2AiYmJniMRonK794dITk6OniMRQpQ3kggVA6lLEEK/5GdQCPGkJBESQgghRKUliZAQQgghKi1JhAQAbdu2ZezYsaX2elOnTiUgIKDUXq80XbhwAZVKRWRkZIlsv7TeuxUrVmBnZ1dmtiOEECVBEiGhFxMmTODPP//UPR46dCi9evXSX0CPUZQPc1dXV+Lj4/Hx8SmRWP773pUlHh4eLFiwIN9YaGgop0+f1k9AQogy7U62mr/OJuk1Brl8XuiFlZUVVlZW+g6j2GVnZ2NiYoKjo2OJvUZ5e+/Mzc11l7gLIcQ9MddSeWvlEWKT0gl/oyU+NW31EofMCBUjRVHIyM7Vy5eiKIWOMz09nZdeegkrKyucnJyYO3dugXWys7N55513qFmzJpaWljRt2pRdu3bplt+bIdmyZQsNGjTAysqKLl26EB8fr1tn165dPPPMM1haWmJnZ0fLli25ePEikP/0ztSpU/n+++/57bffUKlUqFQqdu3aRfv27Rk1alS+uG7cuIGpqSk7dux44L7d2+6yZctwc3PDysqK119/HbVazWeffYajoyM1atTg008/zfe8efPm4evri6WlJa6urrzxxhukpaXp9mPYsGEkJyfr4ps6dSqgnQH55JNPGDp0KLa2trzyyisFTo1Nnz4dZ2dnbty4oXu9Hj160Lp164d2JC/sewd5s2kzZszAwcEBOzs7pk2bRm5uLhMnTsTe3h4XFxeWLVuWb/sqlYrbt2/rxiIjI1GpVFy4cOGBMZ07d46ePXvi4OCAlZUVTZo0Yfv27brlbdu25eLFi7z99tu69wkePJv25ZdfUqdOHUxMTPDy8uLHH3/Mt1ylUvHtt9/Su3dvLCwsqFevHuvXr39gXEKI8kVRFH48cIEeX+zjdEIadhYmpGfl6i0emREqRndy1Hh/tEUvrx09vTMWJoU7nBMnTmTnzp2Eh4fj6OjI+++/z+HDh/N9uA4bNowLFy6watUqnJ2dCQ8Pp0uXLhw7dox69eoBkJGRwZw5c/jxxx8xMDDgxRdfZMKECfzvf/8jNzeXXr168corr7By5Uqys7P5559/HniZ84QJEzh58iQpKSksX74cAHt7e0aMGMGoUaOYO3cupqamAPzvf//D2dmZdu3aPXT/zp07x6ZNm9i8eTPnzp2jX79+xMbG4unpye7du9m/fz/Dhw+nQ4cONGvWDAADAwMWLVqEh4cHsbGxvPHGG7zzzjssWbKEFi1asGDBAj766CNiYmIA8s3IfP7550yePJkPP/zwgfF88MEHbN68mREjRhAeHs5XX33Fnj17iIqKemA38qK8d/fs2LEDFxcX9uzZw19//cXLL7/MgQMHaN26NX///TdhYWGMHDmSTp064erq+tDtPEpaWhpdu3blk08+wczMjO+//57u3bsTExODm5sba9euxd/fn1dffZVXXnnlodsJDw9nzJgxLFiwgI4dO7JhwwaGDRuGi4tLvuM6bdo0PvvsMz7//HMWL17MCy+8wMWLF7G3t3+i+IUQ+ncrPZt3fj3KtugEANp5VWdOf3+qWpnqLSZJhCqZtLQ0vvvuO3744Qc6deoEwPfff4+Li4tunXPnzrFy5UouX76Ms7MzoE1WNm/ezPLly5kxYwagbV731VdfUadOHQBGjRrF9OnTAUhJSSE5OZlu3brpljdo0OCBMVlZWWFubk5WVla+U0p9+/blrbfe4rfffmPAgAEALF++nKFDhz4yKdBoNCxbtgxra2u8vb1p164dMTExbNy4EQMDA7y8vJg9eza7du3SJUL3F4rXqlWLjz/+mNdff50lS5ZgYmKCra0tKpXqgae82rdvz4QJE3SP/zujYmhoyE8//URAQADvvfceixcv5uuvv8bd3f2B8RflvbvH3t6eRYsW6fbvs88+IyMjg/fffx+ASZMmMWvWLP766y8GDhz4yG09jL+/P/7+/rrHn3zyCeHh4axfv55Ro0Zhb2+PoaEh1tbWjzw1OGfOHIYOHcobb7wBwLhx4zh48CBz5szJlwgNHTqUQYMGATBjxgwWL17MP//8Q5cuXZ4ofiGEfh08f4OxqyK5lpKJiaEB74XUZ1hLD733AZNEqBiZGxsSPb2z3l67MM6dO0d2djbNmzfXjdnb2+Pl5aV7fOTIERRFwdPTM99zs7KyqFq1qu6xhYWF7oMatPd5unfPJ3t7e4YOHUrnzp3p1KkTHTt2ZMCAAUW6F5SpqSkvvvgiy5YtY8CAAURGRhIVFcW6dese+TwPDw+sra11jx0cHDA0NMw3++Lg4KCLFWDnzp3MmDGD6OhoUlJSyM3NJTMzk/T0dCwtLR/5ekFBQY/dl9q1azNnzhxee+01QkNDeeGFFx667pO8dw0bNiywf/cXaxsaGlK1atV8+1xU6enpTJs2jQ0bNnD16lVyc3O5c+cOcXFxRdrOyZMnefXVV/ONtWzZkoULF+Yb8/Pz0/3f0tISa2vrp4pfCKEfuWoNi/48w+KdZ1EUqF3dkkUDA/VWE/RfFapGqHfv3lSpUoV+/frlGzcyMiIgIICAgABGjBhRYq+vUqmwMDHSy1dhM+rC1BJpNBoMDQ05fPgwkZGRuq+TJ0/m+7AyNjYusP/3b3/58uUcOHCAFi1aEBYWhqenJwcPHizku6k1YsQItm3bxuXLl1m2bBkdOnR46EzKo+J60Ni9+pyLFy/StWtXfHx8+PXXXzl8+DD/93//BxTulg2PS5Tu2bNnD4aGhly4cIHc3EefDy/qe1fUfb6XNN1/vB63rxMnTuTXX3/l008/Ze/evURGRuLr60t2dvYjn/cg//1+VRSlwNij4hdClA+XbmYQ+vVBFu3QJkGhQa5seKtVmUmCoIIlQqNHj+aHH34oMG5nZ6f7MP/222/1EFnZUbduXYyNjfN9qN66dSvf5c2BgYGo1WoSExOpW7duvq+iXg0VGBjIpEmT2L9/Pz4+Pvz8888PXM/ExER377b7+fr6EhQUxDfffMPPP//M8OHDi/T6hXHo0CFyc3OZO3cuzZo1w9PTk6tXrxYqvsIKCwtj7dq17Nq1i0uXLvHxxx8/9jmFfe+eRPXq1QHyFbc/ru/R3r17GTp0KL1798bX1xdHR8cCpwEL8z41aNCAffv25Rvbv3//Y0//CSHKlz+OxtN10V4OX7yFtakRiwcFMrufX6HrWUtLhUqE2rVrl++UiCjIysqKl19+mYkTJ/Lnn39y/Phxhg4dmu+0iqenJy+88AIvvfQSa9euJTY2ln///ZfZs2ezcePGQr1ObGwskyZN4sCBA1y8eJGtW7dy+vTph37YeXh4cPToUWJiYkhKSso3OzFixAhmzZqFWq2md+/eT/cGPECdOnXIzc1l8eLFnD9/nh9//JGvvvqqQHxpaWn8+eefJCUlkZGRUejtX758mddff53Zs2fTqlUrVqxYwcyZMx86w1PU9+5J1K1bF1dXV6ZOncrp06f5448/Hnj14H+fs3btWt0pyueff77ADI2Hhwd79uzhypUrJCU9uDfIxIkTWbFiBV999RVnzpxh3rx5rF27Nl+dlRCi/MrIzuW9X4/y5s9HSM3MJdDNjo1jnqW7v7O+Q3ugMpMI7dmzh+7du+Ps7IxKpXpgHciSJUuoVasWZmZmNG7cmL179xZq2ykpKTRu3JhWrVqxe/fuYo68/Pn8889p3bo1PXr0oGPHjrRq1YrGjRvnW2f58uW89NJLjB8/Hi8vL3r06MHff/9d6CuOLCwsOHXqFH379sXT05NXX32VUaNG8dprrz1w/VdeeQUvLy+CgoKoXr06f/31l27ZoEGDMDIy4vnnn8fMzOzJd/whAgICmDdvHrNnz8bHx4f//e9/zJw5M986LVq0YOTIkYSGhlK9enU+++yzQm1bURSGDh3KM888o2sF0KlTJ0aNGsWLL76ou0T/fkV9756EsbExK1eu5NSpU/j7+zN79mw++eSTRz5n/vz5VKlShRYtWtC9e3c6d+5Mo0aN8q0zffp0Lly4QJ06dXSzTv/Vq1cvFi5cyOeff07Dhg1ZunQpy5cvp23btsW1e0IIPYm+mkL3xftY9e8lVCoY1a4uv7zWHFd7C32H9lAqpSgNaErQpk2b+Ouvv2jUqBF9+/YlPDw8X6fhsLAwBg8ezJIlS2jZsiVLly7l22+/JTo6Gjc3N916u3bt4osvvmDNmjW6satXr+Ls7Mzx48d57rnnOHbsGDY2NoWKKyUlBVtbW5KTkws8JzMzk9jYWF1yJkrGpUuX8PDw4N9//y3wwSsEyM+iEPqmKAor9l9g5sZTZKs1ONiYMj80gBZ1quktpkd9ft+vzJyoCwkJISQk5KHL582bx8svv6wrdl6wYAFbtmzhyy+/LPDX+3/duwTcx8cHb29vTp8+/dArfbKyssjKytI9TklJKequiGKSk5NDfHw87733Hs2aNZMkSAghyqAbaVm8s+Yof57SXtXZsUENPuvnj72liZ4jK5wyc2rsUbKzszl8+DDBwcH5xoODg9m/f/8jn3vr1i1dYnP58mWio6OpXbv2Q9efOXMmtra2uq8nbT4nnt5ff/2Fu7s7hw8fLlCzI4QQQv/2n00iZOFe/jyViImRAdN6NOSbl4LKTRIEZWhG6FGSkpJQq9U4ODjkG3dwcODatWu6x507d+bIkSOkp6fj4uJCeHg4OTk5vPbaaxgYGKBSqVi4cOEjO9NOmjSJcePG6R6npKRIMqQnbdu2LdKtQ4QQQpSOHLWG+dtO8+XucygK1K1hxeJBgTRwKlzZSVlSLhKhex7Xe2TLlgff3uLYsWOFfg1TU1Pd7RyEEEIIkV/cjQxGr4og8tJtAAY948ZH3bwxNylcY9+yplwkQtWqVcPQ0DDf7A9AYmJigVkiIYQQQpSM3yKv8EH4cdKycrExM2JWXz+6+hb+jgFlUbmoETIxMaFx48Zs27Yt3/i2bdto0aKFnqISQgghKof0rFwmrI5izKpI0rJyCXKvwsYxz5b7JAjK0IxQWloaZ8+e1T2OjY0lMjISe3t73NzcGDduHIMHDyYoKIjmzZvz9ddfExcXx8iRI/UYtRBCCFGxHb+SzFsrI4hNSsdABaPa12N0+7oYGZaLuZTHKjOJ0KFDh/LdefpewfKQIUNYsWIFoaGh3Lhxg+nTpxMfH4+Pjw8bN2587H2nhBBCCFF0Go3Csr9imb35FDlqBSdbM+aHBtCsdtXHP7kcKTOJUGGuEHrjjTd44403SikiIZ6eh4cHY8eOZezYsSX6Om3btiUgIIAFCxaUie0IIcq3pLQsJqyOYlfMdQCCvR2Y3dePKuXosvjCKjOJkCgbhg4dyu3btx94i5NHmTp1KuvWrXvsjTtL24ULF6hVqxYREREEBASU+uv/+++/hb47fWnatWsX7dq149atW9jZ2enG165dW+Cu70KIymXP6euM+yWKpLQsTI0MmNzNmxeauhW4cruikERIiELKyckpcpLwsPttlVWP6rElhKjYsnM1zN0aw9I95wHwdLBi8aBGeDlW7JuZV4xKJ1Eka9aswdfXF3Nzc6pWrUrHjh1JT09n6tSpfP/99/z222+oVCpUKhW7du0C4N1338XT0xMLCwtq167N5MmTdXeIX7FiBdOmTSMqKkr3vBUrVgCQnJzMq6++So0aNbCxsaF9+/ZERUUBEBMTg0ql4tSpU/nimzdvHh4eHrpTpdHR0XTt2hUrKyscHBwYPHhwvjubazQaZs+eTd26dTE1NcXNzY1PP/0UgFq1agEQGBiISqXS3dhTo9Ewffp0XFxcMDU1JSAggM2bN+u2eeHCBVQqFb/88gtt27bFzMyMn3766YHv59SpU3Fzc8PU1BRnZ2dGjx6tW+bh4ZHvNJNKpWLp0qV069YNCwsLGjRowIEDBzh79ixt27bF0tKS5s2bc+7cOd1zhg4dmu++ewBjx4595E1Kf/rpJ4KCgrC2tsbR0ZHnn3+exMRE3b7dq8erUqUKKpWKoUOHAtpTY/efxrt16xYvvfQSVapUwcLCgpCQEM6cOaNbvmLFCuzs7NiyZQsNGjTAysqKLl26EB8f/9DYhBBlz4WkdPp9tV+XBL3YzI31o1pV+CQIJBEqXooC2en6+SpkB+b4+HgGDRrE8OHDOXnyJLt27aJPnz4oisKECRMYMGCA7oMsPj5e157A2tqaFStWEB0dzcKFC/nmm2+YP38+AKGhoYwfP56GDRvqnhcaGoqiKDz33HNcu3aNjRs3cvjwYRo1akSHDh24efMmXl5eNG7cmP/973/5Yvz55595/vnnUalUxMfH06ZNGwICAjh06BCbN28mISGBAQMG6NafNGkSs2fPZvLkyURHR/Pzzz/r+kv9888/AGzfvp34+HjWrl0LwMKFC5k7dy5z5szh6NGjdO7cmR49euT7kAdtAjh69GhOnjxJ586dC7yfa9asYf78+SxdupQzZ86wbt06fH19H3kMPv74Y1566SUiIyOpX78+zz//PK+99hqTJk3i0KFDALo71T+p7OxsPv74Y6Kioli3bh2xsbG6ZMfV1ZVff/0V0Caj8fHxLFy48IHbGTp0KIcOHWL9+vUcOHAARVHo2rWrLgkGyMjIYM6cOfz444/s2bOHuLg4JkyY8FTxCyFKT3jEZZ5btJejl5OxNTfmqxcb80kvX8yMy2eDxKKSU2PFKScDZjjr57Xfvwomj69FiY+PJzc3lz59+uiuuLv/g9vc3JysrCwcHR3zPe/DDz/U/d/Dw4Px48cTFhbGO++8g7m5OVZWVhgZGeV73o4dOzh27BiJiYm6bt1z5sxh3bp1rFmzhldffZUXXniBL774go8//hiA06dPc/jwYX744QcAvvzySxo1asSMGTN02122bBmurq6cPn0aJycnFi5cyBdffMGQIUMAqFOnDq1atQLyTk1VrVo1X2xz5szh3XffZeDAgQDMnj2bnTt3smDBAv7v//5Pt97YsWPp06fPQ9/PuLg4HB0d6dixI8bGxri5ufHMM8888hgMGzZMl8i9++67NG/enMmTJ+sSrTFjxjBs2LBHbuNxhg8frvt/7dq1WbRoEc888wxpaWlYWVnpToHVqFEjX43Q/c6cOcP69ev566+/dAnx//73P1xdXVm3bh39+/cHtKcMv/rqK+rUqQNok7jp06c/VfxCiJKXlpXLR+uOszbiCgDP1LJnQWgAznbmeo6sdMmMUCXj7+9Phw4d8PX1pX///nzzzTfcunXrsc9bs2YNrVq1wtHRESsrKyZPnkxcXNwjn3P48GHS0tKoWrUqVlZWuq/Y2FjdqZ+BAwdy8eJFDh48CGg/aAMCAvD29tZtY+fOnfmeX79+fQDOnTvHyZMnycrKokOHDoV+D1JSUrh69SotW7bMN96yZUtOnjyZbywoKOiR2+rfvz937tyhdu3avPLKK4SHh5Obm/vI5/j5+en+f2/m6v5k1MHBgczMTFJSUgq1Pw8SERFBz549cXd3x9raWnca7XHH7H4nT57EyMiIpk2b6saqVq2Kl5dXvvfJwsJClwQBODk56U7DCSHKpqhLt3lu0V7WRlzBQAXjOnmy8pVmlS4JApkRKl7GFtqZGX29diEYGhqybds29u/fz9atW1m8eDEffPABf//9t66e5r8OHjzIwIEDmTZtGp07d8bW1pZVq1Yxd+7cR76WRqPByclJV2d0v3uzEE5OTrRr146ff/6ZZs2asXLlSl577bV82+jevTuzZ88usA0nJyfOnz9fqP1+kMfduw547BVfrq6uxMTEsG3bNrZv384bb7zB559/zu7dux9aWH3/+L3Xe9CYRqMBwMDAoEBriftPTf1Xeno6wcHBBAcH89NPP1G9enXi4uLo3Lkz2dnZj9yf+z2sncV/36f/7qdKpZKb5QpRRmk0Ct/sPc/nW2LI1SjUtDNnwcAAmnhU3gslJBEqTipVoU5P6ZtKpaJly5a0bNmSjz76CHd3d8LDwxk3bhwmJiao1ep86//111+4u7vzwQcf6MYuXryYb50HPa9Ro0Zcu3YNIyMjPDw8HhrPCy+8wLvvvsugQYM4d+6c7nTVvW38+uuveHh4YGRU8Nu1Xr16mJub8+effzJixIgCy01MtD0v7o/NxsYGZ2dn9u3bR+vWrXXj+/fvf+xprQcxNzenR48e9OjRgzfffJP69etz7NgxGjVqVORtPUj16tU5fvx4vrHIyMiHJlqnTp0iKSmJWbNm4erqCqCrPbrnQe/Lf3l7e5Obm8vff/+tOzV248YNTp8+TYMGDZ54f4QQ+pGYmsn4X6LYe0Z7sUmIjyOz+vhha1G5W2bIqbFK5u+//2bGjBkcOnSIuLg41q5dy/Xr13UfbB4eHhw9epSYmBiSkpLIycmhbt26xMXFsWrVKs6dO8eiRYsIDw/Pt10PDw/dbVGSkpLIysqiY8eONG/enF69erFlyxYuXLjA/v37+fDDD/N9MPfp04eUlBRef/112rVrR82aNXXL3nzzTW7evMmgQYP4559/OH/+PFu3bmX48OGo1WrMzMx49913eeedd/jhhx84d+4cBw8e5LvvvgO0NTDm5ua6Iuvk5GQAJk6cyOzZswkLCyMmJob33nuPyMhIxowZU6T3c8WKFXz33XccP36c8+fP8+OPP2Jubl6sHc/bt2/PoUOH+OGHHzhz5gxTpkwpkBjdz83NDRMTExYvXsz58+dZv369rgbrHnd3d1QqFRs2bOD69eukpaUV2E69evXo2bMnr7zyCvv27SMqKooXX3yRmjVr0rNnz2LbPyFEydsZk0jIgr3sPZOEmbEBM/v4suSFRpU+CQJAEY+UnJysAEpycnKBZXfu3FGio6OVO3fu6CGyJxMdHa107txZqV69umJqaqp4enoqixcv1i1PTExUOnXqpFhZWSmAsnPnTkVRFGXixIlK1apVFSsrKyU0NFSZP3++Ymtrq3teZmam0rdvX8XOzk4BlOXLlyuKoigpKSnKW2+9pTg7OyvGxsaKq6ur8sILLyhxcXH54urfv78CKMuWLSsQ8+nTp5XevXsrdnZ2irm5uVK/fn1l7NixikajURRFUdRqtfLJJ58o7u7uirGxseLm5qbMmDFD9/xvvvlGcXV1VQwMDJQ2bdronjNt2jSlZs2airGxseLv769s2rRJ95zY2FgFUCIiIh75foaHhytNmzZVbGxsFEtLS6VZs2bK9u3bdcvd3d2V+fPn6x4DSnh4+CNfZ+fOnQqg3Lp1Szf20UcfKQ4ODoqtra3y9ttvK6NGjdLti6IoSps2bZQxY8boHv/888+Kh4eHYmpqqjRv3lxZv359gdeZPn264ujoqKhUKmXIkCEP3M7NmzeVwYMHK7a2toq5ubnSuXNn5fTp07rly5cvz/d9cO89Ke1fLeXxZ1GI0pCZk6tM//2E4v7uBsX93Q1K5/m7ldPXUvQdVql41Of3/VSKIifzHyUlJQVbW1uSk5OxsbHJtywzM5PY2Fhq1aqFmZmZniIUQsjPohAFnb+exlsrIzhxVXvhxZDm7kzq2qDSXBb/qM/v+0mNkBBCCFGBKIrCmsOXmbL+BBnZaqpYGPNZP386eTvoO7QySRIhIYQQooJIyczhw/DjrI/SXsHcrLY9C0IDcbSVmdKHkURICCGEqACOxN1izKoILt28g6GBinGdPBnZpg6GBhXzZqnFRRIhIYQQohzTaBS+3H2OedtOo77bG2jRoEAau1fRd2jlgiRCxUDqzYXQL/kZFJVVQkomb4dFsv/cDQC6+TnxaW9fbM3lsvjCkkToKRgaaivvs7OzMTevfG3JhSgrMjIygIJdroWoyP48mcCE1VHcysjB3NiQaT0a0j/IpUCHfPFokgg9BSMjIywsLLh+/TrGxsYYGEh/SiFKk6IoZGRkkJiYiJ2dne6PEyEqsswcNbM2nWLF/gsAeDvZsGhQIHVrWOk3sHJKEqGnoFKpcHJyIjY2tsAtJ4QQpcfOzg5HR0d9hyFEiTubqO0NdDJe2xtoeMtavBvihamR/BHwpCQRekomJibUq1evSDezFEIUH2NjY5kJEhWeoij8cugSU9dHcydHjb2lCXP6+9G+vvQGelqSCBUDAwMD6WYrhBCiRCTfyeH98GP8cTQegJZ1qzJ/QAA1bORzpzhIIiSEEEKUUYcv3mT0ykiu3L6DkYGK8cFevNa6NgbSG6jYSCIkhBBClDFqjcKXu84yf/sZ1BoFN3sLFg0KJMDVTt+hVTiSCAkhhBBlSHzyHd4Oi+Tg+ZsA9Axw5pNePlibSXuIkiCJkBBCCFFGbD1xjXd+PcrtjBwsTAz5uKcPfRrVlN5AJUgSISGEEELPMnPUzNh4kh8OaFux+NS0YfGgRtSqZqnnyCo+SYSEEEIIPTqTkMpbKyM4dS0VgFeercXEzvUxMZImvaVBEiEhhBBCDxRFYeU/l5i+4QSZORqqWZkwd0AAbTyr6zu0SkUSISGEEKKUJWfk8N7ao2w6fg2A1p7Vmdvfn+rWpnqOrPKpUPNuvXv3pkqVKvTr1083dunSJdq2bYu3tzd+fn6sXr1ajxEKIYSo7P69cJOQhXvYdPwaxoYqPujagBVDm0gSpCcqRVEUfQdRXHbu3ElaWhrff/89a9asASA+Pp6EhAQCAgJITEykUaNGxMTEYGlZuAK0lJQUbG1tSU5OxsbGpiTDF0IIUYGpNQpf7DjLwj9Po1HAo6q2N5Cfi52+Q6uQCvv5XaFOjbVr145du3blG3NycsLJyQmAGjVqYG9vz82bNwudCAkhhBBP6+rtO4wNi+SfWG1voD6NajK9pw9WphXqY7hcKjOnxvbs2UP37t1xdnZGpVKxbt26AussWbKEWrVqYWZmRuPGjdm7d2+RXuPQoUNoNBpcXV2LKWohhBDi0TYfv0bIwr38E3sTK1MjFoQGMG9AgCRBZUSZSYTS09Px9/fniy++eODysLAwxo4dywcffEBERATPPvssISEhxMXFFWr7N27c4KWXXuLrr78uzrCFEEKIB8rMUfNB+DFG/nSY5Ds5+LvY8sfoVvQKrKnv0MR9ykw6GhISQkhIyEOXz5s3j5dffpkRI0YAsGDBArZs2cKXX37JzJkzH7ntrKwsevfuzaRJk2jRosVj183KytI9TklJKcJeCCGEEBBzLZW3Vh7hdEIaAK+1qc34Tl7SG6gMKhdHJDs7m8OHDxMcHJxvPDg4mP379z/yuYqiMHToUNq3b8/gwYMf+1ozZ87E1tZW9yWn0YQQQhSWoij8ePAiPb7Yx+mENKpbm/LTy02ZFNJAkqAyqszMCD1KUlISarUaBweHfOMODg5cu3ZN97hz584cOXKE9PR0XFxcCA8PJysri7CwMPz8/HR1Rz/++CO+vr4PfK1JkyYxbtw43eOUlBRJhoQQQjzWrfRs3v31KFujEwBo51WdOf39qWoll8WXZeUiEbrnvzedUxQl39iWLVse+DyNRlPo1zA1NcXUVL5phRBCFN7B8zd4OyyS+ORMTAwNeC+kPsNaesjNUsuBcpEIVatWDUNDw3yzPwCJiYkFZomEEEKI0pKr1rDozzN8sfMsGgVqV7Nk0aBAfGra6js0UUjl4oSliYkJjRs3Ztu2bfnGt23b9tjiZyGEEKIkXL6VQejXB1m0Q5sEDQhy4fe3WkkSVM6UmRmhtLQ0zp49q3scGxtLZGQk9vb2uLm5MW7cOAYPHkxQUBDNmzfn66+/Ji4ujpEjR+oxaiGEEJXRH0fjeW/tUVIzc7E2NeLTPr708HfWd1jiCZSZROjQoUO0a9dO9/hewfKQIUNYsWIFoaGh3Lhxg+nTpxMfH4+Pjw8bN27E3d1dXyELIYSoZDKyc/l4QzQr/7kEQKCbHYsGBuJqb6HnyMSTqlD3GisJcq8xIYQQANFXU3hr5RHOXU9HpYI329ZlTMd6GBuWiyqTSqdS3mtMCCGEKG6KovD9/gvM2HSK7FwNDjamzA8NoEWdavoOTRQDSYSEEEKIh7iZns07a6LYfjIRgI4NavBZP3/sLU30HJkoLpIICSGEEA+w/2wSb/8SSUJKFiZGBnzQtQEvNXeX3kAVjCRCQgghxH1y1BoWbD/Nkl3nUBSoU92SxYMa4e0sdaIVkSRCQgghxF2XbmYwelUEEXG3ARj0jCuTu3ljYSIflxWVHFkhhBACWB91lQ/WHiM1KxcbMyNm9fWjq6+TvsMSJUwSISGEEJVaRnYuU347werDlwEIcq/CgoEBuFSR3kCVgSRCQgghKq3jV5IZvTKC80npGKhgVPt6jG5fFyPpDVRpSCIkhBCi0lEUhWV/XWD2plNkqzU42pixYGAAzWpX1XdoopRJIiSEEKJSSUrLYuLqKHbGXAcg2NuB2X39qCK9gSolSYSEEEJUGnvPXGfcL1FcT83C1MiAD7t582JTN+kNVIlJIiSEEKLCy1FrmLM1hqW7zwPg6WDF4kGN8HK01nNkQt8kERJCCFGhXbyRzuiVEURdTgbghaZufPicN+YmhnqOTJQFkggJIYSosNZFXOHDdcdJy8rF1tyY2X396OLjqO+wRBkiiZAQQogKJy0rl49+O87aI1cAeMbDngUDA3C2M9dzZKKskURICCFEhXL08m1Gr4zgwo0MDFQwpoMno9rXxdBACqJFQZIICSGEqBA0GoXv9sXy2ZZT5KgVnG3NWDgokCYe9voOTZRhkggJIYQo9xJTMxn/SxR7zyQBEOLjyKw+fthaGOs5MlHWSSIkhBCiXNsVk8iE1VEkpWVjZmzAlO4NGdjEVXoDiUKRREgIIUS5lJ2r4fMtp/hmbywA9R2tWTwokHoO0htIFJ4kQkIIIcqd2CRtb6BjV7S9gYY0d2dS1waYGUtvIFE0kggJIYQoNxRF4dcjV/jot+NkZKuxszDm837+dPJ20HdoopySREgIIUS5kJqZw4frjvNb5FUAmtW2Z0FoII62ZnqOTJRnkggJIYQo8yLibjFmVSRxNzMwNFDxdsd6vN5WegOJpyeJkBBCiDJLo1FYuuc8c7fGkKtRqGlnzqJBgTR2r6Lv0EQFIYmQEEKIMikxJZNxv0Sx76y2N9Bzfk7M6O2Lrbn0BhLFRxIhIYQQZc6OUwlMWH2Um+nZmBsbMq1HQ/oHuUhvIFHsJBESQghRZmTlqpm16RTL/7oAgLeTDYsGBVK3hpV+AxMVliRCQgghyoSziWmMXhlBdHwKAMNaevBeSH1MjaQ3kCg5BvoOoDTMmTOHhg0b4uPjw08//aTvcIQQQtxHURR++fcS3RfvIzo+BXtLE5YNDWJK94aSBIkSV6QZoZiYGFauXMnevXu5cOECGRkZVK9encDAQDp37kzfvn0xNTUtqVifyLFjx/j55585fPgwAB06dKBbt27Y2dnpNzAhhBAk38nhg/BjbDgaD0DLulWZPyCAGjbSG0iUjkLNCEVERNCpUyf8/f3Zs2cPTZo0YezYsXz88ce8+OKLKIrCBx98gLOzM7NnzyYrK6uk4y60kydP0qJFC8zMzDAzMyMgIIDNmzfrOywhhKj0Dl+8xXOL9rLhaDxGBire7VKfH4c3lSRIlKpCJUK9evWiZ8+eXL16lR07djBr1izeeustRowYwTvvvMMPP/xAbGwsGzZs4PDhw8ydO7fYAtyzZw/du3fH2dkZlUrFunXrCqyzZMkSatWqhZmZGY0bN2bv3r26ZT4+PuzcuZPbt29z+/ZtduzYwZUrV4otPiGEEEWj1ih8seMMA5Ye4PKtO7jam7N6ZHNeb1sHA2mQKEpZoU6NnTlzBhMTk8eu17x5c5o3b052dvZTB3ZPeno6/v7+DBs2jL59+xZYHhYWxtixY1myZAktW7Zk6dKlhISEEB0djZubG97e3owePZr27dtja2tLkyZNMDKSGnEhhNCHa8mZjA2L4OD5mwD08Hfm094+WJtJbyChHypFUZSn3cjt27dLpeZGpVIRHh5Or169dGNNmzalUaNGfPnll7qxBg0a0KtXL2bOnFlgGyNGjKB3794899xzD3yNrKysfKf2UlJScHV1JTk5GRsbm+LbGSGEqGS2RScwcU0UtzNysDAxZHpPH/o2qim9gUSJSElJwdbW9rGf30W+amz27NmEhYXpHg8YMICqVatSs2ZNoqKinizaJ5Sdnc3hw4cJDg7ONx4cHMz+/ft1jxMTEwFtsfc///xD586dH7rNmTNnYmtrq/tydXUtmeCFEKKSyMxRM+W347zywyFuZ+TgU9OGDW+1ol9jaZAo9K/IidDSpUt1ycG2bdvYtm0bmzZtIiQkhIkTJxZ7gI+SlJSEWq3GwcEh37iDgwPXrl3TPe7Vqxfe3t68+OKLLF++/JGnxiZNmkRycrLu69KlSyUWvxBCVHRnElLp9X9/8f2BiwC88mwt1r7ektrVpUGiKBuKXCwTHx+vS4Q2bNjAgAEDCA4OxsPDg6ZNmxZ7gIXx378oFEXJN3b/7NDjmJqalrkWAEIIUd4oisKqfy8x7fcTZOZoqGZlwpz+/rT1qqHv0ITIp8gzQlWqVNHNkmzevJmOHTsC2m96tVpdvNE9RrVq1TA0NMw3+wPaU2H/nSUSQghROpIzcnjz5yNMWnuMzBwNrT2rs2lMa0mCRJlU5ESoT58+PP/883Tq1IkbN24QEhICQGRkJHXr1i32AB/FxMSExo0bs23btnzj27Zto0WLFqUaixBCCDh04SZdF+1l47FrGBuq+KBrA1YMbUJ1a5lpF2VTkU+NzZ8/Hw8PDy5dusRnn32GlZX2PG98fDxvvPFGsQeYlpbG2bNndY9jY2OJjIzE3t4eNzc3xo0bx+DBgwkKCqJ58+Z8/fXXxMXFMXLkyGKPRQghxINpewOdZeGfp9Eo4FHVgkWDAvFzsdN3aEI8UrFcPl+Sdu3aRbt27QqMDxkyhBUrVgDahoqfffYZ8fHx+Pj4MH/+fFq3bl0sr1/Yy++EEKKyunr7DmPDIvknVtsbqE+jmkzv6YOVqfRsE/pT2M/vJ0qEfvzxR5YuXcr58+c5cOAA7u7uLFiwgFq1atGzZ8+nCryskURICCEebvPxa7z761GS7+RgaWLIp7196RVYU99hCVFyfYS+/PJLxo0bR0hICLdv39YVSNvZ2bFgwYInDlgIIUT5kZmj5sN1xxj502GS7+Tg72LLH6OflSRIlDtFToQWL17MN998wwcffIChoaFuPCgoiGPHjhVrcEIIIcqemGup9PhiHz8djANgZJs6rB7ZAo9qlnqOTIiiK/IJ3NjYWAIDAwuMm5qakp6eXixBCSGEKHsUReGnv+P4ZEM0WbkaqlubMm+AP8/Wq67v0IR4YkVOhGrVqkVkZCTu7u75xjdt2oS3t3exBSaEEKLsuJ2RzTtrjrI1OgGAtl7VmdPfn2pWclm8KN+KnAhNnDiRN998k8zMTBRF4Z9//mHlypXMnDmTb7/9tiRiFEIIoUcHz9/g7bBI4pMzMTZU8V5IA4a18MDAQO4TJsq/IidCw4YNIzc3l3feeYeMjAyef/55atasycKFCxk4cGBJxCiEEEIPctUaFv15hi92nkWjQO1qliwaFIhPTVt9hyZEsXmqPkJJSUloNBpq1Ki4bdPl8nkhRGV0+VYGY1ZFcvjiLQAGBLkwpXtDLKU3kCgnSuzy+fbt23P79m1Ae6+ve0lQSkoK7du3f7JohRBClBl/HI0nZOFeDl+8hbWpEYsGBfJZP39JgkSFVOTv6l27dpGdnV1gPDMzk7179xZLUEIIIUpfRnYuH2+IZuU/2htrB7rZsWhgIK72FnqOTIiSU+hE6OjRo7r/R0dH57vju1qtZvPmzdSsKY20hBCiPIq+msJbK49w7no6KhW80bYOYzt6YmxY5BMHQpQrhU6EAgICUKlUqFSqB54CMzc3Z/HixcUanBBCiJKlKArf77/AjE2nyM7VUMPalAWhAbSoW03foQlRKgqdCMXGxqIoCrVr1+aff/6hevW8BlomJibUqFEjX6dpIYQQZdvN9GzeWRPF9pOJAHSoX4PP+/tjb2mi58iEKD2FToTuNVBMTU3F0lLaqAshRHm2/2wSb/8SSUJKFiZGBnzQtQEvNXdHpZLeQKJyKfLJXwcHB4YPH86+fftKIh4hhBAlKEet4fMtp3jhu79JSMmiTnVL1r3RkiEtPCQJEpVSkROhlStXkpycTIcOHfD09GTWrFlcvXq1JGITQghRjC7dzGDA0gP8385zKAoMesaV399qhbez9EgTldcTN1S8ceMGP/zwAytWrCA6OprOnTszfPhwevTogZFRxek1IQ0VhRAVwe9RV3l/7TFSs3KxNjNiVh8/nvNz0ndYQpSYwn5+P1Vn6XsWL17MxIkTyc7Oplq1aowcOZL33nsPC4vy33tCEiEhRHmWkZ3LlN9OsPrwZQAau1dh4cAAXKqU/9/PQjxKYT+/n3jq5tq1a/zwww8sX76cuLg4+vXrx8svv8zVq1eZNWsWBw8eZOvWrU+6eSGEEE/p+JVkRq+M4HxSOgYqGNWuLqM71MNIegMJoVPkRGjt2rUsX76cLVu24O3tzZtvvsmLL76InZ2dbp2AgAACAwOLM04hhBCFpCgKy/66wOxNp8hWa3C0MWPBwACa1a6q79CEKHOe6O7zAwcO5K+//qJJkyYPXKd27dp88MEHTx2cEEKIoklKy2Li6ih2xlwHINjbgdl9/agivYGEeKAi1whlZGRUiNqfwpIaISFEebH3zHXG/RLF9VRtb6DJ3bx5sambXBYvKqUSqxGysLBArVYTHh7OyZMnUalU1K9fn169elWoq8WEEKK8yFFrmLM1hqW7zwPg6WDFokGB1HeUP96EeJwiZy7Hjx+nR48eJCQk4OXlBcDp06epXr0669evx9fXt9iDFEII8WAXb6QzemUEUZeTAXihqRsfPueNuYnc8kiIwihyIjRixAh8fHw4fPgwVapUAeDWrVsMHTqUV199lQMHDhR7kEIIIQpaF3GFD9cdJy0rF1tzY2b39aOLj6O+wxKiXClyIhQVFcWhQ4d0SRBAlSpV+PTTTx9aPC2EEKL4pGXl8tFvx1l75AoAz3jYs2BgAM525nqOTIjyp8iJkJeXFwkJCTRs2DDfeGJiInXr1i22wIQQQhR09PJtRq+M4MKNDAxUMKaDJ6Pa18XQQAqihXgShUqEUlJSdP+fMWMGo0ePZurUqTRr1gyAgwcPMn36dGbPnl0yUQohRCWn0Sh8ty+Wz7acIket4GxrxsJBgTTxsNd3aEKUa4W6fN7AwCDf5Zf3nnJv7P7HarW6JOLUG7l8Xgihb4mpmYz/JYq9Z5IACPFxZFYfP2wtjPUcmRBlV7FePr9z585iC0wIIUTh7T59nfG/RJKUlo2ZsQFTujdkYBNX6Q0kRDEpVCLUpk2bko6jRMXGxjJ8+HASEhIwNDTk4MGDWFpa6jssIYR4qOxcDZ9vOcU3e2MBqO9ozeJBgdRzsNZzZEJULE/cATEjI4O4uDiys7Pzjfv5+T11UMVt6NChfPLJJzz77LPcvHkTU1NTfYckhBAPFZuk7Q107Iq2N9CQ5u5M6toAM2PpDSREcStyInT9+nWGDRvGpk2bHri8rNUInThxAmNjY5599lkA7O2lsFAIUTYpisLaI1eY/NtxMrLV2FkY83k/fzp5O+g7NCEqLIOiPmHs2LHcunWLgwcPYm5uzubNm/n++++pV68e69evL/YA9+zZQ/fu3XF2dkalUrFu3boC6yxZsoRatWphZmZG48aN2bt3r27ZmTNnsLKyokePHjRq1IgZM2YUe4xCCPG0UjNzeDsskvGro8jIVtOstj2bx7SWJEiIElbkGaEdO3bw22+/0aRJEwwMDHB3d6dTp07Y2Ngwc+ZMnnvuuWINMD09HX9/f4YNG0bfvn0LLA8LC2Ps2LEsWbKEli1bsnTpUkJCQoiOjsbNzY2cnBz27t1LZGQkNWrUoEuXLjRp0oROnToVa5xCCPGkIi9pewPF3czA0EDFuE6ejGxTR3oDCVEKipwIpaenU6NGDUB7mun69et4enri6+vLkSNHij3AkJAQQkJCHrp83rx5vPzyy4wYMQKABQsWsGXLFr788ktmzpyJi4sLTZo0wdXVFYCuXbsSGRn50EQoKyuLrKws3eP7eygJIURx0mgUlu45z9ytMeRqFGrambNoUCCN3as8/slCiGJR5FNjXl5exMTEABAQEMDSpUu5cuUKX331FU5OTsUe4KNkZ2dz+PBhgoOD840HBwezf/9+AJo0aUJCQgK3bt1Co9GwZ88eGjRo8NBtzpw5E1tbW93XvQRKCCGKU2JKJi8t+4fZm0+Rq1Ho5ufExjHPShIkRCkr8ozQ2LFjiY+PB2DKlCl07tyZ//3vf5iYmLBixYriju+RkpKSUKvVODjkP4fu4ODAtWvXADAyMmLGjBm0bt0aRVEIDg6mW7duD93mpEmTGDdunO5xSkqKJENCiGK181Qi41dHcTM9G3NjQ6b1aEj/IBfpDSSEHhQ5EXrhhRd0/w8MDOTChQucOnUKNzc3qlWrVqzBFdZ/f3koipJv7HGn1+5namoql9cLIUpEVq6aWZtOsfyvCwB4O9mwaFAgdWtY6TcwISqxJ+4jdI+FhQWNGjUqjliKrFq1ahgaGupmf+5JTEwsMEskhBD6dO56Gm/9HEF0vLbucFhLD94LqY+pkfQGEkKfipwI3X/a6H4qlQozMzPq1q1Lz549S6Vfj4mJCY0bN2bbtm307t1bN75t2zZ69uxZ4q8vhBCPoygKqw9dZsr6E9zJUWNvacKc/n60ry9/rAlRFhQ5EYqIiODIkSOo1Wq8vLxQFIUzZ85gaGhI/fr1WbJkCePHj2ffvn14e3s/dYBpaWmcPXtW9zg2NpbIyEjs7e1xc3Nj3LhxDB48mKCgIJo3b87XX39NXFwcI0eOfOrXFkKIp5GSmcP7a4+x4ai2rrJl3arMHxBADRszPUcmhLinyInQvdme5cuX6+7mmpKSwssvv0yrVq145ZVXeP7553n77bfZsmXLUwd46NAh2rVrp3t8b0ZqyJAhrFixgtDQUG7cuMH06dOJj4/Hx8eHjRs34u7u/tSvLYQQT+rwxVuMWRXB5Vt3MDJQMT7Yi9da18ZAegMJUaaoFEVRivKEmjVrsm3btgKzPSdOnCA4OJgrV65w5MgRgoODSUpKKtZg9SElJQVbW1uSk5N1iZ8QQjyMWqPw1e5zzNt2GrVGwdXenEUDAwl0k8vihShNhf38LvKMUHJyMomJiQUSoevXr+uaD9rZ2RW4GasQQlR015IzeTsskgPnbwDQw9+ZT3r7YGNmrOfIhBAP80SnxoYPH87cuXNp0qQJKpWKf/75hwkTJtCrVy8A/vnnHzw9PYs7ViGEKLO2RycwcU0UtzJysDAxZHpPH/o2qim9gYQo44p8aiwtLY23336bH374gdzcXEDbtHDIkCHMnz8fS0tLIiMjAW3n6fJOTo0JIR4lM0fNzI0n+f7ARQB8atqwaGAgtatLbyAh9Kmwn99FToTuSUtL4/z58yiKQp06dbCyqpg/9JIICSEe5kxCKm+tjODUtVQARrSqxcQuXtIbSIgyoMRqhO6xsrLCz8/vSZ8uhBDllqIorPr3EtN+P0FmjoZqVibM6e9PW68a+g5NCFFEhUqE+vTpw4oVK7CxsaFPnz6PXHft2rXFEpgQQpRFyRk5TAo/ysZj2o72z9arxtwB/tSwlt5AQpRHhUqEbG1tdQV/tra2JRqQEEKUVYcu3GTMqkiu3Nb2BnqnixcjWklvICHKsyeuEaospEZICKHWKPzfzrMs2H4ajQLuVS1YNDAQf1c7fYcmhHiIEq8REkKIyuDq7TuMDYvkn9ibAPQOrMnHvXywMpVfn0JUBAaFXfHcuXMMHz5c99jNzQ17e3vdV/Xq1YmJiSmRIIUQQh82H79GyMK9/BN7E0sTQ+YN8Gd+aIAkQUJUIIX+aV68eDGOjo66x7du3eKjjz6iRg3tVRJhYWHMnz+fr776qvijFEKIUpSZo+aTP6L56WAcAH4utiwaGIhHNUs9RyaEKG6FToS2b9/O4sWL84317duX2rVrA+Dh4cGIESOKNzohhChlMddSeWvlEU4npAHwWuvajA/2wsSo0BPoQohypNCJ0MWLF6lVq5bu8YgRI/JdQebh4cHly5eLNzohhCgliqLw099xfLIhmqxcDdWsTJk3wJ/WntX1HZoQogQVOhEyMDAgMTERDw8PAObPn59veUJCAsbGcmNBIUT5czsjm3fWHGVrdAIAbTyrM3eAP9WsTPUcmRCipBU6EWrYsCHbt2/nmWeeeeDyLVu24OPjU2yBCSFEafj7/A3GhkUSn5yJsaGKd7vUZ3jLWtIbSIhKotCJ0LBhwxg7diz+/v4899xz+Zb9/vvvzJo1iwULFhR3fEIIUSJy1RoW7TjLFzvOoFGgVjVLFg8KxKemNI0VojIpdCL0yiuvsGPHDrp37079+vXx8vJCpVJx6tQpYmJi6Nu3L6+88kpJxiqEEMXi8q0Mxq6K5NDFWwD0a+zCtB4NsZTL4oWodIrcWXrVqlWsWrWK06dPA1CvXj0GDRrEwIEDSyRAfZPO0kJULBuPxfPer0dJyczFytSIT3v70DOgpr7DEkIUs8J+fsstNh5DEiEhKoY72WqmbzjByn8uAeDvasfigYG4VbXQc2RCiJIgt9gQQoi7oq+mMHpVBGcT01CpYGSbOozr5ImxofQGEqKyk0RICFFhKYrCDwcu8unGk2Tnaqhhbcr80ABa1q2m79CEEGWEJEJCiArpZno276yJYvvJRADa16/B5/38qCq9gYQQ95FESAhR4ew/l8TbYZEkpGRhYmjApK71GdrCA5VKegMJIfKTREgIUWHkqDUs2H6aJbvOoShQu7q2N1BDZ+kNJESZpc4BQ/3dmaLIiVB6ejqzZs3izz//JDExEY1Gk2/5+fPniy04IYQorEs3Mxi9KoKIuNsAhAa5MqWHNxYm8veeEGWORg2xuyEqDM5shdFHwLyKXkIp8m+IESNGsHv3bgYPHoyTk5NMNQsh9O73qKu8v/YYqVm5WJsZMbOPL938nPUdlhDiv64dh6Or4NgaSI3PG4/ZDAGD9BJSkROhTZs28ccff9CyZcuSiEcIIQotIzuXqetP8MuhywA0crNj4cBAXO2lN5AQZUZKPBxbDUfDIOF43riZHfj0Ab+B4Prg+5iWhiInQlWqVMHe3r4kYhFCiEI7fiWZ0asiOH89HZUKRrWry5gO9TCS3kBC6F9WGpzaAFGrtKfAlLtlNAbG4NkZ/AdCvWAw0v9VnEVOhD7++GM++ugjvv/+eyws5K8uIUTpUhSFZX9dYPamU2SrNTjYmLIgNJDmdarqOzQhKjeNGs7v0s78nPwdcjLylrk2Bb9QaNgbLMrWZEqRE6G5c+dy7tw5HBwc8PDwwNg4f6X3kSNHii244pCamkr79u3JyclBrVYzevRouTmsEOVUUloWE1dHsTPmOgAdGzjwWT8/7C1N9ByZEJXYtWPamZ9jayDtWt64fW3taS+//tr/l1FFToR69epVAmGUHAsLC3bv3o2FhQUZGRn4+PjQp08fqlaVvx6FKE/2nUni7V8iuZ6ahYmRAR8+14DBzdzlgg0h9CHlKhz9RfuVeCJv3LwK+PTVJkAuQVAOfj6LnAhNmTKlJOIoMYaGhrpTeJmZmajVauQ+s0KUHzlqDXO3nmbpHm1voHo1rFg0KJAGTnITZCFKVVaq9pRX1CqI3QPc/Sw1NAHPLtq6n7qdwKh8zdCW+arCPXv20L17d5ydnVGpVKxbt67AOkuWLKFWrVqYmZnRuHFj9u7dm2/57du38ff3x8XFhXfeeYdq1eQ+Q0KUBxdvpNPvqwN8tVubBD3f1I31o1pJEiREaVHnwpnt8OsImOMJ617XFj+jgFtz6LYAJpyG0B+h/nPlLgmCQs4I2dvbc/r0aapVq0aVKlUeORV98+bNYgsOtA0c/f39GTZsGH379i2wPCwsjLFjx7JkyRJatmzJ0qVLCQkJITo6Gjc3NwDs7OyIiooiISGBPn360K9fPxwcHB74ellZWWRlZekep6SkFOv+CCEKZ13EFT5cd5y0rFxszIyY3dePEF8nfYclRMWnKHDtqLbZ4fE1kJaQt8y+jnbmx28AVPHQW4jFqVCJ0Pz587G2tgZgwYIFJRlPASEhIYSEhDx0+bx583j55ZcZMWIEoI1vy5YtfPnll8ycOTPfug4ODvj5+bFnzx769+//wO3NnDmTadOmFd8OCCGKJC0rl49+O87aI1cAaOJRhQUDA6lpZ67nyISo4JKvwLFftAnQ9ZN54+b24NtPW/dTs1G5qPspCpVSjgpmVCoV4eHhuoLt7OxsLCwsWL16Nb1799atN2bMGCIjI9m9ezcJCQmYm5tjY2NDSkoKzZs3Z+XKlfj5+T3wNR40I+Tq6kpycjI2NjIdL0RJOnY5mbdWHuHCjQwMVDC6Qz1GtasrvYGEKCmZKdq6n6OrIHYveXU/puAVcrfup6Ne7wX2pFJSUrC1tX3s53ehZoTS09OxtLQs9IsXdf0nlZSUhFqtLnCay8HBgWvXtJfwXb58mZdffhlFUVAUhVGjRj00CQIwNTXF1FT/DZ6EqEw0GoXv9sXy2ZZT5KgVnG3NWDAwkGdqla1+I0JUCOpcOLdDm/yc2gi5d/KWubfU9vvx7gnmdnoLsTQVKhGqW7cub731FkOHDsXZ+cH371EUhe3btzNv3jxat27NpEmTijXQR/lvzZKiKLqxxo0bExkZWWqxCCGK5npqFuNXR7HntLY3UJeGjszq64udRfkruhSizFIUiI/Mq/tJv563rGo98A8F3wFQxV1vIepLoRKhXbt28eGHHzJt2jQCAgIICgrC2dkZMzMzbt26RXR0NAcOHMDY2JhJkybx6quvlnTcAFSrVg1DQ0Pd7M89iYmJDy2GFkKUHbtPX2f8L5EkpWVjamTAR929ef4ZN+kNJERxuX0pr+4nKSZv3KKatt+Pfyg4V7y6n6IoVCLk5eXF6tWruXz5MqtXr2bPnj3s37+fO3fuUK1aNQIDA/nmm2/o2rUrBgaldy7fxMSExo0bs23btnw1Qtu2baNnz56lFocQomiyczXM2RrD13vOA+DlYM3i5wPxdLDWc2RCVACZyRD9mzb5ubgvb9zIDLy6ak991e1QLut+SkKRGiq6uLjw9ttv8/bbb5dUPAWkpaVx9uxZ3ePY2FgiIyOxt7fHzc2NcePGMXjwYIKCgmjevDlff/01cXFxjBw5stRiFEIUXmxSOqNXRnDsSjIALzV35/2uDTAzNtRzZEKUY+ocOPuntu4nZhPkZuYt83j2bt1PDzCz1V+MZVSRO0uXtkOHDtGuXTvd43HjxgEwZMgQVqxYQWhoKDdu3GD69OnEx8fj4+PDxo0bcXevfOc5hSjrfj18mcm/HScjW42dhTGf9fUjuKGjvsMSonxSFLh65G7dz6+QkZS3rJpXXt2Pnav+YiwHytXl8/pQ2MvvhBAPl5qZw+R1x1kXeRWAprXsWTAwACdb6Q0kRJHduphX93PjTN64ZXXw6adNgJwCKnXdDxTz5fNCCPGkIi/dZvTKCOJuZmBooGJsh3q80a4uhgaV+5e0EEVy57a27udoGFz8K2/cyEx7awu/gVCnndT9PAFJhIQQJUKjUfh673nmbIkhV6NQ086cRYMCaOwuvYGEKJTcbDi7/W7dz2ZQ32v2qwKPVtpmhw16gJmcrXgakggJIYpdYkom436JYt9Zbc3Cc75OzOjji625/LUqxCMpClw5rL3D+/Ff4c599++s3uBu3U9/sHXRX4wVTJGvdffw8GD69OnExcWVRDziMfbs2UP37t1xdnZGpVKxbt26B6538uRJevToga2tLdbW1jRr1uyxxywlJYXJkyfTsGFDzM3NqVq1Kk2aNOGzzz7j1q1buvXatm2LSqVCpVJhYGCAg4MD/fv35+LFi8W5q6Kc2nkqkS4L97LvbBJmxgbM7uvLF88HShIkxKPcugC7P4MvguDbDvDvN9okyLIGNHsTXtsDbxyAVm9LElTMipwIjR8/nt9++43atWvTqVMnVq1ale/eXKJkpaen4+/vzxdffPHQdc6dO0erVq2oX78+u3btIioqismTJ2NmZvbQ59y8eZNmzZqxfPlyJkyYwN9//81ff/3FlClTiIyM5Oeff863/iuvvEJ8fDxXrlzht99+49KlS7z44ovFtp+i/MnKVTP992iGrfiXm+nZNHCyYcNbrQhtIg0ShXigO7fg0HJY1gUW+sPOT+HGWTAy1876vPArjDsJXWaAk3+lL34uKU981VhUVBTLli1j5cqV5Obm8vzzzzN8+HAaNWpU3DHqVVm+auy/N6G9Z+DAgRgbG/Pjjz8WelsjR47kp59+IiYmhpo1axZYfv9tS9q2bUtAQAALFizQLf/xxx8ZOXIk6enpT7Qvonw7dz2Nt36OIDo+BYChLTx4L6S+9AYS4r9ys+HsNu2pr9ObQZ19d4EKarW+W/fTHUyluejTKuzn9xO3gfb392fhwoVcuXKFKVOm8O2339KkSRP8/f1ZtmwZclW+fmg0Gv744w88PT3p3LkzNWrUoGnTpg89hXbvOWFhYbz44osPTIKg4P3c7nfz5k1Wr15N06ZNnzZ8Uc4oisIv/16i26J9RMenUMXCmO+GBDG1R0NJgoS4R1Hg0r/wx3iY6wWrnoeT67VJUA1v6DgN3j4BQ9ZDwPOSBJWyJy6WzsnJITw8nOXLl7Nt2zaaNWvGyy+/zNWrV/nggw/Yvn17gdMpouQlJiaSlpbGrFmz+OSTT5g9ezabN2+mT58+7Ny5kzZt2hR4zvXr17l9+zZeXl75xhs3bkxMjPbeNN27d2flypW6ZUuWLOHbb79FURQyMjLw9PRky5YtJbtzokxJyczh/bXH2HA0HoAWdaoyPzQAB5uHn4IVolK5GQtHf9Fe8n7zXN64lSP49tPO/jj4yCkvPStyInTkyBGWL1/OypUrMTQ0ZPDgwcyfP5/69evr1gkODqZ169bFGqgoHI1GA0DPnj11t0IJCAhg//79fPXVVw9MhO7576xPeHg42dnZvPvuu9y5cyffshdeeIEPPvgAgISEBGbMmEFwcDCHDx/G2lr+mqnoDl+8xZhVEVy+dQdDAxXjgz15rXUd6Q0kxJ1bcCJc2+zw0sG8cWML7Skvv1Co3RYMZMa0rChyItSkSRM6derEl19+Sa9evTA2LngliLe3NwMHDiyWAEXRVKtWDSMjI7y9vfONN2jQgH379j3wOdWrV8fOzo5Tp07lG3dzcwPA2tqa27dv51tma2tL3bp1Aahbty7fffcdTk5OhIWFMWLEiGLaG1HWqDUKX+0+x7xtp1FrFFztzVk4MJBGblX0HZoQ+pObBWe2aut+zmzNq/tRGUCtNtqZn/rdwNRKv3GKBypSIqRWq/nuu+/o0aMH9vYPb4pmaWnJ8uXLnzo4UXQmJiY0adJEd0rrntOnTz/0/msGBgYMGDCAn376icmTJz+0TuhRDA21f938d+ZIVBzXkjN5OyySA+dvANDd35lPe/tgYyaXxYtKSFHg0j/aZofH10Lm7bxlDj7amR/f/mDjpLcQReEUKREyNDRk5MiRtGnT5pGJkCg5aWlpnD17Vvc4NjaWyMhI7O3tdTM4EydOJDQ0lNatW9OuXTs2b97M77//zq5dux663RkzZrBr1y6aNm3K9OnTCQoKwtLSkqNHj3LgwAF8fHzyrZ+RkcG1a9cA7amxTz75BDMzM4KDg4t/p4XebY9OYOKaKG5l5GBhYsi0Hg3p19hFLosXlc+Nc3l1P7di88atnbR1P34DwdHn4c8XZU6RL59v0qQJs2bNokOHDiUVU5lS1i6f37VrF+3atSswPmTIEFasWKF7vGzZMmbOnMnly5fx8vJi2rRp9OzZ85HbTk5OZvbs2YSHhxMbG4uBgQH16tWjZ8+ejB07Vpf8tm3blt27d+ueV6VKFfz8/JgyZcoDYxPlV2aOmpkbT/L9AW2zzIbONiweFEjt6jLFLyqRjJvaLs9Hw+Dyv3njxpbauh//UO0pMKn7KVMK+/ld5ERo69atvPvuu3z88cc0btwYS0vLfMvLQrJQnMpaIiREaTmTkMpbKyM4dS0VgBGtajGxixemRvLLXlQCuVnaPj9RYdq6H02OdlxlALXb3a37eQ5MLB+9HaE3JZYIGRjktR66f1r8XsM9tVr9BOGWXZIIicpGURRW/XuJab+fIDNHQ1VLE+YM8KedVw19hyZEyVIUiDuorfs5EQ6ZyXnLHH21p718+4G1o/5iFIVW2M/vIl81tnPnzqcKrLzKzs7O111ZrVajVqsxMDDAyMgo33oAxsbGxbpuTk4OiqJgZGSkS0Y1Gg25ubmoVKp8V++VhXVzc3PRaDQYGhrqCqmLsq6iKOTkaP8CMzExKZF1H/S+F2VdfR77kvo+uZWWxQfhkWw5kYgaA56tV425A/ypYmZIdnZ2qXyflMaxf9rvk9I69vI7oviP5wPXvX0Bw+Or4WgYyu2L5Nz9aDS2ckbl3x/8BqKu5pX3O+Lu88vK90lp/o4oC98nRTn2hVHkROhRfWgqsrlz5/Lhhx/qTgX+9ddf7Ny5k8DAQHr06KFbb86cOeTk5DBmzBjs7OwA+Pfff9myZQu+vr706dNHt+7ChQvJyMjg9ddfp0YN7V/bkZGRbNiwAS8vr3wtCP7v//6P5ORkRowYobuq6/jx44SHh1O7dm0GDx6sW/ebb77h+vXrDBkyBA8PD0B71VhYWBiurq4MHz5ct+6KFSu4evUqgwYNwtPTE9AWYP/00084ODgwcuRI3br/+9//uHjxIv369aNhw4YAXL58meXLl2Nvb89bb72lW/eXX37hzJkz9OzZk4CAAEDb7HHp0qVYW1szbtw43brh4eFER0cTEhLCM888A2i7VX/xxReYmpry3nvv6dbdsGEDUVFRdOzYkZYtWwKQmprK/PnzMTAwYPLkybp1t2zZwqFDh2jTpg1t27YFICsri9mzZwPw4Ycf6n5Y/vzzTw4cOEDz5s11Bd8ajYaZM2cC8O677+ru1bZ37152795NUFAQzz33nO71Zs+ejUaj4e2339b99XHw4EG2b9+Ov79/vluhzJs3j6ysLEaNGkXVqlUBOHz4MJs2bcLb25v+/fvr1l28eDGpqam89tprODpq/xI9duwYv/32G/Xq1eP555/XrfvVV19x8+ZNhg0bpiueP3nyJGvWrMHd3Z2hQ4fq1v3uu+9ISEjgxRdfpE6dOhy6cJOP/7eDRrkn6GpqgU/7PrzybG0MDFQsW7aMS5cuERoaqusZFhcXx/fff0/16tV54403dNtdtWoV58+fp3fv3vj5+QEQHx/Pt99+i62tLWPHjtWtu2bNGmJiYujWrRuNGzcGtA0+v/zySywsLJg4caJu3fXr13Ps2DE6d+5Ms2bNAG1d28KFCzE2Nub999/Xrbtx40YiIiJo166drqdZRkYGc+bMAWDKlCm6dbdv387ff/9Nq1atdLWPOTk5umM/adIk3S/bXbt2sW/fPpo2bUqXLl1027i37oQJE+R3RHn5HZF+gy2/LONQXAZtlAO05QAAWcb2zM4dCsCHYyZhaKw99n9u3Vrpf0cAnD17lpUrV+Ls7Mwrr7yiW/fHH38sk78jCuOJOkvfvn2b7777jpMnT6JSqfD29mb48OHY2to+yeaEEHqkURQWbj/Dwj9P46zKAlOoW8OK19rU0XdoQhS/E+u0Rc9ntoKmNagCtJ2d63TUnvry6ADzF2vXleLnSqHINUKHDh2ic+fOmJub88wzz6AoCocOHeLOnTts3bq1wt509fr161StWrVMTX2Wt+nMsnZ6RE6N5XD1dgYTfz3BPxduAdA7wInJXb2wNjcpF9PecmpMfkc89nhqNCgX95MTGQYn/8Ak+0beug7+aHwGYOjbB0M750Id+7L4fSKnxh68bokVSz/77LPUrVuXb775Rvfm5ubmMmLECM6fP8+ePXuKsrkyT4qlRUW1+fg13v31KMl3crA0MeTjXj70aeSi77CEKB5JZ7Sdno/9Arfj8sZtXMCvv7bhYY0G+otPlLgSS4TMzc2JiIjId28xgOjoaIKCgsjIyHiyiMsoSYRERZOZo+aTP6L56aD2w8HPxZZFAwPxqCaXAYtyLj1J2+8nahVcPZI3bmIN3j21/X7cW8F9Vz+LiqvErhqzsbEhLi6uQCJ06dIludmmEGVczLVU3lp5hNMJaQC81ro244O9MDGSDwZRTuXcgZhN2rqfs9tBk6sdVxlC3Q7amR+vrmBiod84RZlV5EQoNDSUl19+mTlz5tCiRQtUKhX79u1j4sSJDBo0qCRiFEI8JUVR+OnvOD7ZEE1WroZqVqbMG+BPa8/q+g5NiKLTaCBuv3bmJ/o3yErJW+YUoG126NMXrKT3lXi8IidCc+bMQaVS8dJLL5Gbq828jY2Nef3115k1a1axByiEeDq3M7J5Z81RtkYnANDGszpzB/hTzcpUz5EJUUTXT2ubHR5dDcn31f3YuoLfAO3sT3Uv/cUnyqUi1wjdk5GRwblz51AUhbp162JhUTGnHaVGSJRnf5+/wdiwSOKTMzE2VPFul/oMb1kLAwO5WaooJ9Ku373P1yq4GpE3bmpzt+5nILi1kLofUUCJ1QjdY2Fhga+v75M+XQhRgnLVGhbtOMsXO86gUaBWNUsWDwrEp6b0+hLlQM4dOPWH9i7vZ7eDcvfWTQZGULfj3bqfEDA212+cokIociKUmZnJ4sWL2blzJ4mJiWg0mnzLjxw58pBnCiFKw+VbGYxdFcmhi9reQP0auzCtR0MsTZ/47x4hSp5GAxf3aW9yGv0bZKfmLXNulFf3Y1lNfzGKCqnIvxmHDx/Otm3b6NevH88880y+G68KIfRr47F43vv1KCmZuViZGvFpbx96BtTUd1hCPFziqby6n5TLeeN2btqZH79QqFZPf/GJCq/IidAff/zBxo0bdfdwEULo351sNdM3nGDlP5cA8He1Y/HAQNyqVszaPVHOpSXCsTXaBCg+Km/c1BYa9tLO/rg2k7ofUSqKnAjVrFmz3PUL2rBhA+PHj0ej0fDuu+8yYsQIfYckRLE5GZ/CWysjOJuYhkoFI9vUYVwnT4wN5UNElCHZGXfrflbBuZ35637qBWtnfjy7gLGZfuMUlU6RrxrbtGkTixYt4quvvsLd3b2k4io2ubm5eHt7s3PnTmxsbGjUqBF///039vb2hXq+XDUmyipFUfjhwEU+3XiS7FwNNaxNmR8aQMu6UkMhygiNGi7s1db9nFwP2Wl5y2oGaWd+GvYBy6r6i1FUWCV21VhQUBCZmZnUrl0bCwuLfDdHA7h582bRoy1B//zzDw0bNqRmTW2dRNeuXdmyZYs0fxTl2s30bN5ZE8X2k4kAtK9fg8/7+VFVegOJsiAhOq/uJ/Vq3rid+311P3X1F58Q9ylyIjRo0CCuXLnCjBkzcHBwKPFi6T179vD5559z+PBh4uPjCQ8Pp1evXvnWWbJkCZ9//jnx8fE0bNiQBQsW8OyzzwJw9epVXRIE4OLiwpUrV0o0ZiFK0v5zSbwdFklCShYmhgZM6lqfoS085MIFoV+p1/Lqfq4dyxs3s4WGvcFvILg1A/k+FWVMkROh/fv3c+DAAfz9/UsingLS09Px9/dn2LBh9O3bt8DysLAwxo4dy5IlS2jZsiVLly4lJCSE6Oho3NzceNCZv0d9YGRlZZGVlaV7nJKS8tB1hShNOWoNC7afZsmucygK1KluyaJBgTR0lt5AQk+y0+HkBm3yc34XKHfbqRgYa+t+/EOhXmep+xFlWpETofr163Pnzp2SiOWBQkJCCAkJeejyefPm8fLLL+sKoBcsWMCWLVv48ssvmTlzJjVr1sw3A3T58mWaNm360O3NnDmTadOmFd8OCFEMLt3MYPSqCCLibgMQGuTKlB7eWJhIbyBRyjRqiN19t+7nd8hJz1vm0kR72sunL1gUrg5TCH0rcrH01q1bmTZtGp9++im+vr4FaoRKsqBYpVLlOzWWnZ2NhYUFq1evpnfv3rr1xowZQ2RkJLt37yY3N5cGDRqwa9cuXbH0wYMHqVr1wcV5D5oRcnV1lWJpoTe/R13l/bXHSM3KxdrMiJl9fOnm56zvsERlc+24dubn2BpIjc8br+KhPe3lNwCq1tFbeEL8V4kVS3fp0gWADh065BtXFAWVSoVarS7qJp9YUlISarUaBweHfOMODg5cu3YNACMjI+bOnUu7du3QaDS88847D02CAExNTTE1lYJToX8Z2blMXX+CXw5pm8w1crNj4cBAXO2lN5AoJSnxcGw1HA2DhON542Z24NNHmwC5PiN1P6JcK3IitHPnzpKI46n8t+bnXlJ2T48ePejRo0dphyXEEzt+JZnRqyI4fz0dlQpGtavLmA71MJLeQKKkZaXBqQ0QtUp7Cuz+uh/PztpL3usFg5H8wSgqhiInQm3atCmJOJ5ItWrVMDQ01M3+3JOYmFhglkiI8kBRFJb/dYFZm06RrdbgaGPG/NAAmteRPiuiBGnU2mLno/fqfjLylrk20572athb6n5EhfRElZZ79+5l6dKlnD9/ntWrV1OzZk1+/PFHatWqRatWrYo7xocyMTGhcePGbNu2LV+N0LZt2+jZs2epxSFEcbiRlsXENUfZcUrbG6iTtwOf9fWjiqWJniMTFda1Y9qZn2NrIO2+Pyjta+fV/djX0l98QpSCIidCv/76K4MHD+aFF17gyJEjusLi1NRUZsyYwcaNG4s1wLS0NM6ePat7HBsbS2RkJPb29ri5uTFu3DgGDx5MUFAQzZs35+uvvyYuLo6RI0cWaxxClKR9Z5J4+5dIrqdmYWJkwOTnGvBiM3fpDSSKX8pVOPqL9ivxRN64eRXt1V5+A8ElSOp+RKVR5KvGAgMDefvtt3nppZewtrYmKiqK2rVrExkZSZcuXQqcpnpau3btol27dgXGhwwZwooVKwBtQ8XPPvuM+Ph4fHx8mD9/Pq1bty6W15dbbIiSlKPWMHfraZbu0fYGqlfDisXPB1LfUb7XRDHKSr2v389u4O6vfUMT7f29/AdC3U5gJLOPouIo7Od3kRMhCwsLoqOj8fDwyJcInT9/Hm9vbzIzM586+LJEEiFRUi7eSGf0qkiiLt0G4Pmmbkx+zhtzE0P9BiYqBnXu3bqfVdqbnd5f9+PWXNvvp2Ev7UyQEBVQiV0+7+TkxNmzZ/Hw8Mg3vm/fPmrXrl3kQIWojH6LvMIH4cdJy8rFxsyIz/r50cXHSd9hifJOUeDaUW2zw+NrIC0hb5l9He3Mj98Abe8fIQTwBInQa6+9xpgxY1i2bBkqlYqrV69y4MABJkyYwEcffVQSMQpRYaRl5TLltxP8ekTbG6iJRxUWDAykpp25niMT5VryFTj2izYBun4yb9yial7dT81GUvcjxAMUORF65513SE5Opl27dmRmZtK6dWtMTU2ZMGECo0aNKokYhagQjl3W9gaKTUrHQAWjO9RjVLu60htIPJmsVIherz31FbuXvLofU/AKuVv30xEMjR+5GSEquyLXCN2TkZFBdHQ0Go0Gb29vrKysiju2MkFqhMTT0mgUvtsXy2dbTpGjVnC2NWPBwECeqSU9WUQRqXPh/E7tJe+n/oDc++776N5SW/fj3RPM7fQWohBlRYnVCN1jYWFBUFDQkz5diErhemoW41dHsef0dQC6NHRkVl9f7Czk6hxRSIoC8VHaZofH1kB6Yt6yqvW0d3j3HQBV3PUXoxDlWKEToeHDhxdqvWXLlj1xMEJUJLtPX2f8L5EkpWVjamTAR929ef4ZN+kNJAon+fLdfj9hcP1U3rhFVfDpp02AnKXuR4inVehEaMWKFbi7uxMYGMgTnk0TolLIztUwZ2sMX+85D4CXgzWLnw/E08Faz5GJMi8zBaJ/0yY/F/aRr+6nfldt0XPdDlL3I0QxKnQiNHLkSFatWsX58+cZPnw4L774Ivb2UuMgxP1ik9IZvTKCY1eSAXipuTvvd22AmbH0BhIPoc6Bczu0yc+pPyD3vl5s7q20Mz/ePcHMVn8xClGBFalYOisri7Vr17Js2TL279/Pc889x8svv0xwcHCFne6XYmlRWGuPXGbyuuOkZ6uxszDms75+BDd01HdYoixSFLgakVf3k5GUt6yap7bo2W8A2LnpL0YhyrkS6yx9z8WLF1mxYgU//PADOTk5REdHV8grxyQREo+TmpnD5HXHWRd5FYCmtexZMDAAJ1vpDST+43ZcXt1P0um8cYtq4NtPmwA5B0rdjxDFoMSvGlOpVKhUKhRFQaPRPOlmhCjXIi/dZvTKCOJuZmBooGJsh3q80a4uhgbyQSbuykzW1v1EhcHFfXnjRmZQ/zlt3U+ddlL3I4SeFCkRuv/U2L59++jWrRtffPEFXbp0wcBAmsKJykOjUfh673nmbIkhV6NQ086cRYMCaOwudXMCbd3P2e3afj8xm0CddXeBCjxaaZsdNugBZjLLLIS+FToReuONN1i1ahVubm4MGzaMVatWUbVq1ZKMTYgyKTElk3G/RLHvrLau4zlfJ2b08cXWXP6ir9QUBa4c0XZ6Pv4rZNzIW1a9fl7dj62L/mIUQhRQ6BohAwMD3NzcCAwMfGRh9Nq1a4stuLJAaoTE/XaeSmTC6ihupGdjbmzI1B7eDAhyrbAXC4hCuHXxbt3PKrhxNm/csgb49tcmP07+UvcjRCkr9hqhl156SX7Zi0orK1fNZ5tj+G5fLAANnGxYPCiAujWkN1CldOdWXt1P3P68cSNzaNBNW/dTuy0YPnEZphCilBSpoaIQldG562mMXhnBiaspAAxt4cF7IfWlN1Blk5sNZ7dp635ObwZ19t0FKqjV+m7dT3cwleRYiPJE/lwR4iEURWH14ctM+e0Ed3LUVLEwZk5/fzo0cNB3aKK0KApcPnS37mct3LmZt6yGt7bux7c/2NbUX4xCiKciiZAQD5CSmcMH4cf5PUrbG6hFnarMDw3AwcZMz5GJUnEzNq/fz81zeeNWDnfrfkLB0VfqfoSoACQREuI/jsTdYvTKCC7fuoOhgYrxwZ681rqO9Aaq6O7cghPh2rqfSwfzxo0ttKe8/EK1dT8GckpUiIpEEiEh7lJrFL7afY55206j1ii42puzcGAgjdyq6Ds0UVJys+HMVu2pr9Nb8up+VAZQq4227qd+NzCteF3zhRBakggJASSkZPJ2WCT7z2l7v/Twd+aT3j7YmElvoApHUeDyv9qi5xNrtTNB9zj45NX92DjpL0YhRKmRREhUetujE5i4JopbGTlYmBgyrUdD+jV2kXYRFc3N8/fV/ZzPG7dyBL/+2kveHX30F58QQi8kERKVVmaOmlmbTrFi/wUAGjrbsHhQILWry2mQCiPjprbu52gYXPo7b9zYUlv34x+qPQUmdT9CVFqSCIlK6WxiKqN+juDUtVQARrSqxcQuXpgayQdiuZebpa33ORqm/VeTox1XGWiLnf0GapsemljqNUwhRNkgiZCoVBRFIezfS0z9/QSZORqqWpowZ4A/7bxq6Ds08TQURTvjE7VKOwOUeTtvmYOvdubHtz9YO+otRCFE2SSJkKg0ku/k8P7aY/xxLB6AZ+tVY+4Af2pYS2+gcuvGOe3Mz9EwuHUhb9zaSZv4+A8Eh4Z6C08IUfZJIiQqhUMXbjJmVSRXbt/ByEDFxM5evPJsbQykN1D5k3FTe3f3o2Haq7/uMbGCBj20Nzmt1VrqfoQQhSKJkKjQ1BqF/9t5lgXbT6NRwL2qBYsGBuLvaqfv0ERR5GRq7+919Bdt35/7637qtNfW/dTvKnU/Qogik0RIVFjxyXcYuyqSv2O194fqHViTj3v5YGUq3/blgkaj7fActQqi10Fmct4yRz/taS+ffmAt934TQjy5Cv+JkJqaSvv27cnJyUGtVjN69GheeeUVfYclStiWE9d499ej3M7IwdLEkI97+dCnkYu+wxKFkXRW2+n5aBjcjssbt6mZV/dTo4H+4hNCVCgVPhGysLBg9+7dWFhYkJGRgY+PD3369KFq1ar6Dk2UgMwcNZ/+cZIfD14EwM/FlkUDA/GoJqdMyrT0JO3d3Y+ugiuH88ZNrMG7p/aqL/dWYGCgvxiFEBVShU+EDA0NsbCwACAzMxO1Wo2iKHqOSpSEmGupjF4ZQUyCtjfQa61rMz7YCxMj+fAsk3IyIWajdubn7HbQ5GrHVYZQt4P2VhdeXcHEQr9xCiEqNL1/QuzZs4fu3bvj7OyMSqVi3bp1BdZZsmQJtWrVwszMjMaNG7N3794ivcbt27fx9/fHxcWFd955h2rVqhVT9KIsUBSFnw5epMcX+4hJSKWalSk/DH+GSV0bSBJU1mg0cGEf/DYK5tSDNcO0RdCaXHAKgC6zYPwpeGE1+PaTJEgIUeL0PiOUnp6Ov78/w4YNo2/fvgWWh4WFMXbsWJYsWULLli1ZunQpISEhREdH4+bmBkDjxo3Jysoq8NytW7fi7OyMnZ0dUVFRJCQk0KdPH/r164eDw4MLLLOysvJtKyUlpZj2VJSE2xnZvPvrUbacSACgjWd15vT3p7q1qZ4jE/lcP3237mc1JN9X92Prqr3c3S8UqnvpLz4hRKWlUsrQeSKVSkV4eDi9evXSjTVt2pRGjRrx5Zdf6sYaNGhAr169mDlzZpFf4/XXX6d9+/b079//gcunTp3KtGnTCownJydjY2NT5NcTJefv8zcYGxZJfHImxoYq3u1Sn+Eta0lvoLIi7frdfj+r4GpE3ripzd26n4Hg1kLqfoQQJSIlJQVbW9vHfn7rfUboUbKzszl8+DDvvfdevvHg4GD2799fqG0kJCRgbm6OjY0NKSkp7Nmzh9dff/2h60+aNIlx48bpHqekpODq6vpkOyBKRK5aw+IdZ1m84wwaBWpVs2TxoEB8atrqOzSRc0db9xN1t+5HUWvHDYygbse7dT8hYGyu3ziFEOKuMp0IJSUloVarC5zGcnBw4Nq1a4XaxuXLl3n55ZdRFAVFURg1ahR+fn4PXd/U1BRTUzmtUlZduX2Hsasi+PfCLQD6NXZhWo+GWEpvIP3RaODiPm3Rc/R6yLrvdLJzo7v9fvqCpdTmCSHKnnLx6aFS5T/VoShKgbGHady4MZGRkSUQlShtm47F8+6vR0nJzMXK1IhPe/vQM6CmvsOqvBJP5dX9pFzOG7d1u6/ux1N/8QkhRCGU6USoWrVqGBoaFpj9SUxMfGixs6h47mSrmb4hmpX/aIts/V3tWDwwELeqckVRqUtLhGNrtAlQfFTeuKktNOypvdWFW3Op+xFClBtlOhEyMTGhcePGbNu2jd69e+vGt23bRs+ePfUYmSgtJ+NTeGtlBGcT01CpYGSbOozr5ImxoXzQlprsjLt1P6vg3I7/1P100jY79AwBYzP9ximEEE9A74lQWloaZ8+e1T2OjY0lMjISe3t73NzcGDduHIMHDyYoKIjmzZvz9ddfExcXx8iRI/UYtShpiqLw48GLfPLHSbJzNVS3NmX+gABa1ZM6k1Kh0cCFvXl1P9mpectqBmnrfhr2lrofIUS5p/dE6NChQ7Rr1073+N4VW0OGDGHFihWEhoZy48YNpk+fTnx8PD4+PmzcuBF3d3d9hSxK2K30bCauOcr2k9reQO28tL2BqlpJEXuJSzypnfk5thpSruSN27lra378QqFaXf3FJ4QQxaxM9REqiwrbh0AUjwPnbvB2WCTXUjIxMTTgvZD6DGvpUejiePEEUhPg+BptAnTtaN64ma121sdvILg1AzkGQohypEL0ERKVR65aw8I/z/DFzrMoCtSuru0N1NBZegOViOwMOPWHtuj53A5QNNpxA2OoF6yt+6nXWep+hBAVniRCQu8u3cxgzKoIjsTdBiA0yJUpPbyxMJFvz2KlUUPsHjj6C5xcD9lpectcmmhPe/n0BQt7/cUohBClTD5phF79HnWV98OPkZqZi7WpETP6+NLd31nfYVUsCSfu1v2sgdSreeNVPPLqfqrW0Vt4QgihT5IICb3IyM5l2vpowg5dAiDQzY5FAwNxtZfeQMUi9Zq24DkqDBKO5Y2b2WnrfvwHgmtTqfsRQlR6kgiJUnfiajJvrYzg/PV0VCp4s21dxnSsJ72BnlZ2OpzcoK37Ob8rf92PZ2ftzI9nZzCSq++EEOIeSYREqVEUhRX7LzBz4ymy1RocbEyZHxpAizrSi+aJadTapOdomDYJyknPW+byjLbouWEfqfsRQoiHkERIlIobaVlMXHOUHacSAejYwIHP+vlhb2mi58jKqWvH8up+0u67BU2VWtrTXn4DwL62/uITQohyQhIhUeL+OpvE22GRJKZmYWJkwIfPNWBwM3fpDVRUKVfz6n4ST+SNm1fRzvr4D9Re/SXvqxBCFJokQqLE5Kg1zNt2mq92n0NRoG4NKxYPCqSBkzSmLLSsNDj5+926n93A3f6nhiZ3634Gavv+GMnMmhBCPAlJhESJiLuRwVurIoi6dBuAQc+48lG3hpibGOo3sPJAnQuxu7QzP6c2QE5G3jLXZtq6H+9eUvcjhBDFQBIhUex+i7zCB+HHScvKxcbMiFl9/ejq66TvsMo2RdHe3iIqTHu7i7SEvGX2tbUzP34DwL6W/mIUQogKSBIhUWzSs3KZsv4Eaw5fBiDIvQoLBgbgUkV6Az1U8hU49os2Abp+Mm/c3F7b5dkvFFyCpO5HCCFKiCRColgcv6LtDRSblI6BCka1r8fo9nUxkt5ABWWlQvR6bd1P7F7y1f14hWhnf+p2lLofIYQoBZIIiaei0Sgs+yuW2ZtPkaNWcLI1Y0FoAE1rV9V3aGWLOhfO79Re8n7qD8i9k7fMrUVe3Y+5nb4iFEKISkkSIfHErqdmMWF1FLtPXwegc0MHZvf1w85CZjIAbd1PfJS22eGxNZCemLesat27dT/9tff8EkIIoReSCIknsuf0dcb9EkVSWhamRgZM7ubNC03dpDcQQPJl7R3ej4bB9VN54xZVwaeftu6nZiOp+xFCiDJAEiFRJNm5GuZujWHpnvMAeDpYsXhQI7wcrfUcmZ5lpsDJ9dpTXxf2kVf3Ywr1u96t++kAhsZ6DVMIIUR+kgiJQruQlM7oVREcvZwMwIvN3PjwOW/MjCtpbyB1DpzboZ35OfUH5GbmLXNvdbfupyeY2eovRiGEEI8kiZAolLVHLjN53XHSs9XYmhvzWT8/Ojd01HdYpU9R4GpEXt1PRlLesmqe2tNefgPAzk1/MQohhCg0SYTEI6Vl5TJ53XHCI64A8EwtexaEBuBsZ67nyErZ7bi8up+k03njFtXA927dj3Og1P0IIUQ5I4mQeKioS7cZvSqCizcyMFDB2I6evNmuLoYGleTDPjMZon/TNju8uC9v3MgMvLpqb3Jap73U/QghRDkmiZAoQKNR+GbveT7fEkOuRqGmnTkLBwYQ5FEJ7m2lzoGz27UzPzGb8tf9eDyrnfnx7iF1P0IIUUFIIiTySUzNZPwvUew9o6196erryMzefthaVOBZD0WBK0e0nZ6P/woZN/KWVfPSFj37DgA7V/3FKIQQokRIIiR0dsYkMuGXKG6kZ2NmbMDU7g0JbeJacXsD3bqYV/dz40zeuGV18O2vnf1x8pe6HyGEqMAkERJk5ar5bHMM3+2LBaC+ozVfPB9I3RoVsDfQndsQvU5b9xO3P2/cyBzqP6et+6ndDgzlR0MIISoD+W1fyZ2/nsZbKyM4cTUFgKEtPHgvpH7F6g2Um3237mcVxGwGddbdBSqo9ay22WGD7mBmo9cwhRBClD5JhCopRVFYc/gyU9afICNbTRULYz7v509Hbwd9h1Y8FAWuHNZ2ej7+K9y5mbeseoO7dT/9wdZFfzEKIYTQO0mEKqGUzBw+DD/O+qirADSvXZX5oQE42prpObJicDM2r+7n5rm8ccsa2sTHPxQc/aTuRwghBCCJUKVzJO4WY1ZFcOnmHQwNVIzr5MnINnXKd2+gO7fgRLi27ufSwbxxYwuo301b9Fy7rdT9CCGEKKBSfDLExsYyfPhwEhISMDQ05ODBg1haWuo7rFKl0Sh8ufsc87adRq1RcKlizqJBgTRyq6Lv0J5Mbjac2aqt+zm9BdTZdxeooHabu3U/3cC0AhZ8CyGEKDaVIhEaOnQon3zyCc8++yw3b97E1NRU3yGVqoSUTN4Oi2T/OW1/nG5+Tszo44uNWTnrDaQocPlfbd3PibXamaB7ajTMq/uxcdZfjEIIIcqVCp8InThxAmNjY5599lkA7O0rQXfk+/x5MoEJq6O4lZGDubEh03o2pH9jl/LVG+jm+fvqfs7njVs5au/z5T8QHH31F58QQohyy0DfAezZs4fu3bvj7OyMSqVi3bp1BdZZsmQJtWrVwszMjMaNG7N3795Cb//MmTNYWVnRo0cPGjVqxIwZM4ox+rIrM0fN1PUnePn7Q9zKyMHbyYYNo1sxIKicNEjMuAn/fgffBcOiQNg1U5sEGVtqT3sNDodx0dD5U0mChBBCPDG9zwilp6fj7+/PsGHD6Nu3b4HlYWFhjB07liVLltCyZUuWLl1KSEgI0dHRuLm5AdC4cWOysrIKPHfr1q3k5OSwd+9eIiMjqVGjBl26dKFJkyZ06tTpgfFkZWXl21ZKSkox7WnpOZuYylsrIzkZr419eMtavBvihalRGe8NlJulrfuJWqX9917dj8pAW+zsN1Db9NDUSq9hCiGEqDhUiqIo+g7iHpVKRXh4OL169dKNNW3alEaNGvHll1/qxho0aECvXr2YOXPmY7d54MABpk2bxubNmwH4/PPPAZg4ceID1586dSrTpk0rMJ6cnIyNTdluuKcoCmH/XmLq7yfIzNFQ1dKEOf39aVe/hr5DezhFgUt/3637CYfM23nLHHy1dT8+/cDGSW8hCiGEKH9SUlKwtbV97Oe33meEHiU7O5vDhw/z3nvv5RsPDg5m//79D3lWfk2aNCEhIYFbt25ha2vLnj17eO211x66/qRJkxg3bpzucUpKCq6uZf9mm8l3cnh/7TH+OBYPQKu61Zg3wJ8aNmW0N9CNc9qan6NhcOtC3ri1091+PwPBoaHewhNCCFE5lOlEKCkpCbVajYND/m7HDg4OXLt2rVDbMDIyYsaMGbRu3RpFUQgODqZbt24PXd/U1LTcXVV2+OJNRq+M5MrtOxgZqJjQ2YtXn62NQVnrDZRxU9vl+WiY9uqve0ysoEEP8BsAtVqDQRk/hSeEEKLCKNOJ0D3/Le5VFKVIBb8hISGEhIQUd1h6p9YoLNl5lgV/nkGtUXCvasHCgYEEuNrpO7Q8uVlwerO22eGZraDJ0Y6rDKBO+7t1P13BpHL1dRJCCFE2lOlEqFq1ahgaGhaY/UlMTCwwS1TZxCffYeyqSP6O1d5Dq1eAMx/38sG6LPQGUhSIO6Cd+TkRDpnJecsc/bSnvXz6gXXlPoZCCCH0r0wnQiYmJjRu3Jht27bRu3dv3fi2bdvo2bOnHiPTry0nrvHur0e5nZGDpYkhH/fyoU+jMnDz0KSz2k7PR8PgdlzeuE3NvLqfGg30F58QQgjxH3pPhNLS0jh79qzucWxsLJGRkdjb2+Pm5sa4ceMYPHgwQUFBNG/enK+//pq4uDhGjhypx6j1IzNHzad/nOTHgxcB8K1py6JBgdSqpsfTSuk37tb9rNLe7f0eEyvw7qm9z5dHK6n7EUIIUSbpPRE6dOgQ7dq10z2+d8XWkCFDWLFiBaGhody4cYPp06cTHx+Pj48PGzduxN3dXV8h68XphFTe+jmCmIRUAF5tXZsJwV6YGOmhJ2ZOJpzepK37ObsNNLnacZWhtu7HfyB4dQUTi9KPTQghhCiCMtVHqCwqbB+CkqIoCj//E8f036PJytVQzcqUuQP8aeNZvXQD0Wju1v2sghO/QdZ9dT9O/tqiZ99+YFWGexYJIYSoNCpEH6HK7nZGNu/9eozNJ7TF4m08qzOnvz/VrUvx8v7rp+/W/ayG5PvrflzAr782AapRv/TiEUIIIYqRJEJl1N/nbzA2LJL45EyMDVW807k+L7eqVTq9gdKu59X9XI3IGzex1tb9+IeCeysw0Put6oQQQoinIolQGZOr1rB4x1kW7ziDRoFa1SxZNDAQXxfbkn3hnDsQs/Fu3c92UNTacZUh1O2oTX68uoKxecnGIYQQQpQiSYTKkCu37zB2VQT/XrgFQN9GLkzv2RBL0xI6TBoNXPxLO/MTvR6y7rvBrHOg9rSXT1+wKuV6JCGEEKKUSCJURmw6Fs+7vx4lJTMXK1MjPu3tQ8+AmiXzYtdjtDc5PbYaki/ljdu6am9z4TcQqnuWzGsLIYQQZYgkQnp2J1vNx39E8/Pf2kJkf1c7Fg8MxK1qMV96nnYdjq/RJkDxkXnjprbQ8G6/H7cWUvcjhBCiUpFESI9OxqcwemUEZxLTUKlgZJs6jOvkibFhMSUj2Rnaup+jYXD2z7y6HwMjqNtJW/fjGQLGZfQO9UIIIUQJk0RIDxRF4ceDF/nkj5Nk52qobm3K/AEBtKpX7ek3rtHAhb3a5Cd6PWSn5i2r2fhu3U8fsCyG1xJCCCHKOUmE9CArV8OPBy6SnauhnZe2N1BVq6fsDZR4Mq/uJ+VK3ridm/a0l18oVKv3dK8hhBBCVDCSCOmBmbEhi58PZP/ZGwxr6YFK9YS9gVIT8up+rh3NGze1hYa9tLe6cG0mdT9CCCHEQ0gipCf1HW2o7/gEt+zIzoBTf2gveT+3AxSNdtzACOoFa2d+PLtI3Y8QQghRCJIIlQcatbbuJyoMTq6H7LS8ZTWDtDM/DfuAZVX9xSiEEEKUQ5IIlWUJJ+7W/ayB1Kt543bu2uTHdwBUq6u/+IQQQohyThKhsib1mrbgOSoMEo7ljZvZamd9/AeCa1N40roiIYQQQuhIIlQWZKfDyQ3aup/zu+6r+zEGz8536346g1Ep3nVeCCGEqAQkEdIXjRpid9+t+/kdctLzlrk8o2122LAPWNjrL0YhhBCigpNESB80GljSDJJO541VqXW3388AqFpHf7EJIYQQlYgkQvpgYKCd9UlL1HZ59hsIrs9I3Y8QQghRyiQR0pdO06DbPKn7EUIIIfRIEiF9kXt9CSGEEHon914QQgghRKUliZAQQgghKi1JhIQQQghRaUkiJIQQQohKSxIhIYQQQlRakggJIYQQotKSREgIIYQQlZYkQkIIIYSotCQREkIIIUSlJYmQEEIIISotSYSEEEIIUWlJIiSEEEKISksSISGEEEJUWnL3+cdQFAWAlJQUPUcihBBCiMK697l973P8YSQReozU1FQAXF1d9RyJEEIIIYoqNTUVW1vbhy5XKY9LlSo5jUbD1atXsba2RqVS6TucEpOSkoKrqyuXLl3CxsZG3+GUqMq0r1C59lf2teKqTPsr+1o8FEUhNTUVZ2dnDAweXgkkM0KPYWBggIuLi77DKDU2NjYV/gfvnsq0r1C59lf2teKqTPsr+/r0HjUTdI8USwshhBCi0pJESAghhBCVliRCAgBTU1OmTJmCqampvkMpcZVpX6Fy7a/sa8VVmfZX9rV0SbG0EEIIISotmRESQgghRKUliZAQQgghKi1JhIQQQghRaUkiJIQQQohKSxKhSm7q1KmoVKp8X46OjvoOq1js2bOH7t274+zsjEqlYt26dfmWK4rC1KlTcXZ2xtzcnLZt23LixAn9BPuUHrevQ4cOLXCcmzVrpp9gn9LMmTNp0qQJ1tbW1KhRg169ehETE5NvnYpybAuzrxXp2H755Zf4+fnpmus1b96cTZs26ZZXlOMKj9/XinRc/2vmzJmoVCrGjh2rG9PnsZVESNCwYUPi4+N1X8eOHdN3SMUiPT0df39/vvjiiwcu/+yzz5g3bx5ffPEF//77L46OjnTq1El3f7ny5HH7CtClS5d8x3njxo2lGGHx2b17N2+++SYHDx5k27Zt5ObmEhwcTHp6um6dinJsC7OvUHGOrYuLC7NmzeLQoUMcOnSI9u3b07NnT90HYkU5rvD4fYWKc1zv9++///L111/j5+eXb1yvx1YRldqUKVMUf39/fYdR4gAlPDxc91ij0SiOjo7KrFmzdGOZmZmKra2t8tVXX+khwuLz331VFEUZMmSI0rNnT73EU9ISExMVQNm9e7eiKBX72P53XxWlYh9bRVGUKlWqKN9++22FPq733NtXRamYxzU1NVWpV6+esm3bNqVNmzbKmDFjFEXR/8+szAgJzpw5g7OzM7Vq1WLgwIGcP39e3yGVuNjYWK5du0ZwcLBuzNTUlDZt2rB//349RlZydu3aRY0aNfD09OSVV14hMTFR3yEVi+TkZADs7e2Bin1s/7uv91TEY6tWq1m1ahXp6ek0b968Qh/X/+7rPRXtuL755ps899xzdOzYMd+4vo+t3HS1kmvatCk//PADnp6eJCQk8Mknn9CiRQtOnDhB1apV9R1eibl27RoADg4O+cYdHBy4ePGiPkIqUSEhIfTv3x93d3diY2OZPHky7du35/Dhw+W6e62iKIwbN45WrVrh4+MDVNxj+6B9hYp3bI8dO0bz5s3JzMzEysqK8PBwvL29dR+IFem4PmxfoeId11WrVnHkyBH+/fffAsv0/TMriVAlFxISovu/r68vzZs3p06dOnz//feMGzdOj5GVDpVKle+xoigFxiqC0NBQ3f99fHwICgrC3d2dP/74gz59+ugxsqczatQojh49yr59+wosq2jH9mH7WtGOrZeXF5GRkdy+fZtff/2VIUOGsHv3bt3yinRcH7av3t7eFeq4Xrp0iTFjxrB161bMzMweup6+jq2cGhP5WFpa4uvry5kzZ/QdSom6d2Xcvb9E7klMTCzwV0lF5OTkhLu7e7k+zm+99Rbr169n586duLi46MYr4rF92L4+SHk/tiYmJtStW5egoCBmzpyJv78/CxcurJDH9WH7+iDl+bgePnyYxMREGjdujJGREUZGRuzevZtFixZhZGSkO376OraSCIl8srKyOHnyJE5OTvoOpUTVqlULR0dHtm3bphvLzs5m9+7dtGjRQo+RlY4bN25w6dKlcnmcFUVh1KhRrF27lh07dlCrVq18yyvSsX3cvj5IeT62D6IoCllZWRXquD7MvX19kPJ8XDt06MCxY8eIjIzUfQUFBfHCCy8QGRlJ7dq19XtsS7wcW5Rp48ePV3bt2qWcP39eOXjwoNKtWzfF2tpauXDhgr5De2qpqalKRESEEhERoQDKvHnzlIiICOXixYuKoijKrFmzFFtbW2Xt2rXKsWPHlEGDBilOTk5KSkqKniMvukfta2pqqjJ+/Hhl//79SmxsrLJz506lefPmSs2aNcvlvr7++uuKra2tsmvXLiU+Pl73lZGRoVunohzbx+1rRTu2kyZNUvbs2aPExsYqR48eVd5//33FwMBA2bp1q6IoFee4Ksqj97WiHdcHuf+qMUXR77GVRKiSCw0NVZycnBRjY2PF2dlZ6dOnj3LixAl9h1Usdu7cqQAFvoYMGaIoivaSzSlTpiiOjo6Kqamp0rp1a+XYsWP6DfoJPWpfMzIylODgYKV69eqKsbGx4ubmpgwZMkSJi4vTd9hP5EH7CSjLly/XrVNRju3j9rWiHdvhw4cr7u7uiomJiVK9enWlQ4cOuiRIUSrOcVWUR+9rRTuuD/LfREifx1alKIpS8vNOQgghhBBlj9QICSGEEKLSkkRICCGEEJWWJEJCCCGEqLQkERJCCCFEpSWJkBBCCCEqLUmEhBBCCFFpSSIkhBBCiEpLEiEhhBBCVFqSCAkhypQLFy6gUqmIjIzUdyg6p06dolmzZpiZmREQEFBir1OYfd+1axcqlYrbt2+XWBxCVCaSCAkh8hk6dCgqlYpZs2blG1+3bh0qlUpPUenXlClTsLS0JCYmhj///FOvsbRo0YL4+HhsbW0BWLFiBXZ2dnqNSYjyTBIhIUQBZmZmzJ49m1u3buk7lGKTnZ39xM89d+4crVq1wt3dnapVqxZjVEVnYmKCo6NjpU1KhShukggJIQro2LEjjo6OzJw586HrTJ06tcBpogULFuDh4aF7PHToUHr16sWMGTNwcHDAzs6OadOmkZuby8SJE7G3t8fFxYVly5YV2P6pU6do0aIFZmZmNGzYkF27duVbHh0dTdeuXbGyssLBwYHBgweTlJSkW962bVtGjRrFuHHjqFatGp06dXrgfmg0GqZPn46LiwumpqYEBASwefNm3XKVSsXhw4eZPn06KpWKqVOnPnA76enpvPTSS1hZWeHk5MTcuXNp27YtY8eOzbetdevW5XuenZ0dK1asKPS+339qbNeuXQwbNozk5GRUKlW++JYsWUK9evUwMzPDwcGBfv36PTBuISo7SYSEEAUYGhoyY8YMFi9ezOXLl59qWzt27ODq1avs2bOHefPmMXXqVLp160aVKlX4+++/GTlyJCNHjuTSpUv5njdx4kTGjx9PREQELVq0oEePHty4cQOA+Ph42rRpQ0BAAIcOHWLz5s0kJCQwYMCAfNv4/vvvMTIy4q+//mLp0qUPjG/hwoXMnTuXOXPmcPToUTp37kyPHj04c+aM7rUaNmzI+PHjiY+PZ8KECQ/czsSJE9m5cyfh4eFs3bqVXbt2cfjw4Sd6zx617/dr0aIFCxYswMbGhvj4eF18hw4dYvTo0UyfPp2YmBg2b95M69atnygWISq8UrnHvRCi3BgyZIjSs2dPRVEUpVmzZsrw4cMVRVGU8PBw5f5fGVOmTFH8/f3zPXf+/PmKu7t7vm25u7srarVaN+bl5aU8++yzuse5ubmKpaWlsnLlSkVRFCU2NlYBlFmzZunWycnJUVxcXJTZs2criqIokydPVoKDg/O99qVLlxRAiYmJURRFUdq0aaMEBAQ8dn+dnZ2VTz/9NN9YkyZNlDfeeEP32N/fX5kyZcpDt5GamqqYmJgoq1at0o3duHFDMTc3V8aMGaMbA5Tw8PB8z7W1tVWWL19e6H3fuXOnAii3bt1SFEVRli9frtja2ubb5q+//qrY2NgoKSkpj9l7IYTMCAkhHmr27Nl8//33REdHP/E2GjZsiIFB3q8aBwcHfH19dY8NDQ2pWrUqiYmJ+Z7XvHlz3f+NjIwICgri5MmTABw+fJidO3diZWWl+6pfvz6gree5Jygo6JGxpaSkcPXqVVq2bJlvvGXLlrrXKoxz586RnZ2dL+b/b+duQlLpwjiA//uwhZV9KZEQSSkxlIzUJiIGEcNVrYKoaFObIKiFSVBQizZmREEZgYtWYquggnKTm0oiKCRokxHtClq4iTbl+K5u3MnKMS7clzv/Hww4j2eec85sfDhnxurqajQ3N6vO8bvv5q5Gd3c3Ghoa0NjYiOHhYYTDYby8vPxoLET/OhZCRPQlSZLg8XgwMzOT9V1hYSEymYwi9vr6mtVOp9MpzgsKCj6NybKcczy/HhCWZRk9PT1IJBKKI5lMKraASktLc+b8Pe8vmUwmr4eRP96H7/pRc8++ulat8vJyXF5eIhKJoK6uDnNzcxBFka/cE32ChRARfcvv92N/fx/xeFwRN5lMeHx8VPyw/8n//jk7O3v//Pb2houLi/dVn7a2NlxfX8NiscBqtSoOtcUPABgMBpjNZpycnCji8XgcgiCozmO1WqHT6RRjTqVSuLm5UbQzmUx4eHh4P08mk5+u1Hw3949KSkqQTqez4sXFxXC73QgEAri6usL9/T1isZjqORFpRfHfHgAR/b/Z7XYMDQ1hbW1NEXc6nXh6ekIgEEBfXx+i0SgODw9hMBj+SL/BYBA2mw2CIGBlZQWpVAojIyMAgPHxcYRCIQwMDMDn88FoNOL29hbb29sIhUIoKipS3Y/P58P8/DyamprgcDiwtbWFRCKBcDisOkdZWRlGR0fh8/lQU1OD2tpazM7OKrYEAcDlcmF9fR0dHR2QZRnT09NZq2O55v6RxWLB8/Mzjo6OIIoi9Ho9YrEY7u7uIEkSqqqqcHBwAFmWf7xVR/Qv44oQEeW0sLCQtaUjCAI2NjYQDAYhiiLOz8+/fKPqJ/x+PxYXFyGKIo6Pj7G7uwuj0QgAMJvNOD09RTqdhsfjQWtrKyYnJ1FRUZFVfOQyMTEBr9cLr9cLu92OaDSKvb092Gy2vPIsLS1BkiT09vbC7Xajq6sL7e3tijbLy8uor6+HJEkYHBzE1NQU9Hp9XnP/qLOzE2NjY+jv74fJZEIgEEBlZSV2dnbgcrkgCAI2NzcRiUTQ0tKS15yItKAgo3Zzm4iI8uJ0OuFwOLC6uvq3h0JEX+CKEBEREWkWCyEiIiLSLG6NERERkWZxRYiIiIg0i4UQERERaRYLISIiItIsFkJERESkWSyEiIiISLNYCBEREZFmsRAiIiIizWIhRERERJr1Hy8wOOU55U7IAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "nr_of_qubits = np.array(range(2, 41, 2))\n", + "\n", + "\n", + "def log2_memory_needed_to_store_coefficients(\n", + " nr_of_qubits: Union[np.ndarray, int]\n", + ") -> Union[np.ndarray, float]:\n", + " \"\"\"\n", + " Given `nr_of_qubits`, returns memory needed to store a statevector having that number of qubits.\n", + "\n", + " The logarithm is used to avoid excessively large numbers which can not be stored in a\n", + " np.float64.\n", + " \"\"\"\n", + " log2_nr_of_coefficients = nr_of_qubits\n", + "\n", + " # Get nr of bits needed to store complex numbers representing coefficients\n", + " # Note that 2**5 = 128, which is the size of a complex floating point number\n", + " log2_bits = log2_nr_of_coefficients + 5\n", + "\n", + " # Note that 2**3 = 8, which is the nr of bits in a byte\n", + " log2_bytes = log2_bits - 3\n", + "\n", + " log2_gigabytes = log2_bytes - np.log2(1e9)\n", + "\n", + " return log2_gigabytes\n", + "\n", + "\n", + "# Plot memory needed to store density matrix for different nr of qubits\n", + "plt.plot(\n", + " nr_of_qubits,\n", + " 2 ** log2_memory_needed_to_store_coefficients(2 * nr_of_qubits),\n", + " label=\"density matrix simulation\",\n", + ")\n", + "# Plot memory needed to store statevector for different nr of qubits\n", + "plt.plot(\n", + " nr_of_qubits,\n", + " 2 ** log2_memory_needed_to_store_coefficients(nr_of_qubits),\n", + " label=\"statevector simulation\",\n", + ")\n", + "\n", + "plt.legend(loc=\"upper left\")\n", + "plt.yscale(\"log\")\n", + "plt.ylabel(\"Memory (in Gigabytes)\")\n", + "plt.xlabel(\"Number of qubits\")\n", + "plt.axhline(16, c=\"gray\", ls=\":\")\n", + "plt.text(3, 35, \"16 GB\");" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Theory" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Suppose we wish to apply a $k$-qubit unitary matrix (gate) $U$, but our system is affected by a depolarising noise channel. Suppose also that we are interested in computing the expectation value of some operator $O$. This noise channel, in the density matrix formulation, can be written as\n", + "\n", + "\\begin{equation}\n", + "\\sigma = \\mathcal{D}_{p,U}[\\rho] = p_0 U\\rho U^\\dagger + \\sum_{i=1}^{4^k} p_i P_i \\rho P_i^\\dagger, \\tag{1}\n", + "\\end{equation}\n", + "\n", + "The expectation value we wish to compute is $\\text{tr}(\\sigma O)$. Assuming a pure initial state $\\rho_0 = \\ket{\\psi_0}\\bra{\\psi_0}$, this expectation can be expressed as\n", + "\n", + "\\begin{align}\n", + "&\\text{tr}\\,\\left(\\left[p_0 U\\rho_0 U^\\dagger + \\sum_{i=1}^{4^k} p_i P_i \\rho_0 P_i^\\dagger \\right] O \\right ) \\\\\n", + "&= p_0 \\text{tr}\\,\\left(U\\rho_0 U^\\dagger O \\right ) + \\sum_{i=1}^{4^k} p_i \\text{tr} \\, \\left( P_i \\rho_0 P_i^\\dagger O \\right ) \\\\\n", + "&= p_0 \\text{tr}\\,\\left(U\\rho_0 U^\\dagger O \\right ) + \\sum_{i=1}^{4^k} p_i \\text{tr} \\, \\left( P_i \\rho_0 P_i^\\dagger O \\right ) \\\\\n", + "&= p_0 \\langle U\\ket{\\psi_0} \\rangle_O + \\sum_{i=1}^{4^k} p_i \\langle P_i\\ket{\\psi_0} \\rangle_O \\\\\n", + "\\end{align}\n", + "where $P_i \\in \\{I, X, Y, Z\\}^{\\otimes k}$ and $\\sum_{i=0}^{4^k} p_i = 1$.\n", + "\n", + "The above equality means that we can sample expectation values of several pure state simulations according to the probability vector $[p_0,..., p_{4^k}]$ to arrive at the expectation value of the operator $O$ under the noise model, which would otherwise require mixed state simulation. This technique is sometimes called **quantum trajectories** or **Monte Carlo (wavefunction) simulation**.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that, in practice, we will have multiple noisy gates in our circuit. To apply this technique, every time we encounter a noisy gate, we sample and apply a gate according to the probability vector $[p_0,..., p_{4^k}]$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Simulation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For simplicity, we will assume single qubit gates are noiseless and two qubit gates are affected by a depolarising noise channel with probability vector \n", + "\n", + "$$p = [(1-p_0)/4^k, ..., (1-p_0)/4^k],$$\n", + "\n", + "where $p_0 = 0.99$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Constructing 2-qubit Pauli operators" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First let's construct the set of two qubit Pauli combinations $\\{I, X, Y, Z\\}^{\\otimes 2}$" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'I': Array([[1., 0.],\n", + " [0., 1.]], dtype=float32),\n", + " 'X': Array([[0., 1.],\n", + " [1., 0.]], dtype=float32),\n", + " 'Y': Array([[ 0.+0.j, -0.-1.j],\n", + " [ 0.+1.j, 0.+0.j]], dtype=complex64),\n", + " 'Z': Array([[ 1., 0.],\n", + " [ 0., -1.]], dtype=float32)}" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "paulis = {\n", + " \"I\": qujax.gates.I,\n", + " \"X\": qujax.gates.X,\n", + " \"Y\": qujax.gates.Y,\n", + " \"Z\": qujax.gates.Z,\n", + "}\n", + "paulis" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['II', 'IX', 'IY', 'IZ', 'XI', 'XX', 'XY', 'XZ', 'YI', 'YX', 'YY', 'YZ', 'ZI', 'ZX', 'ZY', 'ZZ']\n" + ] + } + ], + "source": [ + "two_qubit_paulis = {\n", + " a + b: jnp.kron(p1, p2).reshape(2, 2, 2, 2)\n", + " for (a, p1) in paulis.items()\n", + " for (b, p2) in paulis.items()\n", + "}\n", + "two_qubit_pauli_strings = list(two_qubit_paulis.keys())\n", + "\n", + "print(list(two_qubit_paulis.keys()))\n", + "\n", + "# The generated k-qubit Paulis have the correct nr. of elements 4**k, where here k=2\n", + "assert len(two_qubit_paulis) == 4**2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Quantum circuit\n", + "\n", + "Let's define a circuit for our experiments." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "noiseless_gate_seq = [\"H\", \"CX\", \"CX\", \"CX\"]\n", + "noiseless_qubit_seq = [[0], [0, 1], [1, 2], [2, 3]]\n", + "noiseless_param_ind_seq = [None] * len(noiseless_qubit_seq)\n", + "\n", + "nr_of_qubits = len(noiseless_gate_seq)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "q0: -----H-------◯-------------------\n", + " | \n", + "q1: -------------CX------◯-----------\n", + " | \n", + "q2: ---------------------CX------◯---\n", + " | \n", + "q3: -----------------------------CX--\n" + ] + } + ], + "source": [ + "qujax.print_circuit(noiseless_gate_seq, noiseless_qubit_seq, noiseless_param_ind_seq);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let's simulate a noisy version of our quantum circuit. To do so, we will use a (currently experimental but largely stable) version of the `get_params_to_statetensor` function that can perform general operations on the circuit and is not restricted to just applying parameterized gates.\n", + "\n", + "For every 2-qubit gate we encounter, we replace it by a `\"ConditionalGate\"` operation, which applies a gate out of several depending on a parameter passed to the circuit. This will allow us to change the 2-qubit gate to one of the several Pauli gates when noise affects the circuit." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "4 CX gates converted\n", + "noisy_operation_seq=['H', 'ConditionalGate', 'ConditionalGate', 'ConditionalGate']\n", + "noisy_operator_metaparameter_seq=[[0], [['CX', 'II', 'IX', 'IY', 'IZ', 'XI', 'XX', 'XY', 'XZ', 'YI', 'YX', 'YY', 'YZ', 'ZI', 'ZX', 'ZY', 'ZZ'], [0, 1]], [['CX', 'II', 'IX', 'IY', 'IZ', 'XI', 'XX', 'XY', 'XZ', 'YI', 'YX', 'YY', 'YZ', 'ZI', 'ZX', 'ZY', 'ZZ'], [1, 2]], [['CX', 'II', 'IX', 'IY', 'IZ', 'XI', 'XX', 'XY', 'XZ', 'YI', 'YX', 'YY', 'YZ', 'ZI', 'ZX', 'ZY', 'ZZ'], [2, 3]]]\n", + "noisy_param_ind_seq=[None, 0, 1, 2]\n" + ] + } + ], + "source": [ + "noisy_operation_seq = []\n", + "noisy_operator_metaparameter_seq = []\n", + "noisy_param_ind_seq = []\n", + "\n", + "nr_of_2_qubit_gates = 0\n", + "for gate, qubit_inds, param_inds in zip(\n", + " noiseless_gate_seq, noiseless_qubit_seq, noiseless_param_ind_seq\n", + "):\n", + " # The only 2-qubit gates in our circuit are CX gates\n", + " if gate == \"CX\":\n", + " noisy_operation_seq.append(\"ConditionalGate\")\n", + " noisy_operator_metaparameter_seq.append(\n", + " [[\"CX\"] + two_qubit_pauli_strings, qubit_inds]\n", + " )\n", + " noisy_param_ind_seq.append(nr_of_2_qubit_gates)\n", + " nr_of_2_qubit_gates += 1\n", + " else:\n", + " noisy_operation_seq.append(gate)\n", + " noisy_operator_metaparameter_seq.append(qubit_inds)\n", + " noisy_param_ind_seq.append(param_inds)\n", + "\n", + "print(f\"{nr_of_2_qubit_gates+1} CX gates converted\")\n", + "print(f\"{noisy_operation_seq=}\")\n", + "print(f\"{noisy_operator_metaparameter_seq=}\")\n", + "print(f\"{noisy_param_ind_seq=}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "gates = get_default_gates()\n", + "# Add two qubit Paulis to set of available gates\n", + "gates |= two_qubit_paulis" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "params_to_statetensor = get_params_to_statetensor_func(\n", + " noisy_operation_seq,\n", + " noisy_operator_metaparameter_seq,\n", + " noisy_param_ind_seq,\n", + " gate_dict=gates,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Monte-Carlo simulation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now sample several runs of our circuit assuming that each $2$-qubit gate has a $p_0$ probability of being applied and all $2$-qubit Paulis, representing noise, have a $(1-p_0)/4^k$ probability of being applied." + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "nr_of_samples = 100000\n", + "seed = 0\n", + "p0 = 0.99\n", + "pj = 1 - 0.99\n", + "k = 2\n", + "nr_of_paulis = 4**k\n", + "\n", + "observables = [[\"Z\"]]\n", + "observable_indices = [[1]]\n", + "coefficients = [1.0]" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "noise_probability_vector = jnp.array([p0] + [pj] * nr_of_paulis)\n", + "gate_samples = choice(\n", + " PRNGKey(seed),\n", + " jnp.arange(nr_of_paulis + 1),\n", + " (nr_of_samples, nr_of_2_qubit_gates),\n", + " p=noise_probability_vector,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Before we proceed, we note that a sample will have no noise with probability ${p_0}^G$, where $G$ is the number of $2$-qubit gates in the circuit. With $p_0 = 0.99$, even for $60$ $2$-qubit gates this is, on average, over half the samples. This means that we can greatly speed up execution by caching the result of running the circuit with no noise and directly replacing it when the sample is noiseless.\n", + "\n", + "To do this, we separate out the noisy samples and run them separately." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of noiseless samples: 63701\n" + ] + } + ], + "source": [ + "is_noisy = jnp.any(gate_samples, axis=1)\n", + "is_noiseless = jnp.logical_not(is_noisy)\n", + "nr_noiseless_samples = jnp.sum(is_noiseless)\n", + "print(f\"Number of noiseless samples: {nr_noiseless_samples}\")\n", + "\n", + "noisy_samples = gate_samples[is_noisy]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, our observable $O$ will be $Z_1$ i.e. we perform a $Z$ Pauli measurement on the first qubit." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "statetensor_to_expectation = get_statetensor_to_expectation_func(\n", + " observables, observable_indices, coefficients\n", + ")\n", + "\n", + "\n", + "def params_to_expectation(params, statetensor_in):\n", + " return statetensor_to_expectation(params_to_statetensor(params, statetensor_in)[0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can batch over the remaining samples in order to greatly speed up execution." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "batched_params_to_expectation = vmap(params_to_expectation, in_axes=(0, None))\n", + "initial_state = all_zeros_statetensor(nr_of_qubits)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "expectation_vector = batched_params_to_expectation(noisy_samples, initial_state)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, all that is left for us to do is to run the noiseless version of the circuit and compute the final expectation value average." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "noiseless_expectation = params_to_expectation(\n", + " jnp.zeros(nr_of_paulis + 1, dtype=int), initial_state\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Estimated expectation: -5.000019e-05\n" + ] + } + ], + "source": [ + "estimated_expectation = (\n", + " jnp.sum(expectation_vector) + nr_noiseless_samples * noiseless_expectation\n", + ") / nr_of_samples\n", + "print(\"Estimated expectation:\", estimated_expectation)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Confirming the result" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can confirm that this value is correct by using density matrix simulation" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [], + "source": [ + "dm_noisy_gate_seq = []\n", + "for gate in noiseless_gate_seq:\n", + " if gate == \"CX\":\n", + " dm_noisy_gate_seq.append(\n", + " [jnp.sqrt(p0) * gates[\"CX\"]]\n", + " + [jnp.sqrt(pj) * g for g in two_qubit_paulis.values()]\n", + " )\n", + " else:\n", + " dm_noisy_gate_seq.append(gate)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "params_to_densitytensor = qujax.get_params_to_densitytensor_func(\n", + " dm_noisy_gate_seq, noiseless_qubit_seq, noiseless_param_ind_seq\n", + ")\n", + "densitytensor_to_expectation = get_densitytensor_to_expectation_func(\n", + " observables, observable_indices, coefficients\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.0\n" + ] + } + ], + "source": [ + "dm_result = densitytensor_to_expectation(params_to_densitytensor())\n", + "print(dm_result)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "qujax", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/qujax/__init__.py b/qujax/__init__.py index 5d75a16..62c6587 100644 --- a/qujax/__init__.py +++ b/qujax/__init__.py @@ -41,6 +41,8 @@ import qujax.typing +import qujax.experimental.statetensor + # pylint: disable=undefined-variable del version del statetensor diff --git a/qujax/densitytensor.py b/qujax/densitytensor.py index ef35102..29249ee 100644 --- a/qujax/densitytensor.py +++ b/qujax/densitytensor.py @@ -18,7 +18,7 @@ from qujax.typing import ( MixedCircuitFunction, KrausOp, - GateFunction, + ParameterizedGateFunction, GateParameterIndices, ) @@ -86,7 +86,7 @@ def kraus( def _to_kraus_operator_seq_funcs( kraus_op: KrausOp, param_inds: Optional[Union[GateParameterIndices, Sequence[GateParameterIndices]]], -) -> Tuple[Sequence[GateFunction], Sequence[jax.Array]]: +) -> Tuple[Sequence[ParameterizedGateFunction], Sequence[jax.Array]]: """ Ensures Kraus operators are a sequence of functions that map (possibly empty) parameters to tensors and that each element of param_inds_seq is a sequence of arrays that correspond to the diff --git a/qujax/experimental/statetensor.py b/qujax/experimental/statetensor.py new file mode 100644 index 0000000..db99f86 --- /dev/null +++ b/qujax/experimental/statetensor.py @@ -0,0 +1,392 @@ +from typing import Any, Callable, Sequence, Tuple, Optional, Mapping, Union + +# Backwards compatibility with Python <3.10 +from typing_extensions import TypeVarTuple, Unpack + +import jax +import jax.numpy as jnp +from jax.typing import ArrayLike + +from qujax import gates +from qujax.typing import Gate, GateFunction + +from qujax.statetensor import apply_gate + +PyTree = Any + +Operation = Union[ + Gate, + str, +] + + +def wrap_parameterised_tensor( + gate_func: Callable, qubit_inds: Sequence[int] +) -> Callable: + """ + Takes a callable representing a parameterised gate and wraps it in a function that takes + the returned jax.Array and applies it to the qubits specified by `qubit_inds`. + + Args: + gate_func: Callable representing parameterised gate. + qubit_inds: Indices gate is to be applied to. + + Returns: + Callable taking in gate parameters, input statetensor and input classical registers, + and returning updated statetensor after applying parameterized gate to specified qubits. + """ + + def unitary_op( + params: Tuple[jax.Array], + statetensor_in: jax.Array, + classical_registers_in: PyTree, + ): + gate_unitary = gate_func(*params[0]) + statetensor = apply_gate(statetensor_in, gate_unitary, qubit_inds) + + return statetensor, classical_registers_in + + return unitary_op + + +def _array_to_callable(arr: jax.Array) -> Callable[[], jax.Array]: + """ + Wraps array `arr` in a callable that takes no parameters and returns `arr`. + """ + + def _no_param_tensor(): + return arr + + return _no_param_tensor + + +def _to_gate_func( + gate: Gate, + tensor_dict: Mapping[str, Union[Callable, jax.Array]], +) -> GateFunction: + """ + Converts a gate specification to a callable function that takes the gate parameters and returns + the corresponding unitary. + + Args: + gate: Gate specification. Can be either a string, a callable or a jax.Array. + + Returns: + Callable taking gate parameters and returning + """ + + if isinstance(gate, str): + gate = tensor_dict[gate] + if isinstance(gate, jax.Array): + gate = _array_to_callable(gate) + if callable(gate): + return gate + else: + raise TypeError( + f"Unsupported gate type - gate must be either a string in qujax.gates, a jax.Array or " + f"callable: {type(gate)}" + ) + + +def parse_op( + op: Operation, + params: Sequence[Any], + gate_dict: Mapping[str, Union[Callable, jax.Array]], + op_dict: Mapping[str, Callable], +) -> Callable: + """ + Parses operation specified by `op`, applying relevant metaparameters and returning a callable + retpresenting the operation to be applied to the circuit. + + Args: + op: Operation specification. Can be: + - A string, in which case we first check whether it is a gate by looking it up in + `tensor_dict` and then check whether it is a more general operation by looking it up + in `op_dict`. + - A jax.Array, which we assume to represent a gate. + - A callable, which we assume to represent a parameterized gate. + params: Operator metaparameters. For gates, these are the qubit indices the gate is to + be applied to. + tensor_dict: Dictionary mapping strings to gates. + op_dict: Dictionary mapping strings to callables that take operation metaparameters and + return a function representing the operation to be applied to the circuit. + + Returns: + A callable encoding the operation to be applied to the circuit. + """ + # Gates + if ( + (isinstance(op, str) and op in gate_dict) + or isinstance(op, jax.Array) + or callable(op) + ): + op = _to_gate_func(op, gate_dict) + return wrap_parameterised_tensor(op, params) + + if isinstance(op, str) and op in op_dict: + return op_dict[op](*params) + + if isinstance(op, str): + raise ValueError(f"String {op} not a known gate or operation") + else: + raise TypeError( + f"Invalid specification for `op`, got type {type(op)} with value {op}" + ) + + +def get_default_gates() -> dict: + """ + Returns dictionary of default gates supported by qujax. + """ + return { + k: v for k, v in gates.__dict__.items() if not k.startswith(("_", "jax", "jnp")) + } + + +def _gate_func_to_unitary( + gate_func: GateFunction, + n_qubits: int, + params: jax.Array, +) -> jax.Array: + """ + Compute tensor representing parameterised unitary for specific parameters. + + Args: + gate_func: Function that maps a (possibly empty) parameter array to a unitary tensor + n_qubts: Number of qubits unitary acts on + params: Parameter vector + + Returns: + Array containing gate unitary in tensor form. + """ + gate_unitary = gate_func(*params) + gate_unitary = gate_unitary.reshape( + (2,) * (2 * n_qubits) + ) # Ensure gate is in tensor form + return gate_unitary + + +Op = Callable[ + [Tuple[jax.Array, ...], jax.Array, jax.Array], Tuple[jax.Array, jax.Array] +] +OpSpecArgs = TypeVarTuple("OpSpecArgs") +OpSpec = Callable[[Unpack[OpSpecArgs]], Op] + + +def get_default_operations( + gate_dict: Mapping[str, Union[Callable, jax.Array]] +) -> Mapping[str, OpSpec]: + """ + Returns dictionary of default operations supported by qujax. Each operation is a function + that takes a set of metaparemeters and returns another function. The returned function + must have three arguments: `op_params`, `statetensor_in` and `classical_registers_in`. + `op_params` holds parameters that are passed when the circuit is executed, while + `statetensor_in` and `classical_registers_in` correspond to the statetensor + and classical registers, respectively, being modified by the circuit. + + Parameters: + `gate_dict`: Dictionary encoding quantum gates that the circuit can use. This + dictionary maps strings to a callable in the case of parameterized gates or to a + jax.Array in the case of unparameterized gates. + """ + op_dict: dict[str, OpSpec] = dict() + + def generic_op(f: Op) -> Op: + """ + Generic operation to be applied to the circuit, passed as a metaparameter `f`. + """ + return f + + def conditional_gate(gates: Sequence[Gate], qubit_inds: Sequence[int]) -> Op: + """ + Operation applying one of the gates in `gates` according to an index passed as a + circuit parameter. + + Args: + gates: gates from which one is selected to be applied + qubit_indices: indices of qubits the selected gate is to be applied to + """ + gate_funcs = [_to_gate_func(g, gate_dict) for g in gates] + + def apply_conditional_gate( + op_params: Union[Tuple[jax.Array], Tuple[jax.Array, jax.Array]], + statetensor_in: jax.Array, + classical_registers_in: jax.Array, + ) -> Tuple[jax.Array, jax.Array]: + """ + Applies a gate specified by an index passed in `op_params` to a statetensor. + + Args: + op_params: gates from which one is selected to be applied + statetensor_in: indices of qubits the selected gate is to be applied to + classical_registers_in: indices of qubits the selected gate is to be applied to + """ + if len(op_params) == 1: + ind, gate_params = op_params[0], jnp.empty((len(gates), 0)) + elif len(op_params) == 2: + ind, gate_params = op_params[0], jnp.array(op_params[1]) + else: + raise ValueError("Invalid number of parameters for ConditionalGate") + + unitaries = jnp.stack( + [ + _gate_func_to_unitary( + gate_funcs[i], len(qubit_inds), gate_params[i] + ) + for i in range(len(gate_funcs)) + ] + ) + + chosen_unitary = unitaries[ind] + + statevector = apply_gate(statetensor_in, chosen_unitary, qubit_inds) + return statevector, classical_registers_in + + return apply_conditional_gate + + op_dict["Generic"] = generic_op + op_dict["ConditionalGate"] = conditional_gate + + return op_dict + + +ParamInds = Optional[ + Union[ + int, + Sequence[int], + Sequence[Sequence[int]], + Mapping[str, int], + Mapping[str, Sequence[int]], + ] +] + + +def get_params( + param_inds: ParamInds, + params: Union[Mapping[str, ArrayLike], ArrayLike], +) -> Tuple[Any, ...]: + """ + Extracts parameters from `params` using indices specified by `param_inds`. + + Args: + param_inds: Indices of parameters. Can be + - None (results in an empty jax.Array) + - an integer, when `params` is an indexable array + - a dictionary, when `params` is also a dictionary + - nested list or tuples of the above + params: Parameters from which a subset is picked. Can be either an array or a dictionary + of arrays + Returns: + Tuple of indexed parameters respeciting the structure of nested lists/tuples of param_inds. + + """ + op_params: Tuple[Any, ...] + if param_inds is None: + op_params = (jnp.array([]),) + elif isinstance(param_inds, int) and isinstance(params, jax.Array): + op_params = (params[param_inds],) + elif isinstance(param_inds, dict) and isinstance(params, dict): + op_params = tuple( + jnp.take(params[k], jnp.array(param_inds[k])) for k in param_inds + ) + elif isinstance(param_inds, (list, tuple)): + if len(param_inds): + if all(isinstance(x, int) for x in param_inds): + op_params = (jnp.take(params, jnp.array(param_inds)),) + else: + op_params = tuple(get_params(p, params) for p in param_inds) + else: + op_params = (jnp.array([]),) + else: + raise TypeError( + f"Invalid specification for parameters: {type(param_inds)=} {type(params)=}." + ) + return op_params + + +def get_params_to_statetensor_func( + op_seq: Sequence[Operation], + op_metaparams_seq: Sequence[Sequence[Any]], + param_pos_seq: Sequence[ParamInds], + op_dict: Optional[Mapping[str, OpSpec]] = None, + gate_dict: Optional[Mapping[str, Union[jax.Array, GateFunction]]] = None, +): + """ + Creates a function that maps circuit parameters to a statetensor. + + Args: + op_seq: Sequence of operations to be executed. + Can be either + - a string specifying a gate in `gate_dict` + - a jax.Array specifying a gate + - a function returning a jax.Array specifying a parameterized gate. + - a string specifying an operation in `op_dict` + op_params_seq: Sequence of operation meta-parameters. Each element corresponds to one + operation in `op_seq`. For gates, this will be the qubit indices the gate is applied to. + param_pos_seq: Sequence of indices specifying the positions of the parameters each gate + or operation takes. + Note that these are parameters of the circuit, and are distinct from the meta-parameters + fixed in `op_params_seq`. + op_dict: Dictionary mapping strings to operations. Each operation is a function + taking metaparameters (which are specified in `op_params_seq`) and returning another + function. This returned function encodes the operation, and takes an array of + parameters, a statetensor and classical registers, and returns the updated statetensor + and classical registers after the operation is applied. + gate_dict: Dictionary mapping strings to gates. Each gate is either a jax.Array or a + function taking a number of parameters and returning a jax.Array. + Defaults to qujax's dictionary of gates. + Returns: + Function that takes a number of parameters, an input statetensor and an input set of + classical registers, and returns the updated statetensor and classical registers + after the specified gates and operations are applied. + """ + if gate_dict is None: + gate_dict = get_default_gates() + if op_dict is None: + op_dict = get_default_operations(gate_dict) + + repeated_ops = set(gate_dict.keys()) & set(op_dict.keys()) + if repeated_ops: + raise ValueError( + f"Operator list and gate list have repeated operation(s): {repeated_ops}" + ) + + parsed_op_seq = [ + parse_op(op, params, gate_dict, op_dict) + for op, params in zip(op_seq, op_metaparams_seq) + ] + + def params_to_statetensor_func( + params: Union[Mapping[str, ArrayLike], ArrayLike], + statetensor_in: jax.Array, + classical_registers_in: Optional[jax.Array] = None, + ) -> Tuple[jax.Array, PyTree]: + """ + Applies parameterised circuit to the quantum state represented by `statetensor_in`. + + Args: + params: Parameters to be passed to the circuit + statetensor_in: Input state in tensor form. + classical_registers_in: Classical registers that can store intermediate results + (e.g. measurements), possibly to later reuse them + Returns: + Resulting quantum state and classical registers after applying the circuit. + + """ + statetensor = statetensor_in + classical_registers = classical_registers_in + for ( + op, + param_pos, + ) in zip( + parsed_op_seq, + param_pos_seq, + ): + op_params = get_params(param_pos, params) + statetensor, classical_registers = op( + op_params, statetensor, classical_registers_in + ) + + return statetensor, classical_registers + + return params_to_statetensor_func diff --git a/qujax/statetensor.py b/qujax/statetensor.py index ee95299..bd0eecf 100644 --- a/qujax/statetensor.py +++ b/qujax/statetensor.py @@ -12,7 +12,12 @@ from qujax import gates from qujax.utils import _arrayify_inds, check_circuit -from qujax.typing import Gate, PureCircuitFunction, GateFunction, GateParameterIndices +from qujax.typing import ( + Gate, + PureCircuitFunction, + ParameterizedGateFunction, + GateParameterIndices, +) def apply_gate( @@ -41,7 +46,7 @@ def apply_gate( def _to_gate_func( gate: Gate, -) -> GateFunction: +) -> ParameterizedGateFunction: """ Ensures a gate_seq element is a function that map (possibly empty) parameters to a unitary tensor. @@ -74,7 +79,7 @@ def _array_to_callable(arr: jax.Array) -> Callable[[], jax.Array]: def _gate_func_to_unitary( - gate_func: GateFunction, + gate_func: ParameterizedGateFunction, qubit_inds: Sequence[int], param_inds: jax.Array, params: jax.Array, diff --git a/qujax/typing.py b/qujax/typing.py index 3d78369..582ccbf 100644 --- a/qujax/typing.py +++ b/qujax/typing.py @@ -1,4 +1,4 @@ -from typing import Union, Optional, Protocol, Callable, Iterable, Sequence +from typing import Union, Optional, Protocol, Callable, Sequence # Backwards compatibility with Python <3.10 from typing_extensions import TypeVarTuple, Unpack @@ -34,7 +34,9 @@ def __call__(self, densitytensor_in: Optional[jax.Array] = None) -> jax.Array: GateArgs = TypeVarTuple("GateArgs") # Function that takes arbitrary nr. of parameters and returns an array representing the gate # Currently Python does not allow us to restrict the type of the arguments using a TypeVarTuple -GateFunction = Callable[[Unpack[GateArgs]], jax.Array] +ParameterizedGateFunction = Callable[[Unpack[GateArgs]], jax.Array] +UnparameterizedGateFunction = Callable[[], jax.Array] +GateFunction = Union[ParameterizedGateFunction, UnparameterizedGateFunction] GateParameterIndices = Optional[Sequence[int]] PureCircuitFunction = Union[PureUnparameterizedCircuit, PureParameterizedCircuit] @@ -42,4 +44,4 @@ def __call__(self, densitytensor_in: Optional[jax.Array] = None) -> jax.Array: Gate = Union[str, jax.Array, GateFunction] -KrausOp = Union[Gate, Iterable[Gate]] +KrausOp = Union[Gate, Sequence[Gate]] diff --git a/tests/test_experimental.py b/tests/test_experimental.py new file mode 100644 index 0000000..d756e88 --- /dev/null +++ b/tests/test_experimental.py @@ -0,0 +1,96 @@ +import jax +import jax.numpy as jnp + +import qujax +from qujax import all_zeros_statetensor, apply_gate +from qujax.experimental.statetensor import get_params_to_statetensor_func + + +def test_get_params_to_statetensor_func(): + ops = ["H", "H", "H", "CX", "Rz", "CY"] + op_params = [[0], [1], [2], [0, 1], [1], [1, 2]] + param_inds = [[], [], [], None, [0], []] + + param_to_st = get_params_to_statetensor_func(ops, op_params, param_inds) + param_to_st = jax.jit(param_to_st) + param = jnp.array(0.1) + st_in = all_zeros_statetensor(3) + st, _ = param_to_st(param, st_in) + + true_sv = jnp.array( + [ + 0.34920055 - 0.05530793j, + 0.34920055 - 0.05530793j, + 0.05530793 - 0.34920055j, + -0.05530793 + 0.34920055j, + 0.34920055 - 0.05530793j, + 0.34920055 - 0.05530793j, + 0.05530793 - 0.34920055j, + -0.05530793 + 0.34920055j, + ], + dtype="complex64", + ) + + assert st.size == true_sv.size + assert jnp.allclose(st.flatten(), true_sv) + + +def test_stochasticity(): + ops = ["ConditionalGate"] + op_params = [[["X", "Y", "Z"], [0]]] + param_inds = [[{"op_ind": 0}]] + + st_in = all_zeros_statetensor(1) + X_apply = apply_gate(st_in, qujax.gates.X, [0]) + Y_apply = apply_gate(st_in, qujax.gates.Y, [0]) + Z_apply = apply_gate(st_in, qujax.gates.Z, [0]) + + param_to_st = get_params_to_statetensor_func(ops, op_params, param_inds) + param_to_st = jax.jit(param_to_st) + + st_in = all_zeros_statetensor(1) + + st_X, _ = param_to_st({"op_ind": 0}, st_in) + st_Y, _ = param_to_st({"op_ind": 1}, st_in) + st_Z, _ = param_to_st({"op_ind": 2}, st_in) + + assert jnp.allclose(X_apply, st_X) + assert jnp.allclose(Y_apply, st_Y) + assert jnp.allclose(Z_apply, st_Z) + + +def test_parameterised_stochasticity(): + ops = ["ConditionalGate"] + op_params = [[["Rx", "Ry", "Rz"], [0]]] + param_inds = [[{"op_ind": 0}, [{"angles": 0}, {"angles": 1}, {"angles": 2}]]] + + st_in = all_zeros_statetensor(1) + params = jnp.array([0.1, 0.2, 0.3]) + + CX_apply = apply_gate(st_in, qujax.gates.Rx(params[0].item()), [0]) + CY_apply = apply_gate(st_in, qujax.gates.Ry(params[1].item()), [0]) + CZ_apply = apply_gate(st_in, qujax.gates.Rz(params[2].item()), [0]) + + param_to_st = get_params_to_statetensor_func(ops, op_params, param_inds) + + st_in = all_zeros_statetensor(1) + + st_CX, _ = param_to_st({"angles": params, "op_ind": 0}, st_in) + st_CY, _ = param_to_st({"angles": params, "op_ind": 1}, st_in) + st_CZ, _ = param_to_st({"angles": params, "op_ind": 2}, st_in) + + assert jnp.allclose(CX_apply, st_CX) + assert jnp.allclose(CY_apply, st_CY) + assert jnp.allclose(CZ_apply, st_CZ) + + batched_op_inds = jnp.array([[0], [1], [2]]) + + batched_param_to_st = jax.vmap( + param_to_st, in_axes=({"angles": None, "op_ind": 0}, None) + ) + + batched_st, _ = batched_param_to_st( + {"angles": params, "op_ind": batched_op_inds}, st_in + ) + + assert jnp.allclose(batched_st, jnp.stack([st_CX, st_CY, st_CZ]))