{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b286a1fb-e2e9-4156-994b-aba2f2aca4e4",
   "metadata": {},
   "source": [
    "# Conserving first integrals via manifold projection\n",
    "\n",
    "As we saw in an [earlier example](<./Outer Solar System.ipynb>), even though heyoka.py is not a [symplectic integrator](https://en.wikipedia.org/wiki/Symplectic_integrator) it is nevertheless capable of optimally preserving conserved quantities in dynamical systems.\n",
    "\n",
    "In order to conserve dynamical invariants, heyoka.py needs to be configured to operate at very low tolerances (i.e., below machine epsilon). High-precision numerical integrations however are computationally expensive, and in some contexts it is preferable, for performance reasons, to maintain a low integration accuracy while at the same time ensuring that the global dynamical invariants are preserved. This is the main reason why symplectic integrators are the tool of choice, for instance, in numerical studies of the long-term stability of extrasolar planetary systems.\n",
    "\n",
    "In this example, we will show how we can enforce the conservation of dynamical invariants in high-tolerance numerical integrations with heyoka.py via a technique known as *manifold projection*.\n",
    "\n",
    "## Setting things up\n",
    "\n",
    "We will re-use here the setup described in the [outer Solar System example](<./Outer Solar System.ipynb>): a gravitational 6-body problem consisting of the Sun, Jupiter, Saturn, Uranus, Neptune and Pluto with realistic initial conditions.\n",
    "\n",
    "Let us begin by defining a few constants and the initial conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fb22cfa7-1e47-4068-af57-fcc78f200717",
   "metadata": {},
   "outputs": [],
   "source": [
    "import heyoka as hy\n",
    "import numpy as np\n",
    "\n",
    "# Masses, from Sun to Pluto.\n",
    "masses = np.array(\n",
    "    [1.00000597682, 1 / 1047.355, 1 / 3501.6, 1 / 22869.0, 1 / 19314.0, 7.4074074e-09]\n",
    ")\n",
    "\n",
    "# The gravitational constant.\n",
    "G = 0.01720209895 * 0.01720209895 * 365 * 365\n",
    "\n",
    "# Initial conditions.\n",
    "ic = [  # Sun.\n",
    "    -4.06428567034226e-3,\n",
    "    -6.08813756435987e-3,\n",
    "    -1.66162304225834e-6,\n",
    "    +6.69048890636161e-6 * 365,\n",
    "    -6.33922479583593e-6 * 365,\n",
    "    -3.13202145590767e-9 * 365,\n",
    "    # Jupiter.\n",
    "    +3.40546614227466e0,\n",
    "    +3.62978190075864e0,\n",
    "    +3.42386261766577e-2,\n",
    "    -5.59797969310664e-3 * 365,\n",
    "    +5.51815399480116e-3 * 365,\n",
    "    -2.66711392865591e-6 * 365,\n",
    "    # Saturn.\n",
    "    +6.60801554403466e0,\n",
    "    +6.38084674585064e0,\n",
    "    -1.36145963724542e-1,\n",
    "    -4.17354020307064e-3 * 365,\n",
    "    +3.99723751748116e-3 * 365,\n",
    "    +1.67206320571441e-5 * 365,\n",
    "    # Uranus.\n",
    "    +1.11636331405597e1,\n",
    "    +1.60373479057256e1,\n",
    "    +3.61783279369958e-1,\n",
    "    -3.25884806151064e-3 * 365,\n",
    "    +2.06438412905916e-3 * 365,\n",
    "    -2.17699042180559e-5 * 365,\n",
    "    # Neptune.\n",
    "    -3.01777243405203e1,\n",
    "    +1.91155314998064e0,\n",
    "    -1.53887595621042e-1,\n",
    "    -2.17471785045538e-4 * 365,\n",
    "    -3.11361111025884e-3 * 365,\n",
    "    +3.58344705491441e-5 * 365,\n",
    "    # Pluto.\n",
    "    -2.13858977531573e1,\n",
    "    +3.20719104739886e1,\n",
    "    +2.49245689556096e0,\n",
    "    -1.76936577252484e-3 * 365,\n",
    "    -2.06720938381724e-3 * 365,\n",
    "    +6.58091931493844e-4 * 365,\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd9258a0-6372-44d4-84ce-8fa8ba242694",
   "metadata": {},
   "source": [
    "See the [outer Solar System example](<./Outer Solar System.ipynb>) for more details on the setup.\n",
    "\n",
    "Next, we setup the dynamical equations of our N-body problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8935f6c9-a178-4f18-a39e-dfb2ba429879",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the ODEs.\n",
    "sys = hy.model.nbody(6, masses=masses, Gconst=G)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c95d9d5-c9c7-469b-b138-164fa6360c78",
   "metadata": {},
   "source": [
    "We also store the 36 state variables (3 Cartesian positions and 3 Cartesian velocities for each of the 6 bodies) in a list for later use:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "95ef6549-b7fc-4ab6-9104-ca7b719609d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "state_vars = [p[0] for p in sys]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f4482bd-efd0-44ec-b102-af4cd084cb4c",
   "metadata": {},
   "source": [
    "Next, we define a [compiled function](./compiled_functions.ipynb) for the rapid computation of the energy of the system from the state vector. Here we are using the ``model.nbody_energy()`` helper, which creates the symbolic formula for the energy constant in the N-body problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e7edc7d0-10de-485b-b782-b2cf8a600b80",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2024-02-12 07:53:01.408] [heyoka] [info] heyoka logger initialised\n"
     ]
    }
   ],
   "source": [
    "# Define a compiled function for the computation\n",
    "# of the energy in an N-body system from its state vector.\n",
    "en_cfunc = hy.cfunc(\n",
    "    [hy.model.nbody_energy(6, masses=masses, Gconst=G)], vars=state_vars\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70c64a69-ae6d-4031-b441-57aa83371036",
   "metadata": {},
   "source": [
    "We can now proceed to the instantiation of the numerical integrator. We choose a high tolerance of $10^{-6}$, so that the integrator will not be able to enforce energy conservation on its own (since we are not operating at sub-epsilon tolerances):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "42a79b43-dbf4-4752-a3b0-82ef39d37229",
   "metadata": {},
   "outputs": [],
   "source": [
    "ta = hy.taylor_adaptive(sys, ic, tol=1e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "358f7349-6b66-45b5-9323-994d163ac797",
   "metadata": {},
   "source": [
    "Let us now compute the initial energy of the system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "feba1bb2-8e22-4572-a2aa-599d37e5afc6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial energy: -0.004286848855986956\n"
     ]
    }
   ],
   "source": [
    "print(f\"Initial energy: {en_cfunc(ta.state)[0]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35760b46-76c3-47e6-9cf0-c5835f3b6338",
   "metadata": {},
   "source": [
    "So far so good (though perhaps not very exciting).\n",
    "\n",
    "The next step is the introduction of a callback to be invoked after every step taken by the numerical integrator. This callback will check the relative energy error, and stop the integration if the error exceeds the $10^{-6}$ threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "107589f8-6449-4d78-a2b4-77824c9073a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "class proj_callback:\n",
    "    def __init__(self, ta):\n",
    "        # Store the initial energy value\n",
    "        # as a data member.\n",
    "        self.E0 = en_cfunc(ta.state)[0]\n",
    "\n",
    "    def __call__(self, ta):\n",
    "        # Compute the relative energy error\n",
    "        # wrt the initial energy value.\n",
    "        rel_err = abs((en_cfunc(ta.state)[0] - self.E0) / self.E0)\n",
    "\n",
    "        # Stop if the error is greater than 1e-6.\n",
    "        if rel_err > 1e-6:\n",
    "            print(f\"Rel energy error: {rel_err}\")\n",
    "            return False\n",
    "        else:\n",
    "            return True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d24ddbf-21db-4769-8931-073ba63ff447",
   "metadata": {},
   "source": [
    "Let us now integrate until the energy error reaches the threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "154f5ff1-2888-4862-8c6f-cc71bc7f4ba5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rel energy error: 1.1480506370240604e-06\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(<taylor_outcome.cb_stop: -4294967301>,\n",
       " 0.7296940138544884,\n",
       " 1.0583722519817735,\n",
       " 36,\n",
       " None,\n",
       " <__main__.proj_callback at 0x7f3ee0402510>)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cb = proj_callback(ta)\n",
    "# Attempt to propagate the state of the\n",
    "# system for 1000 years.\n",
    "ta.propagate_until(1000.0, callback=cb)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef71d8d2-4af5-466e-83e1-d6e5df070d49",
   "metadata": {},
   "source": [
    "We can see from the output that indeed the integration was terminated earlier because of the callback (as indicated by the ``taylor_outcome.cb_stop`` outcome), and that the error threshold was exceeded after 36 steps.\n",
    "\n",
    "## Fixing the energy\n",
    "\n",
    "We can now move to the interesting part, that is, enforcing energy conservation. The basic idea is to alter the current state vector of the system (whose energy has drifted away from the initial value) so that the energy of the updated state vector matches the initial energy of the system.\n",
    "\n",
    "Because there are infinite points in phase space whose energy is equal to the initial value, we impose the additional constraint that the updated state vector should be as close as possible to the original one. In other words, we are projecting the position in phase space of the current state vector onto the hyper-surface (manifold) defined by the conservation of energy.\n",
    "\n",
    "This projection process can be seen as a [constrained optimisation problem](https://en.wikipedia.org/wiki/Constrained_optimization): we want to minimise the distance from the original state vector subject to the constraint that the updated state vector must lie on the manifold defined by the conservation of energy.\n",
    "\n",
    "Let us begin with the definition of the objective function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b586ca8b-c541-4578-8f8a-752e77aa9067",
   "metadata": {},
   "outputs": [],
   "source": [
    "def objfun(x):\n",
    "    return np.sum((x - ta.state) ** 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47545522-c078-422d-8743-3ea846a1a149",
   "metadata": {},
   "source": [
    "As mentioned above, here we are trying to minimise the distance (in phase space) between the original state vector (``ta.state``) and the updated state vector (``x``).\n",
    "\n",
    "Next, we define the constraint:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "56e535a9-8c1f-434b-8960-6326504e2941",
   "metadata": {},
   "outputs": [],
   "source": [
    "def cstr_fun(x):\n",
    "    return en_cfunc(x)[0] - cb.E0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21374e6c-d85d-4ffc-a414-e93fffa5cbee",
   "metadata": {},
   "source": [
    "This equality constraint is satisifed if ``cstr_fun()`` returns zero, that is, if the energy of the updated state vector (``x``) is equal to the initial energy of the system (stored in the callback as ``cb.E0``).\n",
    "\n",
    "We are now ready to run the optimisation via {func}`scipy.optimize.minimize()`. Note that we need an optimiser with support for constraints (here we choose ``'SLSQP'``):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e645b715-b9b9-4638-bfeb-dc837ed0848d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " message: Optimization terminated successfully\n",
      " success: True\n",
      "  status: 0\n",
      "     fun: 1.6028683862604742e-12\n",
      "       x: [ 3.841e-04  4.208e-03 ... -1.215e+00  5.875e-02]\n",
      "     nit: 2\n",
      "     jac: [ 3.113e-08 -4.027e-07 ...  7.715e-11  7.409e-11]\n",
      "    nfev: 76\n",
      "    njev: 2\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "ret = minimize(\n",
    "    objfun,\n",
    "    ta.state,\n",
    "    method=\"SLSQP\",\n",
    "    tol=1e-10,\n",
    "    constraints=[{\"type\": \"eq\", \"fun\": cstr_fun}],\n",
    ")\n",
    "\n",
    "# Store the updated state vector.\n",
    "res = ret.x\n",
    "print(ret)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d827eeae-6bf1-48ab-9de8-6bd700f277fe",
   "metadata": {},
   "source": [
    "We can see how the optimisation terminated successfully.\n",
    "\n",
    "Let us verify that the updated state vector (``res``) conserves energy with a better accuracy than the original state vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "80eeec14-deb2-4bd7-916b-afb8889b45f8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original energy error: 1.1480506370240604e-06\n",
      "Updated  energy error: 1.3963053728083164e-11\n"
     ]
    }
   ],
   "source": [
    "print(f\"Original energy error: {abs((en_cfunc(ta.state)[0] - cb.E0) / cb.E0)}\")\n",
    "print(f\"Updated  energy error: {abs((en_cfunc(res)[0] - cb.E0) / cb.E0)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ae906fd-2926-4e4f-a09c-ed8082f9b96c",
   "metadata": {},
   "source": [
    "## Helping the optimiser out\n",
    "\n",
    "Optimisation algorithms can greatly benefit from the availability of gradients for the objective function and the constraints. Luckily, heyoka.py's expression system lets us compute derivatives without too much pain.\n",
    "\n",
    "Let us begin by defining a symbolic expression representing the square of the distance between the original state vector and the updated state vector. Here, the original state vector will be represented by the symbolic state variables (stored earlier in ``state_vars``) and the updated state vector will be passed in the array of [runtime parameters](<./ODEs with parameters.ipynb>):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8ecc3adc-e466-49e6-af59-72bac1f022f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist2_ex = hy.sum(list((state_vars[i] - hy.par[i]) ** 2 for i in range(36)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "263abf8d-54de-4be7-9f4e-26a85b8a3288",
   "metadata": {},
   "source": [
    "Let us take a look at ``dist2_ex`` just to fix the ideas:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "48d502be-e289-40e1-9c86-d9f3608190ed",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((vx_0 - p3)**2.0000000000000000 + (vx_1 - p9)**2.0000000000000000 + (vx_2 - p15)**2.0000000000000000 + (vx_3 - p21)**2.0000000000000000 + (vx_4 - p27)**2.0000000000000000 + (vx_5 - p33)**2.0000000000000000 + (vy_0 - p4)**2.0000000000000000 + (vy_1 - p10)**2.0000000000000000 + (vy_2 - p16)**2.0000000000000000 + (vy_3 - p22)**2.0000000000000000 + (vy_4 - p28)**2.0000000000000000 + (vy_5 - p34)**2.0000000000000000 + (vz_0 - p5)**2.0000000000000000 + (vz_1 - p11)**2.0000000000000000 + (vz_2 - p17)**2.0000000000000000 + (vz_3 - p23)**2.0000000000000000 + (vz_4 - p29)**2.0000000000000000 + (vz_5 - p35)**2.0000000000000000 + (x_0 - p0)**2.0000000000000000 + (x_1 - p6)**2.0000000000000000 + (x_2 - p12)**2.0000000000000000 + (x_3 - p18)**2.0000000000000000 + (x_4 - p24)**2.0000000000000000 + (x_5 - p30)**2.0000000000000000 + (y_0 - p1)**2.0000000000000000 + (y_1 - p7)**2.0000000000000000 + (y_2 - p13)**2.0000000000000000 + (y_3 - p19)**2.0000000000000000 + (y_4 - p25)**2.0000000000000000 + (y_5 - ..."
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dist2_ex"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fc0a78e-34f9-4fe9-8c79-2074a581486e",
   "metadata": {},
   "source": [
    "We can now compute symbolically the gradient vector of ``dist2_ex``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "a574bdb2-38de-48dd-9fb7-3a9e92e889cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[((2.0000000000000000 * x_0) - (2.0000000000000000 * p0)),\n",
       " ((2.0000000000000000 * y_0) - (2.0000000000000000 * p1)),\n",
       " ((2.0000000000000000 * z_0) - (2.0000000000000000 * p2)),\n",
       " ((2.0000000000000000 * vx_0) - (2.0000000000000000 * p3)),\n",
       " ((2.0000000000000000 * vy_0) - (2.0000000000000000 * p4)),\n",
       " ((2.0000000000000000 * vz_0) - (2.0000000000000000 * p5)),\n",
       " ((2.0000000000000000 * x_1) - (2.0000000000000000 * p6)),\n",
       " ((2.0000000000000000 * y_1) - (2.0000000000000000 * p7)),\n",
       " ((2.0000000000000000 * z_1) - (2.0000000000000000 * p8)),\n",
       " ((2.0000000000000000 * vx_1) - (2.0000000000000000 * p9)),\n",
       " ((2.0000000000000000 * vy_1) - (2.0000000000000000 * p10)),\n",
       " ((2.0000000000000000 * vz_1) - (2.0000000000000000 * p11)),\n",
       " ((2.0000000000000000 * x_2) - (2.0000000000000000 * p12)),\n",
       " ((2.0000000000000000 * y_2) - (2.0000000000000000 * p13)),\n",
       " ((2.0000000000000000 * z_2) - (2.0000000000000000 * p14)),\n",
       " ((2.0000000000000000 * vx_2) - (2.0000000000000000 * p15)),\n",
       " ((2.0000000000000000 * vy_2) - (2.0000000000000000 * p16)),\n",
       " ((2.0000000000000000 * vz_2) - (2.0000000000000000 * p17)),\n",
       " ((2.0000000000000000 * x_3) - (2.0000000000000000 * p18)),\n",
       " ((2.0000000000000000 * y_3) - (2.0000000000000000 * p19)),\n",
       " ((2.0000000000000000 * z_3) - (2.0000000000000000 * p20)),\n",
       " ((2.0000000000000000 * vx_3) - (2.0000000000000000 * p21)),\n",
       " ((2.0000000000000000 * vy_3) - (2.0000000000000000 * p22)),\n",
       " ((2.0000000000000000 * vz_3) - (2.0000000000000000 * p23)),\n",
       " ((2.0000000000000000 * x_4) - (2.0000000000000000 * p24)),\n",
       " ((2.0000000000000000 * y_4) - (2.0000000000000000 * p25)),\n",
       " ((2.0000000000000000 * z_4) - (2.0000000000000000 * p26)),\n",
       " ((2.0000000000000000 * vx_4) - (2.0000000000000000 * p27)),\n",
       " ((2.0000000000000000 * vy_4) - (2.0000000000000000 * p28)),\n",
       " ((2.0000000000000000 * vz_4) - (2.0000000000000000 * p29)),\n",
       " ((2.0000000000000000 * x_5) - (2.0000000000000000 * p30)),\n",
       " ((2.0000000000000000 * y_5) - (2.0000000000000000 * p31)),\n",
       " ((2.0000000000000000 * z_5) - (2.0000000000000000 * p32)),\n",
       " ((2.0000000000000000 * vx_5) - (2.0000000000000000 * p33)),\n",
       " ((2.0000000000000000 * vy_5) - (2.0000000000000000 * p34)),\n",
       " ((2.0000000000000000 * vz_5) - (2.0000000000000000 * p35))]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grad_dist2_ex = list(hy.diff(dist2_ex, state_vars[i]) for i in range(36))\n",
    "grad_dist2_ex"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bfeb689-9a92-43f7-9953-94b4ca4048b2",
   "metadata": {},
   "source": [
    "Shocking, I know.\n",
    "\n",
    "We can now proceed to create a compiled function that will concatenate in a single vector output the value of the distance square and its gradient vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a22e0b7b-fbd9-435d-b875-1f759f220682",
   "metadata": {},
   "outputs": [],
   "source": [
    "objfun_cfunc = hy.cfunc([dist2_ex] + grad_dist2_ex, vars=state_vars)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f69b4e6a-42e6-4d29-9d1a-d2e6155b200f",
   "metadata": {},
   "source": [
    "Finally, we can define an updated objective function that, in addition to returning the distance square, will return also its gradient vector (in the format required by {func}`scipy.optimize.minimize()`, i.e., as a tuple of 2 values):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "6b2785f2-940c-4aa3-914d-0bcdff9133fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "def objfun(x):\n",
    "    ret = objfun_cfunc(x, pars=ta.state)\n",
    "\n",
    "    # Return the distance square and its\n",
    "    # gradient as a tuple of 2 values.\n",
    "    return (ret[0], ret[1:])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "466bc65d-3ca6-4ec8-a76b-fd795228edb4",
   "metadata": {},
   "source": [
    "In a similar fashion, we can define a compiled function for the computation of the gradient of the constraint:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "bd91bf32-f0b9-4027-b87e-d0eb4904a64d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The symbolic expression representing the energy\n",
    "# of the system.\n",
    "en_ex = hy.model.nbody_energy(6, masses=masses, Gconst=G)\n",
    "\n",
    "# The gradient of the energy.\n",
    "grad_cstr = list(hy.diff(en_ex, state_vars[i]) for i in range(36))\n",
    "\n",
    "# The compiled function for the computation of the gradient\n",
    "# of the constraint.\n",
    "grad_cstr_cfunc = hy.cfunc(grad_cstr, vars=state_vars)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0aa7348a-6929-4b87-8053-4a058b277757",
   "metadata": {},
   "source": [
    "We can now put everything together, and invoke again {func}`scipy.optimize.minimize()` supplying the gradients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "bd825b3f-dfc3-4a36-b4da-636583752a46",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " message: Optimization terminated successfully\n",
      " success: True\n",
      "  status: 0\n",
      "     fun: 1.6028685167358091e-12\n",
      "       x: [ 3.841e-04  4.208e-03 ... -1.215e+00  5.875e-02]\n",
      "     nit: 2\n",
      "     jac: [ 3.120e-08 -4.048e-07 ...  2.930e-12 -1.417e-13]\n",
      "    nfev: 4\n",
      "    njev: 2\n"
     ]
    }
   ],
   "source": [
    "ret = minimize(\n",
    "    objfun,\n",
    "    ta.state,\n",
    "    method=\"SLSQP\",\n",
    "    tol=1e-10,\n",
    "    jac=True,\n",
    "    constraints=[{\"type\": \"eq\", \"fun\": cstr_fun, \"jac\": grad_cstr_cfunc}],\n",
    ")\n",
    "\n",
    "res = ret.x\n",
    "print(ret)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a65817a3-4e70-463f-8ecc-b5a8702bd7f9",
   "metadata": {},
   "source": [
    "Like before, the optimisation completed successfully, this time taking far fewer objective function evaluations thanks to the availability of the gradient (which was previously estimated via numerical differencing).\n",
    "\n",
    "We can again verify that the updated state vector (``res``) conserves energy with a better accuracy than the original state vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "1c22fb72-40d4-4584-9689-d2025a1390e0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original energy error: 1.1480506370240604e-06\n",
      "Updated  energy error: 1.3842666877172304e-11\n"
     ]
    }
   ],
   "source": [
    "print(f\"Original energy error: {abs((en_cfunc(ta.state)[0] - cb.E0) / cb.E0)}\")\n",
    "print(f\"Updated  energy error: {abs((en_cfunc(res)[0] - cb.E0) / cb.E0)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d35a9602-3aad-4cb7-a255-caeffa308994",
   "metadata": {},
   "source": [
    "## Packaging it all together\n",
    "\n",
    "As a last step, we now incorporate manifold projection into the step callback, so that we can automatically ensure energy conservation whenever the energy error drifts too far during the numerical integration.\n",
    "\n",
    "Here is the updated version of the callback:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "b5f4670c-0846-45dd-9e7a-694e3660e397",
   "metadata": {},
   "outputs": [],
   "source": [
    "class proj_callback:\n",
    "    def __init__(self, ta):\n",
    "        # Store the initial energy value\n",
    "        # as a data member.\n",
    "        self.E0 = en_cfunc(ta.state)[0]\n",
    "\n",
    "        # Setup a variable to count how many times\n",
    "        # we perform manifold projection.\n",
    "        self.n_proj = 0\n",
    "\n",
    "    def __call__(self, ta):\n",
    "        # Compute the relative energy error\n",
    "        # wrt the initial energy value.\n",
    "        rel_err = abs((en_cfunc(ta.state)[0] - self.E0) / self.E0)\n",
    "\n",
    "        # Fix the energy if the error is greater than 1e-6.\n",
    "        if rel_err > 1e-6:\n",
    "            ret = minimize(\n",
    "                objfun,\n",
    "                ta.state,\n",
    "                method=\"SLSQP\",\n",
    "                tol=1e-10,\n",
    "                jac=True,\n",
    "                constraints=[{\"type\": \"eq\", \"fun\": cstr_fun, \"jac\": grad_cstr_cfunc}],\n",
    "            )\n",
    "            if ret.status != 0:\n",
    "                print(f\"Optimisation returned nonzero status:\\n{ret}\")\n",
    "                return False\n",
    "\n",
    "            # Update the state vector.\n",
    "            ta.state[:] = ret.x\n",
    "\n",
    "            # Update the counter.\n",
    "            self.n_proj += 1\n",
    "\n",
    "        return True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa59d71d-90ba-48a1-a996-16bab3ec6af0",
   "metadata": {},
   "source": [
    "Let us now reset the state of the integrator back to the original initial conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "4c15265c-14eb-468c-8607-7a1257461976",
   "metadata": {},
   "outputs": [],
   "source": [
    "ta.state[:] = ic\n",
    "ta.time = 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2a59596-1dcd-49d8-8ff0-79895f196b55",
   "metadata": {},
   "source": [
    "We are now ready to repeat the numerical integration with automatic manifold projection. We will propagate the state of the system over an equispaced time grid up to $10^5$ years:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "04a71b6e-fca4-4627-9027-e6a42a8be1ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the callback object.\n",
    "cb = proj_callback(ta)\n",
    "\n",
    "# Define the time grid.\n",
    "t_grid = np.linspace(0, 1e5, 1000)\n",
    "\n",
    "# Run the numerical integration.\n",
    "pres_proj = ta.propagate_grid(t_grid, callback=cb)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7de946e-6d52-42cf-818e-6837bbf082b3",
   "metadata": {},
   "source": [
    "Let us also run a numerical integration without projection for comparison:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "069bab3c-ae06-43c7-802e-c5b506910692",
   "metadata": {},
   "outputs": [],
   "source": [
    "ta.state[:] = ic\n",
    "ta.time = 0.0\n",
    "\n",
    "pres_noproj = ta.propagate_grid(t_grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d66be2a1-5a7f-4a24-a62b-f9468a97135d",
   "metadata": {},
   "source": [
    "We can now take a look at the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "7f3640d4-4857-470d-aefe-52160b20ee00",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAJOCAYAAABm7rQwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAC+YElEQVR4nOzdeZwU1b3//3fPKDiDMICgDJssBlFBRCARQWA0V8Etikkw5oohajTGJGqMgl6XqBGMUXEJRsVozC9Rkitmue5fBQGJikQiGJeAQzA4BkdhRocRdKZ/f4zVVFfX3tXd1d2v5+OR66Wnl9PVVafO+ZzPOSeRTCaTAgAAAAAAAPKootAFAAAAAAAAQPkhKAUAAAAAAIC8IygFAAAAAACAvCMoBQAAAAAAgLwjKAUAAAAAAIC8IygFAAAAAACAvCMoBQAAAAAAgLwjKAUAAAAAAIC8263QBShW7e3tevfdd9W1a1clEolCFwcAAAAAACAWksmkPvroI/Xt21cVFc75UASlQnr33Xc1YMCAQhcDAAAAAAAglt555x3179/f8e8EpULq2rWrpI4D3K1btwKXBgAAAAAAIB6am5s1YMCAVOzECUGpkIwpe926dSMoBQAAAAAAYOG13BELnQMAAAAAACDvCEoBAAAAAAAg7whKAQAAAAAAIO9YUyqH2tvbtXPnzkIXAwhs9913V2VlZaGLAQAAAAAoYQSlcmTnzp2qr69Xe3t7oYsChNK9e3f16dPHc2E6AAAAAADCICiVA8lkUg0NDaqsrNSAAQNUUcEsSRSPZDKp7du3a8uWLZKk2traApcIAAAAAFCKCErlwGeffabt27erb9++qq6uLnRxgMCqqqokSVu2bNHee+/NVD4AAAAAQORI4cmBtrY2SVKnTp0KXBIgPCOg+umnnxa4JAAAAACAUkRQKodYiwfFjPMXAAAAAJBLBKWQlaVLlyqRSGjbtm2uzxs0aJDmz5+flzL5MWXKFF1wwQU5/Qy/xwYAAAAAgHJEUAqSpF/+8pfq2rWrPvvss9RjH3/8sXbffXcdccQRac9dvny5EomE3nrrLR1++OFqaGhQTU2NJOn+++9X9+7d81n0UBYvXqxrr702svezC3JZjw0AAAAAANiFoBQkSXV1dfr444/18ssvpx5bvny5+vTpo1WrVmn79u2px5cuXaq+fftq2LBh6tSpk/r06RObqV47d+709byePXuqa9euOS1L3I4NAAAAAABxQlAKkqT9999fffv21dKlS1OPLV26VF/5ylc0dOhQrVy5Mu3xurq61P9vTFFbunSpZs2apaamJiUSCSUSCV199dWp123fvl3f/va31bVrVw0cOFB33323a5mmTJmi888/X+eff766d++uvfbaS//zP/+jZDKZes6gQYN03XXX6Vvf+pZqamp09tlnS5IefvhhHXTQQercubMGDRqkm266KeO9zZlNO3fu1CWXXKJ+/fqpS5cu+tKXvpR2LCTp+eef1+TJk1VdXa0ePXromGOO0datW/Wtb31Lzz33nG699dbU9964caPt9D2vcg0aNEjXX399oOMEAAAAAEAxIigVcw1NrVq5oVENTa05/6wpU6ZoyZIlqX8vWbJEU6ZM0eTJk1OP79y5U3/9619TQSmzww8/XPPnz1e3bt3U0NCghoYGXXzxxam/33TTTRo7dqxeeeUVnXfeefrud7+rN954w7VMv/71r7XbbrvpxRdf1G233aZbbrlFCxcuTHvOjTfeqBEjRmj16tW64oortHr1an3961/XqaeeqrVr1+rqq6/WFVdcofvvv9/xc2bNmqXnn39eDz30kF599VV97Wtf09SpU/XPf/5TkrRmzRodddRROuigg/TXv/5VK1as0AknnKC2tjbdeuutGj9+vM4+++zU9x4wYEDGZ/gtV5jjBAAAAABAsdmt0AWAs0WrNmnO4rVqT0oVCWnu9JGaMW5gzj5vypQpuvDCC/XZZ5+ptbVVr7zyiiZNmqS2tjbddtttkqQXXnhBra2ttkGpTp06qaamRolEQn369Mn4+7HHHqvzzjtPknTppZfqlltu0dKlSzV8+HDHMg0YMEC33HKLEomE9t9/f61du1a33HJLKiNKko488si04Nc3v/lNHXXUUbriiiskScOGDdM//vEP3XjjjfrWt76V8RkbNmzQgw8+qH//+9/q27evJOniiy/WE088ofvuu0/XX3+9fvazn2ns2LFasGBB6nUHHXRQ2nevrq62/d6Gm2++2Ve5whwnAAAAAACKDZlSMdXQ1JoKSElSe1K6bPG6nGZM1dXVqaWlRatWrdLy5cs1bNgw7b333po8ebJWrVqllpYWLV26VAMHDtSQIUMCv//BBx+c+v+NwNWWLVtcX3PYYYelrck0fvx4/fOf/1RbW1vqsbFjx6a95vXXX9eECRPSHpswYULG6wx/+9vflEwmNWzYMO25556p/z333HPasGGDpF2ZUtnwW64wxwkAAAAAgGJDplRM1Te2pAJShrZkUhsbt6u2pionn7nffvupf//+WrJkibZu3arJkydLkvr06aPBgwfr+eef15IlS3TkkUeGev/dd9897d+JRELt7e1Zl7tLly5p/04mkxmLi5vXobJqb29XZWWlVq9ercrKyrS/7bnnnpKkqqrsj7nfcuXqOAEAAAAAECdkSsXU4F5dVGHZtK0ykdCgXtU5/dy6ujotXbpUS5cu1ZQpU1KPT548WU8++aReeOEF26l7hk6dOtlmI4X1wgsvZPz7C1/4QkbwyOzAAw/UihUr0h5buXKlhg0bZvu60aNHq62tTVu2bNF+++2X9j9jOt7BBx+sZ555xvEz/XzvoOUCAAAAAKCUkSkVU7U1VZo7faQuW7xObcmkKhMJXT99RM6ypAx1dXX63ve+p08//TSVKSV1BKW++93v6pNPPnENSg0aNEgff/yxnnnmGY0aNUrV1dWqrg4fSHvnnXd00UUX6ZxzztHf/vY33X777Rk71ln96Ec/0rhx43TttddqxowZ+utf/6o77rgjbT0os2HDhumb3/ymZs6cqZtuukmjR49WY2Ojnn32WY0cOVLHHnus5syZo5EjR+q8887Tueeeq06dOmnJkiX62te+pl69emnQoEF68cUXtXHjRu25557q2bNn1uUCAAAAAMRQ02bpnRel7R/m5v2re0oDviTV9MvN+8cIQakYmzFuoCYN662Njds1qFd1zgNSUkdQqrW1VcOHD9c+++yTenzy5Mn66KOPNHToUNud5QyHH364zj33XM2YMUMffPCBrrrqKl199dWhyzNz5ky1trbqi1/8oiorK/X9739f3/nOd1xfc+ihh+r3v/+9rrzySl177bWqra3VNddcY7vIueG+++7Tddddpx/96EfavHmz9tprL40fP17HHnuspI7A1VNPPaXLLrtMX/ziF1VVVaUvfelL+sY3viGpY2H0M844QwceeKBaW1tVX18fSbkAAAAAoCzlOvAT1jsvSmt/n4cPSkgn3iYdOjMPn1U4iaTbYjtw1NzcrJqaGjU1Nalbt25pf/vkk09UX1+vwYMHa4899ihQCYvflClTdMghh2j+/PmRv/f48eN11FFH6brrrov8vUsF5zEAAABQJuIWAMpb4CfmEhXSBeuKMmPKLWZiRqYUysqOHTu0du1avfbaa/rBD35Q6OIAAAAAKEdxCgIRAIqvZLv04dtFGZTyi6AUysrjjz+umTNn6oQTTtBXv/rVQhcHAAAAQD4QBEIxSlRIPYcUuhQ5RVAKsbV06dLI3/Okk05Sc3Nz5O8LAAAAlK2wAZ/qnlL3faVt/8ptsIggEIpRIiGdcGtJZ0lJBKUAAAAAIN7ilOVjRcAH5Wzk16UBh0X/vtU9pQFfLPmAlERQCgAAAEA5i3PARyLoA+Qq8JONMgoa5RpBKQAAAAC5F8fgDwEfIF2cAkAEfsoCQSkAAACglMUhGETwB7AXlyAQASAUCEEpAAAAIJcKGRQiGASkIwgExApBKQAAAJSXfAaJCAqhnAQN+FT3lLoPlLZtyv31SBAIiCWCUsjK0qVLVVdXp61bt6p79+6Ozxs0aJAuuOACXXDBBXkrW6EkEgk98sgjOumkk3L2GVdffbX++Mc/as2aNTn7DAAACiLXASOCRChmccnysco24NN/TLTlAVA0CEpBkvTLX/5SP/7xj7V161bttlvHafHxxx+rR48eOuyww7R8+fLUc5cvX65JkybpzTff1OGHH66GhgbV1NRIku6//35dcMEF2rZtWyG+RiyCXw0NDerRo0dk72cX5Lr44ov1/e9/P7LPAAAgkFwFjggYoZDiGvCRyPIBULIISkGSVFdXp48//lgvv/yyDjus42a8fPly9enTR6tWrdL27dtVXV0tqSM7qm/fvho2bJgkqU+fPgUrdz59+umn2n333T2fl4/jseeee2rPPffM+ecAAEpA1AEkAkfIVtyCPwR8AKBgCEpBkrT//vurb9++Wrp0aSootXTpUn3lK1/RkiVLtHLlSn35y19OPV5XV5f2/2/dulVr1qzRrFmzJHVk90jSVVddpauvvlqStH37dn3729/WH/7wB/Xo0UP/8z//o+985zupMqxdu1Y//OEP9de//lXV1dU65ZRTdPPNN6eCL1OmTNEhhxyi+fPnp15z0kknqXv37rr//vs1ZcoU/etf/9KFF16oCy+8UJKUTCZtv28ikdCCBQv05z//WUuXLlWfPn30s5/9TF/72tckSRs3btTgwYO1aNEiLViwQC+88ILuvPNOnXHGGbruuut099136/3339cBBxygefPmaerUqWnvbc5s2rx5sy666CI99dRTqqio0MSJE3Xrrbdq0KBBqdf86le/0k033aT169erZ8+eOuWUU3THHXeknnPyySdLkvbdd19t3LgxY/pee3u7a7mM7/Pwww/r9ttv14svvqgvfOEL+uUvf6nx48f7OEMAAHkTZRCJABLMCh0MIvgDALAgKBV3TZulDzdIPYfm/AY+ZcoULVmyRLNnz5YkLVmyRJdccona29u1ZMkSffnLX9bOnTv117/+VbfffnvG6w8//HDNnz9fV155pd58801JSsvmuemmm3Tttdfqsssu0//+7//qu9/9riZNmqThw4dr+/btmjp1qg477DCtWrVKW7Zs0VlnnaXzzz9f999/v6/yL168WKNGjdJ3vvMdnX322Z7Pv+KKKzRv3jzdeuut+s1vfqNvfOMbGjFihA444IDUcy699FLddNNNuu+++9S5c2fdeuutuummm3TXXXdp9OjR+tWvfqUTTzxRr732mr7whS9kfMb27dtVV1enI444QsuWLdNuu+2m6667TlOnTtWrr76qTp066c4779RFF12kefPmadq0aWpqatLzzz8vSVq1apX23ntv3XfffZo6daoqKyttv4vfcl1++eX6+c9/ri984Qu6/PLL9Y1vfEPr169PTdkEAIQUVSCJIFJpK1RQiGAQACCmyrYn+tFHH+nII4/Up59+qra2Nv3gBz/wFcjIq789IP3lh1KyXUpUSCfcKh06M2cfN2XKFF144YX67LPP1NraqldeeUWTJk1SW1ubbrvtNknSCy+8oNbW1lSmlFmnTp1UU1OjRCJhO4Xt2GOP1XnnnSepI9hzyy23aOnSpRo+fLh++9vfqrW1VQ888IC6dOkiSbrjjjt0wgkn6IYbbtA+++zjWf6ePXuqsrJSXbt29TWF7mtf+5rOOussSdK1116rp59+WrfffrsWLFiQes4FF1yg6dOnp/7985//XJdeeqlOPfVUSdINN9ygJUuWaP78+frFL36R8RkPPfSQKioqtHDhwlT22H333afu3btr6dKlOvroo3XdddfpRz/6kX74wx+mXjdu3DhJUu/evSVJ3bt3d/1Ofst18cUX67jjjpMk/eQnP9FBBx2k9evXa/jw4Z7HCwBKDoEk5CtIRFAIAABbZRuUqq6u1nPPPafq6mpt375dI0aM0PTp07XXXnsVumgdmjbvCkhJHf/9ywXS0KNy1qCpq6tTS0uLVq1apa1bt2rYsGHae++9NXnyZJ1++ulqaWnR0qVLNXDgQA0ZMiTw+x988MGp/98IXG3ZskWS9Prrr2vUqFGpgJQkTZgwQe3t7XrzzTd9BaWCsk5bGz9+fMZudmPHjk39/83NzXr33Xc1YcKEtOdMmDBBf//7320/Y/Xq1Vq/fr26du2a9vgnn3yiDRs2aMuWLXr33Xd11FFHhf4eQcpl/g1qa2slSVu2bCEoBaC4RBFMIpAUf7kMGBEkAgAgFso2KFVZWZlauPuTTz5RW1ub4/pDBfHhhl0BKUOyTfrw7Zw1oPbbbz/1799fS5Ys0datWzV58mRJHQt3Dx48WM8//7yWLFmiI488MtT7WxcJTyQSam/v+I7JZDKVSWRlPF5RUZHxG3366aehyuLEWgZzkMzpOW5lb29v15gxY/Tb3/4242+9e/dWRUVFFqUNXi7zb2D8zfgNACAvsg0oEUyKl1wEjggYAQBQNmIblFq2bJluvPFGrV69Wg0NDWkLRxsWLFigG2+8UQ0NDTrooIM0f/58HXHEEb4/Y9u2bZo8ebL++c9/6sYbb1SvXr0i/hZZ6Dm0Y8qeOTCVqJR6Bs9QCqKurk5Lly7V1q1b9eMf/zj1+OTJk/Xkk0/qhRdeSC1mbqdTp05qa2sL/LkHHnigfv3rX6ulpSUVCHr++edVUVGR2uWvd+/eamhoSL2mra1N69atS5tKGOTzX3jhBc2cOTPt36NHj3Z8frdu3dS3b1+tWLFCkyZNSj2+cuVKffGLX7R9zaGHHqpFixZp7733Vrdu3WyfM2jQID3zzDO2UyKljkCS23cKUy4ACIWAUnGLMoBE4AgAAEQgtkGplpYWjRo1SrNmzdIpp5yS8fdFixbpggsu0IIFCzRhwgTdddddmjZtmv7xj39o4MCBkqQxY8Zox44dGa996qmn1LdvX3Xv3l1///vf9Z///EfTp0/XV7/61ZxMEwulpl/HGlJ/uaAjQypRKZ0wP+eNv7q6On3ve9/Tp59+msqUkjqCUt/97nf1ySefOAZPpI4Ay8cff6xnnnlGo0aNUnV1dSojzc03v/lNXXXVVTrjjDN09dVX6/3339f3v/99nX766anf5Mgjj9RFF12kRx99VEOHDtUtt9yibdu2ZXz+smXLdOqpp6pz586ugcY//OEPGjt2rCZOnKjf/va3eumll3Tvvfe6lvPHP/6xrrrqKg0dOlSHHHKI7rvvPq1Zs8Y2E8r4XjfeeKO+8pWv6JprrlH//v21adMmLV68WD/+8Y/Vv39/XX311Tr33HO19957a9q0afroo4/0/PPP6/vf/37qOz3zzDOaMGGCOnfurB49emRdLgBlLGxgiYBS/kUVRCKABAAAYiq2Qalp06Zp2rRpjn+/+eabdeaZZ6YWqp4/f76efPJJ3XnnnZo7d66kjvV8/Nhnn3108MEHa9myZfra175m+5wdO3akBbiam5v9fpXwDp3ZsYbUh293ZEjloTFZV1en1tZWDR8+PC1AN3nyZH300UcaOnSoBgwY4Pj6ww8/XOeee65mzJihDz74QFdddZWuvvpqz8+trq7Wk08+qR/+8IcaN26cqqurdcopp+jmm29OPefb3/62/v73v2vmzJnabbfddOGFF2YEyK655hqdc845Gjp0qHbs2OE6JfMnP/mJHnroIZ133nnq06ePfvvb3+rAAw90LecPfvADNTc360c/+pG2bNmiAw88UH/+859td94zvteyZct06aWXavr06froo4/Ur18/HXXUUanMqTPOOEOffPKJbrnlFl188cXq1auXvvrVr6be46abbtJFF12ke+65R/369dPGjRuzLheAEhAmuERgKfeiCCQRRAIAAGUikYzVQkr2EolE2vS9nTt3qrq6Wn/4wx908sknp573wx/+UGvWrNFzzz3n+Z7/+c9/VFVVpW7duqm5uVnjx4/Xgw8+mLYQtNnVV1+tn/zkJxmPNzU1ZUzL+uSTT1RfX6/Bgwdrjz32CPBNkS/WcyoqO3bs0B577KGnn35aX/7ylyN973zjPAbyhOBS4RFIAgAAiFRzc7NqampsYyZmsc2UctPY2Ki2traMqXb77LOP3nvvPV/v8e9//1tnnnmmksmkksmkzj//fMeAlCTNmTNHF110Uerfzc3NrhlDKD/Nzc1avHixKioq2M0OKEdBg0vVPaX/vCYt/3luy1XKsg0mEUgCAAAoqKIMShmC7IJmNWbMGK1Zs8b3Z3Xu3FmdO3cOUjyUmauuukq/+93vdMMNN6h///6FLg6AbAUJMpG5FFw2ASWCSQAAACWhKINSvXr1UmVlZUZW1JYtW+KzUDliLRezVm+55Rbdcsstkb8vgIgQZIoOASUAAABEoCiDUp06ddKYMWP09NNPp60p9fTTT+srX/lKAUsGAMgbgkzZCRNYIqAEAACACMU2KPXxxx9r/fr1qX/X19drzZo16tmzpwYOHKiLLrpIp59+usaOHavx48fr7rvv1qZNm3TuuecWsNQAgKz5CTYRZNolaHCJwBIAAABiIrZBqZdffll1dXWpfxuLjJ9xxhm6//77NWPGDH3wwQe65ppr1NDQoBEjRuixxx7TvvvuW6giZyiCjQ0BR5y/iJTfrKayDTYlpJFfI7gEAACAspJI0vMMxW17w08//VTr169X3759VVNTU6ASAtn54IMPtGXLFg0bNkyVlZWFLg7iiqymTH4zl6p7St0HSp9ul3oOIbgEAACAkuEWMzGLbaZUMdttt91UXV2t999/X7vvvrsqKioKXSTAt2Qyqe3bt2vLli3q3r07AalyRbBplyBBJjKXAAAAAN8ISuVAIpFQbW2t6uvr9a9//avQxQFC6d69u/r06VPoYiBqBJsIMgEAAAAxQVAqRzp16qQvfOEL2rlzZ6GLAgS2++67kyFVTMp9vSaCTAAAAEBRIiiVQxUVFdpjjz0KXQwAxarcs5q8gk0EmQAAAICiRlAKAPKtnINNfrKaCDYBAAAAZYGgFABEhWCT898JNAEAAACwICgFAF4INjn/nWATAAAAgJAISgEoXwSbnP9OsAkAAABAjhGUAlCavAJOpRhsYr0mAAAAAEWEoBSA4lKO2U1kNQEAAAAoQQSlAMSLW9CJYBMAAAAAlAyCUgDyp1ym1BFsAgAAAABPBKUARKMcAk4EmwAAAAAgMgSlAPhTytPqCDYBAAAAQN4RlALQoVSDTm4BJ4JNAAAAAFAwBKWAcuIUeCrGoBPZTQAAAABQ1AhKAaXILvhUbIEnp6ATwSYAAAAAKAkEpYBiVcxZT0ypAwAAAICyR1AKiLtiy3oi4AQAAAAA8IGgFBAXxRJ8YlodAAAAACACBKWAQrAGoOIUfCLoBAAAAADIA4JSQK7FMQBlF3gi6AQAAAAAyCOCUkDUzEGoQgegrMEnAk8AAAAAgJggKAVkq9BBKLKeAAAAAABFiKAUEEQhp+KR9QQAAAAAKCEEpQAvRiDqzcfzE4Ai+AQAAAAAKAMEpQCrfE7HMwegCD4BAAAAAMoIQSkgX0EoAlAAAAAAAKQQlEJ5yseUPCMIRQAKAAAAAIAMBKVQPnIdiCIIBQAAAACAbwSlUNpyFYhiKh4AAAAAAFkhKIXSk8tA1P7HEoACAAAAACACBKVQOpo2S8tulFbfF837MR0PAAAAAICcISiF4hZlVhRBKAAAAAAA8oagFIpTVFlRTMkDAAAAAKAgCEqhuEQRjCIQBQAAAABAwRGUQvxFMUWPQBQAAAAAALFCUArxlW1WFIEoAAAAAABii6AU4un526Snrwj32jGzpEk/JhAFAAAAAECMEZRCvDRtlv7f1cGn6ZEVBQAAAABAUSEohXgIO1WPrCgAAAAAAIoSQSkUFsEoAAAAAADKEkEpFE7QdaOYogcAAAAAQMkgKIXCeOZaafnP/T135NelL19NIAoAAAAAgBJCUAr5FXQh8/+6Rprww5wWCQAAAAAA5B9BKeRPkOl6rBkFAAAAAEBJIyiF3AuSHUUwCgAAAACAskBQCrnlNzuKdaMAAAAAACgrBKWQO34XMz/iYumoALvwAQAAAACAokdQCrnhNyDFQuYAAAAAAJQlglKInp+AFNP1AAAAAAAoawSlEC0/ASmyowAAAAAAKHsEpRAdr4AU2VEAAAAAAOBzBKUQDa+AFIuZAwAAAAAAk4pCFwAlgIAUAAAAAAAIiKAUsvP8rQSkAAAAAABAYASlEF7TZunpq5z/TkAKAAAAAAA4ICiF8JbdKClp/zcCUgAAAAAAwAVBKYTz/K3S6vvs/0ZACgAAAAAAeCAoheDcpu2NmUVACgAAAAAAeCIoheAcp+0lpEk/zndpAAAAAABAESIohWDcpu391zVSTb/8lgcAAAAAABQlglLwz2va3oQf5Lc8AAAAAACgaBGUgn8v3imm7QEAAAAAgCgQlII/TZullXfY/41pewAAAAAAICCCUvDHaXFzpu0BAAAAAIAQCErBm+Pi5kzbAwAAAAAA4RCUgju3xc0P/z7T9gAAAAAAQCgEpeDuww1yXNz8S+fmuzQAAAAAAKBEEJSCu3dfsX+cxc0BAAAAAEAWCErBWdNm6f9dnfn4ERezuDkAAAAAAMgKQSk4e/FOKdme+fiQKXkvCgAAAAAAKC0EpWCvabO08o7MxxMVUs8h+S8PAAAAAAAoKQSlYM9pgfPx57OWFAAAAAAAyBpBKdjrObQjKypNBTvuAQAAAACASOxW6AIgpjY8IyVNmVKJhHTCrWRJAQAAAACASJAphUxNm6W//FBp0/eSCWnoUQUrEgAAAAAAKC0EpZDpww02u+61Sx++XZDiAAAAAACA0kNQCpnefSXzsUQlu+4BAAAAAIDIEJRCuqbN0v+7OvPxL1/NelIAAAAAACAyBKWQznbqnqS+o/NfFgAAAAAAULIISiEdU/cAAAAAAEAeEJTCLkzdAwAAAAAAeUJQCrswdQ8AAAAAAORJWQeldtttNx1yyCE65JBDdNZZZxW6OIXXc6iUsJwSTN0DAAAAAAA5sFuhC1BI3bt315o1awpdjPio6SedcKv0lwukZFtHQOqE+UzdAwAAAAAAkSvroBRsDD1KOmWhpIQ04IsEpAAAAAAAQE7EdvresmXLdMIJJ6hv375KJBL64x//mPGcBQsWaPDgwdpjjz00ZswYLV++PNBnNDc3a8yYMZo4caKee+65iEpexP72gDR/hPS/s6SHvy1teKbQJQIAAAAAACUqtplSLS0tGjVqlGbNmqVTTjkl4++LFi3SBRdcoAULFmjChAm66667NG3aNP3jH//QwIEDJUljxozRjh07Ml771FNPqW/fvtq4caP69u2rdevW6bjjjtPatWvVrVu3nH+3WGraLP3lh7sWOk+2d0zjG3oU2VIAAAAAACByiWQymSx0IbwkEgk98sgjOumkk1KPfelLX9Khhx6qO++8M/XYAQccoJNOOklz584N/BnTpk3Ttddeq7Fjx/p6fnNzs2pqatTU1FQagaz6ZdKvT8h8/Iz/kwYfkf/yAAAAAACAouQ3ZhLb6Xtudu7cqdWrV+voo49Oe/zoo4/WypUrfb3H1q1bU1lU//73v/WPf/xDQ4Y47zK3Y8cONTc3p/2vpLDzHgAAAAAAyKOiDEo1Njaqra1N++yzT9rj++yzj9577z1f7/H6669r7NixGjVqlI4//njdeuut6tmzp+Pz586dq5qamtT/BgwYkNV3iB1j571EZce/2XkPAAAAAADkUGzXlPIjkUik/TuZTGY85uTwww/X2rVrfX/WnDlzdNFFF6X+3dzcXHqBKXbeAwAAAAAAeVKUQalevXqpsrIyIytqy5YtGdlTUencubM6d+6ck/eOhb89sGuh80RFR9bUoTMLXSoAAAAAAFCiinL6XqdOnTRmzBg9/fTTaY8//fTTOvzwwwtUqiLmtPNe0+aCFgsAAAAAAJSu2GZKffzxx1q/fn3q3/X19VqzZo169uypgQMH6qKLLtLpp5+usWPHavz48br77ru1adMmnXvuuQUsdZH6cMOugJQh2SZ9+DZT+AAAAAAAQE7ENij18ssvq66uLvVvYz2nM844Q/fff79mzJihDz74QNdcc40aGho0YsQIPfbYY9p3330LVeTiZey8Zw5MsfMeAAAAAADIoUQymUwWuhDFqLm5WTU1NWpqalK3bt0KXZzs/e2Bjil7ybZdO++xphQAAAAAAAjIb8wktplSyLNDZ3bsvvfh2x0ZUkzbAwAAAAAAOURQCrvU9CMYBQAAAAAA8qIod98DAAAAAABAcSMohQ5Nm6X6ZR3/BQAAAAAAyDGm7+HzRc5/2LH7XqJCOuFWFjkHAAAAAAA5RaZUuWvavCsgJXX89y8XkDEFAAAAAAByiqBUuftww66AlCHZ1rELHwAAAAAAQI4QlCp3PYd2TNkzS1RKPYcUpjwAAAAAAKAsEJQqdzX9OtaQSlR2/DtRKZ0wv+NxAAAAAACAHGGhc3Qsaj70qI4pez2HEJACAAAAAAA5R1AKHWr6EYwCAAAAAAB5w/Q9AAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAAAAAA5B1BKUhNm6X6ZR3/BQAAAAAAyIPdCl0AFNjfHpD+8kMp2S4lKqQTbpUOnVnoUgEAAAAAgBJHplQ5a9q8KyAldfz3LxeQMQUAAAAAAHKOoFQ5+3DDroCUIdkmffh2YcoDAAAAAADKBkGpctZzaMeUPbNEpdRzSGHKAwAAAAAAygZBqXJW069jDalEZce/E5XSCfM7HgcAAAAAAMghFjovd4fOlIYe1TFlr+cQAlIAAAAAACAvCEqhIxBFMAoAAAAAAOQR0/cAAAAAAACQdwSlAAAAAAAAkHcEpQAAAAAAAJB3BKUAAAAAAACQdwSlAAAAAAAAkHcEpQAAAAAAAJB3BKUAAAAAAACQdwSlAAAAAAAAkHcEpQAAAAAAAJB3BKUAAADgqaGpVSs3NKqhqbXQRQEAACVit0IXAACAOGtoalV9Y4sG9+qi2pqqQhcHAfH7RWPRqk2as3it2pNSRUKaO32kZowbWOhiAQD1fIF5HX9+H3ghKIWSRQW4C8cCxcLtXC3Eeey3I841Fk8EUqLR0NSaOo6S1J6ULlu8TpOG9eZ8B5A14x7apVOlWna2BbqXBqnnuVf7E+Q4eR3/fN2Ho/xt43CexKEM+URQCrEX5qIs1Y6I07FwO0bFeCzKrSJGB6dztaGpVb9aUa97V9RHeh77Gdlz6ohLSr122VvvF901Vg7yEUgpl7qqvrEldRwNbcmkNjZuz+p7l8vxA+DMfO83JCTNnjZc50we6vraIPV8XIJXQd87H/Wk+TOCtGm8jn++BjSi7OvEod8UhzLkG0EpxIJThRvmoizkiG4ubxxOx8LtGBV6dDtXAcUwwbk4iqK8UX/nQh1Dp3N1W+unmvfYGzL3h6M4j/2cZ04d8ftWbNTCFW+rPdnRcJaUKl82ZSu28zfuchVIMUTZaMzXbx/2cwb36qKKhNKOZ2UioepOFVq5oTFUuYul0Z3P69Lus+JUL8SpLEEU8r5WLMerUEEY673fkJQ09/E3pIR0ziTnwJTfej7b4NWkYb0jyeIOWu/lo540f0bQNo3X8Xf6+4MvbtI3vjQwsnar9bed8/Da0O2woP2mXLTDyzEzmaAU8iJMJo+f6Lvde+a6I+IklzcOp2MxvE9X12OU62MRdYaWn4o4THAujqIob9TfuZDH0Olcnfd4ekDK/De389h8bhrvb5ynfm/4dh3xCikVkJIUqmx28nnss21AFUtHyymQMqhXtaSO7/Hyxg+VSCQ0Zt8egb6LUyN4eJ+uGjWgR6ByZvPbB5nyks3n1NZUae70kbps8Tq1JZOqTCR00ui+OnnBylBTW6NqdOfiXAybMZBtmex+H0mxua8V2z3WUKhyR/G5+ZrObg1K+MlQCvPedsfB7t5vdsPjb+jEUX0dv6NXPe/2OeZ7tXE8W3d+ptmL1yppqptmP7xWic8/w/odgmZfBan38pXta/6MoG2aLp0qlUgodbykjnNo/ZaPNKhXte3vI0m3Pbtetz+7XvNOyf56tPtt2yXdt2KjLjvugKzfqy2Z1OqNW9Vzz+yTJ/y0OwrVjy00glLImteN8a5lGzo6lgEzedwuSreGot8bVJRyfeNwOharNm51rbicbgavbt6m8UP3yqpMbsHElzd+GOp4+Gk0hAnOxU0U50tU72EO3BQyw/CDj3fYBoCcGqtu17R1KkBCHQ0t4zwd0LPa1w3f2hGvkHTsyD76v7XvuX6foPVNPkfFvDIrvTo5xZQdZBdIuX76CNXWVGnRqk2a/fDaVAM8IQVqHDs1gk9asFLz8hScsZvy4vSbdulUmfU5NmPcQE0a1lsbG7erulNFKiDl9n5O50sUje5cBBuyyRjIpkxOQU6Z6sRcBu7CTGUOG4TNp0JlHETxuW7nUtT1sDUo4ZShFPRckrzbFU7tVEN7Uq71gls9b+Y0yPRByw7dtWyDbnj8DccyJKW0IJV5Gn+Q3zlovZeP4IRXUFBybtMY52HS8vqkpCv+9Jqu/NNrmnfKSM2dPlJzHl6rdmU+b87icBlNZoN7dUm188wWrnhbsyYOCtQmtmuLJiT94KFXAidPWPltdxSiHxsHBKWQFa8b413Pbei4uX0uSCaP23QBtzRNvzeoKOX6xuF0Mx03qIdrxVVbU6VLpw5P+w0k6WePv+k68mRteNj922mqldON3c/x8KqIwwbnCsmuERfF+ZLte1iDxWdNHFyQY2jtBBojbpWJhC6Zur9ueCLzfKpIyPGatpsKYO1ULj5vvO8bvtERN6bs2QWkjHK3f17uoPVNvkbF3BpQfjJCouzguQW1/Qaq/DzXHEgZ1Ks6VX+ZG4ZS8MaxU0cqmafgjNOUF6ff1K7B7vQ5bsfVuMeu3NDoWW6388VudN24Bt2yHJ2+fxTBhmwzBrIpk1OQ01oIu8yOsEFk4/VrNzel7ttBpjJ7BWGjFiaInW3dGjZwHsXnuq1nGOW57xSUuOHxN3TY4J5q2dmmLp0q9X+vNniu62g93/y0K6xtdis/nXG7et7K+jlGnXj+715xfW87xndIKmn7/ZyyaoIGG9zqSTthzlevoKAkXTJtf9/TLs2M++rzs4/UbaeNtj3WXkFHP2prqnT2EYN19/J6z/f2s1yMuS1a8fn3sLvegmRVBWl3LHvr/fTMM5c2bykhKIXQvEbOGppaNc8SDJG8M3mMCtcpuNSys822cWRO0/Rzg4qS3Y1DiiYjSdp1MzVXaElJb7z3kWcAbmT/moz3c2scWSvmLw7uoZfqt/rKNDGCHHb8NCy8AopO54tXcM4sn1OPnDoEUYyCZPMedsHihcvrC55hmJSUSEq/OG20Dv08rbl79e5pmUpnTRqsWRMGO/52XqN+bcmktu9szzjPLpm6v+obWyTJ9r3NU/bMjHM0m/omm98yyPns1oDy08mJKnjmJ6jtNfofJFPACKSYv4fdKdKelGNnwsraaDQLEpwJ+9u7ned2v6ndU40sgYam1tR7+glQSP46S07ny+3PrNdDqzZlvPb66SMyAmlG2a1lyeZcDDL138rtt1n9r/CDI04DT3I4N/yuvegnCG3mdO0HDcJGLWxmUDZ1azbZSNne493Ob6dASNiOvVOWSXuyI+hoV8/ZnSd255tTu8K6Fp25zf78+ve1YOkG10Eeu2vYWs/bMT5n9catqcwXP6zHx/gttzR/Yjt1zS6rxiij30FzuywkP88Pc76eOXFwKuBo5+B+3dP+3dDUqv979V1fx689KT36aoPGDephe55VJBRJO3PWxMFaaPkO1mvO73IxSUkVSemO00bbBi7dkiecfn+3dodde8H83ERSqYB0KSMohUDMNwKvkbMBPavtG8KmCsiugjZ3DicN6635p45SRSKR6qQ2NLX6StP0c4OKglP6quSdkeSHccyH9+madvNLqqNRsGJ2nVbMrnPsEAdpHNlVzC/Wb0393S3TxGuqld8ov1tA0emGPmpAD183+nyv2+PW0c82m8/rPdwWg7cLFrdL+s7EIbp3RX1BMwyTkv7+TpOOO7ivpOABZqcGtsE498cP3Sv1vq9u3ubaEXfqrF5x3AE69uDatDonjLDnQ9C1LOzS0isTiYyOr2TfyYkqpdxx/TDTgvZuo//ZZsk4nSNunQm7z3drj5sHJNw6mOOH7hXqt3cb3a6QtHX7Ttu/Gd/bnCVgnaZmcDqufjtLTmX83UubMsq7+Lzx2rvbHpow71nbQJq1LGHPRbdrxqmD4ScL0piWYeX3+nCqAyTZPubn/HcKklkDllZuU5ntpuAEDYgEHRgKcr1b39tP3ZqLNc+C1Ol2n+92fv95zbsZ75HNAFJtTZVmT8vMqJcyB1rNrL+7U5/A2q5wWovO+L3GD91L3zxsX8d7frbtuNqaKvXc0zsAbWXO4jYH0M3HyC2rxvgefto0dllICXXUk3bTZcNeI9ZBgNO+NEAPvfSO78COX9c9+roqEtL0Q/tp8d8275q+9vnvF0U700+b2OkY3f7sP23P3Z5dOmtQr2rfyRNuv79Tu8MalHO6juIw8yPXCErBN+uN4NKpw11HzuwCF5J06bThaReWuYI2dw6dRkmDpGnaiSpTxm7tJKtsR26tWUvWjzF3bNxuPH4bR35Git0yTaxTrSok3W7KevHLLaDodEP3utHne20Jr5F8t/L6PUed3sNpHTdjdMspWHzcwX00sn+3tCBwLvldByBIgNmtgW2d9mf895sLX3A9L5w6COaAVLaCBt/C7iRknSJ5/fQR2ry1NeP97To5UQRTJZesC8vznOrPKLJkZh87PGNXR8nf+j1+6knzgIRXACVMZq/1t0j7DpKu/NNrjq+97qSDdOWfXnPNojJ4ZX1Ju4JKdp2lMycO1sLl9a6f0S5p+852X1mORlnCnIte14zxntYpFpdOHa6D+3f37EDaBTmDXB9O54H1MT9TJ92CZHZBaOtznKYyD+/TNSN7JkhAJExAwe/17vTebtdXLtc8c/pcPwvpuwUpb3gi894244v9fZXJyTmTh0oJ7RqckTKCj1bW392pnps1cZBmTRwUaC06p3t+VO04P1PWzJLalTlz6L4d9Zw5gC51HLOffOUgXWGpe+2mc3m1aZwG7B599T3bejbMNWIdjEhKWvTSv3XptOH62eNv+grsWB07so8eX/ueY0bQH195V3/83uF658NWJRJKa2faLQBut5GH8X3tNvdwu9bdMnd/99I7GeU1gkVe9xrzZ37QssM2q+rRVxt03MG1mndK+v3FLijXpVOlY2ZeqSMoBV/sbgQ/e+JNXTptuG547A3bkTNr4KJCHQEpu61d7TqHbqOkftI07US1G8qvVtS7prqGLVNC0tlHDNasiYMlpY+G2n2U1+KD5u/plk1lfC+7TAoro6I2Z5oY72meamVU3EbGS5ScbujWx70y+3K5ZpKfkXy77xF0rR3rezit4+a23pckTR3Rx9dOWl6CBH2zDTA7sWtgO03783NeRBWMMQs7/SBIuY3PcUpLNzeurezWkJCimRrtFlAxc6rfosiSSUg6b8pQHVDbrSNIJ+cUfet39NOhCRpA8fPbW88Z82+xfeenOuuB1WmZtHaSkj75tN13Z8zPlDwjqGQWZDTd/Blux9ValqDnop9rZtKw3hkZyT974k2tmF3n+P5uwbSgUy7szgPrY352lXQLko3ZN3O6u/l93Oq2UQN6aJ7DuRxmsXQ/AQU/17ufgGOQgaqoskKtnxtkIX2789suIClJv3vxHT300jue2bIvb/xQ21o/VY/qThk7f50zaahOHNXXNnhkZbeuo1c95zeg6iaqdpxdWS+Zur8O7t9d1Z0qtGrjVl336OtprzEyZ5y+R7ukHtWdfE/ncuN3wM4ctPGaJil59ynakkkd3K+7Y3/Ba+Dg9MMG6YrjD9TqjVv1zy0f6dZn1me8//ad7Tp+VHq/wG4B8OmH9tMjr2zOOJbWshtJEiP716QF/aycpkg/aMncNZw1cYjvQW/j340fZ07nlDoyxa5/7HXNnT5SK+ccqdUbt2YE5eyOg1Q+60lJBKXgk9ON4OB+3fXI9w53HDmzC1yEXc/BGG04flTwUdJsdoQzs6swnIQZuU1Kunt5vRauqLddIFLa1Wh3m2tv9z1XzK5zXN/KrqHkxJzpZq38zfP1lZDG7Js5opMvfjL7cjn6EOVIvt+1dtzWcXNa76tC0nenDNWdz23IevQxzM5udgHmCmW/xoC5ge3WYfXb+YhynbooguNuG0GY1+pwCiK4Na6lXWtI2I1UBgmg2I1mSruO56OvNmQ0/iX3xeyl9DUw/EyNse4+l5T0i6UbNOfYjoGShqZW3/WD9do2UvbdRjazPX+czhnjt1i5odH2+rYbcbVbg8+O3yl5XgGCIJ9hXYhYlqw+r4CNmbXO8RPMsVsnxavT69SBTMo+uJ5ttrbXvcUrSGZ3/p76xQGasF8vXxmydueynzotSCA96BS8MMEK90WKO7lmjIQRZiF9PwFJg3HfHt6na0ada9vZVebOX+bPsxs48FrX0aueyzbYF/b1dtecW1n37raHrn/sdcfPcSrHmEE9bO8N5raVeWMmp/L5GbCzXnMnj+6nP77yruM0Sac+hZl5WprfwI7da48f1dHXu/3Z9Z6/ldMC4A//bXPGZ9gVvz2p1ECsW3vKrh45c+KgjGMsdVwbsyYOyni9U11uXovRibkvZgTlGppa04KGtv3LZHmsJyURlIJPXguSO42cSekXctD1HKx+8NAratn5mWdKtpnbiG2Y3VDc6nTzaEvYkVvJfYHIxeeN1/ad7Y7v79U489pJL6mOhsc3vjRQv3txU9r3neOQ6WbmZxcvKbcLjrtl9kXZyDR/nt13iWok3xxQcgsYOS2kmJD92hDGmkj1jS36xdINGZ/rdm343ZHRa2c3o6FgXcR/2VvvZ73el5/gSZBdTvy8n5dsph9Yj7m1gWVuhBpZl8cdXOvaiHer2+3qTj9BND+vM76LXYCkQtIj5x2eMU3BmqWakPQdUwfJeny8pj9LHdl0xjS7IEFk67W97K33PV9rPn+C7jIYZlv1ykRCl0zbP6POs67BZw7+GBKSLpm6f8bvHDZAIKUPqDjdJ63HVVKoQJ5TW8Op7G7tBK9Ob22N/ZRhu9dFta6h273FLUhmTCOaNKy3VsyuS+0q+ruX3tFDq+wzbbyyOv3WaX4Wxw8zBc/4zkGDFV6LFBvZlBO/0DuSDXOyXUhf8s4ybUsmU4PExvGbNKy3bWc3KfcdR83HvLpThWu701pGp+eEGazL9vVu15xTWb0+x+3vXtO52iXd/sx6XT99pGv53GaE2F1zf3zl3VT/wG6apF2fwpjO7zW443RczOUKmjlncGq3huHVnrK7v1iPsSTNPna45zkZZm0tc7vabqdKu7dyGtwoRQSl4ItX5RJm4T5r5WH9DLtGlTHKYOzw59VJ9BqxDTJC49agCLt2kuuol+wXnrabT+71nnadTKPTevCA7raZFMcf3FfnH7mfY5qpHb+N01wvOO6W2eeUlhw2SOb1XYIEMpzSi/2OAK/9d5Pt+84cv69+88K/Ms4J85pIQRr0dt/ZaUdGPzu72U2ZCZOpFbSjb2RPpo0e53hUKuz0A+u1O3vacJ0zeWhap8HcCE1qV9aldQTVT+Nakm3d6TYaLznXuW5BSrvyWes5u5H+pKR7l2/UrAmDM9ZRu3Tq8LQ17pzajeaR56BBZPO17fe11sCan/ovmymmM8YNtM0YtJZ3S/MnaRnPSXUE8k88JHOjDruGvTk7z+ke5DWgYndcjX8H4bYzsN3v5NZO8Ntptk4ZdlpQ260udFpw26lOc+tQ2wXJrNOILp06PG1XUbu6OdsMKOPvRiaBNSDldyFi47v6DXZUqCOoGqSDbc1qSaojm7Jb9e6+d1F2+72yWUjfzGsHOesA1q3fOMRX/WcnioEYq2yzRoO83qku8NO28Pocr015jGvari/zu5c2ad9e1TpxVN9QG+I4TYPcvrNd44fu5Ti98DsTh+ie5W8HWjPP7bh4BSv9/FZem9ME5dWesp7TfpeZkZwzr/1yCyouXF6f090JiwFBKaR4deicGnPWlFMnfhrW1s9Y/a+ttqMMxg5/YRbINAQdoXEKIBkNtTBrJxk3HbsdbawLRIZdANepk2l0WhNy3u62tqYjDddN0HWb8rHguFdmn/VzwgbJov4udr+d3QLyTinQdgugStJvXvhXqMBEkOCy3cYGTovqmqfiStGsExHkN3Qb4crFLifma8Tv9Cfj+ZIyph4n9Xm6eqJjmqJTI1XKHEH127h2ej/JfjTeONZudW5bMqn/94//6Mo/v5bWcfIqn1uW6q6FSnetC9Ge/Dw4YF+MNH7WefPL3BExB2kMdoE1P3VGtlNM3YIXUsdv9mHLzoxsSrdr0HhPvxlJl0zd3zaIGTWnKXjWdoO5DH531/TiNGXYuJ4/bMncEdE4xn/++7sZm1NICj14Y7eunjng4nSNWLOq/dzfnM7PVzdvS1sn1MxucfxsNzAY0LNa350yRL9Y2jEV3bgfmteasTJfM+u3fJSxSLWUnk3pxs8AVUaGopQaaLDLTHRitM1adn7muqREWzIpJZ0zRe06vLnMYjdkG+zy+3qn6evmTKVsPsfP3+2m4Ukd51W/7lWu57xTne6Vdeh0TR53cB8tXPF2oDXzJPvphZJ8nSd+jpHdAuAT9+ul5f9szHi+VwCrIqGMZQzchJl14zeI5hR0dgwaHpEeNLRbCL2UEZSCJP8dOnPlErQj77dhbf6MMfvaL3ya9Nn5d8o8CbsjnDWi7ja33i+jQjTS6O1GzKIYSXLqZCbVUXF6rVVlx8+6Tdb1gaJaqNJrBNmuU1Tf2JL6u/l97BreTpkgZrlYPN1uXS67BeSt7+8WDAgbmAjyne12ZHRbVNc8FdfpOq3uVOHnkAUKDv79na2abdnG2Szqdcbs6km3AKDbQrhW5s6SW9ZlWzKpdz5szVhc1MzacPSaTu00ndTtdQl17Apn/ZN5hNeOV5aq3UKl7ZLtgqNGOZIKPjDhh9tGBW6BNa+RXb9B4yCdvYxdGRVsxx+3685cj618+4NUYD0XmbF238eOU7shyt01rcffzzFesf59/WLJrunTRjaHOaAfZsDDHCSzC7jYXSPm39zv/c1xMMVljRW7xfGj2MAg7TOS/teaWfbW+447VnplE0n+70FOWU5GcMAuM9GNNcvYqkLSgJ6ZnX7JvsPrtPFOrjrFuQiAWQeB3DKVvJajiMKsiYN1j80OpO1JaWvLzozfr0LSBy071NDUmqpP7H4jt6xDp3tGy842z0FCt2nwUQTM7Zivi0RC6t+jSps+3K4V6xszjs0jn+/et3rTVt33/MaM95q4X69Am/b4OQf9rAdnZ/axw20HKrx2qgwyQ6WUEJSC681Uso+Eh8kQCZqNYXyuUyaRn86/02eG3REuysWOreW87LgDAmdFeb2n305mUtLtp47WXnt29v3Zbus2mbdaTyp9faAodrXxExA1/1avbt7muFC4U8PbKRPELKodeqw3Rbs1mKzTDv0s4mv9TkEDE3bcvrPdxgaLVm2ybTC7palLHZ2WkxesjHzLcLdR5SA7SRncnmcNgBnfecXsOttppEEbPubOUm1NVca1Z2YOAvph95tI9luFu+00Z3ALsIVZ+0XatQ6D3Sh0RUKOu8MmtGv3wSgbfG73Ra+sXa86I+p7j925FnRwwm2h6ONHddRjYTLDovg+TvwGVaIIVvo5xpdM3d92c4p24wUeZffiFnBxWnPMq+Nkd65az0+v9ZPsNmUI8zv4/d2NQF+Xzrtl7DrntV6on+szyABVbU2Veu4ZzYCW13E230eNXb+2te5Uj+pOGfWf28Y71vuw10YWfmSToW5dz9L4t13byS1TyU8GXLZqa5yn015lyho2HkuqYydYu2Nid77bZR1KzrNb3AYJpfRgk3UafJCAedCAY21NR/bfolWb0tbGNIJ2Rn3wxnsfOV7zCUkr1jdmlNd63dstUO52DvpZJ/GYg/bRE6+9l8pMNU8D9DNwbq7rvGaolKpAQalPP/1URx99tO666y4NGzYsV2VCnjndTM2ZO3478l43VD8Na7sbldsOf5J75Rd1Y97acQ8z0uP0miCj3EEZlaBdgE+SNje5ByysnM6BfjVVrusD2WWcea39YBYkIFpbU6UtzZ+4LhTu1OF1er71/cMEWp1GJp0aAbMfXqtrvnKQvnzgPo4jV9YpM3a+/+Ar2ryttWNqR0hG8MNp7RTjv/WNLdrS/IlrY9+apj68T9e069xvBzbIjmBOZTFPn/DbWHZ7nlMAzPjO44fuFSjbzY75Oy5ataljzRaH54YJBtitHWG3VbjbTnPVnSr06KvvZaxjYUjIe7tjuzrDyFKV7BcqNRqE/bpX2U4BN3YfjJJbkCappG09Y5x31sZykDWEwlj9r62ZdZ6CDU441Z0/eKijnrnhCfvz0a2tEPZ+ajdlz47foEoUx9nufLAeY7fNKRI+A0JunOo9Y3dLpzXHpODBOuv56bZLl3VnMKcFzSW5TsMJUme2y76zH8UyD0EHqKIa0HK6Bs3ZQebBELc2ntvGO+Z7h9tGFkYA3s9gTpilD+x2nHvklc2pAIZsvvfi88Y7ZirlawFpr+m0kjI2nLA7Jk7TEa1ZhwbrNenUD3AKNtlN8fUTMI9qSYykpIrkrkEkSZow71nHesVuNz3rdS+5r5dpdw669ROsm6342W1eyl2SQzELFJTafffdtW7dOiUSCe8no2g4TZ25xzTv2E9H3u8N1a1h7XSjWjG7znGHPz+Vn5/GvLHwcSKRyBhJc+In1dk6ouQ3Mp8LRiVoXYNFkm54LNiokdM54LSGkHXtsG2tn2re58fhhifeUPfq3W1H4qyVeZCAqFdgwJxhYm4sBMnKCxtoNU/p8WoEJCVd8afXdOWfXtPsaZlBK2O64YCe1akpeq9u3paRJZJU+lpEYRjBD+Oct66DEWTevbXOaNnZFmhdG4OfzpNdJ9wsqY5Mv8OG9PTVWPbKMHUKgAXtqFjZTT2zGz21O/Z267t5dSDs6k4/HVXjdQ1NrR1rWDh/JV8LyztdZw1NrTp13AA9tOod2xFKu+mjUU/RNPjZzcsI2Js7UEb9J0U7JcKJUS9aVUgaM8h/9phbB8dtTS+nqblhOjJBdkAKGlTJltu28V7Bm9nHDlf3Ku9p22Z217NToOG2U0enAhRu39tuOrkfdnWyseOi3c5gdgua+zkf/NSZVn7as17LPNitsRM0gOf2fL/BWev7GLzqfjteU8CNheudOvSzH16bCqS6Xb9OQWQ/O/5a77kP/22z43c23nP7znbfu2NGzdyvOHFU31QA2G5XvmTq/6SX33xMoghmzhg3UF0672Y7WGP9fLspvhVSRjvfmiQQdPaNwSnoZgwiOS1DYqz/J9kPUhnlsAberNyyG+3udUlJC5fV6+D+3TVm3+BL3OQyEaEYBZ6+N3PmTN17772aN29eLsqDArC7OR6+314ZC8y5TdGIKt3dLeDgZ9cccwfdT0qxOYXTPO0lIWneKd7zkL1Snb0azGGyF/zwGm0/flRtRlCqXdJ9KzbqsuMO8PUZTueAn05gQ1Nr2k48biNx1src703ZLTPGLsPE6DR+d8pQ3fnchkA3/TCBVqcpPW5r4SSljqwvy+N20w2dskSk9LR1p3PFaSco6zlvXgcjyPQzuzojmwaXW3DQqRNu1ZZMatXGzOCVdfFfr0WLk0raXvNGdoLfjoodu6lnThkZbuvFZLMDZpARPq9MhqQyR6v9ZpLaBZ0v/XxXQvNrcnGvsmOX0ZWU0q6HiqR07UkH6co/vZZWJ0SxhpAfbvViUulTrf1w6+A41WN2U3P9Zk5YNwHwE5AKu45ktvyce3bnjDmo6vc6c7qe3QJjftlNifI6R4yFx+3WMXTaOSzMhihuwa9XN29LTU208tOedVrmwelYB818cHp+0LrZLsvYys991KnjbX69W32elHeGuVub2KuMQTOJze85fuheGbtjOq0zGhXr/cncr3CaRmeVUGYA/8yJg1O7t4a9n9m1052CTXZTfCU51mthZt8YvNqATgu8m9f/c2tD2QXezNw2nXG711kzMLPtm5arwEGpnTt3auHChXr66ac1duxYdenSJe3vN998c2SFQ/5Yp1uc9IuVGc+x7tIR9AbsZ+THq0KydkrCrgdk3ZbbKqmOxq5bh8Ar1Xl4n66h17jIht/RRbtMioUr3tasiYN8l8XpHPBqiHttIe21LbSfTqbbHHCnDJOkpF8+97YunTbccZ2NoIJO6TEaAXZr4RhltOvs2TUGx+zbw/Z3bk92BAKcOhtO55BXlppTNpJ1DRWn7YftOmhnThwkKXxWj1Mn3JrqL3WUb9wg56BqkIWh7Rp8j5x3eMa6D1bmjQ/sprzZTT1zqjed1ouJYtdIvyN8XpkM1kag3w5ZQ1OrbRbkDU+8kbFYsJ+dY6Ni/iy70fB2SZ982m4bjPazhlC25fbqVIYJhDkNRJjPP6epNcZn+cmAtZ4bZ00c7PhdrGtjhVlHMopzxC3Dz9zZ8dpa3qucYbaV9/Mdw9QVdteweRMDP4MPQTKinY7f+KF76cRRfTMWFbf7PL/tWa/jYfzPaQdOK+vvG7ZutssyNn9Xv20Yczb9Q6s2ZQQ+/rzmXc/3MHgFGs0SNgM2QdfNlJx3O5PSF/53W2c0Cnb3J2u/wtrWcWrvGQF8SWntj+9ksdGSU70gZQabnKb4Ol0vTpmHRkBKChZkts6KcVvgXXLeTMAoh1OmlPkYrNzQaDurZdKw3o7noNeAs7VveunU4a47g5ajwEGpdevW6dBDD5UkvfXWW2l/Y1pfcTNujis3NNoGks+aOMS2E+nnYgqyu1+QUe0w6wF5LXZsMDruQT9bcs64sOM0OmQdEY5ynn5tjf02tV7f2Y7dOeDVwHNrlPppjPqZUuB0YzQHBpw+6+B+3W0Xog7D6Tz5/oOvpBa+tS7kOGPcQB02uKe+4hAcNgfN3KYbjh+6l2PaenWnCseRHKdzyO13c8pGqkwkXHf9s7LuRnn38nrd8/l5ahwv8zRZr2m3Tp3w278xWi07P8uoa0YN6OHYWAuyMLR1Suil04a7BqSs17vTlDe7+sKp3nRqTIZdEzBMJ91atoSUWjfDroPst0PmtBaPUx1mrqeCZiKEWbC1tqZKdz23IeNvToFP83ExPzdMwM6Nnw0RwtwDvM4/uwCd+bO8AhV258Y9y+sdpwjf9vm6TdWdKtSysy21k5WVU/0RxbE2Hx+nbBHze0c98OG1rbyf7xi0rvBzDftp5wXNmnVqi9bWdCyebFfXe7Vn/U6HbEsm9eirDTru82yNbM4dP8fbrlxhd5q2e69lb72fCkgltGuKfkNTxzpxdoysUOsgjXkxe7eAeCKZPo3bz7qZxtpkf3zl3bTf1TzAbr32jf9+c+ELrudoFMF/r/uTNRnAulaj+TXWTNqkpHuXb0ytqxiGU71g95j5/Df+7XbNWX8nu7WeggSZgyzwbpTB6bqXdgXeDMZ5LtmvWWWcIytm17lmYhnfyU/f1M/OoOUmcFBqyZIluSgHYsTp5jbr84yFoOwaKXMeds5CCpKBVVsTbD0gr8WOzayZYXaf7ZbqbNfxsHJqHDltCx/lPP1ZEwdnzL12a/SF7Zg5/c0u5d5Ye8vu/DNvkSs5Tymw7txoF3Qwvs9f/p456mccA7fyBzkWtTVVunRqZmAoafqveSFH4/1GDeihGyzbOCdM39Po7HktPG1dYNM4Dk5bA7tNXxs/dC/XgI31VDey0ryyg+yYR9Wso43GNNmTR/fT4r9tdp1267Wui92iupOG9c4IStpNOUnKfmFo85RQqaNOsls3zeA3A8Rt+p9TvWl3HoeZJhnldD9JtvV7kA6wU7anV70dNBMhbACrS6dK2w7cJdP2zwh86vPvkUhmBqjDBOzcArXWutcq7DorXuef03SVVzdvSw1GnDdlqBYs3WAb5HWaovqlQT304sataY8ba2N5TTtzml4zaVjvrDMJnUSRpWjl53q2C7j4KYdTe9DpHPF7Dftp50UxTSnI55kFmQ4pSdc9+rquf+x1241Kgvy+Xr+lW0Az6E7Tdu9lPfeT2jVF3ykT2ljPZ9lb72cEjMyL2V86dbhzlomU0U63O4Z2v+PFx+yf8bu6Xfte52hUwX8/9yfzdWnuy9gdnyh247Syax/YPRYkg7m+sSWjDSVlrvUUJMgcdIF3gzXwZ0ydW3ze+LTprkl1tJH9rDflloll7j/4yYSTcjdVvxgFDkqZ/fvf/1YikVC/fv2iKg9iwM8oVhBOlYnb+kVuAQGzoOsB+Z2TbnT+g2Z1mBtPowb0sA1GSNL36oZq4n69bRtH1m3kzcV1qrzCzNMP8jtHOWpsMN8srKnUJ4/ulxr9Mm7s5jnbTp2Gba2fZqTb2mU8uWXLzfhif9dyhzkWI/vXuP7daTcw880vkVBa0MqpMWP3O5rT1t22Bq5MJLR1+86M8pnPoUnDemv+qaNUkUikyuO0+KR5Md0g/Fyn7cn0RU4l+2m3Xue5cRy9fle/ixY7TVFwunbtGuALl9eHmv7nt94MWsf76cB6BWqtZbN7TpBgWW1NleY5BG3djkGQwFc2ASynDJ6D+3WXZL8OjFOA2q3c5gwNowxO65gYzNfwv7e1RjZV2Wsgwu5+OO+xNzI6YubsDKnjd/jg4x22x3TVv7bq/Lqh+sWSDWmB/j+vedc1OOA2vea2b4wOlB0URNgsRTdh2mx+zyfjvc3HKqmOTr/dTmtBr2HjtzBPdbNeR9lMU7L7PC9+pkPatR/ak/YL/Af5fd0G7P7xblNa29BaLrsOuFuGoN13vPUbh9ieF8Z0civzej5umT/tyY7g1nlT0q9V8/u4tdOta4BZ7yXm+9DLGz90/f3czlGv3z7IAvRB7k/WgSxrXec07ezVzdvSpsbmgt/7oFcbKpt+Zbaba1mDlGdNHJwx3TXIelO1Nd4ZmH4z4aTol3EpVoGDUu3t7bruuut000036eOPP5Ykde3aVT/60Y90+eWXq6Iic0cVFJ8ot6qMav0iqzDrATmNchmjOP26V2V0/r3U1lTpsuMO0KyJgzKOl1MwYuJ+vW1vJH6mFgaZp+9V8bv9zuYR/1yNGhuvt6ZS//GVd7X4vPF658NWff+hVzIaYU4NJ/Ni9cZzV8yuS43EG9yy5X734jt66KV3HNciC3Msgq6rY2bc/Jx47YZnfh+3jIkKSTPG9dedSzOnHF0ybX/XwE0Ui+mahdlRydCelFZv3Kqee+6aCue04K4hqiknkntAza7h4RS0/87EIbp3Rb1tlp9VmKkG2S5WbqyL1nPPlsh2FA3auXYL2jrx6pCYp1AGyTz1s8i/9Tq3WwfGKUDtlaFhBOu91jGxu4ajmqrsxe5+mJTN+njalZ1h7kzYaU9Kw/t0S1tnLynv4IDb9BolM9eEC5tBZhXF7ll2grbZ/JxPxjU8aVjvjOPrtNNa0GvYej5aM42Syn6akpOwO/xaj4eZ3QL/fn5fpzXGzAN2dqz1kV0HPMi0TLtz31gPyFoEu8xdo53htJj9L5ZuSAuYJ5XZTg2anWdwG5z1u1GT2yL8QRf893t/cmrDW5cG2Lb904yg/s8efzPQrtlh+J1S6tWGyqZfGSbwbggy8Oe13pTTwLHTdzK3u3ORoVxqAgelLr/88tTuexMmTFAymdTzzz+vq6++Wp988ol++tOf5qKcKAC3UaWg05eiWr/ILMx6QHYd8bNyOAoXpAHqd2qhnwVBpV0p1X5+H+tzvEb8o4zqO/2O23e269V/b8to+Dk1nORQTuuOH26L4xqcgk1BsyzM14jfdXWC8NoNz4s10+93L71j+7x+Nd6LY0eZXWl3vILEp8zrdUnpU1/tgsFRTjlxC6jZXftOdcSsiYNsA91W2WQx+s0csCtjQspIXZd2Tc/u0nk32zW+vARtuHoFbe2eb3eumjsc5vPGyqn+dttYod3hOg+aVeLUqDUH650CLW6bSKyYXZfz0XYpWLDZCHp6bRbi1JHwCg64Ta8ZM8h+Xbko7ndR15XW9w4ymOZ1PrktQm8OJrpl7Dhdw04ZLV7BRPNgWcvOttR/g679k80Ov24DD24bTIQpi3nAzunzrLuFZTMt0+7ct1sPSHLPhPZaU0ey30lWcs/Oc7q3uQ3OSv4Xtnf67Z3W4PQakPRzf3K6vqxLA6zc0Jjx2nxk2ES5MUGQOsoqbFAryMCflL7Qu7Exj1vmoZ/vZN2N1LwzaJT3gGIXOCj161//WgsXLtSJJ56YemzUqFHq16+fzjvvPIJSJcjauQ7TCQq6fpEfbhWlWyURZRaYlyANUKeGjlcAw+k4GAGpoFkUfkb8K5S5TW1Ybo2Ae2waQuZOg916XtZyWnf8sBshsWN3Q/XbibQG9WZ/vk29n3V1gohqKoj5GNn5wUOv2Abz/CymG5bdOkTWabInje6bMYVPSl+vy+DWiAwz5cSJU2fPbXFdr+mFTnKxPo2f72QsauuYvaLMLZKDfl4u62W7c8u8uKnTpeBWfzudQ24ZemEzwx59tUHXPfp62t+MYL3bOia5mDoWhN155FR/VyTkOGpt/o5JSZu32k9FdgsO1Na4T6/JZTsh2/eOaudIr/PJbRF6K7uMHT+BdCu3YKLb64LUNX520HO7Lt0y7t02mAhaFq9p7HaZSkGCA26BH2v9aNd+d8uE9nOtO2WFSvbZeW73Nq9AoZ+F7e3K7bUGp58BSS9+lwaIIssyis1KotiYIKwwbYOgA3/Wujnb5UvsXm+3pAZCBKU+/PBDDR8+POPx4cOH68MPP4ykUIgPr9TqIKMF2Y4OWivToO9pXfzVz8iw9TNzOU3Gcbe47x2uvbvt4Zoe6rWFqvn389qC1GvEX+poTJz0i5Vpu6CFteyt99NHzz5vaLXsbHPdBXLGuIHq0nm3jB2dzOU9a6JNhp6k08YNTO0uUyHp1C8N0EMvveN5Q/VzztkF9eZ+voDiOZOGhj6P7ETREPC7fpNdMM9rMd1sWd/Pbpps3fC9Hc8BK6dGZNTZC9a1BLx2HoxyBDBXQQZzGe12UrOTiyBZVNeO+dxyWhPNzCvz1Okc8lrkP0xm2HEH1+r6x1637cR4rWOSj46DG+v3/fOad23XXbx02nCN2Tdzs5AKSUlTVCqpjuxQu2n7XsEBr+k1uQyOhnnvhqZW3f7MP/XgS+/42vjE79o3TueTeUqrdcfXpDJ3WvNzHvnJaLELJkpyfZ1dXRN2ep7kfl3aBVysGfd+f1+3sngFA787ZWjGb+/UjrRuFCN5B37MxzHMovNea+q4nTNB721O2bzXnnSQjjpgn0DXmtPub36mfQed4if5b39k206JcrOSXLehohR04M967mcz8OcnAI5dAgelRo0apTvuuEO33XZb2uN33HGHRo0aFVnBUHh2F1M2izjajb6YF7d041SZ+m3Q+1n81eszTx7dT4+8sjln02S8OjVur/ezhWp70t8WpE6BjrtnHqqzHlid1oAxdkELu4aM3ZRF89bAdo0r8y6QY/bt4bol+JhBPTJG+BJS2nbHl36exTSqf3dfN1Svc84pyHPD42/oxFF9QzVanETREHBqzFm/glO6c75vqtZrya7j6sStERx1ZoS1nB27PW623RHN7vl+5Hp00m4gwKmB7iTKIFkuNlyQ/K355mcqdJhzKOqRa6/NEaLoOGQbGDSf69bdQSv0eZ08aaikzMVxnbYWd5q273Vd1dYEm/5ZKItWbdKlD69Ne8ypgxT0OnE6L6z3qkunDtfB/btrUK/qjJ3W7M4ju/PEaQc3adeAlF0w0U/g2O8Oan7rTbdzJ6pr3SmIVN2pIuN3sfrl0rf134ft6xowM+7ldpmrfgI/2S46bz6GQeoep+Oy/v2P9EHLDl+7iiYlXfmn17R7ZUXW2bp+rhHrtO8gAQy/51OQ8866RmK2WdVedWk+Z6EEVaiBv0JnJxebRDJpt1Sfs+eee07HHXecBg4cqPHjxyuRSGjlypV655139Nhjj+mII47IVVljpbm5WTU1NWpqalK3bt0KXZycWLmhUafd82LG43ap1Stm1wW6wII0mhqaWtOmVgT9zIamVh0+91nbKQ3Pzz7S9j3sPtMqzPf2o6GpNZK0/vrGFtvfz+BW/kWrNmXcfAf0rHZ8v7DHwukce/DswzR+6F625bCeJ3c9tyFjpN1cHvN7OI3wGs91O/Z+s+aczjdJuuMbo223kM32PMrmnJEyf+9Lpu6flhFpLqfkf9phVFktfsrvNJXTaJQ7nT/54BYUz/YYuV0j2by3Vx1t/dwZ4/rbrknmVs8aZTRnsDplOWR7H7B7T+v3dVvzze28CTNFur7Re4F4r/f1uu7d6qiw9YW1kxpFtqxdmawdKqepllLu7sVx4XZPkXbdL43nhj0+5t9A8j7ObueRXf0hyXUjlwpJz8/Jvk3mp+x+2hZRcqtPzWUxmJ/T0NRqO8VSSv/tzRqaWh23rDe3d9yOUxRtb2v9Y3fOONVT1jo5bfBS9gPLf39na9qOpkHL7Oc7uV0jdpx+I6/PybZdYN1tzm5dsDBlKydRXAN+6tF8tJMLyW/MJHCm1OTJk/XWW2/pF7/4hd544w0lk0lNnz5d5513nvr2Db71N+LLaTTJnFpdoY7dvoJcSEHTIaOIVNvdM9qTzgut+5nSZN59KmhnxO35YVM67abquY3+ux1Dv+nLft7LjdeIpZ/RDetIu3UEzmvakdtWw4a7lm3QvMffUNJH1lxtTZVmT8vc+rxC0tbtO3MyahL2nDHYHefu1bt7pjt7dfJzkdXiVH6nqZy3fyN9sVA/omwkuG09v631U9+71jkFRScN622bIRI08G9976A76UjSQ6veyTi/L5023PEYOgXrJGWUfUDP6qyuHa/jYfd93HYlDbvGotO6ONZj7Od9jeu+oak1I+vYa+p2mAEP68LUUWTLWr+LXdmj3Fq8GDm1YaRd64WZnxv2OjH+Xt/Yog9bvO9VTueRXf0x5+G1HYFel89vl/Toqw06ziYr0StryO8Oak71Vy7PH6/6dMa4gRrep2taMMX6HLcplnZqa6rUc0/3c8Ep+8d8HoQ9l5yuYes543atm7M+v/9g+r3duIda+w12O5pGmZliLr+f7L0wmcvZtp387jbHjm/evK6RbF+fz3ZyMQgUlPr000919NFH66677mJB8zLgdDEZF8y8zztTNzzR0fke2b/G104oUcwVD1KZuu2y4/Qefhb2TGjX7lNBOyNRVz52NyHrOhtWXsfQ2nioranSpdOGa95jb2Qcy7A3Nz8Vvp8OlNeigebOm9e5ZO10WjOx2pNKW1zbrsNuDZRJHQ3uK//0Wsa5GJeGgfU4ezXa3c7nfC3AbWY3jc9usVAvfq5TPxks5oxFp6C4Eeg0/u10jIJOJfZz/J2ydS6dOlwVFZl1n5+ddKzrrJinYdkdI7tgnXmreXPZF583PvR9wO/5aFfnmQVdY9Er2GdlHGPJ/1QLu/N10rDeoaZuW4+Z3TopdqK8vnO9tXgxcmrDSJlB37X/bsp4jt/rxJoFF+ReZa3zMoKucvgCFtc9+rquf+x123PUbq0+uzX7opieFyU/bV6vYEqYzrGf4+B2LYVtezsFJYf36Zq2xp6fa90IrvkdWM7XgttOn5WQUvevsGvn+q37g6yZFpflF7JRqIyibO83Tq8vRDs57gIFpXbffXetW7dOiUTC+8koCU4ZMzdYOlPWrBC3tP6gNw27m/ElU/dXfWNL6u9uamvcd9nx85kVkqaO6KMnX/uP7e5TXpWJ00067LbpVk6NHvM6G9luQbpo1aaO3/3zf1unRbm9l9vNJKoOhp8GZtBRi0unDtc8m0V4rew67OdMGqrDBvdMG/1MquO4uW0THydBRsK9thAPOloZtAGS7YiW0/eyXqd+prVZzyG7DmVC8jWia1cmr6Co1/F3y9axW3Ra8tewnzSst+afOkoViUTGotFWTsG6pOyPy/ad7aF/36jOxyBrLNqtHWeX7WVmHGO7dXf8nhuXLV6nW79xSFYBJGtwQvKOJ0SVjZCPrcXDClInRbFJisGuDSNJcyxB34am1tQgodkl07yz2e026PB7r/KTpV0hZeymaO7Am3kFd7O9z+ebnzZvtgEkO36Pg9MxDXscnYIiJy1YqXkB17SSgg0s5/O3d/qsbNqzfo9JmDXTnHabKwaFzijK9n5j93rWm8oUePrezJkzde+992revHm5KA9iyHox+ZnalpRzWn+Ym4b5Zvzq5m2+p7xYX++0y47Ta7a1fprKCHvitfdSC316TQOzcrpJZ7Ntuplbg8b4/cYP3Sv0FqR2I/wJSXecNtrzWAaZhpIPQUYtzEE4N04ddrvRz6Sk208NPqUsTrxuptmOVoZtgGQb4PS6Tv1kxthlLM4+Nj3DMJGQZlveS7I/Rn6nEpvrHrfj75WtY8dPHR30N3PqaNh1VI2yjx+6V6jfN1e7VDptX1/dqcI2WGSX7WV+nbFwrnVBa6fyOl2HSmZuEmF9jtO9yi444UdU2Qj5zHQIIsj5HSaz0Stg5acN41RXHNyvu+f3s3utESBOqGOZBqdd/uzum6d+aaAeemlTWkBLkm0H3m69JLfgbj7uBVHyM7CabQDJSa6yPNzYZetJHeeS+X4ZJKMtyMByPn97p88K+5l+jomfHd2C7DaXD9kE5f1mFBXb2kxxvdcVUuCg1M6dO7Vw4UI9/fTTGjt2rLp06ZL295tvvjmywiGe/ExtMzhVHmFuGsZzvrnwhVDpjrU19rvsuC0Ga80I+9kTb2rF7DoN6lUdqDJxO2ZRpGzmqkFjcOqQ9ezS2TNDKo7pqX5HLew6nZJ07Mg+enLdfzyDqk43naBTyqwKffP1uplmM1qZ7TmTTYDT6zr12n3ULWNx5ZwjMzqUbut2+SmTwVr3uB1/P+tgmB0/slaXH39A5Ne5U/bH2UcM1l5dOztmdYb5faMYPXfaDeq8KUN159K30963ZWeb7Xlgl+11ydT9UzuaSdLhc5+1/Xy7bBe3+sVr/Z0gW7FbVUj6xpcG6kFL0CGKusj6WxnrVkr+d+uNWtApNUEyG72mQVszhNx2Csymk+NWzyTV0fY58ZC+vu+bv3txU9pubcZzF583PmO6nd16SZK0Yv37unPphoLcC5yEvff6GVjNVTAlF1keTpyy9QxhpyQGHVjOx0Cn+VwIsli42znk55j4ybCJU1A2SB1nJ+gukcWyNlM+s/qKReCg1Lp163TooYdKkt566620vzGtrzxYLyQvTqOyYW4aUac7ulVkbp81qFe1zpw4WPeuqHdslFsrW7djFkXKZi5vQmEbu8WUnup3Yf9Lpw3XOZOH+t6pL+qbTiFvvkG+l7kRqUTHmk9RNUByxes6dcqMMa4Dr4xFa4fSzzVr13A5aXRf/fGVd13PKaf3DjKoIEmPrWvQ5ccf4PqcsL+Z+Rz569sf6MGXNunuzxdkNW8/H8Xvnm39aHdutEtasHRDRlnd1q5zy/ZauaHRMTPJLtvFrVFr7QT7nbrtdH5Yp2vPGDdQ5x+5X07uN9Ys5XmPv9Gx/pqirfP8BhiCnN9BMhsl57XDwmQIZdPJCdtG8Qpm3bt8o/bas3NGEMbcia+tqdKlUzM3B7lzyQbXQYB8y+bea5xrXTpVpq81aZPlEre2URBe5382UxK9grL5FPZc8PM6r2MS1Zpp+RjcdAvo+63jvL5vXAe//bCuk9eys00NTa2xL3euBApKtbW16eqrr9bIkSPVs2fPXJUJRcCpwWsnynTEKNMdvSoyp896dfO2VLaWeSTQXIk43XiMDph1m94KSdWdKgJ/B6tcNWjCNnazWXQ135y+44xxA22nPToda7vf3m53tDDC3HyjaniE+V7mRod5bZpsGiC55nadWoOUdlk8Qa8T42/WNfLMv5tdI/XiY/b3FRS1juDaldHI1vnLq5v1uxffSXt+e9J5l1Ljsz74eEfo36y2pkpjBintWJuzUqOsz7KpHxuaWjWgZ7XunnmoznpgtW0Grd/zwKkcYTblcOvAGJ8TZOq2U4DAbrq20/fIts6xZimbj0fYDkc2uyb66RQZ7+036Pvq5m1KKmkb7Fq9cavvNSit3yub4KtX3Wd3DvoJZvnZ0GFk/5qM19oNAlTI+VoIw++56rTeoJ/z0LpGm/UomYOUxTT9yO7YuZ3/2WTwFzo73FqWMEGQIK9zOyZRZNjka3DTKaBvV8e5HYtsd4mM0/ljVVtTFXqacqkJFJSqrKzUMccco9dff52gVJmxu6DtGrzVnSr06KvvaeGKtyNP6zc+M6rME6+KzKnzZh7lMkYCjfR0yfvGc/yoKrXs/CxjxP3kBStjXREFbexms+hqobitD+CnzE6//YrZdYHSu534XQDZEFXDI8z3clubxq2TFeU1HpbddeoVpDQEvU7sfiNJtr+bXcDBz/v5HYkd1KtaD730ju/gUsZuXZ93IIP+ZnHPqPTbqTSXNewUdae1UyTn6Wt+6qcgATmnAIHXdG0pmjrHK9si6LkRdNdEK6e2QH1ji/7893czsoCszz36oH30+Lr30t7zZ4+/6bijpHVBcMl+DUqnY51N8NWo+zZ/Hhj0047zGnDzc227ZSqb1+NLqmOwI4p2UpBz1Wma4n0rNuqy45wzSf2u0fabFzbqiXXvpZVl0rDese1Iu517blOUw3yPuE3NCnu/ivI+l03wOZ+ZRU7XtV0d5xaczWaXyLidP1bFnOkVtcDT90aOHKm3335bgwcP9n4yioZbFDnoQtWjBvQIvcNDkAU/s808CbPbiZ+bit/53sP7dE3bma0YKqIgjV2njkVlIhHr9NRsGvT3fT6d0yyqzvWiVZs02+cCyFK0Nzq3c9r4e5dOlWrZ2Za6dr06lm4L/cdlPYSwQUq3v5vrOClz+s6ch9emNdiC/G7ZjsQGCQjadbYqkv42P7AqdHacGz+dSrcskrCdDfPaKcveel8T5j2bFlQZ2b8mZ51Vp/WzwmwFH6bO8co2CnJuBN010U+7w7wmkJk5WG/OIq1vbMkISrntKDlm3x6ea1AO79M1Zx0ZY6ddIwjrtMi5mVMg/5Kp+/va0MGp7pk0rHdG1lwU3zPoueq0ePfCFW9r1sRBjmXxM51Tkh5bu+v8aE+q417/eZA/bh1pr2MX5f07jh32sDMAor7PhW2rZhMci2pX5M1bWzOea52JYj3vnb6vW7sljuePVdwH5fIpcFDqpz/9qS6++GJde+21GjNmTMZC5926dYuscMgPr0XowlzQYSpLu9FMp4Z3NoED83uEWRw8iq19Jfud2UqpInLqWFz36Ou6/rHXY9XIikJDU6vuWV6f8bjb1Bs/72kEfOYsXms7rccpaBDljc7PVFZzmYxRXj/TWJzqlCiu8ShEWQ5rHXfWxMG2o+/WH9rv7+b1mwfd5ctYC8zvZ/nNprGKQ3acE6dOpXFu56KsRgdfsr8HG1PrctVZNX4Pc8aWnwyVqOoc6/mQkFId9Ciy8LzWhnMrl6SMOs/M+L7jh+7lq91gt8ZYQ1Nr2pqVdp/xzOtbctKRsQvCWhc5twbW3aYPLnvr/bTjnHC5Z1lfK0n/9+q7OfmeQc5Vt8W725Pu05udArzWoKhVMvV/sutI52Lakp9jF9V9M24d9mxmABTqPmc9B8IGx6LaFVmSJszL3NDj3ClDXNdbC/IZQZIECi3Og3L5FjgoNXXqVEnSiSeemLaweTKZVCKRUFtbW3SlQ855BZ3ydUEXouEtBc/I8HNTqa2p0qXThnumvpd6RWQ9VmZxHK3IVn1ji20WxVkTh2Sdsm43ZUiSbjt1tI4f1df29UHOL6+Gq1NKvle2gFPH0ipujYRcsKvjFn6+qLe1w2JNbc9mBy3jtUEalH7WN4i6/opLdpyV0/e020UsF9wyLXJZj04a1jstcOMnQyXKc8KuIxPm3PCzgUWQzmGQhZwbmlr18sYPlUgkdOm04b52lLTW+6d9aUDGlNqEpNufXe/62WF5ZftayydlrhNofB+jzjO/XSLZcW45MV5r/pxcfM8g56rbb+5VFsd75xOZ9043Ye6R1t/q7CMGa9bEwVnXFflsu8atnex0LthtRGEn3/c5p/t+0OBYlLsiO+0A3KO6U1b9TbtAaNzOHztxHpTLt8BBqSVLluSiHCgQr6BTri9oozP8YcvOgjS8peAjOl43Fb+p714pp3FdS0DyP/pmHKtHX23QdY++nva3UgtEOI2Izpo4KPB7+Z0yNKBnles6M35udH6DFX6mshqM39auYxlkEd1S4pSx8Z2JQ3Tvivq030hSJDtomd8vmy3tzc/N5c6SccmOM3M6pqMG2GeQRc1rKluu6tEwA1JBsmL8sJ4PYd7H6ffzWhvOid+FnI3p1sbTEpJmT3PfUdKu3l/00r/TAloVnz9u/XijrSE5rz3mxSvb12udQGud4lTneZ2v1s8xi6rDFqQj6PSbu2Uqm9m1GbtX756xo+ojf9vsmEEV9B5p91vdvbxeC1fUZz3Im89OdNw67FH0ifJ1n3O7lwcNjuUj837coMxpy9m2DeN2/jiJ66BcvgUOSk2ePDkX5UCBeFWwubygMxbJlfMikHELYDjdVPykvpvZVURxX5TPqXxOgaramiodd3Ctrn/s9ZIORER5rfiZMnTS6L46ecFK1/PE60YXdPTLet47XbMViY7dJI3OkXkxdKcFxONybeeKU8dmr66dbHcxzHYHLfNr7UYmg2xpbzzXLoMqqp0l46yQDUa3jFMpPhkKYbJi8iWbDSys9zWnzBdzsKmhqTUtICV11JM3PPGGnp99pONnOl17B/frnrrOPmjZofN/90rGa5OS5j3+RsdOdwrXdvDK9nXKcDCX1VynhO3AO93/rjjuAB17cK2k8IE3M7/XtfU3r5B0ls2uy26s55rdZ88cv2/aOqMGv8EvM6djGNUgbz63so9Th71YghySv82csgnGh733uA305OLYxun8cRPHQbl8CxyUkqTly5frrrvu0ttvv60//OEP6tevn37zm99o8ODBmjhxYtRlRA75qWBzcUHbBW8Sylx3wVAsAYwwownmiijui/I5lW9b66dp2WGzpw3XOZOHpl5XTDfybOR6EX5jylB1p4pUQEpKX/zWvNi45H6jy2b0q7amSrOnZW4dL0lTR/RxDZgVSyMhSrU1Vbp0aubx+tnjb+rEUX0zdjHMpoFifW2QBqXTc6s7VeR0Z0kvuc4edXp/8+P5+J52rItsh5l2FlTQOjtsVky+hLmenAZgvOovp+BOezL4+kPGdWqUv6Gp1bGd5JW55MUr2zfoAvRh7/tOx+HYg2vTAuNuU9L81hd+z4tc3LOsnz1qQA/NyzL4ZXD7raIa5K2tyd9W9nHqsBdL+yUfgaSw393pGObq2Mbp/IGzwEGphx9+WKeffrq++c1v6m9/+5t27NghSfroo490/fXX67HHHou8kMgtP5VA1Be0XQM2Ken7dftp/z5d9e9traEb3oWc+pbtTSDui/I5lW+eZXecuY+/ISWkcybtCkwVy43cKsxuI1E09tymDDllvhijrH4bh9mer+dMHioltGtLdEnfnTJUdz63wTOwWoqNBK9zZWT/mozH8nF9B2lQOj23ZWdbweqmXGeP3rVsQ6oOM79/nLJWjetl/NC9Qk07CyNIne0U0KjuVJGz8uWS1wCRW/01uFcX2yxSr00v/FynXplzZkGvT6/Pt/7dzwL0Ye77TuWQ5GtKmtfmPWHbh/m4Z0XVTjKO4ZyH12ZMCYxqkDfug6i5VAztl3wFkrIpXz76mygeiWTS465mMXr0aF144YWaOXOmunbtqr///e8aMmSI1qxZo6lTp+q9997zfpMS0NzcrJqaGjU1NbHjYAgNTa2pLa6tzLt3Ban8Gppa9asV9akda8J2IrINai1atcl2/Qq/n209LpWJhFbMrrMtS9iyZvM6a/mcdpKpSMh1qkIY+Q442jVuJw3rnbcyNDS12l4DbtePwe28McvmfLUrZ31ji06758WM5zx49mEZmSZxXzstCD8BjKDXd9Sczic/zy1U2XP9uXc9tyEje83ISjRnI0b9uaXIXJcYCh3MC2vlhkbf9ZidjDWlEtI8n8fBz3Xa0NRquz6fWdjz1evzzX+Xwi1AH6YcTr+JtOu7SnKsL/KV1RMnDU2tum/FRi1c8XZq6n+Ye7ydbK+RbPltP/h5Xim1RayC3PfjpJR/E6tS/65+YyaBM6XefPNNTZo0KePxbt26adu2bUHfrmDefPNNzZgxI+3fDz74oE466aTCFaqMuI32hZkWYm0Amt8nyKhNFCPj2YwmBBnZCFvWbL6jXfkumbp/ah0LM6epCmEr33xnLdiNAs5+eK0Sn2cD5KMMbiNJ1jUurIFBv6PkUYx+WcvpJ/sqTlko2fI7Yhz1yKW1DF7XVZARSOtzc1l2N7nMHm1oatU8m+mnbcmkVm3cGuus1TiaMW6ghvfpmrYuTrFmT2SbRWrUq6s3blUiIR26b4/Q157Tc6zr8/nJXIri8+3qhlywfo6fKWlJJW2v29Ubt5ZlVk9tTZUuO+4AzZo4KPLAhNPv8ermbTkPSvltP/h5Xim1RewUY+ZRqf8mZuX0Xb0EDkrV1tZq/fr1GjRoUNrjK1as0JAhQ6IqV87tv//+WrNmjSTp448/1qBBg/Rf//VfhS1UmTEabQ++uEm3WbY3DtL4t1tgNZv3iaLREuYmYHQqJw3r7bl4cNiyRvEdnYIYdtkGboEIu7WngpR7zsNrc9qgdJpiGpfOlnWhUbusDr8dKLfzNcz0Ra/gRaml/a/+l/8ARi6msearUVOIKbi53AHWae2fhKTOu1fEfivpOGrZ2ZaxUHMxBvOiCMIagaNcstvhtBizIvzwOyXN7rqVTfAk6vMyztkOuQhM1Na4r5OYq2Pgtz3op51Ram0RN3E+P83K7Tcpl+/qR+Cg1DnnnKMf/vCH+tWvfqVEIqF3331Xf/3rX3XxxRfryiuvzEUZc+7Pf/6zjjrqKHXp0qXQRYm1XFRoy956X7dbAlKS9y4/5nK4bU1fLOs5Be1Uhi1rVN/R2sCxri3kJxDhtPaU33K3S7pvxUZddtwBvssdhNeirlLhO1vm3yFIB8rvtRw22OFnIeBiy0JxOmZGlqaVW90TZQch342afI+6Bg0OBLlPuV3jV/zxtY5dYSPIPCknuQwi5luQIGwhO3z5ylyKA+M3sZuS5nYvHLNv9NvNS7t+97Wbm3atrZiDgYG4BhQKsU6i3/agn3ZGMbZF7HidH8WUjVMqv4kf5fRd/QgclLrkkkvU1NSkuro6ffLJJ5o0aZI6d+6siy++WOeff35kBVu2bJluvPFGrV69Wg0NDXrkkUcyptYtWLBAN954oxoaGnTQQQdp/vz5OuKIIwJ/1u9//3vNnDkzopKXplxUaE4ZTm7b3zqt8WPXsQi6jW6hGtNhOpVhy5rL73jOpKGuC/A6BQ9vePwNz1E1p4VjF654W7MmDspZ59s6RS5pKUOFpA9adgTaCjlXjUu/HSi/13K2wQ634EWxdVydjlmYOixqfho12Z5zhe4QRX1uG6zXeOLzx5Om/1YkpTtOGx1o+lU5K9Q0z1zxE4Qtpg5fKXCakuaVbR71eWn+3c2s98oo1imN6/lViHu53/agn7IVW1vEjtf5UWzZOKXwm/hVTt/Vj1Dbovz0pz9VY2OjXnrpJb3wwgt6//33de2110ZasJaWFo0aNUp33HGH7d8XLVqkCy64QJdffrleeeUVHXHEEZo2bZo2bdqUes6YMWM0YsSIjP+9++67qec0Nzfr+eef17HHHhtp+UuJU4XW0NSa1fs6BSluO3V0oE6y1FEJVyY6uhQVkr4zabCen31koBu30Zg23icOa6Y4qa2p0qXThqvi816U37Lm+jvW1nTsDmX3fkZDwspYe8rrfc8+YnCo12ZjxriBWjG7Tg+efZien3Ok5p2y69gZjaLzf/eKJsx7VotWbVJDU6tWbmh0vDYWrdqkCfOe1Wn3vJh6TZTcjr8U7FoOc14GKWchrrUw3I5Z0DosF4xGjZm5UZPtOZfrc9avKM9tM/M1fvtpozPXxpPUs0vnVOfS7fpGB/MxXTG7LjYd6FzIVfsoTDn8npu5PI/zeY2Y6wRrPbXsrfcz6osoz0vr725l3CuzrT/jcn45KcS93G970E/ZiqktYsfP+ZHLtlwuRPWbFMP9utjPv6gFzpQyVFdXa+zYsVGWJc20adM0bdo0x7/ffPPNOvPMM3XWWWdJkubPn68nn3xSd955p+bOnStJWr16tefn/OlPf9IxxxyjPfbYw/V5O3bs0I4dO1L/bm5u9vM1SkKu0gudIsRjBvUIXI6o1jopljVTFq3alEoVT0i6ZOr+vhtYdt8xH5kQtTVVmj0tc/0Bv6MCsyYO1sLPd1b089qovpN5pNy8eK1516P2pPci6HEYrQpyLTtt8R7VCE4hrjVDkHPD7ZgFrcNywS0zJdtzLg7nrF/Z3KeMa7yhqdWxLo5ztkIc5XuaZ6EEPe9yca8Ncm7m8jwu1DUSpJ6K6rx0WzZC6qg3qjtVZF1/FsP0nkLcy/22B/2Uzes5DU2tennjh0okEhoTIGM2H+1qP+dHMWbjZHtOxfF+7XQ+FLItHDehMqUKbefOnVq9erWOPvrotMePPvporVy5MtB7/f73v0/bhc/J3LlzVVNTk/rfgAEDAn1OMfMaiQ8raITYqxxeI+lByhXF+wT5vCDHwW5tpp898Wag0QC3EcZcZkKcM3mo5hwbPMPLKLPf4xTF6KTTCEttTZV67tkpoyGQlDIan35Gq1Zv3Jq30Zwg17JxvM1PT6pjHbio5Ptak4KfG27HLC6jXE4ZANmOkBbTCGsU9ymn31NSrLMVUDhBzrtc3GuDZNLkMuumkBk9hain7H53g1FvtOxsy7pcuWp/R51FEud2s5+yOT1n0apNOnzus/r+g2t0/u9e0eFz/V23+WpX+zk/4tJOCSrsORXH7EKv86EQbeE4Cp0pVUiNjY1qa2vTPvvsk/b4Pvvso/fee8/3+zQ1Nemll17Sww8/7PncOXPm6KKLLkr9u7m5uWwCU7lcIyJIhLjU1qowC3IcgoyceY3UFCITwrr2lCSt3NDoazTJz3HK9ju5jbAYx7NLp8rAi6DbjVYlpFTGVT5Gc4JeQ5OG9U4t9Cx1BKXiminjR5hzw+uYxWWUyy4DINsR0mIaYY3q/mD3e67c0Bj7bIVyUOi1zez4Pe9yda8N0h7IZdZNITN6ClFP2f3ul0zdXwf3756WgZ5tuXLR7o1jFkkY1h2IW3a2BVrf00tDU6tmP5y+ZmRSHQMUbtdt0Gs9m3rN7/kRl3ZKPsQtu7CYMs4LrSiDUoZEIj08nEwmMx5zU1NTo//85z++ntu5c2d17tw5UPlKSS4rtCDp1IWqWPM1vc3Pe6/9d1PGY3YNHT8Nj0JV3sZ3DdM48jpO2Xwnt5vHsrfeTyvryaP76Y+vvOu4CLrTaJV14XTzZ815eK2G9+mqUQNyN/0rVwHQYhD2+3gds7hOVcq2Q1NsAwFR3R+sv2eYTm8cAyjFLM4daT/nXb6XQbA7N3MZvClkALtQ9ZSf+0KuAuV2/NQ5pdZBrq2pymibRVU31De2ZKwxKO1au8rpeAW51qOo1/yeH3Ftp0R9r4zbYFqptaNzKXBQqqWlRV26dMlFWXzr1auXKisrM7KitmzZkpE9heiEqdBy0TDPd8VqvmkkJM2eNlznTB6at883a2hq1Q1PvJHx+CXT9k87Jn4bHoWsvO3KOOfhterSebdA8/atsvlOblPsrGX94yvvavF547V9Z7te3bxN8x7b9bskHHZfMzcePmjZofN/90ra39slnbRgpeblIWPKz/GN2809W9l8H7tjVgyBh2wDNX6zE70yMvO1+18u7g9BO5dxDqDEid/fNa4daWv53critEZf0N1brYKcm7kM3hQ6gF2oAUuv3z1XgXKDcQ6u3dyUWmfUrc4ptQ5yLusGp13+KhLua2v6bWdEWfa4Bpy85OJeWei6yKrU2tG5FDgotc8+++jrX/+6vv3tb2vixIm5KJOnTp06acyYMXr66ad18sknpx5/+umn9ZWvfKUgZUKmuDTMs+kQ2a3fNPfxN6RExzS0fHNaXLMykUhr2PpteBSy8rYrY7s6drPL5nwJ8p2s54bTzUM2U/Xakklt39muQb2q9c2FL6Q1XBLJjqlvTuVzSu2XOqbKxaHDJeXn/MhnYCfK7xOX+s2PbBusbq/3Og7ZHKeGplb9akW97v18QdtCHucg2QpxDKDETZDzIo4d6aDntbXuMe/eajfYFaReDBL4yGXwptBThOLaMc9VucznoJlbnVNqHeRc1g21NVWad8rItCl8ic+vdbf39tvOiFO9VogBtlzeKwtdF5nFLUgWZ4GDUg8++KDuv/9+HXXUUdp333317W9/WzNnzlTfvn0jLdjHH3+s9evXp/5dX1+vNWvWqGfPnho4cKAuuuginX766Ro7dqzGjx+vu+++W5s2bdK5554baTkQTlwa5tl2HJ2CQDc8/oZOHNU375WKXYNCkq579HVd/9jrqe8XpOERReUd5obm9F2k7M8XP9/J6dywu3mM2beH4/F0Cq55NSyMG9Wch9eq3fK3Qne4zHJ5cy9EYCeq8z0O9VuheR2HbI7TolWbMtbzKPRx9tO5jFNHI66Cnhdx60iHPa+ddm+1DnblYlp72Of6ESRjDNGwnoNWTnVOMXaQ3dqXua4bzNdsIiEd6pHFb5R10rDeWjG7zrWdEZd6rVADbLm+V8apLopTkCzOAgelTjjhBJ1wwgn64IMP9MADD+j+++/XFVdcoWOOOUbf/va3deKJJ2q33bJfqurll19WXV1d6t/GIuNnnHGG7r//fs2YMUMffPCBrrnmGjU0NGjEiBF67LHHtO+++2b92cg+ah6HhnkUHUen9F2vOeW5Ym1QWMtk/n5BGh7ZVN5hb2hu30XK/nxx+05u54bTzcPteIZtWMwYN1DD+3TVSQtWynoIXt28TeOH7hX8i+dALm7uhQzsZPt94lC/xYHXcQh7nIxzw66/5ff1hZpWGZeORpwFPS/i1pHO5vqvralSzz2dB7sOG9yzqALexZQxWkqcBkwNbnWOXRsnrlPRvc6vfNQNtTVVOn6U9/tlmz1ZiHqtkO2wcrtXxilIFleho0d77bWXLrzwQl144YW6/fbb9eMf/1iPPfaYevXqpXPPPVezZ89WdXX4E2vKlClK2nRUzc477zydd955oT8D9qJoZHhVNvm4AUbRcaytqdLsacM7RjFNCllxGg2KR19t0HWPvp72N/P3y0dkPtsbmtPIseT/GIc5l7zODbubh9PxzLZhMWpAD82emnmO/ezxNwuSjZcvxRzYKbfGlBOv4xD2OLl1uLxeX+hOchw6GnEX5rzI9f3M77poXTpV6oOPd2S9q6XTYNeqjVuLpl4kY7Rw3DLN/dQ55jZOVHVm1O16v+eXW92Qr2BbttmThcqgKWQ7rFjvlXEN4JaC0EGp9957Tw888IDuu+8+bdq0SV/96ld15pln6t1339W8efP0wgsv6KmnnoqyrMiDqBoZbpVNvjoNUXUcz5k8VEootYhkHCrO2poqHXdwra5/7HXX75fryHxUgb/jR1WpZedngW9OYc+lsOeG0/HMtmExsn9NxmNx7YhEpZgDO8XamIqa13EIe5ycOlwVDhsIGOLSSS50RyPuwp4XUdzPzMGllp1tGtyri+fuXXZr9yTUsb5MMkSbwG2wa9wg56niUYqiY1XMAwvFzu4aumTq/jq4f3ffdU5DU6te3vhhJHVmLgJbQc4vu7ohnwMU2WZPFup6KXQ7rNjulYUe9Cp1gYNSixcv1n333acnn3xSBx54oL73ve/pv//7v9W9e/fUcw455BCNHj06ynIiT6JsZEwa1lvzTx2likQiNQ87n52GKDuO50waqhNH9Y1Vxen0/SRp5YbGvETxo7yhmbOmlJDG7NvD9fnZnEu5CCpk07AodMOgEIo9sFNsjalc8ToOYY6T9dyokHTWpMGaNWGw6+vj1EkmVd9dIdYydAouSbuylrzWRZPp+RVJ6Y7TRnuuM2PHabBr1IAeOa8Xo+pYleN9K06yuYacFkmXgteZUbXrreflpVOHhz6/8j1AUazXQhzaYcVyr4zLoFcpCxyUmjVrlk499VQ9//zzGjdunO1zhgwZossvvzzrwiH/oqpYnRo9+e40RNlxjGPFaf1+y956XxPmPZu3KH7UNzSvEWuzbM+lKLa6j0ocGgaFUOyBnTjWCYXgdRzCHCevc8Pu2izWjkG5yub6CRpYcQsuWXmti2Zol9SzS+fQ38FpsCuX9WLU29CX430rTsJcQ16LpAetM6No19udlz974k1dOm24fvb4m4HPr3z3NYr5Wij2dli+xGnQq1QFDko1NDR4rhVVVVWlq666KnShUDhRVKxujZ58dxrKYe6v0SgpVBQ/qhtaIXZkcmvQ5TtNt1wbBgR24MTp3HC6Nou5YyCVx/0qCmHudV4LQ5t5rYtm97ywnM7xXNWLUXesyvW+Vcy81uwLWmdG0RZzOi8P7tfdcxe7XJUpqGK+FmiHeWPQK/cCB6U+++wzNTc3ZzyeSCTUuXNnderUKZKCoXCyrVjdGj3jh+6Vt05Duc39LfSChdl+Rpx2ZCpUgI+GAeDO69os1o5Bud2vshHmXucWXDLWh7JbM9J6nzEUW8DTkIuOFfet4mJ3DlRIuj3kVNQo2mJu52WY86tQAxRcC6Wr2Ae9ikHgoFT37t2VSCQc/96/f39961vf0lVXXaWKioqsCofCyeX6OMWwK1wxKvYofpx2ZCJNF4gnP9dmsXUMyvF+lY0w9wqv4JLfddGqO1Vo+872ogp4mtGxgtM5cNzBfUO/Z7ZtsVycl2HKRLYq3FjvBS0729TQ1Mq5EpHAQan7779fl19+ub71rW/pi1/8opLJpFatWqVf//rX+p//+R+9//77+vnPf67OnTvrsssuy0WZEXN+bi657jSUY1Ahqpt6oW7KYcufi3Op2AN85SBX5ymN0ngr1mvT7bxyu18Zf+d83CXsvcIruBT1umhxVazZhIhOLs6BbK+RfJbJrj4mW7X4FKK9VltTFWj9W/iXSCaTNsnMzo466iidc845+vrXv572+O9//3vdddddeuaZZ/Sb3/xGP/3pT/XGG284vEvxa25uVk1NjZqamtStW7dCFyeWGppaC9boaWhqTS34bahMJLRidl1RNcDCVLjZHPc43JRzed4EOZ6LVm3K6PRw04mHXJ2ncTj/4a3Yrk2v88rpfnXJtP1Tu7NxPmYqZBsDxcN6389XR5YBjniyq48nDetdEn2GclKo9lqp9C/zyW/MJHBQqrq6Wn//+9/1hS98Ie3xf/7znxo1apS2b9+u+vp6HXTQQdq+fXu40hcBglLxV2wdF6t8V7ilXtGGOZ50euInV+dpqZ//paZYrk2/55X1fnXJ1P11wxNvcD4CWbDe908e3U+PvLI55+0qBjjiyak+vvUbh+j8372S8fwHzz5M44fulccSwo9CttdWbmjUafe8mPE454ozvzGTwNP3+vfvr3vvvVfz5s1Le/zee+/VgAEDJEkffPCBevToEfStgUgVc4q63Rojcx5eqy6dd9OYEAtR+hH3KY/ZjDqGXbOllKZsGIp99DZX52ncz3+kK5Zr0+95Zb1fcT4C2bG77z/8t82pv+dq7TbWiIsvp3pVnwcPs5kWXqi2VbG36cIo5P2xWJcQKAaBg1I///nP9bWvfU2PP/64xo0bp0QioVWrVumNN97Q//7v/0qSVq1apRkzZkReWCCoYum4WNlVuO2Szv/dKzkbdYtzRZvtqGMpdfCyaYCUwuhtrs7TOJ//CCcOjfUg55X1fsX5CIRnd9+3ykU7oJTaG6XGqT4eM6hHVmuyFqptVQptujAK2V5js4jcCbw93oknnqi33npLxx57rD788EM1NjZq2rRpeuONN3T88cdLkr773e/q5ptvjrywQK40NLVq5YZGNTS12v4734wK144x6hZ12YyKtvLz3TXjUtE6jToG+f52x7MYO3iLVm3ShHnP6rR7XtSEec9q0apNvl8bxXGMg2VvvS/zpPNEQpGcp3E9/xFONteKIYr7QNjzivOxMAp970d03NpRhly0A0qlvVGK3OrVGeMGasXsOj149mFaMbvOd3CnUG2rUmnThVHo+2PYcwXuAmVKffrppzr66KN11113ae7cubkqE5BXhVpzwI3T9tWGXI26xXHKYxSjjqUwspHtlIBSGL01joH5aySS0qRhvSN5/zie/wguiukzUY5Ahz2vjNet3rhVSkhj9mVZhFwq16yDUmV33z9pdF/98ZV3c9oOKIX2Rilzq4/DzK4oVNuqFNp02Sh0e61YZ+LEWaCg1O67765169YpkfAYegCKRKHWHPDD3CH5wUOv5CVNNQ7TXayiStMt9A0sW9k2QEpheprTtNYoG2E0NIpfttdKLtaECXteRbn1dBzr91wIu2st6wCVHrv7/sXH7J/zdkCxtzdKXZT3+UK1rUqhTZct2mulJfD0vZkzZ+ree+/NRVmAvAuy5kAh1NZU6fhRffOSphrFdJdciDJNt7amSuOH7lWUN7FspwQUOt05LPN0GqZFwI9szxO3oFY+RTk9I671e9TCfs+4/OaInvW+X1tTldpMIJdTnYq5vQH/CtW2KtY2HeAk8ELnO3fu1MKFC/X0009r7Nix6tKlS9rfWUsKxaKhqVUffLwjY6TBKg6d3lyPusV9lJhRx2imBBTbcbSbTsO0CHjJ9lqJywh0VNMz4l6/RyWb7xmX3xy5xzRNRK1QbSu7zy2XjFiUnsBBqXXr1unQQw+VJL311ltpf2NaX3kohQrP3ChJqGOx5GRSeVtzIIxcpqkWw9x00nSjafgUy3F06mCumF2nFbPriiawVs4Kea/I5lqJy5owUQVKiqF+j0I23zMuvzlyKx8B2lJoIyO4QrWtzJ9LwBXFLHBQasmSJbkoB4pEKVR41kZJUlJFUrrjtNE6dN8eeVtzIE4YJQ4v3w3QYgkqufFzzNw6mEyJiL843CuyuVbikFUYVaCkXOr3bL9nHH5z5FauA7RxqPdQnsolIxalK3BQyrB+/Xpt2LBBkyZNUlVVlZLJJJlSJa4YKrywnd12ST27dE5bcyAu3ykfGCUOhwZocH6PWbl0pEtRMdwr/IjDfSCq7MhyqN+j+J5x+M2RO7m8r5RKvYfiVC4ZsShdgYNSH3zwgb7+9a9ryZIlSiQS+uc//6khQ4borLPOUvfu3XXTTTflopyIgbhXeHR2s8MocTA0QIMLcszKpSNdiuJ+ryg2UQRKyqV+L5fviXByeV+h3kMh0bdBsQsclLrwwgu1++67a9OmTTrggANSj8+YMUMXXnghQakSFucKj85uNBgl9o8GaHBBjxkdzOIU53tFOSuX+r1cvifCCXpf8TtFn3oPhUTfBsUucFDqqaee0pNPPqn+/funPf6FL3xB//rXvyIrGOInzhUenV3kGw3Q4MIcMzqYxSfO9woA8HtfCTJFn3oPhZaLvg0L9yNfAgelWlpaVF2d2YFobGxU586dIykU4iuuwRw6u8g3GqDBcczKR1zvFQDgR5gp+tR7KLQo+zasm4p8ChyUmjRpkh544AFde+21kqREIqH29nbdeOONqquri7yAyL2gUfA4BnPo7KIQaIAGxzHLXrGMXMbxXgEAfoSdok+9h1LAuqnIt8BBqRtvvFFTpkzRyy+/rJ07d+qSSy7Ra6+9pg8//FDPP/98LsqIHCqlKDidXRQCDdDgCnXM8hHMyfVnlFKdDQBxxRR9lDPWTUW+VQR9wYEHHqhXX31VX/ziF/Vf//Vfamlp0fTp0/XKK69o6NChuSgjcsQpCt7Q1FrYgmWhtqZK44fuRYUJIM2iVZs0Yd6zOu2eFzVh3rNatGpT0X1GKdbZiE5DU6tWbmjkfEDRieO5a2TgVyYSkkQGPsqKEZQ1IyiLXAqcKSVJffr00U9+8pOoy4I8IwoOoBzkIw09H59BnQ0nZNChWMX53CUDH+WKZVGQb6GCUtu2bdNLL72kLVu2qL29Pe1vM2fOjKRgyD1SkwH4VSzrGNnJRzAnH59BnQ07rP2BYlUM5y5T9FGuCMoinwIHpf7yl7/om9/8plpaWtS1a1clErty+xKJBEGpIkIUHIAfcR7J9iMfwZx8fAZ1NuyQQYdiVSrnbjEP2gBuCMoiXwIHpX70ox/p29/+tq6//npVVzM6W+zsouDcXKPDsUSxK4aRbC/5CObkK2DEyCWsyKBDsSqFc7fYB20AIA4SyWQy6f20Xbp06aK1a9dqyJAhuSpTUWhublZNTY2amprUrVu3QhcnMtxco+N0LAlUoZis3NCo0+55MePxB88+TOOH7lWAEoXX0NSa82BOPj4DsFq0alNGQJR7N4pBMZ+7DU2tmjDv2Yyg2orZddT/ACD/MZPAmVLHHHOMXn755bIPSpWiUsiIiAunY7mt9VPd8PgbBP1QNEphJNuQjzR0Ut1RCGTQoVgV87lbKtMP4R8Dy0BuBA5KHXfccfrxj3+sf/zjHxo5cqR23333tL+feOKJkRUO+cXNNTpOx3Le428oSdAPLuLW4GEdI6A4EBB1F7e6FbsU67lbSoM28MZsEiB3Agelzj77bEnSNddck/G3RCKhtra27EuFguDmGh27Y1khEfSDq7g2eIp5JBsA4lq3orgxaFM+mE0C5FZF0Be0t7c7/o+AVHEzbq6Vn++oyM01PLtjeem04apIpD+PoB8MTg2ehqbWwhbsc7U1VRo/dC/qAwBFJe51K4rbjHEDtWJ2nR48+zCtmF3nGuxsaGrVyg2NnHtFyG02CYDsBc6UMvvkk0+0xx57RFUWxAAZEdGxO5bdq3dnRC2AcppuwfRZIFM51QHIDepW5Jqf6Ydk6xU3ZpMAuRU4KNXW1qbrr79ev/zlL/Wf//xHb731loYMGaIrrrhCgwYN0plnnpmLciKPinVufxxZjyVBP//KrQFHgwdIV251QFwVe2CQuhWFxtSv4sdUTSC3Ak/f++lPf6r7779fP/vZz9SpU6fU4yNHjtTChQsjLRyKE+nJ7pgG5a0cp1swfRbYpRzrgDhatGqTJsx7Vqfd86ImzHtWi1ZtKnSRAqNuRaEx9as0BJmqCSCYwJlSDzzwgO6++24dddRROvfcc1OPH3zwwXrjjTciLRyKDyPbiEK5Trcgky6YYs/ggLNyrQPipJSyO6hbUUhk65UOZpMAuRE4KLV582btt99+GY+3t7fr008/jaRQKE6l1IBFYZVzA44Gjz8EwEtbOdcBcVFqgUHqVhQKU78AwF3g6XsHHXSQli9fnvH4H/7wB40ePTqSQqE4kZ4MN0GmdTLdAm6Y2hWduE63pg4oPCMwaEZgEAiHqV8A4CxwptRVV12l008/XZs3b1Z7e7sWL16sN998Uw888ID+7//+LxdlRJFgZBtOwmS1MN0CTkotg6NQ4p5tRh1QWGR3ANEiWw8A7CWSyWTS+2npnnzySV1//fVavXq12tvbdeihh+rKK6/U0UcfnYsyxlJzc7NqamrU1NSkbt26Fbo4sbFo1aaMBmycOjnIv4amVk2Y92xGsHLF7DoaZwiFcyp7HEP41dDUSmAQAAAE5jdmEjhTSpKOOeYYHXPMMaELh9LFyDasyGpB1MjgyB7XJfwiuwMAAORSqKAU4IYGLMyY1olcIACeHa5LAAAAxEHghc4BIAgWLEau1NZUafzQvTiXQuC6RFBxXRQfAAAUt1BrSoE1pYCgWJcEiB+uS/gR90XxAQDRamhqVX1jiwb36kL7AKHldE0pAAiKaZ1A/HBdwktDU2sqICV1TPm8bPE6TRrWm3MHAEoQAxHIN6bvAQAQE0yRQty4LYoPAMi9fLYNnAYiaJcglyLLlPrTn/6kpqYmzZw5M6q3BACgbORjZJJ0fATFovgAUDj5zlpid14UQmSZUpdeeqlmzZoV1dsBAFA28jEyuWjVJk2Y96xOu+dFTZj3rBat2hTZe6N0sSg+ABRGIbKWjIEIMwYikGuRZUq98cYbUb0VAABlJdcjk+W2LhAZYdGaMW6gJg3rzaL4AJBHhchaMgYiLlu8Tm3JJAMRyAsWOgcAoMByPUWqnNLxWaA1N1gUHwDyq1DTpxmIyA4DY8H5Ckq9+uqrvt/w4IMPDl0YAADKUa5HJstlXaByywgDAJSuQmYtMRARDgNj4fgKSh1yyCFKJBJKJpO2fzf+lkgk1NbWFmkBAQAoB7kcmSyXdPxyyggDAJQ+spaKBwNj4fkKStXX1+e6HAAAlL1cjkyWQ8O2XDLCAADlg6yl4sDAWHi+glL77rtvrssBAAByrNQbtuWSEQYAAOKFgbHwKsK86De/+Y0mTJigvn376l//+pckaf78+frTn/4UaeEAAACCmDFuoFbMrtODZx+mFbPrWMsBAADknDEwVplISBIDYwEE3n3vzjvv1JVXXqkLLrhAP/3pT1NrSHXv3l3z58/XV77ylcgLCQAA4FepZ4QBAID4KYelEnIhcKbU7bffrnvuuUeXX365KisrU4+PHTtWa9eujbRwAAAAAAAAxaC2pkrjh+5FQCqAwEGp+vp6jR49OuPxzp07q6WlJZJCAQAAAAAAoLQFDkoNHjxYa9asyXj88ccf14EHHhhFmQAAAAAAAFDiAq8p9eMf/1jf+9739MknnyiZTOqll17Sgw8+qLlz52rhwoW5KCMAAAAAAABKTOCg1KxZs/TZZ5/pkksu0fbt23XaaaepX79+uvXWW3XqqafmoowAAAAAAAAoMYlkMpkM++LGxka1t7dr7733liRt3rxZ/fr1i6xwcdbc3Kyamho1NTWpW7duhS4OAAAAAABALPiNmQReU8qsV69e2nvvvfXee+/p+9//vvbbb79s3g4AAAAAAABlwndQatu2bfrmN7+p3r17q2/fvrrtttvU3t6uK6+8UkOGDNELL7ygX/3qV7ksKwAAAAAAAEqE7zWlLrvsMi1btkxnnHGGnnjiCV144YV64okn9Mknn+jxxx/X5MmTc1lOAAAAAAAAlBDfQalHH31U9913n7785S/rvPPO03777adhw4Zp/vz5OSweAAAAAAAASpHv6XvvvvuuDjzwQEnSkCFDtMcee+iss87KWcEAAAAAAABQunwHpdrb27X77run/l1ZWakuXbrkpFAAAAAAAAAobb6n7yWTSX3rW99S586dJUmffPKJzj333IzA1OLFi6MtIQAE1NDUqvrGFg3u1UW1NVWFLg4AAAAAwIbvoNQZZ5yR9u///u//jrwwAJCtRas2ac7itWpPShUJae70kZoxbmChiwUAAAAAsEgkk8lkoQtRjJqbm1VTU6OmpiZ169at0MUBoI4MqQnznlW7qVarTCS0YnYdGVMAAAAAkCd+Yya+15QCgLirb2xJC0hJUlsyqY2N2wtTIAAAAACAI4JSQAE1NLVq5YZGNTS1FrooJWFwry6qSKQ/VplIaFCv6sIUCAAAAADgiKAUUCCLVm3ShHnP6rR7XtSEec9q0apNhS5S0autqdLc6SNVmeiITFUmErp++gim7gEAAABADLGmVEisKYVssPZRbjU0tWpj43YN6lXN8QQAAACAPPMbM/G9+x6A6LitfUQQJXu1NVUcRwAAAACIOabvAQXA2kcAAAAAgHJHUAooANY+AgAAAACUO6bvAQUyY9xATRrWm7WPAAAAAABliaAUUECsfQQAAFA4DU2tqm9s0eBeXWiTAWWAaz5+CEoBAAAAAdGxKX6LVm3SnMVr1Z6UKhLS3OkjNWPcwEIXC4gM9VQ6rvl4IigFAAAABEDHpvg1NLWmfkNJak9Kly1ep0nDetN5R0mgnkrHNR9fZbvQ+c9//nMddNBBGjFihP6//+//K3RxAAAAUAScOjYNTa2FLZiLhqZWrdzQGOsy5lt9Y0vqNzS0JZPa2Li9MAUCIlSM9VSucc3HV1lmSq1du1a/+93vtHr1aknSUUcdpeOPP17du3cvbMEAAAAQa24dmziOtpMtYW9wry6qSCjtt6xMJDSoV3XhCgVEpNjqqXzgmo+vssyUev3113X44Ydrjz320B577KFDDjlETzzxRKGLBQAAgJgzOjZmce3YkC3hrLamSnOnj1RlouPHrEwkdP30EWXbYUdpKaZ6Kl+45uMrlkGpZcuW6YQTTlDfvn2VSCT0xz/+MeM5CxYs0ODBg7XHHntozJgxWr58ue/3HzFihJYsWaJt27Zp27ZtevbZZ7V58+YIvwEAAABKUTF1bJiu4m7GuIFaMbtOD559mFbMriODDCWjmOqpfOKaj6dYTt9raWnRqFGjNGvWLJ1yyikZf1+0aJEuuOACLViwQBMmTNBdd92ladOm6R//+IcGDuw4scaMGaMdO3ZkvPapp57SgQceqB/84Ac68sgjVVNTo3Hjxmm33WJ5KAAAABAzM8YN1KRhvbWxcbsG9aqObUeP6SreamuqYvv7Adkolnoq37jm4yeRTCaT3k8rnEQioUceeUQnnXRS6rEvfelLOvTQQ3XnnXemHjvggAN00kknae7cuYE/46yzztLJJ5+s4447zvE5O3bsSAtyNTc3a8CAAWpqalK3bt0CfyYAAACQa4tWbdJli9epLZlMZUuQHQAAyLXm5mbV1NR4xkyKLj1o586dWr16tWbPnp32+NFHH62VK1f6fp8tW7Zo77331ptvvqmXXnpJv/zlL12fP3fuXP3kJz8JVWYAAMpdQ1Or6htbNLhXF0YogTwq5WwJ6hUAKH5FF5RqbGxUW1ub9tlnn7TH99lnH7333nu+3+ekk07Stm3b1KVLF913332e0/fmzJmjiy66KPVvI1MKAAC4Y/cvoLBKcboK9QoAlIaiC0oZEon07QSSyWTGY26CZFVJUufOndW5c+dArwEAoNw57f41aVjvkuskA8gP6hUApaTcsz6LLijVq1cvVVZWZmRFbdmyJSN7CgAAFJbb7l/l2PACkD3qFQClgqxPqaLQBQiqU6dOGjNmjJ5++um0x59++mkdfvjhBSoVAACwY+z+ZcbuXwCyQb0CoBQ4ZX02NLUWtmB5Fsug1Mcff6w1a9ZozZo1kqT6+nqtWbNGmzZtkiRddNFFWrhwoX71q1/p9ddf14UXXqhNmzbp3HPPLWCpAQCAVW1NleZOH6nKz6fYG7t/kc0AICzqFQClwC3rs5zEcvreyy+/rLq6utS/jQXGzzjjDN1///2aMWOGPvjgA11zzTVqaGjQiBEj9Nhjj2nfffctVJEBAICDUt79C0BhUK8AKHZG1qc5MFWOWZ+JZDKZ9H4arJqbm1VTU6OmpiZ169at0MUBAAAAAABFZNGqTbps8Tq1JZOprM9SWVPKb8wklplSAAAAAAAApYysT4JSAAAAAAAABVFbU1WWwShDLBc6BwAAAAAAQGkjKAUAAAAAAIC8IygFAAAAAACAvCMoBQAAAAAAgLwjKAUAAAAAAIC8IygFAAAAAACAvCMoBQAAAAAAgLwjKAUAAAAAAIC8IygFAAAAAACAvCMoBQAAAAAAgLwjKAUAAAAAAIC8IygFAAAAAACAvCMoBQAAAAAAgLwjKAUAAAAAAIC8IygFAAAAAACAvCMoBQAAAAAAgLwjKAUAAAAAAIC8IygFAAAAAACAvCMoBQAAAAAAgLwjKAUAAAAAAIC8IygFAABQ5hqaWrVyQ6MamloLXRQAAFBGdit0AQAAAFA4i1Zt0pzFa9WelCoS0tzpIzVj3MBCFwsAAJQBMqUAAADKVENTayogJUntSemyxevImAIAAHlBUAoAAKBM1Te2pAJShrZkUhsbtxemQAAAoKwQlAIAAChTg3t1UUUi/bHKREKDelUXpkAAAKCsEJQCAAAoU7U1VZo7faQqEx2RqcpEQtdPH6HamqoClwwAAJQDFjoHAAAoYzPGDdSkYb21sXG7BvWqJiAFAADyhqAUAABAmautqSIYBQAA8o7pewAAAAAAAMg7glIAAAAAgLLS0NSqlRsa1dDUWuiiAGWN6XsAAAAAgLKxaNUmzVm8Vu1JqSIhzZ0+UjPGDSx0sYCyRKYUAAAAAKAsNDS1pgJSktSelC5bvI6MKaBACEoBAAAAAMpCfWNLKiBlaEsmtbFxe2EKBJQ5glIAAAAAgLIwuFcXVSTSH6tMJDSoV3VhCgSUOYJSAAAAAICyUFtTpbnTR6oy0RGZqkwkdP30EaqtqSpwyYDyxELnAAAAAICyMWPcQE0a1lsbG7drUK9qAlJAARGUAgAAAACUldqaKoJRQAwwfQ8AAAAAAAB5R1AKAAAAAIAi19DUqpUbGtXQ1FroogC+MX0PAAAAAIAitmjVJs1ZvFbtSakiIc2dPlIzxg0sdLEAT2RKAQAAAABQpBqaWlMBKUlqT0qXLV5HxhSKAkEpAAAAAACKVH1jSyogZWhLJrWxcXthCgQEQFAKAAAAAIAiNbhXF1Uk0h+rTCQ0qFd1YQoEBEBQCgAAAACAIlVbU6W500eqMtERmapMJHT99BGqrakqcMkAbyx0DgAAAABAEZsxbqAmDeutjY3bNahXNQEpFA2CUgAAAAAAFLnamiqCUSg6TN8DAAAAAABA3hGUAgAAAAAAQN4RlAIAAAAAAEDeEZQCAAAAAABA3hGUAgAAAAAAQN4RlAIAAAAAAEDeEZQCAAAAAABA3hGUAgAAAAAAQN4RlAIAAAAAAEDeEZQCAAAAAABA3hGUAgAAAAAAQN4RlAIAAAAAAEDeEZQCAAAAAABA3hGUAgAAAAAAQN4RlAIAAAAAAEDeEZQCAADwqaGpVSs3NKqhqbXQRQEAACh6uxW6AAAAAMVg0apNmrN4rdqTUkVCmjt9pGaMG1joYhVcQ1Or6htbNLhXF9XWVBW6OAAAoIgQlAIAAPDQ0NSaCkhJUntSumzxOk0a1rusAzEE6gAAQDaYvgcAAOChvrElFZAytCWT2ti4vTAFigGnQB1TGwEAgF8EpQAAADwM7tVFFYn0xyoTCQ3qVV2YAsUAgToAAJAtglIAAAAeamuqNHf6SFUmOiJTlYmErp8+oqyn7hGoAwAA2WJNKQAAAB9mjBuoScN6a2Pjdg3qVV3WASlpV6DussXr1JZMEqgDAACBEZQCAADwqbamiqCLCYE6AACQDYJSAAAACI1AHQAACIs1pQAAAAAAAJB3BKUAAAAAAEBJamhq1coNjWpoai10UWCD6XsAAAAAAKDkLFq1SXMWr1V7UqpISHOnj9SMcQMLXSyYkCkFAAAAAEAekLWTPw1NramAlCS1J6XLFq/j2McMmVIAAAAAAOQYWTv5Vd/YkgpIGdqSSW1s3M4GHTFCphQAAAAAADlE1k7+De7VRRWJ9McqEwkN6lVdmALBFkEpAACAIsLUDwAoPm5ZO8iN2poqzZ0+UpWJjshUZSKh66ePIEsqZpi+BwAAUCSY+gEAxcnI2jEHpsjayY2GplbVN7ZocK8umjFuoCYN662Njds1qFc1AakYIigFAABQBJymfkwa1ptGNgDEnJG1c9nidWpLJsnayRGnwRuOc3wRlAIAACgCLNgKAMWNrJ3cYvCmOBGUAgAAKAJM/QBQjMxTqQgMdGRMcRxyg8Gb4sRC5wAAAEWABVsLj0XmgWAWrdqkCfOe1Wn3vKgJ857VolWbCl0khFQM9V+Uu+0Vw/ctFYlkMpn0fhqsmpubVVNTo6amJnXr1q3QxQEAAGWioamVqR8FwCLzpY1snug1NLVqwrxnM7I7V8yu4xgXmVzVf7m47hat2pSxbpffshrlWbu5STc8/gb1fZb8xkzKYvreySefrKVLl+qoo47S//7v//r+GwCgfNFBQVwx9SP/WKektBFwzA2mUpWGXNV/ubruwq7bZS6PGfV97pXF9L0f/OAHeuCBBwL/DQBQnophugFp5UD+uHWuUdycOtzUrdmLcioVCicX9V+ur7vamiqNH7qX7yCStTxW1Pe5VRZBqbq6OnXt2jXw3wAA5acYOijFEDQDSgmd69JFwDF3WAevNOSi/ovbdWdXHjPq+9wqeFBq2bJlOuGEE9S3b18lEgn98Y9/zHjOggULNHjwYO2xxx4aM2aMli9fnv+ClglG3nfhWADlKW4NJatiCJoBpYbOdeki4JhbM8YN1IrZdXrw7MO0YnYd0yKLUC7qv7hdd3blMVDf517B15RqaWnRqFGjNGvWLJ1yyikZf1+0aJEuuOACLViwQBMmTNBdd92ladOm6R//+IcGDuyo1MaMGaMdO3ZkvPapp55S3759c/4dSgXz6XfhWADly2iYWBdmjUsHhTU6gMIIu04J4s3ocFsXRub3jQ7r4BW/qOu/uF13duW5ZOr+Orh/d+r7PCh4UGratGmaNm2a499vvvlmnXnmmTrrrLMkSfPnz9eTTz6pO++8U3PnzpUkrV69Oufl3LFjR1rgq7m5OeefmU8s4LkLxwIob3FrKFnFPWgGlDI616XJrcPNphdAh6jrv7gF+uNWnnJS8KCUm507d2r16tWaPXt22uNHH320Vq5cmdeyzJ07Vz/5yU/y+pn5xMj7LhwLAHFumMQ9aAYAxciuw03mPJBbcQv0x6085SLWQanGxka1tbVpn332SXt8n3320Xvvvef7fY455hj97W9/U0tLi/r3769HHnlE48aN8/yb2Zw5c3TRRRel/t3c3KwBAwaE/Gbxw8j7LhwLAFK8GyZxDpoBQCkgcx7Fgmw+FLtYB6UMiUT6qmPJZDLjMTdPPvlkqL+Zde7cWZ07d/b9mcWGkfddOBYAikGcg2YAUOzInEcxIJsPpSDWQalevXqpsrIyIytqy5YtGdlTyB4j77twLAAAAMoXmfOIuzhn85G9hSAqCl0AN506ddKYMWP09NNPpz3+9NNP6/DDDy9QqUpbbU2Vxg/di8pDHAsAAIByZWTOV34+O4PMecSNWzZfIS1atUkT5j2r0+55URPmPatFqzYVtDyIv4JnSn388cdav3596t/19fVas2aNevbsqYEDB+qiiy7S6aefrrFjx2r8+PG6++67tWnTJp177rkFLDUAAACAUkbmPOIsjtl8cc7eQnwVPCj18ssvq66uLvVvYzHxM844Q/fff79mzJihDz74QNdcc40aGho0YsQIPfbYY9p3330LVWQAAAAAZYD1+xBXcVwHl7XYEEYimUwmvZ8Gq+bmZtXU1KipqUndunUrdHEAAAAAAGWmoak1Ntl8DU2tmjDv2YzsrRWz6wpeNuSf35hJrNeUAgAAAAAA9uK0Di5rsSGMgk/fAwAAAAAAxY+12BAUQSkAAAAAABAJ1mJDEEzfAwAAAAAAQN4RlAIAAAAAAEDeEZQCAAAAAABA3hGUAgAAAAAAQN4RlAIAAAAAAEDeEZQCAAAAAABA3hGUAgAAAAAAQN4RlAIAAAAAAEDeEZQCAAAAAABA3hGUAgAAAAAAQN4RlAIAAAAAAEDeEZQCAAAAAABA3hGUAgAAAAAAQN4RlAIAAAAAoEw1NLVq5YZGNTS1FrooKEO7FboAAAAAAAAg/xat2qQ5i9eqPSlVJKS500dqxriBhS4WygiZUgAAAAAAlJmGptZUQEqS2pPSZYvXkTGFvCIoBQAAAABAmalvbEkFpAxtyaQ2Nm4vTIFQlghKAQAAAABQZgb36qKKRPpjlYmEBvWqLkyBUJYISgEAAAAAUGZqa6o0d/pIVSY6IlOViYSunz5CtTVVBS4ZygkLnQMAAAAAUIZmjBuoScN6a2Pjdg3qVU1ACnlHUAoAAAAAgDJVW1NFMAoFw/Q9AAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAABBIQ1OrVm5oVENTa6GLgiK2W6ELAAAAAAAAiseiVZs0Z/FatSelioQ0d/pIzRg3sNDFQhEiUwoAAAAAAPjS0NSaCkhJUntSumzxOjKmEApBKQAAAAAA4Et9Y0sqIGVoSya1sXF7YQqEokZQCgAAAAAA+DK4VxdVJNIfq0wkNKhXdWEKhKJGUAoAAAAAAPhSW1OludNHqjLREZmqTCR0/fQRqq2pKnDJUIxY6BwAAAAAAPg2Y9xATRrWWxsbt2tQr2oCUgiNoBQAAAAAAAiktqaKYBSyxvQ9AAAAAAAA5B1BKQAAAAAAAOQdQSkAAAAAAADkHUEpAAAA4P9v796Dor7u/4+/FoEFuawCBURUwHohQSMBUzW2xKgxE41x4qSJV6jTWNsQoZkm0VxJUsXOpLRjEnMxKamNCTSjzRim04pVQSuRFKSCNJILokENTUUgRUDY8/0jPz+/bgRjjNkl7vMxszPuOe/d897PfN6zw9vz+SwAAHA7mlIAAAAAAABwO5pSAAAAAAAAcDuaUgAAAAAAAHA7mlIAAAAAAABwO5pSAAAAAAAAcDuaUgAAAAAAAHA7mlIAAAAAAABwO5pSAAAAAAAAcDuaUgAAAAAAAHA7mlIAAAAAAABwO5pSAAAAAAAAcDuaUgAAAAAAAHA7mlIAAAAAAABwO5pSAAAAAAAAcDtfTyfwbWWMkSS1trZ6OBMAAAAAAID+41yv5FzvpC80pS5RW1ubJGnYsGEezgQAAAAAAKD/aWtrk8Ph6HPeZr6sbYVeOZ1OHT9+XCEhIbLZbJ5O55K1trZq2LBhOnbsmEJDQz2dDtCvUB9A36gPoG/UB3Bh1AjQtyulPowxamtrU0xMjHx8+r5zFDulLpGPj49iY2M9ncZlExoa+q0+4YFvEvUB9I36APpGfQAXRo0AfbsS6uNCO6TO4UbnAAAAAAAAcDuaUgAAAAAAAHA7mlJezm636/HHH5fdbvd0KkC/Q30AfaM+gL5RH8CFUSNA37ytPrjROQAAAAAAANyOnVIAAAAAAABwO5pSAAAAAAAAcDuaUgAAAAAAAHA7mlJebMOGDYqPj1dAQIBSUlK0Z88eT6cEfC25ubmaOHGiQkJCFBkZqXnz5unw4cMuMcYY5eTkKCYmRoGBgbrhhht06NAhl5jOzk7de++9ioiIUFBQkObOnauPP/7YJaa5uVlLliyRw+GQw+HQkiVLdPr0aZeYo0eP6tZbb1VQUJAiIiK0cuVKdXV1fSOfHfiqcnNzZbPZlJ2dbY1RH/BmjY2NWrx4scLDwzVw4EBNmDBBFRUV1jz1AW/W3d2tRx55RPHx8QoMDFRCQoKefPJJOZ1OK4YagbcoLS3VrbfeqpiYGNlsNr311lsu8/2tFqqrq5WWlqbAwEANHTpUTz75pPrVrcUNvFJBQYHx8/MzGzduNLW1tSYrK8sEBQWZhoYGT6cGXLJZs2aZ/Px8U1NTY6qqqszs2bPN8OHDzWeffWbFrFu3zoSEhJgtW7aY6upqc+edd5ohQ4aY1tZWK2bFihVm6NChpri42FRWVppp06aZa665xnR3d1sxN998s0lKSjL79u0z+/btM0lJSWbOnDnWfHd3t0lKSjLTpk0zlZWVpri42MTExJjMzEz3HAzgAsrLy01cXJwZP368ycrKssapD3irU6dOmREjRpiMjAyzf/9+U19fb3bs2GE++OADK4b6gDf75S9/acLDw01RUZGpr683b775pgkODja//e1vrRhqBN7iz3/+s3n44YfNli1bjCTzpz/9yWW+P9VCS0uLiYqKMnfddZeprq42W7ZsMSEhIebpp5/+5g7QV0RTyktdd911ZsWKFS5jY8eONatWrfJQRsDl19TUZCSZkpISY4wxTqfTREdHm3Xr1lkxHR0dxuFwmBdeeMEYY8zp06eNn5+fKSgosGIaGxuNj4+P+ctf/mKMMaa2ttZIMu+8844VU1ZWZiSZ9957zxjz+ZeVj4+PaWxstGLeeOMNY7fbTUtLyzf3oYEv0dbWZkaNGmWKi4tNWlqa1ZSiPuDNHnzwQTN16tQ+56kPeLvZs2ebZcuWuYzdfvvtZvHixcYYagTe64tNqf5WCxs2bDAOh8N0dHRYMbm5uSYmJsY4nc7LeCQuHZfveaGuri5VVFTopptuchm/6aabtG/fPg9lBVx+LS0tkqSwsDBJUn19vU6ePOly7tvtdqWlpVnnfkVFhc6ePesSExMTo6SkJCumrKxMDodD3/ve96yYSZMmyeFwuMQkJSUpJibGipk1a5Y6OztdLgcB3O2ee+7R7NmzNWPGDJdx6gPebNu2bUpNTdUdd9yhyMhIJScna+PGjdY89QFvN3XqVP3tb39TXV2dJOmf//yn9u7dq1tuuUUSNQKc099qoaysTGlpabLb7S4xx48f15EjRy7/AbgEvp5OAO736aefqqenR1FRUS7jUVFROnnypIeyAi4vY4zuu+8+TZ06VUlJSZJknd+9nfsNDQ1WjL+/vwYPHnxezLnXnzx5UpGRkeetGRkZ6RLzxXUGDx4sf39/6gweU1BQoMrKSr377rvnzVEf8GYfffSRnn/+ed1333166KGHVF5erpUrV8put2vp0qXUB7zegw8+qJaWFo0dO1YDBgxQT0+P1qxZowULFkjiOwQ4p7/VwsmTJxUXF3feOufm4uPjL+VjXlY0pbyYzWZzeW6MOW8M+LbKzMzUwYMHtXfv3vPmLuXc/2JMb/GXEgO4y7Fjx5SVlaXt27crICCgzzjqA97I6XQqNTVVa9eulSQlJyfr0KFDev7557V06VIrjvqAtyosLNRrr72m119/XVdffbWqqqqUnZ2tmJgYpaenW3HUCPC5/lQLveXS12s9gcv3vFBERIQGDBhw3v8kNDU1nddpBb6N7r33Xm3btk27du1SbGysNR4dHS1JFzz3o6Oj1dXVpebm5gvGfPLJJ+et++9//9sl5ovrNDc36+zZs9QZPKKiokJNTU1KSUmRr6+vfH19VVJSovXr18vX19flf83+F/UBbzBkyBBdddVVLmOJiYk6evSoJL4/gPvvv1+rVq3SXXfdpXHjxmnJkiX6+c9/rtzcXEnUCHBOf6uF3mKampoknb+by1NoSnkhf39/paSkqLi42GW8uLhYU6ZM8VBWwNdnjFFmZqa2bt2qnTt3nrcdNT4+XtHR0S7nfldXl0pKSqxzPyUlRX5+fi4xJ06cUE1NjRUzefJktbS0qLy83IrZv3+/WlpaXGJqamp04sQJK2b79u2y2+1KSUm5/B8e+BLTp09XdXW1qqqqrEdqaqoWLVqkqqoqJSQkUB/wWtdff70OHz7sMlZXV6cRI0ZI4vsDaG9vl4+P65+OAwYMkNPplESNAOf0t1qYPHmySktL1dXV5RITExNz3mV9HuO+e6qjPykoKDB+fn7mlVdeMbW1tSY7O9sEBQWZI0eOeDo14JL99Kc/NQ6Hw+zevducOHHCerS3t1sx69atMw6Hw2zdutVUV1ebBQsW9PoTrbGxsWbHjh2msrLS3Hjjjb3+ROv48eNNWVmZKSsrM+PGjev1J1qnT59uKisrzY4dO0xsbCw/V4x+5X9/fc8Y6gPeq7y83Pj6+po1a9aY999/32zevNkMHDjQvPbaa1YM9QFvlp6eboYOHWqKiopMfX292bp1q4mIiDAPPPCAFUONwFu0tbWZAwcOmAMHDhhJJi8vzxw4cMA0NDQYY/pXLZw+fdpERUWZBQsWmOrqarN161YTGhpqnn76aTccqYtDU8qLPffcc2bEiBHG39/fXHvttaakpMTTKQFfi6ReH/n5+VaM0+k0jz/+uImOjjZ2u9384Ac/MNXV1S7vc+bMGZOZmWnCwsJMYGCgmTNnjjl69KhLzH/+8x+zaNEiExISYkJCQsyiRYtMc3OzS0xDQ4OZPXu2CQwMNGFhYSYzM9Pl51gBT/tiU4r6gDd7++23TVJSkrHb7Wbs2LHmpZdecpmnPuDNWltbTVZWlhk+fLgJCAgwCQkJ5uGHHzadnZ1WDDUCb7Fr165e/+ZIT083xvS/Wjh48KD5/ve/b+x2u4mOjjY5OTnG6XRe9uNyqWzG/L+7XAEAAAAAAABuwj2lAAAAAAAA4HY0pQAAAAAAAOB2NKUAAAAAAADgdjSlAAAAAAAA4HY0pQAAAAAAAOB2NKUAAAAAAADgdjSlAAAAAAAA4HY0pQAAAAAAAOB2NKUAAAC+hpycHE2YMMFj6z/66KNavny5x9a/WM8++6zmzp3r6TQAAEA/YjPGGE8nAQAA0B/ZbLYLzqenp+vZZ59VZ2enwsPD3ZTV//fJJ59o1KhROnjwoOLi4ty+/lfR2dmpuLg4vfnmm5o6daqn0wEAAP2Ar6cTAAAA6K9OnDhh/buwsFCPPfaYDh8+bI0FBgYqODhYwcHBnkhPr7zyiiZPnuzxhlRPT49sNpt8fPrehG+327Vw4UI988wzNKUAAIAkLt8DAADoU3R0tPVwOByy2WznjX3x8r2MjAzNmzdPa9euVVRUlAYNGqQnnnhC3d3duv/++xUWFqbY2Fj97ne/c1mrsbFRd955pwYPHqzw8HDddtttOnLkyAXzKygocLkkbtOmTQoPD1dnZ6dL3Pz587V06VLr+dtvv62UlBQFBAQoISHByu+cvLw8jRs3TkFBQRo2bJh+9rOf6bPPPrPmX331VQ0aNEhFRUW66qqrZLfb1dDQoN27d+u6665TUFCQBg0apOuvv14NDQ3W6+bOnau33npLZ86cuajjDwAArmw0pQAAAC6znTt36vjx4yotLVVeXp5ycnI0Z84cDR48WPv379eKFSu0YsUKHTt2TJLU3t6uadOmKTg4WKWlpdq7d6+Cg4N18803q6urq9c1mpubVVNTo9TUVGvsjjvuUE9Pj7Zt22aNffrppyoqKtKPfvQjSdJf//pXLV68WCtXrlRtba1efPFFvfrqq1qzZo31Gh8fH61fv141NTX6/e9/r507d+qBBx5wWb+9vV25ubl6+eWXdejQIYWFhWnevHlKS0vTwYMHVVZWpuXLl7tcApmamqqzZ8+qvLz86x9kAADwrUdTCgAA4DILCwvT+vXrNWbMGC1btkxjxoxRe3u7HnroIY0aNUqrV6+Wv7+//v73v0v6fMeTj4+PXn75ZY0bN06JiYnKz8/X0aNHtXv37l7XaGhokDFGMTEx1lhgYKAWLlyo/Px8a2zz5s2KjY3VDTfcIElas2aNVq1apfT0dCUkJGjmzJl66qmn9OKLL1qvyc7O1rRp0xQfH68bb7xRTz31lP74xz+6rH/27Flt2LBBU6ZM0ZgxY9TT06OWlhbNmTNHI0eOVGJiotLT0zV8+HDrNed2UH3ZDjAAAOAduKcUAADAZXb11Ve73F8pKipKSUlJ1vMBAwYoPDxcTU1NkqSKigp98MEHCgkJcXmfjo4Offjhh72uce4SuICAAJfxu+++WxMnTlRjY6OGDh2q/Px8ZWRkWDuWKioq9O6777rsjOrp6VFHR4fa29s1cOBA7dq1S2vXrlVtba1aW1vV3d2tjo4O/fe//1VQUJAkyd/fX+PHj7feIywsTBkZGZo1a5ZmzpypGTNm6Ic//KGGDBnikl9gYKDa29sv7kACAIArGjulAAAALjM/Pz+X5zabrdcxp9MpSXI6nUpJSVFVVZXLo66uTgsXLux1jYiICEmfX8b3v5KTk3XNNddo06ZNqqysVHV1tTIyMqx5p9OpJ554wmWd6upqvf/++woICFBDQ4NuueUWJSUlacuWLaqoqNBzzz0n6fPdUecEBgae9+uE+fn5Kisr05QpU1RYWKjRo0frnXfecYk5deqUvvOd73zZIQQAAF6AnVIAAAAedu2116qwsFCRkZEKDQ29qNeMHDlSoaGhqq2t1ejRo13mfvzjH+s3v/mNGhsbNWPGDA0bNsxlrcOHD+u73/1ur+/7j3/8Q93d3fr1r39t7fb64qV7F5KcnKzk5GStXr1akydP1uuvv65JkyZJkj788EN1dHQoOTn5ot8PAABcudgpBQAA4GGLFi1SRESEbrvtNu3Zs0f19fUqKSlRVlaWPv74415f4+PjoxkzZmjv3r29vl9jY6M2btyoZcuWucw99thj2rRpk3JycnTo0CH961//UmFhoR555BFJnze7uru79cwzz+ijjz7SH/7wB73wwgtf+hnq6+u1evVqlZWVqaGhQdu3b1ddXZ0SExOtmD179ighIUEjR478KocHAABcoWhKAQAAeNjAgQNVWlqq4cOH6/bbb1diYqKWLVumM2fOXHDn1PLly1VQUGBdBnhOaGio5s+fr+DgYM2bN89lbtasWSoqKlJxcbEmTpyoSZMmKS8vTyNGjJAkTZgwQXl5efrVr36lpKQkbd68Wbm5uRf1Gd577z3Nnz9fo0eP1vLly5WZmamf/OQnVswbb7yhu++++yscGQAAcCWzGWOMp5MAAADAV2eM0aRJk5Sdna0FCxa4zM2cOVOJiYlav369h7JzVVNTo+nTp6uurk4Oh8PT6QAAgH6AnVIAAADfUjabTS+99JK6u7utsVOnTqmgoEA7d+7UPffc48HsXB0/flybNm2iIQUAACzslAIAALiCxMXFqbm5WY8++qh+8YtfeDodAACAPtGUAgAAAAAAgNtx+R4AAAAAAADcjqYUAAAAAAAA3I6mFAAAAAAAANyOphQAAAAAAADcjqYUAAAAAAAA3I6mFAAAAAAAANyOphQAAAAAAADcjqYUAAAAAAAA3I6mFAAAAAAAANzu/wDA721B7ESEwQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1200x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib.pylab import plt\n",
    "\n",
    "fig = plt.figure(figsize=(12, 6))\n",
    "\n",
    "# Use log scale on the y axis.\n",
    "plt.yscale(\"log\")\n",
    "\n",
    "# Compute the energy error over the time grid.\n",
    "en_proj_hist = en_cfunc(np.ascontiguousarray(pres_proj[-1].transpose()))\n",
    "err_proj_hist = abs((cb.E0 - en_proj_hist[0]) / cb.E0)\n",
    "err_noproj_hist = en_cfunc(np.ascontiguousarray(pres_noproj[-1].transpose()))\n",
    "err_noproj_hist = abs((cb.E0 - err_noproj_hist[0]) / cb.E0)\n",
    "\n",
    "plt.plot(t_grid, err_proj_hist, marker=\".\", linestyle=\"None\", label=\"With projection\")\n",
    "plt.plot(\n",
    "    t_grid, err_noproj_hist, marker=\".\", linestyle=\"None\", label=\"Without projection\"\n",
    ")\n",
    "\n",
    "plt.xlabel(\"Time (years)\")\n",
    "plt.ylabel(\"Rel. energy error\")\n",
    "plt.legend()\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "234dbb3f-d2b6-4179-90d9-2152a5fee824",
   "metadata": {},
   "source": [
    "We can indeed see how the energy error with projection is kept under the $10^{-6}$ threshold.\n",
    "\n",
    "We can also take a look at the number of manifold projections performed throughout the integration and compare it to the total number of steps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "dd78519a-8fa4-4fab-b3c1-582ab54b0ae9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of projections: 3765\n",
      "Number of steps: 108881\n"
     ]
    }
   ],
   "source": [
    "print(f\"Number of projections: {cb.n_proj}\")\n",
    "print(f\"Number of steps: {pres_proj[3]}\")"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}