QFit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting functions with a quantum circuit\n",
    "___________________________________\n",
    "\n",
    "_This notebook is part of the Q-Hack 2021 presentation (you can find the [slides here](https://albacl.github.io/files/QHack2021.pdf) and the [recording here](https://www.twitch.tv/videos/920118091))._\n",
    "\n",
    "In this notebook, I will present simple circuits to encode the data into a quantum circuit utilizing the angles in the rotational gates. With that, we can check if the quantum circuit can approximate a given function. The explicit circuits presented are for only one qubit, but its generalization to multiple-qubits follows.\n",
    "\n",
    "This notebook is based on the references:\n",
    "\n",
    "- _One qubit as a Universal Approximant_, A. Pérez-Salinas, D. López-Núñez, A. García-Sáez, P. Forn-Díaz, J. I. Latorre, [arXiv:2102.04032 [quant-ph]](https://arxiv.org/abs/2102.04032) (2021).\n",
    "- _The effect of data encoding on the expressive power of variational quantum machine learning models_, M. Schuld, R. Sweke, J. J. Meyer, [arXiv:2008.08605 [quant-ph]](https://arxiv.org/abs/2008.08605) (2020).\n",
    "- _Data re-uploading for a universal quantum classifier_, A. Pérez-Salinas, A. Cervera-Lierta, E. Gil-Fuster, J. I. Latorre, [Quantum 4, 226 (2020)](https://quantum-journal.org/papers/q-2020-02-06-226/).\n",
    "\n",
    "You will need to install [Tequila](https://github.com/aspuru-guzik-group/tequila) package to run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tequila as tq\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum Circuit\n",
    "\n",
    "We will encode the data into the rotational gates of a quantum circuit.\n",
    "As an example, for one qubit, we select the following encoding:\n",
    "$$ U\\left(\\vec{\\theta},x\\right) = R_{y}\\left(\\theta^{(1)}\\right)R_{z}\\left(\\theta^{(2)}x + \\theta^{(3)}\\right), $$\n",
    "where $x$ is a given data point (in this case, we consider 1D data), $\\vec{\\theta}=\\left(\\theta^{(1)}, \\theta^{(2)}, \\theta^{(3)}\\right)$ and $R_{k}\\left(\\theta\\right)=e^{i\\frac{\\theta}{2}\\sigma_{k}}$ with $k=x,y,z$.\n",
    "\n",
    "The strategy will be the following: the total quantum circuit will be a concatenation of $U$ gates that feed the circuit with the data and some parameters to be optimized, i.e.\n",
    "$$ U_{feed}\\left(\\vec{\\theta}_{1},\\cdots,\\vec{\\theta}_{L};x\\right) = U\\left(\\vec{\\theta}_{L},x\\right)U\\left(\\vec{\\theta}_{L-1},x\\right)\\cdots U\\left(\\vec{\\theta}_{1},x\\right). $$\n",
    "We need to establish the number of _layers_ $L$ of our circuit and we will have to find the optimal values of the $3L$ parameters. We expect that the larger the number of layers, the more sophisticated will be the encoding and the higher the accuracy of the function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def qcircuit(x, layers, param):\n",
    "    qc = tq.QCircuit()\n",
    "    p = 0 # all parameters are stored in a single list\n",
    "    for l in range(layers): # add layers to the circuit\n",
    "        qc += tq.gates.Ry(param[p],0) + tq.gates.Rz(x*param[p+1]+param[p+2],0) # feed the troll!\n",
    "        p = p + 3\n",
    "    return qc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning the function\n",
    "\n",
    "Once we have the data encoded into the quantum circuit (a.k.a. _feature map_), we have to design a cost function that will allow us to optimize the circuit parameters.\n",
    "\n",
    "Since we are dealing with only one qubit, a simple approach is to measure it on the computational basis and depending on the expectation value of $Z$ operator, we compare it with the value of the target function to learn. This function has to be rescaled between -1 and 1 accordingly. Thus, we construct the cost function by taking a set of data training points $\\left(x_{1},\\cdots,x_{M}\\right)$, running out $U_feed$ circuit with each of them and computing the expectation value of $Z$ with the resultant state:\n",
    "$$\\mathcal{L}oss\\left(\\vec{\\theta}_{1},\\cdots,\\vec{\\theta}_{L}\\right)=\\sum_{i=1}^{M}\\left(|\\langle 0|U_{feed}^{\\dagger}(x_{i})ZU_{feed}(x_{i})|0\\rangle|^2 - f(x_{i})\\right)^2, $$\n",
    "were I ommited the $\\theta$ dependency in $U_{feed}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function\n",
    "def f(x):\n",
    "    func = np.sin(10*x)/(10*x)\n",
    "    return func\n",
    "\n",
    "# Cost function\n",
    "def loss(xval,layers,param):\n",
    "    obj = 0\n",
    "    for i in range(len(xval)):\n",
    "        qc = qcircuit(xval[i],layers, param)\n",
    "        expval = tq.ExpectationValue(qc, H=tq.paulis.Z(0))\n",
    "        obj += (expval-f(xval[i]))**2\n",
    "    return obj/len(xval)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are ready for the circuit optimization! To compare how the accuracy changes with the number of layers, we will optimize for up to `laymax` layers. We have also to create the training data set (in this example, we will take random points between -1 and 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1  layers\n",
      "Loss function =  0.13260757536486664\n",
      "2  layers\n",
      "Loss function =  0.04352773925273443\n",
      "3  layers\n",
      "Loss function =  0.00035727648165960176\n",
      "4  layers\n",
      "Loss function =  0.0001931700809334706\n"
     ]
    }
   ],
   "source": [
    "training_set = 25 # number of training points\n",
    "laymax = 4\n",
    "\n",
    "# generate the training data\n",
    "xval=[]\n",
    "for i in range(training_set):\n",
    "    rand = random.uniform(-1,1)\n",
    "    xval.append(rand)\n",
    "\n",
    "# Optimization details    \n",
    "grad = '2-point' # numerical gradient (= None: analytical gradient)\n",
    "mthd = 'BFGS' # minimization method\n",
    "backend = 'qulacs' # you can choose another, Qiskit, pyquil, cirq,... check the documentation\n",
    "mthd_opt = {'eps':1.e-4} # method options (that's the stepsize for the gradients)\n",
    "\n",
    "train_param = [[]]*(laymax+1)\n",
    "for l in range(1,laymax+1):\n",
    "    \n",
    "    # Generate the needed variables\n",
    "    param = [tq.Variable(name='th{}'.format(i)) for i in range(0,3*l)]\n",
    "    # Initialize them at random\n",
    "    param0 = {key : random.uniform(0, 2*np.pi) for key in param}\n",
    "    \n",
    "    # Optimize!\n",
    "    train = tq.minimize(objective=loss(xval,l,param), initial_values = param0, method = mthd, \n",
    "                        gradient = grad, method_options = mthd_opt, backend = backend, \n",
    "                        silent=True) # silent=False to see the optimization\n",
    "  \n",
    "    print(\"{}\".format(l),\" layers\")\n",
    "    print(\"Loss function = \", train.energy)\n",
    "    \n",
    "    train_param[l] = train.angles # saving the parameters \n",
    "    # You can plot the history plot of the optimization! both for the energy and angles. Check the documentation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test and results\n",
    "\n",
    "Once we have the optimized parameters, we can proceed to test the results by running new circuits with other data values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABao0lEQVR4nO2deXgURfrHPzWTTBLCHSaAAhIFBAIhHIGEQwIoggdhQVy88EDwxHM9EK/9oeKKrojiCggIgriiKKgcrsTIFTAcQe5bBBVnDHJKSDLz/v7oGRjC5J4jR32ep5+Z7q6u+nZP97xd9Va9pUQEjUaj0WgKwhRsARqNRqMp32hDodFoNJpC0YZCo9FoNIWiDYVGo9FoCkUbCo1Go9EUSkiwBfiDevXqSdOmTYMtQ6PRaCoM69ev/0NErN72VUpD0bRpU9atWxdsGRqNRlNhUEodKGifbnrSaDQaTaFoQ6HRaDSaQtGGQqPRaDSFUil9FBqNpuqRm5vLoUOHyM7ODraUck14eDiNGjUiNDS02MdoQ6HRaCoFhw4dokaNGjRt2hSlVLDllEtEhKysLA4dOkRMTEyxjwtq05NSarpSyqaU2lLAfqWUmqiU2qOU+lEp1SHQGjUaTcUgOzubqKgobSQKQSlFVFRUiWtdwfZRfAD0K2R/f6C5axkJ/CcAmjQaTQVFG4miKc01CmrTk4gsV0o1LSRJCjBLjFjoa5RStZVSDUXkt8Ao1Gh8w4kTJ1gwfTobli+nSZs2XHXjjcTGxgZblkZTLIJdoyiKi4GDHuuHXNsuQCk1Uim1Tim1zm63B0ScRlMUTqeTV558kiYNGnDbI4/w5vz5PPp//0ebNm1ISUlh7969wZao8RF33XUX0dHRtGnTpsA0L774Iq+//noAVfmG8m4ovNWRvM60JCJTRKSTiHSyWr2OQtdoAsqZM2e44ZprGDN+PBcDo5s0Yerll/Nm88sYdU0yaWlpdO7cmVWrVgVbqsYH3HHHHSxZsiSoGvLy8vySb3k3FIeAxh7rjYBfg6RFoyk2v2/YwDVduvD50qX8PdrK6EsuoU31SKqZTTSMCOW5u9uzfv16oqKiuPLKK8nIyAi25KpJejqMG2d8lpErrriCunXrFjv91KlTSUhIoF27dgwePJi//vqLEydOEBMTQ25uLgDHjx+nadOm5ObmsnfvXvr160fHjh3p0aMHO3bsAAwD9dhjj9GrVy+eeuopvv/+e+Lj44mPj6d9+/acOHGizOdW3g3FQmCYq/dTInBM+yc05R17Zib3XnMNqZs2cftFDRlYPwqTCUxmoVn7P+l922Gs3YfQrFkzVq1aRYMGDUhJSeHXX/U7UEBJT4c+feC554xPHxiLkjBo0CAyMjLYtGkTrVq1Ytq0adSoUYPk5GS+/vprAD7++GMGDx5MaGgoI0eO5O2332b9+vW8/vrr3H///Wfz2rVrF99++y1vvPEGr7/+OpMmTSIzM5MVK1YQERFRZq1BdWYrpeYCyUA9pdQh4AUgFEBE3gMWAdcAe4C/gDuDo1SjKT4LZs7ki99/p1ft2vStU4dL448TWfMM0U3PYL3iZogZBtYkAKxWKwsXLiQpKYnbbruNb7/9VvfcCRRpaZCTAw6H8ZmWBklJASt+y5YtPPvssxw9epSTJ09y9dVXA3D33Xfz2muvMXDgQGbMmMHUqVM5efIkq1evZsiQIWePP3PmzNnvQ4YMwWw2A9CtWzcee+wxbrnlFgYNGkSjRo3KrDXYvZ5uKmK/AA8ESI5GU2Z+Sk/n+RnTaGixcHvD+pjNTmKG3I31UgtEJ581EJ60bduWN954g3vvvZfpY29g+H3/8JpO42OSk8FiMYyExWKsB5A77riDL774gnbt2vHBBx+QlpYGGH/0P/30E99//z0Oh4M2bdpw/PhxateuTWZmpte8IiMjz35/+umnufbaa1m0aBGJiYl8++23tGzZskxay3vTk0ZTYbBnZjJqYAq/HTvByEYNaNXhGL1vPWgYidjRhf75jxjYhitamXj8X/Oxfd4b7IFtBqmSJCXBsmUwdqzxGcDaBBhdphs2bEhubi5z5sw5b9+wYcO46aabuPNOoxGlZs2axMTEMG/ePMAYYb1p0yav+e7du5e2bdvy1FNP0alTp7O+jLKgDYVG4yMy509jif0PetauzeURkUTWysPaxGnUJIrA9MdyJt8FJ7Ph5flnwJbmd70aDOMwerRPjMRNN91EUlISO3fupFGjRkybNq3Q9GPHjqVLly5cddVVF7zx33LLLfz555/cdNO5Rpc5c+Ywbdo02rVrR2xsLAsWLPCa74QJE2jTpg3t2rUjIiKC/v37l/nclNG6U7no1KmT6ImLNIHEnpnJnTcP4Jsdh/h388uwhofQ+6H6WAeMK14zkj0d+8xruG9GNgt3ZrNr7Wc07TjI/8IrEdu3b6dVq1bBluETPv30UxYsWMCHH37ol/y9XSul1HoR6eQtvQ4KqNGUEXtmJh/ePJRF2w9ybb26dO58iph2p7EOmF5sX4P9lwhSP7qEno6/WOjcw+j/m8ncBdpQVEVGjRrF4sWLWbRoUbClnEUbCo2mjOz/7xS+/O0wIUpxXVQUkTHNsd7+bIkc0raMDBy5DqLMofSuW5dPv/qK8YcO+aTHiqZi8fbbbwdbwgVoH4VGUwbsmZmsX7CMlceO07N2LWpZzES3v7zEvZaiExIwWywos5lrGzTACbz11lv+Ea3RlBBtKDSaMmBL+4yl9iwcIvSPiuLSdiewdh9S9IH5sMbH03vaNOJGjWLo7NkMGTKEyZMnc+zYMT+o1mhKhjYUGk1Z+CuT1D+P0qFGdS6uFkpM/56lHgNhjY8ndsQIrPHxPP7445w4cYLZs2f7WLBGU3K0odBoSok9bTZTJm/ghMPBlXXr0OHqP7Be9aBP8k5ISKBDhw5MnjyZytgzUVOx0IZCoykltpWL+TbrKNbQUNpGRpIT3tmnI6pHjhzJ5s2bWTv3fj0Ar5xz8OBBevXqRatWrYiNjS3Qv6TDjGs0VYyTTduz7dRf9K5Tm5AQiO5daESaEnNz38uoHg5TpkyG1D7aWJRjQkJCeOONN9i+fTtr1qxh0qRJbNu2LeA6qmqYcY2m3PL1ruOYTIp77+hO77fGYE2+1af51/grgxu7KOatFf46rUdr+wNfRRlv2LAhHTp0AKBGjRq0atWKX375pdBjKlKYcUSk0i0dO3YUjcafOJ1OiYmJkb59+/qvENtq+e45iwAy96FQEdtq/5VVCdi2bVuJ0q9eLRIRIWI2G5+rfXR59+/fL40bN5Zjx45dsO+FF16Q8ePHi4jIH3/8cXb7mDFjZOLEiSIicscdd8jnn38uIiKTJ0+Wxx57TEREevfuLbt27RIRkTVr1kivXr1EROT222+Xa6+9VvLy8kRE5LrrrpOVK1eKiMiJEyckNzf3Ah3erhWwTgr4T9U1Co2mFKxevZr9+/dz662+rUWchzWJKx5IpVH9Wny4rZOOKOtjvEUZLysnT55k8ODBTJgwgZo1axaadsuWLfTo0YO2bdsyZ84ctm7dChhhxmfMmAHAjBkzuPPOO88LMx4fH88999zDb7+dm5rHW5jxiRMncvToUUJCyj6uWhsKjaYUTH3teSJCzXSvd6boxGXAVL8bt9xxL0u/+wGbzebXsqoa7ijjZrNvoozn5uYyePDgs/NAFMUdd9zBO++8w+bNm3nhhRfIzs4GvIcZdzqdZ8OMu5ft27efzSt/mPH333+f06dPk5iYqKPHajTB4LdlM5n/9ffER0Tyw9NvYE/z31gHe2YmXcLCcDgczJ8/32/lVEV8GWVcRBg+fDitWrXiscceK9YxOsy4RlOJWTxnFiccDrrUrInTobCtXOyXcuyZmaQOH85fn33GRWFhfORqjtD4Dl9FGV+1ahUffvghqampZx3JRQX1q0hhxnVQQI2mhKzIUoQpRbsakZjMQnT3sj+I3rBlZODIyUGJ0LlmTb5ctw673Y7VavVLeZrS071792INjHzxxRfPfr/vvvu47777vKZbuXIlN9xwA7Vr1z67LSYmhiVLllyQ9oMPPjhv3R9BBbWh0GhKgMPhYPHaLVx9RQIJVzYjunt/n3eLdeMOFOjMzSUxKoov7HYWLFjA3Xff7ZfyNOUDHWZco6ngrFq1it9//51bJk4k9sYb/VqWO1CgLSODPp06MXXIED799FNtKCo55THMuDYUGk0JWLBgARaLxSftvsXBGh+PNT4egJSUFN555x1OnjxJ9erVA1K+RgPama3RFBsRYeHChfTu3ZsaNWoEvPxrr72WnJwcls2fBFvH6ZAemoChDYVGU0x27tzJnj17uP7664NSfvfu3alRvRpfz3wGNo2BZb20sdAEhKAaCqVUP6XUTqXUHqXU017211JKfamU2qSU2qqUujMYOjUagC+//BIgaIbCYrHQN6EBizY6jR42zjOwf1ZQtGiqFkEzFEopMzAJ6A+0Bm5SSrXOl+wBYJuItAOSgTeUUpaACtVoXCxcuJD4+HgaN24cNA3XdmvCL3/CpgNBk6DxQnZ2Np07dz47xuGFF17wmk6HGS85nYE9IrJPRHKAj4GUfGkEqKGUUkB14Ajgnzi6Gk0hHDt2jPT0dK655pqg6uj/98cB+DoTMFkgZlhQ9WgMwsLCSE1NZdOmTWRmZrJkyRLWrFkTcB2VMcz4xcBBj/VDrm2evAO0An4FNgMPi4gzMPI0mnOkpqbicDi4+uqrg6qjQZvr6BTfkq93NYE+aTpQYFmxp/ukY4BS6mxPtNzcXHJzczHebwumIoUZD6ah8HYV8w9tvBrIBC4C4oF3lFJeQzIqpUYqpdYppdbZ7XZf6tRUdezpfPPxq1SPrEZSWWM9+IBrU/7Omh8P8sfOr7UzuyzY040JoTY955OJoRwOB/Hx8URHR3PVVVfRpUuXQtMPGjSIjIwMNm3aRKtWrZg2bRo1atQgOTmZr7/+GoCPP/6YwYMHExoaysiRI3n77bdZv349r7/+Ovfff//ZvHbt2sW3337LG2+8weuvv86kSZPIzMxkxYoVRERElOm8ILiG4hDg2djbCKPm4MmdwHxXuPQ9wH6gJV4QkSki0klEOukQBxqfYU+HZb34Ju0HerfMJvToumAr4tqujRERlnz0ip75rizY0sCRAzjAmVPmiaHMZjOZmZkcOnSIH374gS1bthSaXocZLx4ZQHOlVIzLQT0UWJgvzc9AHwClVH3gcmBfQFVqqjb7Z7F2A+yzQVIjc7noZdSx4WHq14KvMwUc2eVCU4UkOhnMFlBmw98TneyTbGvXrk1ycrLXuEye6DDjxUBE8oAHgaXAduATEdmqlLpXKXWvK9lYoKtSajOwDHhKRP4IjmJNVcS+5ziT59QGIHx7Y+x7jgdXELA3w0y7GpEs2ahwOgX2zdC1itJgTYLeyyBurPFZBn+P3W7n6NGjAJw+fZpvv/32goiw+dFhxouJiCwSkRYicpmIvOza9p6IvOf6/quI9BWRtiLSRkT8F/hfo/HC/s21+PHEKayhoVjNFvZvrhVUPbvnzSPjjQ+51FmLo9nCl0urg+Tp+bRLizUJYkeXuVPAb7/9Rq9evYiLiyMhIYGrrrqK6667rtBjdJhxjaaSkBdaj62n/qJbrVooZYLw+kHVc/CbbwBo42pqWLw+hJRrTGCJCqasKk9cXBwbN24sMp0OM67RVEKymjYl2+mkbfXqmCwWYlLyD/UJLI379uXw6tXUCQ3lIouF7dmnQByw4RGo3VZ3l60E6DDjGk0F44f9+zGbzfz98ce5rGfPs5Fcg0XzIUMAo2bRq66TLzJ2kJMHFuXqtaMNRYWnPIYZ10EBNZpCWLp0KYmJiSQ+/HDQjYSb5kOG0HvqVAbdNYxTZyBjn8mnvXY0mvxoQ6HRFMCxY8dYv349ffr0KTJtejqMGwdTphif6QHohJSccg9KKZb90avMvXY0msLQTU8aTQGsXLkSp9NJr169LtiXng5paRAVBYsXw5dfgtMJImAyQWgo3HknDHOFYkpLg+Rk8NXAbntmJraMDOJatmTZ5lye10ZC40e0odBoCiAtLQ2LxXJBKIb0dOjTB86cMYxDfpxOY9/kyTBjhmE8HA4ICTlnPMpiMOyZmaQOH44jJ4cmWVks2b2bv/76i2rVqpU+U42mEHTTk0ZTAGlpaSQmJl4QKyctDXJyvBsJAHcsOBEjXW6uYSjcxqNPn7I1TdkyMnC4BLSOiCA3L4+VK1eWPkONz3A4HLRv377AMRQ6zLhGU4k4fvw4GzZsIDk5mSlT4OqrDf8DGE1IFovRxOSJ2QxPPgn33GPsN5uNz9DQC41HWlrptUUnJGByxe9pERGBWSm+nDu39BlqfMZbb71Fq1atgla+v8KM66YnjcYLbv/EiRM9ueceY5trrBsjR8KyZed8FO5xVp5NSsOGnfNLAMyaBdOnGzULi+Xc9tJgjY/n0oED2TNvHuEmE80jIkj97rvSZ1iFcft6ohMSytyr7dChQ3z99deMGTOGf//730Wmnzp1KlOmTCEnJ4dmzZrx4Ycf4nA4iIuLY9euXYSGhnL8+HHi4uLYvXs3P//8Mw888AB2u51q1aoxdepUWrZsyR133EHdunXZuHEjHTp0YMCAATz88MOAEf58+fLlZZ7jXRsKjcYLbv/Ejz8mnrf9s88MQ5GUVLifIf/+pKTzjUdZndoxKSnsX7gQZ24usTVr8vnBgxw/fpyaNb1G4dd4wdPXY7ZY6D1tWpmMxSOPPMJrr71W7PkfBg0axIgRIwB49tlnmTZtGqNGjTobZnzgwIEXhBl/7733aN68OWvXruX+++8nNTUVOBdm3Gw2c/311zNp0iS6devGyZMnCQ8PL/U5udFNTxqNF77//nu6dOnCjTee7yAePLj0eSYlwejRvun5ZI2Pp/e0aVw2eDC9+vbF6XSyevXqsmdchfD09Thzc7FlZJQ6r6+++oro6Gg6duxY7GN0mHGNpgJz/Phx1q9fT3JyMiNHGg7ovn2Nz5Ej/VeueyxGSRzd+xcuJGLdOszA0k8/9Zu2ykh0QgJmiwVlNmMKDSU6IaHUea1atYqFCxfStGlThg4dSmpqKrfeemuhx1SkMOO66UmjceEeGxEWtgqHw0HPnj0Bwzj400C4y+7Tx3B0WyyGD6Somof7jTgcaBoRwYrly/0rspLhrpX5wkcxbtw4xo0bBxjNlq+//jqzZxce7Dp/mPGLLz43E7Q7zPhzzz0HnB9mfMiQIYgIP/74I+3atbsgX3eY8bZt25Kens6OHTuKDHleFLpGodFw7o/6uefgqafSCAkJDei0p+4utw5H8XtFeb4Rt6pRg80//XT2rVRTPKzx8cSOGBGU8CwVKcy4Esk/TXXFp1OnTrJuXfCnrNRUDNLT4cUX4dtv3WMjErnkklB++mlFQDWUtEYB53rtbDh1imGPPsry5cvp0aOH/wWXQ7Zv3x7Urqm+5NNPP2XBggV8+OGHfsnf27VSSq0XkU7e0uumJ02VZsoUeOAB401eBJQ6hcg6evZ8KqA6kpLOdbktSa8oa3w81vh4GmRlwaOPsmLFiiprKCoL5THMuG560lRZ0tPhwQchL89tJKBDh7WAg5tuCvyfrbdeUcV1cEdFRdGmTRuWaz9Fheftt99mz549tGjRIthSzqJrFJoqS1qaUZNwYzZDhw6r2LBBBdQ/URAlbY7q0aMHs2fPJi8vzyddIjUaN7pGoamyJCdDWNi5aK+TJsGBAytp06YNtWoFd25sKLmDu0ePHpw4cYJNmzYFQp6mCqENhaZKkZ4O991nLGC8pb/0Enz/PQwfsJL0Vd/TvWOz4Ip04Y4p5Y4ZVVTYD7dvYsWKwDnhNVUDbSg0VYb0dOjRA957z1jc00yMHg1JzdLZPPNKTpw6Q7fqX4E9ADMPFYHbwT12bPF6QTVq1IiYmBhtKDQ+RxsKTZXhtdfO90mcOePRnGNLY9X2HAC6N3cY80+XA0rq4O7RowcrVqygMnZ7r8p88cUXbNu2LWjla0OhqTLs2nX+ulIezTnRyazcbeLiutAkuvzOP52ebtSExowxPvMbi27dumG329m7d29wBGr8QpU2FEqpfkqpnUqpPUqppwtIk6yUylRKbVVKfR9ojZrKQXr6hYbiiSc83tStSazYU5u2FzXhjwZTy+3807NmGTUhEeNz1qzz93ft2hVABwgMErNnz6Zz585nA/etXbuWuLg4srOzOXXqFLGxsWzZsoWTJ0/Sp08fOnToQNu2bc8bZT1r1izi4uJo164dt912G6tXr2bhwoU88cQTxMfHB+UlIGh96JRSZmAScBVwCMhQSi0UkW0eaWoD7wL9RORnpVR0UMRqKjxpacafq5uBA+Ff/zq3vnHxYn6xZdHHbCH18TfpPa1NUMI6lJTDh89fb926NTVr1iQ9PZ1h7gm7qyCPPPIImZmZPs0zPj6eCRMmFLh/+/bt/Pe//2XVqlWEhoZy//33s3PnTgYMGMCzzz7L6dOnufXWW2nTpg15eXl8/vnn1KxZkz/++IPExEQGDBjAtm3bePnll1m1ahX16tXjyJEj1K1blwEDBnDddddxww03+PScikswO1t3BvaIyD4ApdTHQArgWb+6GZgvIj8DiIgt4Co1FRp3oL+oKKPnkHtMwpNPnp/u2/nzAbg8PPxsyOnyaCiGDYNp04zpVQG+/NIYXe4OWmgymUhMTNQ1iiCwbNky1q9fT4IrCu3p06eJjo7m+eefJyEhgfDwcCZOnAiAiPDMM8+wfPlyTCYTv/zyC7///jupqanccMMN1KtXD4C6desG7Xw8CaahuBg46LF+COiSL00LIFQplQbUAN4SkXyVbQOl1EhgJECTJk18LlZT8cg/YG3CBMjK8h4iY/vx44SbTDSJjCxzyGl/kpQEw4cbvbbAcM67u/q6jUVSUhJjx46t0hMZFfbm7y9EhNtvv/1sFFk3hw8f5uTJk+Tm5pKdnU1kZCRz5szBbrezfv16QkNDadq0KdnZ2YgIyj1vbjkimD4Kb1cjf1eNEKAjcC1wNfCcUsrruHYRmSIinUSkk9Vq9a1STYUk/4C1rKwCJg6yp7Nh0wq6dGhN+4ceKvNMZ/5m2DDwHHjtdBrxqtyO7a5du+J0Ovnhhx+CI7CK0qdPHz799FNsNqPh48iRIxw4cICRI0cyduxYbrnlFp56yoghduzYMaKjowkNDeW7777jwIEDZ/P45JNPyMrKOpsHQI0aNYo9c54/CKahOAQ09lhvBPzqJc0SETklIn8Ay4ELA7BrNF4o1oA1ezrHv+rNj7t+o1fTHcQOLP++iaQkYxS5yePpdTrPdfXt0qULSinSSzIDkqbMtG7dmpdeeom+ffsSFxfHVVddxcyZMwkJCeHmm2/m6aefJiMjg9TUVG655RbWrVtHp06dmDNnztkw47GxsYwZM4aePXvSrl07HnvsMQCGDh3K+PHjad++fXB6tIlIUBaM2sI+IAawAJuA2HxpWgHLXGmrAVuANkXl3bFjR9FoRERWrxZ55RXj0ytbXpFvR5sEkKVPm0S2vBJQfWVh8mSRkBARk0kkIuL8c2zTpo3069cveOKCwLZt24ItocLg7VoB66SA/9Sg+ShEJE8p9SCwFDAD00Vkq1LqXtf+90Rku1JqCfAj4ATeF5EtwdKsqXgkJRUxojk6mTX7TICTzs3L7/gJb4wcCW3beg9NnpSUxLx583A6nZhMeriUpmwENcSkiCwCFuXb9l6+9fHA+EDq0lQhrEks/7kNl1r3kXvppHI7fqIgCjKEXbt2ZerUqezYsYPWrVsHXpimUqFfNTRVGtvGjaRnbKFxronUx9/E7uO+98HCHSa9qvkpRIcuKZLSXCNtKDRVmvWfvs+JvDyaRUScHT9RGWjRogV169atUuMpwsPDycrK0saiEESErKwswsPDS3Scnt1EU3Wxp7N/3wcANK8WjinEVG7HT5QUpYzJl6pSjaJRo0YcOnQIu90ebCnlmvDwcBo1alSiY7Sh0FRdbGlsyzpN9TDo3/cEDa8aWO67xpaErl278vXXX58NA1HZCQ0NJSYmJtgyKiW66UlTdYlOJn0PdGkGcT1PYe0+JNiKfIrbT7F27dogK9FUdLSh0FRZ/opsx6YDisvqxWKvP6XC9XgqioSEBEwmU5XyU2j8gzYUmirLso8/xuF0UmP/iUrV48lN9erVadeuXZXyU2j8gzYUmgpNYbO9FUXaImMITzOPiLGVjaSkJNauXYvDc2o/jaaEaEOhqbC4o8M+95zxWVJjsePPP6kfFkatsLByHTG2LHTt2pWTJ0+yZYsOaKApPdpQaCos+aPDnp3/uhiICBu2baPHlVcSN2pUuY8YW1rcDm3tp9CUBW0oNBWWYkWHLYCDBw9y+PBhkvv3J3bEiEppJABiYmKoX7++9lNoyoQeR6GpkLhnritsMqLCWLNmDQCJiYn+kFduUEqRmJh49nw1mtKgDYWmwpF/5rply0pmJLCnk/7VRCIiwoiLi/ObzvJCYmIiCxYsICsri6ioqGDL0VRAdNOTpsJRFt8E9nRI7cOa9FV0uiSX0KPr/KSy/KAH3mnKijYUmgpHWXwT2NI4k32GDT9BYjMBW5pfNJYnOnXqhMlk0n4KTanRTU+aCkdSktHc5G3CniKxRJH5syInDxJbhFaoiYpKS2RkJHFxcdpPoSk12lBoKiRFzlznDXs6bHiENbudACT+bWylC9tREElJScyePRuHw4HZbA62HE0FQzc9aaoOtjRwnCFtk4kG1U2E2g8HW1HASExM5MSJE2zfvj3YUjQVEG0oNFUHSxT2Q2Gs3K64xBRJ6v8tq3TxnQrC3Q1YNz9pSoM2FJqqw58b2botnD9y82geUQ1nXl6ljO/kjebNm1O3bl1tKDSlQhsKTZVi259nAGhWyWa0K4r8A+/KEkxRU/XQhkJTZbD/1Y0l6eGYgaYRYXQYdUulDd3hjcTERLZt28b//neMPn1gzBjo3h2eeirYyjTlHW0oNFUG2+7T7D6VTdOICCwmCzmOqjVKOSkpCRFh7twfyM4GEXA64bXXYMqUYKvTlGeKZSiUUtFKqb8ppR5QSt2llOqslCqzkVFK9VNK7VRK7VFKPV1IugSllEMpdUNZy9RULHzZRFL3Ugv7Tp+meUQEmKpOs5ObhIQElFKEhFzop3jrrSAI0lQYCv2zV0r1UkotBb4G+gMNgdbAs8BmpdQ/lVI1S1OwUsoMTHLl2xq4SSnVuoB0/wKWlqYcTcWlrPNNnIc9nR9mPUKOCM0iwpG8PI7u3u0zrRWBWrVq0bp1aw4dSqdHj/P37dih/RWagimqVnANMEJEEkRkpIg8KyL/EJEBQDtgI3BVKcvuDOwRkX0ikgN8DKR4STcK+AywlbIcTQWlTDGd8mNLI3W9AqBZRDUADn7zTZk1VjSSkpJYs2YN48YJpnxPf5mur6ZSU6ihEJEnROTnAvblicgXIvJZKcu+GDjosX7Ite0sSqmLgb8B7xWVmVJqpFJqnVJqnd1uL6UkTXmiTDGd8hOdzAH5i1pmM/VCjYAEjfv29YXMCkViYiJ//vkn9UyfkP7+OLpfno7JBGFhZby+mkpNsUJ4KKU+BB4UkWOu9abANBHpU4aylZdtkm99AvCUiDiU8pbc40CRKcAUgE6dOuXPR1PBKOt8ExdgTWJbdjQJrRUNO7Sl8bVDaD5kiG/EViDcA+/SZ93G7d0dfP+8mSV/vEOdhJFlu76aSk1xYz2tBNYqpR7DeOt/Ani8jGUfAhp7rDcCfs2XphPwsctI1AOuUUrlicgXZSxbU44p83wTXsjKymL3voNcN2gQbR95pkp1i/WkVatW1KwexppdZ7i9O5hwck29+6FZW0BbCo13itVzSUQmA3cDC4D/A64QkS/LWHYG0FwpFaOUsgBDgYX5yo0RkaYi0hT4FLhfG4nKj099Ey6WfvQRAJEbN5I6fHiVCd2RH5PJRJeOcazZ47nVAftnBUuSpgJQ3O6xtwHTgWHAB8AipVS7shQsInnAgxi9mbYDn4jIVqXUvUqpe8uSt6Zi41PfhIvlS5diAmLCwnDm5laZ0B3eSLyiHz8ehFPZ52/Xo7U1BVHcpqfBQHcRsQFzlVKfYxiM9mUpXEQWAYvybfPquBaRO8pSlqbiUKb5Jgpgm91Ok4gIIkJDMYWGVrkxFJ4kJSXhdELGPkh2dUjf+2d7+vwNzpwBkwkmTYKRI4OrU1N+KJahEJGB+dZ/UEp18YsijYZSzjfhDXs6jt9Sydy2hcHXXUfcVVcRnZBQZX0UAJ07dwZgzR5Fcmuj38ex/Rs5c8YYqe10woMPQtu2PvoNNBWeogbcPauUquttn4jkKKV6K6Wu8480jaaMuObH3rboeU6c/IvePVsSO2JElTYSAFFRUbS4rAlr9p7rSdiu5nS6Nj/X5uRw6HEVmnMU5aPYDHyplFqmlBqvlHpSKfW8UupDpdRm4HpAz9iuKZ/Y0sCRQ/oeY0a7pEv/Cq6eckRit2TS94Yhro7kZuXgjafSCA1Fj6vQXEBRTU83iEg3pdSTGCOjGwLHgdnASBE57W+BGk2piU4Gs4X03dnUqyFc1nFQsBWVG5KSkpg1axY/HQkjpl4emCx0vi6Z7783ahJRUedqFLr5SVOUoeiolLoEuAXolW9fBKANhab8Yk2C3stIHz2ApK4tUdFdg62o3HB2xjvLaGLiLIZRtSaRZDX2JydDbi6EhhoGQxuLqk1RTU/vAUuAlsA6j2W961OjKdccMV/Ozv1/kHhF/2BLKVe0adOGyMhI0rdmQexow6i6mDXLGL8iYnzO0kMsqjxFxXqaKCKtgOkicqnHEiMilwZIo0ZTar5xDbRrZbUGWUn5IiQkhISEBD01qqZYFHdk9n3+FqLR+Bp7ZibzXn4ZBRyfNKnKjsYuiMTERDZu3Mjp0+e3IA8bZjizlTI+hw0LkkBNuUHPcKeptNgyMth14gRNwsOxOJ1VejS2NxITE8nLy2PDhg3nbU9Kgu++g5dfNj61f0JT3JHZGk2FI6pDB/aePk232rWr/Gjs/NgzM6m9xwj4tGbNGrp163befp8NeNRUCnSNQlNpsVssnHY66Z2SQu9p06r8QDs39sxMUocP55cPPiDaYuH7JUuCLUlTztE1Ck2lJd0V3W7IM89gbdYsyGrKD7aMDBw5OeB0cllEBGt/+KHIY9zzg/gq9pamYqFrFJpKy5o1a6hXrx6XXXZZsKWUK6ITEjCFGO+IzSMisB0/zqalBU9J79O5yzUVEm0oNJWW9PR0EhMTKWp2xKqGNT6eht27A9AsIgKARdOmFZjeH/ODaCoW2lBoKiVHjhxhx44dJOl2Eu+4jOcl4eGEKkXmnj0FJvXH/CCaioX2UWgqJWvXGrEqtaHwTkRUFAAhShETHs6OP/4oMK0/5gfRVCy0odBUStLT0zGZTCToLrFeiUlJYd8XX+DMzaVZZCTfHj5Mzi/LsRxddTbukye6u2zVRjc9aSol6enpxMXFUb169WBLKZdY4+Pp+MwzKLOZZuHh5OTmkvbatbDpOUjtY8zlodG40IZCU+lwOBysXbv2bIRUjXdyjh5FnE6auxza320QwAGObNivIwFqzqENhSYgpKfDuHGB6Vq5fft2Tpw4of0TRRCdkIDZYiEqPJy6oaHsPHnGtUdg3wxdq9CcRfsoNH7H3Q8/J8foNbNsmX/bu1etWgVA1656/onCsMbH03vaNGwZGfT49FPW/7gaOAUISJ4xQ6BVG1uNrlFoAkCg++GvXLmSBg0a6IF2xcAaH0/siBH0ue46fj58kp//tIAyg8k1mZFGg65RaAKAux++u0bh7374K1eupHv37nqgXQlwBwVcFTKaJh4z3mk0EOQahVKqn1Jqp1Jqj1LqaS/7b1FK/ehaViul2gVDp6ZsuPvhjx3r/2anQ4cO8dNPP10QDVVTOO4eYqs22y+Y8U6jCVqNQillBiYBVwGHgAyl1EIR2eaRbD/QU0T+VEr1B6YAXQKv1s/Y0432YHdV3/29EjysnsHkRo/2f3lu/0R3V4gKTfEICQkhKSmJlStXBluKphwSzKanzsAeEdkHoJT6GEgBzhoKEVntkX4N0CigCgOBPd3ot+7IAZMZUODMA7MFei+r0MYi0E5sMAxFZGQk8TqkeInp3r07L774IseOHaNWrVrBlqMphEBH8w2mobgYOOixfojCawvDgcUF7VRKjQRGAjRp0sQX+vyHZw3ClmYYCRzgdLoSCDhzjL7sFbh24c2J7e+beuXKlSQmJhISot1vJaV79+6ICOnp6fTrWKtC33uVmWC8gAXzafLmaRSvCZXqhWEoCmxPEJEpGE1TdOrUyWs+5QLPGoTZAh0mYP+1BrafQoluegZrozNG10QVAvumg9NRYWsXgXZiHz9+nE2bNvHcc8/5t6BKSpcuXTCbzaz85mP6Hfvk3D1aAe+9ysysWZCdDSKBewELpqE4BDT2WG8E/Jo/kVIqDngf6C8iWQHS5j9sadgPmLAdqE1002yw7CD1o0tw5OZiXh1KhweHkvPrZqIvBSv/xahp5FTIPu2BDia3Zs0anE6n9k+UksjISNq3b8/KFcuh0xnACY4zFfLeq6ykp8P06YaRAAgJCUw032AaigyguVIqBvgFGArc7JlAKdUEmA/cJiK7Ai/RR+yeAgc/g8aDsdsbkzqnMQ6HwrxSiLkmG0euw3gmcx2se2su4nRiDjXT4cq65PwF0U1zsV6ZHOyzKBWBDCa3cuVKzGYzXbpUvv4OgaJ79+68959J5OQ5sYQAOMESFWxZFQPPJmU/Gda0NKMpF4xI8XfeWcl9FCKSp5R6EFgKmIHpIrJVKXWva/97wPNAFPCuq098noh0CpbmUrF7CmTcY3w//A22A3fjcJpBBKfTDOH1MVssOHNzQSnE6QSnE0curFtSHxEn5tWh9L4yAquVgNyMFZWVK1cSHx9PjRo1gi2lwtK9e3cmTJjAhv2KxOYCmCCneBX5Kjtdqj3d8CcGoKk4f3PusGE+L8IrQfX4icgiYFG+be95fL8buDvQunzKwc/YvaE2B3fUoHHLE0S32o3ZEoYzNxdTaCgxKSnEpKRgy8jAUrs2G1599ZzRcDgMv3aeE1tGBtaLT5/v3yiHbcfB+rPIzc1lzZo1jBw5MnCFVjLsmZnU278fgJW7zSS2kGKP0A6GgzWouFsJ6sRjXzYF2/4QLBHVyDltJrppNtb8zXU+esEL1twgumuIn9m9rQUZi43OXYf3VyehaVt6TxuOLSOD6IQErK5unO7P2s2bX2A0TKGhRCckgG2xRw+p8ue3COafxcaNGzl9+rT2T5QSe2YmqcOH4zhzhvoWC0u2XsY/Rt9W7D82zx5u2dmGw7XSGgqPVgL7uhWkzrkER965vjnmEKF3QmOs7g1nO7CcAWWCTpOgeelfaIIxN4g2FH7m4KZsPDt4HdyUTfMH488ahvxY4+MvMBpnDYr9tFGTcOaUy1g8wegO62b58uWAHmhXWmwZGTjOnAERWlarxtrNu/j9zNXUt3Yo1vHJyYZj1eEwHK3TpxvNIpXSWHi0EphDBYdDYTzjAiicTjO23aexJrvS29KwH1Ds/zEaFMQc+gfWO9qWq5e8otBBAf1M4759C10vDHfAtrNGxZqEvf4Utu79O/b6U4xtW8eVm3DQwZxbOS0tjZYtW9KgQYPAFVqJiE5IQJnNALSqVo2TDgerFiwo9vFJSYZj1R1ey+Hwf/DHgGNPh63j2L2xDhmLG3J4f3V+2VUDZTK53gUVmEyYLBajBcB9mL0xy2Y3Yc/GOuzZUIdlH16MfeHocvPcFgddo/AzzYcMAeDgN9/QuG/fs+ulwZ6ZSerjb+LIycE8fyu9bz6A9aLjRrTPTu+UqTrrC4LVfpqXl8eKFSu4+eabi06s8Yo1Pp5Ozz7LupdeoqVrVsAdp0+XKI9hw2DmTDhzxjAYUZWps5TH+KeDaY2Bamd31WndhsZ9+mCpXZuco0fPa1IGsO0+jdN57p3c6YDNn++n7ZGrsQ58PejPbXHQhiIANB8ypEwGwo0tIwNHTg44nThzBdv+EKwXOUGcsO4BqB386mzA20/t6WSmzub48eMkB7IKUwlpPmTI2ebOJv/8J+v27CnR8UlJMGECPPigUaN45BFo27aSND95jH+qU/8vDu87ZyguGzy40Oc7OiEBU6gFZ07O2W2H90fy+4FIOtlG0/zJ4D+3RaENRTnA6XSyZs0avv76a9avX8/evXs5deoUERERNG7cmC5dutCvXz9ad+x4tiutKcREdNMz5zKRPNj2GvT8PHgnEmhcb3lpC7MB6BkbGWRBFR+3j6xPejoLFizA6XRiMhW/hTory4hE43QG3k/lT84b/2QWrNf14MuVm9njdPLayy9je/hhQkJCiI6OJjY2lt69e3PttddSr14943rOmMH+BQs4sjmDI9v3AQpxCuuWRFO76zysfyvfF0kbiiCSnZ3N1KlTmThxInv27MFsNhMXF0fHjh2pUaMGp06dYv/+/bz55pu89tprNG3alPtuuom+jRtzcVISHFzAD59MBhFi4o5h5QujR0YFqMr6BFecrO+2CS0vggZqM3BdsFVVCpKTk5kxYwabN2+mXbviR/cPdNiWQGHbfRqH08yh06f57I8/WP/6VBxOJ40bN6Zdu3Z06tQJh8PBL7/8wuLFi5k1axYWi4UhQ4bw9NNP08ZlgO2ZmXw77Daj6zsKcYJtxwmsW8eV77FRIlLplo4dO0p5xul0yieffCKNGzcWQJKSkmTmzJny559/ek1/6tQp+eijjyQpKUkAadGihXw5Y4bMjY+XOa1by5zWrWRu25Zi+1eEyOLOgT2ZYGJbLb++WksiQ5XcnmgRsa0OtqJKw08//SSATJgwocTHrl4t8sorxmdlYf/q1dLPahUFEmEyyQO33y7bt2/3mtbhcEhGRoaMGjVKatSoIUopueOOO8Rut4uIyK5PPpGP2raRObGt5eP4ONkwvL4su7ax7PpH3aDew8A6KeA/VYmU3/h5paVTp06ybt260h1c1IixKVPgs89g8GCjAbY4nluPY47ccAN3DRzIghUraN+iBePffZc+1aqdywe8f09KQkRY9MYbPPDyy/x89Ch/q1ePv1mtmJTRNa9Z+z/pfI0dGr4La7LO6UpPNzq2g+Fx3Lz53Dm4B6jlP2/3+tGjkJl54fnm01bg9XOvR0UZ7RLF9XIXdVx6OvZ585j65QLG7NnHQ02a8OyCBRd2O/bUA/D007B9O7RqBa++amxzX5v27b2X5fmb579eUVGwcSMcPgwNGpzLw1O3+1oVdQ283Xvetk2ZAtOmQXg41K1rlFvSvqjp6ca12LEDWra84FrYHQ46zJpF25gYFm3fXryRlO77zH0tPDUV99y85QcXnl9J9IDxu3j+TgVdLy/5bt26lb/97W/s3r2bm5o354VHH+Vys/nC/wEvv++RI0d49dVXmfDmm9QMCeH9nj0ZOGgQ9u+/x3b0KDmND7N9xblOAwl9Q2gefVPxz8uHfZCVUuuloMgXBVmQiryUukaxerVIRISI2Wx85n8lmjxZxOgmbiwhIQWn9XLMepCmtWpJKMjrSklueLix312mxSISFmZ8Dwsz1j3zd+k7bjLJjaGhAkjH6tVlRsvLZU7rVvJRm5ZiGxch8jfzueMmTzbycms2m88/h8mTLzxvtyalzk8bGlqwNm/Xz52PyWQcbzIVfq3y/w4FHefav6VePbkpOloAebd5c9ny7LMF/55hYefycy8mk3FO+bd5lpX/N/e8Xvnzcy/u62YynbtWRV0Db/eet2359bgXi6X4r/CrV194H5hMxv0MYouIkLktW0rPWrUk0mSSw3ffLRJrEUlRxqe3clavPv8+A2O9oPMo6lnLn5/n+RV1bEF6irpeXvL96quvJDIyUhpERcl37ns+/z0TGlr477t6tWwJCZFOxkALeRIk13Xs4m6XyJzWrc62DCy+4pLCn5PVqw3t+a+xD6CQGoUeR+GJtxFjnnz22fnreXkFp813zHSgK5B36hQrTCYeFyEkN9fY7y4zN/f88nNzz8/fpa+G08nHubm8FBrKhpMnefXAz/zlcCJOxbb0KNjiOHecO3837ohinvryn7f7mPy1TU89+bV5u37ufNzzbHh6OIvzOxR0nGt/9KlT7Dh1iossFuqazUSfOuU9H7ees/N9cC7f3NwLt3mWlf8397xe+fNz475u7rxyc4u+Bt7uPW/b8utxk5tb/IELnpHl3Didxv0M7K9VC6dStIqM5JTTyf9WLoHHc2CwGJ+rZ3nP0/M+g8LPo6hnLX9+nudX1LEF6fHE2/XKl+9HEyYwcOBAWrZsyfoRI0h2OM49p/nzKuz3TUsj1uFgJXAP8BrQF7ADEafO/x3M4QKXFvKcpKWdX35xnicfoA2FJ0WNGBs8+Pz1y82QoqCluUCvnQwaxFiMyTR6ABvuvZcuYWHnyhg8+FyZoaHnlx8aer4WD33KYmEM8H+1arHn9GnG//wz2U4nv+6ujv33mhfm78Y1qOq8c8p/3u5j8vd28dSTX5u365c/H5OpeB5Odz4FHZecDC3N1LnyNLtPn6K7yUTvw4ex5u+imF9P/vMxmYzzyL/Ns6z8v7nn9SqoN1B+3aGhRV8Db/eet2359bgJDS2+5zg5+cL7wGQyhlZ70LKa0QV0g0WMbi9mjM/WBeTpeZ9B4edR1LOWPz/P8yvOyE5vejzxdr1c9xUpiv82hFvnzaN79+6kpqZy0YAB5z+n+fMq7PdNTobQUMKA94APgHSMyXVq7vjDNUjReLnI+q0a9l4RBf+nuPI6S4B6DOheT54UNWLM3T792WcwJB77n+9h+8lCdLccrM0uzE5EeHznTt4EhjVsyLTnniPkvvvg5pvPL6Ogtn/P724tnvqAe955B/vy5bx96BATDh7iH5dcgu2+G7E2/hlaD4ZeI438i/JR5D9vt6YCfBTpm6uT9lkWyYOjSEpqW/D1K6T9tli/g7fjmgFjFBt3w8kvYdAN3bHe/38X5ptfD5TcR+H5m3u7Xr7yURR07+Xf5t5eFh9FUhKsWFGgjyLm+HH2/fgjVqWIDg1lR80GEPKHETomxAKdvIQsTUqC774r2EdRnHMrKD84P6/ijOzMf3xxfBTNwH57OPOXh/Pgr8fp3qUdX3/9NdWqVbuwTM9nqKh7PCnJ2O/Scnv79jRbtIhrFy8mpVYt3rBWI9f2F0b4D9j8az3a3n7M639K/rwCFieloDapirwEoteTbf6j8nFcS5kT20o+jmsptvmPnrc/Ly9P7rrrLgHkoYceEofD4Tctuz75REY0aiSA9I6qK7+/Vktkjlnk4wi/9KIoThOxX9nyisgcs4z7u9Hm+1va6AALqPzYNm6ULVOmyM0pKVK7dm3J/XW5cd1LeT9VhJ5QtvmPyr9bxEg1k0mahofJ3tn3+7W8jRs3SnR0tNSrU0devbyFzIlt7fJXeP9P8TdoH0XxSU+HceOMz8KwHYg0goGJwulQ2A6cG+wlItx3331Mnz6d559/ngkTJpRo0FJJaT5kCC9/+SV39+9PatYR3l95mq2ramP/2WSMNfAxxWki9ivRydh/rcHnK0OJbWiiQevrAyyg8uOOMzbgpps4evQo6w6EQuzoUvXzd0cVfu4547OoZytY7N1hZvxPv2BWikcbN+a0va5fy4uPj2fFihVEVK/O+N9tHGnonvBT4XAoflxbfgaQakPhQUlu6OjkwZgtFpRJGUHAko22YxHh8ccfZ+rUqYwZM4Z//vOfKHekND9ijY9n8ldf0a9rB57/Mof5S2qQOrsxdnvjog8uIcEM/gdg/yWCRR82Yv2vuVzqiML+S0RgBVQhrrzySpRSfPPNN6XOI+gvFsXA6XTywpIfsOXk8EiTxtSvHnn2mfYnLVq0YNmyZWAO58nl6/nljCAocp0WMnP8X35x0YbCg5Lc0Nb4eHpP/4C4hx6m9/QPzvbff/HFF3nzzTd56KGHGDt2bCBkn8VkMvHckNuoHRLK2wd/4WSuwra7ZIHdCsNd2wKjuXbs2OBMUGPLyGDLsRM4BNpGVseWkRFYAVWIqKgoOnXqxNKlS0udR7BfLIrDSy+9xDfLl/Py6NH8ffQz5z3T/uaPP5pz8uT/OJ0jPLPvGLN+u43xv31A1yGBKb9YFNQmVZGX0vooitP2vmvaOFk2tI/smjbugn1vv/22AHLXXXf51SdRGLaNG+WfLVqICaRL7dry+4YNPsk36H4JD2wbN0o/q1UsSsms+HixbdwYPDFVgDFjxojZbC4wckBxKM8+irVr14rZbJZbb71VnE5nwMt/5RX3sJZ1AjUlMrKVLFlyJOA6KMRHEfQ/dX8sZXFmF3ZD75o27qyzaU7rVucZi/nz54tSSlJSUiQvL6/U5fsC28aN8tigQQLIe++955M8z93Mxucrr/gk21LTPCZGusfGaiMRAL7//nsBZP78+cGW4nNOnjwpzZs3lyZNmsjRo0eDouGWW4wxmkqJWCzfSWioRa644go5ffp0QHVoQ+Ejlg3pet4oymVDuoqIyKpVqyQ8PFwSExPl1KlTfim7pDgcDrnyyiulevXq8tP6z8rUY0WkfNUoDhw4IID8+9//Dp6IKkROTo5Ur15d7rnnnmBL8Tn333+/AJKamhqU8p98Us4b5H3LLSJz584VQG688caAtkwUZii0j6IENO7kngVXzq7v3LmT66+/nsaNG/Pll18afa7LASaTialTpyJOB/fcNQTJfNaYeKWUs2q5u5EHyy/hidux2rcEswVqSk9oaCi9e/dm6dKlxttlJWHp0qW8++67PPbYY/Tq1SsoGubPP3997VoYOnQo48eP55NPPuGJJ54Iiq4LKMiCVOTFb+MobKtl1xP1ZNm1TWTXE/Xkt81fStOmTSU6Olr27t3rnzLLyNvPXC+AvNQ7Umz/qiay9t5gSyozQ4YMkYsvvjgo7clVlXFPPy2ArFmwINhSfMLJkyelSZMm0qpVq4A38XiSv0bx5JPGdqfTKQ899JAA8uabbwZEC7rpyYfYVotseUVO7PtWOnbsKNWqVZMffvjBf+WVkcPLZsrl1SKkmskk/2nVTGzja5aoCWryZJG+fY3P8kBeXp7UqVNH7rzzzmBLqTLYNm6UCa1bCyB3Nm7sM79QMB3c//jHPwSQFStWBL7wfDz5pEizZueMhJu8vDwZNGiQKKVk3rx5ftdRbg0F0A/YCewBnvayXwETXft/BDoUJ19/j8zOycmRfv36idlslq+++sqvZZWVLVOmyOvNLpNQpaRzjRqy9taGhr+iGHgLnBps1qxZI4DMnTs32FKqDFumTJHZsbFiDQ2VjjVqyJYpU8qcZzB9XpmZmWI2m2XEiBGBKzQfxTWSf/31l3Tr1k3CwsJk+fLlftVUmKEImo9CKWUGJgH9McKM3aSUyh9urD/Q3LWMBP4TUJFeEBHuuecelixZwn/+8x+uvfbaYEsqlOiEBC6OjCSlXj1+OHGCL1aGFHsQnrfAqcFm0aJFKKW48sorgy2lyhCdkEBIWBhxNWqw5dQpasbFlTnPYA3Cczgc3HPPPdStW5dX3fGtAsyUKXDFFfDss0UP7I2IiGDBggU0bdqUAQMGsHXr1sAJ9SCYzuzOwB4R2SciOcDHQEq+NCnALJfBWwPUVko1DLRQT8aMGcOMGTN44YUXGDFiRDClFAtrfDyX/m0Q19aLon5oKB/8+ju/bD9R5HHp6ZDfL19Q4NJAsmDBArp160a9evWCLaXKYI2Pp/e0adxwyy2ccTrZnJVV5jyDNQjvjX8MZ+3atfxzxGDq1vVviA5vpKfDgw8aEd2dTjhzpmgjGRUVxZIlS4iIiKBfv34cPHgwIFrPo6Cqhr8X4AbgfY/124B38qX5Cujusb4M6FRU3v5qepo4caIAMnLkyArlSLVt3Cgfd+ggTzZtKoA8+9BDhab3bBYIDRXp3Ll8NDvt379fABk/fnywpVRJsrOzpXr16jJy5Eif5BdoH8XOBf+R6maTxEZWk7ltLxfbdx8GpmAPXnnl/PmuQkKKf/6bNm2SmjVrSqtWrSQrK8vn2iiPTU8Y/of85O97V5w0RkKlRiql1iml1tnt9jKLy88nn3zCww8/THz8QPbvf5epU/0fv8lXuN8Ihz3zDP2Sk/n3++8X+lbi2SzgdMLAgeeiaweThQsXApCSkr/iqQkEYWFh9OvXj4ULF+IsaNKmEpCUBKNHB66r9XPPvswph5Nb6zdAnCZsKxcHpmAPkpMhLOzcVCiTJhX//OPi4liwYAF79+yhT+fOHAhkdMWCLIi/FyAJWOqxPhoYnS/NZOAmj/WdQMOi8vZ1jWLZsmVisVikWbMeAn+VK+duSdm/f7+Eh1vkxqvbFtj7qTwNrvOkV69e0rp162DLqNLMnj1bAElPTw+2lBKxYsarYgK5sk5tmdO6lcxt2zLgNQp3DWry5NLXpGwbN8ojMTGiQDrWqiW/ZmT4TB/lsdcTxqRJ+4AYwAJsAmLzpbkWWIxRs0gEfihO3r40FCtXrpTq1atLbGys9Op15LxeQH37+qyYwGFbLU9cFSaALLivWqHGojzF5snKyhKz2SyjR+u5J4LJkSNHxGw2y9NPPx1sKcXG6XRKYuO6Us1kkvdatJA5rVvJ2vv7BFTD5MnnptYuy8vXlilTZE6bWLmjQQMB5Oq2bX0WMqgwQxG0picRyQMeBJYC24FPRGSrUupepdS9rmSLMIzJHmAqcH8gNa5Zs4b+/ftz0UUX8b///Y+hQ+uct788OHdLiv1/7xD76yVEhYTw0AwHh5dO9Jou0M0CRbFo0SIcDodudgoyderUoWfPnixYsCDYUorN3FefYM3BIwyyWqkRYsZkFmKG3BGw8tPT4YEHzk2tXRwHdkFEN4/ApBxcVbcOf4+2snTzZm5KScGRfw50HxPUqVBFZBGGMfDc9p7HdwEeCKQme2YmtowMDkVEcOMDD1C/fn1SU1Np2LBhgbNiViT2r/qdUDFzU/36vPPLL0yclskrtwZbVdEsWLCAhg0bkpCQEGwpVZ6UlBQefvhhdu3aRYsWLYItp1BycnIY88YUGlosXOXq5XRpj4uwJgfupk9LMwyEG3MB02EXB6v1IJe2O8aeDbUZUM+KAJ98/TWRd9/NtGnT/DZBmo715IE9M5PU4cP56JVX+Nvtt1OnZk2+++47Lr744rNpRo6EpUuNz+LOhleuqGk82Ik1a3J5RATv/fATR48eLb/nYk8ne8P/sXjR1wwYMMCvMwVqioe7VlcRahWTJk3ip6wTDLsomlATmEOEmCF3BlSDpwM7JATeeacMNfXoZGLiTmAOEZQSBjWI4om/J/LBBx9w++23k5ub60vp5yioTaoiL6X1UWyZMkUebNxIQpSSRmFh8u2//lVg2vLq8C0K28aNMrddW5kT21pead5MlFIydOhj5fNcbKvFNr62vHttTQFk0dw3gq1I4yI+Pl66devm1zLK6iez2WxSs3p16R4bK7umjZMtY28OSpdYER/7/HZNFttrNWTL/fXENr62OH9fJS+//LIAcv3115c6dhXl0Zntz6W0hmLsA7cIIC2rRcj7rZsXelOVt/kZSoJt40bZMmWK2DZulLvvvltMphAxmXaUu3OxzX9UPo5rKT1r15IIk0kO/ndUsCVpXPzzn/8UpZQcOnTIL/n74kVs2ODBYgJ5rdll8nH7dpVr7hJXzDnPziiTJk2SgQMHSm5ubqmyLMxQ6Hq8i6ysLCbO/YJr25iZc2sEA+74Fau14LEGFWF6x4KwxscTO2IE1vh4XnrpJSIiqgGPlbtzsR2IJDtXyDh+gk41anDs11rBlqRx8fe//x0RYd68eX7J33MsT3Y2zJpVsuN//PFHZs+fz1V163CxJQxnTg62tHIQg8ZXWJMgdrTx6eL+++9n/vz5hIT43vWsDYWLqKgo1iz5gAXPhNIh+RjWJk6ITi4wfXman6Es1K9fn3/+83mczkXceuvicnUu0cmD2Zx9hr+cTrrVqxuQye41xePyyy+nffv2zJ071y/5Jycb7flgdEafPr34/jMR4ZFHHqFmRAiDo+sBgsksRF9yyi9ayxNK+WcgsDYUHlyacAPmK1Mhbiz0XnaetfZGeetCWlzyO65HjRpF8+bNWbPmUTp18pMzrBRY4+PZ27wFdapXZ9Tc/wZssntN8Rg6dCg//PADe/fu9XneSUlw553g/t9zOIrfpfSLL77gu+++42+16xJpCkGZoMPVf2DtPsTnOqsK2lDkx0uVrjKRnm68rY0ZY3ymp4PFYuHFOwewc+dO/vVoYHuEFIQ9bTY/PHcjS79PY+htt3GR7hZb7hg6dCjAuVqFPR22jiv1LIr5GTYMwsNL1rybnZ3N448/TvPoWvSqXQdQIJAT3tnvz3S57TnoA4I6jkITeGbNMtp+wficNQuanZmNae5XxFWPZNx/5jL4yq60GhjQsY3nYU+bTerDL7Psj2OcPpPDdbHWog/SBJwmTZrQs2dPZs6cyZiRvVHfXQmOHDBbilUjLwp3825ammEkilNznzBhAvv37+fT158mb9YXOB2uZqfeN5VJS1Gkpxshw3NyDKNWnppwfYGuUWiwrVyM02ni1voNyHY6eeGVN4Oux+FQfP/nURpaLDQ5sjuoejQFc+edd7Jnzx5WLvnAMBI4wJkDtjSf5F+S5t3ffvuNl156iZSUFAY/Po7eb40hbmh7er81xu8D7II1v0ag0IaiijFsmDH4Rynjc9gwiO7eH7NZaBRu4ep6dfh03V4yMzODpjG6e39seWfYdfo0yXVrUb/HNUHToimcG264gerVqzNjya9GTQKTsVii/Fqut2aeZ555htzcXF5//XUArMm3EvvsHL8bifR0+Plno4msvPUc9BkF9ZutyIu/p0KtiHgO+PE2+Mf23YeyZezNsmvhexIVFSU9e/YM6pwbD918nZiUkh/nTQyaBk3xuOuuuyQyMlKObZgg8lGIyByTyMcRJZqbvSR4G2Pxww8/CCBPPfWUX8osjpawMJF77y1HA1ZLCHrAXXAJdiTWyZONCVKKG7nyP//5jwABmdDdG9nZ2WK1WiUlJSUo5WtKxtq1awWQd8ZcLzLHLDIHkTlKZO29fikv/2DXl192SlJSktSvX1+OHz/ulzLzM3myET164MCKO/A2P9pQBJHVq0UsFhGljM9AG4vVq43wxu7Q6CZT0TdzXl6exMXFySWXXCJ//fVXYIR68OGHHwog33zzTcDL1pQcp9MpHTt2lBaXXCSb748W278iDGMxN8wvtYr8NYrnnpslgEyfPt3nZXlj8uRzzxMYz1e5C39TCgozFNpH4WfcvYxEzvUyCiRpaYaDzY3JVHT7qdls5sUHH+TAgQP83z/+4U95Xpk0aRItWrSgT58+AS9bU3KUUgy7/np2HfiVuYsiSJ1zCfZDESB5PnNqe+I52PXzz4/w3nuPk5iYyO233+7zsrzxWb4B3u3bV/yBt0Whu8dWQtLTz3UpdEeuPHPGcLQVJ3KlPTOTnP9MIqFmDd587z2GDhhAu6uvDoByWLt2LWvWrGHChAk6UmwFolu9elQ3m1nyRxZtIhtj+7k61ktOFhrdoCwkJRnLiBFPceTIEb799lu/3y/u5yo+Hr755tz24cMr5pQDJaKgqkZFXspb01NYmNH0FBZWsDPZl+Xld/SVtLwtbz4rc2JbyZvNmolFKUlq1FB+37DB92K9MHDgQKlTp46cOHEiIOVpfINt40YZ7Jp17bXmzWXXhDv95sx2s3z5cgHkySef9Gs5Ihc+V08+afgoKuJ0yAVBIU1PukbhZ5KS4Lvvzr3hg+8H5kyZAtOmwUUXQYMGF/bnLmmYkehLTmE2C9FhoQyJtjLn0G/838CBvLBggV/DaGzfvp0vvviC5557jurVq/utHI3vscbH88RLL/HViBF8abfTZNYmaveMwOqnsZI5OTncc889XHLJJTz//PNe03jWrMv6jOUfJ1G7tjEvTVVBG4oA4K4mg9H329vAnNLc0Onp8PTTsHz5uW3uvtxQ+v7c1u5D6H3bDNYvrUM/Z11+OH6Cmb/8wrWzZ9Pfj4bilVdeISIiglGjRvmtDI3/qO500qtuXf6XlcWgEyfY78cXi2dGjWL79u18NHEikZGRF+x3j5T2bHItS/OQO1q0+wWv0o2TKAJtKAJM/hsuKur8GsaECZCVde5GzG9A3G9JUVHw0EPGg+CJwwEjRkCTJmV4k7ImYb19EXX3P8eRX39l5EUX8cy+fbz25Xz6jR/vlwiVmzdvZs6cOfzjH//A6q/XUI1fiU5I4PoGDVh25Aif2mw0+OILYlJSfG4svv7gA/49ZQo9atfG9MEH2Hv0uKCMtDTj2XA6jeX++2HxYmNfgwbGQFNvz0ZBtZDShBOpVBTUJlWRl/Lko/CGp8/As0+4yXSuq11YmNGdNr+vwd1OGhJi+D08u+mBsd1Xvg/bnJtlbtuWMqd1Kxla32q0P999t18mgLn++uulVq1akpWV5fO8NYFj7YsvyvX16okCebVZM9kyZYpP8z927Jg0qldPrKGhMvXyy+WjNm28lpG/W3j+xVtX9Yo6a6WvQHePLV94xq/xnADJZDJqBO5mqdzc85uoPNtJnU4jvRul4IorjGYoX73tWJvVpM+tB2iXbOP5m3NoXi2C56dP56Nbb8XuwxAfy5Yt48svv+Spp56ibt26PstXE3hiUlJIuegiIkwmPvz9d6ydOvk0/1GjRvHrkSwebNKQyBATJlMe0c0jLkiXlGQ0NxU0h09u7oXxmCp7vKayoA1FkPHsEz5pktGV1R0vJjT0/NgxnkYlLAzefRfuvddYVq2C77/3cZU4ZhjWJk5iux3BmRPOA40aYQLe2rePHZ9+ytapU8tsMLKzs7nvvvu47LLLePTRR30iWxM8rPHxDJg5k8eGDmXriRN8u327scMHIcjnzp3LrFmzeHZkL+4bmUVcTxu9bz1Y4EyUI0caL0733ms8S56Ehl7oZ6jIs1b6G2XUOCoXnTp1knXr1gVbRqnwbCOFgn0UAWsntaeDLQ27vTGpj73OD1lH+PfPB+lRpw73XHQRIWFh9J42rdTt0M8//zxjx47lm2++4aqrrvKtdk3QcDgcJCUlceDAAbZ+/wH1Ng0uUwjyDRs20L17dzp27Ejqf18hdMXVRpRaU/HyS083BrsePlw6H0VVQCm1XkS8VgG1odAUD3s69pnXYNsfylsrzzD1x+PcGB3NwPr1iRs1itgRIwo/PDMTW0YG0QkJZ43K8uXL6dWrF7feeiszZ84MwEloAsmmTZvo3LkzfZMuZVr/n7H/HEF002ys14wxJgcrJocOHaJr164AZGRkUL9+/bMvMEQnV9pJxgJNYYYiKL2elFJ1gf8CTYGfgBtF5M98aRoDs4AGgBOYIiJvBVap5iy2NKwXncB6kYOxTaqx9edafGKzUSs8nCs9Zp/LbxDsmZnsX7CAfV98gTMvD7PFQu9p03A0aMDNN9/MZZddxjvvvBPEE9P4i3bt2vH666/z0EMP8dhOK9dE1cO8Uuid0Jji9GuzZ2ay8aOPuGfGDP48fZoVK1caRgIM46ANRMAIlo/iaWCZiDQHlrnW85MHPC4irYBE4AGlVOsAatR4Ep1sNBsoM/WbCp9MG0WP5g15/+BB5q1eDRgPdurw4WyaOJHU4cPZPW8eqcOHs2fePJw5OeB04szNZe/339O/f3/+/PNP/vvf/1KjRo3gnpvGbzz44IP07dCBOYftrDp6HKczBNvu00UeZ8/MZN4ttzB84kR+ycri0QbRXBwAvRrvBGscRQqQ7Po+E0gDnvJMICK/Ab+5vp9QSm0HLga2BUyl5hzWJKMt2JYGligu3vAIS585w9B3FA888AA//fQTN8fE4PAwCAe/+cZYdzdvKsVx4MFp09iyfTtfffUV7du3D+ZZafyMUor3332Xq6+8ksm//IKEhJytgdrTZmNbuZjo7v0vmFxo2bx5PL97N6ccDp5o0oSWYeHY0j7za2QATcEEy1DUdxkCROQ3pVR0YYmVUk2B9sDaQtKMBEaCMZevxg+4q/tbx4EjhwiLk88eMfHQV10YP348qW3bMig3lxiLBVNoKI379sW+YQPO3FyUCWwx1RmfvosjJ04zb948rg5QoEFNcGncpQtfLlrELSNH8p8dOzC9/z6P9otj3Zg3cTgU5nkb6P2WMSPdqVOneOONN3jptX9RSymevaQJMRERxrzXl5wK9qlUWfzmzFZKfYvhX8jPGGCmiNT2SPuniNQpIJ/qwPfAyyIyvzhla2e2n7GnQ2qf83qd/Df1Z0aNGoXdbqdbbCyDBg2ic9++HFm/hPXfLGDZzu2s2uugWQPFf2dPp0OfO4J9FpoAk5ubyxNPPMHEiROpHWEhMSKS2GqR1AwxU7drE/bV78ScOXP4/fffGXJ1G15K3MvRvbVAIKbdSax3fKP9En6kMGd2UEZOAzuBhq7vDYGdBaQLBZYCj5Uk//I+MrtSYFstsuWV8yKEHjt2TF58YrjENKojwHlLs/rI+JuR7Jkm4zhNlWX9+vXSr1sHCVHqvHvEYrFISkqKrFy50rivPo4wplX9KFRkVyUK01pOoZCR2UHpHquUGg9kicirSqmngboi8mS+NArDf3FERB4pSf66RhEkXDUNyTvDblso+0x9MB9eQpMoJy0aGu3VmMNL1Y9eU/nY+9UU1iz8FOclbbj8qqG0bNmSmjVrnkugu8AGlHI3jkIpFQV8AjQBfgaGiMgRpdRFwPsico1SqjuwAtiM0T0W4BkRWVRU/tpQBImt42DTc4ADlBkuGwH7ZxpNVMoMl94FMcP0Q6/RlEPK3TgKEckCLpjnUkR+Ba5xfV8J+D5MqcZ/uLvQun0XMcOMRb8VajQVGh1mXOM7PLvQehoGbSA0mgqNNhQa36JHzGo0lQ4dPVaj0Wg0haINhUaj0WgKRRsKjUaj0RSKNhQajUajKRRtKDQajUZTKNpQaDQajaZQKuUMd0opO3CglIfXA/7woRxfoXWVnPKqTesqGVpXySmNtktExOucUpXSUJQFpdS6goaxBxOtq+SUV21aV8nQukqOr7XppieNRqPRFIo2FBqNRqMpFG0oLmRKsAUUgNZVcsqrNq2rZGhdJcen2rSPQqPRaDSFomsUGo1GoykUbSg0Go1GUyhV0lAopYYopbYqpZxKqQK7kCml+imldiql9rimbHVvr6uU+p9Sarfrs46PdBWZr1LqcqVUpsdyXCn1iGvfi0qpXzz2XRMoXa50PymlNrvKXlfS4/2hSynVWCn1nVJqu+s3f9hjn0+vV0H3i8d+pZSa6Nr/o1KqQ3GP9bOuW1x6flRKrVZKtfPY5/U3DaC2ZKXUMY/f6PniHutnXU94aNqilHIopeq69vntmimlpiulbEqpLQXs9889VtBk2pV5AVoBlwNpQKcC0piBvcClgAXYBLR27XsNeNr1/WngXz7SVaJ8XRoPYwyUAXgR+IcfrlexdAE/AfXKel6+1AU0BDq4vtcAdnn8jj67XoXdLx5prgEWY8zcmAisLe6xftbVFajj+t7frauw3zSA2pKBr0pzrD915Ut/PZAaoGt2BdAB2FLAfr/cY1WyRiEi20VkZxHJOgN7RGSfiOQAHwMprn0pwEzX95nAQB9JK2m+fYC9IlLaUejFpaznG7TrJSK/icgG1/cTwHbgYh+V70lh94un3llisAaorZRqWMxj/aZLRFaLyJ+u1TVAIx+VXWZtfjrW13nfBMz1UdmFIiLLgSOFJPHLPVYlDUUxuRg46LF+iHN/MPVF5Dcw/oiAaB+VWdJ8h3LhDfqgq8o53VdNPCXQJcA3Sqn1SqmRpTjeX7oAUEo1BdoDaz02++p6FXa/FJWmOMf6U5cnwzHeSN0U9JsGUluSUmqTUmqxUiq2hMf6UxdKqWpAP+Azj83+vGZF4Zd7rNJOhaqU+hZo4GXXGBFZUJwsvGwrc1/iwnSVMB8LMAAY7bH5P8BYDJ1jgTeAuwKoq5uI/KqUigb+p5Ta4XoDKjU+vF7VMR7mR0TkuGtzqa+XtyK8bMt/vxSUxi/3WhFlXphQqV4YhqK7x2af/6Yl1LYBo2n1pMuH9AXQvJjH+lOXm+uBVSLi+Zbvz2tWFH65xyqtoRCRK8uYxSGgscd6I+BX1/fflVINReQ3V7XO5gtdSqmS5Nsf2CAiv3vkffa7Umoq8FUgdYnIr65Pm1Lqc4zq7nKCfL2UUqEYRmKOiMz3yLvU18sLhd0vRaWxFONYf+pCKRUHvA/0F5Es9/ZCftOAaPMw6ojIIqXUu0qpesU51p+6PLigVu/na1YUfrnHdNNTwWQAzZVSMa6396HAQte+hcDtru+3A8WpoRSHkuR7Qbuo68/Szd8Arz0j/KFLKRWplKrh/g709Sg/aNdLKaWAacB2Efl3vn2+vF6F3S+eeoe5eqYkAsdcTWbFOdZvupRSTYD5wG0isstje2G/aaC0NXD9hiilOmP8Z2UV51h/6nLpqQX0xOO+C8A1Kwr/3GP+8MyX9wXjT+EQcAb4HVjq2n4RsMgj3TUYvWT2YjRZubdHAcuA3a7Puj7S5TVfL7qqYTwstfId/yGwGfjRdRM0DJQujN4Um1zL1vJyvTCaUcR1TTJdyzX+uF7e7hfgXuBe13cFTHLt34xHj7uC7jUfXaeidL0P/OlxfdYV9ZsGUNuDrrI3YTjau5aHa+ZavwP4ON9xfr1mGC+HvwG5GP9hwwNxj+kQHhqNRqMpFN30pNFoNJpC0YZCo9FoNIWiDYVGo9FoCkUbCo1Go9EUijYUGo1GoykUbSg0Go1GUyjaUGg0Go2mULSh0Gj8jFIqwRV4MNw1cnerUqpNsHVpNMVFD7jTaAKAUuolIByIAA6JyLggS9Joio02FBpNAHDF18kAsjHCUDiCLEmjKTa66UmjCQx1geoYs+yFB1mLRlMidI1CowkASqmFGLOKxWAEH3wwyJI0mmJTaeej0GjKC0qpYUCeiHyklDIDq5VSvUUkNdjaNJrioGsUGo1GoykU7aPQaDQaTaFoQ6HRaDSaQtGGQqPRaDSFog2FRqPRaApFGwqNRqPRFIo2FBqNRqMpFG0oNBqNRlMo/w9LtQAyJLFKdwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "test_set = 100 # number of test points\n",
    "\n",
    "col = ['red','blue','orange','brown']\n",
    "\n",
    "# generate the test data\n",
    "xtest=[]\n",
    "for i in range(test_set):\n",
    "    rand = random.uniform(-1,1)\n",
    "    xtest.append(rand)   \n",
    "    \n",
    "for l in range(1,laymax+1):\n",
    "    \n",
    "    param = [tq.Variable(name='th{}'.format(i)) for i in range(0,3*l)]\n",
    "    ytrain=[]\n",
    "    for i in range(test_set):\n",
    "        \n",
    "        qc = qcircuit(xtest[i], l, param) # construct again the QFit\n",
    "        expval = tq.ExpectationValue(qc, H=tq.paulis.Z(0)) # compute the expectation value\n",
    "        res_test = tq.simulate(expval,variables=train_param[l]) # use the optimized angles to check the result!\n",
    "        \n",
    "        ytrain.append(res_test)\n",
    "        \n",
    "    plt.plot(xtest, ytrain, '.', color=col[l-1], label = \"{} layers\".format(l))\n",
    "    \n",
    "xfunc = np.arange(-1., 1., 0.01)\n",
    "plt.plot(xfunc,f(xfunc), '-', color='black', label =\"exact\") \n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"f(x)\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions and some open questions\n",
    "\n",
    "As you can see, one qubit has the capability of learning an arbitrary function. The data re-uploading strategy used (repeating the same encoding multiple times) increases the accuracy as more layers are considered. This result comes from the fact that a qubit can represent the Fourier series, which in turn can represent arbitrary continuous functions. It can also be proved the universality of a single-qubit as a function approximator by means of the Universal Approximation Theorem.\n",
    "\n",
    "- These are the results for only one qubit, what happens if we consider more and entanglement between layers? Do we require more or fewer parameters to learn the function with similar accuracy?\n",
    "- We used a particular encoding. Are there encodings that involve less number of parameters without compromising the accuracy?\n",
    "- We only tested a 1d model. What about higher-dimensional functions?\n",
    "- What about using quantum data instead of classical data that we introduce \"manually\" to the circuit?\n",
    "- To construct the cost function, we measure the state on the computational basis. What about measuring it on another basis? (hint: by using a final rotation without the data point that projects into a different basis)\n",
    "- Can we design more clever cost functions?\n",
    "\n",
    "**I hope you enjoyed this small tutorial and do not hesitate to contact me with questions or ideas!**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:tqenv2] *",
   "language": "python",
   "name": "conda-env-tqenv2-py"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}