diff --git a/algorithms/discrete_poisson_solver/discrete_poisson_solver.ipynb b/algorithms/discrete_poisson_solver/discrete_poisson_solver.ipynb new file mode 100644 index 00000000..6ef8f9ba --- /dev/null +++ b/algorithms/discrete_poisson_solver/discrete_poisson_solver.ipynb @@ -0,0 +1,720 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Quantum Algorithm for Solving the Poisson Equation\n", + "\n", + "The Poisson equation is a partial differential equation that appears in various research areas, such as physics and engineering. It describes the distribution of a potential field, $u$, under some source term $b$:\n", + "$$\n", + "\\nabla^2 u = b,\n", + "$$\n", + "where the $\\nabla^2$ is the Laplacian (second derivatives) operator. One approach for numerically solving the Poisson equation is by moving from the continuos description to a discrete one, namely, by using the finite difference method that casts the problem into a set of linear equations. Then, the solution can be obtained by a linear solver." + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "In this notebook we treat the Poisson equation on a rectangular geometry, $L_x\\times L_y$, with Dirichlet boundary condition on the $x$ axis and a Neumann boundary condition on the $y$ axis:\n", + "$$\n", + "u(0) = u(L_x)=f_0,\\,\\,\\,\\, \\partial_y u|_{y=0} = \\partial_y u|_{y=L_y} = g_0.\n", + "$$\n", + "Furthermore, we will assume that $f_0=g_0=0$. The discritezation of space, including the treatment of the above boundary conditions, is given in Figure 1. The resulting linear equation reads:\n", + "$$\n", + "\\mathcal{L}\\cdot \\vec{u} = \\vec{b}, \\,\\,\\,\\,\\,\\,\\, \\mathcal{L} = \\mathcal{L}_{xx} \\otimes I_y + I_x \\otimes \\mathcal{L}_{yy},\n", + "$$\n", + "where $\\mathcal{L}_{xx}$ and $\\mathcal{L}_{yy}$ are the Laplacian operators [[1](#CST)]:\n", + "$$\n", + "\\mathcal{L}_{xx} = \\frac{1}{\\Delta x^2}\n", + "\\begin{pmatrix}\n", + "3 & -1 & 0 & \\cdots & 0 \\\\\n", + "-1 & 2 & -1 & \\cdots & 0 \\\\\n", + "0 & -1 & 2 & \\cdots & 0 \\\\\n", + "\\vdots & \\vdots & \\vdots & \\ddots & -1 \\\\\n", + "0 & 0 & \\cdots & -1 & 3 \\\\\n", + "\\end{pmatrix},\\,\\,\\,\n", + "\\mathcal{L}_{yy} = \\frac{1}{\\Delta y^2}\n", + "\\begin{pmatrix}\n", + "1 & -1 & 0 & \\cdots & 0 \\\\\n", + "-1 & 2 & -1 & \\cdots & 0 \\\\\n", + "0 & -1 & 2 & \\cdots & 0 \\\\\n", + "\\vdots & \\vdots & \\vdots & \\ddots & -1 \\\\\n", + "0 & 0 & \\cdots & -1 & 1 \\\\\n", + "\\end{pmatrix}\n", + "$$\n", + "and $\\Delta x$ and $\\Delta y$ are the discretization of the $x$ and $y$ axes, respectively. The above square matrices, which are of dimensions $N_x$ and $N_y$ respectively, represent the solution at the inner part of our geometry." + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "
\n", + "\n", + "
Figure 1. A schematic description of discrediting the Poisson equation. The area in which we would like to solve the problem is designated by the blue line. The linear equations are written only for the inner filled grid points. The unfilled ghost points are used to impose the boundary conditions. In the example above, a Dirichlet boundary condition is given by setting $u_{0,k}=-u_{-1,k}$, and $u_{3,k}=-u_{4,k}$, whereas a Neumann boundary condition reads $u_{j,0}=u_{j,-1}$, and $u_{j,3}=u_{j,4}$.
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "***\n", + "In this notebook we solve the Poisson problem with the [HHL](https://github.com/Classiq/classiq-library/blob/main/algorithms/hhl/hhl/hhl.ipynb) quantum linear solver. We utilize similar ideas appearing in Ref.[[2](#PoissonQuantum)], where a quantum cosine and a quantum sine transforms [[3](#QCST)] are performed towards achieving scalable implementation.\n", + "***" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "## How to Build the Algorithm with Classiq" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "The HHL algorithm essentially applies a matrix inversion. Here we treat the Laplacian matrix, which can be diagonalized by a quantum sine and cosine transforms, thus, the matrix we need to invert is a diagonal one. The main four quantum blocks of the algorithm are thus (see Figure 2):\n", + "1. Prepare the amplitudes of the source term on a quantum variable.\n", + "2. Perform QST and QCT at the beginning of the computation. This is done by applying the QST to the x qubits and the QCT to the y qubits.\n", + "3. Perform matrix inversion for a diagonal matrix.\n", + "4. Uncompute the QST and QCT at the end of the computation." + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "
\n", + "\n", + "
Figure 2. The quantum circuit that for solving the Poisson equation.
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "7", + "metadata": {}, + "source": [ + "Below we define several classical and quantum functions for constructing our quantum linear solver." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import sympy\n", + "\n", + "from classiq import *" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "### Calculating the Diagonal Hamiltonian" + ] + }, + { + "cell_type": "markdown", + "id": "10", + "metadata": {}, + "source": [ + "The eigenvalues of the Poisson equation with Dirichlet boundary conditions in the x direction and Neumann boundary conditions in the y direction are given by:\n", + "$$\n", + "\\lambda_{k,j} \\equiv \\lambda_{x,k} +\\lambda_{y,j}\n", + "$$\n", + "where\n", + "$$\n", + "\\lambda_{x,k} = \\frac{4}{\\Delta x^2} \\sin^2\\left(\\frac{\\pi}{2N_x} (k+1)\\right),\\,\\,\\,\\,\\,\\,\\,\\,\n", + "\\lambda_{y,j} = \\frac{4}{\\Delta y^2} \\sin^2\\left(\\frac{\\pi}{2N_y} j\\right),\n", + "$$\n", + "and $k = 0, 1, \\ldots, N_y-1$ and $j = 0, 1, \\ldots, N_x-1$.\n", + "\n", + "The HHL algorithm requires the application of an Hamiltonian simulation $e^{iH}$. In this notebook this quantum block is implemented using the Suzuki-Trotter built-in funcion. We start with defining a function that gets $N_x$ and $N_y$, and returns the corresponding Hamiltonian. The decomposition of the diagonal matrix to the Pauli basis is done using the Walsh Hadamard transform." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "11", + "metadata": {}, + "outputs": [], + "source": [ + "def get_poisson_dirichletx_neumanny_ham(nx, ny):\n", + " dx2 = ny / nx\n", + " dy2 = nx / ny\n", + " eigenvalues_x = dx2 * 4 * np.sin(np.pi / 2 * np.arange(1, nx + 1) / nx) ** 2\n", + " eigenvalues_y = dy2 * 4 * np.sin(np.pi / 2 * np.arange(0, ny) / ny) ** 2\n", + " num_qubits_x = int(np.log2(nx))\n", + " num_qubits_y = int(np.log2(ny))\n", + "\n", + " # Decompose the eigenvalues of the Poisson equation into Pauli terms and construct the corresponding Hamiltonian\n", + " pauli_coefficients_x = sympy.fwht(eigenvalues_x / nx)\n", + " pauli_coefficients_y = sympy.fwht(eigenvalues_y / ny)\n", + "\n", + " def convert_bitstring_to_pauli_representation(\n", + " bitstring: int, num_qubits: int\n", + " ) -> list[Pauli]:\n", + " return [\n", + " Pauli.Z if (bitstring >> s) & 1 else Pauli.I for s in range(num_qubits)\n", + " ][::-1]\n", + "\n", + " hamiltonian_y = [\n", + " PauliTerm(\n", + " pauli=convert_bitstring_to_pauli_representation(i, num_qubits_y)\n", + " + [Pauli.I] * num_qubits_x,\n", + " coefficient=pauli_coefficients_y[i],\n", + " )\n", + " for i in range(ny)\n", + " ]\n", + " hamiltonian_x = [\n", + " PauliTerm(\n", + " pauli=[Pauli.I] * num_qubits_y\n", + " + convert_bitstring_to_pauli_representation(i, num_qubits_x),\n", + " coefficient=pauli_coefficients_x[i],\n", + " )\n", + " for i in range(nx)\n", + " ]\n", + "\n", + " hamiltonian = hamiltonian_x + hamiltonian_y\n", + " return hamiltonian" + ] + }, + { + "cell_type": "markdown", + "id": "12", + "metadata": {}, + "source": [ + "### Hamiltonian Evolution for QPE\n", + "\n", + "The HHL is based on a [QPE](https://github.com/Classiq/classiq-library/blob/main/tutorials/high_level_modeling_flexible_qpe/high_level_modeling_flexible_qpe.ipynb) applied on $e^{iHt}$. For this, we need to define a function implementing $\\left(e^{iHt}\\right)^p$ for an integer power $p$. Since in our case the Hamiltonian is diagonal, an exact implementation is given by the fisrt order Suzuki-Trotter formula, where in addition $\\left(e^{iHt}\\right)^p = e^{ipHt}$." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "13", + "metadata": {}, + "outputs": [], + "source": [ + "@qfunc\n", + "def powered_hamiltonian_evolution(\n", + " hamiltonian: CArray[PauliTerm], # the hamiltonian H\n", + " scaling: CReal, # the scaling factor t\n", + " p: CInt, # the power\n", + " qba: QArray[QBit],\n", + "):\n", + " suzuki_trotter(\n", + " pauli_operator=hamiltonian,\n", + " evolution_coefficient=p * (-2 * np.pi * scaling),\n", + " order=1,\n", + " repetitions=1,\n", + " qbv=qba,\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "14", + "metadata": {}, + "source": [ + "### Sine and Cosine Transforms\n", + "\n", + "For the system treated in this notebook, diagonalization of the $2D$ Laplacian is given by appying a quantum sine (cosine) transform on the $x$ ($y$) dimension. We define a quantum function for implementing this, using the open library functions:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "15", + "metadata": {}, + "outputs": [], + "source": [ + "@qfunc\n", + "def qsct_2d(xy_variable: QArray[QNum, 2]):\n", + " qst_type2(xy_variable[0])\n", + " qct_type2(xy_variable[1])" + ] + }, + { + "cell_type": "markdown", + "id": "16", + "metadata": {}, + "source": [ + "### Matrix Inversion" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "Finally, we define a matrix inversion using a QPE routine, similar to the example given in the basic [HHL notebook](https://github.com/Classiq/classiq-library/blob/main/algorithms/hhl/hhl/hhl.ipynb)." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "18", + "metadata": {}, + "outputs": [], + "source": [ + "@qfunc\n", + "def inverse_amplitude_load(prefactor: CReal, phase: QNum, ind: QBit):\n", + " ind *= prefactor / phase" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "@qfunc\n", + "def matrix_inversion_HHL(\n", + " prefactor: CReal,\n", + " my_unitary: QCallable[CInt, QArray],\n", + " state: QArray[QBit],\n", + " phase: QNum,\n", + " indicator: Output[QBit],\n", + "):\n", + " allocate(1, indicator)\n", + "\n", + " within_apply(\n", + " compute=lambda: qpe_flexible(\n", + " unitary_with_power=lambda power: my_unitary(power, state),\n", + " phase=phase,\n", + " ),\n", + " action=lambda: inverse_amplitude_load(prefactor, phase, indicator),\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "## Example: non-separable source term" + ] + }, + { + "cell_type": "markdown", + "id": "21", + "metadata": {}, + "source": [ + "We solve an example with a square grid of $N_x,\\, N_y=2^3$. For the source term we take a non-seperable $2^{N_x+N_y}$ vector that represents the function\n", + "$$\n", + "b = F\\left[xy(x-L_x)(y-L_y)\\right].\n", + "$$\n", + "We choose $F(z)=\\tanh(z)^2$, which satisfies the boundary conditions." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "22", + "metadata": {}, + "outputs": [], + "source": [ + "# Set discretization of X and Y axes.\n", + "NUM_QUBITS_X = 3\n", + "NUM_QUBITS_Y = 3\n", + "nx = 2**NUM_QUBITS_X\n", + "ny = 2**NUM_QUBITS_Y\n", + "\n", + "# Set the source term as function of x and y such that b(x,y) = exp(-1/xy(x-Lx)(y-Ly))\n", + "xgrid = (np.arange(nx) + 0.5) / nx\n", + "ygrid = (np.arange(ny) + 0.5) / ny\n", + "zgrid = np.kron(ygrid, xgrid) * np.kron(\n", + " (ygrid - 1), (xgrid - 1)\n", + ") # in classiq the variable order is [y,x]\n", + "b_vector = np.tanh(zgrid) ** 2\n", + "\n", + "# Normalize the source term\n", + "b_vector = b_vector / np.linalg.norm(b_vector)" + ] + }, + { + "cell_type": "markdown", + "id": "23", + "metadata": {}, + "source": [ + "We can plot the source term:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "24", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAG2CAYAAACkgiamAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMKUlEQVR4nO3df1hUdd4//iegDJqAmvHTSdJa0fUHBsKime19kWx5V+zVD9ZacbnLLhVcc/bjKqZgWqKucU+brNxS9MMfN1616u1VhhnJtn6l2GDZtVL64Q8QG8RtYxSTIWa+f3gxOjIDM8PMOe9zzvNxXXMV45mZ9zlz5ryf5/V+z5kAm81mAxEREZEgAuVuABEREdH1GE6IiIhIKAwnREREJBSGEyIiIhIKwwkREREJheGEiIiIhMJwQkREREJhOCEiIiKhMJwQERGRUBhOiIiISCgMJ0RERBpWXFyMuLg4hISEICUlBTU1NS6X3bNnD5KSkjB06FDcdNNNSEhIwPbt2+3/3tnZieXLl2PixIm46aabEBMTg6ysLJw7d86jNjGcEBERadTu3bthMBhQUFCAuro6TJ48Genp6Th//rzT5YcPH45nn30W1dXV+Oc//4ns7GxkZ2fj4MGDAIDLly+jrq4Oq1evRl1dHfbs2YOGhgY8+OCDHrUrgD/8R0REpE0pKSmYOnUqtmzZAgCwWq3Q6/VYvHgxVqxY4dZz3HnnnZg9ezbWrVvn9N//9re/ITk5GWfOnMGtt97q1nMOcK/56mG1WnHu3DmEhoYiICBA7uYQEZHAbDYbLl68iJiYGAQG+m+w4cqVK7BYLP1+HpvN1qNv0+l00Ol0PZa1WCyora1FXl6e/b7AwECkpaWhurrardf68MMP0dDQgI0bN7pcrq2tDQEBARg6dKjb66G5cHLu3Dno9Xq5m0FERArS1NSEkSNH+uW5r1y5gltvvQmtrdZ+P9eQIUNw6dIlh/sKCgqwZs2aHsteuHABXV1diIyMdLg/MjISJ06ccPkabW1tiI2NRUdHB4KCgvCnP/0J9957r9Nlr1y5guXLl2POnDkICwtzez00F05CQ0MBXN3RPNlQD9/xO381iYiIJPLnr170aHmz2Qy9Xm/vO/zBYrGgtdWKqk8iMGSI9xX9S5dsuCflfI/+zVnVpD9CQ0NRX1+PS5cuobKyEgaDAaNHj8Y999zjsFxnZycee+wx2Gw2bN261aPX0Fw46S53hYWFeZbiAoP91SQiIpKIJ8f960kxDWDIkAAMCe3P0NHVyou7/duIESMQFBSElpYWh/tbWloQFRXl8nGBgYG4/fbbAQAJCQk4fvw4CgsLHcJJdzA5c+YMPvzwQ4+3O7+tQ0REpEHBwcFITExEZWWl/T6r1YrKykqkpqa6/TxWqxUdHR32v7uDyVdffYUPPvgAN998s8dt01zlhIiIiK4yGAyYN28ekpKSkJycDKPRiPb2dmRnZwMAsrKyEBsbi8LCQgBAYWEhkpKSMGbMGHR0dODAgQPYvn27fdims7MTjzzyCOrq6vDOO++gq6sLJpMJwNWvIQcHuzcKwXBCRESkUZmZmWhtbUV+fj5MJhMSEhJQUVFhnyTb2Njo8C2l9vZ2LFq0CGfPnsWgQYMQHx+PHTt2IDMzEwDQ3NyM/fv3A7g65HO9w4cP95iX4ormrnNiNpsRHh6OtrY2j8bA7ovO8WOriIhICu99W+zR8t72Gd68xqefR/Zrzsmli1Yk/bTFr22VCuecEBERkVAYToiIiEgoDCdEREQkFIYTIiIiEgrDCREREQlF9nBSXFyMuLg4hISEICUlBTU1NS6X7ezsxNq1azFmzBiEhIRg8uTJqKiokLC1RERE5G+yXudk9+7dMBgMKCkpQUpKCoxGI9LT09HQ0ICIiIgey69atQo7duxAaWkp4uPjcfDgQfzyl7/E0aNHMWXKFBnWgMi3LPH++WExJQg+cVbuJhCRIGS9zklKSgqmTp2KLVu2ALh6CVy9Xo/FixdjxYoVPZaPiYnBs88+i5yca9ccefjhhzFo0CDs2LHDrdfkdU5ITloOH/3F8EK+wOucKINslROLxYLa2lrk5eXZ7wsMDERaWhqqq6udPqajowMhISEO9w0aNAhHjhxx+TodHR0O1/w3m839bDlR7xhA/KOv7crwQqQesoWTCxcuoKury36J3G6RkZE4ceKE08ekp6ejqKgId999N8aMGYPKykrs2bMHXV1dLl+nsLAQzz33nE/bTtrG8CGm3t4XBhciZVHUb+u89NJLmD9/PuLj4xEQEIAxY8YgOzsbZWVlLh+Tl5cHg8Fg/9tsNkOv10vRXFIwBhB1YdWFSFlkCycjRoxAUFAQWlpaHO5vaWlBVFSU08fccsst2LdvH65cuYJ//etfiImJwYoVKzB69GiXr6PT6aDT6XzadlI+hg+6HqsuRGKRLZwEBwcjMTERlZWVyMjIAHB1QmxlZSVyc3N7fWxISAhiY2PR2dmJP//5z3jsscckaDEpCcMH+QqrLkTSk3VYx2AwYN68eUhKSkJycjKMRiPa29uRnZ0NAMjKykJsbCwKCwsBAJ988gmam5uRkJCA5uZmrFmzBlarFb///e/lXA0SBAMJyeHG/Y5hhaj/ZA0nmZmZaG1tRX5+PkwmExISElBRUWGfJNvY2IjAwGtfq7py5QpWrVqFkydPYsiQIbj//vuxfft2DB06VKY1IBEwlJBILPEjGVCI+knW65zIgdc5UQ+GEhIdQ4p4eJ0TZVDUt3WIAIYSUo7ufZUhhcgzsv+2DpEnGExIiSzxI7nvEnmAlRNSBB7YSQ1YSSFyDysnJDSecZIacZ8m6h0rJyQkHrxJ7VhFIXKN4YSEwlBCWsOQQtQTwwkJgaGEtI4hhegazjkh2TGYEF3DeVZErJyQjHgAJnKNlRTSMlZOSHI8MyRyHz8rpEWsnJBkeJAl8g6rKKQ1DCfkdwwlRL7BkEJawXBCfsVg0j9tY3RyN8Gvwr/pkLsJisSQQmrHcEJ+wVDiHbWHkRtdv74MKp5jSCG1Yjghn2Io8YzWwkhvbtwWDCvus8SPZEAhVWE4IZ9gKHEPw4j7GFY8wyoKqQnDCfULQ0nfGEh8g2HFPQwppAYMJ+QVhhLXGEakwfkqvWNIISVjOCGPMZg4YhiRH6sqrjGkkBLxCrHkNl7Z9aq2MTqHG4mH71FP/PySK8XFxYiLi0NISAhSUlJQU1PjctnS0lLMmDEDw4YNw7Bhw5CWltZj+UuXLiE3NxcjR47EoEGDMH78eJSUlHjUJoYT6hMPamBHp3AMK9do/bNMjnbv3g2DwYCCggLU1dVh8uTJSE9Px/nz550uX1VVhTlz5uDw4cOorq6GXq/HrFmz0NzcbF/GYDCgoqICO3bswPHjx/HMM88gNzcX+/fvd7tdATabzdbvtVMQs9mM8PBwtLW1ISwszO3H3Red48dWiUnLBzGtd2Bao9VhIC0O9bz3bbFHy3vbZ3jzGp9+Hokhod7XDC5dtCLppy0etTUlJQVTp07Fli1bAABWqxV6vR6LFy/GihUr+nx8V1cXhg0bhi1btiArKwsAMGHCBGRmZmL16tX25RITE3Hffffh+eefd6tdrJxQD1qslPDMWtu0+t5r8bOuBWaz2eHW0eE8fFssFtTW1iItLc1+X2BgINLS0lBdXe3Wa12+fBmdnZ0YPny4/b5p06Zh//79aG5uhs1mw+HDh/Hll19i1qxZbq8DJ8SSnZYOUlrrhMh9Wpxcy0mzYth3cTJCbAO9fvyVS50A3oder3e4v6CgAGvWrOmx/IULF9DV1YXIyEiH+yMjI3HixAm3XnP58uWIiYlxCDgvv/wynn76aYwcORIDBgxAYGAgSktLcffdd7u9LgwnBED9wYRhhLylpbDCkKIOTU1NDsM6Op1/jn8bNmxAeXk5qqqqEBISYr//5Zdfxscff4z9+/dj1KhR+Oijj5CTk9MjxPSG4UTj1BxKGEjIH7QQVhhSlC0sLMytOScjRoxAUFAQWlpaHO5vaWlBVFRUr4/dvHkzNmzYgA8++ACTJk2y3//DDz9g5cqV2Lt3L2bPng0AmDRpEurr67F582a3wwnnnGiY2oIJ542QHNS8z6ntGEGOgoODkZiYiMrKSvt9VqsVlZWVSE1Ndfm4TZs2Yd26daioqEBSUpLDv3V2dqKzsxOBgY7xIigoCFar1e22sXKiUWo76KixYyDl6d4P1VRN4Y8KqpvBYMC8efOQlJSE5ORkGI1GtLe3Izs7GwCQlZWF2NhYFBYWAgA2btyI/Px87Nq1C3FxcTCZTACAIUOGYMiQIQgLC8PMmTOxbNkyDBo0CKNGjcJf/vIXvPnmmygqKnK7XQwnGqSmYMJQQiJqG6NjQCFFyMzMRGtrK/Lz82EymZCQkICKigr7JNnGxkaHKsjWrVthsVjwyCOPODzP9ZNuy8vLkZeXhyeeeALfffcdRo0ahRdeeAELFixwu10MJxqjlmDCUEKiU1sVhQFFvXJzc5Gbm+v036qqqhz+Pn36dJ/PFxUVhddee61fbeKcEw1hMCGSHvdXIs/JHk48uaY/ABiNRowdOxaDBg2CXq/H0qVLceXKFYlaq1wMJkTyUct+q5bjCIlP1nDi6TX9d+3ahRUrVqCgoADHjx/Hq6++it27d2PlypUSt5zkoJYDPGkT918i98kaToqKijB//nxkZ2fbf7Vw8ODBKCsrc7r80aNHMX36dDz++OOIi4vDrFmzMGfOnD6rLVqnhrMdHthJDdSwH6vheELik21CbPc1/fPy8uz39XVN/2nTpmHHjh2oqalBcnIyTp48iQMHDmDu3LkuX6ejo8PhdwXMZrPvVoIkoYYDujcujgqQuwl+F3pGU787CkAd3+Th5FjyN9nCiTfX9H/88cdx4cIF3HXXXbDZbPjxxx+xYMGCXod1CgsL8dxzz/m07Uqi9LMcrQUTLQSS612/vloKKmoIKET+JPuEWE9UVVVh/fr1+NOf/oS6ujrs2bMH7777LtatW+fyMXl5eWhra7PfmpqaJGwx9YdWgsnFUQH2m5ZpbTsoff9W+okPiU22yok31/RfvXo15s6di6eeegoAMHHiRLS3t+Ppp5/Gs88+2+NyucDVHzzy148eiU7JBw+lH7j7opUO2FtaqaiwgkLknGyVE2+u6X/58mWn1+sHAJtNvQcwrVFrMNFaZcBX1L7dlLy/K/kEiMQm6xViPb2m/wMPPICioiJMmTIFKSkp+Prrr7F69Wo88MAD9pBCVyn1oKHkA7Uzau1Q5aLWigorKESOZA0nnl7Tf9WqVQgICMCqVavQ3NyMW265BQ888ABeeOEFuVaBfEhNwYShxP/UFlSUGlD4zR3yhwCbxsZDzGYzwsPD0dbWhrCwMLcfd190jh9b5VtKrJqoIZgwkIhB6UFFiQFFSeHkvW+LPVre2z7Dm9dY9fEshAwZ6PXzXLnUied/9r5f2yoV/vCfyjCYSIuBRDxKr6gosYLC6gn5GsMJyUqJwYSBRDmUGlSUGFCIfInhREWUVjVRUjBhIFE+pQUVpQUUVk/IlxR1ETZSDyUEE7V/hVXLlPK+KuFzcj2lnSCRuFg5UQklHRREPuAqocMi31FCNUVpFRQiX2A4IUmJGEwYSAgQO6goKaBweId8geFEBZRSNREpmDCQUG9EDCpKCihE/cVwQpIQIZgwkJA3RAoqSgkorJ5QfzGcKJwSqiZyBhMGEvIlEYKKUgIKUX8wnCgYg4lzDCQkBTmDihICCqsn1B/8KjH5jdTBRClfDyX1kWPfE2GolMhfWDlRKNGrJnIEEzXquNUidxP8RtcYLHcTfO7iqABJqyiiV1BYPSFvsXJCPsdg0n8dt1pUHUwA9a4jKyiORD+RIjExnCiQyB92BpP+UWuH3ZvudVbTejOgEPUPwwn5DIOJ99TWOXtLTduBAeUakU+oSEwMJwoj6oecwcRzaqwY+IpatgsDCpF3GE6o3xhMPKOWjlcKathWDChXiXpiRWJiOFEQET/cDCbuU0NHKxelbzsGFCLPMJyQ1xhM3KP0jlUkSh4KY0AR8wSLxMRwohCifagZTPqm1E5UKZS4fRlQiNzDcKIADCbKCSZKPrNXKqVtb60HFNGOZyQmhhPyCIOJc0rrINVISe+B1gMKUV8YTgQn0lkGg0lPSuoQtUIp74mWA4pIxzUSE8MJuYXBxJFSOkAtU8J7xIBCIiguLkZcXBxCQkKQkpKCmpoal8uWlpZixowZGDZsGIYNG4a0tLRel1+wYAECAgJgNBo9ahPDicBE+fAymFyjhA6PHIk+D0jLAYXkt3v3bhgMBhQUFKCurg6TJ09Geno6zp8/73T5qqoqzJkzB4cPH0Z1dTX0ej1mzZqF5ubmHsvu3bsXH3/8MWJiYjxuF8MJ9YrBRPzOjdwn6vuo1YAiygmYlhUVFWH+/PnIzs7G+PHjUVJSgsGDB6OsrMzp8jt37sSiRYuQkJCA+Ph4vPLKK7BaraisrHRYrrm5GYsXL8bOnTsxcOBAj9vFcCIoET60Wg8monZk1H8ivrdaDSgkH4vFgtraWqSlpdnvCwwMRFpaGqqrq916jsuXL6OzsxPDhw+332e1WjF37lwsW7YMP/3pT71qG8MJOaXlYCJix0X+Idp7rcWAIsKJmNqYzWaHW0dHh9PlLly4gK6uLkRGRjrcHxkZCZPJ5NZrLV++HDExMQ4BZ+PGjRgwYAB++9vfer0OA7x+JPmN3B9WrQYTkTopklb3e69rDJa5JVc/D6FnbJK9XtsYHcK/cd55kbQOt/wEAy55f/z9sb0DwPvQ6/UO9xcUFGDNmjX9a5wTGzZsQHl5OaqqqhASEgIAqK2txUsvvYS6ujoEBHh/bGc4IQdaDCYMJdTt+n1BzqCitYBiiR+J4BNnZXt9tWlqakJYWJj9b53O+XF9xIgRCAoKQktLi8P9LS0tiIqK6vU1Nm/ejA0bNuCDDz7ApEmT7Pf/9a9/xfnz53Hrrbfa7+vq6sLvfvc7GI1GnD592q11YDgRjJxVE60FE6WEkriRrXI3wedOn71F7ib0Se5qitYCCvlOWFiYQzhxJTg4GImJiaisrERGRgYA2Ce35ubmunzcpk2b8MILL+DgwYNISkpy+Le5c+c6DPEAQHp6OubOnYvs7Gy310GIOSeefMf6nnvuQUBAQI/b7NmzJWyxfzCY+J+SvnkTN7JVlcEEuLZuSlg/OfcXLc1BkXs4W6sMBgNKS0vxxhtv4Pjx41i4cCHa29vtQSIrKwt5eXn25Tdu3IjVq1ejrKwMcXFxMJlMMJlMuHTpEgDg5ptvxoQJExxuAwcORFRUFMaOHet2u2SvnHR/x7qkpAQpKSkwGo1IT09HQ0MDIiIieiy/Z88eWCzXDhT/+te/MHnyZDz66KNSNltVtBBMlBBGuimhw/al7vUVvZoiVyWFFRTyp8zMTLS2tiI/Px8mkwkJCQmoqKiwT5JtbGxEYOC1OsbWrVthsVjwyCOPODyPr+e1BNhsNun2eidSUlIwdepUbNmyBcDVkpJer8fixYuxYsWKPh9vNBqRn5+Pb7/9FjfddFOfy5vNZoSHh6Otrc2tsle3+6Jz3F7WG3KdNag9mDCUKI/oIaWb1CFFyoACQLaA4u+5J+99W+zR8t72Gd68xvT/y8WAm/o3Ifb/e2iLX9sqFVmHdXzxHetXX30Vv/rVr1wGk46Ojh5fq6Kr1BxMlDJ0A6h7+MYbStkeUg8RamWIh8M7BMgcTvr7Heuamhp89tlneOqpp1wuU1hYiPDwcPvtxq9YiUALH0apDqxKCSVKmnMhFyVtH6n2O60EFCIhJsR669VXX8XEiRORnJzscpm8vDy0tbXZb01NTRK2UFxSHnSkDCaiU1KHKwolbTM1BhQ5aOGEjXon64TY/nzHur29HeXl5Vi7dm2vy+l0Opff8RaB2j+EDCZXKaVzFZmSJs76ez6KlJNkOUGW5CBr5eT671h36/6OdWpqaq+Pfeutt9DR0YFf//rX/m6m6khVNWEwUdZZv1IoYZuygtJ/aj9xo97JPqzj6Xesu7366qvIyMjAzTffLHWTfUaODx+DiTSU0IEqnejbWE0BhXNPSGqyX+fE0+9YA0BDQwOOHDmC999/X44mkyBEDCYid5ZqpZThHn+R+jooUuJl7bVL9nACALm5uS4vlVtVVdXjvrFjx0Lmy7P0m5pLllKczTGY0I3iRrYKF1CkmH8iFc49ISnJPqxD0mFp1n8YTMTA94FIHRhOiPqJHaJYRHs/1DT3hEgqDCfkU1ob0hGtI6Sr+L74B6uvJBWGEyIvsQMUm0jvj0iBmkgJGE40Qi1nPKIc5EXq+Mg1Lb1Pah3aUfOXB8g1hhPyGbUeHG+kpQ5PDUR5v0QJ1kRKwHAiA54JeEeEg7soHR15hu+b76ilCktiYzjRAB5MfIMdnLKJ8P75O2BrpXpJ6sdwQj7h74Oi3FUTETo26j++j0TKwHBC1Ad2aOrC95NIfAwnJDw5qybsyNRJzvdVDUM7HComf2M4UTkpDiJqHedmMFE3vr/KwS8RaA/DCZET7Li0Qa73WQ3VEyJ/YjghockxpMNgoi18v4nEM0DuBmiN2sqTajtDE62jujfqhNxN8JtDpni5m2AXN7IVp8/eIulrdtxqga4xWNLX9KW2MTqEf9MhdzNIpVg5UTGlT1qTumrCYCIt0dZPtPe/v9R24kDawnBCXlPTwU+0jkm0jttfRFtPqfcDua/fQyQqhhMSkpQHbQYTeYm2vqLtD0RaxHBCmiZaRyRaRy0V0dZbyv3Cn0Hc39VNpQ8dk7gYTlTK3wcNfx70pKqaMJiIRbT1F23/0Dq1fZmAesdwQpokWscjWscsF9G2g2j7CZFWMJyQUKSomojW4YjWIctNtO0hxf7CoR0iRwwnElJLWVLJ39JhMFEG0baLaPsNkdoxnKgQz2ScE62DEa0DFo1o28ff+w+/Vkx0DcMJCcOfB2cGE2USbTuJth+5S8nVTtImhhPyiBIPcqJ1KKJ1uKITbXv5c39SavWE1VryNYYTlVHqQcJfB2UGE3UQbbuJtl9phVrm7VHfGE5ItUTrQETrYJVGtO3nr/3LX0FdiVVP0i6GE3Kbvw5u/jgYM5iok2jbUbT9jMgbxcXFiIuLQ0hICFJSUlBTU+Ny2dLSUsyYMQPDhg3DsGHDkJaW1mN5m82G/Px8REdHY9CgQUhLS8NXX33lUZsYTiTCcqR0ROswROtQlU607emP/U2Jc0+UOqSsdbt374bBYEBBQQHq6uowefJkpKen4/z5806Xr6qqwpw5c3D48GFUV1dDr9dj1qxZaG5uti+zadMm/PGPf0RJSQk++eQT3HTTTUhPT8eVK1fcbhfDiYoo8eDg64Mwg4k2iLZdRdvvXOHQDt2oqKgI8+fPR3Z2NsaPH4+SkhIMHjwYZWVlTpffuXMnFi1ahISEBMTHx+OVV16B1WpFZWUlgKtVE6PRiFWrVuGhhx7CpEmT8Oabb+LcuXPYt2+f2+1iOCG38KDmOdE6ULURbfsqJaCQ+pnNZodbR0eH0+UsFgtqa2uRlpZmvy8wMBBpaWmorq5267UuX76Mzs5ODB8+HABw6tQpmEwmh+cMDw9HSkqK288JCBBOPBnrAoDvv/8eOTk5iI6Ohk6nw09+8hMcOHBAotaSL6m5aiJax6lWat7OShzaof5pPDcCp8/e4vWt8dwIAIBer0d4eLj9VlhY6PT1Lly4gK6uLkRGRjrcHxkZCZPJ5Fably9fjpiYGHsY6X5cf54TAAa4vaQfdI91lZSUICUlBUajEenp6WhoaEBERESP5S0WC+69915ERETg7bffRmxsLM6cOYOhQ4dK33giF9TcYYro3qgTOGSKl7sZAK4G5NNnb5G7Gb26OCoAoWdsfnnutjE6hH/j/CydpNPU1ISwsDD73zqdf4b8N2zYgPLyclRVVSEkJMSnzy1r5cTTsa6ysjJ899132LdvH6ZPn464uDjMnDkTkydPlrjl4vHnfBMlDOmIUjVhMJGHWrc7qyc98csFfQsLC3O4uQonI0aMQFBQEFpaWhzub2lpQVRUVK+vsXnzZmzYsAHvv/8+Jk2aZL+/+3HePOf1ZAsn3ox17d+/H6mpqcjJyUFkZCQmTJiA9evXo6ury+XrdHR09Bh/I/mp8aCr1g5SKUTZ/qIEZaK+BAcHIzEx0T6ZFYB9cmtqaqrLx23atAnr1q1DRUUFkpKSHP7ttttuQ1RUlMNzms1mfPLJJ70+541kCyfejHWdPHkSb7/9Nrq6unDgwAGsXr0aL774Ip5//nmXr1NYWOgw9qbX6326HiQ/EToDUTpGrVPj++CPIO/PaqgSvzWoZQaDAaWlpXjjjTdw/PhxLFy4EO3t7cjOzgYAZGVlIS8vz778xo0bsXr1apSVlSEuLg4mkwkmkwmXLl0CAAQEBOCZZ57B888/j/379+PYsWPIyspCTEwMMjIy3G6XrHNOPGW1WhEREYFt27YhKCgIiYmJaG5uxh/+8AcUFBQ4fUxeXh4MBoP9b7PZLHlAUXIZ0h8HMbVVTdTYISqZCHNQlDD3hAgAMjMz0draivz8fJhMJiQkJKCiosJeOGhsbERg4LU6xtatW2GxWPDII484PE9BQQHWrFkDAPj973+P9vZ2PP300/j+++9x1113oaKiwqN5KbKFE2/GuqKjozFw4EAEBQXZ7xs3bhxMJhMsFguCg4N7PEan0/ltMpAotHymInfVhMFETCIEFF/quNUCXWPP4xuRL+Tm5iI3N9fpv1VVVTn8ffr06T6fLyAgAGvXrsXatWu9bpNswzrejHVNnz4dX3/9NaxWq/2+L7/8EtHR0U6DCYlHTVUTBhOxyf3+yB2c+6KEie6kXbJ+W8fTsa6FCxfiu+++w5IlS/Dll1/i3Xffxfr165GTkyPXKqia6AcvOQ/+cnd85B41vU9KCvZaruaSb8g658TTsS69Xo+DBw9i6dKlmDRpEmJjY7FkyRIsX75crlUgDyjp4NobNXV4WqC2IR4iLZB9QqwnY10AkJqaio8//tjPrVIOrZ6hyFU1YTBRJrkCiugTY/15QTZ/ssSPRPCJs3I3g/xI9svXk5h8PaSjhqoJg4myqeH9U8PniMgdsldO1E7JXyMWlRxVEzk7tkfC6mR7bX9623yn5K8pRwVF9OqJv/BS9tQfrJwomFaHdKTGYOIfcq2b0isovq6eiD7xnbSJ4YR6EHlIR+qqCYOJf2kloIj+tWIi0TCcELnAYCINrQQUX+LcE1I7hhPyK6UeRBlMpKWFgCJy9cRfQzsceiZvMZwolL8+9CKPP0t1cGcwkYcWAoovKTX4+wq/bKBuDCfkN0o8eDKYyEvtAUWL1RMibzCckCJIcVBnMNE2pVZQiNSI4YSIiIiEwnBCREREQmE4ISIiIqEwnBAREZFQGE6IiIhIKAwnREREJBSGEyIiIhIKwwkREREJheGEiIiIhMJwQkREREJhOCEiIiKhMJwQERGRUBhOiIiISCgMJ0RERCQUhhMiIiISCsMJERERCYXhhIiIiITCcEJERERCYTghIiIioTCcEBERkVAYToiIiEgoDCdEREQkFCHCSXFxMeLi4hASEoKUlBTU1NS4XPb1119HQECAwy0kJETC1hIREamHJ33w559/jocffhhxcXEICAiA0Wh0ulxzczN+/etf4+abb8agQYMwceJEfPrpp263yeNwMm/ePHz00UeePsyl3bt3w2AwoKCgAHV1dZg8eTLS09Nx/vx5l48JCwvDt99+a7+dOXPGZ+0hIiLSCk/74MuXL2P06NHYsGEDoqKinC7z73//G9OnT8fAgQPx3nvv4YsvvsCLL76IYcOGud0uj8NJW1sb0tLScMcdd2D9+vVobm729CkcFBUVYf78+cjOzsb48eNRUlKCwYMHo6yszOVjAgICEBUVZb9FRkb2qw1ERERa5GkfPHXqVPzhD3/Ar371K+h0OqfLbNy4EXq9Hq+99hqSk5Nx2223YdasWRgzZozb7fI4nOzbtw/Nzc1YuHAhdu/ejbi4ONx33314++230dnZ6dFzWSwW1NbWIi0t7VqDAgORlpaG6upql4+7dOkSRo0aBb1ej4ceegiff/65y2U7OjpgNpsdbkRERGp1Y5/X0dHhdDlv++C+7N+/H0lJSXj00UcRERGBKVOmoLS01KPn8GrOyS233AKDwYB//OMf+OSTT3D77bdj7ty5iImJwdKlS/HVV1+59TwXLlxAV1dXj8pHZGQkTCaT08eMHTsWZWVl+L//+z/s2LEDVqsV06ZNw9mzZ50uX1hYiPDwcPtNr9d7trJEREQSCG4Khq7R+1twUzAAQK/XO/R7hYWFTl/Pmz7YHSdPnsTWrVtxxx134ODBg1i4cCF++9vf4o033nD7OQZ4/eoAvv32Wxw6dAiHDh1CUFAQ7r//fhw7dgzjx4/Hpk2bsHTp0v48vVOpqalITU21/z1t2jSMGzcO//M//4N169b1WD4vLw8Gg8H+t9lsZkAhIiLVampqQlhYmP1vV8Mv/mK1WpGUlIT169cDAKZMmYLPPvsMJSUlmDdvnlvP4XE46ezsxP79+/Haa6/h/fffx6RJk/DMM8/g8ccft2+MvXv34r/+67/6DCcjRoxAUFAQWlpaHO5vaWlxOdHmRgMHDsSUKVPw9ddfO/13nU4n+RtDREQkl7CwMIdw4oov+mBnoqOjMX78eIf7xo0bhz//+c9uP4fHwzrR0dGYP38+Ro0ahZqaGnz66adYsGCBw4b4+c9/jqFDh/b5XMHBwUhMTERlZaX9PqvVisrKSofqSG+6urpw7NgxREdHe7oqREREmuWLPtiZ6dOno6GhweG+L7/8EqNGjXL7OTyunPz3f/83Hn300V6vLTJ06FCcOnXKreczGAyYN28ekpKSkJycDKPRiPb2dmRnZwMAsrKyEBsbax8zW7t2LX72s5/h9ttvx/fff48//OEPOHPmDJ566ilPV4WIiEjTPO2DLRYLvvjiC/v/Nzc3o76+HkOGDMHtt98OAFi6dCmmTZuG9evX47HHHkNNTQ22bduGbdu2ud0uj8PJ3LlzPX1IrzIzM9Ha2or8/HyYTCYkJCSgoqLCPkGnsbERgYHXCjz//ve/MX/+fJhMJgwbNgyJiYk4evRojxISERER9c7TPvjcuXOYMmWK/e/Nmzdj8+bNmDlzJqqqqgBc/brx3r17kZeXh7Vr1+K2226D0WjEE0884Xa7Amw2m803q6gMZrMZ4eHhaGtrc2tMrtt90TlevZ4lfqRXj+tL2xj/zKO5OCrAZ8/VcavFZ88VN7LVZ8/lyr1RJ/z+Gq48ElYn22uL5m3znbK99iFTvN9f4/TZW3z2XLrGYJ89FwCEnvF9dxD+jfOvsfpC8Ann39LszXvfFnu0vLd9hjevMWblegT144rnXVeu4Jv1K/3aVqkIcfl6IiIiom4MJ0RERCQUhhMiIiISCsMJERERCYXhhIiIiITCcEJERERCYTghIiIioTCcEBERkVAYToiIiEgoDCdEREQkFIYTIiIiEgrDCREREQmF4YSIiIiEwnBCREREQmE4ISIiIqEwnBAREZFQGE6IiIhIKAwnREREJBSGEyIiIhIKwwkREREJheGEFOH02Vv8/hqHTPF+fw1X3jbfKdtri0TO7SDn+y+C0DM2uZtAZMdwQn6jawyWuwkeY0CRjxaCiRQhm0gNGE4UKvybDr88r8hnT1Id2BlQpKeFYOJrSgz/vhR84qzcTSA/Yjghv1LqAZQBRTpaCSYiV038dVLir5MoUj+GE1IUKQ/wDCj+p5Vg4mtKDf1E7mI4oR58fRal5AMpA4r/aCmYiFw1IRIRw4mCabVkKvWBngHF97QUTHzN12Ff5HlmpF0D5G6A2gWfOAtL/Ei5m0H9dMgUj3ujTsjy2moNKHKQI5hotWqi1ZMn8g1WTsgpDu30pPQzbq1Tw/unhs8RkTsYThROq2cncp2NqqGD0yK53jfRqyZKHdLh14jVj+GEJKOWsz4GFGVRy/ulls8PkTuECCfFxcWIi4tDSEgIUlJSUFNT49bjysvLERAQgIyMDP82UKNEP6uS86xULR2e2sn5PoleNfEnrVZ0yXdkDye7d++GwWBAQUEB6urqMHnyZKSnp+P8+fO9Pu706dP4f//v/2HGjBkStZR8QU1nfwwoYuP70zvRTz5I22QPJ0VFRZg/fz6ys7Mxfvx4lJSUYPDgwSgrK3P5mK6uLjzxxBN47rnnMHr0aAlbKyYtn6XIfXbKDlBMantf1BTqidwhazixWCyora1FWlqa/b7AwECkpaWhurra5ePWrl2LiIgIPPnkk32+RkdHB8xms8NNakqevOWPsyu1HWjV1hEqnQjvh9yhWU5aPlki35E1nFy4cAFdXV2IjIx0uD8yMhImk8npY44cOYJXX30VpaWlbr1GYWEhwsPD7Te9Xt/vdpNYROgIROgQSZ3vgz/CPId06HqezPv8/PPP8fDDDyMuLg4BAQEwGo09liksLMTUqVMRGhqKiIgIZGRkoKGhwaM2yT6s44mLFy9i7ty5KC0txYgRI9x6TF5eHtra2uy3pqYmP7eS3KG26gmgzo5RSUTZ/iKEZSJ3eTrv8/Llyxg9ejQ2bNiAqKgop8v85S9/QU5ODj7++GMcOnQInZ2dmDVrFtrb291ul6xXiB0xYgSCgoLQ0tLicH9LS4vTlf7mm29w+vRpPPDAA/b7rFYrAGDAgAFoaGjAmDFjHB6j0+mg0+n80HqxhH/TgbYx/lnP0DM2XBwV4Jfn9pXTZ29B3MhWuZsh65VktUyUYOJragzx/aXkYXIRXT/vEwBKSkrw7rvvoqysDCtWrOix/NSpUzF16lQAcPrvAFBRUeHw9+uvv46IiAjU1tbi7rvvdqtdslZOgoODkZiYiMrKSvt9VqsVlZWVSE1N7bF8fHw8jh07hvr6evvtwQcfxM9//nPU19dzyEZh1HrgVWtHKSqRtrcSqib+HNLhfBMx3DjPsqPD+fvi7bxPT7W1tQEAhg8f7vZjZP9tHYPBgHnz5iEpKQnJyckwGo1ob2+3p7isrCzExsaisLAQISEhmDBhgsPjhw4dCgA97iftEaV6ArCCIhWRgomvqTW8k2uhTTYEBXsfHrssVx9744l6QUEB1qxZ02P53uZ9njjhm+OX1WrFM888g+nTp3vUT8seTjIzM9Ha2or8/HyYTCYkJCSgoqLCvrEaGxsRGKioqTGqpIShHdEwoPiXaMFECVUT0oampiaEhYXZ/5ZzakNOTg4+++wzHDlyxKPHyR5OACA3Nxe5ublO/62qqqrXx77++uu+b5BC+XPeib/oGoPRcavFZ88nUvUEYEDxF7UHE39VTfgtHW0ICwtzCCeueDrv01O5ubl455138NFHH2HkyJEePZYlCYlwEpd0RDuDFa0jVTrRtqdo+5tcON9EeTyd9+kum82G3Nxc7N27Fx9++CFuu+02j5+D4YTc5q+zLn+cJYrWYYjWoSqVaNtRtP2MyFMGgwGlpaV44403cPz4cSxcuLDHvM+8vDz78haLxf6FFIvFgubmZtTX1+Prr7+2L5OTk4MdO3Zg165dCA0Nhclkgslkwg8//OB2u4QY1iHfUeLQjr9wiEddtBJMOKTjGivQvufpvM9z585hypQp9r83b96MzZs3Y+bMmfZpGFu3bgUA3HPPPQ6v9dprr+E3v/mNW+1iOCEh+HruSTcGFHXQSjAhkoMn8z7j4uJgs/UedPv6d3dwWIc8osSzL9E6EtE6WtGJtr38uT8p9evDnG9CvsZwQsLw54GZAUWZRNtOou1H7lLiSQVpG8OJCvEsxjnROhbROl7RiLZ9/L3/KLVqQuQPDCcSUstkLn+ehfn7AM2AogyibRfR9hsitWM4Ic0RraMRrSOWm2jbQ4r9xZ+h3N9DOqzUkj8wnJBXlD6GzYAiJtG2g2j7CZFWMJyolJLPZqQaexet4xGtY5aaaOsv1f7BuSbuUcuwOLmH4YQ0jQFFDKKtt2j7hbc4pENKxXBCXlPyxNjridYRidZR+5to6yva/kCkRQwnRBCvQxKtw/YX0dZT6v1AyRNhifyJ4UTFlF5ylXosngFFWqKtn2jvP5GW8bd1JBZ84iws8SPlbobPhJ6x4eKoALmb4TMi/hYP+Z8cwUTpE2GVfvJDYmPlhIQmxwGcZ9Daosb3m0M6pHQMJ0ROqLHDop7kep+VXjWRGr9GrD0MJyonRenV32dpch3IGVDUje8vkbgYToh6wQ5MneR8X/0dtqUY0uF8E/I3hhOiPjCgqAvfTyLxMZxoAId2+o8dmjrI/T7KvR8TKQXDiQw4uUuZ5O7YqH+08P7xWzqkFgwnpBginHVqoYNTI75vvsP5JiQFhhPyGa2ctbGjUxZR3i8RwjWRUjCcaIRaznZEOcCL0uFR77T0Pqn15IDD4NrEcELkJS11fEok0vsjSqjuL7Wc5JD4GE7Ip6Q4exPpQC9SB0jX8H0hUjaGE6J+YkcoFtHeDynCtFqHdEi7GE40RE0lWZGqJ4B4HaJW8X0gUgeGE5moeZKXVs/i2DHKS8TtL1qI7g81ndyQ+IQIJ8XFxYiLi0NISAhSUlJQU1Pjctk9e/YgKSkJQ4cOxU033YSEhARs375dwtaSO7Q296Tb6bO3CNlJqpmo21yq/VPNJwNqPomj3skeTnbv3g2DwYCCggLU1dVh8uTJSE9Px/nz550uP3z4cDz77LOorq7GP//5T2RnZyM7OxsHDx6UuOX9J8cHT8qzH60GFEDcDlNNRN7GagsmrJqQ1GQPJ0VFRZg/fz6ys7Mxfvx4lJSUYPDgwSgrK3O6/D333INf/vKXGDduHMaMGYMlS5Zg0qRJOHLkiMQtVy4GFOl0d6CidqJKJPr2VFswkQurJtomazixWCyora1FWlqa/b7AwECkpaWhurq6z8fbbDZUVlaioaEBd999t9NlOjo6YDabHW4i0cIHUOsBpZvonarolLD91BhMWDUhOcgaTi5cuICuri5ERkY63B8ZGQmTyeTycW1tbRgyZAiCg4Mxe/ZsvPzyy7j33nudLltYWIjw8HD7Ta/X+3QdlErqAw4DyjVK6GRFopTtpcZgIhctnLRR72Qf1vFGaGgo6uvr8be//Q0vvPACDAYDqqqqnC6bl5eHtrY2+62pqUnaxrpBrg+iWgMKQ4ryKWk4TMp9TupgwqoJyWWAnC8+YsQIBAUFoaWlxeH+lpYWREVFuXxcYGAgbr/9dgBAQkICjh8/jsLCQtxzzz09ltXpdNDpdD5tt5qEf9OBtjHSbZ/QMzZcHBXg99fp7iw6brX4/bX6q7sDjhvZKnNL5KeEMNJN6hCslWDCqgkBMldOgoODkZiYiMrKSvt9VqsVlZWVSE1Ndft5rFYrOjqUnfDl/ECqsYLSjZUUZVDSusuxTzGYkNbIPqxjMBhQWlqKN954A8ePH8fChQvR3t6O7OxsAEBWVhby8vLsyxcWFuLQoUM4efIkjh8/jhdffBHbt2/Hr3/9a7lWQRXUHFCAax2KEoKKkjrq/lLSusq1/2glmJB8PLnWGAC89dZbiI+PR0hICCZOnIgDBw44/PulS5eQm5uLkSNHYtCgQfZv4npC1mEdAMjMzERrayvy8/NhMpmQkJCAiooK+yTZxsZGBAZey1Dt7e1YtGgRzp49i0GDBiE+Ph47duxAZmamXKvgM8EnzsISP1K211frEM+NlDLko5ROW83kDrNaCiasmsij+1pjJSUlSElJgdFoRHp6OhoaGhAREdFj+aNHj2LOnDkoLCzEf/7nf2LXrl3IyMhAXV0dJkyYAOBq0eHDDz/Ejh07EBcXh/fffx+LFi1CTEwMHnzwQbfaFWCz2dQ/9fs6ZrMZ4eHhaGtrQ1hYmNuPuy86x4+tciRnQAEgaUABIEtAuZ7oIYWkJ3coAbQVTADpwsl73xZ7tLy3fYY3r5Ew9wUEBYd4/Txdliuo3/6sR21NSUnB1KlTsWXLFgBXp0no9XosXrwYK1as6LF8ZmYm2tvb8c4779jv+9nPfoaEhAR7dWTChAnIzMzE6tWr7cskJibivvvuw/PPP+9Wu2Qf1iHxqH2I50ZKGe4h/xNlX2Awof648dperuZkenOtserqaoflASA9Pd1h+WnTpmH//v1obm6GzWbD4cOH8eWXX2LWrFlur4PswzrUk9zDO4B2hniup5ThHvI9EQJJN60FE7om7FQHBgzw/jj4449X38sbr+dVUFCANWvW9Fi+t2uNnThxwulrmEymPq9N9vLLL+Ppp5/GyJEjMWDAAAQGBqK0tNTlxVKdYTghl7QYUADHjopBRd1ECiWANoMJqya+19TU5DCsI/XlNF5++WV8/PHH2L9/P0aNGoWPPvoIOTk5iImJ6VF1cYXhRFAiVE8A7QaUbqymqJNooQTQZjAh/wgLC3Nrzok31xqLiorqdfkffvgBK1euxN69ezF79mwAwKRJk1BfX4/Nmze7HU4454T6pLU5KM6IMheBvCfy18m1GkxYNZGXN9caS01NdVgeAA4dOmRfvrOzE52dnQ7fsgWAoKAgWK1Wt9vGyonARKmeAKygdGMlRXlEDCPXYzAhORkMBsybNw9JSUlITk6G0Wjsca2x2NhYFBYWAgCWLFmCmTNn4sUXX8Ts2bNRXl6OTz/9FNu2bQNwtWozc+ZMLFu2DIMGDcKoUaPwl7/8BW+++SaKiorcbhfDCbmNAeUahhTxiR5KAO0GExKHp9camzZtGnbt2oVVq1Zh5cqVuOOOO7Bv3z77NU4AoLy8HHl5eXjiiSfw3XffYdSoUXjhhRewYMECt9vF65y4ScrrnNxIlOpJN61dB8VdDCpiUEIoAbQdTOSsmoh8nZO778rHgAHeX+fkxx+v4KMja/3aVqlwzokCiFb+5BwU50Sdz6AVStr+Wg4mRO5gOCGvMKC4pqROUulEnuTqitaDiWgnWyQmhhOFEPEDzYDSO6V1mkqi1G2r9WBC5C6GE+oXBpS+KbUjFZGStyWDiZgnWSQmhhMFEfWDzYDiHiUOQYhC6duNwYTIMwwn5BMMKJ5RemcrFTVsJwaTq0Q9uSIxMZwojMgfcAYUz6mh8/U1NVWYGEyIvMNwQj7FgOIdtXTG/aG2bcBgco3IJ1UkJoYTBRL9g86A4j21ddDuUOM6M5hcI/rxisTEcEJ+wYDSP2oa2nBFrevHYELUf/xtHYUS6UcBXeFv8fiGGjtwtWIwccSqCXmLlRMFU8IHX44KitqqKCQ+OfY70YMJUX+wckJ+J3UFBXA8g1VjNYXkJ2cIVkIwUcLJE4mL4UThlDC8A8gTULoxqJCviFCVU0IwIeovhhOSjJwBpRuDCnlKhEDSTSnBhFUT6i+GExVQSvUEECOgdGNQIVdECiTdlBJMiHyB4YQkJ1JA6dbdGTGkaJeIgaSbkoIJqybkCwwnKqGk6gkgZkABWE3RGpEDSTcGE9IihhOSjagBpRuDijopIZB0U1IwIfIlhhMVUVr1BBA/oHRjUFE2JQWSbkoLJqyakC8xnKgMA4r/MagogxIDSTelBRMiX2M4ISEoLaB0Y1ARi5IDSTclBhNWTcjXGE5USInVE0C5AaUbg4o81BBIuikxmBD5A8MJCUXpAaUbg4p/qSmQdFNqMGHVhPxBiB/+Ky4uRlxcHEJCQpCSkoKamhqXy5aWlmLGjBkYNmwYhg0bhrS0tF6X1yolHzCUepB2pftH4dTYoUpNrdtRbfs8UX/JHk52794Ng8GAgoIC1NXVYfLkyUhPT8f58+edLl9VVYU5c+bg8OHDqK6uhl6vx6xZs9Dc3Cxxy8mf1HqwZlDxnNq3mZL3dSWfBJHYZA8nRUVFmD9/PrKzszF+/HiUlJRg8ODBKCsrc7r8zp07sWjRIiQkJCA+Ph6vvPIKrFYrKisrJW65+JR+4FDyQdsdau90+0Mr20bt+ziRt2Sdc2KxWFBbW4u8vDz7fYGBgUhLS0N1dbVbz3H58mV0dnZi+PDhTv+9o6MDHR3XDgBms7l/jSZJqWUOSl84R0Wd80h6o/RgovSTHxKbrJWTCxcuoKurC5GRkQ73R0ZGwmQyufUcy5cvR0xMDNLS0pz+e2FhIcLDw+03vV7f73YriRoOIEo/iHtKK1UDQFvrej2l79NqOK6Q2BT9bZ0NGzagvLwcVVVVCAkJcbpMXl4eDAaD/W+z2ay5gKIGWqmg3EhrnbYWKD2YEElB1nAyYsQIBAUFoaWlxeH+lpYWREVF9frYzZs3Y8OGDfjggw8wadIkl8vpdDrodNrr1K6n1Oue3EirAYXUQw3BhFUTkoKswzrBwcFITEx0mMzaPbk1NTXV5eM2bdqEdevWoaKiAklJSVI0VfHUckBRw8GdtIn7LpH7ZB/WMRgMmDdvHpKSkpCcnAyj0Yj29nZkZ2cDALKyshAbG4vCwkIAwMaNG5Gfn49du3YhLi7OPjdlyJAhGDJkiGzrQdJhBYWURi3BRC0nOSQ+2cNJZmYmWltbkZ+fD5PJhISEBFRUVNgnyTY2NiIw8FqBZ+vWrbBYLHjkkUccnqegoABr1qyRsumKo5bhHYABhZRDLcGESEqyhxMAyM3NRW5urtN/q6qqcvj79OnT/m+QijGgEElDbaGEVROSkuwXYSPpqekgE/5Nh+o6AVI2Ne6TajpmUE+e/IQMALz11luIj49HSEgIJk6ciAMHDrhcdsGCBQgICIDRaPSoTQwnGqW2g40aOwRSHjXug2o7VpAjT39C5ujRo5gzZw6efPJJ/P3vf0dGRgYyMjLw2Wef9Vh27969+PjjjxETE+NxuxhONEyNB53ukKLGToLEc/3+prZ9LvjEWVUeI8iRpz8h89JLL+EXv/gFli1bhnHjxmHdunW48847sWXLFoflmpubsXjxYuzcuRMDBw70uF0MJxqn5gOQmjsOkocW9ik1HxPIUfdPyFx/hfW+fkKmurq6xxXZ09PTHZa3Wq2YO3culi1bhp/+9KdetU2ICbEkv+6DkVomyzpzY2fCybTkDrWGEGcYStThxt+Qc3Ux0t5+QubEiRNOn9tkMvX5kzMbN27EgAED8Nvf/tbbVWA4IUdaCCndGFbIGS2FkW4MJWII/vIcBgQGe/34QKsFAHr8RIuUl9qora3FSy+9hLq6OgQEeP8jpgwn5JSavnLsrus7JQYV7dBiGOnGUKJOTU1NCAsLs//t6idcvPkJmaioqF6X/+tf/4rz58/j1ltvtf97V1cXfve738FoNLp9ORDOOSGXtDz2rIW5BVrF91bbn20tCAsLc7i5Cife/IRMamqqw/IAcOjQIfvyc+fOxT//+U/U19fbbzExMVi2bBkOHjzo9jqwckJ90tJQjyscAlIurQYQZxhI6Eae/oTMkiVLMHPmTLz44ouYPXs2ysvL8emnn2Lbtm0AgJtvvhk333yzw2sMHDgQUVFRGDt2rNvtYjghtzGkXMOwIjYGEkcMJeSKpz8hM23aNOzatQurVq3CypUrcccdd2Dfvn2YMGGCT9sVYLPZbD59RsGZzWaEh4ejra3NYUyuL/dF5/ixVcrEkOIcg4r0GEZcYzBx9N63xR4t722f4c1rpEU81a8JsT9aLfjg/Ct+batUWDkhr7GS4hyrKv7HMNI3hhJSMoYT6jctfrPHEwwr/ccw4j6GElIDhhPyCVZR3Mew4h4GEs8wlJCaMJyQTzGkeI5h5SqGEe8wlJAaMZyQXzCkeE8rF4NjGOk/BhNSK4YT8iuGlP5hB07OMJSQ2vEKsSQJHkyJ+o9XdiWtYOWEJMMqCpF3GEhIaxhOSHIMKUTuYSghrWI4IdkwpBA5x1BCWsc5JyQ7jqMTXcPPAhErJyQQVlJIyxhKiK5h5YSEw4M0aQkrh0Q9sXJCQmIVhdSOgYTINYYTEhpDCqkNQwlR3zisQ4rA0jepAfdhIvewckKKwkoKKRFDCZFnWDkhReLBnpSAFT8i77ByQorFKgqJioGEqH8YTkjxGFJIFAwlRL4h+7BOcXEx4uLiEBISgpSUFNTU1Lhc9vPPP8fDDz+MuLg4BAQEwGg0StdQEl53Cf36G5E/cX8j8g9ZKye7d++GwWBASUkJUlJSYDQakZ6ejoaGBkRERPRY/vLlyxg9ejQeffRRLF26VIYWk9L01mGw0kJ9YeAgkoes4aSoqAjz589HdnY2AKCkpATvvvsuysrKsGLFih7LT506FVOnTgUAp/9O5Im+Oh6GF21gACESj2zhxGKxoLa2Fnl5efb7AgMDkZaWhurqap+9TkdHBzo6Oux/m81mnz03qRurLurA8EGkPLKFkwsXLqCrqwuRkZEO90dGRuLEiRM+e53CwkI899xzPns+IoBVF5EwfBCpj+q/rZOXlweDwWD/22w2Q6/Xy9gi0gJWXXyLAYRIW2QLJyNGjEBQUBBaWloc7m9paUFUVJTPXken00Gn0/ns+Yj6i1WXnhg+iOh6soWT4OBgJCYmorKyEhkZGQAAq9WKyspK5ObmytUsItmxoyYirZN1WMdgMGDevHlISkpCcnIyjEYj2tvb7d/eycrKQmxsLAoLCwFcnUT7xRdf2P+/ubkZ9fX1GDJkCG6//XbZ1oOIiIh8R9ZwkpmZidbWVuTn58NkMiEhIQEVFRX2SbKNjY0IDLx2nbhz585hypQp9r83b96MzZs3Y+bMmaiqqpK6+UREROQHsk+Izc3NdTmMc2PgiIuLg81mk6BVREREJBfZL19PREREdD2GEyIiIhIKwwkREREJheGEiIiIhMJwQkREREJhOCEiItKw4uJixMXFISQkBCkpKaipqel1+bfeegvx8fEICQnBxIkTceDAAYd/t9lsyM/PR3R0NAYNGoS0tDR89dVXHrWJ4YSIiEijdu/eDYPBgIKCAtTV1WHy5MlIT0/H+fPnnS5/9OhRzJkzB08++ST+/ve/IyMjAxkZGfjss8/sy2zatAl//OMfUVJSgk8++QQ33XQT0tPTceXKFbfbxXBCRESkUUVFRZg/fz6ys7Mxfvx4lJSUYPDgwSgrK3O6/EsvvYRf/OIXWLZsGcaNG4d169bhzjvvxJYtWwBcrZoYjUasWrUKDz30ECZNmoQ333wT586dw759+9xul+wXYZNa90XczGazR4/70WrxR3OIiEhCnh77u5eX4gKgP9osgLWfj0fPdXT1A7gWiwW1tbXIy8uz3xcYGIi0tDRUV1c7fY3q6moYDAaH+9LT0+3B49SpUzCZTEhLS7P/e3h4OFJSUlBdXY1f/epXbq2L5sLJxYsXAQB6vV7mlhARkdTCw1/x6nEXL15EeHi4j1tzVXBwMKKiolBlerPfzzVkyJAe/VtBQQHWrFnTY9kLFy6gq6vL/pMx3SIjI3HixAmnz28ymZwubzKZ7P/efZ+rZdyhuXASExODpqYmhIaGIiAgQO7mCMdsNkOv16OpqQlhYWFyN0fxuD19j9vUt7g9e2ez2XDx4kXExMT47TVCQkJw6tQpWCz9r9DbbLYefZuzqonoNBdOAgMDMXLkSLmbIbywsDAeqHyI29P3uE19i9vTNX9VTK4XEhKCkJAQv7/O9UaMGIGgoCC0tLQ43N/S0oKoqCinj4mKiup1+e7/trS0IDo62mGZhIQEt9vGCbFEREQaFBwcjMTERFRWVtrvs1qtqKysRGpqqtPHpKamOiwPAIcOHbIvf9tttyEqKsphGbPZjE8++cTlczqjucoJERERXWUwGDBv3jwkJSUhOTkZRqMR7e3tyM7OBgBkZWUhNjYWhYWFAIAlS5Zg5syZePHFFzF79myUl5fj008/xbZt2wAAAQEBeOaZZ/D888/jjjvuwG233YbVq1cjJiYGGRkZbreL4YQc6HQ6FBQUKHKMUkTcnr7Hbepb3J7alpmZidbWVuTn58NkMiEhIQEVFRX2Ca2NjY0IDLw2yDJt2jTs2rULq1atwsqVK3HHHXdg3759mDBhgn2Z3//+92hvb8fTTz+N77//HnfddRcqKio8GrYKsEnx/SgiIiIiN3HOCREREQmF4YSIiIiEwnBCREREQmE4ISIiIqEwnGiQJz+PXVpaihkzZmDYsGEYNmwY0tLS+vw5ba3x9OfGu5WXlyMgIMCjr9dphafb9Pvvv0dOTg6io6Oh0+nwk5/8pMfPuGuZp9vTaDRi7NixGDRoEPR6PZYuXerRL8oS9ZuNNKW8vNwWHBxsKysrs33++ee2+fPn24YOHWpraWlxuvzjjz9uKy4utv3973+3HT9+3Pab3/zGFh4ebjt79qzELReTp9uz26lTp2yxsbG2GTNm2B566CFpGqsQnm7Tjo4OW1JSku3++++3HTlyxHbq1ClbVVWVrb6+XuKWi8nT7blz506bTqez7dy503bq1CnbwYMHbdHR0balS5dK3HLSMoYTjUlOTrbl5OTY/+7q6rLFxMTYCgsL3Xr8jz/+aAsNDbW98cYb/mqionizPX/88UfbtGnTbK+88opt3rx5DCc38HSbbt261TZ69GibxWKRqomK4un2zMnJsf3Hf/yHw30Gg8E2ffp0v7aT6Hoc1tGQ7p/Hvv6nrPv6eewbXb58GZ2dnRg+fLi/mqkY3m7PtWvXIiIiAk8++aQUzVQUb7bp/v37kZqaipycHERGRmLChAlYv349urq6pGq2sLzZntOmTUNtba196OfkyZM4cOAA7r//fknaTATwCrGa4s3PY99o+fLliImJcTjYaZU32/PIkSN49dVXUV9fL0ELlcebbXry5El8+OGHeOKJJ3DgwAF8/fXXWLRoETo7O1FQUCBFs4XlzfZ8/PHHceHCBdx1112w2Wz48ccfsWDBAqxcuVKKJhMB4IRY8sCGDRtQXl6OvXv3Sv7rmWpw8eJFzJ07F6WlpRgxYoTczVENq9WKiIgIbNu2DYmJicjMzMSzzz6LkpISuZumSFVVVVi/fj3+9Kc/oa6uDnv27MG7776LdevWyd000hBWTjTEm5/H7rZ582Zs2LABH3zwASZNmuTPZiqGp9vzm2++wenTp/HAAw/Y77NarQCAAQMGoKGhAWPGjPFvowXnzT4aHR2NgQMHIigoyH7fuHHjYDKZYLFYEBwc7Nc2i8yb7bl69WrMnTsXTz31FABg4sSJ9t9JefbZZx1+Z4XIX7iXaYg3P48NAJs2bcK6detQUVGBpKQkKZqqCJ5uz/j4eBw7dgz19fX224MPPoif//znqK+vh16vl7L5QvJmH50+fTq+/vpre9ADgC+//BLR0dGaDiaAd9vz8uXLPQJId/Cz8afYSCpyz8glaZWXl9t0Op3t9ddft33xxRe2p59+2jZ06FCbyWSy2Ww229y5c20rVqywL79hwwZbcHCw7e2337Z9++239tvFixflWgWheLo9b8Rv6/Tk6TZtbGy0hYaG2nJzc20NDQ22d955xxYREWF7/vnn5VoFoXi6PQsKCmyhoaG2//3f/7WdPHnS9v7779vGjBlje+yxx+RaBdIgDutojKc/j71161ZYLBY88sgjDs9TUFCANWvWSNl0IXm6Palvnm5TvV6PgwcPYunSpZg0aRJiY2OxZMkSLF++XK5VEIqn23PVqlUICAjAqlWr0NzcjFtuuQUPPPAAXnjhBblWgTQowGZjnY6IiIjEwVM6IiIiEgrDCREREQmF4YSIiIiEwnBCREREQmE4ISIiIqEwnBAREZFQGE6IiIhIKAwnREREJBSGEyIiIhIKwwkREREJheGESMNaW1sRFRWF9evX2+87evQogoODHX7JlohISvxtHSKNO3DgADIyMnD06FGMHTsWCQkJeOihh1BUVCR304hIoxhOiAg5OTn44IMPkJSUhGPHjuFvf/sbdDqd3M0iIo1iOCEi/PDDD5gwYQKamppQW1uLiRMnyt0kItIwzjkhInzzzTc4d+4crFYrTp8+LXdziEjjWDkh0jiLxYLk5GQkJCRg7NixMBqNOHbsGCIiIuRuGhFpFMMJkcYtW7YMb7/9Nv7xj39gyJAhmDlzJsLDw/HOO+/I3TQi0igO6xBpWFVVFYxGI7Zv346wsDAEBgZi+/bt+Otf/4qtW7fK3Twi0ihWToiIiEgorJwQERGRUBhOiIiISCgMJ0RERCQUhhMiIiISCsMJERERCYXhhIiIiITCcEJERERCYTghIiIioTCcEBERkVAYToiIiEgoDCdEREQkFIYTIiIiEsr/D3b1Yuvw2TMPAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "b_matrix = b_vector.reshape((nx, ny))\n", + "xmesh, ymesh = np.meshgrid(xgrid, ygrid)\n", + "plt.contourf(xmesh, ymesh, b_matrix.transpose())\n", + "plt.xlabel(\"x\")\n", + "plt.ylabel(\"y\")\n", + "ax = plt.gca()\n", + "ax.axis(\"equal\")\n", + "plt.colorbar()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "25", + "metadata": {}, + "source": [ + "Next, we calculate the Hamiltonian for the specific discretization: " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "26", + "metadata": {}, + "outputs": [], + "source": [ + "hamiltonian = get_poisson_dirichletx_neumanny_ham(nx, ny)" + ] + }, + { + "cell_type": "markdown", + "id": "27", + "metadata": {}, + "source": [ + "Now, we need to define the resolution (QPE size) for the quantum solver. The required number of QPE phase qubits depends on the condition number $\\kappa$ of the inverted matrix. In the case of the Laplacian matrix, this parameter is of the order of the matrix dimension $O(2^{N_x+N_y})$, which eliminates the exponential advantage. However, in our example we have a smooth source term, and high modes are expected to have minor effects on the solution. The highest mode of $\\mathcal{L}$ is $\\lambda_{\\max}=8$, whereas the smallest one is $\\lambda_{x,0}\\sim 4\\pi^2/(2N_x)^2$. Let us assume that the highest mode participating in the solution is $\\sim 2^6 \\lambda_{x,0}$, and take a QPE size of 6 qubits." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "28", + "metadata": {}, + "outputs": [], + "source": [ + "# number of qubits for the QPE\n", + "QPE_SIZE = 6\n", + "MAX_EIG = 4 * (np.pi / (2 * nx)) ** 2 * 2**6\n", + "# parameters for the amplitude preparation\n", + "PREFACTOR = 2**-QPE_SIZE" + ] + }, + { + "cell_type": "markdown", + "id": "29", + "metadata": {}, + "source": [ + "We build the model and synthesize it:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "30", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Opening: https://platform.classiq.io/circuit/506deb07-15d3-455f-883f-aa52f269fcda?version=0.43.1\n" + ] + } + ], + "source": [ + "@qfunc\n", + "def main(\n", + " x_variable: Output[QNum[NUM_QUBITS_X, False, 0]],\n", + " y_variable: Output[QNum[NUM_QUBITS_Y, False, 0]],\n", + " phase: Output[QNum],\n", + " indicator: Output[QBit],\n", + "):\n", + "\n", + " xy_variable = QArray(\"xy_variable\", QNum[NUM_QUBITS_X, False, 0], 2)\n", + " prepare_amplitudes(b_vector.tolist(), 0.0, xy_variable)\n", + "\n", + " allocate(QPE_SIZE, phase)\n", + "\n", + " within_apply(\n", + " compute=lambda: qsct_2d(xy_variable),\n", + " action=lambda: matrix_inversion_HHL(\n", + " prefactor=PREFACTOR,\n", + " my_unitary=lambda p, target: powered_hamiltonian_evolution(\n", + " hamiltonian=hamiltonian,\n", + " scaling=1 / MAX_EIG,\n", + " p=p,\n", + " qba=target,\n", + " ),\n", + " state=xy_variable,\n", + " phase=phase,\n", + " indicator=indicator,\n", + " ),\n", + " )\n", + "\n", + " bind(xy_variable, [x_variable, y_variable])\n", + "\n", + "\n", + "qmod = create_model(main, constraints=Constraints(max_width=18))\n", + "write_qmod(qmod, \"discrete_poisson_solver\", decimal_precision=12)\n", + "qprog = synthesize(qmod)\n", + "show(qprog)" + ] + }, + { + "cell_type": "markdown", + "id": "31", + "metadata": {}, + "source": [ + "### Executing and Plotting the Result" + ] + }, + { + "cell_type": "markdown", + "id": "32", + "metadata": {}, + "source": [ + "We run the quantum program on a statevector simulator to retrieve the full solution." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "33", + "metadata": {}, + "outputs": [], + "source": [ + "from classiq.execution import (\n", + " ClassiqBackendPreferences,\n", + " ClassiqSimulatorBackendNames,\n", + " ExecutionDetails,\n", + " ExecutionPreferences,\n", + " set_quantum_program_execution_preferences,\n", + ")\n", + "\n", + "execution_preferences = ExecutionPreferences(\n", + " backend_preferences=ClassiqBackendPreferences(backend_name=\"simulator_statevector\")\n", + ")\n", + "qprog = set_quantum_program_execution_preferences(qprog, execution_preferences)\n", + "result = execute(qprog).result()[0].value" + ] + }, + { + "cell_type": "markdown", + "id": "34", + "metadata": {}, + "source": [ + "We define a postprocess function that gets the quantum solution out of the execution and returns solution" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "35", + "metadata": {}, + "outputs": [], + "source": [ + "import itertools\n", + "\n", + "\n", + "# Retrieve the solution of the Poisson equation from the results of the quantum program.\n", + "def extract_result_from_statevector_simulation(\n", + " result: ExecutionDetails,\n", + " normalization: float,\n", + " indicator_name: str,\n", + " phase_name: str,\n", + " x_name: str,\n", + " y_name: str,\n", + ") -> np.ndarray:\n", + " state_matrix = np.zeros((nx, ny), dtype=complex)\n", + "\n", + " # Filter only the successful states.\n", + " for state in result.parsed_state_vector:\n", + " if state.state[indicator_name] == 1.0 and state.state[phase_name] == 0.0:\n", + " state_matrix[\n", + " int(state.state[x_name]), int(state.state[y_name])\n", + " ] += state.amplitude\n", + "\n", + " # Normalize the solution of the Poisson equation.\n", + " global_phase = np.angle(state_matrix[0, 0])\n", + " result_matrix = np.real(state_matrix / np.exp(1j * global_phase)) / normalization\n", + "\n", + " return result_matrix" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "36", + "metadata": {}, + "outputs": [], + "source": [ + "result_matrix = extract_result_from_statevector_simulation(\n", + " result, PREFACTOR / MAX_EIG, \"indicator\", \"phase\", \"x_variable\", \"y_variable\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "37", + "metadata": {}, + "source": [ + "The resulting statevector is given up to a sign. The sign in the center of the solution is expected to be positive. We can correct accordingly:" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "38", + "metadata": {}, + "outputs": [], + "source": [ + "if result_matrix[nx // 2, ny // 2] < 0:\n", + " result_matrix = -result_matrix" + ] + }, + { + "cell_type": "markdown", + "id": "39", + "metadata": {}, + "source": [ + "Finally, we print the result and compare it to the classical solution" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "40", + "metadata": {}, + "outputs": [], + "source": [ + "## getting the classical solution\n", + "Hx = (\n", + " 2 * np.diag(np.ones(nx))\n", + " - np.diag(np.ones(nx - 1), 1)\n", + " - np.diag(np.ones(nx - 1), -1)\n", + ")\n", + "Hx[0, 0] = Hx[nx - 1, nx - 1] = 3\n", + "\n", + "Hy = (\n", + " 2 * np.diag(np.ones(ny))\n", + " - np.diag(np.ones(ny - 1), 1)\n", + " - np.diag(np.ones(ny - 1), -1)\n", + ")\n", + "Hy[0, 0] = Hy[ny - 1, ny - 1] = 1\n", + "\n", + "# we need to flip the order of x and y\n", + "b_classical = (b_matrix.T).reshape(nx * ny)\n", + "\n", + "H = np.kron(Hx, np.identity(ny)) + np.kron(np.identity(nx), Hy)\n", + "classical_result_vector = np.linalg.solve(H, b_classical)\n", + "classical_result_matrix = classical_result_vector.reshape((nx, ny))" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "41", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9UAAAHwCAYAAACotEyzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAAByIElEQVR4nO3de3gU5f3//1cOZgNCAjSQk4FwEBCRxILkGxGVNhgPxdpPbRGrwQhUkCiaWhWFBDwQFaVYRakURG0tVGqp/cIPimlTpaZSofl4BOUYQBJASxJiSSSb3x98WVlyYGeyszuz+3xcVy7NZnb23tllXvO+5557Ipqbm5sFAAAAAAAMiwx2AwAAAAAAcCqKagAAAAAATKKoBgAAAADAJIpqAAAAAABMoqgGAAAAAMAkimoAAAAAAEyiqAYAAAAAwCSKagAAAAAATKKoBgAAAADAJIpqAAAAAABMoqgGADjeW2+9pXHjxiklJUURERFavXr1GZ/z29/+VhkZGercubOSk5N166236osvvrC+sQAAwBSjef/6669r7Nix6tmzp+Li4pSdna3169e3WG7RokVKT09XbGyssrKytGnTJkPtoqgGADhefX29MjIytGjRIp+W/8c//qG8vDxNmjRJH330kV577TVt2rRJU6ZMsbilAADALKN5/9Zbb2ns2LFau3atNm/erDFjxmjcuHH697//7Vlm5cqVKiwsVHFxsbZs2aKMjAzl5ubq4MGDPrcrorm5udnwuwEAwKYiIiL0xz/+Udddd12byzz55JN6/vnntWPHDs9jzzzzjB5//HHt27cvAK0EAAAd4Uvet+b888/X+PHjVVRUJEnKysrSRRddpGeffVaS5Ha7lZaWpjvuuEP333+/T+vkTDUAIOxkZ2dr7969Wrt2rZqbm1VdXa1Vq1bp6quvDnbTAACARdxut+rq6tSjRw9JUmNjozZv3qycnBzPMpGRkcrJyVF5ebnP6432e0sBAGHj2LFjamxstGTdzc3NioiI8HrM5XLJ5XJ1eN2jRo3Sb3/7W40fP17Hjh3T8ePHNW7cOJ+HkwEAEC6cmvWtefLJJ3X06FH9+Mc/liQdPnxYTU1NSkxM9FouMTFRW7du9Xm9FNUAAFOOHTum3r3P1qFDbkvW36VLFx09etTrseLiYs2ZM6fD6/744481Y8YMFRUVKTc3VwcOHNDPf/5zTZ06VUuXLu3w+gEACAVOzvrTvfrqq5o7d67+9Kc/qVevXn5dN0U1AMCUxsZGHTrkVtm7vdSlS8SZn2DA0aPNujzroPbu3au4uDjP4/7quS4pKdGoUaP085//XJI0bNgwnX322Ro9erQeeeQRJScn++V1AABwMidn/alWrFihyZMn67XXXvMa6p2QkKCoqChVV1d7LV9dXa2kpCSf109RDQDokC5dItSlq7+n6DjRIx4XF+cVtP7y1VdfKTraOwKjoqIknRiKBgAAvuHErD/pd7/7nW699VatWLFC11xzjdffYmJiNHz4cJWWlnomPHO73SotLVVBQYHPr0FRDQBwvKNHj2r79u2e33ft2qWKigr16NFDvXv31syZM7V//369/PLLkqRx48ZpypQpev755z3Dv++66y6NHDlSKSkpwXobAACgHUbz/tVXX9XEiRP19NNPKysrS1VVVZKkTp06KT4+XpJUWFioiRMnasSIERo5cqQWLlyo+vp65efn+9wuimoAgOO99957GjNmjOf3wsJCSdLEiRO1fPlyHThwQJWVlZ6/33LLLaqrq9Ozzz6rn/3sZ+rWrZu+853v6PHHHw942wEAgG+M5v0LL7yg48ePa/r06Zo+fbrn8ZPLS9L48eN16NAhFRUVqaqqSpmZmVq3bl2Lycvaw32qAQCm1NbWKj4+Xu99lOj3IWFH69wacX61ampqLB0SBgAA2kbW+4b7VAMAAAAAYBJFNQAAAAAAJlFUAwAAAABgEkU1AAAAAAAmUVQDAAAAAGASRTUAAAAAACZRVAMAAAAAYBJFNQAAAAAAJlFUAwAAAABgEkU1AAAAAAAmUVQDAAAAAGASRTUAAAAAACZRVAMAAAAAYBJFNQAAAAAAJlFUAwAAAABgEkU1AAAAAAAmUVQDAAAAAGASRTUAAAAAACZRVAMAAAAAYBJFNQAAAAAAJlFUAwAAAABgEkU1AAAAAAAmUVQDAAAAAGASRTUAAAAAACZRVAMAAAAAYBJFNQAAAAAAJlFUAwAAAABgEkU1AAAAAAAmUVQDAAAAAGASRTUAAAAAACZRVAMAAAAAYBJFNQAAAAAAJlFUAwAAAABgEkU1AAAAQkZ6erpuueWWoL3+LbfcovT0dMvWf/nll+vyyy+3bP0AjKOohqN99NFHuummm5SamiqXy6WUlBTddNNN+vjjj4PdNI+PP/5Yc+bM0e7du4PdFAAAHG3Hjh267bbb1K9fP8XGxiouLk6jRo3S008/rf/+97/Bbh6AMBUd7AYAZr3++uuaMGGCevTooUmTJqlv377avXu3li5dqlWrVmnlypX6/ve/H+xm6uOPP9bcuXN1+eWXW9pzDQBAKFuzZo1+9KMfyeVyKS8vT0OHDlVjY6M2btyon//85/roo4/0wgsvBLuZWrJkidxud7CbASCAKKrhSDt27NDNN9+sfv366a233lLPnj09f5sxY4ZGjx6tm266Se+//7769u0bxJYCAICO2rVrl2644Qb16dNHf/3rX5WcnOz52/Tp07V9+3atWbMmiC38xllnnRXsJgAIMIZ/w5Hmz5+vr776Si+88IJXQS1JCQkJ+tWvfqWjR49q/vz5ktq+vmnOnDmKiIjweuzFF1/Ud77zHfXq1Usul0tDhgzR888/3+K56enp+t73vqeNGzdq5MiRio2NVb9+/fTyyy97llm+fLl+9KMfSZLGjBmjiIgIRUREqKysTJIUERGhOXPmtLruU68HW758uSIiIrRx40bdeeed6tmzp7p166bbbrtNjY2NOnLkiPLy8tS9e3d1795d9957r5qbm33ZlAAA2N4TTzyho0ePaunSpV4F9UkDBgzQjBkzWn3ul19+qXvuuUcXXHCBunTpori4OF111VX63//93xbLPvPMMzr//PPVuXNnde/eXSNGjNCrr77q+XtdXZ3uuusupaeny+VyqVevXho7dqy2bNniWaa1Yw63262nn35aF1xwgWJjY9WzZ09deeWVeu+99zzL+Hr8AcB+OFMNR/rzn/+s9PR0jR49utW/X3rppUpPT9ef//xnPffcc4bW/fzzz+v888/Xtddeq+joaP35z3/W7bffLrfbrenTp3stu337dl1//fWaNGmSJk6cqGXLlumWW27R8OHDdf755+vSSy/VnXfeqV/+8pd64IEHdN5550mS579G3XHHHUpKStLcuXP1z3/+Uy+88IK6deumd955R71799a8efO0du1azZ8/X0OHDlVeXp6p1wEAwE7+/Oc/q1+/frr44osNP3fnzp1avXq1fvSjH6lv376qrq7Wr371K1122WX6+OOPlZKSIunEsO0777xT119/vWbMmKFjx47p/fff17vvvqsbb7xRkjR16lStWrVKBQUFGjJkiL744gtt3LhRn3zyib797W+32YZJkyZp+fLluuqqqzR58mQdP35cb7/9tv75z39qxIgRkowdfwCwF4pqOE5NTY0+//zzM14vPWzYML3xxhuqq6sztP6///3v6tSpk+f3goICXXnllVqwYEGLUNu2bZveeustT3H/4x//WGlpaXrxxRf15JNPql+/fho9erR++ctfauzYsR2erTMxMVFr165VRESEbr/9dm3fvl3z58/Xbbfd5unN/ulPf6r09HQtW7aMohoA4Hi1tbXav3+/6XlSLrjgAn366aeKjPxmgObNN9+swYMHa+nSpZo9e7akE9dsn3/++XrttdfaXNeaNWs0ZcoUPfXUU57H7r333nZf/29/+5uWL1+uO++8U08//bTn8Z/97Gdeo8qMHH8AsBeGf8NxThbJXbt2bXe5k383WlSfGmg1NTU6fPiwLrvsMu3cuVM1NTVeyw4ZMsTrbHnPnj01aNAg7dy509Br+mrSpElew9WzsrLU3NysSZMmeR6LiorSiBEjLGsDAACBVFtbK+nMud8Wl8vlKaibmpr0xRdfqEuXLho0aJDXsO1u3bpp3759+te//tXmurp166Z3331Xn3/+uc+v/4c//EEREREqLi5u8bdTM93I8QcAe6GohuP4WizX1dUpIiJCCQkJhtb/j3/8Qzk5OTr77LPVrVs39ezZUw888IAktQi13r17t3h+9+7d9Z///MfQa/rq9NeLj4+XJKWlpbV43Ko2AAAQSHFxcZKMd5Kf5Ha79Ytf/ELnnnuuXC6XEhIS1LNnT73//vteuX7fffepS5cuGjlypM4991xNnz5d//jHP7zW9cQTT+jDDz9UWlqaRo4cqTlz5pyxE3vHjh1KSUlRjx492l3OyPEHAHuhqIbjxMfHKyUlRe+//367y73//vs655xzFBMT02IyspOampq8ft+xY4e++93v6vDhw1qwYIHWrFmjDRs26O6775akFrfIiIqKanW9HZ0k7PR2nen1WnucicoAAKEgLi5OKSkp+vDDD009f968eSosLNSll16q3/zmN1q/fr02bNig888/3yvXzzvvPG3btk0rVqzQJZdcoj/84Q+65JJLvM4w//jHP9bOnTv1zDPPKCUlRfPnz9f555+v/+//+/869B6NHn8AsBeKajjSuHHjtGvXLm3cuLHVv7/99tvavXu3Z+bt7t2768iRIy2W27Nnj9fvf/7zn9XQ0KA33nhDt912m66++mrl5OR4Dckyqq2Cvq12NTY26sCBA6ZfDwCAUPO9731PO3bsUHl5ueHnrlq1SmPGjNHSpUt1ww036IorrlBOTk6rxwVnn322xo8frxdffFGVlZW65ppr9Oijj+rYsWOeZZKTk3X77bdr9erV2rVrl771rW/p0UcfbfP1+/fvr88//1xffvllm8tYcfwBIHAoquFI99xzjzp37qzbbrtNX3zxhdffvvzyS02dOlVxcXEqKCiQdCLQampqvM5uHzhwQH/84x+9nnvyjO+pZ3lramr04osvmm7r2WefLUmthnf//v311ltveT32wgsvtHmmGgCAcHTvvffq7LPP1uTJk1VdXd3i7zt27PCaBOxUUVFRLUZvvfbaa9q/f7/XY6cfT8TExGjIkCFqbm7W119/raamphbDsHv16qWUlBQ1NDS02fYf/vCHam5u1ty5c1v87WS7rDj+ABA4zP4NRxowYIBefvllTZgwQRdccIEmTZqkvn37avfu3Vq6dKn+85//aMWKFerbt68k6YYbbtB9992nH/zgB7rzzjv11Vdf6fnnn9fAgQO9Jim54oorFBMTo3Hjxum2227T0aNHtWTJEvXq1cv02ePMzExFRUXp8ccfV01NjVwul+c+lJMnT9bUqVP1wx/+UGPHjtX//u//av369YavAwcAIJT1799fr776qsaPH6/zzjtPeXl5Gjp0qBobG/XOO+/otdde0y233NLqc7/3ve/poYceUn5+vi6++GJ98MEH+u1vf6t+/fp5LXfFFVcoKSlJo0aNUmJioj755BM9++yzuuaaa9S1a1cdOXJE55xzjq6//nplZGSoS5cuevPNN/Wvf/3Lazbw040ZM0Y333yzfvnLX+qzzz7TlVdeKbfbrbfffltjxoxRQUGBJccfAAKHohqO9cMf/lBbtmxRSUmJfv3rX+vgwYNyu92KjY3V5s2bNWTIEM+y3/rWt/THP/5RhYWFuvfee9W3b1+VlJTos88+8yqqBw0apFWrVmnWrFm65557lJSUpGnTpqlnz5669dZbTbUzKSlJixcvVklJiSZNmqSmpib97W9/U69evTRlyhTt2rVLS5cu1bp16zR69Ght2LBB3/3udzu8fQAACCXXXnut3n//fc2fP19/+tOf9Pzzz8vlcmnYsGF66qmnNGXKlFaf98ADD6i+vl6vvvqqVq5cqW9/+9tas2aN7r//fq/lbrvtNv32t7/VggULdPToUZ1zzjm68847NWvWLElS586ddfvtt+svf/mLXn/9dbndbg0YMEDPPfecpk2b1m7bX3zxRQ0bNkxLly7Vz3/+c8XHx2vEiBGe+25bcfwBIHAimpnNCCHk5Zdf1i233KKbbrpJL7/8crCbA4S02tpaxcfH672PEtWlq3+vJjpa59aI86tVU1Pjmfm3PW+99Zbmz5+vzZs3ey7tuO6669p9TkNDgx566CH95je/UVVVlZKTk1VUVMQBLAAA/4+dsl4yl/dlZWUqLCzURx99pLS0NM2aNctrZMucOXNaXJ4xaNAgbd261ef3wplqhJS8vDwdOHBA999/v8455xzNmzcv2E0CEAD19fXKyMjQrbfeqv/5n//x6Tk//vGPVV1draVLl2rAgAE6cOAAM+wCAGBjRvN+165duuaaazR16lT99re/VWlpqSZPnqzk5GTl5uZ6ljv//PP15ptven6PjjZWJlNUI+Tcd999uu+++4LdDAABdNVVV+mqq67yefl169bp73//u3bu3Om5d2x6erpFrQMAAP5gNO8XL16svn37euY9OO+887Rx40b94he/8Cqqo6OjlZSUZLpdzP4NALCt2tpar5/2Ztg14o033tCIESP0xBNPKDU1VQMHDtQ999yj//73v35ZPwAA8I1VWS9J5eXlysnJ8XosNze3xe35PvvsM6WkpKhfv376yU9+osrKSkOvw5lqAECHrK7LUGzzWX5d57GjX0v6i9LS0rweLy4u1pw5czq8/p07d2rjxo2KjY3VH//4Rx0+fFi33367vvjiC25hAwDAaZyY9ZJUVVWlxMREr8cSExNVW1ur//73v+rUqZOysrK0fPlyDRo0SAcOHNDcuXM1evRoffjhh+ratatPr0NRDQCwrb1793pNXuJyufyyXrfbrYiICP32t79VfHy8JGnBggW6/vrr9dxzz6lTp05+eR0AANA+q7LeV6cOJx82bJiysrLUp08f/f73v9ekSZN8WgdFNQDAtuLi4nyeEdSI5ORkpaamegpq6cR1Vs3Nzdq3b5/OPfdcv78mAABoyaqsl07c2ra6utrrserqasXFxbXZgd6tWzcNHDhQ27dv9/l1HFFUu91uff755+ratasiIiKC3RwAcJzm5mbV1dUpJSVFkZFMpzFq1Ci99tprOnr0qLp06SJJ+vTTTxUZGalzzjknyK0LT2Q9AHQMWd9Sdna21q5d6/XYhg0blJ2d3eZzjh49qh07dujmm2/2+XUcUVR//vnnLcbaAwCM27t3b0gWjUePHvXqUd61a5cqKirUo0cP9e7dWzNnztT+/fs996+/8cYb9fDDDys/P19z587V4cOH9fOf/1y33norQ7+DhKwHAP8I1ayXjOf91KlT9eyzz+ree+/Vrbfeqr/+9a/6/e9/rzVr1njWcc8992jcuHHq06ePPv/8cxUXFysqKkoTJkzwuV2OKKpPXiB+eUKeoiNjOry+xnOTO7yOttT2i7Vs3aerOyc4PfmNvRuD8rpwppjKjv+bNaPrvuaAvVbczmOWrTvmswN+Wc9xd6PKDr/s84QbTvPee+9pzJgxnt8LCwslSRMnTtTy5ct14MABr5k8u3Tpog0bNuiOO+7QiBEj9K1vfUs//vGP9cgjjwS87Tjh5Hfz4v9zn6KjrbuerrZv4K7Vq0sL7hn3xjRn53XvlMPBboIjVH6eEOwmdEjM3uAcJ5yq697AHDPE7fLfrNatOX68Qe/88/GQzXrJeN737dtXa9as0d13362nn35a55xzjn7961973U5r3759mjBhgr744gv17NlTl1xyif75z3+qZ8+ePrfLEUX1yWFg0ZExfimq3dHWFb5RMYErqqNiAx/WDX0aFanAvUc43/FBkmtP4AMzKiZwRXW0hXtSf+zzThWqw2ovv/xyNTe3/ZkvX768xWODBw/Whg0bLGwVjPBkfbRL0ZbmdGCK6ro+EYoKyCu1LbKTc4d/pp9zSFJgJytyqn7n1mn3Pt8P/u3m+EDJFaQO+JMCdcwQHR2YDA7VrJfM5f3ll1+uf//7320+Z8WKFR1ul3P3tjZUMyCAZ6l7B6egBgAA9tfAqLKwcqITAkCwUFQDsFwwOmSC0fEEAOg4CkQATkNRDZ9wlhodxXcIAHAmFNTmOXnbBXtkRV0fOuLRMRTVANBBgbz0A4D9BfsAPdgFCoLHyYV1OKjpzzwBoYqi2oECPayVM4wAAMBqFIQAnMoRs3/7U+Og1GA3AUEQKkHt5Nk9pRMdNIGcCbyud4S6VgZuFnAA9sZZIvsKlZy2g/RzDjn+eAFwGs5UI+QR1ACAcOHEod/ktP85cZs68bsLnGSqqF60aJHS09MVGxurrKwsbdq0qc1lv/76az300EPq37+/YmNjlZGRoXXr1plusF2F6jWVTh/67cRQaU+ovR8A9kXWmxPs66kBmMO/XXSE4aJ65cqVKiwsVHFxsbZs2aKMjAzl5ubq4MGDrS4/a9Ys/epXv9Izzzyjjz/+WFOnTtUPfvCDdm/AjbZxmyA4vbB2ekcNEA7IegSK0zPNzti2QOAYLqoXLFigKVOmKD8/X0OGDNHixYvVuXNnLVu2rNXlX3nlFT3wwAO6+uqr1a9fP02bNk1XX321nnrqqQ43HmhPKIdJKL83f3N6RxTzQCAYyHoEAllmPadtY4aAw6kMFdWNjY3avHmzcnJyvllBZKRycnJUXl7e6nMaGhoUG+s9NLpTp07auHGjieYikJx8RtFpIQLnC9VLQBB+yHrnoiBBazgmAqxnqKg+fPiwmpqalJiY6PV4YmKiqqqqWn1Obm6uFixYoM8++0xut1sbNmzQ66+/rgMHDrT5Og0NDaqtrfX6AXwVLuHh5Pfp5A4bINSR9QgEJ2cYQhfXVcMsy2f/fvrpp3Xuuedq8ODBiomJUUFBgfLz8xUZ2fZLl5SUKD4+3vOTlpZmdTMdwenDWOF/HJQAsAOy/gQOyH1DdgWek7Z5qI+44NZ+oclQUZ2QkKCoqChVV1d7PV5dXa2kpKRWn9OzZ0+tXr1a9fX12rNnj7Zu3aouXbqoX79+bb7OzJkzVVNT4/nZu3evkWbCD5x6JtFJoeEv4fiejaJDCvAdWe9MTilEyKzgYdsD1jFUVMfExGj48OEqLS31POZ2u1VaWqrs7Ox2nxsbG6vU1FQdP35cf/jDH/T973+/zWVdLpfi4uK8fuyMayntgbBwFqd23AChLlSznrNDAACrGB7+XVhYqCVLluill17SJ598omnTpqm+vl75+fmSpLy8PM2cOdOz/LvvvqvXX39dO3fu1Ntvv60rr7xSbrdb9957r//ehY+cPIsuZ9rQHjoUAPiTk7Me9kVWBR+fAWCNaKNPGD9+vA4dOqSioiJVVVUpMzNT69at80xoUllZ6XUN1bFjxzRr1izt3LlTXbp00dVXX61XXnlF3bp189ubgH858QwiIXFiG+ze1zPYzQhrNQNiFb/9WLCbAXQYWe8sThj6TU7bhxOOFxp6N8pVGROU167rE6Gue5qD8tpwLsNFtSQVFBSooKCg1b+VlZV5/X7ZZZfp448/NvMygE8I6m84IShP1dCnUa49gQnNut4R6lpJSAK+IuuNYZIyOInTjhcAu7N89m/AShTUAADYG1kNINRRVDtAIK+nduLQb3jj4CX0OHk+CADhjUyyL7t/Nk64rAE4iaIajmX3MAgmJ20bOnIAwDwKD3SEk44XADujqO4gbqcVHIQAzGAWfQAIHLIaaB23+As9FNXwcMoZQ0LaN2wnAECwkEHOYefPKlgjMZh4EEZRVNscZ9bQEXYOylM5pUPHF4xeAezH6rNCwToAt+vQb6dkD77BZwZ0DEU1HIWdvnFsMwAAAMA6FNWQ5IwzhRSH8AdGfwCAdchq5+KzA8wLm6KaW9I4Gzv6jnHC9nNCxw4A2IUdh347IWvQPjt+hnb8rgOnC5ui2gpWXzvJGTX4kx2DEgAA2AvHC4BxFNWw/RlCdu7+w7Y8gQ4rAPAv8gWhhhnAYQRFNWyNkA4vdu/gAYDTBePA227DYcnq0GO3z9Ru33ngdBTVQJixW1CGIisuDWFeCAB2RKaELj5ba1l9qz8EFkW1TQVqeKqdzwyyM7cO2xYAAADwD4pq2BJFn/Xsuo0D1dHDddVAeOBskHXsmiPwHzt9xgwBh51RVMN27LQDBwDATuxSWJDV4YPPGjgzimqTrL6dViDYeeg3AoOgBAAAZxKuxwvMAA5fUVTbUDgPSw3XnXYw2XGb0+EDwAnC9YDbjrmB8GCXkRrA6cKiqGbWXGcgpBFo4dyBBcB57FBQkNXhi88eaFtYFNVoyW5nAtlRBxfb3/9C4RIRADgVWQG+A0DrKKoBSLJfUNqt4wcAAABoDUW1zYTjcFS7FXPhjM8CAOwp2EO/yQecFOzvQrD/LfgTt/wLHRTVYchOZwCDvWMGnNSRxfwQAIKBrMbpwuk7Ea4TEsIYimoTuFbSP8Jph+wkdvpc7NQBBMB5rDwLxIE2wp2djheAYKOoBtACQQkAkMgD2FMoDQFHaKCotpFADEO1y5k/Qtr++Iw6jlEtAPwhWAUEOYAz4TsCnEBRjYBjBwwjAtER5KTrqgEAsBOO6wCKagDtICgBIDyx/4fdMQQcdmKqqF60aJHS09MVGxurrKwsbdq0qd3lFy5cqEGDBqlTp05KS0vT3XffrWPHjplqcKgKl6HfhLTz8JnBCd566y2NGzdOKSkpioiI0OrVq31+7j/+8Q9FR0crMzPTsvY5EVlvD8EoHNjvw6hQ/84wMaG9GMmnr7/+Wg899JD69++v2NhYZWRkaN26dR1aZ2sMF9UrV65UYWGhiouLtWXLFmVkZCg3N1cHDx5sdflXX31V999/v4qLi/XJJ59o6dKlWrlypR544AGjL20Kt6Cxj1Df4YayYH92dugQgr3V19crIyNDixYtMvS8I0eOKC8vT9/97nctapkzOS3rAy2UD7CDvb+Hc/HdQSAYzadZs2bpV7/6lZ555hl9/PHHmjp1qn7wgx/o3//+t+l1tsZwUb1gwQJNmTJF+fn5GjJkiBYvXqzOnTtr2bJlrS7/zjvvaNSoUbrxxhuVnp6uK664QhMmTDBc/dsFEw+Zw44Wdsd11c521VVX6ZFHHtEPfvADQ8+bOnWqbrzxRmVnZ1vUMmcK96wPV2Q1OorvkHFW3vovFBnNp1deeUUPPPCArr76avXr10/Tpk3T1Vdfraeeesr0OltjqKhubGzU5s2blZOT880KIiOVk5Oj8vLyVp9z8cUXa/PmzZ5g3blzp9auXaurr77ayEujg4J5po8dbGjgc0SoefHFF7Vz504VFxcHuym2QtbbB9eMAu3j30h4MZNPDQ0Nio31PinaqVMnbdy40fQ6WxNt5I0cPnxYTU1NSkxM9Ho8MTFRW7dubfU5N954ow4fPqxLLrlEzc3NOn78uKZOndrukLCGhgY1NDR4fq+trTXSTMfhDBmcIv2cQ9q9r2dQXruhT6Nce2KC8todUTMgVvHb/XddaeOgVMVs2++39dnd6ft/l8sll6vjvfqfffaZ7r//fr399tuKjjYUhSEvVLKesz/hZ2xS699Pp9lQNTjYTeiQYB4rwJl8zXoz+ZSbm6sFCxbo0ksvVf/+/VVaWqrXX39dTU1NptfZGsuPJMrKyjRv3jw999xzysrK0vbt2zVjxgw9/PDDmj17dqvPKSkp0dy5c61uGgKAs5uhh7DE6f5WPVDRR/1bwByvb5D0F6WlpXk9XlxcrDlz5nRo3U1NTbrxxhs1d+5cDRw4sEPrwglkvfM5Pa9DpaCWTrwXpxfWCD1Oy/qTnn76aU2ZMkWDBw9WRESE+vfvr/z8fENDu31hqKhOSEhQVFSUqqurvR6vrq5WUlJSq8+ZPXu2br75Zk2ePFmSdMEFF6i+vl4//elP9eCDDyoysuUI9JkzZ6qwsNDze21tbYuNDd8Fa+i30wMabQvVwrqud4S6VjYHuxk4xd69exUXF+f53R9nqevq6vTee+/p3//+twoKCiRJbrdbzc3Nio6O1l/+8hd95zvf6fDrOBVZD6cJpYL6JApr3zX0bpSr0vqRbHV9ItR1D8cIVvA1683kU8+ePbV69WodO3ZMX3zxhVJSUnT//ferX79+ptfZGkPXVMfExGj48OEqLS31POZ2u1VaWtrmJC9fffVVizCNioqSJDU3t/7FdLlciouL8/qBs1BQwwrMAh5+Ts8CfxTVcXFx+uCDD1RRUeH5mTp1qgYNGqSKigplZWX5oeXORda3L1AzfwfyWlEnZ3YoFtQnOfm9Ofk7hcDzNevN5NNJsbGxSk1N1fHjx/WHP/xB3//+9zu8zlMZHv5dWFioiRMnasSIERo5cqQWLlyo+vp65efnS5Ly8vKUmpqqkpISSdK4ceO0YMECXXjhhZ4hYbNnz9a4ceM8gesUVsz8zfXUcKpQPVsNZzp69Ki2b9/u+X3Xrl2qqKhQjx491Lt3b82cOVP79+/Xyy+/rMjISA0dOtTr+b169VJsbGyLx8NVOGc9nMPJRaevOGMNeDOaT++++67279+vzMxM7d+/X3PmzJHb7da9997r8zp9YbioHj9+vA4dOqSioiJVVVUpMzNT69at81zcXVlZ6dVbPWvWLEVERGjWrFnav3+/evbsqXHjxunRRx81+tIwIRhn9uidDB+hWFgzBNyZ3nvvPY0ZM8bz+8lhxRMnTtTy5ct14MABVVZWBqt5jkPWA/bh1MI6FI8REHxG8+nYsWOaNWuWdu7cqS5duujqq6/WK6+8om7duvm8Tl9ENLc1LstGamtrFR8fr5xekxUdaeyaicZBqX5rhxPPVAe6qKagDk+BDk2rZwH3d1Htz9m/JZma/fu4u1FvHvy1ampq/DbM9uS+edSfChR9tv8nL/nH95/1a3thbye/T5deUqTo6I7nrVWzfwdi+DdDv88sHM5Sn8qJRbUU2OODQFxXbdU11fE7Gs68kA+OHz+mtzY+RNYHgaFrqp3GnwU1zsypwQxYzYoOOQAIV+FWUEvh+Z7DCbcAdL6QLqrtLpTOUlNQhzc+fwDhIFCTlKFt4VxcOvG9h9rxAfsAtIWiGoBfBDI4re4wYgJBwNmcfNaHod9tc2JR6W9sg7YF8t8OcDrDE5UBp3NCKIdKCNn9miomJQEAwFpOnbgMCGWcqfaRv6+JDJWh3xTUOJ0TvhNOx3wRADrKaftqsty5nPZdA8ygqIZpTthJhloIh9r76Yhg3C4OABB4ZF9LbJPWMQQcwUJRDTiME4LUCR0uZ+Lv0STMAA7gTCgIWnJC5gWLk7ZNKBwXnMRkZWgNRXUICsQZPCfsHJ0UNkY54b054TsCAL4KpQNpp+yfnZB1wcY2AuyBojoInD6zsBPCOBxCxgnv0ervCkPAASA0OSHj7IJtBQQfRTUMoaAGAMD/AjX02wk5DuOccOwTqO+eUy+jcPKtABHCRbU/Z8d10rWQVp65c0IQOyFU/MkJ79cJ35u2OH1UCQA4kROyDQBOFbJFNRAunHDwYWVhzRBwAKfibI+zOSHT7IptBwQPRXWAOfXMlxPONoZzmDjhvTvhOwQAoczu+2EnZJnd2X0bhsoQ8FCauBD+QVEdQqw6Y2f3EJbsHyKBwDawhj87wpx0KQmAb1h9AO3Ua0D9iQzzH7YlEHgU1XA8wuMbdt8WVnXQhOMQcH/OGwEAwWT37HIiO29TJ5ysAYyiqIaj2Tk0AACwCwqZ8BPux0iMAEEgUVQHkJXXU4fj0O9wD4u2sF0AwDnC/cCfzAIQCiiqz4BrIOFEdj5IceIQcKdOMAgAkn07yO2cVaHCrtvYrt9JwCyKajiSXUPCTthGAAC7IqMCh21tDSsmMOSWgM5FUR0Cwm3oN+HgO7aV/TD6BbCOFQekVs78Ha5Dv8mmwAvXbR6u/8YQeCFZVNtxVlyGj/pHuIZCqHHiEHAAcCK7dpAj8Ox2DMV3E6EkJItqdJwdd3R2CwOnYLv5Bx1jANBxZBKAUERR7XCcmYMv7HgQY8eOGwCAdeyYReHGbp9BII4FGAKOQKCobgfXPtqH3ULAicJlG9LRBMAJAnGgb6fOy3DJICfgs/AfK+dcgLNQVAeA04aN2imEJXb+/sS2DC12nD8CCCUcMPsH2WM/fCaAf1FUO1g4nJFjpx/a7NaBAwChwi77V3Lcvuzy2djluwp0BEU1bMsuO/tQEw7b1aoOJ3+NOuHSEiC8cY0n7CIcjgkk/s3BehTV8EJvYXgIlxAFAAQPWQMYV9PfFewmwARTRfWiRYuUnp6u2NhYZWVladOmTW0ue/nllysiIqLFzzXXXGO60Qj9od8EsfXsso3pyAHsyWlZz4HoN+ywX7VLxuDM7PBZ2eE7C3SE4aJ65cqVKiwsVHFxsbZs2aKMjAzl5ubq4MGDrS7/+uuv68CBA56fDz/8UFFRUfrRj37U4cY7gdMmKbMDO+zcw0Uob+tQ73gCrETWoyNCOVtCVTh8ZlYNAWdCQ0gmiuoFCxZoypQpys/P15AhQ7R48WJ17txZy5Yta3X5Hj16KCkpyfOzYcMGde7c2bKg9ddsuOF4zaMdegnDYaduN2xzAKeze9YHglUHyqF+bSeZ4lzB/uzscBwKmGWoqG5sbNTmzZuVk5PzzQoiI5WTk6Py8nKf1rF06VLdcMMNOvvss9tcpqGhQbW1tV4/+EaonoEL9s4cweOkIGX0CUIdWe9swdyfkuPOx2cImGOoqD58+LCampqUmJjo9XhiYqKqqqrO+PxNmzbpww8/1OTJk9tdrqSkRPHx8Z6ftLQ0I82ECcEuatiJB1eobn87d0CF42gYOANZD4S3UD0mkEJ/pAiCJ6Czfy9dulQXXHCBRo4c2e5yM2fOVE1Njedn7969AWohEL6CHaLB7tgB4B9kfdusPqDnLDWcjmMBOJWhojohIUFRUVGqrq72ery6ulpJSUntPre+vl4rVqzQpEmTzvg6LpdLcXFxXj9OZMUwUTufeTOLILYPPgvn8dc8EsBJZD2MIjtCD58pYIyhojomJkbDhw9XaWmp5zG3263S0lJlZ2e3+9zXXntNDQ0Nuummm8y1FJahZxunCrXPJBQ7ogArkfUwItQyA98I1mfrxLPVzAAOw8O/CwsLtWTJEr300kv65JNPNG3aNNXX1ys/P1+SlJeXp5kzZ7Z43tKlS3XdddfpW9/6VsdbbTGudQwMghinc0qQMlkZQl04ZH17rDhADsWh3+R46AvFz9gJ11XX9HcFuwkwKNroE8aPH69Dhw6pqKhIVVVVyszM1Lp16zwTmlRWVioy0rtW37ZtmzZu3Ki//OUv/mk1AEuNTdqqDVWDg90MAEHitKznABQAEEyGi2pJKigoUEFBQat/Kysra/HYoEGD1NzcbOalcAorhrEG68xgKPZ8AkAoIevRHnI8fNDRDpxZQGf/BuAcHDABAIBgcMrlYMBJFNUIOIo1tMeKIGWyMgChjAIEVuPYDWgfRTWANhGi1mJSRABORDbA6ZwwWRmchaLaIswODAAAJA7gERroTAHaRlEdxrj9BnzBZwYAOIlMAICWQqqobhyUGuwmWIZrQgEAABAunDZXgBX3t4dzhFRRDXujdxu+ckKQcokHEJqcdmAcyP0lOQ6+A0DrKKoBnFEohCijPQAAAGCF6GA3AOEhFIoyAPb11ltvaf78+dq8ebMOHDigP/7xj7ruuuvaXP7111/X888/r4qKCjU0NOj888/XnDlzlJubG7hGAw7ihBy/Pm5LsJvgF6tqvx3sJrRrbNJWbagaHOxmIIwtWrRI8+fPV1VVlTIyMvTMM89o5MiRrS57+eWX6+9//3uLx6+++mqtWbNGknTLLbfopZde8vp7bm6u1q1b53ObOFN9mnC5xY0ThtfCXpxwQIXwVV9fr4yMDC1atMin5d966y2NHTtWa9eu1ebNmzVmzBiNGzdO//73vy1uKQC0L1Q6B+yOWfmdaeXKlSosLFRxcbG2bNmijIwM5ebm6uDBg60u//rrr+vAgQOenw8//FBRUVH60Y9+5LXclVde6bXc7373O0Pt4ky1Azh92CrFGMxIP+eQdu/rGexmwCGuuuoqXXXVVT4vv3DhQq/f582bpz/96U/685//rAsvvNDPrUM448A9MChEAytQZ6s5FsDpFixYoClTpig/P1+StHjxYq1Zs0bLli3T/fff32L5Hj16eP2+YsUKde7cuUVR7XK5lJSUZLpdnKm2ABMYIVTRQYJQ5Xa7VVdX1yJ8ATsL1Kgz9v2BRycB0FJjY6M2b96snJwcz2ORkZHKyclReXm5T+tYunSpbrjhBp199tlej5eVlalXr14aNGiQpk2bpi+++MJQ2yiqAcDhQvl2grW1tV4/DQ0NlrzOk08+qaNHj+rHP/6xJesHYB0K0OCgs8VaNf1dwW5CwPia9YcPH1ZTU5MSExO9Hk9MTFRVVdUZX2fTpk368MMPNXnyZK/Hr7zySr388ssqLS3V448/rr///e+66qqr1NTU5PN7YPh3GOL2G+gIJ09Q0tCnUa49McFuRsip/DxBkZ38Ox+F+7/HJElpaWlejxcXF2vOnDl+fa1XX31Vc+fO1Z/+9Cf16tXLr+sGnI4cD57r47bYftIyhA+nZ7104iz1BRdc0GJSsxtuuMHz/xdccIGGDRum/v37q6ysTN/97nd9WjdFNQAEUc2AWMVvPxbsZtjW3r17FRcX5/nd5fJvz/2KFSs0efJkvfbaa17DyQA4A2epAefzNesTEhIUFRWl6upqr8erq6vPeD10fX29VqxYoYceeuiM7enXr58SEhK0fft2n4tqhn/DMvRuo6PsPks98ydYLy4uzuvHn0X17373O+Xn5+t3v/udrrnmGr+tF4EVTkMkA40cDz67dxoE4jti5bGAvycSrOvDcYEZvmZ9TEyMhg8frtLSUs9jbrdbpaWlys7Obvc1XnvtNTU0NOimm246Y3v27dunL774QsnJyT6/B85U25zTZ/62C7uHkq/sMgzMyUPAEZqOHj2q7du3e37ftWuXKioq1KNHD/Xu3VszZ87U/v379fLLL0s6MeR74sSJevrpp5WVleW5FqtTp06Kj48PyntA6LFy5m+7dzoGQqhkOwDfFRYWauLEiRoxYoRGjhyphQsXqr6+3jMbeF5enlJTU1VSUuL1vKVLl+q6667Tt771La/Hjx49qrlz5+qHP/yhkpKStGPHDt17770aMGCAcnNzfW4XZ6phCTv1bhO6QOh77733dOGFF3puh1VYWKgLL7xQRUVFkqQDBw6osrLSs/wLL7yg48ePa/r06UpOTvb8zJgxIyjthz1wlukbdsrxcGf34xi+Kwik8ePH68knn1RRUZEyMzNVUVGhdevWeSYvq6ys1IEDB7yes23bNm3cuFGTJk1qsb6oqCi9//77uvbaazVw4EBNmjRJw4cP19tvv21odBxnqsMMPdvOZqdJS5x6tprJykLT5Zdfrubm5jb/vnz5cq/fy8rKrG0QAEvZvdAEYJ2CggIVFBS0+rfW8n3QoEFtHiN06tRJ69ev73CbOFMNv7NTjyWhCwCAMXbKcZxg9+MZvjMIdxTVfsbERbCanYLV6ROUAACcx045iMDgWAB2R1F9ipoB/r33GoKL0AUAhCIKjPBk9+MazlYjnFFU25gTZ/5mhxoYdgpWPnMACB1236fbKf8QGqycpR/hg6I6jIRTzzahCwAIJg7UYQWObwB7CpmiunFQarCbEPbs3rsdaghW8/w5CoR5FACECrvnOLlnf3b/DgFWCZmiGjiJ0A08q0M0nEZZAABaItu/wbYA7IeiGn5Bz2RwEKyhgUkSAfiKTkaEK6d89+v6MIItHJkqqhctWqT09HTFxsYqKytLmzZtanf5I0eOaPr06UpOTpbL5dLAgQO1du1aUw0G2kORGTx0rAChJdyyPtwPhO28DyfbW7LzNrHzdwmwiuGieuXKlSosLFRxcbG2bNmijIwM5ebm6uDBg60u39jYqLFjx2r37t1atWqVtm3bpiVLlig1lWug2+Pvmb+t7N1j5xlcdg5WAM5E1gMIJ3adWLCmvyvYTYCPDBfVCxYs0JQpU5Sfn68hQ4Zo8eLF6ty5s5YtW9bq8suWLdOXX36p1atXa9SoUUpPT9dll12mjIyMDjfebpiwKLgoLoPPyg4Wpwz7AkIBWd8xdj1AdyKyvW1sG8A+DBXVjY2N2rx5s3Jycr5ZQWSkcnJyVF5e3upz3njjDWVnZ2v69OlKTEzU0KFDNW/ePDU1NXWs5QA8CFbjnHgfeCAQyHoAAIwxVFQfPnxYTU1NSkxM9Ho8MTFRVVVVrT5n586dWrVqlZqamrR27VrNnj1bTz31lB555JE2X6ehoUG1tbVeP7Anuwz9pqhEuOO2gvAXJ2V9uA2NDLcRO2T7mdl1G1l1fBhu/wbgHJbP/u12u9WrVy+98MILGj58uMaPH68HH3xQixcvbvM5JSUlio+P9/ykpaVZ3cyQxg4oPNglWO3S0QIgcMh6Z2O/DQAdY6ioTkhIUFRUlKqrq70er66uVlJSUqvPSU5O1sCBAxUVFeV57LzzzlNVVZUaG1sffjlz5kzV1NR4fvbu3WukmQgQu4SwXYpJhC/mU0AoIethB2S779hWQPAZKqpjYmI0fPhwlZaWeh5zu90qLS1VdnZ2q88ZNWqUtm/fLrfb7Xns008/VXJysmJiYlp9jsvlUlxcnNeP1ex0n1iu9YRZdglWhn0BzhXKWQ8AgBUMD/8uLCzUkiVL9NJLL+mTTz7RtGnTVF9fr/z8fElSXl6eZs6c6Vl+2rRp+vLLLzVjxgx9+umnWrNmjebNm6fp06f7710g4DhLDQChi6xHMJHtxtlxm9nlWNFXzNqPjog2+oTx48fr0KFDKioqUlVVlTIzM7Vu3TrPhCaVlZWKjPymVk9LS9P69et19913a9iwYUpNTdWMGTN03333+e9dAPC4Pm6LVtV+O9jN0NikrdpQNTjYzWhXQ59Gufa0fhYNCGdkvXlWHZhbNVLHaYUPkH7OIe3e1zPYzWhXXZ8Idd3THOxmIIAMF9WSVFBQoIKCglb/VlZW1uKx7Oxs/fOf/zTzUuggK0LYLgFsx15ZAAgVZD2CgWw3zy6d6kA4snz273DBREWwE7sclNilA8YJ7DSvAxBu6vqQ4QAA8yiqbYZJynxjl6IRgcdkZQAQmsj2jmMbAsFBUQ1DOPPoHAQrAKA9ZDqsxncM4YKiGo5DsegsBCoAdEw4jNAh2/0nHLalVf8mmAEcZoVEUd04KDXYTbAlf+9wKI6cJxyCtaP8dckF8yoA4IAcAMJTSBTVCB8Uic7k7w6ZcDhrAwDhgmz3P7YpEFgU1fAJZ6mdi2AFAADBwjFkx9T0dwW7CfABRbWNMPN3+ygOnY1QBQD7sMs+mWy3DtsWCByKanF/WAAAAMBJuBQMdkJRDQAAAACASRTVAAAAAACYRFENAAAAAIBJFNUAAAAAAJhEUQ0AAAAAgEkU1QAAAAAAmERRDQAAAACASRTVAAAAAACYRFENAAAAAIBJFNUAAAAAAJhEUQ0AAAAAgEkU1QAAAAAAmERRDQAAAACASRTVAAAAAACYRFENAAAAAIBJFNUAAAAAAJhEUQ0AAAAAgEkU1QAAAAAAmERRDQAAAACASaaK6kWLFik9PV2xsbHKysrSpk2b2lx2+fLlioiI8PqJjY013WAAAE731ltvady4cUpJSVFERIRWr159xueUlZXp29/+tlwulwYMGKDly5db3k4nIesBAHZkJJ8k6ciRI5o+fbqSk5Plcrk0cOBArV27tkPrPJ3honrlypUqLCxUcXGxtmzZooyMDOXm5urgwYNtPicuLk4HDhzw/OzZs8foywIA0Kb6+nplZGRo0aJFPi2/a9cuXXPNNRozZowqKip01113afLkyVq/fr3FLXUGsh4AYEdG86mxsVFjx47V7t27tWrVKm3btk1LlixRamqq6XW2xnBRvWDBAk2ZMkX5+fkaMmSIFi9erM6dO2vZsmVtPiciIkJJSUmen8TERKMvCwBAm6666io98sgj+sEPfuDT8osXL1bfvn311FNP6bzzzlNBQYGuv/56/eIXv7C4pc5A1gMA7MhoPi1btkxffvmlVq9erVGjRik9PV2XXXaZMjIyTK+zNYaK6sbGRm3evFk5OTnfrCAyUjk5OSovL2/zeUePHlWfPn2Ulpam73//+/roo4/afZ2GhgbV1tZ6/QAA4C/l5eVeWSZJubm57WZZuAhU1gMAYISZfHrjjTeUnZ2t6dOnKzExUUOHDtW8efPU1NRkep2tMVRUHz58WE1NTS16nxMTE1VVVdXqcwYNGqRly5bpT3/6k37zm9/I7Xbr4osv1r59+9p8nZKSEsXHx3t+0tLSjDQTABAiTu9gbWho8Mt6q6qqWs2y2tpa/fe///XLazhVoLKeDnQAgOR71pvJp507d2rVqlVqamrS2rVrNXv2bD311FN65JFHTK+zNdE+L2lSdna2srOzPb9ffPHFOu+88/SrX/1KDz/8cKvPmTlzpgoLCz2/19bWUlgDgE3FVMYoKjbGr+tsOuaWpBb7/uLiYs2ZM8evr4WOM5P1JSUlmjt3bqCaCADogJi9zsx6t9utXr166YUXXlBUVJSGDx+u/fv3a/78+SouLvbLa0gGi+qEhARFRUWpurra6/Hq6molJSX5tI6zzjpLF154obZv397mMi6XSy6Xy0jTAAAhaO/evYqLi/P87q9sSEpKajXL4uLi1KlTJ7+8hlMFKuvpQAcASL5nvZl8Sk5O1llnnaWoqCjPY+edd56qqqrU2Njol8yTDA7/jomJ0fDhw1VaWup5zO12q7S01KuHuj1NTU364IMPlJycbOSlAQBhKC4uzuvHX0V1dna2V5ZJ0oYNG3zOslAWqKx3uVwtPl8AQPjxNevN5NOoUaO0fft2ud1uz2OffvqpkpOTFRMT45fMk0zM/l1YWKglS5bopZde0ieffKJp06apvr5e+fn5kqS8vDzNnDnTs/xDDz2kv/zlL9q5c6e2bNmim266SXv27NHkyZONvjQAAK06evSoKioqVFFRIenELbMqKipUWVkp6cRZ0by8PM/yU6dO1c6dO3Xvvfdq69ateu655/T73/9ed999dzCabztkPQDAjozm07Rp0/Tll19qxowZ+vTTT7VmzRrNmzdP06dP93mdvjB8TfX48eN16NAhFRUVqaqqSpmZmVq3bp3n4u7KykpFRn5Tq//nP//RlClTVFVVpe7du2v48OF65513NGTIEKMvDQBAq9577z2NGTPG8/vJYcUTJ07U8uXLdeDAAU+BLUl9+/bVmjVrdPfdd+vpp5/WOeeco1//+tfKzc0NeNvtiKwHANiR0XxKS0vT+vXrdffdd2vYsGFKTU3VjBkzdN999/m8Tl9ENDc3N/vvbVqjtrZW8fHxyuk1WdGRLS+QbxyU2sqzfFczILZDz5ekut4RHV5HQ5/GDq/jVOnnHPLbusYmbfXbusy6Pm5LsJvgWKtqvx3sJkiSNlQN9tu6du/r6bd1ufb4Z+KNrpUd253Gbz/WoefHbNvf5t+Ouxv15sFfq6amxm/DbE/um/vPnKeo2I7vR0/VdOyYdpQ84Nf2wt5Ofp8uvaRI0dG+fZ9q+vvncoC6Ph3PcElq6O3fHD/Jn3l+kh1yXSLbrRaK+X8qfx4LnOSq9NMxwR7/lVjxO3y788Xx48f01saHrMn6ByzK+nmhkfWGh38DAAAAAIATKKoBAAAAADCJohoAAAAAAJMoqgEAAAAAMImiGgAAAAAAkyiqAQAAAAAwiaIaAAAAAACTKKoBAAAAADCJohoAAAAAAJMoqgEAAAAAMImiGgAAAAAAkyiqAQAAAAAwiaIaAAAAAACTKKoBAAAAADCJohoAAAAAAJMoqgEAAAAAMImiGgAAAAAAkyiqAQAAAAAwiaJaUvz2Y8FuAgAAAAAf7d7XM9hNADwoqm3EtScm2E2wtVW13w52E9ABG6oGB7sJAID/xy77ZLLdOmxbIHAoquETu4QvjCNUAQBAsHAM2THxOxqC3QT4gKIajkKB6Ez+DlSGfAFA6CDb/Y9tCgRWSBTVMdv2B7sJtuTvwoOeRuchVM/MX5dddK1s9st6ADiXq5LLuAAgHIVEUY3wQqHoLHTGAEDHhMPoHLLdf8JhW1r1b4KOMZhFUQ1DKJCcIxxCFQBgHpkOq/EdQ7igqLYZZgD3DQVj+AqHMzYAEI7I9o5jGwLBQVHtJ1xPCTuxS6jSQ+27+O3Hgt0EIGx13UOGAwDMo6gOcVac1bNLoWSXwhEAAPgH2W4e2w4IHlNF9aJFi5Senq7Y2FhlZWVp06ZNPj1vxYoVioiI0HXXXWfmZQH4wC6hapfOl/ZwuQXQNrLeHKsmOrLq0hcn7KuBUznhMjBGv4Qfw0X1ypUrVVhYqOLiYm3ZskUZGRnKzc3VwYMH233e7t27dc8992j06NGmGwv7sEsI26WABIBQQtYjmMh24+y4zexyrOgrZv5GRxguqhcsWKApU6YoPz9fQ4YM0eLFi9W5c2ctW7aszec0NTXpJz/5iebOnat+/fp1qMFWsdP1jJw9g1l2CVWrgtQJvdNAKAjVrAcAwAqGiurGxkZt3rxZOTk536wgMlI5OTkqLy9v83kPPfSQevXqpUmTJvn0Og0NDaqtrfX6gf3YpQfSLoUkwhcTFSKUkPWwA7Ldd2wrIPgMFdWHDx9WU1OTEhMTvR5PTExUVVVVq8/ZuHGjli5dqiVLlvj8OiUlJYqPj/f8pKWlGWkmTsPZvfBgl1C1S2cLAHPI+vDDfhsAOsbS2b/r6up08803a8mSJUpISPD5eTNnzlRNTY3nZ+/evRa2Eh1hlyC2S0EJBEvMtv3BbgLCVDCzPn5Hg+HnOFm4dZKT7Wdm123EZWAIN9FGFk5ISFBUVJSqq6u9Hq+urlZSUlKL5Xfs2KHdu3dr3LhxnsfcbveJF46O1rZt29S/f/8Wz3O5XHK5XEaaBoQ1u4aqnTF3AdA6sh4AAGMMnamOiYnR8OHDVVpa6nnM7XartLRU2dnZLZYfPHiwPvjgA1VUVHh+rr32Wo0ZM0YVFRUhN9SL6yqDi8Iy+KwcuUDvNBAYZH3HMYuw/5DtbWPbAPZh6Ey1JBUWFmrixIkaMWKERo4cqYULF6q+vl75+fmSpLy8PKWmpqqkpESxsbEaOnSo1/O7desmSS0ehzfXnhg19Gn02/p27+up9HMO+W19p9pQNVhjk7Zasm6cGaEKwN/I+vBDliOc2bUjLNwucXEyw0X1+PHjdejQIRUVFamqqkqZmZlat26dZ0KTyspKRUZaeqk20KZVtd/W9XFbgt2MsGSX6+sBdFw4Zn3XPc2q6xMR7GagFWR7S3buUOd4AOHIcFEtSQUFBSooKGj1b2VlZe0+d/ny5WZeEjZHD3dw2DlU4bv47ceC3QSgBbLenqwceWZnFNZwymVgXfdwOWg4Cq1uZkAUmsFgda+0U4IUAACrcZwD2E/IFNXcTib4GO4TWISqef6c+ZsJCgGECrvnOLlnf3b/DgFWCZmiGmcWTmf7CF4AQDDZdeIjOBvHN4A9UVTbmBPvo0sPZWDYKVT5zAEgdNh9n26n/ENooAMM/kBRfQomCwotBC8AIBSF08gzfMPuxzV275ABrERR7WdcXwmr2SlUAxGgHDwCAE5lpxxEYHAsALujqIbf2amnkuAFAMAYO+U4TrD78QzfGYQ7iuowQ0+fs9kpVJ0aoE6cqwAA4M1OeQgAFNWwhJ0KLoIXCB+LFi1Senq6YmNjlZWVpU2bNrW7/MKFCzVo0CB16tRJaWlpuvvuu3XsGPNrhKuue7iE6yQ75Xi4s/txDN8VBJrRrD9pxYoVioiI0HXXXef1+C233KKIiAivnyuvvNJQm6INLY2Ac+2JUUOfxmA3w/HsHkhOQ4DCjlauXKnCwkItXrxYWVlZWrhwoXJzc7Vt2zb16tWrxfKvvvqq7r//fi1btkwXX3yxPv30U0+wLliwIAjvAKHGVRmjht7WZPjufT2Vfs4hS9btFKtqv63r47YEuxkAAsho1p+0e/du3XPPPRo9enSrf7/yyiv14osven53uVyG2sWZaliGwgsdZffLFZiY0F4WLFigKVOmKD8/X0OGDNHixYvVuXNnLVu2rNXl33nnHY0aNUo33nij0tPTdcUVV2jChAk+93jDHuJ3NAS7CSGLHA8+u58UcPqEpf6+nRajXaxnNOslqampST/5yU80d+5c9evXr9VlXC6XkpKSPD/du3c31C6KagAIIm7l177a2lqvn4aG1guoxsZGbd68WTk5OZ7HIiMjlZOTo/Ly8lafc/HFF2vz5s2eInrnzp1au3atrr76av+/EQCWsHvRCeDMrMx6SXrooYfUq1cvTZo0qc1lysrK1KtXLw0aNEjTpk3TF198Yeg9MPw7DAVyyNiGqsEam7Q1IK+FwHDymQsmKbNG133Niorxb+98U+OJ9aWlpXk9XlxcrDlz5rRY/vDhw2pqalJiYqLX44mJidq6tfV90I033qjDhw/rkksuUXNzs44fP66pU6fqgQce8M+bAEIAOR48dBjATrrudWbWb9y4UUuXLlVFRUWb7bjyyiv1P//zP+rbt6927NihBx54QFdddZXKy8sVFRXl03uhqAYAh4vZtj/YTbDM3r17FRcX5/nd6DVO7SkrK9O8efP03HPPKSsrS9u3b9eMGTP08MMPa/bs2X57HQDW4trq4HByJ7sThNOlLVZlfV1dnW6++WYtWbJECQkJbS53ww03eP7/ggsu0LBhw9S/f3+VlZXpu9/9rk+vRVFtga6VzarrHRHsZgB+R4Ai0OLi4ryCti0JCQmKiopSdXW11+PV1dVKSkpq9TmzZ8/WzTffrMmTJ0s6EaT19fX66U9/qgcffFCRkVwhBXsL1MgzzlYHHmepEU6syvodO3Zo9+7dGjdunOcxt9stSYqOjta2bdvUv3//Fs/r16+fEhIStH37dp+Lao4YHMDpQ1YpxGCG3Scpg73ExMRo+PDhKi0t9TzmdrtVWlqq7OzsVp/z1VdftSicTw7zam5mshn4h78nQkLrKEIDK1DHdhwL4FRGs37w4MH64IMPVFFR4fm59tprNWbMGFVUVLQYdn7Svn379MUXXyg5OdnntnGm+jTx24+pZkBssJthOW7FAaPoHIHdFRYWauLEiRoxYoRGjhyphQsXqr6+Xvn5+ZKkvLw8paamqqSkRJI0btw4LViwQBdeeKFn+Pfs2bM1btw4n6+hAgB/o4MgMOjwciYjWR8bG6uhQ4d6Pb9bt26S5Hn86NGjmjt3rn74wx8qKSlJO3bs0L333qsBAwYoNzfX53ZRVCMgGDoGwGrjx4/XoUOHVFRUpKqqKmVmZmrdunWeCU0qKyu9zkzPmjVLERERmjVrlvbv36+ePXtq3LhxevTRR4P1FgDbckKOU4wGBp3sCCajWX8mUVFRev/99/XSSy/pyJEjSklJ0RVXXKGHH37Y0LXdFNUAzigUAtTpl1HANwUFBSooKGj1b2VlZV6/R0dHq7i4WMXFxQFoGQAA8AcjWX+65cuXe/3eqVMnrV+/vsNt4ppqBEwoFGYIDCdcQ9W1kmtugVDUdY+z/m0Hcn9JjoPvANC6kCqqQ/m2MpxlAwAAQLhwQgf7qZzWIQf/CqmiGsYEY2dFD6fz8JkBAE4iEwCgJYpqizA0FAAASMwyjNBAhwrQNopqAG0iQK0Vv/1YsJsAAIaRDXA6OrrgbxTVCDjCGO2x4rIE5iQAEMqcdu0pnIdjN6B9FNUAWkWAAgCAYKCjCE5DUe0gVpxtC9ZOi4INAADnIsfDB581cGamiupFixYpPT1dsbGxysrK0qZNm9pc9vXXX9eIESPUrVs3nX322crMzNQrr7xiusEArEeAAnBS1sfvaAjYawEAcDrDRfXKlStVWFio4uJibdmyRRkZGcrNzdXBgwdbXb5Hjx568MEHVV5ervfff1/5+fnKz8/X+vXrO9x4qzB5UGBQuOF0Thnuxez+CHXhkPXtseJ+s1ZPjMRtMmGFUPyMnTBJGR2FzmO4qF6wYIGmTJmi/Px8DRkyRIsXL1bnzp21bNmyVpe//PLL9YMf/EDnnXee+vfvrxkzZmjYsGHauHFjhxsP/whmIROKO2unC7XPhEnKAOPIevgq1DID3wjWZ+uUDvZTWdERB2cxVFQ3NjZq8+bNysnJ+WYFkZHKyclReXn5GZ/f3Nys0tJSbdu2TZdeeqnx1jqMFWezQrFAIJDtg8/CeWK27Q92ExBiyHoYRXaEHj5TwJhoIwsfPnxYTU1NSkxM9Ho8MTFRW7dubfN5NTU1Sk1NVUNDg6KiovTcc89p7NixbS7f0NCghoZvhj3U1tYaaSYAE4IdoE7smQZCEVlvHVdljBp6N1q2/t37eir9nEOWrb89G6oGa2xS298PwBccC8CpAjL7d9euXVVRUaF//etfevTRR1VYWKiysrI2ly8pKVF8fLznJy0tLRDNDGvB3okFu6ALd6G6/e08soO5GxBqyHogNITqMYHkjOup4UyGiuqEhARFRUWpurra6/Hq6molJSW1/SKRkRowYIAyMzP1s5/9TNdff71KSkraXH7mzJmqqanx/Ozdu9dIM0OenQuFjgjlnTjaF+xOHSOYpAyhjqx3NuZJQUfwGQLmGCqqY2JiNHz4cJWWlnoec7vdKi0tVXZ2ts/rcbvdXkO+TudyuRQXF+f143Mb/XR9YTieRbJDYcPOPPDY5gBO5YSsDwSrJh4K9TNlZIpzBfuzs8NxKGCWoWuqJamwsFATJ07UiBEjNHLkSC1cuFD19fXKz8+XJOXl5Sk1NdXTO11SUqIRI0aof//+amho0Nq1a/XKK6/o+eef9+87samulc2q6x0R7GY4CtdlBU6wA9RKoTqiAwgEsh4dQY47TygfD5xkVYcWM39DMlFUjx8/XocOHVJRUZGqqqqUmZmpdevWeSY0qaysVGTkNyfA6+vrdfvtt2vfvn3q1KmTBg8erN/85jcaP368/95FGHLtiVFDH+smOwk2Atl6dglQeqYB+3Fi1sfvaFBNf1fAXs/Ogjlh2UnkuHPY4XiAYwE4neGiWpIKCgpUUFDQ6t9On5TkkUce0SOPPGLmZRAEdghiWM8OAQrA3sh6dBSFNWBc/I62L5uBfQVk9m/ADAo/a4TDdrVq6Le/JikLxzkbAHwj1K+rhnOEwzGBxL85WI+i2sHC4ZrRcNnZhyuGewGANeyyfyXH7csun41dvqtAR1BUB4DTbsFjt52bXXb6oYBtGVr8dbcDAK1jAiL/IHvsh88E8C+K6nYwRNM+2Pl3XLhsw3AYwQHA+QIxHNVOneThkkFOwGfhP3S84SSKaoejgIAv7BigdjrYAwBYz45ZFG7s9hkE4liA66kRCBTVaJUdCx67BYFTsN38w2mXcQCAHZFJAEJRSBbVdrzOkANy/yCMQ4NVnTaM3AAAb3bsJEdw2O0Yiu8mQklIFtXhxqpCwq47O7uFgp2xreyHuRoA61hxf1crr5kM12GpZFPghes2D9d/Ywg8imo4UriGgxFsIwCAXZFRgcO2toYVHW5WdAwiMCiqz4CzSnAiOweoE4d+c/kGACdj5Fn4sus2tut3EjCLojqArDwwD7ch4JJ9gyLY2C4A4BzhPjyVzAIQCiiq4WiEMQAAZ2bnTnJYI9yPkcK9wwqBRVENxwv30DiV3beFE4d+25Ud73IAAGbYPbucyM7blA4ehCKK6hASjkPAT7JzeAQK28Aa/rxsgzkaAGeycgZwiTNqEhnmT2xLIPAoqgPMqRMeUVjbmxPeuxO+QwAQyuy+H3ZCltmd3bdhoL6DVndUWd3RBuehqAYczu4BKlkbouE49BtA27gljbM5IdPsim0HBE/IFtX+vN7QSUM2rSww7N7DLYVfoDjh/Trhe9MWp44sAQAnc0K2AcCpQraohjWcUCARxgAApwnUddVOyHEY54Rjn1AZ+m0VRtk4G0V1EDj97JcTAtkJ4dJRTniPVn9XGPoNAKHJCRlnF2wrIPgoqkNQIAoNCuvgcsJ7c8J3BAB8FUoTEzll/+yErAs2thFgDxTVgMM4IUCdcsDWHn+PKHHS3AwAgsOpw1at5ITMCxYnbZtQOC44KZQ62OA/FNUwzQk7SCcFji9C7f10BEO/ASA8kH0tsU1aR8cUgoWi2kf+Pstk9XXVgSo4KKxxOid8J5zOn3c3ABCenLavJsudy2nfNcCM6GA3AM63e19PpZ9zKNjNaBdhHBgEJwAA1uKYBrAfzlQD8ItAFtRWj8Rw+gz9QLhz8q1pAjl81WkdoRSTbIP2MPQbwURRHUShMgRccl4ww7/4/AGEAyYoCr5wLiqd+N5D7fiAfQDaEtJFNdcdBlao7TgBf2HmbwDwHycWlx0Vju85nDh5dA1OCOmiGoFHYR1+Av2ZM+s3gFDHEPAzo8jEqRj6jWAzVVQvWrRI6enpio2NVVZWljZt2tTmskuWLNHo0aPVvXt3de/eXTk5Oe0uD/8KRgHi1ICGcaH4WXM9tbMZySdJOnLkiKZPn67k5GS5XC4NHDhQa9euDVBr7Y2sB+zBqR0IoXiMAHswkk+vv/66RowYoW7duunss89WZmamXnnlFa9lmpubVVRUpOTkZHXq1Ek5OTn67LPPDLXJcFG9cuVKFRYWqri4WFu2bFFGRoZyc3N18ODBVpcvKyvThAkT9Le//U3l5eVKS0vTFVdcof37nTc024ohnBzAw6kIS9iN0XxqbGzU2LFjtXv3bq1atUrbtm3TkiVLlJqaGuCW2084Zz2cw6nFphHh8B4BI4zmU48ePfTggw+qvLxc77//vvLz85Wfn6/169d7lnniiSf0y1/+UosXL9a7776rs88+W7m5uTp2zPfaL6K5udlQVZeVlaWLLrpIzz77rCTJ7XYrLS1Nd9xxh+6///4zPr+pqUndu3fXs88+q7y8PJ9es7a2VvHx8crpNVnRkcbOvDYO8u/BUc2AWL+uT5Lqekf4fZ2naujTaOn622L322yhY4JRVAdi5IUVHV3+7pAzM1/EcXej3jz4a9XU1CguLs4v7Ti5b868+VFFxfh339jUeEwVrzxoqL1G82nx4sWaP3++tm7dqrPOOsuv7Xe6YGb9pZcUKTq649+nmv6uDq+jLXV9rM3tkxp6By6/nZzZY5O2BrsJlnByQR3wS8MCNPzbqonK/HVN9fHjx/TWxofI+jP49re/rWuuuUYPP/ywmpublZKSop/97Ge65557JEk1NTVKTEzU8uXLdcMNN/i0TkNnqhsbG7V582bl5OR8s4LISOXk5Ki8vNyndXz11Vf6+uuv1aNHjzaXaWhoUG1trdcPzAvWNaicyQxdofrZMnLEfk7PgoaG1g88zOTTG2+8oezsbE2fPl2JiYkaOnSo5s2bp6amJkvei1OQ9XAaJxefbQnF92QVpxfUsDbrT9Xc3KzS0lJt27ZNl156qSRp165dqqqq8lpnfHy8srKyfM48SYr2eUlJhw8fVlNTkxITE70eT0xM1NatvvUS3nfffUpJSfFq+OlKSko0d+5cI02DTe3e19PRvd9oKVQLapgXt/OYog2lyZkdP37i7H5aWprX48XFxZozZ06L5c3k086dO/XXv/5VP/nJT7R27Vpt375dt99+u77++msVFxf75404EFkfnpye1xuqBofMGWsKathR3K4GRUf7d5TO8eMnimcrs146ceY5NTVVDQ0NioqK0nPPPaexY8dKkqqqqjzrOH2dJ//mCz8fBrXvscce04oVK1RWVqbY2LaHD8ycOVOFhYWe32tra1ts7FDStbLZ8iHggD8Es6B26qzfdhj67WR79+71GhLmcvlvSK/b7VavXr30wgsvKCoqSsOHD9f+/fs1f/78sC6qO8ouWR+/o8HSIeCwH4pRe6DzHUZZmfWS1LVrV1VUVOjo0aMqLS1VYWGh+vXrp8svv9xvr2GoqE5ISFBUVJSqq6u9Hq+urlZSUlK7z33yySf12GOP6c0339SwYcPaXdblcvl9Y4Y7156YoF1b7fTeb5xASCIY4uLifLrOykw+JScn66yzzlJUVJTnsfPOO09VVVVqbGxUTIwzO3I6iqy3D1dlTECvqwachltphQYrs146MUR8wIABkqTMzEx98sknKikp0eWXX+55XnV1tZKTk73WmZmZ6fN7MHRNdUxMjIYPH67S0lLPY263W6WlpcrOzm7zeU888YQefvhhrVu3TiNGjDDykrZjxQzg4YCCDHbH9dTOZiafRo0ape3bt8vtdnse+/TTT5WcnBy2BbVE1oczshodxXfIOH9NUhYOzObT6dxut+e67b59+yopKclrnbW1tXr33XcNrdPwLbUKCwu1ZMkSvfTSS/rkk080bdo01dfXKz8/X5KUl5enmTNnepZ//PHHNXv2bC1btkzp6emqqqpSVVWVjh49avSlTQm3oZJ2xo7WuYL92Tl16DcCy2g+TZs2TV9++aVmzJihTz/9VGvWrNG8efM0ffr0YL0F23Ba1gdaKE9YFOz9PZyL7w4CwWg+lZSUaMOGDdq5c6c++eQTPfXUU3rllVd00003SZIiIiJ011136ZFHHtEbb7yhDz74QHl5eUpJSdF1113nc7sMX1M9fvx4HTp0SEVFRaqqqlJmZqbWrVvnubi7srJSkZHf1OrPP/+8Ghsbdf3113utp60L0MNVIK6rDuYQ8JMYCu48hCScwmg+paWlaf369br77rs1bNgwpaamasaMGbrvvvuC9RZsg6y3j2AMASerYVSoHyuEckea0xjNp/r6et1+++3at2+fOnXqpMGDB+s3v/mNxo8f71nm3nvvVX19vX7605/qyJEjuuSSS7Ru3bp25wU5neH7VAdDR+5TLTnjXtWS9ferloJ3z+pTEdTOYZeQdOr9qSX7TFRm5X2qLx012y/3FT7V8ePH9NY/HvZre2Fv/r5PtRQa96qWAnu/6lOR1/BVMI4XAnk9tZVFtT+Hf1t5n2p/7ptPsqK9wWJ4+DfQUXYp1OAMTi6oAQAIdRzXARTVthKIA3u7XJvKDtj++Iw6jokNAfhDsGY4JgdwJnxHgBMoqgG0QEgCACTyAPbErbRgNxTVJnD2yT8Ianuy0+dil5EVAJzJylvVMHERwp2djheAYKOoDkN2KlTYISPYnHQ9NbcIBBAMZDVOF07fCTrQ4AuKaptx0gG+v4TTjtnu+CwAwJ6CPdyVfMBJwf4uBPvfgj9ZOZoGgUVRDUBS8EPydHYaUQEAAAC0haI6TNmtYLFbQRdu2P7+x9wLAEINWQG+A0DrwqKo5jpEZ2BHjUALx8stADiXHYa9ktXhi88eaFtYFNVOE84H+uywA8+O29xuIykAoDXhOoGRHXMD4cEOHUtAayiqTQqFoZ0ULuDACAAAnEm4Hi+Ea8cZjKOohu2E644bAIAzscuZOrI6fPBZA2dGUQ1bYgduPbtu40CNoAjnyyyAcMIta6xj1xyB/9jpM7ZLhxLQGopqmwrUAb+dh4DbaUceati2AAAAgH9QVANhhoLaelbMucBdDADYEZkSuvhsrcUomtBCUQ1bY4ceXuw8cgIAWhOMiYzsNgyWrA49dvtM7fadB05HUQ3bFzJ227E7GdvyBK6nBgD/Il8Qapj5G0ZQVHeA1bfV4sAf/sQBDwAAOBOOFwDjwqao5npEZ2MH3zFO2H52HzEBAHZix+GwTsgatM+On6Edv+vA6cKmqEb7nFDQ2HFHD+dhBAgAWIesdi4+O8A8imo4Cjt849hmAAAAgHUoqm2Os2roCKcU1E4YKeErq+daAGCc1beuCdaERnYdFuuU7ME3+MyAjqGohodTCht2/L5hOwEAgoUMcg47f1bB6jhi5m8YRVHdQZyVCg47BwDsi5EfABA4ZDXQOqtHzyDwKKrhWIR125y0bZwyQgIA7MiuQ8DhDE46XgDsjKLaAQJ5do0Cx/kIyNDDLQEBOBWZZF92/2zoMIKTUFTD0eweCAAAhDuyGkCoo6iG4xHW33DatgjkyAiupwZgJSY2gpM47XgBsDtTRfWiRYuUnp6u2NhYZWVladOmTW0u+9FHH+mHP/yh0tPTFRERoYULF5ptKwLEiUPACQe2gR0wcSFCCVnvHE4YJktG2YcTPotgfqfpIIMZhovqlStXqrCwUMXFxdqyZYsyMjKUm5urgwcPtrr8V199pX79+umxxx5TUlJShxvcEU6+LpGzbGiPEwISgHM4OethX2RV8PEZANYwXFQvWLBAU6ZMUX5+voYMGaLFixerc+fOWrZsWavLX3TRRZo/f75uuOEGuVyuDjfYjjg7ZQ8EhbM4cUQEEC5CMeu5hQ0AwCqGiurGxkZt3rxZOTk536wgMlI5OTkqLy/3W6MaGhpUW1vr9YPAcmrBE46FdTi+Z6MY6QH4jqx3JicMAZfIrGBi2wPWMVRUHz58WE1NTUpMTPR6PDExUVVVVX5rVElJieLj4z0/aWlpflu3k1EY4HQEJAB/I+vN41pM35Bdgeekbe6UDiKzGDUTmmw5+/fMmTNVU1Pj+dm7d2+wmwQHcVJwdIST36dTR0IA8B+yPrw5OcMQuugYg1nRRhZOSEhQVFSUqqurvR6vrq7268QkLpfLttdkhRPXnhg19GkMdjNM2b2vp9LPORTsZiCMMLcCQgVZ71yuyhg19HZmbsM6dGAA1jN0pjomJkbDhw9XaWmp5zG3263S0lJlZ2f7vXFAR4RyiITye/M3p1824eS7FsCZyHoECllmPadt41Af+o3QZehMtSQVFhZq4sSJGjFihEaOHKmFCxeqvr5e+fn5kqS8vDylpqaqpKRE0okJTz7++GPP/+/fv18VFRXq0qWLBgwY4Me3Eh66VjarrndEsJuBIHJaQJ6Ood+A/ZH1CBRGllnH6ccLgJMYLqrHjx+vQ4cOqaioSFVVVcrMzNS6des8E5pUVlYqMvKbE+Cff/65LrzwQs/vTz75pJ588klddtllKisr6/g7sIn47cdUMyA22M3wOycPAZdCL6wJSACBQNab13VPs+r60PkNOA3XU6MjDBfVklRQUKCCgoJW/3Z6eKanp6u5mS8pgifUCmsACASy3pmceF01Oe1/TuyEZ+g3nMxUUe1kMdv2q3FQarCbgQBzYriEokAP/Xb69dQA/Ct+R4Nq+jM5mh1RWPsPxzxA4NnyllpoX6ALBa6BBQAAVqMYBOBUFNUA0EHcTgvAqYJ9bSbDaMMXHRP2Fr+jIdhNgEUoquETzlajo/gOAQDOhKLQPCdvu2B3BAW7IwzOR1ENwHLBKKi5nhoAnMnJxSGA8ERR7UeBHAIajIKBM40AADhDsM/8IbDoiACCi6IahlBYwyi+MwDCEcNJO4Yi0XdO31Z0ACEUOOKWWifvfXnc7Z/7Lh4/bt0Z5aYA3hqy6VhE4F7sFO7/Ouv+lwiemMoYNSnwk3h13despgC+nlX7lJjPDui4n9Z1cv8Z6vcSXrRokebPn6+qqiplZGTomWee0ciRI8/4vBUrVmjChAn6/ve/r9WrV1vfULTgyfrj1k7k09QYuH8Dwcrpk5ye1zs/66reKYeD3Qxbq/w8QQpCzvpT0zF3UF+/697AHTNYvX87uX6y/hsfffSRioqKtHnzZu3Zs0e/+MUvdNddd3ktM2fOHM2dO9frsUGDBmnr1q0+t8kRRXVdXZ0kqezwy/5Z4UH/rKZV/7Bw3QDQQXV1dYqPjw92MyyxcuVKFRYWavHixcrKytLChQuVm5urbdu2qVevXm0+b/fu3brnnns0evToALYWpzuZ9e/883FrX2ijtauHf+0NdgMAByLrv/HVV1+pX79++tGPfqS77767zfWef/75evPNNz2/R0cbK5MdUVSnpKRo79696tq1qyIigtvr6y+1tbVKS0vT3r17FRcXF+zmOALbzDi2mTmhuN2am5tVV1enlJSUYDfFMgsWLNCUKVOUn58vSVq8eLHWrFmjZcuW6f7772/1OU1NTfrJT36iuXPn6u2339aRI0cC2GKciqyHxDYzi+1mXChuM7K+ZdZfdNFFuuiiiySpzWMB6UQRnZSUZLpdjiiqIyMjdc455wS7GZaIi4sLmX/IgcI2M45tZk6obbdQ7bWWpMbGRm3evFkzZ870PBYZGamcnByVl5e3+byHHnpIvXr10qRJk/T2228HoqloA1mPU7HNzGG7GRdq24ysN+ezzz5TSkqKYmNjlZ2drZKSEvXu3dvn5zuiqAYAhKfa2lqv310ul1wuV4vlDh8+rKamJiUmJno9npiY2OY1URs3btTSpUtVUVHht/YCAABjrMx6X2RlZWn58uUaNGiQDhw4oLlz52r06NH68MMP1bVrV5/WQVENAOiQmM8OKDrSv7O3Rv6/idXS0tK8Hi8uLtacOXM6vP66ujrdfPPNWrJkiRISEjq8PgAAQlnMp587Lut9ddVVV3n+f9iwYcrKylKfPn30+9//XpMmTfJpHRTVQeJyuVRcXNxqLwxaxzYzjm1mDtvNPk6/1q2tzyQhIUFRUVGqrq72ery6urrVa6R27Nih3bt3a9y4cZ7H3O4TM9BGR0dr27Zt6t+/vz/eAsIY+xLj2GbmsN2MY5vZh1VZb1a3bt00cOBAbd++3efnRDSH+pzrAABL1NbWKj4+Xjm9Jvu99/q4u1FvHvy1ampqfL7WLSsrSyNHjtQzzzwj6USR3Lt3bxUUFLSYnOTYsWMtwnLWrFmqq6vT008/rYEDByomhnunAgDCm5Oz/nTp6em66667WtxS63RHjx5V7969NWfOHN15550+tYsz1QCAkFBYWKiJEydqxIgRGjlypBYuXKj6+nrPDKF5eXlKTU1VSUmJYmNjNXToUK/nd+vWTZJaPA4AAOzBSNZLJyY3+/jjjz3/v3//flVUVKhLly4aMGCAJOmee+7RuHHj1KdPH33++ecqLi5WVFSUJkyY4HO7KKoBACFh/PjxOnTokIqKilRVVaXMzEytW7fOM6FJZWWlIiMjg9xKAABgltGs//zzz3XhhRd6fn/yySf15JNP6rLLLlNZWZkkad++fZowYYK++OIL9ezZU5dccon++c9/qmfPnj63i+HfAABT7DYkDAAA+BdZ7xu67C20aNEipaenKzY2VllZWdq0aVObyy5ZskSjR49W9+7d1b17d+Xk5LS7fKgyss1OtWLFCkVEROi6666ztoE2ZHSbHTlyRNOnT1dycrJcLpcGDhyotWvXBqi19mF0uy1cuFCDBg1Sp06dlJaWprvvvlvHjh0LUGsB2BVZbxxZbw55bxxZj0ChqLbIypUrVVhYqOLiYm3ZskUZGRnKzc3VwYMHW12+rKxMEyZM0N/+9jeVl5crLS1NV1xxhfbv3x/glgeP0W120u7du3XPPfdo9OjRAWqpfRjdZo2NjRo7dqx2796tVatWadu2bVqyZIlSU1MD3PLgMrrdXn31Vd1///0qLi7WJ598oqVLl2rlypV64IEHAtxyAHZC1htH1ptD3htH1iOQGP5tkaysLF100UV69tlnJZ2YmS4tLU133HHHGWemk6SmpiZ1795dzz77rPLy8qxuri2Y2WZNTU269NJLdeutt+rtt9/WkSNHtHr16gC2OriMbrPFixdr/vz52rp1q84666xAN9c2jG63goICffLJJyotLfU89rOf/UzvvvuuNm7cGLB22w1DwhDuyHrjyHpzyHvjyHr/IOt9w5lqCzQ2Nmrz5s3KycnxPBYZGamcnByVl5f7tI6vvvpKX3/9tXr06GFVM23F7DZ76KGH1KtXL59vzB5KzGyzN954Q9nZ2Zo+fboSExM1dOhQzZs3T01NTYFqdtCZ2W4XX3yxNm/e7Bk2tnPnTq1du1ZXX311QNoMwH7IeuPIenPIe+PIegQas39b4PDhw2pqavLMQndSYmKitm7d6tM67rvvPqWkpHjtDEKZmW22ceNGLV26VBUVFQFoof2Y2WY7d+7UX//6V/3kJz/R2rVrtX37dt1+++36+uuvVVxcHIhmB52Z7XbjjTfq8OHDuuSSS9Tc3Kzjx49r6tSpDAkDwhhZbxxZbw55bxxZj0DjTLUNPfbYY1qxYoX++Mc/KjY2NtjNsaW6ujrdfPPNWrJkiRISEoLdHMdwu93q1auXXnjhBQ0fPlzjx4/Xgw8+qMWLFwe7abZWVlamefPm6bnnntOWLVv0+uuva82aNXr44YeD3TQADkXWnxlZbx55bxxZj47gTLUFEhISFBUVperqaq/Hq6urlZSU1O5zn3zyST322GN68803NWzYMCubaStGt9mOHTu0e/dujRs3zvOY2+2WJEVHR2vbtm3q37+/tY0OMjPfs+TkZJ111lmKioryPHbeeeepqqpKjY2Nionx77UydmRmu82ePVs333yzJk+eLEm64IILVF9fr5/+9Kd68MEHufcxEIbIeuPIenPIe+PIegQa3w4LxMTEaPjw4V4THbjdbpWWlio7O7vN5z3xxBN6+OGHtW7dOo0YMSIQTbUNo9ts8ODB+uCDD1RRUeH5ufbaazVmzBhVVFQoLS0tkM0PCjPfs1GjRmn79u2egxJJ+vTTT5WcnBzyAXuSme321VdftQjTkwcqzPUIhCey3jiy3hzy3jiyHoHGmWqLFBYWauLEiRoxYoRGjhyphQsXqr6+Xvn5+ZKkvLw8paamqqSkRJL0+OOPq6ioSK+++qrS09NVVVUlSerSpYu6dOkStPcRSEa2WWxsrIYOHer1/G7duklSi8dDmdHv2bRp0/Tss89qxowZuuOOO/TZZ59p3rx5uvPOO4P5NgLO6HYbN26cFixYoAsvvFBZWVnavn27Zs+erXHjxnmdBQAQXsh648h6c8h748h6BBJFtUXGjx+vQ4cOqaioSFVVVcrMzNS6des8EyZUVlZ69YY9//zzamxs1PXXX++1nuLiYs2ZMyeQTQ8ao9sMxrdZWlqa1q9fr7vvvlvDhg1TamqqZsyYofvuuy9YbyEojG63WbNmKSIiQrNmzdL+/fvVs2dPjRs3To8++miw3gIAGyDrjSPrzSHvjSPrEUjcpxoAYAr3rgQAILSR9b6hKxAAAAAAAJMoqgEAAAAAMImiGgAAAAAAkyiqAQAAAAAwiaIaAAAAAACTKKoBAAAAADCJohoAAAAAAJMoqgEAAAAAMImiGgAAAAAAkyiqAQAAAAAwiaIaAAAAAACTKKoBAAAAADCJohoAAAAAAJMoqgEAAAAAMImiGgAAAAAAkyiqAQAAAAAwiaIaAAAAAACTKKoBAAAAADCJohoAAAAAAJMoqgEAAAAAMImiGgAAAAAAkyiqAQAAAAAwiaIaAAAAAACTKKoBAAAAADCJohoAAAAAAJMoqgEAAAAAMImiGgAAAAAAkyiqAQAAAAAwiaIaAAAAAACTKKoBAAAAADCJohoAAAAAAJMoqgEAAAAAMImiGgAAAAAAkyiqAQAAAAAwiaIaAAAAAACTKKoBAAAAADCJohoAEDIWLVqk9PR0xcbGKisrS5s2bWpz2SVLlmj06NHq3r27unfvrpycnHaXBwAAwWck6yXptdde0+DBgxUbG6sLLrhAa9eu9fp7c3OzioqKlJycrE6dOiknJ0efffaZoTZRVAMAQsLKlStVWFio4uJibdmyRRkZGcrNzdXBgwdbXb6srEwTJkzQ3/72N5WXlystLU1XXHGF9u/fH+CWAwAAXxjN+nfeeUcTJkzQpEmT9O9//1vXXXedrrvuOn344YeeZZ544gn98pe/1OLFi/Xuu+/q7LPPVm5uro4dO+ZzuyKam5ubO/zuAABhp7a2VvHx8crpNVnRkTF+Xfdxd6PePPhr1dTUKC4uzqfnZGVl6aKLLtKzzz4rSXK73UpLS9Mdd9yh+++//4zPb2pqUvfu3fXss88qLy+vQ+0HACAUOD3rx48fr/r6ev3f//t/PY/9n//zf5SZmanFixerublZKSkp+tnPfqZ77rlHklRTU6PExEQtX75cN9xwg0/t4kw1AMDxGhsbtXnzZuXk5Hgei4yMVE5OjsrLy31ax1dffaWvv/5aPXr0sKqZAADAJDNZX15e7rW8JOXm5nqW37Vrl6qqqryWiY+PV1ZWls/HD5IUbeSNAABwuuPNjZLbgnXqRA/5qVwul1wuV4vlDx8+rKamJiUmJno9npiYqK1bt/r0mvfdd59SUlJahC8AAOHOqVlfVVXV6vJVVVWev598rK1lfEFRDQAwJSYmRklJSSqretmS9Xfp0kVpaWlejxUXF2vOnDl+f63HHntMK1asUFlZmWJjY/2+fgAAnCiUst5KFNUAAFNiY2O1a9cuNTY2WrL+5uZmRUREeD3WWs+1JCUkJCgqKkrV1dVej1dXVyspKand13nyySf12GOP6c0339SwYcM61mgAAEKI07M+KSmp3eVP/re6ulrJycley2RmZvr8PiiqAQCmxcbG2uLMbkxMjIYPH67S0lJdd911kk5MXlJaWqqCgoI2n/fEE0/o0Ucf1fr16zVixIgAtRYAAOdwctZnZ2ertLRUd911l+exDRs2KDs7W5LUt29fJSUlqbS01FNE19bW6t1339W0adN8bhtFNQAgJBQWFmrixIkaMWKERo4cqYULF6q+vl75+fmSpLy8PKWmpqqkpESS9Pjjj6uoqEivvvqq0tPTPddOdenSRV26dAna+wAAAK0zmvUzZszQZZddpqeeekrXXHONVqxYoffee08vvPCCJCkiIkJ33XWXHnnkEZ177rnq27evZs+erZSUFE/h7guKagBASBg/frwOHTqkoqIiVVVVKTMzU+vWrfNMPlJZWanIyG9uevH888+rsbFR119/vdd6nHgtFwAA4cBo1l988cV69dVXNWvWLD3wwAM699xztXr1ag0dOtSzzL333qv6+nr99Kc/1ZEjR3TJJZdo3bp1hs7Oc59qAAAAAABM4j7VAAAAAACYRFENAAAAAIBJFNUAAAAAAJhEUQ0AAAAAgEkU1QAAAAAAmERRDQAAAACASRTVAAAAAACYRFENAAAAAIBJFNUAAAAAAJhEUQ0AAAAAgEkU1QAAAAAAmERRDQAAAACASf8/Wry9UhQzZBkAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "xmesh, ymesh = np.meshgrid(xgrid, ygrid)\n", + "\n", + "fig, axs = plt.subplots(1, 2, figsize=(12, 6))\n", + "contour0 = axs[0].contourf(xmesh, ymesh, result_matrix.transpose())\n", + "axs[0].axis(\"equal\")\n", + "axs[0].axis(\"square\")\n", + "axs[0].title.set_text(\"Quantum\")\n", + "fig.colorbar(contour0, ax=axs[0])\n", + "\n", + "contour1 = axs[1].contourf(xmesh, ymesh, classical_result_matrix.transpose())\n", + "axs[1].axis(\"equal\")\n", + "axs[1].axis(\"square\")\n", + "axs[1].title.set_text(\"Classical\")\n", + "fig.colorbar(contour1, ax=axs[1])\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "42", + "metadata": { + "collapsed": false, + "jupyter": { + "outputs_hidden": false + } + }, + "source": [ + "## References\n", + "\n", + "[1]: [Strang, Gilbert, 1999 SIAM Review 41 135. The Discrete Cosine Transform.](https://doi.org/10.1137/S0036144598336745)\n", + "\n", + "[2]: [Yudong Cao et al 2013 New J. Phys. 15 013021. Quantum algorithm and circuit design solving the Poisson equation.](https://iopscience.iop.org/article/10.1088/1367-2630/15/1/013021/pdf)\n", + "\n", + "[3]: [Klappenecker, A., & Rotteler M., \"Discrete Cosine Transforms on Quantum Computers\".](https://arxiv.org/abs/quant-ph/0111038)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/algorithms/discrete_poisson_solver/discrete_poisson_solver.qmod b/algorithms/discrete_poisson_solver/discrete_poisson_solver.qmod new file mode 100644 index 00000000..31d94879 --- /dev/null +++ b/algorithms/discrete_poisson_solver/discrete_poisson_solver.qmod @@ -0,0 +1,280 @@ +qfunc qsct_2d(xy_variable: qnum[2]) { + qst_type2(xy_variable[0]); + qct_type2(xy_variable[1]); +} + +qfunc powered_hamiltonian_evolution(qba: qbit[]) { + suzuki_trotter(qba); +} + +qfunc inverse_amplitude_load(phase: qnum, ind: qbit) { + ind *= prefactor / phase; +} + +qfunc matrix_inversion_HHL(arg1: qbit[])>(state: qbit[], phase: qnum, output indicator: qbit) { + allocate<1>(indicator); + within { + qpe_flexible() { + my_unitary(state); + }>(phase); + } apply { + inverse_amplitude_load(phase, indicator); + } +} + +qfunc main(output x_variable: qnum<3, False, 0>, output y_variable: qnum<3, False, 0>, output phase: qnum, output indicator: qbit) { + xy_variable: qnum<3, False, 0>[2]; + prepare_amplitudes<[ + 0.000929993206, + 0.006286469531, + 0.012502019419, + 0.01640293531, + 0.01640293531, + 0.012502019419, + 0.006286469531, + 0.000929993206, + 0.006286469531, + 0.042483535101, + 0.084462252894, + 0.110795378713, + 0.110795378713, + 0.084462252894, + 0.042483535101, + 0.006286469531, + 0.012502019419, + 0.084462252894, + 0.167861987331, + 0.220148535874, + 0.220148535874, + 0.167861987331, + 0.084462252894, + 0.012502019419, + 0.01640293531, + 0.110795378713, + 0.220148535874, + 0.288681749865, + 0.288681749865, + 0.220148535874, + 0.110795378713, + 0.01640293531, + 0.01640293531, + 0.110795378713, + 0.220148535874, + 0.288681749865, + 0.288681749865, + 0.220148535874, + 0.110795378713, + 0.01640293531, + 0.012502019419, + 0.084462252894, + 0.167861987331, + 0.220148535874, + 0.220148535874, + 0.167861987331, + 0.084462252894, + 0.012502019419, + 0.006286469531, + 0.042483535101, + 0.084462252894, + 0.110795378713, + 0.110795378713, + 0.084462252894, + 0.042483535101, + 0.006286469531, + 0.000929993206, + 0.006286469531, + 0.012502019419, + 0.01640293531, + 0.01640293531, + 0.012502019419, + 0.006286469531, + 0.000929993206 + ], 0.0>(xy_variable); + allocate<6>(phase); + within { + qsct_2d(xy_variable); + } apply { + matrix_inversion_HHL<0.015625, lambda

(target) { + powered_hamiltonian_evolution<[ + PauliTerm { + pauli=[ + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I + ], + coefficient=2.25 + }, + PauliTerm { + pauli=[ + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::Z + ], + coefficient=-0.25 + }, + PauliTerm { + pauli=[ + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::Z, + Pauli::I + ], + coefficient=-0.520598050073 + }, + PauliTerm { + pauli=[ + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::Z, + Pauli::Z + ], + coefficient=-0.020598050073 + }, + PauliTerm { + pauli=[ + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::Z, + Pauli::I, + Pauli::I + ], + coefficient=-1.256834873031 + }, + PauliTerm { + pauli=[ + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::Z, + Pauli::I, + Pauli::Z + ], + coefficient=-0.049728091845 + }, + PauliTerm { + pauli=[ + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::Z, + Pauli::Z, + Pauli::I + ], + coefficient=-0.103553390593 + }, + PauliTerm { + pauli=[ + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::Z, + Pauli::Z, + Pauli::Z + ], + coefficient=0.103553390593 + }, + PauliTerm { + pauli=[ + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I + ], + coefficient=1.75 + }, + PauliTerm { + pauli=[ + Pauli::I, + Pauli::I, + Pauli::Z, + Pauli::I, + Pauli::I, + Pauli::I + ], + coefficient=-0.25 + }, + PauliTerm { + pauli=[ + Pauli::I, + Pauli::Z, + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I + ], + coefficient=-0.520598050073 + }, + PauliTerm { + pauli=[ + Pauli::I, + Pauli::Z, + Pauli::Z, + Pauli::I, + Pauli::I, + Pauli::I + ], + coefficient=0.020598050073 + }, + PauliTerm { + pauli=[ + Pauli::Z, + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I + ], + coefficient=-1.256834873031 + }, + PauliTerm { + pauli=[ + Pauli::Z, + Pauli::I, + Pauli::Z, + Pauli::I, + Pauli::I, + Pauli::I + ], + coefficient=0.049728091845 + }, + PauliTerm { + pauli=[ + Pauli::Z, + Pauli::Z, + Pauli::I, + Pauli::I, + Pauli::I, + Pauli::I + ], + coefficient=0.103553390593 + }, + PauliTerm { + pauli=[ + Pauli::Z, + Pauli::Z, + Pauli::Z, + Pauli::I, + Pauli::I, + Pauli::I + ], + coefficient=0.103553390593 + } + ], 0.101321183642, p>(target); + }>(xy_variable, phase, indicator); + } + xy_variable -> {x_variable, y_variable}; +} + diff --git a/algorithms/discrete_poisson_solver/discrete_poisson_solver.synthesis_options.json b/algorithms/discrete_poisson_solver/discrete_poisson_solver.synthesis_options.json new file mode 100644 index 00000000..f28ff899 --- /dev/null +++ b/algorithms/discrete_poisson_solver/discrete_poisson_solver.synthesis_options.json @@ -0,0 +1,5 @@ +{ + "constraints": { + "max_width": 18 + } +} \ No newline at end of file