diff --git a/docs/algorithms/arnoldi.rst b/docs/algorithms/arnoldi.rst index cf35a10..ec6f429 100644 --- a/docs/algorithms/arnoldi.rst +++ b/docs/algorithms/arnoldi.rst @@ -38,7 +38,7 @@ based on previous estimates, using the formula :math:`|{\xi_{k+1}}\rangle=(1-\ga memory factor :math:`\gamma=-0.75`. -The rule in :ref:`alg_descent` is a particular case of the Arnoldi iteration with a Krylov basis with $L=2$. +The rule in :ref:`alg_descent` is a particular case of the Arnoldi iteration with a Krylov basis with :math:`L=2`. Ref. :cite:t:`GarciaMolina2024` presents this algorithm and its implementation for global optimization problems. It is also suitable for evolution problems. diff --git a/docs/algorithms/generated/seemps.evolution.euler.euler.rst b/docs/algorithms/generated/seemps.evolution.euler.euler.rst new file mode 100644 index 0000000..6c18a8c --- /dev/null +++ b/docs/algorithms/generated/seemps.evolution.euler.euler.rst @@ -0,0 +1,11 @@ + + +seemps.evolution.euler.euler +============================ + +.. currentmodule:: seemps.evolution.euler + + + +.. autofunction:: seemps.evolution.euler.euler + diff --git a/docs/algorithms/generated/seemps.evolution.euler.euler2.rst b/docs/algorithms/generated/seemps.evolution.euler.euler2.rst new file mode 100644 index 0000000..4a583b4 --- /dev/null +++ b/docs/algorithms/generated/seemps.evolution.euler.euler2.rst @@ -0,0 +1,11 @@ + + +seemps.evolution.euler.euler2 +============================= + +.. currentmodule:: seemps.evolution.euler + + + +.. autofunction:: seemps.evolution.euler.euler2 + diff --git a/docs/algorithms/generated/seemps.evolution.runge_kutta.runge_kutta.rst b/docs/algorithms/generated/seemps.evolution.runge_kutta.runge_kutta.rst new file mode 100644 index 0000000..34639ee --- /dev/null +++ b/docs/algorithms/generated/seemps.evolution.runge_kutta.runge_kutta.rst @@ -0,0 +1,11 @@ + + +seemps.evolution.runge\_kutta.runge\_kutta +========================================== + +.. currentmodule:: seemps.evolution.runge_kutta + + + +.. autofunction:: seemps.evolution.runge_kutta.runge_kutta + diff --git a/docs/algorithms/generated/seemps.evolution.runge_kutta.runge_kutta_fehlberg.rst b/docs/algorithms/generated/seemps.evolution.runge_kutta.runge_kutta_fehlberg.rst new file mode 100644 index 0000000..d52296c --- /dev/null +++ b/docs/algorithms/generated/seemps.evolution.runge_kutta.runge_kutta_fehlberg.rst @@ -0,0 +1,11 @@ + + +seemps.evolution.runge\_kutta.runge\_kutta\_fehlberg +==================================================== + +.. currentmodule:: seemps.evolution.runge_kutta + + + +.. autofunction:: seemps.evolution.runge_kutta.runge_kutta_fehlberg + diff --git a/docs/algorithms/index.rst b/docs/algorithms/index.rst index e21f1a4..b468c35 100644 --- a/docs/algorithms/index.rst +++ b/docs/algorithms/index.rst @@ -9,6 +9,9 @@ Index of algorithms tensor_split mps_simplification gradient_descent + cgs arnoldi dmrg + runge_kutta + crank_nicolson tebd_evolution diff --git a/docs/algorithms/runge_kutta.rst b/docs/algorithms/runge_kutta.rst new file mode 100644 index 0000000..bca568a --- /dev/null +++ b/docs/algorithms/runge_kutta.rst @@ -0,0 +1,70 @@ +.. _mps_runge_kutta: + +******************* +Runge-Kutta methods +******************* + +Runge-Kutta methods use a Taylor expansion of the state to approximate + +.. math:: + \psi_{k+1} = \psi_k + \sum_{p}\frac{1}{p!}(-\Delta \beta H)^p\frac{\partial^{(p)} \psi}{\partial \beta^{(p)}}. + +:math:`\Delta \beta` can either be real or imaginary time, and hence these +methods are suitable for both evolution and optimization problems. When using +MPS, these methods perform a global optimization of the MPS at each step, +i.e., they update all tensors simultaneously. + +The order of the expansion `p` determines the truncation error of the method, which +is :math:`O(\Delta \beta ^{p+1})`, and also the cost of the method, since +a higher order implies more operations. Thus, it is important to consider +trade-off in cost and accuracy to choose the most suitable method for each application. + +The SeeMPS library considers four methods. + +1. Euler method +---------------- + +This is an explicit, first-order Taylor approximation of the evolution, with an error :math:`\mathcal{O}(\Delta\beta^2)` +and a simple update with a fixed time-step :math:`\beta_k = k \Delta\beta`. + +.. math:: + \psi_0 &= \psi(\beta_0), \\ + \psi_{k+1} &= \psi_k - \Delta\beta H \psi_k, \quad \text{for } k=0,1,\dots,N-1. + +2. Improved Euler or Heun method +--------------------------------- + +This is a second-order, fixed-step explicit method that uses two matrix-vector multiplications and two linear combinations of +vectors to achieve an error :math:`\mathcal{O}(\Delta\beta^3)`. + +.. math:: + \psi_{k+1} = \psi_k - \frac{\Delta\beta}{2} \left[v_1 + H(\psi_k - \Delta\beta v_1)\right], \\ + \text{with } v_1 = H \psi_k. + +3. Fourth-order Runge-Kutta method +----------------------------------- + +This algorithm achieves an error :math:`\mathcal{O}(\Delta\beta^5)` using four matrix-vector multiplications and four linear combinations of vectors. + +.. math:: + \psi_{k+1} &= \psi_k + \frac{\Delta\beta}{6}(v_1 + 2v_2 + 2v_3 + v_4), \\ + v_1 &= -H \psi_k, \\ + v_2 &= -H\left(\psi_k + \frac{\Delta\beta}{2}v_1\right), \\ + v_3 &= -H\left(\psi_k + \frac{\Delta\beta}{2}v_2\right), \\ + v_4 &= -H\left(\psi_k + \Delta\beta v_3\right). + +4. Runge-Kutta-Fehlberg method +------------------------------- +The Runge-Kutta-Fehlberg algorithm is an adaptative step-size solver that combines a fifth-order accurate integrator +:math:`O(\Delta\beta^5)` with a sixth-order error estimator :math:`O(\Delta\beta^6)`. This combination dynamically adjusts the step size :math:`\Delta\beta` +to maintain the integration error within a specified tolerance. The method requires an initial estimate of the step size, which can be obtained from a simpler +method. Each iteration involves six matrix-vector multiplications and six linear combinations, and it may repeat the evolution steps if the proposed step size +is deemed unsuitable. + +.. autosummary:: + :toctree: generated/ + + ~seemps.evolution.euler.euler + ~seemps.evolution.euler.euler2 + ~seemps.evolution.runge_kutta.runge_kutta + ~seemps.evolution.runge_kutta.runge_kutta_fehlberg diff --git a/docs/seemps_analysis_differentiation.rst b/docs/seemps_analysis_differentiation.rst index 1a6ba7d..c02652e 100644 --- a/docs/seemps_analysis_differentiation.rst +++ b/docs/seemps_analysis_differentiation.rst @@ -40,6 +40,8 @@ leading to This approximation is improved with a smoother formula to avoid noise resilience following `Holoborodko `_. +An example on how to use these functions is shown in `Differentiation.ipynb `_. + .. autosummary:: :toctree: generated/ diff --git a/docs/seemps_examples.rst b/docs/seemps_examples.rst index a8f4574..c324d28 100644 --- a/docs/seemps_examples.rst +++ b/docs/seemps_examples.rst @@ -11,6 +11,8 @@ In the directory `examples` of the library, we have created various Jupyter note - `Computation of a spin-1/2 model ground state using DMRG `_. - `Solution of a Hamiltonian PDE `_. - `Function interpolation `_. +- `Function differentiation `_. + Further examples of use appear in the datasets and associated codes to the following publications: diff --git a/examples/Differentiation.ipynb b/examples/Differentiation.ipynb new file mode 100644 index 0000000..f230a94 --- /dev/null +++ b/examples/Differentiation.ipynb @@ -0,0 +1,192 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Function differentiation\n", + "\n", + "This notebook shows how to differentiate functions represented as MPS using the SeeMPS library. We focus on the finite difference method, since the QFT is not yet optimally implemented to efficiently approximate derivatives. More information about these methods can be found in https://arxiv.org/abs/2303.09430." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, let us import the necessary set of tools." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from seemps.state import MPS\n", + "from seemps.analysis.finite_differences import smooth_finite_differences_mpo" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Gaussian function\n", + "\\begin{equation}\n", + " f(x) = \\frac{1}{\\mathcal{N}}e^{-x^2/2},\n", + "\\end{equation}\n", + "where $\\mathcal{N}$ is normalization constant, acts as benchmark for the differentiation techniques." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us define the functions to compute the error." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "def remove_phase(v):\n", + " \"\"\"\n", + " Removes the phase from a complex vector `v` by normalizing it with the phase \n", + " of its maximum magnitude element.\n", + " \n", + " Parameters\n", + " ----------\n", + " v : numpy.ndarray\n", + " A complex vector or array-like object.\n", + "\n", + " Returns\n", + " -------\n", + " numpy.ndarray\n", + " A complex vector of the same shape as `v` with its phase removed.\n", + " \n", + " Notes\n", + " -----\n", + " The phase of the vector is adjusted by multiplying it by the conjugate of \n", + " the element with the maximum magnitude, normalized by the absolute value of \n", + " that element. This operation effectively rotates the vector to align its \n", + " largest element with the positive real axis.\n", + " \"\"\"\n", + " shape = v.shape\n", + " v = v.flatten()\n", + " k = np.argmax(np.abs(v))\n", + " return (v * np.abs(v[k]) / v[k]).reshape(shape)\n", + "\n", + "def norm2_difference(a, b):\n", + " \"\"\"\n", + " Calculates the norm-2 difference between two vectors or MPS objects after \n", + " removing their phases and normalizing them.\n", + "\n", + " Parameters\n", + " ----------\n", + " a : numpy.ndarray or MPS\n", + " The first vector or MPS object.\n", + " b : numpy.ndarray or MPS\n", + " The second vector or MPS object.\n", + "\n", + " Returns\n", + " -------\n", + " float\n", + " The norm-2 (Euclidean) difference between the two vectors after phase \n", + " removal and normalization.\n", + "\n", + " Notes\n", + " -----\n", + " If `a` or `b` are instances of the `MPS` class, they are first converted to \n", + " vector form using the `to_vector()` method. The vectors are then phase-aligned \n", + " using the `remove_phase` function and normalized to unit length. The difference \n", + " is calculated using the Euclidean norm of the difference vector.\n", + " \"\"\"\n", + " if isinstance(a, MPS):\n", + " a = a.to_vector()\n", + " if isinstance(b, MPS):\n", + " b = b.to_vector()\n", + " a = remove_phase(a)\n", + " b = remove_phase(b)\n", + " a /= np.linalg.norm(a)\n", + " b /= np.linalg.norm(b)\n", + " return np.linalg.norm(a - b)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us study the error scaling with the number of qubits $n$, which is equivalent to studying how it varies with $\\Delta x$. The finite difference function has a default minimum value of $\\Delta x$ to avoid round-off error. This value can be modified according to the problem." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "a, b = -10, 10\n", + "n_qubits = range(6,23)\n", + "f = lambda x: np.exp(- x ** 2)\n", + "df = lambda x: - 2 * x * f(x)\n", + "errors = np.empty(shape=(len(n_qubits)))\n", + "\n", + "for i, n in enumerate(n_qubits):\n", + " x, dx = np.linspace(a, b, num=2 ** n, endpoint=False, retstep=True)\n", + " f_mps = MPS.from_vector(f(x), [2]*n, normalize=False)\n", + " K = smooth_finite_differences_mpo(n, order=1, filter=9, dx=dx, periodic=True, tol=1e-3)\n", + " errors[i] = norm2_difference((K @ f_mps), df(x))" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABCAElEQVR4nO3deXxU1cH/8e9kDyEJaxIiCSAKGJZEIURcQaOIffBBasW2WlxqHzWgmJaK+lNpreJKqTKPPLbWpWpLbQtuFdCoLIKsRsQgi4Y1JCGGZLJOkpn7+wOZEiGEhEzOLJ/365U/cnOY+V6D4Ztzz7nXZlmWJQAAAPi9ENMBAAAA0DEodgAAAAGCYgcAABAgKHYAAAABgmIHAAAQICh2AAAAAYJiBwAAECAodgAAAAEizHQAX+d2u1VUVKTY2FjZbDbTcQAAQJCxLEtVVVVKTk5WSMiJ5+Qodq0oKipSSkqK6RgAACDI7d27V3379j3hGIpdK2JjYyUd/o8ZFxdnOA0AAAg2DodDKSkpnk5yIhS7Vhy5/BoXF0exAwAAxpzMkjA2T7TAbrcrLS1NmZmZpqMAAACcFJtlWZbpEL7M4XAoPj5elZWVzNgBAIBO15YuwowdAABAgKDYAQAABAiKHQAAQICg2AEAAAQIih0AAECAoNgBAAAECIodAABAgKDYAQAABAiKHQAAQICg2AEAAAQIih0AAECAoNgBAAAECIpdC+x2u9LS0pSZmdkp79fkcnfK+wAAgMBlsyzLMh3ClzkcDsXHx6uyslJxcXEd/vqWZenVT3frlTW79c87zlNcVHiHvwcAAPBfbekizNgZVtfo0vMrv9GO0mr9de0e03EAAIAfCzMdINh1iQjT76/N0OZ9lbrxvP6m4wAAAD9GsfMBo/r30Kj+PUzHAAAAfo5LsT7G5ba06LN9crtZ+ggAANqGGTsfYlmWbnxxnVbuKNOhmkbdfMEA05EAAIAfYcbOh9hsNl2elqiukWHqHsPuWAAA0DbM2PmY68/tp8uHJikxLsp0FAAA4GeYsfMxNputWaljrR0AADhZFDsftr2kSlfZV+nTb741HQUAAPgBip0Pe3n1Lm3Z79Aj724VDwgBAACtYY2dD7vvyrPktqS7LztTNpvNdBwAAODjKHY+LCYyTHMmDzcdAwAA+AkuxfqRz/YcUlFFnekYAADAR1Hs/MQbG/bqmgVr9Ks3PmenLAAAOC6KnZ8Y2a+7IkJD1KtrpJxNbtNxAACAD2KNnZ84vXdXLZ1xkVJ7djEdBQAA+KigmLG7+uqr1b17d11zzTWmo5yS75c6boECAACOFhTF7q677tIrr7xiOkaHqW906eF3CvTYe1+ZjgIAAHxIUBS7sWPHKjY21nSMDrN+V7leWFWo51d+o68PVpuOAwAAfITxYrdixQpNnDhRycnJstlsWrx48TFj7Ha7+vfvr6ioKGVlZWndunWdH9SHXHhmb90xdqBemDpKA3t3NR0HAAD4COObJ2pqapSenq6bb75ZkydPPubrCxcuVG5urhYsWKCsrCzNmzdP48eP17Zt25SQkCBJysjIUFNT0zF/dtmyZUpOTvb6OZjw6yuGmI4AAAB8jPFiN2HCBE2YMKHFr8+dO1e33nqrbrrpJknSggUL9O677+rPf/6zZs2aJUnKz8/vsDxOp1NOp9PzucPh6LDX9pZqZ5O2FTs0sl8P01EAAIBBxi/FnkhDQ4M2btyo7Oxsz7GQkBBlZ2drzZo1XnnPOXPmKD4+3vORkpLilffpKHvLazXhDyt044vrtZ+nUgAAENR8utiVlZXJ5XIpMTGx2fHExEQVFxef9OtkZ2frRz/6kf7973+rb9++JyyF9957ryorKz0fe/fubXf+ztAnPkq9ukYqLipc5dUNpuMAAACDjF+K7QwffPDBSY+NjIxUZGSkF9N0rLDQENl/co66RoUpLircdBwAAGCQTxe7Xr16KTQ0VCUlJc2Ol5SUKCkpyVAq35PcLdp0BAAA4AN8+lJsRESERo4cqby8PM8xt9utvLw8jRkzxqvvbbfblZaWpszMTK++T0f7aFuprv/TWtU3ukxHAQAAncx4sauurlZ+fr5nZ2thYaHy8/O1Z88eSVJubq7++Mc/6uWXX9bWrVt1++23q6amxrNL1ltycnJUUFCg9evXe/V9OlJtQ5N+/Y/NWrWzTC+sKjQdBwAAdDLjl2I3bNigcePGeT7Pzc2VJE2dOlUvvfSSpkyZooMHD+rBBx9UcXGxMjIytGTJkmM2VEDqEhGmJ344Qsu3H9QtFwwwHQcAAHQym8WT5E/I4XAoPj5elZWViouLMx0HAAAEmbZ0EeOXYn2Vv66x+74Nu8pFdwcAIDgwY9cKf52xsyxLv/7HZr2xcZ/mXpuuyef0NR0JAAC0AzN2kM1mU7+eXRRik0ocztb/AAAA8HvGN0/Ae24fe4bGDk7QsNPiTUcBAACdgBm7FgTCGrvQEBulDgCAIMIau1b46xq77ztY5dQDi7do2iVnUPYAAPAjrLHDMZ5c+pWWfFmse/65mV2yAAAEKNbYBYlZE87Sgcp63f+Ds2Sz2UzHAQAAXkCxCxI9YiL0l1uyTMcAAABexKXYILXvUK0O1TSYjgEAADoQxa4FgbArtiXvF5RowryVum/RF6y3AwAggFDsWpCTk6OCggKtX7/edJQO1yc+SnWNLpU46lXT4DIdBwAAdBDW2AWhYafF62+/OFcZKd0UFkq3BwAgUFDsgtSo/j1MRwAAAB2M6Zog53Zb+suaXXpt7W7TUQAAwClixi7ILf2yWA+8+aWiw0M1dnCCTusWbToSAABoJ4pdC+x2u+x2u1yuwN5cMH5oki4dkqCLBvVWn7go03EAAMAp4FmxrQiUZ8WeiGVZPI0CAAAfxbNi0SZHlzqX21JFLTcuBgDAH1Hs4LGrrEY/WrBad7y2iRsXAwDgh1hjh2YKDjgUFhKirw9W64yEWNNxAABAG1Ds4NG/V4z+cN3ZGnZaPLtjAQDwQxQ7NDN+aJLpCAAAoJ1YY4cWbSuu0vsFJaZjAACAk0Sxa4HdbldaWpoyMzNNRzFi055D+q9nVyp3Yb6KKupMxwEAACeBYteCnJwcFRQUaP369aajGJHet5uGnRavrNN7KiyUe9wBAOAPWGOH4woNsenlm0crNjKMmxcDAOAnmLFDi+KiwpuVOu5tBwCAb6PYoVXOJpeeXrZN//OXjZQ7AAB8GJdi0aq95XX6v+XfqMHl1qfflGvMwJ6mIwEAgOOg2KFVZyR01f/7r7PUq2skpQ4AAB9GscNJ+dmY/qYjAACAVrDGDm1W3+jSlv2VpmMAAIDvodi1INhvUNySfYdqdeUzK3X9C2t1sMppOg4AADgKxa4FwX6D4pYkxkUpKixUEaEh2nuo1nQcAABwFNbYoU3CQ0P03PXnqFt0hOK7hJuOAwAAjkKxQ5v16xljOgIAADgOLsXilHz6zbd65N0C0zEAAICYscMpKK6s1w0vrFWjy9Ko/j00fmiS6UgAAAQ1ih3aLSk+SrddPFBl1Q06jxsXAwBgHMUOpyT3skGy2WymYwAAALHGDqfo+6WuorbBUBIAAECxQ4eocTZp1j8367Lfr6DcAQBgCMUOHSI0xKYNuw/pYJVTH20rNR0HAICgxBo7dIio8FDNvTZddQ0uZZ3ORgoAAEyg2KHDjOjbzXQEAACCGpdi4RUVtQ3627o9pmMAABBUmLFrgd1ul91ul8vlMh3F71Q7m3TFvJUqdtQrMT5K4wYnmI4EAEBQYMauBTk5OSooKND69etNR/E7XSPD9F8j+mhg7xh17xJhOg4AAEHDZlmWZTqEL3M4HIqPj1dlZaXi4uJMx/Eb9Y2HZzqjwkMNJwEAwL+1pYtwKRZe8f1CZ1kWT6gAAMDLuBQLr7IsS+9sLtJV8z9RVX2j6TgAAAQ0ih28ytnk1hNLtumL/ZV6YVWh6TgAAAQ0LsXCq6LCQ/XENSO0+utvdcfYM0zHAQAgoFHs4HXnnt5T5/I0CgAAvI5LsehUlmXpsz2HTMcAACAgUezQaZpcbv385Q26+n9Xa/XXZabjAAAQcCh26DRhoSFKjI9SRGiIdpXVmo4DAEDAYY0dOtW9E4boxvP6a1BirOkoAAAEHIodOlVsVLhio8JNxwAAICBxKRbG7DtUq1//43PVNbhMRwEAICAwYwcj3G5LN764XjtLq9W9S4TuvfIs05EAAPB7zNjBiJAQm+67cogy+3fXlMwU03EAAAgIzNjBmEuGJGrc4ATZbDbTUQAACAjM2MGoo0tdtbPJYBIAAPxfwBe7vXv3auzYsUpLS9OIESP0xhtvmI6E77EsS39eVajz5uRp6wGH6TgAAPitgC92YWFhmjdvngoKCrRs2TLNmDFDNTU1pmPhe9YWfitHfZP+tm6P6SgAAPitgF9j16dPH/Xp00eSlJSUpF69eqm8vFwxMTGGk+EIm82mhycN09jBCZoyio0UAAC0l/EZuxUrVmjixIlKTk6WzWbT4sWLjxljt9vVv39/RUVFKSsrS+vWrWvXe23cuFEul0spKZQHX5MQG6Ufj05VSAgbKQAAaC/jxa6mpkbp6emy2+3H/frChQuVm5urhx56SJs2bVJ6errGjx+v0tJSz5iMjAwNGzbsmI+ioiLPmPLycv3sZz/T888/7/Vzwqlpcrn1r0375HZbpqMAAOBXbJZl+cy/njabTYsWLdKkSZM8x7KyspSZman58+dLktxut1JSUjR9+nTNmjXrpF7X6XTqsssu06233qobbrih1bFOp9PzucPhUEpKiiorKxUXF9f2k0KbWJal619Yq092fqvfXDVUU8/rbzoSAABGORwOxcfHn1QXMT5jdyINDQ3auHGjsrOzPcdCQkKUnZ2tNWvWnNRrWJalG2+8UZdcckmrpU6S5syZo/j4eM8Hl207l81m0xVDkxQbFabYqIBfAgoAQIfy6WJXVlYml8ulxMTEZscTExNVXFx8Uq/xySefaOHChVq8eLEyMjKUkZGhL774osXx9957ryorKz0fe/fuPaVzQNv9NKuf8n55sSaf09d0FAAA/ErAT4lccMEFcrvdJz0+MjJSkZGRXkyE1oSE2JQQG2U6BgAAfsenZ+x69eql0NBQlZSUNDteUlKipKQkQ6nQmb4sqtRP/vipShz1pqMAAODzfLrYRUREaOTIkcrLy/Mcc7vdysvL05gxY7z63na7XWlpacrMzPTq+6BllmXpoTe/1Oqvv9Vj731lOg4AAD7P+KXY6upq7dy50/N5YWGh8vPz1aNHD6Wmpio3N1dTp07VqFGjNHr0aM2bN081NTW66aabvJorJydHOTk5np0o6Hw2m02PXD1cz364Q/deOcR0HAAAfJ7x2518/PHHGjdu3DHHp06dqpdeekmSNH/+fD355JMqLi5WRkaGnnnmGWVlZXVKvrZsMQYAAOhobekixoudr6PY+Za95bVK6dHFdAwAADpNwNzHziTW2Pmep5Zu09inPtYHBSWtDwYAIAhR7FqQk5OjgoICrV+/3nQUfKfR5ZbLbemTr8tMRwEAwCcZ3zwBnKy7Lxukcwf21LjBCaajAADgk5ixg9+ICg+l1AEAcAIUO/ilugaXnsnboRpnk+koAAD4DC7FtsBut8tut8vlcpmOguP4+Svr9cnOb1Ve06DZVw01HQcAAJ/AjF0L2Dzh2/7nooHqEx+liwf1Nh0FAACfwYwd/NJFg3rro1+NVVR4qOkoAAD4DGbs4LeOLnUuN/fZBgCAYge/t2pHmbLnLteW/ZWmowAAYBTFrgU8ecJ//HX9HhWW1WjeBztMRwEAwCieFdsKnhXr+8qqnXp+xTe669IzFRPJslEAQGBpSxfhX0H4vV5dI3XflWeZjgEAgHFcikXAWVdYzmYKAEBQotghoMx+60td+39r9OInhaajAADQ6Sh2CCiDk2Jls0nlNQ2mowAA0OlYY9cCHinmn67LTFFGSjed1YeNLgCA4MOu2FawKxYAAJjUli7CpVgErPKaBt3zj80qqqgzHQUAgE7BpVgErF//43N9sLVU39Y49aep3GgaABD4KHYIWLMmDFFplVPTLznTdBQAADoFxQ4B64yEWL2Zc75sNpvpKAAAdArW2CGgHV3qHPWNBpMAAOB9FDsEhTfz9+uiJz7Ski0HTEcBAMBrKHYtsNvtSktLU2Ymi+4DwfaSKlXUNuqv6/aajgIAgNdwH7tWcB+7wFDf6NLC9Xv1k6xUhYfy+wwAwH+0pYuweQJBISo8VFPP6286BgAAXsXUBYKOZVla/Nl+VTubTEcBAKBDUewQdO5btEUzFubr8fe+Mh0FAIAORbFD0Jk4oo+iwkPUt3u06SgAAHQo1tgh6Jx3Ri+tuucS9eoaaToKAAAdihk7BCVKHQAgEFHsENT2fFurG15Yq837KkxHAQDglFHsENT+kLdDK3eU6aG3vhS3dAQA+DvW2LXAbrfLbrfL5XKZjgIvuu/KIaptaNI9Vwxp9lxZAAD8EU+eaAVPngAAACa1pYtwKRY4yt7yWrnc/K4DAPBPFDvgO39Zs0vZc5frxU8KTUcBAKBdKHbAd0JDQuRscuvTb8rZSAEA8EtsngC+8+PRKeodG6nssxLYSAEA8EsUO+A7NptNl6Ulmo4BAEC7cSkWOI5Gl1sLln+t4sp601EAADhpFDvgOP7foi167L2v9P8Wf8F6OwCA36DYAcdxy4UD1Ds2UlcO72M6CgAAJ401dsBxDEqM1ap7xikyLNR0FAAAThozdkALji513LQYAOAPKHZAKzbvq9APnlmpvK0lpqMAAHBCFDugFe9uPqCviqv09LLtbKQAAPg01ti1wG63y263y+VymY4Cw2ZkD1Kjy1LOuIHcuBgA4NNsFlMQJ+RwOBQfH6/KykrFxcWZjgMAAIJMW7oIl2KBNtq8r0J1DczkAgB8D8UOaIPnPv5ak+yfaO7720xHAQDgGBQ7oA0GJ3WV25LKaxrZSAEA8DlsngDa4JIhiXpn+gUadlq86SgAAByDGTugjSh1AABfRbED2qnG2aTfvP2lth5wmI4CAIAkih3Qbo/+e6te/GSXZv1zs9w8cgwA4AModkA73XnpmUrrE6dfXj5YISHcuBgAYB6bJ4B2SoyL0rt3XsDTKAAAPoMZO+AUHF3qapxNXJIFABhFsQM6wIrtB3XZ3OV6fd0e01EAAEGszcWusbFRYWFh2rJlizfyAH7p64PVKqqs12tr9zBrBwAwps1r7MLDw5WamiqXi2dlAkf8bEx/udyWfpKVykYKAIAx7boUe//99+u+++5TeXl5R+cB/FJoiE0/v/B0dYlgPxIAwJx2/Ss0f/587dy5U8nJyerXr59iYmKafX3Tpk0dEg7wVx9+VaIRfbupV9dI01EAAEGkXcVu0qRJHRwDCBy/f3+7/pC3Qz8Y0Uf2n5xjOg4AIIi0q9g99NBDHZ3DayoqKpSdna2mpiY1NTXprrvu0q233mo6FgLYZWmJem751+rbLVout6VQ1twBADqJzbKsdm/h27hxo7Zu3SpJGjp0qM4+++wOC9ZRXC6XnE6nunTpopqaGg0bNkwbNmxQz549T+rPOxwOxcfHq7KyUnFxcV5Oi0BR6qhXQlyU6RgAgADQli7Srhm70tJSXXfddfr444/VrVs3SYdnxsaNG6e//e1v6t27d3te1itCQ0PVpUsXSZLT6ZRlWTqFLgucFEodAMCEdu2KnT59uqqqqvTll1+qvLxc5eXl2rJlixwOh+688842vdaKFSs0ceJEJScny2azafHixceMsdvt6t+/v6KiopSVlaV169a16T0qKiqUnp6uvn37aubMmerVq1eb/jzQXgernMp5fZNW7SgzHQUAEATaVeyWLFmi//3f/9VZZ53lOZaWlia73a733nuvTa9VU1Oj9PR02e3243594cKFys3N1UMPPaRNmzYpPT1d48ePV2lpqWdMRkaGhg0bdsxHUVGRJKlbt276/PPPVVhYqNdff10lJSXtOGug7f648hu9u/mAHnhzi1zcuBgA4GXtuhTrdrsVHh5+zPHw8HC53e42vdaECRM0YcKEFr8+d+5c3XrrrbrpppskSQsWLNC7776rP//5z5o1a5YkKT8//6TeKzExUenp6Vq5cqWuueaa445xOp1yOp2ezx0Ox0meCXCsOy89U98crNGM7DPZRAEA8Lp2zdhdcskluuuuuzwzYpK0f/9+3X333br00ks7LFxDQ4M2btyo7Oxsz7GQkBBlZ2drzZo1J/UaJSUlqqqqkiRVVlZqxYoVGjx4cIvj58yZo/j4eM9HSkrKqZ0EglrXyDD9aeooDTst3nQUAEAQaFexmz9/vhwOh/r376+BAwdq4MCBGjBggBwOh5599tkOC1dWViaXy6XExMRmxxMTE1VcXHxSr7F7925deOGFSk9P14UXXqjp06dr+PDhLY6/9957VVlZ6fnYu3fvKZ0DcLTSqnrVN/I4PgCAd7TrUmxKSoo2bdqkDz74QF999ZUk6ayzzmo2s+YrRo8efdKXaiUpMjJSkZE8LQAd753NRbrvX19o6nn99cvLW541BgCgvdpc7BobGxUdHa38/Hxddtlluuyyy7yRS5LUq1cvhYaGHrPZoaSkRElJSV57X+nwTly73S6Xi9kVdIxQm02O+iat2lmmuy49U2Gh7ZowBwCgRW3+lyU8PFypqamdUngiIiI0cuRI5eXleY653W7l5eVpzJgxXn3vnJwcFRQUaP369V59HwSPCcP7aMH15+iN/xlDqQMAeEW7/nW5//77dd9996m8vPyUA1RXVys/P99zubSwsFD5+fnas2ePJCk3N1d//OMf9fLLL2vr1q26/fbbVVNT49klC/iTK4b1odQBALymXY8UO/vss7Vz5041NjaqX79+iomJafb1TZs2nfRrffzxxxo3btwxx6dOnaqXXnpJ0uHNGk8++aSKi4uVkZGhZ555RllZWW2N3S48UgzeYFmWFq7fq8wBPTSwd1fTcQAAPszrjxSbNGlSe/7YcY0dO7bVR3xNmzZN06ZN67D3BEx7etl2zf9opzL7d9fCX4xRCPe4AwB0gDYXu6amJtlsNt18883q27evNzL5BDZPwJuuG52iv63fo/FDk8TzKAAAHaVdl2JjY2P1xRdfqH///l6I5Fu4FAtvqW90KSo81HQMAICPa0sXafeTJ5YvX96ucAAOO7rUud1Wq0sSAABoTbvW2E2YMEGzZs3SF198oZEjRx6zeeKqq67qkHBAMCgsq9E9/9ysKaNS9MORgbu8AQDgfe0qdnfccYckae7cucd8zWazBcS6NNbYobMs2VKsdYXl2ldeq6sykhXO7VAAAO3UrjV2wYQ1dvC2Rpdbj7y7VbdcMEApPbqYjgMA8DFeW2N35ZVXqrKy0vP5Y489poqKCs/n3377rdLS0tqWFghy4aEhmn3VUEodAOCUtanYLV26VE6n0/P5o48+2uzpE01NTdq2bVvHpQOC0M7SalXWNpqOAQDwQ20qdt+/astVXKBjLVy/R1f+YaV+926B6SgAAD/EKu0W2O12paWlKTMz03QUBJGBvbuq0e3WwWqnGprcpuMAAPxMmzZPhIaGqri4WL1795Z0+EbFmzdv1oABAyRJJSUlSk5ODqidpGyeQGf7fG+FRvSNl83GY8YAAF58VqxlWbrxxhsVGRkpSaqvr9dtt93muY/d0evvALRPeko30xEAAH6qTcVu6tSpzT6//vrrjxnzs5/97NQSAZAkNTS5tWD51zr/jJ4a2a+H6TgAAD/AfexawaVYmDJ32TY98+FOnd47Ru/ddaEiw3iuLAAEI68/KxaA991ywekakhSruy49UxE8jQIAcBLa9UgxAN4X3yVc/77zQoWEsIkCAHBymAZoAbc7gS84utQ5m1xqcnELFABAy1hj1wrW2MEXfLbnkH79j826+pzTdMfYM0zHAQB0ItbYAQHmm4M12lFardc+3SNnU+DcJxIA0LFYYwf4gcnnnKbymgZdM7Ivu2MBAC2i2AF+wGaz6daLTjcdAwDg47gUC/ihDbvKte9QrekYAAAfQ7ED/Mwra3bpR/+3Rvct2iL2PgEAjkaxawG3O4GvOv+MXooIDVFibKScTdz+BADwH9zupBXc7gS+qKiiTsndok3HAAB0Am53AgQ4Sh0A4HgodoAfq6pv1P9b/IWWbDlgOgoAwAdQ7AA/9sqa3Xr10z164M0vVd/IjYsBINhxHzvAj91ywQBt3H1IP79ggKLCuXExAAQ7ih3gx6LCQ/XnG9m5DQA4jEuxQABx1DeqxtlkOgYAwBCKHRAgVu44qMvmLteTS7eZjgIAMIRiBwQIy5JKHE6t2HFQdQ1spACAYMQauxbY7XbZ7Xa5XPwDCf9w0aDesv/kHF16VgIbKQAgSPHkiVbw5AkAAGAST54AoCVbivVlUaXpGACATkSxAwLQC6sKddurGzXzjc1qdLlNxwEAdBKKHRCAJqb3Ua+uERo3pLdYbAEAwYPNE0AASoiN0sczx6lrJP+LA0AwYcYOCFBHlzrLssQ+KQAIfBQ7IMCVOur1P3/ZqFfW7DYdBQDgZRQ7IMAtKyjRsoISzX1/O48bA4AAxwIcIMD9ZHSqthVX6SdZqYphzR0ABDR+ygMBLiTEpocnDTMdAwDQCbgUCwSZA5V1Kqt2mo4BAPACih0QRJZ+WazL567QA4u3mI4CAPACih0QRPp2j1Zdo0sHKutVzUYKAAg4rLFrgd1ul91ul8vlMh0F6DBDk+O18H/OVUZKd4WG2EzHAQB0MJvFXUtPyOFwKD4+XpWVlYqLizMdBwAABJm2dBEuxQJByu229Jc1u7RqR5npKACADkKxA4LUi6t36YE3v9Q9/9zMjYsBIEBQ7IAgdV1migYnxuoXF52u6PBQ03EAAB2AzRNAkIqJDNO/77qQTRQAEECYsQOC2NGlzu221NDkNpgGAHCqKHYAtLO0Wj/6vzX6Q95201EAAKeAYgdAO0urtXH3Ib366R456htNxwEAtBNr7ADoimFJmjVhiK5KT1ZcVLjpOACAdqLYAZAk3XbxQNMRAACniEuxAI6xvaRKhWU1pmMAANqIYgegmbc+L9IPnlmpmW98LrebJw4CgD+h2AFo5pzUbooIDVG3LuGqaeCJFADgT1hjB6CZvt276N93XajUHl1ks3HzYgDwJxQ7AMfo1zPGdAQAQDtwKRZAi5xNLs1dtk2LPttnOgoA4CQwYwegRW9s2KdnPtyp+OhwXTI4UfFduMcdAPiyoJmxq62tVb9+/fSrX/3KdBTAb1yXmaJxg3trzuThlDoA8ANBM2P3yCOP6NxzzzUdA/ArYaEhevGm0aZjAABOUlDM2O3YsUNfffWVJkyYYDoK4NfqG12q4lmyAOCzjBe7FStWaOLEiUpOTpbNZtPixYuPGWO329W/f39FRUUpKytL69ata9N7/OpXv9KcOXM6KDEQnPL3VujKZ1bqN28XmI4CAGiB8UuxNTU1Sk9P180336zJkycf8/WFCxcqNzdXCxYsUFZWlubNm6fx48dr27ZtSkhIkCRlZGSoqenYG6kuW7ZM69ev16BBgzRo0CCtXr3a6+cDBCqX263CshpV1zfpUE2DusdEmI4EAPgem2VZPvPMIJvNpkWLFmnSpEmeY1lZWcrMzNT8+fMlSW63WykpKZo+fbpmzZrV6mvee++9evXVVxUaGqrq6mo1Njbql7/8pR588MHjjnc6nXI6nZ7PHQ6HUlJSVFlZqbi4uFM7QcDPvfV5kS4+szcbKQCgEzkcDsXHx59UFzF+KfZEGhoatHHjRmVnZ3uOhYSEKDs7W2vWrDmp15gzZ4727t2rXbt26amnntKtt97aYqk7Mj4+Pt7zkZKScsrnAQSKq9KTKXUA4MN8utiVlZXJ5XIpMTGx2fHExEQVFxd75T3vvfdeVVZWej727t3rlfcB/N2n33yrz/YcMh0DAHAU42vsOtONN97Y6pjIyEhFRkZ6PwzgxxZ/tl8zFubr9N4x+vedFyoqPNR0JACAfHzGrlevXgoNDVVJSUmz4yUlJUpKSjKUCsC4wQlKiI3Uuaf3lMvtM8t0ASDo+XSxi4iI0MiRI5WXl+c55na7lZeXpzFjxnj1ve12u9LS0pSZmenV9wH8UXyXcL2fe7EevXq4YiKDauIfAHya8Z/I1dXV2rlzp+fzwsJC5efnq0ePHkpNTVVubq6mTp2qUaNGafTo0Zo3b55qamp00003eTVXTk6OcnJyPDtRADQXH918E4VlWbLZbIbSAAAkHyh2GzZs0Lhx4zyf5+bmSpKmTp2ql156SVOmTNHBgwf14IMPqri4WBkZGVqyZMkxGyoAmFFZ16g5/96q1J5ddMfYM0zHAYCg5lP3sfNFbbl3DBCM3szfr7v+lq/IsBCtnnWJenZl8xEAdKS2dBHjM3a+ym63y263y+VymY4C+LSr0pO1fle5rko/jVIHAIYxY9cKZuwAAIBJAfPkCQD+p6K2QSWOetMxACAoUewAdJg1X3+r7Lkr9Ks3PhcXAwCg81HsAHSYhLhIOeobVVRRp29rGkzHAYCgw+aJFrB5Ami7gb276i83j1Z6SjceMwYABrB5ohVsngAAACaxeQKAcZZl6Z3NRcrbWtL6YABAh+BSLACvWJy/X3cv/Fy9ukbqg9zu6tYlwnQkAAh4zNgB8IoJw/poSFKsfpqVqi4R/A4JAJ2Bn7YtYPMEcGqiwkP19vQLFB7K748A0FnYPNEKNk8AHcOyLDW5LYoeALQRmycA+JQDlXX6+csb9Mi7W01HAYCARrED4HU7S6uV91Wp/rpuj0p53BgAeA1r7AB43YVn9tasCUN06ZAEJcRFmY4DAAGLYgegU9x28UDTEQAg4HEptgV2u11paWnKzMw0HQUIOPsO1WpHSZXpGAAQcNgV2wp2xQIda8X2g7r91Y3q1zNGb047n12yANAKdsUC8FlpyXGKCAtRl4hQVdQ2mo4DAAGFNXYAOlWvrpH6x+3naUDPGIWE2EzHAYCAQrED0OkG9u5qOgIABCQuxQIwxuW29OdVhfrLp7tNRwGAgMCMHQBjln1ZrN++U6Co8BBdMiRBp3WLNh0JAPwaxQ6AMVcMS9JlaYm6eFBv9eHGxQBwyih2LbDb7bLb7XK5XKajAAHLZrPp+RtGymZjEwUAdATuY9cK7mMHdB6X21K1s0nx0eGmowCAz+A+dgD8zs7Sak1+brXuXpgvft8EgPbhUiwAn2BZlrYWORQZFqK95XVK7dnFdCQA8DsUOwA+4czEWM27LkPnpHZXUjwbKQCgPSh2AHzGlcP7mI4AAH6NNXYAfFJBkUOrd5aZjgEAfoViB8DnrNpRpqvmr9JdC/NVWdtoOg4A+A0uxQLwOaP6d1dqzy4alBCrJrfbdBwA8BsUOwA+Jyo8VP+6/Tx16xJhOgoA+BUuxbbAbrcrLS1NmZmZpqMAQen7pY572wFA63jyRCt48gRgVn2jS/M+2KFGl1sP/Fea6TgA0Ona0kW4FAvAp322p0ILln8tm0368ehUnZHQ1XQkAPBZFDsAPm3MwJ66fexAZaR0o9QBQCsodgB83j1XDDEdAQD8ApsnAPiVugaX9pbXmo4BAD6JYgfAb2w94NCEP6zQL/6yUY0u7m8HAN9HsQPgNxJiI1VZ16iK2gbtYdYOAI7BGjsAfqNn10i9cGOmzkjoqriocNNxAMDnUOwA+JVzUrubjgAAPotLsQD81ic7y/TGhr2mYwCAz2DGDoBfWvvNt/rpn9YqKjxEmf17qH+vGNORAMA4ih0AvzR6QA9dcEYvnd47Rr1iI03HAQCfQLED4JdsNptevClT4aGsKAGAI/iJ2AK73a60tDRlZmaajgKgBd8vdQ1N3NsOQHCzWZZlmQ7hyxwOh+Lj41VZWam4uDjTcQAcR0Vtgx5+Z6sO1TbohamjZLPZTEcCgA7Tli7CpVgAfq+s2qm3Py9So9utL/ZXakTfbqYjAYARFDsAfu+MhFg9PGmozkjoSqkDENQodgACwpTMVNMRAMA4Nk8ACDjlNQ36bM8h0zEAoNNR7AAElIIihy6bu1y/+MtGVdY2mo4DAJ2KYgcgoJzeO0bduoSrR5cIldU4TccBgE7FGjsAASUqPFQv3TRaiXFRigjjd1cAwYViByDgpPToYjoCABjBr7MAApZlWVr82X49vWyb6SgA0CmYsQMQsL4scmjGwnxJ0iVDEnR2anezgQDAyyh2AALWsNPidcO5/ZQYF6lhp8WbjgMAXkexAxDQHp40zHQEAOg0rLEDEDTcbktf7Ks0HQMAvIZiByAo1DW4dMvL6/XD51br870VpuMAgFdQ7AAEhajwEIWHhshmk/YdqjMdBwC8gjV2AIKCzWbTU9ema195ndKS40zHAQCvCIpi179/f8XFxSkkJETdu3fXRx99ZDoSAAPiosKVlhzu+dzlthQaYjOYCAA6VlAUO0lavXq1unbtajoGAB+xt7xWd7y2SdMuOUPjhyaZjgMAHYI1dgCC0l/X7dEX+yv16L+3qtHlNh0HADqE8WK3YsUKTZw4UcnJybLZbFq8ePExY+x2u/r376+oqChlZWVp3bp1bXoPm82miy++WJmZmXrttdc6KDkAf3b3ZYM0dUw/vX7ruQoPNf6jEAA6hPFLsTU1NUpPT9fNN9+syZMnH/P1hQsXKjc3VwsWLFBWVpbmzZun8ePHa9u2bUpISJAkZWRkqKmp6Zg/u2zZMiUnJ2vVqlU67bTTdODAAWVnZ2v48OEaMWKE188NgO8KDw3Rb/6bmxcDCCw2y7Is0yGOsNlsWrRokSZNmuQ5lpWVpczMTM2fP1+S5Ha7lZKSounTp2vWrFltfo+ZM2dq6NChuvHGG4/7dafTKafT6fnc4XAoJSVFlZWViotjJx0QqDbvq9BXxVW6dlSK6SgA0IzD4VB8fPxJdRGfvv7Q0NCgjRs3Kjs723MsJCRE2dnZWrNmzUm9Rk1NjaqqqiRJ1dXV+vDDDzV06NAWx8+ZM0fx8fGej5QUfsgDgW5HSZWueW6N7v3XF9qwq9x0HABoN+OXYk+krKxMLpdLiYmJzY4nJibqq6++OqnXKCkp0dVXXy1JcrlcuvXWW5WZmdni+HvvvVe5ubmez4/M2AEIXGckdNWE4Umqcbo0KCnWdBwAaDefLnYd4fTTT9fnn39+0uMjIyMVGRnpxUQAfI3NZtMT14xQeEiIQrivHQA/5tOXYnv16qXQ0FCVlJQ0O15SUqKkJO/ed8putystLe2Es3sAAkdkWGizUre9pMpgGgBoH58udhERERo5cqTy8vI8x9xut/Ly8jRmzBivvndOTo4KCgq0fv16r74PAN8z9/3tuvz3K/Rm/n7TUQCgTYxfiq2urtbOnTs9nxcWFio/P189evRQamqqcnNzNXXqVI0aNUqjR4/WvHnzVFNTo5tuuslgagCBzO0+fLOAgiKH/jvjNMNpAODkGS92GzZs0Lhx4zyfH9m4MHXqVL300kuaMmWKDh48qAcffFDFxcXKyMjQkiVLjtlQAQAd5e7LBmn0gB66aFBv01EAoE186j52vsRut8tut8vlcmn79u3cxw4IYkd+TNpsbKwA0Pnach87il0r2vIfE0DgqWtw6f7FX+j0XjGadsmZpuMACEJt6SLGL8UCgC/7aFup/rVpv8JDbbr6nL46rVu06UgA0CKKHQCcwJXD++i2iwfqokG9KHUAfB7FDgBaMWvCENMRAOCk+PR97EziBsUAjqes2qkXVhWajgEAx8XmiVaweQLAEXUNLmXPXa79FXV6bPJwXTc61XQkAEGgLV2EGTsAOEnREaH6SVaqTu8Vo5H9upuOAwDHYMauFczYATiaZVmqbXApJpIlygA6BzN2AOAlNputWanbW14rl5vfjwH4BopdC9g8AaA17xeU6Mo/rNTTy7aZjgIAkih2LcrJyVFBQYHWr19vOgoAH1XX6FKVs0kbdh1So8ttOg4AcB87AGivq9KTFRkWokuGJCg8lN+TAZhHsQOAUzB+aFKzzy3Lks1mM5QGQLDjV0wA6ACWZenl1buU+/fPxc0GAJhCsQOADvBNWY0efqdAiz7brw+/KjUdB0CQ4lIsAHSAgb276sGJaWpocuuSIQmm4wAIUhS7FtjtdtntdrlcLtNRAPiJn43pbzoCgCDHkydawZMnALRHk8ut19bu0Y9HpyoijFUvANqvLV2EGTsA8II7XtukZQUlKiyr0eyrhpqOAyBI8GskAHjBj0alKDYyTFkDepxwnNttqb7R1ewGx/WNLm0vqdL2kqpmYz/bc0jvbj6g3d/WeI59W+3U/368U39c8U2zsa+t3a3chfn6eFtps7G/eGWDcl7b1Gzsn1cV6oYX1mrxZ/s9xxz1jfrRgtW65rnVch/1yLQ/rfxG//XsSr28elezvNlzl+vSpz9WXcN/lq88v+JrXfD4h/r9+9ubvd8Fj3+o8x/7UN9WOz3HXl69S+c+mqdH3i1oNnbcUx9r9CMfaN+hWs+xv63bo1G/e1/3Lfqi2djsuct1zsPva8dR/90WfbZPZ/92me7622fNxv7gmZXK+O0ybd5X4Tm2ZMsBZfx2mX7xyoZmY3/43Gql/2aZ1n7zrefYR9tKNWL2Ul3/p7XNxl7/p7UaMXtps//un37zrUbMXqofPre62difv7xBI2Yv1XtfHPAc+3xvhYbPXqoJf1jZbOy01zdp+Oyl+temfZ5j20uqNHz2Uo176uNmY2e+8bmGz16q19bu9hzbW16r4bOX6txH85qNffDNLRo+e6n+tPI/f38OVjk1fPZSpf9mWbOxc/69VcNnL9X8D3d4jlU7mzR89lINn71U9Y3/+d7//v3tGj57qZ5c+pXnmMttecZW1DZ4ji9Y/rWGz16q377d/Hs/6nfva/jspTpQWec59vLqXRo+e6nu/Vfz7/0Fj3+o4bOX6puD1Z5jf1+/V8NnL9XdC/Objb1s7nINn71UXxZVeo699XmRhs9eqtv+srHZ2Kvmr9Lw2Uu1cXe559j7BSUaPnuppv55nXwRM3YA4AWXpSVqxa/HqXtMhCQpf2+Fpv91k5Ljo7Xwf8Z4xk19cZ1W7ijT76ek6+qz+0o6/A/2VfM/UXJ8lFbfe6ln7HMff61lBSX63aRh6tczRpJ0qLZBTyzZpvjocN160emesRt2HdKiz/ZrSJ9YjR18eDOHs8mtZQUlx1wa3nmwWit3lGlUv/+UUJfL0vpdhyRJR6/XKa6s15b9Dl1wRn3z1yg9/A+q66jVPVX1Tdp3qE6VdY3Nxu6vqJNlSUc/Yrfa2aRiR/0xY0sd9appcMl91IM96hpdKqtuUFV9U7OxFbUNKq9paJahocmtQ7WNqnE2H+uob1RFbaOajgrR4LJUUduo6u+Nra5vUmVdY7NnAje5LDnqm44d62ySo75JTa7/jHW5vxv7vby1DYfHNh71ui7LUlV90zF56xpcqqpvavYLgPu7sZFhoc3HNn43tunYse7vPde4/ruxDUe9rqXDY0O+dzvGI2OdTc2fsvL974N0+O9aVX2TnI3HH3v0IrCG78bWNzVf0+6ob1JDk/v4Yxubj612Nh0+v6PHug6PPfqXjWZjj4rW9N3Y2hZe13W8sQ3HnrcvYI1dK1hjB6AjbNx9SD98brVSe3TRil+P8xy/5aX1yvuqVI//cLimZKZKknaUVGnK858qITZSS2Zc5Bn79LJtWvtNuW6+YICuGHb4xshl1U499t5X6hoZ1uySb97WEu0srda5p/dUeko3SVKNs0lv5hcpLMSmazNTPGM/23NIu76t0Vl94jQk6fDPuYYmt/K2lshmky5PS1LId//K7yyt0t5DdUrt0UUDe3eVdHjWcW1huWw2aVS/7gr77ikcByrrVOJwqlfXCPXt3sXzfkdmyc7qE+d5YkdpVb1KHU516xLebOxXxQ653dLAhBhPgSmvaVBpVb1io8J1Wrdoz9ivD1bL7baU0qOLosIPj62sa9TBqnp1iQhT8lFjd5XVqMntVt/u/xlbVd+oEodT0RGhzV53b3mtGlxu9YmPUpeIMM9/y2JHvSLDQprl3V9RJ2ejS4lxUYqJPDy2rsGlA5V1Cg8NUUqP/4wtqqhTfaNLvWMjFRsVLulweSqqOHZscWW9ahua1Cs2UnHfjXU2uVRUUa9Qm02pPf8z9kgZ7hETofjocM/3c39FnWyS+veK8Yw9WOVUtbNJ3buEq1uXw7+ENLnc2nvo8CzZgKPGllU7VVXfpG7R4Z5fWNxuS7vLD8+m9uvRxfP3pLymQZV1jYqLClPPrpGSDt/rcde3h8em9uii0O/GHqppUEVdo2KjwtTru7FHvkeWpL7doz1/TyprG1Ve26CYyFAlxEZ5xu7+tkZuSzqtW7TnF5fKukaV1zSoS0SoEuP+M3Zvea2a3Jb6xEc1+96XVTcoOjxUSfEnHlvjbFJplVORYSHN/k55U1u6CMWuFRQ7AB2h2tmk7SVViokI0+CkWM/xqvrDM1TR4aGeQgQAR2PzRAfgdicAOlLXyDCdk9r9mONHZmoAoCMwY9cKZuwAAIBJbekizPsDAAAECIodAABAgKDYAQAABAiKHQAAQICg2AEAAAQIih0AAECAoNgBAAAECIodAABAgKDYtcButystLU2ZmZmmowAAAJwUnjzRCp48AQAATOLJEwAAAEGIYgcAABAgKHYAAAABgmIHAAAQIMJMB/B1R/aWOBwOw0kAAEAwOtJBTma/K8WuFVVVVZKklJQUw0kAAEAwq6qqUnx8/AnHcLuTVrjdbhUVFSk2NlY2m80r7+FwOJSSkqK9e/cG1S1VOG/OOxgE63lLwXvunDfn3dEsy1JVVZWSk5MVEnLiVXTM2LUiJCREffv27ZT3iouLC6r/GY7gvIML5x18gvXcOe/g4u3zbm2m7gg2TwAAAAQIih0AAECAoNj5gMjISD300EOKjIw0HaVTcd6cdzAI1vOWgvfcOW/O2yQ2TwAAAAQIZuwAAAACBMUOAAAgQFDsAAAAAgTFzqD9+/fr+uuvV8+ePRUdHa3hw4drw4YNpmN5lcvl0gMPPKABAwYoOjpaAwcO1MMPP3xSj0nxNytWrNDEiROVnJwsm82mxYsXN/u6ZVl68MEH1adPH0VHRys7O1s7duwwE7YDnei8Gxsbdc8992j48OGKiYlRcnKyfvazn6moqMhc4A7S2vf7aLfddptsNpvmzZvXafm85WTOe+vWrbrqqqsUHx+vmJgYZWZmas+ePZ0ftgO1dt7V1dWaNm2a+vbtq+joaKWlpWnBggVmwnagOXPmKDMzU7GxsUpISNCkSZO0bdu2ZmPq6+uVk5Ojnj17qmvXrvrhD3+okpISQ4k7RmvnXV5erunTp2vw4MGKjo5Wamqq7rzzTlVWVnZ6VoqdIYcOHdL555+v8PBwvffeeyooKNDTTz+t7t27m47mVY8//riee+45zZ8/X1u3btXjjz+uJ554Qs8++6zpaB2upqZG6enpstvtx/36E088oWeeeUYLFizQ2rVrFRMTo/Hjx6u+vr6Tk3asE513bW2tNm3apAceeECbNm3Sv/71L23btk1XXXWVgaQdq7Xv9xGLFi3Sp59+quTk5E5K5l2tnffXX3+tCy64QEOGDNHHH3+szZs364EHHlBUVFQnJ+1YrZ13bm6ulixZoldffVVbt27VjBkzNG3aNL311ludnLRjLV++XDk5Ofr000/1/vvvq7GxUZdffrlqamo8Y+6++269/fbbeuONN7R8+XIVFRVp8uTJBlOfutbOu6ioSEVFRXrqqae0ZcsWvfTSS1qyZIluueWWzg9rwYh77rnHuuCCC0zH6HQ/+MEPrJtvvrnZscmTJ1s//elPDSXqHJKsRYsWeT53u91WUlKS9eSTT3qOVVRUWJGRkdZf//pXAwm94/vnfTzr1q2zJFm7d+/unFCdoKXz3rdvn3XaaadZW7Zssfr162f9/ve/7/Rs3nS8854yZYp1/fXXmwnUSY533kOHDrV++9vfNjt2zjnnWPfff38nJvO+0tJSS5K1fPlyy7IO/xwLDw+33njjDc+YrVu3WpKsNWvWmIrZ4b5/3sfz97//3YqIiLAaGxs7MZllMWNnyFtvvaVRo0bpRz/6kRISEnT22Wfrj3/8o+lYXnfeeecpLy9P27dvlyR9/vnnWrVqlSZMmGA4WecqLCxUcXGxsrOzPcfi4+OVlZWlNWvWGEzW+SorK2Wz2dStWzfTUbzK7Xbrhhtu0MyZMzV06FDTcTqF2+3Wu+++q0GDBmn8+PFKSEhQVlbWCS9TB4rzzjtPb731lvbv3y/LsvTRRx9p+/btuvzyy01H61BHLjX26NFDkrRx40Y1NjY2+9k2ZMgQpaamBtTPtu+fd0tj4uLiFBbWuU9vpdgZ8s033+i5557TmWeeqaVLl+r222/XnXfeqZdfftl0NK+aNWuWrrvuOg0ZMkTh4eE6++yzNWPGDP30pz81Ha1TFRcXS5ISExObHU9MTPR8LRjU19frnnvu0Y9//OOAf7bk448/rrCwMN15552mo3Sa0tJSVVdX67HHHtMVV1yhZcuW6eqrr9bkyZO1fPly0/G86tlnn1VaWpr69u2riIgIXXHFFbLb7broootMR+swbrdbM2bM0Pnnn69hw4ZJOvyzLSIi4phf1ALpZ9vxzvv7ysrK9PDDD+sXv/hFJ6eTOrdGwsPtdmvUqFF69NFHJUlnn322tmzZogULFmjq1KmG03nP3//+d7322mt6/fXXNXToUOXn52vGjBlKTk4O6PPGsRobG3XttdfKsiw999xzpuN41caNG/WHP/xBmzZtks1mMx2n07jdbknSf//3f+vuu++WJGVkZGj16tVasGCBLr74YpPxvOrZZ5/Vp59+qrfeekv9+vXTihUrlJOTo+Tk5GazWf4sJydHW7Zs0apVq0xH6VStnbfD4dAPfvADpaWlafbs2Z0bTszYGdOnTx+lpaU1O3bWWWf5/U6x1sycOdMzazd8+HDdcMMNuvvuuzVnzhzT0TpVUlKSJB2zU6ykpMTztUB2pNTt3r1b77//fsDP1q1cuVKlpaVKTU1VWFiYwsLCtHv3bv3yl79U//79Tcfzml69eiksLCzoftbV1dXpvvvu09y5czVx4kSNGDFC06ZN05QpU/TUU0+Zjtchpk2bpnfeeUcfffSR+vbt6zmelJSkhoYGVVRUNBsfKD/bWjrvI6qqqnTFFVcoNjZWixYtUnh4eKdnpNgZcv755x+zRXz79u3q16+foUSdo7a2ViEhzf/ahYaGen6zDxYDBgxQUlKS8vLyPMccDofWrl2rMWPGGEzmfUdK3Y4dO/TBBx+oZ8+epiN53Q033KDNmzcrPz/f85GcnKyZM2dq6dKlpuN5TUREhDIzM4PuZ11jY6MaGxsD8medZVmaNm2aFi1apA8//FADBgxo9vWRI0cqPDy82c+2bdu2ac+ePX79s62185YO/wy//PLLFRERobfeesvYzm8uxRpy991367zzztOjjz6qa6+9VuvWrdPzzz+v559/3nQ0r5o4caIeeeQRpaamaujQofrss880d+5c3Xzzzaajdbjq6mrt3LnT83lhYaHy8/PVo0cPpaamasaMGfrd736nM888UwMGDNADDzyg5ORkTZo0yVzoDnCi8+7Tp4+uueYabdq0Se+8845cLpdn3U2PHj0UERFhKvYpa+37/f0CGx4erqSkJA0ePLizo3ao1s575syZmjJlii666CKNGzdOS5Ys0dtvv62PP/7YXOgO0Np5X3zxxZo5c6aio6PVr18/LV++XK+88ormzp1rMPWpy8nJ0euvv64333xTsbGxnv9/4+PjFR0drfj4eN1yyy3Kzc1Vjx49FBcXp+nTp2vMmDE699xzDadvv9bO+0ipq62t1auvviqHwyGHwyFJ6t27t0JDQzsvbKfuwUUzb7/9tjVs2DArMjLSGjJkiPX888+bjuR1DofDuuuuu6zU1FQrKirKOv30063777/fcjqdpqN1uI8++siSdMzH1KlTLcs6fMuTBx54wEpMTLQiIyOtSy+91Nq2bZvZ0B3gROddWFh43K9Jsj766CPT0U9Ja9/v7wuU252czHm/8MIL1hlnnGFFRUVZ6enp1uLFi80F7iCtnfeBAwesG2+80UpOTraioqKswYMHW08//bTldrvNBj9FLf3/++KLL3rG1NXVWXfccYfVvXt3q0uXLtbVV19tHThwwFzoDtDaebf090GSVVhY2KlZbd8FBgAAgJ9jjR0AAECAoNgBAAAECIodAABAgKDYAQAABAiKHQAAQICg2AEAAAQIih0AAECAoNgBAAAECIodAABAgKDYAQAABAiKHQAAQIAIMx0AAALR2LFjNWLECEVFRelPf/qTIiIidNttt2n27NmmowEIYMzYAYCXvPzyy4qJidHatWv1xBNP6Le//a3ef/9907EABDCbZVmW6RAAEGjGjh0rl8ullStXeo6NHj1al1xyiR577DGDyQAEMmbsAMBLRowY0ezzPn36qLS01FAaAMGAYgcAXhIeHt7sc5vNJrfbbSgNgGBAsQMAAAgQFDsAAIAAQbEDAAAIEOyKBQAACBDM2AEAAAQIih0AAECAoNgBAAAECIodAABAgKDYAQAABAiKHQAAQICg2AEAAAQIih0AAECAoNgBAAAECIodAABAgKDYAQAABAiKHQAAQID4/+4f3Ml/emfxAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fig, ax = plt.subplots()\n", + "ax.set_xlabel('n')\n", + "ax.set_ylabel('Error')\n", + "ax.set_xticks(n_qubits[::2], n_qubits[::2])\n", + "ax.plot(n_qubits, errors, linestyle='dotted')\n", + "ax.set_yscale('log')\n", + "plt.tight_layout()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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": 2 +}