{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Speed test of stochastic solvers \n", "## Based on development-smesolve-milstein notebook\n", "\n", "Denis V. Vasilyev\n", "\n", "30 september 2013\n", "\n", "\n", "Modified by Eric Giguere, March 2018 \n", "Include new solvers from the pull request #815 for qutip" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Populating the interactive namespace from numpy and matplotlib\n" ] } ], "source": [ "%pylab inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from qutip import *\n", "from numpy import log2, cos, sin\n", "from scipy.integrate import odeint\n", "from qutip.cy.spmatfuncs import cy_expect_rho_vec, spmv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Single stochastic operator" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "th = 0.1 # Interaction parameter\n", "alpha = cos(th)\n", "beta = sin(th)\n", "gamma = 1\n", "\n", "# Exact steady state solution for Vc\n", "Vc = (alpha*beta - gamma + sqrt((gamma-alpha*beta)**2 + 4*gamma*alpha**2))/(4*alpha**2)\n", "\n", "#********* Model ************\n", "NN = 200\n", "tlist = linspace(0,50,NN)\n", "Nsub = 200\n", "N = 20\n", "Id = qeye(N)\n", "a = destroy(N)\n", "s = 0.5*((alpha+beta)*a + (alpha-beta)*a.dag())\n", "x = (a + a.dag())/sqrt(2)\n", "H = Id\n", "c_op = [sqrt(gamma)*a]\n", "sc_op = [s]\n", "e_op = [x, x*x]\n", "rho0 = fock_dm(N,0) #initial vacuum state" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Solution of the differential equation for the variance Vc\n", "y0 = 0.5\n", "def func(y, t):\n", " return -(gamma - alpha*beta)*y - 2*alpha*alpha*y*y + 0.5*gamma\n", "y = odeint(func, y0, tlist)\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on function stochastic_solvers in module qutip.stochastic:\n", "\n", "stochastic_solvers()\n", " Available solvers for ssesolve and smesolve\n", " euler-maruyama:\n", " A simple generalization of the Euler method for ordinary\n", " differential equations to stochastic differential equations.\n", " Only solver which could take non-commuting sc_ops. *not tested*\n", " -Order 0.5\n", " -Code: 'euler-maruyama', 'euler', 0.5\n", " \n", " milstein, Order 1.0 strong Taylor scheme:\n", " Better approximate numerical solution to stochastic\n", " differential equations.\n", " -Order strong 1.0\n", " -Code: 'milstein', 1.0\n", " Numerical Solution of Stochastic Differential Equations\n", " Chapter 10.3 Eq. (3.1), By Peter E. Kloeden, Eckhard Platen\n", " \n", " milstein-imp, Order 1.0 implicit strong Taylor scheme:\n", " Implicit milstein scheme for the numerical simulation of stiff\n", " stochastic differential equations.\n", " -Order strong 1.0\n", " -Code: 'milstein-imp'\n", " Numerical Solution of Stochastic Differential Equations\n", " Chapter 12.2 Eq. (2.9), By Peter E. Kloeden, Eckhard Platen\n", " \n", " predictor-corrector:\n", " Generalization of the trapezoidal method to stochastic\n", " differential equations. More stable than explicit methods.\n", " -Order strong 0.5, weak 1.0\n", " Only the stochastic part is corrected.\n", " (alpha = 0, eta = 1/2)\n", " -Code: 'pred-corr', 'predictor-corrector', 'pc-euler'\n", " Both the deterministic and stochastic part corrected.\n", " (alpha = 1/2, eta = 1/2)\n", " -Code: 'pc-euler-imp', 'pc-euler-2', 'pred-corr-2'\n", " Numerical Solution of Stochastic Differential Equations\n", " Chapter 15.5 Eq. (5.4), By Peter E. Kloeden, Eckhard Platen\n", " \n", " platen:\n", " Explicit scheme, create the milstein using finite difference instead of\n", " derivatives. Also contain some higher order terms, thus converge better\n", " than milstein while staying strong order 1.0.\n", " Do not require derivatives, therefore usable for\n", " :func:qutip.stochastic.general_stochastic\n", " -Order strong 1.0, weak 2.0\n", " -Code: 'platen', 'platen1', 'explicit1'\n", " The Theory of Open Quantum Systems\n", " Chapter 7 Eq. (7.47), H.-P Breuer, F. Petruccione\n", " \n", " rouchon:\n", " Scheme keeping the positivity of the density matrix. (smesolve only)\n", " -Order strong 1.0?\n", " -Code: 'rouchon', 'Rouchon'\n", " Efficient Quantum Filtering for Quantum Feedback Control\n", " Pierre Rouchon, Jason F. Ralph\n", " arXiv:1410.5345 [quant-ph]\n", " \n", " taylor1.5, Order 1.5 strong Taylor scheme:\n", " Solver with more terms of the Ito-Taylor expansion.\n", " Default solver for smesolve and ssesolve.\n", " -Order strong 1.5\n", " -Code: 'taylor1.5', 'taylor15', 1.5, None\n", " Numerical Solution of Stochastic Differential Equations\n", " Chapter 10.4 Eq. (4.6), By Peter E. Kloeden, Eckhard Platen\n", " \n", " taylor1.5-imp, Order 1.5 implicit strong Taylor scheme:\n", " implicit Taylor 1.5 (alpha = 1/2, beta = doesn't matter)\n", " -Order strong 1.5\n", " -Code: 'taylor1.5-imp', 'taylor15-imp'\n", " Numerical Solution of Stochastic Differential Equations\n", " Chapter 12.2 Eq. (2.18), By Peter E. Kloeden, Eckhard Platen\n", " \n", " explicit1.5, Explicit Order 1.5 Strong Schemes:\n", " Reproduce the order 1.5 strong Taylor scheme using finite difference\n", " instead of derivatives. Slower than taylor15 but usable by\n", " :func:qutip.stochastic.general_stochastic\n", " -Order strong 1.5\n", " -Code: 'explicit1.5', 'explicit15', 'platen15'\n", " Numerical Solution of Stochastic Differential Equations\n", " Chapter 11.2 Eq. (2.13), By Peter E. Kloeden, Eckhard Platen\n", " \n", " taylor2.0, Order 2 strong Taylor scheme:\n", " Solver with more terms of the Stratonovich expansion.\n", " -Order strong 2.0\n", " -Code: 'taylor2.0', 'taylor20', 2.0\n", " Numerical Solution of Stochastic Differential Equations\n", " Chapter 10.5 Eq. (5.2), By Peter E. Kloeden, Eckhard Platen\n", " \n", " ---All solvers, except taylor2.0, are usable in both smesolve and ssesolve\n", " and for both heterodyne and homodyne. taylor2.0 only work for 1 stochastic\n", " operator not dependent of time with the homodyne method.\n", " The :func:qutip.stochastic.general_stochastic only accept derivatives\n", " free solvers: ['euler', 'platen', 'explicit1.5'].\n", " \n", " Available solver for photocurrent_sesolve and photocurrent_mesolve:\n", " Photocurrent use ordinary differential equations between\n", " stochastic \"jump/collapse\".\n", " euler:\n", " Euler method for ordinary differential equations.\n", " Default solver\n", " -Order 1.0\n", " -Code: 'euler'\n", " \n", " predictorâ€“corrector:\n", " predictorâ€“corrector method (PECE) for ordinary differential equations.\n", " -Order 2.0\n", " -Code: 'pred-corr'\n", "\n" ] } ], "source": [ "# list solver\n", "help(stochastic_solvers)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total run time: 1.09s\n" ] } ], "source": [ "# Euler-Maruyama\n", "sol_euler = smesolve(H, rho0, tlist, c_op, sc_op, e_op,\n", " nsubsteps=Nsub, method='homodyne', solver='euler-maruyama',\n", " options=Odeoptions(store_states=True, average_states=False))" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total run time: 1.55s\n" ] } ], "source": [ "# Milstein\n", "sol_mil = smesolve(H, rho0, tlist, c_op, sc_op, e_op,\n", " nsubsteps=Nsub, method='homodyne', solver='milstein',\n", " options=Odeoptions(store_states=True, average_states=False))" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total run time: 2.42s\n" ] } ], "source": [ "# predictor-corrector\n", "sol_pc_euler = smesolve(H, rho0, tlist, c_op, sc_op, e_op,\n", " nsubsteps=Nsub, method='homodyne', solver='pc-euler',\n", " options=Odeoptions(store_states=True, average_states=False))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total run time: 3.77s\n" ] } ], "source": [ "# predictor-corrector-2\n", "sol_pc_euler2 = smesolve(H, rho0, tlist, c_op, sc_op, e_op,\n", " nsubsteps=Nsub, method='homodyne', solver='pred-corr-2',\n", " options=Odeoptions(store_states=True, average_states=False))" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total run time: 7.01s\n" ] } ], "source": [ "# Default: taylor1.5\n", "sol_taylor15 = smesolve(H, rho0, tlist, c_op, sc_op, e_op,\n", " nsubsteps=Nsub, method='homodyne', solver='taylor1.5',\n", " options=Odeoptions(store_states=True, average_states=False))" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total run time: 8.04s\n" ] } ], "source": [ "# Taylor 2.0\n", "sol_taylor20 = smesolve(H, rho0, tlist, c_op, sc_op, e_op,\n", " nsubsteps=Nsub, method='homodyne', solver='taylor2.0',\n", " options=Odeoptions(store_states=True, average_states=False))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Variance $V_{\\mathrm{c}}$ as a function of time" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd8VFX6+PHPmUmH0DsBEnpCEgIkVAEVaYpREEXE/kXWgvWrX3GLuiq76vpTF5eVdRfBgoKKBVdULISmlACREkooAUJNAgnpycw8vz9mEgdIyBASkjjP+/WaF3PPnHvuucNknjn33PtcIyIopZRSltrugFJKqbpBA4JSSilAA4JSSikXDQhKKaUADQhKKaVcNCAopZQCNCAopZRy0YCglFIK0ICglFLKxae2O3AhWrRoIaGhobXdDaWUqlc2btyYISItK6tXrwJCaGgoiYmJtd0NpZSqV4wxBzypp4eMlFJKARoQlFJKuWhAUEopBdSzOQSlakNJSQlpaWkUFhbWdleUOq+AgABCQkLw9fWt0voaEJSqRFpaGsHBwYSGhmKMqe3uKFUuESEzM5O0tDTCwsKq1IZHh4yMMWOMMbuMMXuMMTPKef1OY0y6MSbJ9Zjq9todxpgU1+MOt/J+xpitrjZnGf1LU3VUYWEhzZs312Cg6jRjDM2bN7+okWylAcEYYwVmA2OBCGCyMSainKqLRCTG9fiPa91mwDPAAKA/8Iwxpqmr/pvAPUA312NMlfdCqRqmwUDVBxf7OfVkhNAf2CMi+0SkGFgIXOdh+6OB70TkpIicAr4Dxhhj2gKNRGStOO/h+S5wfRX675FJgwbzu9Gjaqp5pZT6TfAkILQHDrktp7nKznaDMWaLMeYTY0yHStZt73peWZsYY6YZYxKNMYnp6ekedPdcK5OTWZe8u0rrKlUXWK1WYmJiyh4vvvjieevPnz+f6dOnX6Leqd+K6ppU/hL4UESKjDG/A94BrqyOhkXkLeAtgNjYWKlKGxaLBbs4qqM7StWKwMBAkpKSaqx9m82Gj4+eY+LtPBkhHAY6uC2HuMrKiEimiBS5Fv8D9Ktk3cOu5xW2WZ2sxmB3aEBQvz2hoaFkZGQAkJiYyOWXX35OnfT0dG644Qbi4uKIi4tjzZo1ADz77LPcdtttDBkyhNtuu+1SdlvVUZ78JNgAdDPGhOH80r4ZuMW9gjGmrYgcdS3GAztcz78F/uI2kTwKeEpEThpjThtjBgLrgNuBNy5uVypmtVhxaEBQ1eDPX24n+cjpam0zol0jnrm213nrFBQUEBMTU7b81FNPMWnSJI/af/jhh3n00Ue57LLLOHjwIKNHj2bHDuefaHJyMqtXryYwMLDqO6B+MyoNCCJiM8ZMx/nlbgXeFpHtxpjngEQRWQI8ZIyJB2zASeBO17onjTHP4wwqAM+JyEnX8/uB+UAg8LXrUSMsFqOHjFS9djGHjL7//nuSk5PLlk+fPk1ubi4A8fHxGgxUGY8OGorIUmDpWWVPuz1/CniqgnXfBt4upzwRiLyQzlaV1Viw6QhBVYPKfslfaj4+PmWj34rOP3c4HKxdu5aAgIBzXmvQoEGN9k/VL16Ry8hiseghI/WbFBoaysaNGwFYvHhxuXVGjRrFG2/8ekS2JienVf3mFQHBqmcZqXqudA6h9DFjhjNhwDPPPMPDDz9MbGwsVqu13HVnzZpFYmIi0dHRREREMGfOnEvZdVWPGOd1YfVDbGysVOUGOVEhHcnMzeFI1qka6JX6rduxYwfh4eG13Q2lPFLe59UYs1FEYitb1ytGCD4WCw4dISil1Hl5RUCwGCsOR/0ZCSmlVG3wioBg1dNOlVKqUl4SECx6pbJSSlXCKwKCxRidQ1BKqUp4RUCwWqzYdQ5BKaXOy0sCgp5lpOo3Ywy33npr2bLNZqNly5aMGzcOgCVLlpSlxH722Wd55ZVXKmwrISGBn376qdJturepvINX5Lu1WoyOEFS91qBBA7Zt20ZBQQGBgYF89913tG//6y1E4uPjiY+P96ithIQEGjZsyODBg89b70LaVL8NXjVC0FNPVX129dVX89VXXwHw4YcfMnny5LLXKrohzqxZs4iIiCA6Opqbb76Z1NRU5syZw2uvvUZMTAyrVq2qMD22e5t33nknDz30EIMHD6Zz58588sknl2CP1aXmFSMEi8WC3SHkFxbSMEgzO6qL8PUMOLa1ettsEwVjKz80c/PNN/Pcc88xbtw4tmzZwt13382qVavOu86LL77I/v378ff3JysriyZNmnDvvffSsGFDHn/8cQBuueWWCtNjuzt69CirV69m586dxMfHM3HixKrtr6qzvCIg+FicA6Hc7NMaEFS9FR0dTWpqKh9++CFXX321x+tMmTKF66+/nuuvL/+25edLj+3u+uuvx2KxEBERwfHjx6u2E6pO84qAYHEFhJysk7Rp27qWe6PqNQ9+ydek+Ph4Hn/8cRISEsjMzKy0/ldffcXKlSv58ssvmTlzJlu3nju6OV96bHf+/v5lz+tTDjTlOa+ZQwA4nVX5H5BSddndd9/NM888Q1RUVKV1HQ4Hhw4d4oorruCll14iOzub3NxcgoODycnJKaun6bFVKe8ICK60wHlZ1XvrQ6UutZCQEB566CGP6trtdm699VaioqLo06cPDz30EE2aNOHaa6/ls88+K5tU1vTYqpRH6a+NMWOAv+O8heZ/RKTccbMx5gbgEyBORBKNMVOAJ9yqRAN9RSTJGJMAtAUKXK+NEpET5+tHVdNf3zxoCIvW/sTSRR8w9qbJla+glBtNf63qk4tJf13pHIIxxgrMBkYCacAGY8wSEUk+q14w8DCwrrRMRBYAC1yvRwGfi4j7eHSK61aaNcrHxzlCyM3VEYJSSlXEk0NG/YE9IrJPRIqBhcB15dR7HngJKP/GrjDZte4lVzqHUJyXXxubV0qpesGTgNAeOOS2nOYqK2OM6Qt0EJGvztPOJODDs8rmGWOSjDF/MsYYTzpcFVbXCKEgL6+mNqGUUvXeRU8qG2MswKvA/56nzgAgX0S2uRVPEZEoYKjrcVsF604zxiQaYxLT09Or1EcfH+eRseLCigYvSimlPAkIh4EObsshrrJSwUAkkGCMSQUGAkuMMe4TGDdz1uhARA67/s0BPsB5aOocIvKWiMSKSGzLli096O65fKzO3Sws0ICglFIV8SQgbAC6GWPCjDF+OL/cl5S+KCLZItJCREJFJBRYC8SXTha7RhA34TZ/YIzxMca0cD33BcYB7qOHalU6qVyiIwSllKpQpQFBRGzAdOBbYAfwkYhsN8Y8Z4zxJBXiMOCQiOxzK/MHvjXGbAGScI44/n3BvfeQtfSQUXFxTW1Cqd+M1NRUIiMja7sbqhZ4lLpCRJYCS88qe7qCupeftZyA8zCSe1ke0O8C+nlR/Px8ASgp0oCgVHWz2Wxl83SqfvOKK5VLP6w2HSGoeio1NZWePXsyZcoUwsPDmThxIvn5+WzYsIHBgwfTu3dv+vfvf0ZKilJ79+5lzJgx9OvXj6FDh7Jz507AmdLaPY11w4YNz1nXbrfzxBNPEBcXR3R0NP/6178A5z0Vhg4dSnx8PBERETW01+pS84qw7uvvB0BJia2We6Lqu5fWv8TOkzurtc2ezXryZP8nK623a9cu5s6dy5AhQ7j77rv5xz/+wZw5c1i0aBFxcXGcPn2awMBzs/lOmzaNOXPm0K1bN9atW8f999/Pjz/+6FHf5s6dS+PGjdmwYQNFRUUMGTKEUaNGAbBp0ya2bdtGWFjYhe2wqrO8IiD4+zsPGdlsGhBU/dWhQweGDBkCwK233srMmTNp27YtcXFxADRq1OicdXJzc/npp5+48cYby8qKioo83uayZcvYsmVL2UgiOzublJQU/Pz86N+/vwaD3xivCAi+fs4Rgq2kpJZ7ouo7T37J15Szr91s1KgRheWcOXfXXXexefNm2rVrx8KFC2nSpEm5GUx9fHxwOJz3Gnc4HOWedCEivPHGG4wePfqM8oSEBBo0aHAxu6PqIK+YQ/B35Xm32ey13BOlqu7gwYP8/PPPAHzwwQcMHDiQo0ePsmHDBgBycnKw2WzMmzePpKQkli5dSqNGjQgLC+Pjjz8GnF/wv/zyCwChoaFs3LgRgCVLllBSzg+m0aNH8+abb5a9tnv3bvL0iv/fLK8ICH6ugGDXgKDqsR49ejB79mzCw8M5deoUDz74IIsWLeLBBx+kd+/ejBw5stwRw4IFC5g7dy69e/emV69efPHFFwDcc889rFixgt69e/Pzzz+X+4t/6tSpRERE0LdvXyIjI/nd736nh15/wzxKf11XVDX99bsvv8AdT/6Ju64cwds/fF8DPVO/ZXUh/XVqairjxo1j27Yau35T/UZcTPprrxghBDYIApyn0CmllCqfVwSEINdQ2G531HJPlKqa0NBQHR2oGucVAaFBsPOCG7tDRwhKKVURrwgIwY2CAR0hKKXU+XhFQAhq1BgAmwYEpZSqkFcEhIaNmwOUXYSjlFLqXF4REBo0agqAXQOCUpWqjvTXCxYsIDo6mqioKAYPHlx2MZyq27widUXDJs0AHSEoVRPKS38dFhbGihUraNq0KV9//TXTpk1j3bp1tdRD5SmvGCH4BQVijE4qq/qrvqW/Hjx4ME2bOkfmAwcOJC0trVreB1WzvGKEAGA1FuyiAUFdnGN/+QtFO6o3/bV/eE/a/P73ldarr+mv586dy9ixYz3anqpdXhMQLBajcwiqXquP6a+XL1/O3LlzWb16tcfbVLXHo4BgjBkD/B2wAv8RkRcrqHcD8AkQJyKJxphQnPdh3uWqslZE7nXV7QfMBwJx3p7zYanBxEpWY8HuqD95m1Td5Mkv+ZpSl9Nfz549m3//23lb9KVLl9KuXTu2bNnC1KlT+frrr2nevPmF77C65CqdQzDGWIHZwFggAphsjDnnoKExJhh4GDh75miviMS4Hve6lb8J3AN0cz3GVG0XPGO1WHRSWdVrdTn99QMPPEBSUhJJSUm0a9eOgwcPMmHCBN577z26d+9e/W+GqhGeTCr3B/aIyD4RKQYWAteVU+954CXg3J8sZzHGtAUaicha16jgXeB6z7t94SzG6ByCqtfqU/rr5557jszMTO6//35iYmKIja000aaqC0TkvA9gIs7DRKXLtwH/OKtOX2Cx63kCEOt6HgrkAZuBFcBQV3ks8L3b+kOB/1aw/WlAIpDYsWNHqaqmQUHSP7RrlddX3is5Obm2uyD79++XXr161XY3VD1Q3ucVSJRKvutF5OInlY0xFuBV4M5yXj4KdBSRTNecwefGmF4X0r6IvAW8Bc77IVS1nzpCUEqp8/MkIBwGOrgth7jKSgUDkUCCa9KrDbDEGBMvIolAEYCIbDTG7AW6u9YPOU+b1U7nEFR9pumv1aXgyRzCBqCbMSbMGOMH3AwsKX1RRLJFpIWIhIpIKLAWiBfnWUYtXZPSGGM645w83iciR4HTxpiBxhlFbge+qN5dO5PFYsFej+4Op5RSl1qlIwQRsRljpgPf4jzt9G0R2W6MeQ7ncakl51l9GPCcMaYEcAD3ishJ12v38+tpp1+7HjXGedqp3g9BKaUq4tEcgogsxXmtgHvZ0xXUvdzt+WJgcQX1EnEearokrBYLDh0hKKVUhbwilxGAxVj0SmWllDoPrwkIVk1doeqprKws/vnPf1Zp3aqmsh4zZgxNmjRh3LhxFdaZP38+LVu2JCYmhpiYGP7zn/9UqY+q7vCagGAxFhx62qmqhy4mIFyo0ovOnnjiCd57771K60+aNKnsCuWpU6fWdPdUDfOagGC16CEjVT/NmDGDvXv3EhMTw6OPPsqIESPo27cvUVFRZVcdP/3007z++utl6/zhD3/g73//+xntFBYWctdddxEVFUWfPn1Yvnw54PylHx8fz5VXXsmIESMAGDFiBMHBwZdoD1Vd4TXZTnVSWVWHVR/tJuNQbrW22aJDQ4beVHG+nxdffJFt27aRlJSEzWYjPz+fRo0akZGRwcCBA4mPj+fuu+9mwoQJPPLIIzgcDhYuXMj69evPuD/C7NmzMcawdetWdu7cyahRo9i9ezfgTGW9ZcsWmjVrdkF9X7x4MStXrqR79+689tprdOjQofKVVJ3lNSMEnVRWvwUiwu9//3uio6O56qqrOHz4MMePHyc0NJTmzZuzefNmli1bRp8+fc7JMLp69WpuvfVWAHr27EmnTp3KAsLIkSMvOBhce+21pKamsmXLFkaOHMkdd9xRPTupao1XjRA0IKiLdb5f8pfCggULSE9PZ+PGjfj6+hIaGlqW0G7q1KnMnz+fY8eOcffdd19Qu+UltquMe8CZOnUq//d//3fBbai6xatGCDqprOqj4ODgskM/2dnZtGrVCl9fX5YvX86BAwfK6o0fP55vvvmGDRs2nHP/AoChQ4eyYMECwJnG+uDBg/To0aPK/Tp69GjZ8yVLlhAeHl7ltlTd4DUjBB/Xaacics6NRpSqy5o3b86QIUOIjIwkLi6OnTt3EhUVRWxsLD179iyr5+fnxxVXXEGTJk2wWq3ntHP//fdz3333ERUVhY+PD/Pnz8ff37/cbZbeezk3N5eQkBDmzp3L6NGjefrpp4mNjSU+Pp5Zs2axZMkSfHx8aNasGfPnz6+pt0BdIkbq0URrbGysJCYmVmndQV26s+NoGhk5ufhYvWZgpKrBjh076sWvX4fDQd++ffn444/p1q1bbXdH1ZLyPq/GmI0iUulNKbzmm7F0DsGmt9FUv0HJycl07dqVESNGaDBQVeY1h4xKTzstsTsI8D13OK1UfRYREcG+fftquxuqnvOaEULpaac2u44QlFKqPF4TEHwsFuzioERPPVVKqXJ5TUCwWAwOh+gIQSmlKuA1AcFqteIQoaCgqLa7opRSdZL3BASL89qDnJysWu6JUhfmUqe/TkpKYtCgQfTq1Yvo6GgWLVpUbr2ioiImTZpE165dGTBgAKmpqVXqo6o7PAoIxpgxxphdxpg9xpgZ56l3gzFGjDGxruWRxpiNxpitrn+vdKub4GozyfVodfG7UzGrxXlmUd7JjJrcjFLV7lKnvw4KCuLdd99l+/btfPPNNzzyyCNkZZ37Q2ru3Lk0bdqUPXv28Oijj/Lkk09ekj6qmlNpQDDGWIHZwFggAphsjIkop14w8DCwzq04A7hWRKKAO4CzE6xPEZEY1+NEFffBI6UjhNxsHSGo+uVSp7/u3r172bUM7dq1o1WrVqSnp5/Try+++KIsod3EiRP54YcfqE8XuqpzeXIdQn9gj4jsAzDGLASuA5LPqvc88BLwRGmBiGx2e307EGiM8ReRS34gv/RS/nw9ZKQuwvL5b3HiQPWe79+qU2euuHNaha/XZvrr9evXU1xcTJcuXc7p1+HDh8vSXfv4+NC4cWMyMzNp0aJFdbwtqhZ4EhDaA4fcltOAAe4VjDF9gQ4i8pUx5gnKdwOw6axgMM8YYwcWAy9IDf68KE1XkXc6p5KaStVdpemvV65cicViKTf99fHjx8vSX7sHhNWrV/Pggw8CnqW/Pnr0KLfddhvvvPMOFovXTDd6tYu+UtkYYwFeBe48T51eOEcPo9yKp4jIYdehpsXAbcC75aw7DZgG0LFjxyr308c1QijIPV3lNpQ63y/5S+FSpb8+ffo011xzDTNnzmTgwIHlrtO+fXsOHTpESEgINpuN7Ozsc+7BoOoXT8L+YcD9NkghrrJSwUAkkGCMSQUGAkvcJpZDgM+A20Vkb+lKInLY9W8O8AHOQ1PnEJG3RCRWRGJbtmzp6X6dw+rj3NXCvPwqt6FUbbjU6a+Li4sZP348t99+OxMnTqywX/Hx8bzzzjsAfPLJJ1x55ZWaSbie82SEsAHoZowJwxkIbgZuKX1RRLKBsoOGxpgE4HERSTTGNAG+AmaIyBq3Oj5AExHJMMb4AuOA76thfyrk6xohFBUU1ORmlKp2lzr99UcffcTKlSvJzMwsS2k9f/58YmJizkh//T//8z/cdtttdO3alWbNmrFw4cIaew/UJSIilT6Aq4HdwF7gD66y54D4cuomALGu538E8oAkt0croAGwEdiCc7L574C1sn7069dPqureMaMFkOcf/98qt6G8U3Jycm13wSN2u1169+4tu3fvru2uqFpU3ucVSBQPvus9mkMQkaXA0rPKnq6g7uVuz18AXqig2X6ebLu6+LkynBYXFl/KzSp1SSQnJzNu3DjGjx+v6a9VlXlN+mtfH18ASgr1kJH67dH016o6eM25ZEFBzmOlxQU6qayUUuXxmoAQ2CAIgOKCwlruiVJK1U1eExAaBjvPsy4p0mynSilVHq8JCMGNGgFQXFRSyz1RSqm6yWsCQuPmzsvyS4r1LCOl3CUkJPDTTz95XD81NZUPPvig2rb/7LPP8sorr1xUG57uw4Xua1U0bNiw0jqvv/46+fm/zmdeffXV5WaUvdS8JiA0dSXcKimx1XJPlKpbajsgVIe6FBA8cXZAWLp0KU2aNKnFHjl5T0Bo0w6AEpseMlL1z/vvv0///v2JiYnhd7/7HXa7nQMHDtCtWzcyMjJwOBwMHTqUZcuWAXD99dfTr18/evXqxVtvvVXWzjfffEPfvn3p3bs3I0aMIDU1lTlz5vDaa68RExPDqlWrztjuihUriImJISYmhj59+pCTk8OMGTNYtWoVMTExvPbaa9jtdp544gni4uKIjo7mX//6FwC5ubnlpuoGmDlzJt27d+eyyy5j165dAOzdu5e+ffuW1UlJSTljudSsWbOIiIggOjqam2++udx9+PLLLxkwYAB9+vThqquu4vjx4+XWS09P54YbbiAuLo64uDjWrFlzzva2b99e9t5HR0eTkpICwKuvvkpkZCSRkZFnpB4vlZCQwLhx48qWp0+fzvz585k1axZHjhzhiiuu4IorrgAgNDSUjIyMCttNTU0lPDyce+65h169ejFq1CgKaiLrgidXr9WVx8VcqXw4ZYcAcm1M1dtQ3unsKz+HDx8u8+bNExGR4uJiGT58uLz33nsiIpKXlyfDhw+XhQsXiohIVlaWDB8+XBYvXiwiIunp6TJ8+HBZsmSJiIgcPXrUo+2PGzdOiouLRUTkvvvuk3feeUdERP7973/LxIkT5eWXX5Zp06aVrZOZmSkiIvn5+dKrVy/JyMiQEydOSEhIiOzbt++MOs8884z87W9/K3fb48aNk9WrV4uISE5OjpSUlMjy5cvlmmuuKavzr3/9S55//nkRESksLJR+/frJvn37pKSkRLKzs8v2u0uXLuJwOCQxMVEiIyMlLy9PsrOzpUuXLmXbv/zyy2Xz5s0iIvLUU0/JrFmzzulT27ZtpbCwUERETp06Ve4+nDx5UhwOR9l79Nhjj5Vbb/LkybJq1SoRETlw4ID07NnznO1Nnz5d3n//fRERKSoqkvz8/LJ9yM3NlZycHImIiJBNmzaJiEiDBg1ERM55nx544IGyz02nTp0kPT297LXS5Yra3b9/v1it1rL35sYbbyz7zJ2txq9U/i1o2SEMgBK7vZZ7otSF+eGHH9i4cSNxcXEAFBQU0KqV8waDU6dO5eOPP2bOnDkkJSWVrTNr1iw+++wzAA4dOkRKSgrp6ekMGzaMsDDn38LZ6a7LM2TIEB577DGmTJnChAkTCAkJOafOsmXL2LJlC5988gngTMCXkpJCSEhIuam6V61axfjx4wkKcp4KHh8fX9bW1KlTmTdvHq+++iqLFi1i/fr152wvOjqaKVOmcP3113P99deX2++0tDQmTZrE0aNHKS4uLtvns33//fckJ/96a5fTp0+Tm5t7xjzAoEGDmDlzJmlpaUyYMIFu3bqxevVqxo8fX5YldsKECaxatYo+ffpU9paeV0XtxsfHExYWRkxMDAD9+vWrkVuWek1A8PX3x8diocSmcwjq4iQkJJQ99/X1PWM5KCjojOXGjRufsdyiRYszltu0aVPp9kSEO+64g7/+9a/nvJafn09aWhrgPEQTHBxMQkIC33//PT///DNBQUFcfvnlZSmyL9SMGTO45pprWLp0KUOGDOHbb78tt39vvPHGORlW58+fX2Gq7orccMMN/PnPf+bKK6+kX79+5abT/uqrr1i5ciVffvklM2fOZOvWrefUefDBB3nssceIj48nISGBZ599ttztORwO1q5dS0BAQIV9uuWWWxgwYABfffUVV199ddkhscr4+PjgcDjKlqv6f1DKPRGh1WqtkUNGXjOHAM6MpyV2R+UVlapDRowYwSeffMKJE867zJ48ebIs7fWTTz7JlClTeO6557jnnnsA5y/0pk2bEhQUxM6dO1m7di0AAwcOZOXKlezfv7+sHTgzvfbZ9u7dS1RUFE8++WRZptWz648ePZo333yTkhLn/Nzu3bvJy8urMFX3sGHD+PzzzykoKCAnJ4cvv/yyrK2AgABGjx7Nfffdx1133XVOfxwOB4cOHeKKK67gpZdeIjs7uywQuvcpOzub9u3bA5Sl6C5vX0eNGsUbb7xRtuw+yiq1b98+OnfuzEMPPcR1113Hli1bGDp0KJ9//jn5+fnk5eXx2WefMXTo0DPW69SpE8nJyRQVFZGVlcUPP/xQYT9KedJuTfK+gKAjBFXPRERE8MILLzBq1Ciio6MZOXIkR48eZcWKFWzYsKEsKPj5+TFv3jzGjBmDzWYjPDycGTNmlN3gpmXLlrz11ltMmDCB3r17M2nSJACuvfZaPvvss3InlV9//XUiIyOJjo7G19eXsWPHEh0djdVqpXfv3rz22mtMnTqViIgI+vbtS2RkJL/73e+w2WxMmTKFxMREoqKiePfdd8tSdfft25dJkybRu3dvxo4dW3YorNSUKVOwWCyMGjWKs9ntdm699day+0I/9NBDNGnS5Jx9ePbZZ7nxxhvp16/fGbf0PLverFmzSExMJDo6moiICObMmXPONj/66CMiIyOJiYlh27Zt3H777fTt25c777yT/v37M2DAAKZOnXrO4aIOHTpw0003ERkZyU033XTG69OmTWPMmDFlk8qlPGm3JhmpRzfFjo2NlcTExCqv36JhQ0JbtCYxdW/llZVy2bFjB+Hh4bXdDa/xyisNbjzUAAAgAElEQVSvkJ2dzfPPP1/bXamXyvu8GmM2ikhsZet6zRwCOG+jqZPKStVd48ePZ+/evfz444+13RWv5FUBwdfqg81uw2Z34GP1qqNlStULpWdGqdrhVd+KvlYfSux2Cm06sayUUmfzKCAYY8YYY3YZY/YYY2acp94NxhgxxsS6lT3lWm+XMWa0W7lHbVYnP6uVYruNwhI9bKSUUmerNCAYY6zAbGAsEAFMNsZElFMvGHgYWOdWFgHcDPQCxgD/NMZYPW2zuvn6OEcIBcUaEJRS6myejBD6A3tEZJ+IFAMLgevKqfc88BLgfvXFdcBCESkSkf3AHld7nrZZrXwtVkrsNopsGhCUUupsngSE9sAht+U0V1kZY0xfoIOIfOXhupW2WRN8fXyw2e0UlugcglKlqisDaGpqKpGRkZXW+8tf/nLG8uDBgy9626p6XPSksjHGArwK/O/Fd6fc9qcZYxKNMYnp6ekX1ZZf6aSyziEoVeZSp4Q+OyDUhXTUysmTgHAY6OC2HOIqKxUMRAIJxphUYCCwxDWxXNG6lbVZRkTeEpFYEYlt2bKlB92tmK+P8zqEAg0Iqp6pS+mvRYQnnniCyMhIoqKiWLRo0Tn9nT9/PtOnTy9bHjduHAkJCcyYMYOCggJiYmKYMmUK8OsNZSpqNyEhgcsvv5yJEyfSs2dPpkyZQn26oLY+8eQ6hA1AN2NMGM4v7ZuBW0pfFJFsoOzacGNMAvC4iCQaYwqAD4wxrwLtgG7AesCcr82aUjapXKTpK1TVPPLII+Xmu7kYMTEx5ebTL7Vjxw4WLVrEmjVr8PX15f7772fBggXcfvvtPPnkk9x3333079+fiIiIsnQPb7/9Ns2aNaOgoIC4uDhuuOEGHA4H99xzDytXriQsLIyTJ0/SrFkz7r33Xho2bMjjjz9+zrZfeeUVZs+ezZAhQ8jNzSUgIIBPP/2UpKQkfvnlFzIyMoiLi2PYsGEe7euLL77IP/7xj3Lfw/O1u3nzZrZv3067du0YMmQIa9as4bLLLvNom8pzlY4QRMQGTAe+BXYAH4nIdmPMc8aY+ErW3Q58BCQD3wAPiIi9ojYvblcq5+/jg0OErJMnanpTSlUb9/TXMTEx/PDDD+zbtw9wpos+ffo0c+bMOeM2lLNmzaJ3794MHDiwLP312rVrq5z+etasWWRlZeHj48Pq1auZPHkyVquV1q1bM3z4cDZs2HDR+3m+dvv3709ISAgWi4WYmJgaSf2sPLxSWUSWAkvPKnu6grqXn7U8E5jpSZs1zc/XubunjxwEel/KTavfiPP9kq8pdT39dXlqOvWzTZNU1givulLZ3xUQ8tKP1XJPlPJcXUt/PXToUBYtWoTdbic9PZ2VK1fSv3//M9YLDQ0lKSmpLF21+41ufH19y1Jlu/OkXVWzvCqXkZ+fLwD5eshI1SPu6a8dDge+vr7Mnj2b1NRUNmzYwJo1a7BarSxevJh58+Zxyy23MGfOHMLDw+nRo0e56a8dDgetWrXiu+++49prr2XixIl88cUXvPHGG2fk33/99ddZvnw5FouFXr16MXbsWPz8/Pj555/p3bs3xhhefvll2rRpc8ZhnCFDhhAWFkZERATh4eFn3Bt52rRpREdH07dvXxYsWFBWPn78+HLb3blzZ82/yQrwsvTX940ZxZxvv+OZx57g2f/3cjX2TP2WafprVZ9cTPprrzpkFBTgPA5ZlHu6lnuilFJ1j3cFhCDnfVOL8vJruSdKKVX3eFVACAgKBKC4QAOCujD16dCq8l4X+zn1qoDQMLgBAMWFRbXcE1WfBAQEkJmZqUFB1WkiQmZmJgEBAVVuw6vOMmrYqBEAtuLiWu6Jqk9CQkJIS0vjYnNpKVXTAgICCAkJqfL6XhUQGjdtCkBJiQYE5TlfX9+yq3uV+i3zqkNGTVo4Uy6VlOhVjkopdTavCghNW7cCNCAopVR5vCogNG/XEQCb3jFNKaXO4VUBoW2nLgAUa2IspZQ6h1cFBP+GwVgtFkrsGhCUUupsXhUQAHytFkp0hKCUUufwwoBgpcRmp9jmqLyyUkp5ES8MCD7YHHbyi3WUoJRS7rwwIFgpttvIK9YzjZRSyp1HAcEYM8YYs8sYs8cYM6Oc1+81xmw1xiQZY1YbYyJc5VNcZaUPhzEmxvVagqvN0tdaVe+ulc/HYqXEbie/SEcISinlrtLUFcYYKzAbGAmkARuMMUtEJNmt2gciMsdVPx54FRgjIguABa7yKOBzEUlyW2+KiFT9jjdV4OfjS7GthFwNCEopdQZPRgj9gT0isk9EioGFwHXuFUTE/Y4zDYDy0kJOdq1bq4L8/CgoLiZfDxkppdQZPAkI7YFDbstprrIzGGMeMMbsBV4GHiqnnUnAh2eVzXMdLvqTMcZ42OeLEujnT0FJMXk6QlBKqTNU26SyiMwWkS7Ak8Af3V8zxgwA8kVkm1vxFBGJAoa6HreV164xZpoxJtEYk1gd6Ycb+PvrCEEppcrhSUA4DHRwWw5xlVVkIXD9WWU3c9boQEQOu/7NAT7AeWjqHCLylojEikhsy5YtPeju+TXw96ewpIRTp7Mvui2llPot8SQgbAC6GWPCjDF+OL/cl7hXMMZ0c1u8Bkhxe80C3ITb/IExxscY08L13BcYB7iPHmpMgwB/BMjak1xpXaWU8iaVnmUkIjZjzHTgW8AKvC0i240xzwGJIrIEmG6MuQooAU4Bd7g1MQw4JCL73Mr8gW9dwcAKfA/8u1r2qBLBDZy3l8s5sONSbE4ppeoNj+6YJiJLgaVnlT3t9vzh86ybAAw8qywP6HchHa0ujYMbApB77EhtbF4ppeosr7tSuVmzxgDkZ2XWck+UUqpu8bqA0KK1c2K6MCe3lnuilFJ1i9cFhHYdnXdNK8wrqOWeKKVU3eJ1ASEsPAKA/EINCEop5c7rAkLn3s657MKi4lruiVJK1S1eFxAaNGuBv48PBRoQlFLqDF4XEAAC/fzILyqq7W4opVSd4p0BwdeP/BINCEop5c47A4IrBbZSSqlfeWVACHKlwC62OWq7K0opVWd4Z0Dw9aOguIgCTYGtlFJlvDIgOG+SU0Jesd4kRymlSnllQAgK8KeopITMrFO13RWllKozvDIglN4TYf+WTbXdFaWUqjO8MiA0DHTeE+Hw7u213BOllKo7vDIgNG7UAIDjqQdruSdKKVV3eGVAaNmiKQCnjp+o5Z4opVTd4VFAMMaMMcbsMsbsMcbMKOf1e40xW40xScaY1caYCFd5qDGmwFWeZIyZ47ZOP9c6e4wxs4wxpvp26/zCI3sCcDI941JtUiml6rxKA4IxxgrMBsYCEcDk0i98Nx+ISJSIxAAvA6+6vbZXRGJcj3vdyt8E7gG6uR5jLmI/LsigsfEAZGXnXKpNKqVUnefJCKE/sEdE9olIMbAQuM69goicdltsAMj5GjTGtAUaichaERHgXeD6C+r5RWjbM5IgPz9O5epd05RSqpQnAaE9cMhtOc1VdgZjzAPGmL04RwgPub0UZozZbIxZYYwZ6tZmWmVt1hhjaBLUgKz8vEu2SaWUquuqbVJZRGaLSBfgSeCPruKjQEcR6QM8BnxgjGl0Ie0aY6YZYxKNMYnp6enV1V2aBjXgVH4ehSWavkIppcCzgHAY6OC2HOIqq8hCXId/RKRIRDJdzzcCe4HurvVDPGlTRN4SkVgRiW3ZsqUH3fVMU9cI4Xh2YbW1qZRS9ZknAWED0M0YE2aM8QNuBpa4VzDGdHNbvAZIcZW3dE1KY4zpjHPyeJ+IHAVOG2MGus4uuh344qL35gI0Cw4mv7iYrb9svpSbVUqpOqvSgCAiNmA68C2wA/hIRLYbY54zxsS7qk03xmw3xiThPDR0h6t8GLDFVf4JcK+InHS9dj/wH2APzpHD19W1U55o3awxAMnLl17KzSqlVJ3l40klEVkKLD2r7Gm35w9XsN5iYHEFryUCkR73tJp17NAWgCO7UmqrC0opVad45ZXKAJF9owHIOKEXpymlFHhxQBg8znnZQ1b26UpqKqWUd/DagNAqrBsN/f314jSllHLx2oAA0DgwiFP5edjsem9lpZTy6oDQomEjTuRkk5qpVywrpZRXB4TQli3JzM1l5XfLarsrSilV67w6IET1DAMg6fOFtdwTpZSqfV4dEK67fTIAh/brndOUUsqrA0Ls6OtpHBjIQb1RjlJKeXdAwBg6NGvBoVPpnMwrru3eKKVUrfLugACEtmpFek4Oy7/5tra7opRStcrrA0JE904ArPnw3VruiVJK1S6vDwh3Pjwdq8Xwy7Ydtd0Vpcp1IDOPvCJbbXdDeQGvDwjhg66gZ9v2bDy4j7yTmbXdHaXOUGSzM27Wap5dsr3COjmFJZewR+q3zOsDAsCAXuFkFxTwyiMPnVF+6GQ+dofUUq/UpbL54CmWbj2Kow7+X29MPUVOkY0lvxwhO//cL/4vkg4T89x37DymSRrVxdOAANz2h6fxtVr54edNZWWHTuZz5f9L4O8/eHa/BJG692XyW5JXZCMzt6ja2xURHl6YxP0LNjHhzZ9Iz6n+bVyMVXsyMAaKbA4+T3LeZfZwVgF/+GwrW9KyeP6/ydgdwvfJx2u5p+q3QAMCMGTIYKJDOrE+dQ/L3n4TgI83plFiF95evZ/UjDzueTeRVSnp56xbWGLn6S+2EfvC9xzT+zOfV16Rjac+3crwvy0v99fu+Tz04WYmvbW22vu0fv9JDp7MZ0Lf9mxJy+Ldn1OrfRsXY1VKOnGhzYhq35gP1x9ERJi/Zj8L1h0k/h9rOJlXTJtGAazYfe5nUzk5HML/fvQLHyUequ2u1HkaEABfq4W46yZidzh46s9/w2Z38EniIbq2akhukY1rZq3iu+Tj/OPHPees+z/vbODdnw+QmVfM19uOlpV/nHiITzelXcrdOENmbhH7M+pW0r7J/17Lh+sPciAzn2XJxzxeL+V4Dj/sPMGeE7mkncq/4O2W2B2s23fm/NBjHyUx9Z1E3lt7gIb+PrxwfSSXdWvJp5sO15lDR5m5RWw/cpph3VowuX9Hdh7LYe2+kyzdeoy40KaMjWzD/47qwY2xIWw6mEV2Qd2eSyixO2ols/AnG9NYvCmND9drRoLKeBQQjDFjjDG7jDF7jDEzynn9XmPMVmNMkjFmtTEmwlU+0hiz0fXaRmPMlW7rJLjaTHI9WlXfbl24kXc8yIjIaDYd3M+Do8dwJLuQR67qxphebcgvsTO4S3PW7T/JkayCsnWOZBWwZk8mj17Vne6tG/LNNueX3LHsQv7w2TYe//gX1u6rnYnqpz7dylWvrmD28j3V8gUnIhd1WCwzt4gtadk8Pqo77ZsE8vU2zwPC22v2YzHO5z/vvfD389+r9jHprbVsO5wNwM5jp/l002G+33Gc/245yjVRbQny82FivxAOZxXU2v/Z2ValZCACl3VryYS+7WnWwI/ff7aVw1kFTIrryJu39uOBK7oyvHtL7A5hzZ7au+J++c4TzF29v9zXCkvsXPn/Euj2h68Z+vLySzqSPpVXzF+/3oExsDUtm4Ji+yXbdn1UaUAwxliB2cBYIAKYXPqF7+YDEYkSkRjgZeBVV3kGcK2IRAF3AO+dtd4UEYlxPU5czI5crAl9Qxj4yPN0at6Ct378nssSXuSyTg145abe/PfBy/jrhCgAlvxypGydhF3OYfrVUW0Y3asNG1JPkplbxL9X7cMuQrsmgTyyMIns/BJK7A4+XH+QDA+Pg5d++W5Jy+LpL7bx7JLtfLopzaMPdIndwZo9GTQK8OFv3+7i7TXl/6FeiPve38RjH/1S5fV3Hc8BoHeHJlwd1YZVKelkF5QgIjz+8S+8XcGXSWZuEZ9uOsxNsR1o1sCPn/dlkph6kvke7pPN7uD9nw8AkLDL+RH7z6r9BPpaee66XrRo6M9tg5zXooyKaE1wgA+fbKy9kV2p04Ul/O3bXYS1aEBU+8YE+Fq5Y1Ao+zPy8LUaRka0Lqsb06EJwQE+rNhVfYeN/rvlCBsPnPK4/hs/pvD8f5PLDaZf/nKEfel53DGoE6fyi/nTF9s8/nFxsWdQzVuzn6yCEp4Y3QObQ0g6lHXG63aHkJVffMlGhdn5JTzwwSae/mLbJdnehfJkhNAf2CMi+0SkGFgIXOdeQUTcT3FoAIirfLOIlH6DbgcCjTH+F9/tmvHH26/m6X/Oo3PLVixYt5q47j34/firKTq4hg7NAunTsQmLNhxi+a4TFNscJOw6QfsmgXRt1ZDRvdrgEHhl2W4+WHeQ63q3459T+nLsdCHvrU3l001pPPXpVia++ROHTlZ82OOjxEOMeX0lPf/0DbfNXceEf/7Ex4lpLN6YxmMf/cJlL/3IwczzHzb55VAWecV2/jI+iqHdWjB7+Z6L+sMSEdbszWBVSnqVRwm7jzkDQo82wYyNakuJXfh2+zE+2ZjGJxvT+HLLkXLX+8/q/RTbHUwd2plBnZuzKiWDe9/fyLNfJnt0SO675OMcyS4kwNfCit3pnDhdyBdJh7kxNoTbB4Wy4Q8jiGzfGIAAXysT+4XwWdJhVqfUbn6rZ7/YzrHThbx6U2+sruHRbYM6EeBr4bKuLWgc6FtW18dqYWDn5mxIPVkt2/52+zGmf7CZm9/6mcUb08gvPv81EIUldra6Rl9/+nwbxbZfDwuJCO/8nEq3Vg15Nr4Xj17Vne+Sj5eNps9nwboDRP95GbfNXccu1+fnQtgdwicb0xjarSVT+nfCGEh0e4/yimxcN3s1Mc99R+Sz35Jy/MK3cSF+2pPBtf9YzVdbjrJow6Ez3qe6wlT2B26MmQiMEZGpruXbgAEiMv2seg8AjwF+wJUiklJOO/eKyFWu5QSgOWAHFgMvSDmdMcZMA6YBdOzYsd+BAweqsJtw+eWXc+edd3LnnXdSUlLCyJEjmTp1Krfeeiv5+flcffXV3HfffUyaNInD+1Lo2zeO0/n5FJaUYDEGq8VC60aNadu4KX4+PiQfSaNL67a0bdIUu8POhn176NkuhNaNm5JXVMDG/XvpExZK46BgTufls3H/HqI6dKJ5cCNOF+STdGAfsWGdaR7ckJO5uazft4+4zp0J9A8kKzeXLQdS6RvamWYNGpJVkMum1P0M6NodP19/MnNOsz3tAHFdutG8YQOOZJ1iy4EDDAsPp4G/P6npmWw/dJArIiII9PVj3/F0dh8/zIjwXvj7+5OamU5yWhpXRUXh5+PLvhPH2Xn4MKOie+NjsbDn+HF2Hz3C6OgYrBbD7qNH2X30KFf0iibI34fdRw6z/0Q6Y6J7IxbYefgwBzMzubJXJA4g5dhRjp46yVVRUSCw7dAhjmdlMahnOEG+FrYcPEhGTi4Du/XACGxPO0ROQQHDw3tggKQDB8krKmJojx4UlDjYcvAANnsJA7t0o8gOmw7uxyEO+nTuSpCvhQ0pKYAwoHMYBuGnvfvxsVqJ69KVwhIH6/fsJsDXl14dQvGxGH7evYtGQQH06+QcGazYuYvGgYHEdOyIAD8m76BJgwbEhHYEICE5mZaNGtGrQwgG+GHbdto2aUqvkBAQWLZtCx2bNye8bXsAvtnyC2EtW9KjbRtE4OstW+jWujVdW7fG4RC+3baF7q3b0LVVa2x2O8uSt9GzbVvCWrWh2F7CD9u206V1G8JatsRut7F8RzKR7UPo0Lw5eUXFrNy1g+iQDrRv3oy8wiJW7tpFZEgHWgQ3xmYrZnXKbvp06kSbJk3JLsjnp5TdxIV1pnXjxpzKy+OnlN3079KFFsGNyMjJZd3eFAZ06UZwYCAZOTlsSt1Hv7CuNG3QgOPZWWw9eID+XbrSMrghx7Kz2XwglWE9etIwMIgjJ0+SdCCV2M7daBwUyKHMTHYdTWN4eDiBfv4czMhge1oaV/SMICjAnwOZ6Ww/lMawHhEENwh0fvaOHGZEVBQ+Viv7jh0j5dgxRkRFUWSHgxnH2X/8OJdH9ybIx0rKkSOknjjB6N4xgJR99kZFRYNA8uE0jmadYkRUNA6HkHTgAKfycrgqKoq8Yjs70g6Sk5/H8PAIiuwOtqQdJK+wkJjwHgRarSTvTyWvsIjLIsIxwMY9eygsLmFIjx4AbNizB5vdzqBu3UFg/d4UEBjQpQsI/Lx3D1arhbguXcHAml27CPDxpXfHTtjsDtbt3UODgAAiQjoS4Gth9a6dNAoMJKZjqPOzuCOZpg0bEt3R+dn7cft2WjVqRGTHjvx98cc0at6iSt9/xpiNIhJbacXSY8MVPYCJwH/clm8D/nGe+rcA75xV1gvYC3RxK2vv+jcYWAbcXllf+vXrJ1U1fPhwmTdvnoiIFBcXy/Dhw+W9994TEZG8vDwZPny4LFy4UEREsrKyZPjw4fLBe+/Li9OmylURUdLAz19aNAyWpkFB4mOxCM5RkD70oQ99XJLHsX17qvz9ByRW9v0qIh6NEAYBz4rIaNfyUwAi8tcK6luAUyLS2LUcAvwI3CUiaypY504gVs4adZwtNjZWEhMTz9vfS8Lh4IPVu/h67U7+PKYT9uIcTp8+TW5eFtgEq1ixFxZid9goscD76w9TIHDXoE6A4a0Ve+nSPJj43u2wABYMBli0/jB5BUXcNawzJcZBEXYQwSLgsNsxdkHEwZKko/j5+nAqv4QgX18m9WuPlJTw4+5M9mQUkGdzENu5BcPDW+GwGE4V2vjX6r10bBbIxF6tMEUlYCwU223kFRbRLMCXI6cLWJd6iq6tgwlv1wirxcr76w/SvGEgzRsG8NP+TCwCAzo2ZkivEE4W2vnv5kOM7NoCxMGCdYfp1jYYP6th98GTjO/TgY4tg8HHB2Ox8Mr3KfRuE8xV4e1AwPj6YAyAg5STxcxbe4SHRkUR1qIBxVJMYVEhr36zlY7Bhht7t8QENQaLD6fSj9IYG6aohKXJJzhl9+XWyyPBxx+sVozFgMPB1v0n+CYplUGdGjCsextsVh/e/mErVoswZUhXAgMbgNUHY7WAxYIBEAERHLZiTpzO48vEw2QXQHRYR0aGt8LX2BGHDQfOfz/bfIDOLYKI6tCUQpsdq68FY6wIsP1YDl9tOU6jAH/aNQ5k14k8+oc1Y0xkGxwiWIwF5985lP7NS2Exh04W8P66A1we0YbBXVuC1QIY1wO3fwWxl4BDEGwU2YXXlu3msm6tyMkvJjU9h+4tGnIyv5gAf8PQ7o1c77dzc8b8+rkDw+LEw+SVOLhjUBhWq8FgMBaDwQdjgRPZxby//iAdmgYxsW9rikpysThKsDiEb3dkcLpAuHVwVyxWC7uOF7DneC6tAvxo2sSP5k39CQrwB5vzM2wPCORfP2whJNDG2MhOBPr542PxA9d7UlBs443vUxjWpQmxHZqw70QOX289woTYENo2CwYDYkCMce2Ica4rguD8G7HbhTdXpBDRugmjeramsKSAtOxiPks6iTirMqxzc+LaBGDLzmHZxgP4BPpzZXQ7wIax+IDVuD4jVrBaMFYft3+tzu1arOBrxbiOvgsCNgfYbEiJDWOx8G1KBusO5vDoyG408jGIvYS5K/bSLMiH8X07gHH9h5T+54iAOLBg8DM+2B12ihzFDJ94Jz5+flX6yvJ0hODjQVsbgG7GmDDgMHAzzlGA+8a6ya+HiK4BUlzlTYCvgBnuwcAY4wM0EZEMY4wvMA743oO+1A0WC7cMC+eWYeEeVY++pgCbXejQLAiADUHb+WDdQW7q159th7OZ/1MqN/QNYanvHv7nyjAuG3v+dsfc7vz39e938/cfUvhz/AhaNvTniZnf4xNqwSHC+DviiAppXLbOyZB9vPDVDpr7duDWoZ3wtVqY/sEm9mXksezRYcz7fBtrAzJxZMHwVs6zWlZuS8LPWBjQuhm5PvkE+VnZ2ziAP06I45Z/r+Nnn5ZYAlrRoqEf+zscZeEfRuBntXDZS8tZ17Ipd97aD3CejbV9049MvjqSKwd2Omd/OmfmM3Pfcujci56RbfGxGn45lMWO1YYHJvfhquh25b4P+1bu5S9Ld/KXUVfStnFgWXlhiZ3f/y2B1rHh/PG+wVhcx+DzehymaZAfw7q39Oj/7cqJRbyybBfvrT9EqybdeXBEt7LX9pzI4dNNK/E7ZeHVkb350+fbaNHQnzdv7UdYiwbMfHUFgb2iWTh9CD5WC48tSuLT7cfo2z6Sp7/Yzlu3xTKoS/Mztrc6JYO/LthIYJdQ7p9+OQ38Pfnz/NXLh1aQgGF3YQ6dOgWx/lQB7VsGciAzn04hPXjgiq7lrrcqJZ3vtq1n5uRIRg849/+nVF74fv78ZTJXde/N+D4hgPMc/8eeW8bV/dsy9IZoAIZ40NedjfvxyrLdLNlluKxrC343vDODuzgPh6zYnc7mbet54pYBDO7Sgq7ZBTzz1x8ZG9eLoYNDKSyx88wX2xkR3opBXZrz5S9HGd2rNc0b/jo9mXzkNBu3ruL2G2IYGtO+rLzPoSy+3naUiLaNuM6t/KuFm1mz/yR/fXgEALOX72HviVxenRTjwd7AtsPZvL1mP3+4OvyMfqTnFPHA35YzZkQb4m/8ta1v2MI324/x+o0jyz6fdUGlnzgRsRljpgPfAlbgbRHZbox5DucwZAkw3RhzFVACnMJ5RhHAdKAr8LQx5mlX2SggD/jWFQysOIPBv6txv+oU9y8rgMn9O/LOz6nc7LrQKqRpYNkV0SPDW5+9eoXGRrbl9e9TWLb9OEO7tSAjt5iZ4yOZUs4f9Z2DQ0nNzOOjDWks3OC8QKdxoC9+Vgv/98kWNh44xf+N6YGf1cILX+1g44FTNAnyJSu/hFUpGYzu1fr/t3fnYVXX+QLH359zWI7AAWTfZZNNRVBBETdocUnNsSZtMcuespua1b312PJUM7dl6rl3Jm/btFlzm7HRKcv0WmmpuUyj4ULkhisKCriAAols3/vH+YGgIosgefi+noeH8/txlu/nnN/5fb6/7wZWiyNrdhfzyQXJtGIAAA94SURBVOYj/HDgJAmB7qzeXYyT2cQtA0Nwt9g6OscnBrJw82HKKquxWhwbOgRjA6yXjCO4Zw+cHEwcOFHB3Qs2cbKiij5B7rg6mbn+Mu/HyBg/Xlqxm3W5x5mSEoZSimOnK1n04xEKz1TypylJTb5sjU8AreHt5szLkxPZdcw2D6JxQliz2zaix9nBxOyF2/B3d+ZURRUT39jAdfH+HDxRwdt3DsDBbKs53jc8giXbCnh0kW2k1uKsI00SQnFZJTM++pEIH1fenz6ozckAbKONFmfl42gWFs1Mw9fNGRF4+O/b+eOqXIZGeZMc1pNfqmp4dukO9h8v59MHh7I8+xjuFtuw28u5Oy2cJVsLeGP1voaEsD2/lDOVNQwK92pTWR8aFU1sgDtb8kr4bGs+d7y3icUz00iN8GJrXgkmgcQQTwAC3C14uzo1dFyv2lnEoqwjLMo6gquTmYqqWorLKnnk+piG588psI0m6hfs0eR1k0I9SQr1vKg80X5ufLH9KOXnajhztpr53+6lqraOmSOjmj1u61XX1vHY4u3kFpWz8+gZFt4/BC9XJ9buKebxT3+iplbxwMjIJo9JifBiUdYRcovLiAtwb9N715laNQ9BKbVCKRWjlIpSSr1o7HvWSAYopeYqpfoo2/DRDKXUDmP/C0opV3V+aGmSUqpYKVWhlBqolEo0HjdXKdVtBgjHBlhZ+cgIPrw3haWz0vn+8QzuGhJGYogHyWE9W/08Mf5uRPq6siz7/BDBgb0u/XgHs4kXJvVjw7wM5k9N4tnxCSyfM4y703qxJa8EJwcTU1PCmJEewdAob8rP1fD0uHiCPW3JLD7QnX7BHpysqOKpz3NIjfDi4/tSsTiaqKqtY1qjmv+E/kFU1dSxcodtOYWth0uM8l76i2U2CeHeLny7s4js/NPkl5zlmx1FjO4bgMXRfNn4G8/SXbDxEEP/sJr53+1lZIzvRTXw9sqI9SM7v7TJ0hlr9hQT62/lv3/bn/4hHnxy/xD+7+HhpEV6syz7aMPIs3p9gjzIjPMj3NuFGxL8+XZnEedqzh/yX+UUUlVbxxt3JDdcSbZVUqjts5/QPwh/dwsmkyAivPSbvvi6OfPMFz9TfKaS37z5Tz7dks+2w6VsP1LKhn0nSIvyxtmh+fcabJ/TbYNC2H+8omFEzvvrD2C1ODC6T+srMgAmk23o7LyxcXz/+Cj83Z159evdKKXYdqSUGH8rbkZSFBH6Bns0zCNZuv0o/u7OPJwZzdBoH/zdndl9rOkIoez801gtDoR7u7aqPNF+tmNzf3E5r6/eh0LhZDa1OJlNKcV76w+QW1TOzJGRHDxRwYTXN/DSil3c+9GPeLk48cWs9ItO+inhts9qa17ppZ62y+iZyl2kt7+VjFg/+od6YjYJL0zqx5ezhzUMMWwNEeGWASFsOniKJVsLsDo70Nvv8rUZP6uFm5OCmTEsglAvFx4YEYmbswOTk20Tn0wm4bWpSTw/IYFJycGM7x8IQFyAO+nRPlgtDswcEcn/zkjF282ZOZm9mZwcTELQ+QN+QJgnIT17sDT7KBXnavj4X3lkxvk1GSp5oShfNw6cqMAkMG9sHE5mW4JqKf6RMb6szz1BZXVtw4l4/tQk3rxzQKvfx5ZkxPmiFKwzli4pq6zmx0OnGBXny419Alg6exiRvm4EeFj44J4UFt4/mHenDbyoKeDPdw1k1WMjuXNwGGXnalife35o67Lso8QFWOndTNJsjVGxvvQNduehUU2bhqwWR566KZ4dR88w+rV1HD71C2/ckYzZJHy48SAFpWdJj27d6JXRfQIQgRU5hRw4Xs5XPxcybUgvrJbmP9uWuDg5MCezN1l5JazaWcT2wyUXVYz6BXuwt7icojOVfJ9bzMT+QTx2Yyzv3T2I5NCeDfNc6uXkn6ZfsEerm2N6+7sB8PWOQv6RdYQ7UsMY0zfgsnN/Nh88xeCXvuPVr/eQGefHvDFxLJqZhtkkvLvuAKMTAvhiVnqT70a9MC8XPHo4Nlz11DtbVcvktzbyn8t3Nix5XlunyO3kIbH1dEK4xt06MASzSdiw7wRJYZ5tSihgaxZZ+egInp/Yp2Gfn9XCPekROJpN3DW4F5lxfqRFehPt50bO86N5clx8Q819Vkb0Re2s9YlqXe5xpi/YTOkv1czJvHT7db1IX1tNLj3ahwdHRpH93I2kRrTcDDExKYiyczX89V95bD9SyqSkIG5OCm6oXXaEvkEe+Lg5NTQTfbOjiOpaRUbspSfXD43yIdLX7aL9Tg4mHM0m0o15BEu25VNXpygoPUtWXgkT+l+6r6S1gjx7sHzOcKL9Ln7tCYmBDI7wovxcDe9MG8j4xCAG9erJ8p9sy620NiH4uVtI6eXFl9kF/G7ZThzNJu5Nj7iicgNMSQkl0seVmX/dwpnKGpLDmjbr9A12p7ZO8eii7VTXqibNf7EBVg6drGg4cZ+rqWV34ZmGJqfW6OXlgqNZeHvtflydHZid2ZvbU8M4U1nDvCU/UX6J/0fx5+/3o4BXbunH67cnIyIkhXqyYu5wPrw3hbfuHEAPp0tfdYkIiSEe5BSUUlZZzW3v/MCWvBLW7T3O1sOlfLDhIGPnr2dvURmP/yObm9/Y2K5lW9pKJ4RrnL+7peHE1FxzUUuCPHs02zQT6uXCgntS8HBpWw1wdmY0E/oHkZVXwvDePi02hUX62E5i9SfF5r5IF0qL9CbUqwevfrMHgBsSAlp4RNuZTMKoWD+++vkY0z7YxBOfZhPl69ru99vRbOLWgSGsyClk3P+sZ9bfbKvsTmim87wjiAjvTx/EN4+MaOhUr++fCXC3EOnTuqYVgLH9Ath/vIIN+07w9Lh4fK1XPtfU0Wxi8YNpTE8LJ8zLhWEXJKiUcC/CvV34535b31WfRrXu+EArSsHe4jKUUuw6VkZ1rSIxxOPCl2mWg9lEhPEePHOTLaYhkV7Mva43y7KPMuWdH5rMZj59tpr1e48zKSmIKSlhTfp83JwdyIj1a/HqpF+wB3sKy/gqp5DNB0/x9tp9rNxRhEcPRxbeP5hfqmoYM389S7YV8NCoKEJ6tq8psS06rhqldZm7hoTx7a4i0iI7ps28IziaTbw2JYnUCC9GtWJUz/UJ/szOiG7zSdFkEqYMCuW/VuYS5uVCjP/FteOO8MSYWJwdTKzdc5zfDgzlmfHxOJrbX596elw8iSEevLvuAGerark7rRdh3p37hbdaHJs07WTG+/Hiil0MjfZGpPVXlpOTQ8gtKmdqSij9L9FB214+bs48P7FPk6vVet5uzqx9PIPK6lrMRt9IvVijfX5LXgkPfryFYmMJ8ws7lFsytm8g8YEVDZ3rIsKjN8QQ6GFh3pIcfjhwkqFR3lRW1/HtTttV4th+ge0Nl8QQD6prFW+ttS2auXp3Ma5ODtyQ4M/QKB8++7ehPLpoO2P6BvDAiKh2v05btDgP4dfkVzMP4Vco72QFvVrZgWZvCk9XMuyV1cwYFsFT41o3FFizdYi+/f1+MmL9iA/89Yx0aavaOkWf577G4mim9Jdqbk8NJcC9Bw9fF92mRNecyupahrz8HelRPrg6m/ky+yi+VmdqaxUb52W2+zUKSs+S/ofVAIyI8TWWhrH1NY3p27FXuh05D0G7BnTXZAAQ4GFh2ZxhhLVzdE53JSIXdUBfi8wmIcbfyk/5p0mP9ublyYkd+vwWRzO3DAhpWM01LsDK7sIyZo6IvKKEE+RhG057sqKKGenhCLDp4ElGxLRveYqOoBOCZheu5RquduVijYTQ3OS7K3V7aigLNh5kVIwvH0xPYf/x8nYPD64nIvQL8SDrUAlpUd7EBlg5WnoWF6euOy3rhKBp2jVvWlovwrxcOq0fLdrPyoqHhxPh44rJJFc0PLixJ8fGc7zsHM4OZgI9elw0ifVq0wlB07RrXmKIZ5uGmbZHZ1yFxgZYW5wJfTXpYaeapmkaoBOCpmmaZtAJQdM0TQN0QtA0TdMMOiFomqZpgE4ImqZpmkEnBE3TNA3QCUHTNE0zXFOL24nIcSCvnQ/3AU60eC/7omPuHrpbzN0tXrjymHsppVpcdviaSghXQkSyWrPanz3RMXcP3S3m7hYvXL2YdZORpmmaBuiEoGmaphm6U0J4t6sL0AV0zN1Dd4u5u8ULVynmbtOHoGmapl1ed7pC0DRN0y7D7hOCiIwRkT0isk9E5nV1eTqLiCwQkWIR+bnRPi8RWSUie43fPbuyjB1JREJFZI2I7BSRHSIy19hvzzFbRGSziGQbMf/O2B8hIpuMY3yRiDh1dVk7moiYRWSbiCw3tu06ZhE5JCI5IrJdRLKMfZ1+bNt1QhARM/AmMBZIAG4XkYSuLVWn+QgYc8G+ecB3SqnewHfGtr2oAf5dKZUADAFmGZ+tPcd8DshUSvUHkoAxIjIEeAX4k1IqGigB7uvCMnaWucCuRtvdIeYMpVRSo+GmnX5s23VCAFKBfUqpA0qpKuDvwM1dXKZOoZRaB5y6YPfNwF+M238BJl3VQnUipdQxpdRW43YZtpNFMPYds1JKlRubjsaPAjKBT439dhUzgIiEADcB7xvbgp3H3IxOP7btPSEEA0cabecb+7oLf6XUMeN2IeDflYXpLCISDiQDm7DzmI2mk+1AMbAK2A+UKqVqjLvY4zH+GvAEUGdse2P/MStgpYhsEZEHjH2dfmzr/6ncTSillIjY3ZAyEXEDPgMeUUqdsVUebewxZqVULZAkIp7A50BcFxepU4nIeKBYKbVFREZ1dXmuomFKqQIR8QNWicjuxn/srGPb3q8QCoDQRtshxr7uokhEAgGM38VdXJ4OJSKO2JLB35RSS4zddh1zPaVUKbAGSAM8RaS+cmdvx3g6MFFEDmFr8s0E5mPfMaOUKjB+F2NL/KlchWPb3hPCj0BvY0SCEzAV+LKLy3Q1fQlMN25PB5Z2YVk6lNGO/AGwSyn1x0Z/sueYfY0rA0SkB3ADtr6TNcCtxt3sKmal1JNKqRClVDi27+9qpdSd2HHMIuIqItb628CNwM9chWPb7iemicg4bG2QZmCBUurFLi5SpxCRT4BR2FZFLAKeA74AFgNh2FaJvU0pdWHH8zVJRIYB64EczrctP4WtH8FeY07E1ploxlaZW6yU+r2IRGKrPXsB24C7lFLnuq6kncNoMvoPpdR4e47ZiO1zY9MBWKiUelFEvOnkY9vuE4KmaZrWOvbeZKRpmqa1kk4ImqZpGqATgqZpmmbQCUHTNE0DdELQNE3TDDohaJqmaYBOCJqmaZpBJwRN0zQNgP8HuYTVEP5Nz+EAAAAASUVORK5CYII=\n", "text/plain": [ "