{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-00",
   "metadata": {},
   "source": [
    "# L1A — The one-dimensional suite\n",
    "\n",
    "Mixed finite element methods — a crash course, Lecture 1.\n",
    "\n",
    "This notebook contains the live demonstrations of Lecture 1 and the associated\n",
    "self-study material. Setting: the one-dimensional mixed Laplacian on $(0,1)$\n",
    "with $\\kappa=1$ and $u(0)=u(1)=0$, discretized with the pairs\n",
    "$P_1$–$P_0$, $P_1$–$P_1$ and $P_2$–$P_0$. Signs and notation follow the\n",
    "course blueprint: $b(\\tau,v)=(\\tau',v)$ and $G(v)=-(f,v)$. The boundary\n",
    "condition is natural in the mixed formulation; the code therefore contains no\n",
    "`DirichletBC`.\n",
    "\n",
    "**Cell labels.** Each code cell carries a `Cell N` label in its first line;\n",
    "the slides refer to these labels. Tags: `[LECTURE]` cells are run during the\n",
    "lecture; `[ADD-BACK n]` cells are the optional additions of the timing plan;\n",
    "`[SELF-STUDY Sk]` cells contain material omitted from the lecture and used by\n",
    "Exercise Sheet 1.\n",
    "\n",
    "**Execution.** Run the cells from top to bottom. On Google Colab, the first\n",
    "cell installs Firedrake through FEM on Colab (several minutes on first run).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cell-01",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Firedrake available.\n"
     ]
    }
   ],
   "source": [
    "# Cell 1 [SETUP] -- Firedrake via FEM on Colab (no effect on a local installation)\n",
    "try:\n",
    "    import firedrake\n",
    "except ImportError:\n",
    "    !wget \"https://fem-on-colab.github.io/releases/firedrake-install-release-real.sh\" -O \"/tmp/firedrake-install.sh\"\n",
    "    !bash \"/tmp/firedrake-install.sh\"\n",
    "    import firedrake\n",
    "print(\"Firedrake available.\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cell-02",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "setup complete\n"
     ]
    }
   ],
   "source": [
    "# Cell 2 [LECTURE] -- imports, data, and a solver for the mixed pair\n",
    "from firedrake import *\n",
    "import numpy as np\n",
    "import scipy.sparse as sp\n",
    "import scipy.linalg as sla\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Data: UFL expressions (for forms and errors) and NumPy versions (for plots).\n",
    "CASES = {\n",
    "    \"f=1\": dict(\n",
    "        f=lambda x: Constant(1.0),\n",
    "        sig=lambda x: 0.5 - x,          u=lambda x: 0.5 * x * (1 - x),\n",
    "        sig_np=lambda x: 0.5 - x,       u_np=lambda x: 0.5 * x * (1 - x)),\n",
    "    \"f=pi^2 sin(pi x)\": dict(\n",
    "        f=lambda x: pi**2 * sin(pi * x),\n",
    "        sig=lambda x: pi * cos(pi * x), u=lambda x: sin(pi * x),\n",
    "        sig_np=lambda x: np.pi * np.cos(np.pi * x),\n",
    "        u_np=lambda x: np.sin(np.pi * x)),\n",
    "}\n",
    "\n",
    "def build(n, Vdeg_fam, Qdeg_fam, case):\n",
    "    \"\"\"Assemble the mixed problem on a uniform mesh with n elements.\"\"\"\n",
    "    mesh = UnitIntervalMesh(n)\n",
    "    V = FunctionSpace(mesh, *Vdeg_fam)\n",
    "    Q = FunctionSpace(mesh, *Qdeg_fam)\n",
    "    W = V * Q\n",
    "    sigma, u = TrialFunctions(W)\n",
    "    tau, v = TestFunctions(W)\n",
    "    a = (sigma * tau + tau.dx(0) * u + sigma.dx(0) * v) * dx\n",
    "    x, = SpatialCoordinate(mesh)\n",
    "    L = -CASES[case][\"f\"](x) * v * dx          # data sign: G(v) = -(f, v)\n",
    "    return mesh, W, a, L\n",
    "\n",
    "def solve_pair(n, Vdeg_fam, Qdeg_fam, case, solver=\"mumps\"):\n",
    "    mesh, W, a, L = build(n, Vdeg_fam, Qdeg_fam, case)\n",
    "    w = Function(W)\n",
    "    solve(a == L, w, solver_parameters={\n",
    "        \"mat_type\": \"aij\", \"ksp_type\": \"preonly\",\n",
    "        \"pc_type\": \"lu\", \"pc_factor_mat_solver_type\": solver})\n",
    "    return mesh, w\n",
    "\n",
    "def rel_L2(err_expr, ref_expr):\n",
    "    return sqrt(assemble(err_expr**2 * dx)) / sqrt(assemble(ref_expr**2 * dx))\n",
    "\n",
    "def eval_at(f, points):\n",
    "    \"\"\"Evaluate a Function at physical points (input ordering preserved).\n",
    "\n",
    "    Recent Firedrake provides the PointEvaluator class and deprecates\n",
    "    Function.at; older installations fall back to Function.at.\n",
    "    \"\"\"\n",
    "    pts = np.asarray(points, dtype=float)\n",
    "    if pts.ndim == 1:                        # 1D mesh: (n,) -> (n, 1)\n",
    "        pts = pts[:, None]\n",
    "    try:\n",
    "        from firedrake import PointEvaluator\n",
    "    except ImportError:\n",
    "        seq = pts[:, 0].tolist() if pts.shape[1] == 1 else pts.tolist()\n",
    "        return np.asarray(f.at(seq))\n",
    "    pe = PointEvaluator(f.function_space().mesh(), pts)\n",
    "    return np.asarray(pe.evaluate(f))\n",
    "\n",
    "print(\"setup complete\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-03",
   "metadata": {},
   "source": [
    "## The lowest-order pair $P_1$–$P_0$\n",
    "\n",
    "$V_h$ = continuous $P_1$ ($\\dim = n+1$), $Q_h$ = piecewise constants\n",
    "($\\dim = n$). The next cell reproduces the plot and the convergence table of\n",
    "the lecture.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cell-04",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nn= 3\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAEmCAYAAAD/UpNPAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAW8VJREFUeJzt3XdYk+f+P/B3EhI2yJC9ZDpZKoh7UUSROKutttrW2p5qW62tbU9P1eo5p8P2aHftst9WWxdogooL9ygOZIgIEpS9t+wk9+8Pfj42BRQ0IWA+r+vK1cP9rDs5kefNcy8eY4yBEEIIITqLr+0KEEIIIUS7KAwQQgghOo7CACGEEKLjKAwQQgghOo7CACGEEKLjKAwQQgghOo7CACGEEKLjKAwQQgghOk5P2xV4EKVSiYKCApiamoLH42m7OoQQQkivwRhDbW0tHBwcwOd3/Pd/jw8DBQUFcHZ21nY1CCGEkF4rNzcXTk5OHW7v8WHA1NQUQOsbMTMz03JtCCGEkN6jpqYGzs7O3L20Iz0+DNxtGjAzM6MwQAghhDyEBzWzUwdCQgghRMdRGCCEEEJ0HIUBQgghRMdRGCCEEEJ0nM6FgejoaPj5+cHQ0BB+fn6Ijo7WdpUIIYQQrdKpMBAdHY3Zs2cjJSUFjY2NSElJwezZs7Fnzx4wxrRdPUIIIUQrevzQQnX64IMPwOPxuBs/Yww8Hg8rV65ERUUFeDweBAIBRCIRhEIhRCKRyv82MDCAkZERDA0NVf6rr69PsyMSQgjptTQeBrKysnDw4EGUlZXB29sbs2bNgoGBgaYv266MjIw2TwAYYygsKkQLWiBkQsjlcsjl8i6dVyAQwMTEBGZmZjA1NeVeZmZmMDc3h1AoVOfbIIQQQtSKxzT4fHzz5s3YunUrJk6cCDMzM0gkEtTU1ODcuXOwt7fv1Dlqampgbm6O6urqR550yM/PDykpKaqBgAcYOBlgwPoB8FJ6ob+iPxyYA3hQ31/6pqam6NOnDywsLFReFBIIIYRoUmfvoRoNA9evX8eAAQO4R+gNDQ1wcnLCO++8g7feeqtT51BnGLjbZ+BuU8Hd/7q/7g6jACNuPwcjB4Q6hGKCzQRYCCzQ3NyM5uZmNDY2or6+Hg0NDdx/GxoaHqq/AY/HQ58+fWBtbY2+ffvC2toa1tbW0NPTqZYbQgghGtQjwsDf1dbWwtnZGRs3bsSLL77YqWPUGQaA1kCwfv16pKenw8fHB2vXrsXkaZNx5PYRSGVSJJQkqOw/zHYYxJ5ihLqGwlho3OZ8SqUS9fX1qK2tbfOqrq5GfX19p+vG4/FgYWEBOzs77mViYvLI75kQQrrDxo0bkZubiy+++ELbVSH/X48JA0VFRfj0009RV1eHM2fOIDQ0FJ988kmHj8ibmprQ1NTE/Xx3kQV1hYEHyanJgVQmRYwsBgV1BVy5oZ4hJrtMhthTjOF2w8HndW4gRlNTEyorK1FVVYXKykpUVlaioqKi0yHBxMSECwaOjo4wMzOjzoqEkB5pxYoVyMzMxP79+7VdFfL/9ZgwUF5ejq1bt6KqqgpSqRQWFhbYvXs3bGxs2t1/3bp1+OCDD9qUd1cYuEvJlLhcdBkSmQRHs4+iQd7AbbM3tkeEewTEnmK4mrk+1Pnr6+tRVlaG0tJS7r+dCQgmJiZwcnKCk5MTHBwctNYZkxDS89y4cQMbN25Eamoq7Ozs8NRTT2HevHkAgLNnz+K1117Dtm3bMHDgQABAXFwcVq9ejR07dsDLywvr16/n5l6xtrZGSEgI3n77bZUnlDU1Nfjf//6HkydPQl9fH08//TQWLVqEH3/8Ee+++y6am5vRr18/AMCXX36JMWPGdFjfS5cuISkpCeXl5SrNrc8//3yH9wjSNT0mDPxVU1MT/P39MWHCBHzzzTcd7qPNJwPtqW+px9Hso5DKpLhYdFFlm39ff0R6RiLMLQxmokerX319PYqLi1FYWIji4mKUlZU9sD9C37594ezsDFdXV1hbW9NTA0J0VHJyMsaNG4eVK1di0qRJyMnJwVtvvYU333wTK1asAAAsWLAAqampiI+Px507d+Dn54dFixbhP//5D4DWpeLLy8sBtD7V/fDDD6Gvr48jR44AAOrq6hASEgJ9fX289957MDY2xvbt2zF79mwEBQXhjTfeUGkmcHd3b/f3tkKhwLx583D69GmIxWKUlJRAKpWib9++8PDwwNGjR6mJVE16ZBgAWr+M+fn5OHnyZKf2V3efgUeVfycfMbIYSGVS5NbmcuX6An1MdJ6ISM9IhNiHQMAXPPK1mpubUVJSgqKiIhQUFKC4uPi+4cDIyAiurq5wdXWFo6MjBIJHrwMhuo4xBqVSqbXr8/n8ToX8iIgIeHl5YdOmTVzZ7t27sWzZMpSUlABo/X3q5+eHmTNn4vbt28jLy8P58+c77LhcWloKGxsbyGQyuLu7Y/PmzdiwYQOysrJgbm7O7dfc3AyRSNTpZoL169fj+++/x+XLl2FnZwcAmDt3LvLz83H+/PkHvlfSeZ29h2q06/qpU6cwbtw47ueKigqcPn0ac+bM0eRlNcrRxBEv+72Ml3xfQmJpIiSZEhy+fRh3Wu4g9nYsYm/HwsbQBtM8pkHsIYZHH4+HvpZIJOKaBIDWf3CFhYXIz89HXl4eqqqqVPavr69HWloa0tLSIBQK4eTkBHd3d7i4uNAwRkIeklKpxE8//aS167/wwgudCvanTp1Camoqzpw5A8YYGGOor69HaWkpKisrYWFhATMzM2zbtg3jxo2DgYEBrl69qhIEysrKsHnzZly8eBFlZWVQKpXg8XjIysqCu7s7zpw5g7Fjx6oEAaD1d1VnKZVKfPPNN9iwYQMXBABg2LBhuHTpUqfPQ9RLo2Hgm2++wdtvvw0/Pz80NDQgNjYWgwYNwvvvv6/Jy3YLHo+HAJsABNgE4O2gt3E85zikMikuFFxASUMJtl7biq3XtmKw1WBEekZiar+pMNc3f/CJ70MkEnF/+QPAnTt3kJeXh+zsbOTl5UGhUHD7trS04NatW7h16xb09PTg4uICDw8PODs70/BFQh5D9fX1eOGFFxAREdFmm6mpKfe/7z5lMDExgYWFhcp+U6ZMgaWlJZYvXw57e3vo6elh+PDhaGxsBNA6PNzKyuqR6nnlyhUUFxcjNDRUpTwvLw8uLi6PdG7y8DR6V9i5cydSUlIQHx8PoVCI1157DcOGDdPkJbXCUM8Q09ynYZr7NBTXFWN/1n5IZBLcqr6Fa+XXcK38GjZe2ojxzuMR6RGJUY6jIOQ/+l/qJiYm6N+/P/r37w+5XI78/HxkZ2cjOzsbDQ33OjzK5XJkZWUhKysLQqEQrq6u8PT0hJOTE/h8nVqegpDHlpeXFwoKCuDv79/hPrW1tVi4cCGWLVuGP//8E0uWLMG+ffsAtN6Mr1y5gszMTHh4tD7RzMjIUPkjw8fHB3FxcR2en8/nP7CfU35+PgCoBJHm5mZIpVKubwPpft3eZ6Crelqfgc5ijOFa2TVIZBLE3opFTXMNt83SwBLT3FubEXwsfTRy7ZKSEi4A1NXVtbufoaEhvLy84O3tDUtLS7XXg5DHQW/pM/Dtt99ixYoViI6OxrRp0wC0Tge/a9cuvPPOOwCAZ599FklJSbh48SJyc3MREBCAzz77DEuXLkVtbS0sLS3xyy+/YMGCBaivr8fcuXNx8OBBxMTEICIiAqmpqQgMDMSGDRuwevVqAMDJkychEokwcuRIfPLJJ9i2bRuSkpI6rHNSUhL8/f3x+++/46mnnoJSqcTrr7+O2NhYJCQk9Krf871Bj+1A2FW9NQz8VbOiGSdzT0Iik+Bc/jko2L2k3d+yPyI9WpsRrAwf7fFbexhjKC4u5oJBR8MXra2t4e3tDQ8PDxgaGqq9HoQQzfvkk0/w3//+F4aGhuDz+TAyMsLHH3+MWbNmYefOnVi8eDEuXbqEwYMHAwB+/vlnvPbaa0hISIC3tze++uorvPnmm3B2dkZxcTHmzJmD7du3Iyoqimt+2L9/P/7xj3+grq4ORkZG8Pb2xvbt22Fvb4+cnByMHj0aAGBpadnh0MIFCxZgz549GDFiBPLy8sDn87F3716uXkR9KAz0UGUNZTiQdQASmQQ3K29y5Xo8PYx2Gg2xhxjjnMZBKFB/hz+lUomioiJkZmZCJpOhpaWlzT58Ph/9+vXDgAEDYG9vT0MVCell7jYLmpqaqqwBk5mZCT6fD3d3d5X9U1JSYG1tze1bW1uL3Nxc2Nvbw8LCAsnJyXBzc1P5/csYg0wmg6mpKWxtbdtcPycnBzU1NR0OLQRa5z3IysqCo6Mjxo0bR32ZNITCQA/HGMONihuQyqQ4kHUAlU2V3DZzfXNM7TcVYg8xBloN1MgNWS6X4/bt27h58yby8vLabefr06cPBgwYAC8vL5rciBBCeiEKA71Ii6IFZ/LPQJIpwem805Cze0soe/bxRKRHJCLcI9DXqK9Grl9XV4ebN28iIyOjzXBFoHWJZnd3dwwaNIhmBSOEkF6EwkAvVdlYiYO3DkKSKUFaRRpXzufxMdJhJMQeYkxwmQB9gb7ar80YQ1FREdLS0pCVldVupylbW1sMHjwY/fr1o5EIhBDSw1EYeAzcrLzJLZpU3ljOlZsKTTGl3xREekTCr6+fRpoRGhsbkZ6ejrS0NNTU1LTZbmxsjEGDBqF///7UhEAIIT0UhYHHiFwpx/mC85DKpDiecxwtynsd/9zM3BDpEYnpHtNhZ2x3n7M8HMYY8vPzkZqaiuzs7DbbBQIBvL294evr22ZWMkIIIdpFYeAxVd1UjcO3D0OSKUFyWTJXzgMPwfbBiPSIxCSXSTASGqn92jU1Nbh27RrS09PbjETg8Xjo168f/P39YW1trfZrE0II6ToKAzogqzoLMbIYxMhiUFxfzJUb6RnhCbcnIPYQI9A2EHyeetv2m5ubkZGRgWvXrrXbhODo6Ag/Pz84OjrS0ERCCNEiCgM6RKFUIL4oHlKZFHHZcWhUNHLbHE0cIfYQY7rHdDiZOqn1ukqlEjk5OUhMTORWRfsra2trBAYGwtXVlUIBIYRoAYUBHXWn+Q6OZB+BJFOChJIElW1DbYdC7CHGE25PwFhorLZrMsZQWFiIpKQk5ObmttluZWWFoUOHUigghJBuRmGAILcmF9Ks1tEI+XfyuXJDPUNMcpkEsacYQXZBam1GKC8vR1JSEmQyWZuJjCgUENL7ffvtt+Dz+XjppZe0XZXH3oYNGzBgwADMmTPnoc9BYYBwlEyJK8VXIMmU4Ej2ETTI761oaGdsh+nu0xHpEQk3cze1XbOmpgZXr15FRkYGhQJCHhP5+fkYMmQIEhMTVZYbnjdvHsrLW4c/6+vrw8fHBytXroSzs/NDXUcsFnMLrBkYGGDgwIFYuXKlyvTKuiA+Ph4zZszAzZs3YWJi8lDnoDBA2lXfUo9jOccgzZQiviheZZtfXz9EekRiSr8pMBOp57OmUEDI4+Ott95CdnY2du3axZWVl5fD2toaH330EYYOHYqqqiqsW7cODQ0NSEtLg0gk6tI1srOz4ebmhs8//xwDBw5EeXk5/vWvf0EkEiElJUXnJjsbOnQonnvuOSxfvvyhjqcwQB6o4E4BYmQxkMqkyKnN4cpFfBEmukyE2FOMEPsQCPiCR77W/UKBra0tgoODYWen/nkSCNE1JSUlSEhIQEFBgcosoqNGjcKAAQMe+rxKpRK2trb47rvvMHv2bK78yJEjCAsLQ05ODvckICYmBpGRkUhKSoKvr2+XrhMdHY3Zs2ejvLycW1r9jz/+wNNPPw2ZTNZmoaWu2rJlC3JyciAWi/HDDz+gpKQETzzxBJYtW/ZI5/27NWvWwNbWFp6enti2bRvq6urw1FNPYe7cudw+8+fPx3PPPYfCwkLs378fQqEQK1euRFBQELfPxx9/jJ07dyIhIaG9yzxQZ++htEyUDnMwccBLfi9hqe9SJJYmQpIpweHbh3Gn5Q4O3T6EQ7cPoa9hX0S4RyDSIxKeFp4PfS0zMzOMGzcOAQEBbUJBcXExpFIpXFxcEBQUxP0CIIR0zbfffovVq1dj5MiREIlEiI2NhUAggLOzM3766adHCgPJyckoKyvDiBEjVMoTEhJgYWGh0iTQ0NDaFPkwy6EnJCTA2dlZ5ffAo5zv7yQSCdLT05GSkoJnnnkGeXl5eO211+Dp6YmwsDCVfXfv3o0tW7bc93zr16/HyJEj25Tv2LEDADBy5EjMmjUL8fHxmD9/Pvz8/ODt7Y3q6mrs3LkTqampmD59Op599ln89ttvmDZtGoqKiiAQCLjj3333XZSXl8PKSv3L3N9FTwaIikZ5I07knoAkU4ILhRegZPf+shhkNQiRHpGY2m8q+hj0eaTr1NTU4MqVK7h582abbd7e3hg2bNhDt5ERok5ypRxlDWVauba1oTX0+J37m+3s2bOYOHEiJBIJwsPDAbR2QPvkk09QUVEBoVCI27dvY9euXVi9enWX67Jnzx7Mnz8fLS0tKs16c+fORVlZGU6cOAEAqK6uxpQpU9DS0oLLly93+Trh4eHQ09NDTEwMAKCiogKTJk2Cubk5Tp482eXz/Z2joyMCAgKwf/9+rszPzw+LFy/GypUrVfbNyclBRkbGfc/n5+eHvn1VF5Grq6uDmZkZli9fjs8//5wrNzY2xo4dOzB9+nScOXMGY8eOxa+//opnnnkGAJCWloaBAweioKCA6x9xt9kkPj5e5YlBZ9GTAfJQDPQMEN4vHOH9wlFcV4z9WfshlUmRVZ2F1PJUpJanYuPljRjvNB6RHpEY7TQaQr6wy9cxMzPDhAkT4Ovri0uXLiEn514zRUZGBmQyGQYOHIiAgABa+4BoVVlDGUL3hGrl2kfnHO30NOObNm3C3LlzuSAAAGPHjsWaNWtQVlYGe3t7nD9/vt0A3hkNDQ3Q19dv07/nypUraGhowOTJk9HY2IjU1FSMGzcOX331FQBg48aNSEpKAgD8+uuvD2zzT0hIgEgkwuTJk9HQ0IDU1FSEhYXhiy++AAC0tLTgq6++gkwmw9SpUzF16tROv4eysjIUFBTgt99+UynPz8+Hk1PbeVhcXFxUOkp2VkpKCpRKJVatWsWVVVVVob6+nrtOUlIS7OzssGDBApV6iEQilXBx92lIfX19l+vRFRQGSIdsjW3xwpAX8Pzg55Fanop9mfsQeysWNc01OJZzDMdyjsHSwBJT+02F2FOM/pb9u3wNKysrTJkyBYWFhbh48SKKi1tnUlQoFEhJScGNGzcQGBiIwYMHc4/NCCGqFAoFDh061OaRdlFREfT19blH7klJSfD0fLjmvr59+6K+vh6NjY1cQK+srMStW7ewYcMGjBgxAmZmZnB3d1eZknzo0KGwt7fHkiVL8Msvv3Bh4Pvvv1fpiPjrr79CoVCgpKQEn376Kfz8/GBubg4PDw+VJoOlS5eipKQE4eHhWLZsGb755huVAHQ/SUlJ0NPTU3msn5eXh/Lycvj5+bXZ/2GbCZKSkuDq6qoSJBITE6Gnp4eBAwdy+4waNUolHCUmJmLQoEHQ07t3a747SkPTy8dTGCAPxOPxMNh6MAZbD8bq4atxKu8UJJkSnM0/i4rGCmxL24ZtadvgY+GDSI9ITHOfBivDrrVt2dvbIzIyEtnZ2bh48SKqqqoAtP4VEB8fj7S0NIwYMYJGHpBuZ21ojaNzjmrt2p1RVlaG+vr6NjcMiUSCKVOmQF+/dcnzpKQkLFiwAOvXr0dpaSlWr17d6eF/w4YNA4/HQ3JyMve4+sqVKwCAhQsXws3Nrd3jJk6cCAB4+eWXVconTJig0hnQwsIChw8fBgA8++yzbR69A61/HUdHR6OgoADGxsaws7PD999/36Uw8PeVVhMTE2FsbNxuSAoODoaFhcV9z+nl5dXudQICAlTKEhMT0b9/f5X/L8RicZt9/h5KkpOTYW5uDh8fn/u/uUdEYYB0iUggQqhrKEJdQ1HWUIaDWQchkUmQUZmB9Mp0bLy8Ef+78j+McRyDSM9IjHMaB5Ggc0OLeDwe3Nzc4OLigoyMDFy5coUba1xTU4MjR47AwcEBISEhGu1IQ8hf6fH1NLIiqDpZWFjA0NAQJ0+exBNPPAEA2L9/P6Kjo3Hq1Cluv6SkJIhEIkyZMgU3btzA2rVr8fPPP3fqGtbW1ggJCcGxY8e4MJCQkABLS8sOg8D9eHl5tbmRJiQkwMXFpd0gALQ+RndwcICxcesMqoMHD8atW7e47R999BFyc3Px9ddft3t8RzfpIUOGtNt88bDNBElJSW06IyYmJsLf3x9A68iMa9euYd26dW32efHFF1XK4uLiMG3aNI0/Ge2WAZt1dXWora3tjkuRbmRtaI1nBz2LqMgo7J6+GwsHLISFvgUUTIGTeSfxxsk3MHH3RPznz//gWtm1NkMKO8Ln89G/f3/MmzcPQ4cOVflHUFBQgOjoaJw+fZrrYUyIrhOJRHjnnXfw0UcfYerUqQgLC8O8efOwZcsWBAcHA2gdclhTU4OtW7filVdeweLFi7v8lO21117DL7/8wv1bvnLlCoYOHaq29/Gg84lEIjQ1NXE/NzU1qcxjcObMGZSVddzZs6Mw0F4TwcNijCElJeW+18nIyEBDQ4PKPk1NTUhPT1epS0NDA3bv3v3Qcwx0hUZHE0gkEvz3v/9FWloaAMDBwQGbNm3q9CMdgEYT9DYtyhaczTsLiUyCU3mnIFfKuW3u5u4Qe4oR4R4BG6POt3/duXMHFy9eRGZmpkq5UCik/gSE/MWZM2dw+fJlmJubIzw8XGXGvqNHj+Kzzz7DoUOHAADr1q2Dvb19l6YVZowhODgYq1atwrx583Dp0iWYmJh0asiiiYkJqqqqVNrD/+7PP/+ElZVVu4/eAUAul8Pe3h7x8fFwd3fHZ599hrS0NPz4449oamritnV0/KlTpzBw4ECVJw+XLl1C3759H+rpRkd1PHnyJEaMGKEyIurUqVMYMGAAbGxsUFFRgcTERK4JBWi98Z87dw6jRo3iOg1u3LgRFy5cQHR09EPXp0dMOvSPf/wDzz//PJf0NmzYgI8++gjJyckd/p/1dxQGeq/KxkrE3oqFRCbB9fLrXDmfx0eIQwjEHmJMcJ4AA73OjRYoKSnB+fPn26yQaGZmhlGjRj301KeE6IJPP/0UTU1NeO+99wAAkyZNwhdffIFBgwZ16Tx5eXmoqKjo9GRCv/76K44cOYKdO3fiySefxLhx47B06dIu1/+uL774Ah9//DECAgKQkJCAuLg4DBgwAHV1dbhx44Zan1Ro26VLl+Dm5tZhs0ln9Igw8HcKhQKGhob45ptvsGTJkk4dQ2Hg8ZBZmQmpTIqYrBiVMdumQlOE9QuD2EMMv75+D3xsyRiDTCZDfHw815/grn79+iEkJITmJyCkHZcuXYKdnR0Xmvfu3YsZM2ZovEPupUuXkJ6ezv3s4eGBkJCQRzpnamoqsrKyEBwcrPFe9r1djwwDmZmZ8PLywoEDBzo9NpTCwONFrpTjQsEFSGVSHM85jmZlM7fN1cwVkR6RmO4+HfYm91+QRC6XIykpCYmJiVAoFFy5np4eAgMDMWTIEGo6IITovB4XBuRyOZ544gmUlZXhypUrEArbn6imqalJpYNITU0NnJ2dKQw8hqqbqnH49mFIZVIklSZx5TzwEGQXBLGnGJNcJsFIaNThOWpra3H+/HlkZ2erlPfp0wejRo2Co6OjxupPCCE9XY8KA0qlEosXL8bRo0dx5syZ+056sW7dOnzwwQdtyikMPN5uVd/iFk0qri/myo30jPCE2xOI9IjEUNuh4PPaHwCTk5ODc+fOtRm14uHhgREjRnBDkQghRJf0mDDAGMMLL7yA2NhYnDhxAv3733+WOnoyoNsUSgUuFl2EVCbFsexjaFQ0ctscTRxbmxE8psPZtG1nQblcjsTERCQlJak0HQiFQgwfPhwDBw7UueVPCSG6rUeEAcYYlixZggMHDuDEiRMPtWIW9RnQXXea7+Bo9lFIZBJcKb6isi3QJhAzPGcg1DUUJiLVDoM1NTU4d+4ccnNzVcr79u2LsWPH0oRFhBCd0SPCwMsvv4zt27dj165dKkGgT58+6NOnT6fOQWGAAEBubS7XjJB/J58rNxAYYJLrJIg9xAiyC4KA39ppkDGG7OxsnD9/Hnfu3OH25/F48PPzQ2Bg4H3HOxNCyOOgR4QBHx8flUf+d61YsQIrVqzo1DkoDJC/UjIlrhRfgVQmxZHbR1Avv7eSl62RLSI9IhHpEQk3czcA4JZRvXZNdQZEMzMzjB49ut2Vyggh5HHRI8KAOlAYIB2pb6lHXE4cJDIJLhZeBMO9r7JvX1+IPcQIcwuDub45SktLcfr0aW4FsLu8vLwQEhJCyyQT8pg4dOgQlEpll5Y2fpxRGCA6pfBOIWKyWpsRsmvuDTMU8UWY4DIBYg8xgu2CkZaahsuXL6t0MNTX18fIkSPh6elJKyIS0o5PP/0UTz755EMt2tPdFi9eDLlcjm3btmm7Kj0ChQGikxhjSCpNgkQmweFbh1Hbcm+oobWhNSLcIzDJbhIKkguQl5encqyLiwvGjBlDwxAJ+Rsej4ejR49i8uTJ2q7KA1EYUEVhgOi8RnkjTuSegEQmwYWCC1AyJbdtoNVAjDQbCf1b+uA33RtuKBKJEBISAm9vb3pKQHqN27dvIyoqCq+99prKhG4ZGRnYv38/Xn/9dW5GzuTkZJw7dw4mJiYYM2aMygI9O3fuhJOTE2xsbHDixAlYWFigtLQUy5Ytw9KlS+Hl5QULCwu88MILAFrnkDl69Chu3rwJR0dHTJ48GaampgCACxcu4NKlS3jllVe4zrrZ2dnYvXs3Fi5cCDs71WWhy8vLsXXrVjz//POwtLTkygsLC7F9+3YsXboUhoaG+PzzzwG0zjbq5uaG0NBQlQD/9zCwZcsWjB49WmUNht27d8POzg5jxozhyoqKinDkyBE0NDTA39+fW+3xrtTUVFy4cAGGhoYYP378I01oVlZWhsbGe8OmhUIhbG1tH/p899PZeygNuiaPLQM9A4T3C8d3k7/D0TlHsXLoSniYewAArpdfx4+3fsQW/hac6nMKWfwsKKBAc3MzTp06hUOHDqmMQiCkq6Kjo+Hn5wdDQ0P4+fk90spzD2JpaYl//etfiI2NVSnftGkTpFIpFwRef/11hIaG4sqVKzh8+DD8/f3x66+/cvt/++23eOWVVzBt2jQkJCSgsrISxcWtk4BVVlaiqKgIpaWlAICKigqMGDEC77zzDq5fv46vv/4agwYNwo0bNwC09sf56KOPsGHDBgCta9M89dRTOH36dJsgAAAWFhb43//+hx07dqiU//zzz9iyZQvMzMzAGENRURGKiopw69YtfPzxxxg4cCBycnI6/Gw2bNiAS5cuqZRt2bIFBw4c4H6OiopC//79ERMTg4SEBMyZMwfPPfcct/3zzz/HyJEjcebMGRw7dgyhoaHc6o8A8NNPPyEuLq7DOgBAXV0dXn/9dZiZmcHDwwPOzs7ca9GiRfc9tluwHq66upoBYNXV1dquCnkMKJVKdq30Gvv3hX+zUX+MYoN/Gcy9gn4JYgt+XMD+/f2/2ZYtW9jPP//Mrl+/zpRKpbarTXqZqKgoBoDxeDyV/0ZFRWnsmnPnzmVPPvkk93NzczOzsrJiP/zwA2OMsd27dzN7e3tWWlrK7RMbG8tMTExYZWUlY4yxcePGMTs7O1ZVVaVybgDs6NGjKmWLFi1is2fPVvn3sWzZMhYWFsb9fOTIEaanp8fOnj3L1qxZw+zs7FhJSUmH72HlypVs5MiRKmUDBgxga9eu7fCY2bNns6VLl6rUa8GCBdzPjo6ObOvWrSrHTJo0ib399tuMMcaKioqYkZERi4uL47aXlZUxKysrFhMTwxhjzNPTk/3444/c9rq6OpaUlMT97OHhwZYtW9ZhHRljbOHChWzw4MEsIyODMcbYgQMHGAD23XffcZ/JwYMHWV1d3X3P01WdvYfSQGuiU3g8HgZZD8Ig60F4a/hbOJV3CtJMKc7kn0E9q0eSXhKSkAQrpRUGKAeg+kw1tm3bBqlUiszMTHh7e2Pt2rWYNWuWtt8K6cE++OAD8Hg8bjgrYww8Hg/r16/X2Hdn4cKFmD9/Pmpra2Fqaso93ZozZw4AYMeOHXB0dMS2bdvAGANjDAqFAnV1dUhOTsbYsWMBAJGRkTA3N7/vtRQKBfbs2YM5c+bg888/585XX1+PCxcucPuFhobitddew5w5c1BWVob9+/ffdznehQsXYtOmTbh16xb69euHq1evIi0tDQsXLuT2aWlpwbFjx3Dr1i3U19dzM48+rJiYGPB4PCQnJyMpqXWNFMYYrKyscP78eURERMDGxgZHjx5FREQEbG1tYWRkpLKE85IlS+47u+7Vq1exbds2pKSkwMvLCwAwdepUuLi4oLS0lPtMFi1a1GaytO5CYYDoLJFAhFDXUIS6hqKsoQwHsw5CKpMivTId5fxynOWfxcGrB5HzdQ7AA8CAlJQUzJ49G1FRURQISIcyMjJU5rUAWm8wf13KV93Cw8NhZGSE6OhoLFq0CNu3b0dERAQ3wVt+fj4YY206zr7xxhsqs3Le72Z9V1VVFerq6lBfX69yPktLS7z44otQKpXc1N/Lli3Dpk2bMHz4cISFhd33vIGBgRgwYAB+//13vPfee9i+fTuCg4O59WyKioowatQoGBsbIzg4GGZmZqirq0NFRUWnPqP25OfnQ19fv83nMn36dAwdOhQAsHXrVrz55pvo168fPD09MWPGDLz55ptcG/w777xz32tER0fD398fgwcPVilvaWnhllzPz8+HoaEh9PX1H/q9PAoKA4SgdaTBs4OexbODnsWNihvYm7EX0ptS3JTc5IIA0D1/4ZHez9vbGykpKSqBgMfjwcfHR2PXFAqFmDt3LrZv345Zs2ZBKpXi999/57ZbWVlBIBDg008/feRrmZmZQSgUIjQ0FC+++GKH+zHG8OKLLyIoKAhJSUnYvn07FixYcN9zL1iwANu3b8e7776LHTt2qNxof/jhB1hYWODixYtc2Fi3bh1u3brV4fn09PRUhhIDre33d1lZWaGxsREff/xxh8uee3t7QyqVoqmpCWfOnMEbb7yBjIyMNv0bOpKSktJmgb6MjAwUFhZyIzSSkpLg4eGBoqIinD9/HsHBwd266ip1ICTkb/pb9se7I97F6QWnoSxRAn8bb8MYQ1pamnYqR3qFtWvXcsERANdksHbtWo1ed8GCBTh+/Di++eYbGBoaqky8M2PGDMTGxuLatWsqx/w9tLTH1NRU5QYqFAoxdepUfPnllyq94oHW0Qp3ffrpp0hKSsLevXvxySefYNmyZW2WG2/vPdy4cQP/+9//UFxcjHnz5nHb7vaMvxsEGhsbVQJPe1xcXJCSksL9nJeXxzUHAMC0adPQ3NyMr776SuW48vJy7mnB3WYIfX19TJ48GXPnzuU6SgIP7kBoaGjYppPj+++/j9GjR3NPC5KSklBfX4/Fixfj22+/RVBQ0AP/f1ErtfZU0ADqQEi0ydfXl+v8xb14YObO5uzA0QOssbFR21UkPVRUVBTz8/NjBgYGzM/Pj0VHR2v8mkqlkrm5uTGRSMReeukllW0KhYI99dRTzNTUlP3jH/9g77//Pps+fToLCAhgCoWCMdbagfC9995rc95p06axYcOGsQ8//JDrSJebm8u8vb2Zj48Pe+edd9iqVavY0KFDuY55V69eZSKRSOV9h4eHszFjxnDX68ioUaOYSCRi06ZNUym/cuUKE4lEbMGCBWzNmjVs8ODBzMnJiXl4eHD7/L0D4Y4dO5hIJGLLli1j7777Luvfvz+zs7Pj6skYY99//z0TCoVs9uzZbN26dWzRokXM09OTJSYmMsYYGzlyJAsPD2dr1qxhK1euZGZmZuzrr7/mjn9QB0KpVMoAsBdffJH99NNPbNq0acze3p5lZ2dz+8ybN48tX76c+9nNzU0tv186ew+leQYIuY/o6GjMnj37Xmew/99k4PKqC5wDnDFTOBPzxs+jNQ5IjxEVFYULFy5g8eLFbdqoAeDMmTM4deoUGGMIDAxEeHg495f2t99+i379+mHKlCkqx9TW1mL79u3Izs6Gubk59+i+qakJ+/btQ0pKCqysrDBhwgT4+/sDaP1rubGxEcuWLePOU1xcjE8//RSLFi1qt253HTt2DIcOHcLMmTMxatQolW1paWnYu3cvmpqaEBISAnNzc5w6dYqr086dO6FUKvHUU09xx5w7dw7Hjx+HhYUFIiIicOTIETg5Oak8OZHJZJBKpaioqICXlxdmzJjB3XMYYzh8+DAuXrwIQ0NDhIaGcu8TAD766CP0798fM2bM6PA9RUVFYfv27WhoaEBQUBCWL1+u0j9jwIABiI6OxoABA1BZWYlRo0bh+vXrHZ6vs2jSIULUJDo6GuvXr0d6ejq8vLwwKHwQbgy4ATlPDgETYLR8NOb5zENISAithEgI6bKGhgb069cPhYWF4PF42L9/P/bt24cff/zxkc/d2Xso/eYi5AFmzZql0lmQMYa4xDj8O+nfKOeV45TwFPJv5mNG/gyETwyHjY2NFmtLCOltysvLsXLlSq6PiVwux/z587u1DvRkgJCHVFJZgtWHVuNK8xUAgBkzQ7g8HFP8pyAwMJB79EoIIdpC0xETomE2FjbYOn8rXnZ9GUImRA2vBrv1duPnpJ8hkUhQU1Oj7SoSQkinUBgg5BHweDwsG78MP034CTZ8Gyh5SpzRO4NfKn/B9qjt7U4+QwghPQ2FAULUIMA1ANInpRjXZxwAIEuQhd94v2HHqR04fvw4mpqatFxDQgjpGIUBQtTEWN8YX4m/wrtD3oUIItTyahEljMLu27uxZ88eFBUVabuKhBDSLgoDhKjZ04FP44+pf8BRzxFKnhLn9M5hR9MO7IrZhcuXL0OpVGq7ioQQooLCACEa4N3XG/vm7cMU29bJW24LbuMP4R/Yf3U/pFIpdS4khPQoFAYI0RADPQNsnLIRG4ZvgAEMcId3B9HCaBwsP4g9UXtw8+ZNbVeREEIAUBggRONmDJyBPeI9cDNwA+MxXNC7gChE4cCJAzh+/Diam5u1XUVCiI7TeBjIzc3FmjVrEBERgZMnT2r6coT0SK59XBE9JxozXWYCAHL4Odgh2oGTspPYu3cvysrKtFxDQogu02gY+L//+z+MGTMGfD4fBw4c4JaDJEQXCQVCrJ+wHp+O/hSGPEPU8eqwT7gPR+4cQfS+aKSmptKcBIQQrdBoGAgLC4NMJsO6des0eRlCepUwjzDsnbkXXsZeYDyGi3oXsVewF0fOHcHRo0dpTgJCSLfTaBiws7ODQCDQ5CUI6ZUcTR2xc9ZOzHdvXYwkj5+HP0R/4HT2aURFRaG4uFjLNSSE6JIet2phU1OTyl9GNASLPK6EfCHeG/MeQpxD8M8z/0Sdsg4SoQT5DfmoldYieHgw/Pz8uJXMCCFEU3rcaIIPP/wQ5ubm3MvZ2VnbVSJEoya6TcTemXsxwHQAwAMu611GtF40jl88jtjYWDQ0NGi7ioSQx1yPCwPvvvsuqquruVdubq62q0SIxtmb2OP3Gb/jWa9nAQAF/AL8IfoDZwvOIioqCgUFBVquISHkcdbjmgn09fWhr6+v7WoQ0u30+Hp4a+RbCHEOweqTq1GrrEWMMAb5Tfmo3V+L4YHDERgYCD6/x2V4QkgvR79VCOlhRjuPhmS2BIPNBwMAEvQSsFe4F6cSTuHAgQOoq6vTcg0JIY8bjYaBpKQkREREICIiAgDw2WefISIiAl9++aUmL0tIr9fXqC+2RW7D8/2fBw88FPIL8YfoD5wrPoeoqCias4MQolY8psFZTsrKyvDnn3+2KXdxcYGvr2+nzlFTUwNzc3NUV1fDzMxM3VUkpMf7M/9PrDqxCjWK1pE1/nJ/jFSMRNDQIAQGBtJoA0JIhzp7D9VoGFAHCgOEAOUN5Vh5dCWuVl4FANgqbRHWEoZBToMwYcIEGBoaarmGhJCeqLP3UOozQEgvYGVohV+m/4KlA5eCBx6K+cWtaxsUnER0dDSKioq0XUVCSC9GYYCQXoLP4+PV4a/ix9Af0UfQB828ZsQKY3Gw6SD2xuxFcnIyrW1ACHkoFAYI6WWCHIIgnSPFUIuhAIBkQTJ26+3GoT8P0doGhJCHQmGAkF7IwsACP0//Gf8Y+A/wwUcpvxQ7RDtwJOcIoqOjaUlkQkiXUBggpJfi8/h4Zfgr+Cn0J1gILNDCa8Fh4WFIG6SIkkTh+vXr1GxACOkUCgOE9HLDHIZBOkeKYX2GAQCuCa7hD/4fiDkXgxMnTqClpUXLNSSE9HQUBgh5DPQx6IOfI3/GsoHLwAcf5fxy7BTuxP6s/di7dy8qKyu1XUVCSA9GYYCQxwSPx8PLw1/GL0/8AkuBJVp4LTgqPIqoO1HYtXcXMjIytF1FQkgPRWGAkMdMgH0AYubGIKhPEADguuA6fuf9jqhTUThz5gwUCoWWa0gI6WkoDBDyGDLTN8OPkT/i1YGvQgABKvgV2CXchaj0KEgkEty5c0fbVSSE9CAUBgh5TPF4PCwdvhT/98T/wUpgBTlPjjhhHH6v+h1/RP1Bix0RQjgUBgh5zPnZ+yFmbgyCzYMBAOmCdPyq/BW/xv6KhIQEGn5ICKEwQIguMNU3xQ/iH/Cqz6vQY3qo4ldhl3AXfkn4BYcOHaJZCwnRcRQGCNERPB4PS0csxc+TfoY1zxoKngInhCfwQ+EP+CPqD5q1kBAdRmGAEB0T4BwAyVwJgoxbRxvcFNzEj00/4gfJD0hPT9dy7Qgh2kBhgBAdZGZohp/m/IRX3F+BHtNDNb8aOwQ78MWZL3Dq1CnI5XJtV5EQ0o0oDBCiw/4x5h/4fuz3sIY1lDwlTgtPY5NsE3ZJd6G2tlbb1SOEdBMKA4TouOHuw7F39l4MNxgOAJAJZPi65mt8Hf01cnNztVw7Qkh3oDBACEEfkz746cmfsNRpKYRMiBpeDf5gf+DDIx/i8uXLNPyQkMcchQFCCIDW0QavTnoVX4d8jb7oCyVPiTN6Z7A+eT2iD0ajsbFR21UkhGgIhQFCiIoQnxDsEu/CMGHrksi3BLfwaemn+DLqSxp+SMhjqlvCQFZWFk6cOEHtj4T0EtZ9rPHDkz/g+b7PQ8iEuMO7g99afsP70vdxPe26tqtHCFEzjYYBuVyOBQsWwNfXF++88w58fHzw+uuva/KShBA10dPTw8qpK/E////BhtmA8RjOCs5i9YXV2H98Pw0/JOQxotEw8OWXXyI2NhbJycmIj4/H2bNn8d1332Hnzp2avCwhRI3G+4/Hb+G/YRi/tdkgW5CNf+f8G19Gf4mamhot144Qog48psFuwv7+/hgxYgS+++47riwyMhItLS2IjY3t1Dlqampgbm6O6upqmJmZaaqqhJAHaGpqwhexX2Bn1U408ZrAYzyM4o3C2+Pfhpurm7arRwhpR2fvoRp7MiCXy5GamoqAgACV8oCAACQlJXV4XFNTE2pqalRehBDt09fXx5viN/Gf/v+BrdK2tdkAZ/HK8VcQ92cclEqltqtICHlIGgsDd+7cgVwuh6WlpUq5lZUVKisrOzzuww8/hLm5OfdydnbWVBUJIV3E4/EQNiIMP0z6AcPQ2myQy8/FP2/8E19Jv6Lhh4T0UhoLAyKRCADQ0NCgUl5fX89ta8+7776L6upq7kUjEAjpefq59MMXs77A04ZPw4AZoJ5Xjx+qfsDK3StRWFyo7eoRQrpIY2HAyMgINjY2bW7meXl5cHNz6/A4fX19mJmZqbwIIT2PqakpVs9ejX+5/gv2SnuAB5xVnsXzB5/H2cSzNGshIb2IRkcThIWFQSKRcL8UWlpaEBMTgylTpmjysoSQbiIQCCCeIMbmkM0YpmxtNsjj5+GNxDfwbey3NPyQkF5Co6MJbt68ieHDhyMiIgLTpk3D9u3bkZCQgKtXr8LW1rZT56DRBIT0DuXl5fju8HfY17wPjbxGgAGjRaPx7/B/w8rCStvVI0QnaX00AQB4eXnh0qVLsLCwwI4dO9C/f39cunSp00GAENJ7WFlZ4c3Zb+Itm7fgoHRobTZoOYsFkgW4nH5Z29UjhNyHRp8MqAM9GSCkd2GMISExAd9c/QYX+RcBHmDADPCSy0t4fvzz4PNpSRRCukuPeDJACNE9PB4PQwOGYn3YeszlzYURM0IjrxGf536OV3e/itq6Wm1XkRDyNxQGCCEa4ejoiFVzVuFV81fhpHQCAJxuPI05e+Yg+XaylmtHCPkrCgOEEI0xNjbGU5FP4Z9e/0SwPBg8xkMBCvD8yefxy9lfaPghIT0EhQFCiEYJBAKMGT0G74x7B3OUc2DMjNHEa8Jnss+wPHo56hvrtV1FQnQehQFCSLfw9PTE8hnLsdRgKVyULgCA03dOY+bumbief13LtSNEt1EYIIR0G0tLSyyctRCvOb6GEHlIa7OBsgDPHH0Gv8X/pu3qEaKzKAwQQrqVSCTCE6FPYPnw5Zgtnw0TZoJmXjM+ufEJXpW8ioaWhgefhBCiVhQGCCHdjsfjwdfXF0unLcVz/OfgpnADAJysOgnxTjFuFN/QbgUJ0TEUBgghWmNvb4+FcxZiidUSjJaPBp/xUagoxNOHnsbvCb9ru3qE6AwKA4QQrTIyMsL06dOxaNAizG6ZDVNmiha04MOUD/Ha/teo2YCQbkBhgBCidXw+HyNGjMDCSQvxDHsG7gp3AMCJ8hMQ7xIjvSxdyzUk5PFGYYAQ0mO4u7vjqZlPYYHJAoxtGdvabCAvxFMHnsIfyX9ou3qEPLYoDBBCepQ+ffpg5syZmNVvFua2zIUZM0MLWvDfq//F67Gvo76FJikiRN0oDBBCehyhUIiJEydiRsgMPNXyFDwVngCA4yXHMWP3DGo2IETNKAwQQnokHo+HwYMHY27kXMwWzcb4lvEQMAEKW1qbDbYnb6e1DQhREwoDhJAezdbWFrNnz8YU+ymY2zIXfZR90IIWfHT1I7x66FXUtdRpu4qE9HoUBgghPZ6hoSHCw8MxJXAK5rXMg7fCGwBwquQUxLvFuF5GaxsQ8igoDBBCegU+n4+hQ4di5rSZEOuJMbFlIvSYHopbirHgwAL8mvQrNRsQ8pAoDBBCehVHR0fMmTMHk20n48mWJ2GhtIAccmxM3IhlsctQ21yr7SoS0utQGCCE9DpGRkaYNm0aJvtPxryWeRigGAAAOFN6BuI9YlwrvablGhLSu1AYIIT0Snw+H8OHD0dkeCQi9CIwuWUy9JgeSltKsfDgQmxN3ErNBoR0EoUBQkiv5uzsjFmzZmGCzQTMa5kHK6UVFFDgf0n/w8uxL6O6qVrbVSSkx9N4GDhx4gSefPJJeHp6QiKRaPpyhBAdZGJigoiICEzwnYAnW57EIMUgAMD50vOYETUDicWJ2q0gIT2cRsPApk2b8MEHH2DWrFmQyWSoraWOPYQQzeDz+QgODsa0sGkIF4QjrCUMQiZEWUsZFh1ahB+u/kDNBoR0gMc0+K+jubkZIpGo9UI8Hn777TcsXLiwS+eoqamBubk5qqurYWZmpolqEkIeM3fu3MGxY8eQXpqOQ3qHUMYvAwAEWQXhs8mfoY9BH+1WkJBu0tl7qEafDNwNAoQQ0p1MTEwwffp0jB08FnNb5mKIYggA4GL5RYijxLhceFnLNSSkZ9HTdgX+rqmpCU1NTdzPNTU1WqwNIaS3EggECAkJgb29PYxPGcOpxQlxenGokFfghSMvYOnApfjHsH+Az6N+1IR0KQxs2rQJX3/99X332b59O4KDgx+6Qh9++CE++OCDhz6eEEL+ys3NDdbW1ugT1wd9i/vikPAQSvgl+O76d4gviMfmsM2wNLDUdjUJ0aou9RmoqKhARUXFffdxdHSEoaFh2wt1ss9Ae08GnJ2dqc8AIeSRKJVKJCQk4FLCJZwXnEeiXiIAwFxgjk/Hf4oRTiO0W0FCNKCzfQa69GTA0tISlpaaTdD6+vrQ19fX6DUIIbqHz+dj2LBhcHBwgOlxUzg2OuKY3jFUK6qxNG4pFnsvxuvBr0PAF2i7qoR0O2osI4ToFAcHB8yZMwfjncZjfvN82CntwMCwNWMrFu5biNL6Um1XkZBup9EwcP78eXh6esLT0xMA8Oabb8LT0xP//Oc/NXlZQgi5LwMDA4SFhSEsJAxzFHMQKA8EAFyrvQZxlBins09ruYaEdC+NzjPQ0NCA/Pz8NuXm5ubo27dvp85B8wwQQjSprKwMcXFxSKxNxDG9Y2jkNQIAnu73NFaPXk3NBqRX6+w9VKNhQB0oDBBCNK25uRnnzp3D1ZtXcUh4CIX8QgCAj5EPvgr/CnYmdlquISEPp0dMOkQIIb2BSCTChAkTEDE+Ak+yJzFMPgxgQHp9OsTRYhzOOKztKhKiUT1u0iFCCNEWb29v2NrawjLOEo7ljjgiPIJ61OPNC2/i7O2zWDNpDYQCobarSYja0ZMBQgj5C3Nzc4jFYkz3nY75zfPhqHQEAOwr3IdZO2fhdsVt7VaQEA2gMEAIIX8jEAgQHByM+RHzsVC0EEHyIIABt1tuY27MXOy5ukfbVSREragDISGE3EdTUxPOnTuHE7ITrc0GvHoAwOQ+k/HfKf+FoX7bGVcJ6SmoAyEhhKiBvr4+Jk6ciEUTFuFZPAtnpTMA4FjVMczcNRPXcq5puYaEPDoKA4QQ0gmenp54dvazeMnqJYyQjwCP8ZCvzMfi44vx06mfoFQqtV1FQh4ahQFCCOkkU1NTTI+YjmVDl2G2YjaMmTGaeE3YfHszXt71Msoqy7RdRUIeCoUBQgjpAj6fD39/fywTL8NLhi/BVeEKALjQdAHzJPMQdyUOPbwrFiFtUBgghJCHYG1tjWdmP4O3vd/GSPlI8BgPJbwSvJXyFha+vRBDhgyBoaEh/Pz8EB0dre3qEnJfNJqAEEIeUUFBAbYd34a98r3Iv5KP3K9yAR4ABvB4PDDGEBUVhVmzZmm7qkTH0GgCQgjpJg4ODljx5AqscV2Dqn1VXBAAAMYYeDwe1q9fr80qEnJfFAYIIUQNRCIRwieEo7m0mQsCdzHGkHI9BT+l/ISS+hLtVJCQ+6AwQAghauTj4wMej6dayANEdiJsTtiM0D2hePnYy4i9FYtGeaN2KknI39BCRYQQokZr167F7Nmzub4Cd/87YvoINLAG1KEO5/LP4Vz+OZgKTRHWLwxiDzH8+vq1DRGEdBPqQEgIIWoWHR2N9evXIz09HV5eXpg1axYcHByghBI5/Bzc4N/ALcEtyCHnjnE1c0WkRySmu0+HvYm9FmtPHiedvYdSGCCEkG4gk8lw7tw5NDa2Ng00ohG3hLeQbZqNm3U3uf144CHIPghiDzEmuUyCkdBIW1UmjwEKA4QQ0sM0Njbi/PnzyMzMVCkX9BWg3LYcR/KPoLi+mCs30jPCE25PINIjEkNth4LPo25epGsoDBBCSA+Vk5ODM2fOoK6ujisTCAQICAxAo00jYrJicCz7GBoV9zoYOpo4tjYjeEyHs6mzNqpNeiEKA4QQ0oM1Nzfj4sWLuH79ukq5lZUVxo4dC0NzQxzNPgqJTIIrxVdU9gm0CcQMzxkIdQ2FicikO6tNehkKA4QQ0gsUFhbi9OnTqK6uVikfNGgQhg8fDpFIhNzaXMTIYiCVSZF/J5/bx0BggEmukyD2ECPILggCvqC7q096uB4TBnJzc3Hu3Dm0tLRg2LBhGDBgQJeOpzBACHncyeVyJCQkICkpSWWRIyMjI4SEhMDd3R08Hg9KpkRCcQIkMgmO3D6Cenk9t6+tkS0iPSIR6REJN3M3LbwL0hP1iDDw8ssv4+jRoxg+fDgEAgH27duHF198EZs3b+70OSgMEEJ0RXl5Oc6cOYOSEtVZCp2dnTFq1CiV34H1LfWIy4mDRCbBxcKLYH+Z9tC3ry/EHmKEuYXBXN+82+pPep4eEQYkEgkiIiIgELQ+ujp79izGjBmDkydPYty4cZ06B4UBQoguYYwhLS0NFy9eRHNzM1cuEAgQEBAAPz8/7nfqXYV3ChGT1dqMkF2TzZWL+CJMcJkAsYcYIQ4h0OPTPHO6pkeEgfbo6+vjyy+/xNKlSzu1P4UBQoguqq+vx59//tlmGGKfPn0wevRoODg4tDmGMYak0iRIZBIcvnUYtS213DZrQ2tEuEcg0iMSXhZeGq8/6Rl6ZBg4ePAgpk2bhitXriAwMLDdfZqamtDU1MT9XFNTA2dnZwoDhBCdlJ+fj7Nnz7bpYOju7o4RI0bAxKT90QSN8kaczD2JfbJ9uFBwAUqm5LYNtBqISI9ITO03FRYGFpqsPtEyjYSBCxcuID4+/r77zJkzB05OTm3Kc3NzERwcjCeeeAK//PJLh8evW7cOH3zwQZtyCgOEEF0ll8uRlJSEq1evQqm8d1O/23Tg6+sLPb2OmwBK6ktwIOsAJJkSyKplXLkeXw/jnMYh0iMSY5zGQMgXavR9kO6nkTBw8OBBHDly5L77LF++HJ6eniplRUVFGDduHNzd3bFv3z7o6+t3eDw9GSCEkPZVVVXhwoULyM3NVSk3NTVFSEgIXF1d77vYEWMM18uvQyKT4OCtg6huuve0wdLAElP7TYXYU4z+lv019h5I9+oxzQRFRUWYMGECXF1dsW/fPhgYGHTpeOozQAgh9zDGkJOTg/Pnz6O2tlZlm5OTE0aOHIk+ffo88DzNimaczjsNSaYEZ/LPQMEU3DZvC29EekRimvs0WBtaq/stkG7UI8JAcXExJkyYABcXl4cKAgCFAUIIaY9cLkdKSgquXr0Kufze6oc8Hg9DhgxBQEDAfZ/C/lV5QzkO3joISaYE6ZXpXLmAJ8Box9GI9IjEeOfxEAlEan8fRLN6RBjw9/eHTCbDe++9pxIERowYgREjRnTqHBQGCCGkY3fu3EF8fDxkMplKub6+PoYOHYoBAwa0GYp4P+kV6ZDIJDiQdQAVjRVcuZnIDOH9wiH2EGOw9eD7NkeQnqNHhIG33noLLS0tbcqnTJmCKVOmdOocFAYIIeTBCgsLce7cOVRUVKiUm5ubIygoCG5ubl26gbcoW3Au/xykMilO5J6AXHnv6YO7uTsiPSIR4R4BW2Nbtb0Hon49IgyoA4UBQgjpHKVSiRs3buDy5ctobGxU2WZra4uQkBDY2Nh0+bxVjVWIvR0LSaYEqeWpXDmfx0eIfQgiPSIx0WUiDPS63hRMNIvCACGE6Kjm5mYkJSUhOTkZCoVCZZu7uzuCgoIe+veprEoGiUyC/bL9KG0o5cpNhCYIcwuD2FMM/77+1IzQQ1AYIIQQHXfnzh1cvnwZGRkZKuV8Ph+DBg2Cv78/DA0NH+rccqUcfxb+CWmmFHE5cWhW3ps62cXUBZEekZjuMR0OJm1nSiTdh8IAIYQQAEBZWRni4+ORn5+vUq6np4chQ4bA19e30yMP2lPTXIPDtw9DmilFYmmiyrYguyCIPcWY7DIZRkKjh74GeTgUBgghhHAYY8jNzUV8fDwqKytVtolEIvj6+mLIkCEQCh9tFsLb1bchlUkRkxWDoroirtxQzxChrqGY4TkDQ22Hgs/jP9J1SOdQGCCEENKGUqlEZmYmrly50mbSIgMDA/j7+2PgwIH3nd64U9dhSlwsughpphTHco6hQd7AbXM0ccR0j+mIdI+Es5nzI12H3B+FAUIIIR1SKBRIT09HQkIC6uvrVbYZGRkhMDAQPj4+XZqjoCN1LXU4cvsIpDIpLhdfVtkWaBMIsacYT7g+ARNR+4sukYdHYYAQQsgDyeVyXL9+HYmJiW2GIxobG8PX1xf9+/d/5OaDu/Jq8xAji4FEJkH+nXt9GAwEBpjoMhFiTzGC7YIh4D96CCEUBgghhHRBc3Mzrl27huTkZDQ3N6tsMzAwwJAhQzBo0CCIROqZkljJlEgoToBUJsWR7COoa6njttkY2WC6+3REekbC3dxdLdfTVRQGCCGEdFlTUxOSk5Nx7dq1NjPICoVCDBo0CEOGDHnoIYntaZA3IC4nDpJMCeIL48Fw77bka+0LsacYYW5hMNc3V9s1dQWFAUIIIQ+tubkZqampSElJadN8oKenh/79+8PX1xcmJupt5y+qK0KMLAZSmRS3a25z5SK+COOdx0PsKcZIh5HQ4z9aB0ddQWGAEELII5PL5UhLS0NycjLq6upUtvF4PHh4eMDX1xfW1upd6pgxhuSyZEgyJTh06xBqW+6NfLA2tMa0ftMQ6RkJbwtvtV73cUNhgBBCiNooFApkZGQgMTGxzZBEALC3t4evry9cXFzUPhVxk6IJJ3JOQCKT4HzBeSiZkts2wHIAxJ5iTO03FRYGFmq97uOAwgAhhBC1UyqVkMlkSEpKarNCItC6SuLgwYPh5eWlts6Gf1VaX4oDWQcgkUmQWZXJlevx9TDWcSwiPSMx1nEshAL1jH7o7SgMEEII0RjGGAoKCpCcnIzc3Nw224VCIby8vDBw4EBYWlpq5PrXK65DkinBwVsHUd1UzW2z0LfAVPepiPSIxADLATq9aBKFAUIIId2isrISKSkpuHnzZptVEoHWJoSBAweiX79+4PPVPw1xi6IFp/NOY59sH87mnYWcybltXhZeEHuIMc19GqwN1duvoTegMEAIIaRbNTQ0IC0tDWlpaW06GwKAoaEhfHx84OPjA3NzzQwTLG8oR+ytWEhkEtyouMGVC3gCjHIchUiPSIx3Hg99wcMvzNSbUBgghBCiFUqlEtnZ2bh+/XqblRLvsrOzg4+PD9zd3dU2u+HfpVekQyqTYn/WflQ03uvfYCoyxdR+rc0IQ6yHPNbNCBQGCCGEaF1VVRWuX7+OjIyMNjMbAq19C9zd3eHt7Q07OzuN3JhblC04n38eEpkEJ3NPokV5bzKlfub9EOkRienu02FrbKv2a2sbhQFCCCE9hlwuR1ZWFtLT01FYWNjuPiYmJvDw8ICnpycsLS01Egyqm6oReysWUpkUKWUpXDkPPIQ4hCDSIxITXSbCUE99MyxqE4UBQgghPVJ1dTUyMjKQkZHRbt8CALCwsOCCgaZ+92dVZUEik2C/bD9KGkq4cmOhMaa4TUGkRyQCbAJ6dTMChQFCCCE9mlKpRH5+PtLT05Gdnd3uSAQAsLKygpubG/r16wcLCwu135wVSgX+LPwTEpkEx3OOo0nRxG1zNnVGpEckIj0i4WDioNbrdgcKA4QQQnqN5uZmZGdnIzMzE3l5eejo1mRmZgY3Nze4ubnBxsZG7UMVa5trcfj2YUhlUlwtuaqybbjdcIg9xAh1DYWR0Eit19WUHhEG8vPz8d///hdSqRRlZWXw8vLCqlWrsGjRok6fg8IAIYTolsbGRmRlZSEzMxNFRUUd7mdgYAAnJyc4OTnB2dlZrSspAkB2TTakMiliZDEorLvXz8FQzxChrqEQe4gxzG4Y+Dz1z52gLj0iDHzwwQdwc3PDlClTYGZmhh07dmDJkiXYv38/wsPDO3UOCgOEEKK76uvrkZ2djdu3byM/Px9KpbLDfa2treHs7AwnJyfY2NhAIBCopQ5KpsSlokuQyqQ4mn0UDfIGbpuDsQOme0xHpEckXMxc1HI9deoRYaA9Li4ueOGFF7B27dpO7U9hgBBCCNDalJCTk4Pbt28jNzcXLS0tHe4rEAhgY2MDBwcH2Nvbw8bGBnp6j77scV1LHY5mH4VUJsWloksq2wJsAiD2EOMJtydgKjJ95GupQ48KA3K5HHV1ddizZw9WrFiB06dPIyAgoFPHUhgghBDydwqFAsXFxcjNzUVeXh7Ky8vvu79AIEDfvn1hY2PDvYyNjR+pM2L+nXxIZVJIM6XIu5PHlesL9DHRZSJmeMxAsH0wBHz1PKF4GBoJA0ql8r6PaIDWD/yvH25qair8/PygUChgYGCAH374AQsXLuzw+KamJjQ13evJWVNTA2dnZwoDhBBCOlRXV4e8vDzk5uaisLAQDQ0NDzzGyMiICwjW1tawtLSEkZFRlwMCYwwJJQmQyqQ4fPsw6lruDZe0MbJBhHsExB5iuPdx7/L7elQaCQNr167Ff/7zn/vuExcXh3HjxrUpv3PnDnbt2oWXX34Zu3btwowZM9o9ft26dfjggw/alFMYIIQQ0hmMMVRXV6OgoACFhYUoKCjoVDgAAH19fVhaWnIvCwsLmJmZwdDQsFMhoUHegLicOEgzpfiz8E8w3LvFDrEegkiPSIT3C4e5vmbWZvi7HtVM8FczZ85EU1MTDh482O52ejJACCFEne6Gg+LiYpSUlKC0tBTl5eUdDl9sj1AohJmZGczMzGBubg4zMzMYGxvD2NgYRkZG0NfXbxMWiuqKsD9rPySZEtyuuX3vXHwhxjuPh9hDjFGOo6DH10N0dDQ++OADZGRkwNvbG2vXrsWsWbMe+b332DAwdepU8Pl87N+/v1P7U58BQggh6iaXy1FWVsaFg4qKClRVVXUpIPyVQCCAkZERjI2NYWBgAH19fe4lEomQ05KDk6UncbrkNO7I73DHWepbwinLCb+/9zt4PB4YY9x/o6KiHjkQdPYe+uhdK+9j5syZeP311+Hn54eGhgZs374dhw8fxu7duzV5WUIIIeS+9PT0YGdnBzs7O65MoVCgqqoKFRUV3Ku6uhq1tbUPDAkKhQK1tbWora3tcB8PeMAVrrjFv4U0fhpy+DmoaKpA/LfxAA/cNe4GgvXr16vl6UBnaDQMrF69Gv/+978RHx8PoVCIQYMG4cCBA5gyZYomL0sIIYR0mUAggJWVFaysrFTKlUolamtrUVNTg+rqatTU1KC2thZ1dXWoq6vrdH8EANCDHryUXvBSeqEOdUgXpCO1KBX4W9ZgjCE9PV0db6uT9dKgkJAQHDhwQJOXIIQQQjSKz+fD3Nwc5ubmcHZ2brNdqVSivr4e9fX1qKur4/q+/f3V0tIChULBvYwURrBUWOKA3QHk5eepBAIejwcfH59ue48aDQOEEELI447P58PExAQmJiYPdbypqSlmz57dps9AZyfnU4eeO6EyIYQQogNmzZqFqKgo+Pr6wsDAAL6+voiOjsbMmTO7rQ60aiEhhBDymOrsPZSeDBBCCCE6jsIAIYQQouMoDBBCCCE6jsIAIYQQouMoDBBCCCE6rsfPM3B3sENNTY2Wa0IIIYT0LnfvnQ8aONjjw8DdeZ7bm/WJEEIIIQ9WW1sLc/OOl03u8fMMKJVKFBQUwNTUtFNrSXfG3WWRc3Nzae4CNaDPU/3oM1Uv+jzVjz5T9dPEZ8oYQ21tLRwcHMDnd9wzoMc/GeDz+XByctLIue+uTU3Ugz5P9aPPVL3o81Q/+kzVT92f6f2eCNxFHQgJIYQQHUdhgBBCCNFxOhkG9PX1sXbtWujr62u7Ko8F+jzVjz5T9aLPU/3oM1U/bX6mPb4DISGEEEI0SyefDBBCCCHkHgoDhBBCiI6jMEAIIYTouMcyDBw/fhxPPvkkxo8fjxUrVqCkpEQjx+iK+vp6rF+/HhMnTsT06dPx+++/P/CYY8eO4YUXXsCkSZPw/PPPIz4+vhtq2ntcvHgRCxYswLhx4/Dyyy8jOzu708dKJBKMGDECGzZs0GANe5eWlhZ89tlnCA0NRXh4OLZs2fLA6VcBICUlBS+99BImTJiAV199FUVFRd1Q297h+vXreO655zBu3Dg899xzuH79+gOP2bdvH+bPn4/x48dj/vz52Lt3bzfUtHdobm7GH3/8gSlTpiA0NLRTxyiVSnz77bfcMZs2bYJcLtdI/R67MHD48GGEhYVh0KBBeOutt5CSkoLRo0ejvr5ercfoklmzZmHnzp1Yvnw5pk6diiVLlmDz5s0d7v/ZZ5/h448/xqhRo/Dee+/B3t4eI0eOhFQq7b5K92CXL1/G2LFjYWtri3feeQclJSUICQlBaWnpA4/NycnB8uXLUVJSAplM1g217R2WLFmCzz//HEuWLMH8+fPx7rvv4t13373vMbGxsRg+fDhMTEywZs0aBAYGYvHixd1T4R4uMzMTI0eOBI/HwzvvvAMAGDly5H2/c99//z3mzZuHESNGYN26dQgKCsK8efPw/fffd1e1e7QRI0ZAIpHA0dERly5d6tQxb731FtasWYMFCxbghRdewKeffoqXX35ZMxVkj5nAwEC2aNEi7ufq6mpmZGTEvvrqK7UeoyuOHz/OALCUlBSu7MMPP2Tm5uassbGx3WNqa2vblM2dO5dNmDBBY/XsTaZOncrCwsK4n5ubm5m9vT3717/+dd/j5HI5GzVqFPv+++/ZuHHjVL6zuuz69esMADt69ChXtnXrViYUCllpaWm7xzQ0NDBbW1u2YsWKNuWEseeff575+/szpVLJGGNMqVSyIUOGsCVLlnR4THh4OJszZ45K2ezZs1l4eLhG69pbVFVVMcYY+/LLL5m5ufkD9y8sLGQCgYD98ccfXFlMTAzj8Xjs5s2baq/fY/VkoLKyEgkJCZg+fTpXZmZmhvHjx+PYsWNqO0aXxMXFwdXVFYMHD+bKxGIxqqurcfny5XaPMTExabesublZY/XsLZRKJU6cOKHyfRMKhQgPD3/g923t2rWwtrbGiy++qOlq9ipxcXEwNjbGhAkTuDKxWIyWlhacOnWq3WOOHj2K4uJivPLKKyrlBgYGGq1rbxEXF4eIiAhuPRgej4fp06ff9zs6bNgwpKSkcIvL1dTUICUlBUFBQd1S556uM1MC/9XJkyehUCgQERHBlYWFhUEkEiEuLk7d1Xu8mglycnIAAA4ODirlDg4OHbbJPswxuiQ7OxuOjo4qZXc/q85+PpmZmdi5cydmzJih7ur1OqWlpWhoaOjy9+348eP45Zdf8OOPP2q6ir1OdnY2bG1tIRAIuDILCwsYGhp2+Jlev34dZmZmqKysxJw5cxAaGopVq1ahsLCwu6rdYzHGkJub2+539O7vy/asXbsWYrEYzs7OCAgIgIuLC2bNmoU1a9ZousqPpezsbJiZman8cSUUCtG3b1+N3Jt6/EJFXdHS0gIAbWZvMjQ05Lap4xhd0tLS0u5nc3fbg5SXlyMyMhIjRozAihUrNFHFXuVhvm+lpaV45pln8PPPP8Pa2lrjdext2vuOAq1/5Xf0mTY2NqK5uRnPP/883n//fZibm2Pz5s3cX7eWlpaarnaPpVAooFQq2/2OKpVKKBQKleB11759+/Ddd9/h/fffR2BgIBISErBhwwYEBQVh5syZ3VX9x0ZH32tN3ZseqzBw9x9wRUWFSnl5eTmsrKzUdowusbS0xI0bN1TKysvLAeCBn09FRQVCQ0NhZWUFiUQCPb3H6uv2UCwsLMDj8br0fTty5AgqKiqwdu1arF27FkDrX7apqakYMWIEDh8+3OVHkI8TS0vLNp+nQqFAVVXVff/dNzY24vPPP8ekSZMAAGPGjIG1tTX27NmDpUuXarzePZWenh5MTU3b/Y6am5u3GwQA4M0338SLL76IN998EwAwceJEFBUVYdWqVRQGHoKlpSUqKyvBGOOaawDN3Zseq2YCNzc3WFpatumpefHiRQQEBKjtGF0SGBiIGzdu4M6dO1zZ3WGC/v7+HR5XWVmJ0NBQGBkZITY2tt1+BLrI2NgY3t7ebb5v8fHxHX7fwsLCcOLECWzevJl7eXp6IigoCJs3b4axsXF3VL3HCgwMRHFxscoj7EuXLoEx1uFnOmzYMACqzYPGxsYwNTVFZWWlZivcCwQGBnbpOwq0hv/2mhT/HipI5wQGBkIulyMxMZErk8lkqKio0My9Se1dErXsjTfeYK6urqyoqIgxxti2bdsYn89nSUlJ3D5r165lixcv7tIxuqqiooL16dOHvffee4yx1t7WISEhKr3hy8vLWXBwMNebu6qqig0dOpSNGjWq3ZEFuu6TTz5hVlZWLDMzkzHG2JEjRxifz1fpDf/FF1+w6dOnd3gOGk1wT2NjI3NxceF6usvlchYREcH8/f1V9gsODmY7duxgjLX2jvf19WXLly/neszv2rWL8fl8dvHixe59Az3Qr7/+ygwNDdmVK1cYY4xdvnyZGRgYsN9++43bZ/v27WzUqFHcz3c/84qKCsZY6+8Of39/FhER0b2V7+HuN5ogLCyMfffdd4yx1u/o4MGD2axZs5hCoWCMMfbMM8+wfv36sebmZrXX67ELA/X19UwsFjNDQ0Pm6enJjI2N2Q8//KCyz4IFC9jQoUO7dIwuO3r0KLOxsWEuLi7M3NycDR8+nOXn53PbCwsLGQBuCMxbb73FALDBgwez4OBg7jV16lRtvYUeRS6Xs8WLFzN9fX3m7e3NDAwM2EcffaSyz6pVq5ijo2OH56AwoOrSpUvM1dWV2dvbM2trazZgwACWnp6usg8AtmnTJu7njIwMNnjwYObo6Mh8fHyYubk527JlSzfXvOdatWoVE4lEzNvbm4lEIrZq1SqV7Rs3bmQCgYD7OTc3l40dO5aZmJgwX19fZmJiwsaNG8dyc3O7u+o90urVq1lwcDBzc3NjAoGA+72YkZHB7WNlZcX94cUYY2lpaczHx4f17duX2dnZMTc3N5aQkKCR+j22qxbm5eWhtLQU3t7ebR6jymQyNDQ0qAyXe9Axuq6lpQVpaWkwMjKCp6dnm21XrlyBl5cXrKyskJ2d3W6vbJFIhMDAwO6qco9XXFyMgoICuLu7t2nzz8nJQXl5eYePA69fvw59fX14eHh0R1V7BYVCgbS0NOjp6cHHx0elnRUA/vzzT7i5ucHOzk6lPCMjAwqFAu7u7rQc79+Ul5cjOzsbrq6ubdqpCwsLkZOTg+DgYJXykpISFBQUwMHBATY2Nt1Z3R4tIyOj3SYTX19fGBkZAQCuXLkCW1tbODk5cdsZY7hx4waUSiUGDBgAPl8zrfuPbRgghBBCSOc8Vh0ICSGEENJ1FAYIIYQQHUdhgBBCCNFxFAYIIYQQHUdhgBBCCNFxFAYIIYQQHUdhgBBCCNFxFAYIIYQQHUdhgBBCCNFxFAYIIV22Zs0arFixQqXszTffxLJly6BUKrVTKULIQ6PpiAkhXRYfH4+QkBBcvnwZgYGB+Ne//oX/+7//w/nz5+Hs7Kzt6hFCuojCACHkocyYMQNyuRyRkZF4++23cebMmTaLfxFCegcKA4SQh5Kamgo/Pz8YGBjg4MGDGDt2rLarRAh5SNRngBDyUGQyGQDAx8eHggAhvRyFAUJIl124cAFPP/00Nm7ciNTUVEgkEm1XiRDyCKiZgBDSJenp6Rg1ahTee+89rFy5EsuXL8fJkyeRnJwMPp/+viCkN6IwQAjptKKiIoSEhGDWrFn47LPPAAAFBQXw8PDAli1b8Oyzz2q5hoSQh0FhgBDSabdu3UJubi7GjBkDHo/HlScmJkJPT49GExDSS1EYIIQQQnQcNfARQgghOo7CACGEEKLjKAwQQgghOo7CACGEEKLjKAwQQgghOo7CACGEEKLjKAwQQgghOo7CACGEEKLjKAwQQgghOo7CACGEEKLjKAwQQgghOo7CACGEEKLj/h8effqHa1RITgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 600x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    n  rel sig err   rate    rel u err   rate  max mid err   rate\n",
      "    8    1.404e-02           1.138e-01           1.872e-02\n",
      "   16    3.517e-03   2.00    5.674e-02   1.00    4.784e-03   1.97\n",
      "   32    8.797e-04   2.00    2.835e-02   1.00    1.203e-03   1.99\n",
      "   64    2.200e-04   2.00    1.417e-02   1.00    3.011e-04   2.00\n",
      "  128    5.499e-05   2.00    7.085e-03   1.00    7.529e-05   2.00\n",
      "\n",
      "f = 1:  ||sigma - sigma_h||_0 = 1.8463195827281382e-16\n"
     ]
    }
   ],
   "source": [
    "# Cell 4 [LECTURE] -- P1-P0: flux plot (n = 8) and measured convergence rates\n",
    "case = \"f=pi^2 sin(pi x)\"\n",
    "\n",
    "nn = 3\n",
    "mesh, w = solve_pair(nn, (\"CG\", 1), (\"DG\", 0), case)\n",
    "sig_h = w.subfunctions[0]\n",
    "xs = np.linspace(0, 1, 200)\n",
    "nodes = np.linspace(0, 1, nn+1)\n",
    "print('nn=',nn)\n",
    "plt.figure(figsize=(6, 3))\n",
    "plt.plot(xs, CASES[case][\"sig_np\"](xs), lw=2.2, color=\"0.6\",\n",
    "         label=\"exact $\\\\sigma$\")\n",
    "plt.plot(xs, eval_at(sig_h, xs), lw=1.6, color=\"tab:green\",\n",
    "         label=\"$\\\\sigma_h$  ($P_1$-$P_0$, $n=nn$)\")\n",
    "plt.plot(nodes, eval_at(sig_h, nodes), \"o\", color=\"k\", ms=4,\n",
    "         label=\"vertex values: $\\\\sigma_h$\")\n",
    "plt.legend(frameon=False); plt.xlabel(\"$x$\"); plt.show()\n",
    "\n",
    "print(f\"{'n':>5} {'rel sig err':>12} {'rate':>6} {'rel u err':>12} \"\n",
    "      f\"{'rate':>6} {'max mid err':>12} {'rate':>6}\")\n",
    "prev = None\n",
    "for n in (8, 16, 32, 64, 128):\n",
    "    mesh, w = solve_pair(n, (\"CG\", 1), (\"DG\", 0), case)\n",
    "    sig_h, u_h = w.subfunctions\n",
    "    x, = SpatialCoordinate(mesh)\n",
    "    es = rel_L2(sig_h - CASES[case][\"sig\"](x), CASES[case][\"sig\"](x))\n",
    "    eu = rel_L2(u_h - CASES[case][\"u\"](x), CASES[case][\"u\"](x))\n",
    "    mids = (np.arange(n) + 0.5) / n\n",
    "    em = np.max(np.abs(eval_at(u_h, mids) - CASES[case][\"u_np\"](mids)))\n",
    "    if prev:\n",
    "        r = tuple(np.log2(p / c) for p, c in zip(prev, (es, eu, em)))\n",
    "        print(f\"{n:>5} {es:12.3e} {r[0]:6.2f} {eu:12.3e} {r[1]:6.2f} \"\n",
    "              f\"{em:12.3e} {r[2]:6.2f}\")\n",
    "    else:\n",
    "        print(f\"{n:>5} {es:12.3e} {'':>6} {eu:12.3e} {'':>6} {em:12.3e}\")\n",
    "    prev = (es, eu, em)\n",
    "\n",
    "# For f = 1 the discrete flux is exact:\n",
    "mesh, w = solve_pair(8, (\"CG\", 1), (\"DG\", 0), \"f=1\")\n",
    "x, = SpatialCoordinate(mesh)\n",
    "print(\"\\nf = 1:  ||sigma - sigma_h||_0 =\",\n",
    "      float(sqrt(assemble((w.subfunctions[0] - (0.5 - x))**2 * dx))))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-05",
   "metadata": {},
   "source": [
    "## Modification 1: $P_1$–$P_1$\n",
    "\n",
    "Both fields in continuous $P_1$. The dimension count is satisfied with\n",
    "equality; the kernel of $B^T$ is spanned by the alternating (checkerboard)\n",
    "mode, and the system is singular for every $n$ (Lecture 1; Exercise 1).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "cell-06",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ConvergenceError raised; message:\n",
      "\n",
      "Nonlinear solve failed to converge after 0 nonlinear iterations.\n",
      "Reason:\n",
      "   DIVERGED_LINEAR_SOLVE\n"
     ]
    }
   ],
   "source": [
    "# Cell 6 [LECTURE] -- P1-P1: the direct factorization fails\n",
    "n = 64\n",
    "mesh, W, a, L = build(n, (\"CG\", 1), (\"CG\", 1), \"f=1\")\n",
    "w = Function(W)\n",
    "try:\n",
    "    solve(a == L, w, solver_parameters={\n",
    "        \"mat_type\": \"aij\", \"ksp_type\": \"preonly\",\n",
    "        \"pc_type\": \"lu\", \"pc_factor_mat_solver_type\": \"petsc\"})\n",
    "    print(\"solver returned without an error (unexpected)\")\n",
    "except Exception as exc:\n",
    "    print(type(exc).__name__, \"raised; message:\\n\")\n",
    "    print(exc)\n",
    "# The exact wording of the error depends on the PETSc build; a zero-pivot\n",
    "# message in the LU factorization is expected.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-07",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 7 [LECTURE] -- P1-P1: dense SVD and the null vector\n",
    "A_mat = assemble(a, mat_type=\"aij\")\n",
    "ai, aj, av = A_mat.petscmat.getValuesCSR()\n",
    "M = sp.csr_matrix((av, aj, ai)).toarray()\n",
    "sv = sla.svd(M, compute_uv=False)\n",
    "print(\"three smallest singular values:\", sv[-3:])\n",
    "\n",
    "U_, S_, Vt = sla.svd(M)\n",
    "wnull = Function(W)\n",
    "with wnull.dat.vec_wo as vv:\n",
    "    vv.array[:] = Vt[-1]\n",
    "signull, qnull = wnull.subfunctions\n",
    "print(\"norm of the V-part of the null vector:\",\n",
    "      np.linalg.norm(signull.dat.data))\n",
    "\n",
    "nodes = np.linspace(0, 1, n + 1)\n",
    "qvals = eval_at(qnull, nodes)\n",
    "qvals = qvals / qvals[0]\n",
    "plt.figure(figsize=(6, 2.6))\n",
    "plt.axhline(0, color=\"0.85\")\n",
    "plt.plot(nodes, qvals, \"o-\", ms=4, color=\"tab:red\")\n",
    "plt.title(\"Q-part of the null vector: the alternating mode\")\n",
    "plt.xlabel(\"$x$\"); plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-08",
   "metadata": {},
   "source": [
    "## Modification 2: $P_2$–$P_0$\n",
    "\n",
    "$V_h$ = continuous $P_2$, $Q_h$ = piecewise constants. The system is\n",
    "invertible for every $n$, and the method does not converge (Lecture 1;\n",
    "Exercise 2).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cell-09",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P2-P0 solve at n = 8 completed (invertible system).\n"
     ]
    }
   ],
   "source": [
    "# Cell 9 [LECTURE] -- P2-P0: solve at n = 8\n",
    "mesh, w = solve_pair(8, (\"CG\", 2), (\"DG\", 0), \"f=1\")\n",
    "sig_h, u_h = w.subfunctions\n",
    "print(\"P2-P0 solve at n = 8 completed (invertible system).\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cell-10",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAE1CAYAAAAvcQAkAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlcelbwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAnCNJREFUeJzsnXdYk2fbh88EEvZSQAQHuFBk495t3VoH1ro67Hprd+2wu2r3suvt+Nq3wy61rlpHtVZbV90MRREUEFEQEZUNCeP5/kjylMiQkZBg7/M4OFqecd9XEPL8ck2FJEkSAoFAIBAIBK0EpaUNEAgEAoFAIGgMQrwIBAKBQCBoVQjxIhAIBAKBoFUhxItAIBAIBIJWhRAvAoFAIBAIWhVCvAgEAoFAIGhVCPEiEAgEAoGgVSHEi0AgEAgEglaFraUNMAdVVVVkZWXh4uKCQqGwtDkCgUAgEAgagCRJFBYW4uvri1JZt3/luhQvWVlZdOzY0dJmCAQCgUAgaAJnz56lQ4cOdZ6/LsWLi4sLoHvxrq6uFrZGIBAIBAJBQygoKKBjx47yc7wurkvxYggVubq6CvEiEAgEAkEr41opHyJhVyAQCAQCQatCiBeBQCAQCAStCiFeBAKBQCAQtCqEeBEIBAKBQNCqEOJFIBAIBAJBq0KIlwaydu1aQkJCsLe3JywsjLVr11raJIFAIBAI/pUI8dIA1q5dy7Rp0zh+7BgajYaEhASmTZsmBIxAIBAIBBZAiJcGsHjxYhQKBZL+e0mSUCgUvPDCC1RVVVnUNoFAIBA0josXL3LHHXfI79/33HMPhw4datJaycnJTJ8+3ZTmtXrmz59PYmKiWfcQ4qUBnDx5EkmSjI5JkkRKSgqrVq3i5MmTQsQIBAJBK+HVV18lICBAnp2zfft2Lly40KS1rly5wubNm01pXqsnJCSEZ5991qx7CPHSAHr06FGj259CocDHx4f8/Hx27NjBypUrhYgRCAQtwiuvvMLq1asttv/u3buZOHGi/JWcnFzrdR999JEcXt++fTsPPfQQc+bM4b333qOgoKAlTZYpKCjg22+/5d5777XI/i3FM8880yxRtX37du65555az/3222/MmDGDsrKyWs/PnDmTv/76i5SUlCbvfy2EeGkACxculENFAAp0npdJkybJ1xQUFMgiJjk5WYgYgUBgNg4ePGjWB8O16Nq1K/PmzeP2229n06ZNXLlypdbrPvjgAzp37szLL7/MO++8Q1BQEOPGjWPlypUMGzYMjUbTwpbDL7/8Qq9eva774b1///03Z86cafL9zs7OLF26lNLSUqPjpaWlPPjgg0RERGBvb1/rvY6OjowaNYply5Y1ef9rcV3ONjI10dHRrFmzhpcff5xT587RxcWFN5YuZcyYMRw5coTExEQqKysBnYjZuXMnsbGxRERE0KNHj3rHegsEgtbN5cuX+fTTT0lMTKR9+/Y8+OCDdOvWjdLSUu69915uu+02xo0bB8CuXbv45JNP+OqrryguLua+++4DwMnJiYiICB599FEcHR3ltRMTE1m6dClnz55l6NChzJs3j6+//ppDhw6RmprKnj17aNeuHV9//XWttu3bt49PP/2UlJQUo0/JTzzxBHfccUeTX7Ovry++vr7k5eXVeU1CQgLl5eVERkbSqVMnXnnlFfnc6NGjadeuHYcOHWLIkCFNtsPAI488Qv/+/bntttsAiI+P55VXXpHFU3V27dpF3759a6xRXFzMq6++SmJiIkFBQTz55JNG/xarVq3ijz/+QKvVMnbsWGbOnFmnPWvWrGHz5s1UVlYyatQoZs+eDUBOTg533303S5Ys4euvvyYjI4OBAwfy0EMPYWurexyfO3eODz/8kOzsbAYNGoSHhwdpaWm88MILVFZW8r///Y89e/bg6OjIrFmzuOGGG2rs/+GHH5KYmMhnn33Gxo0b6dKlCx9//DGSJLF06VJ27tyJWq1m0qRJTJw4sdbXEBISAsCxY8eMfl5vvPEGjo6OPPnkk3W+foD+/fuzdevWeq9pDuKp2kCio6PZ+/XXxPcI5NdeQUydOhVHR0cGDhzIrFmzCA0NxcbGRr6+sLCQXbt28fPPP5OUlCSLG4FA0HgkSaKysrJFv67Oc6uNnJwcgoODSU9PZ8qUKTg5OdG3b19SU1NxcHAgOjqa2267jYyMDHJycpg5cyZjx47F1dUVNzc35s2bx7x585gyZQo7d+5kwoQJ8tpbt26lX79+spc3MTGR1157jcGDB9O5c2f69+/PvHnzmDNnTq22/fHHH4wZM4bhw4fz7bffMnbsWI4cOcJ//vMfhg8fLl936NAhoxDQ1V8PPPBAk/7NNmzYwMSJE1EoFHh5eRmdy83NBaBNmzZNWvtqhgwZwvz58yksLCQ5OZlx48Yxffr0GsIFIDU1lQ4dOtQ4/uCDD6LRaJg4cSJbt25l9OjR8u/AQw89xMKFC+nbty833HADixYtYtGiRbXasnDhQh5++GEiIyMZNGgQzz77LE888QQAJSUlbNq0iZkzZ9K1a1fGjh3Le++9x9tvvw3o8mf69OnDxYsXmThxIrGxsdx///1yMvGbb77Jhx9+yMiRIxk0aBBvvfUWu3fvrmHDDTfcQPv27WXBa0govvfee3njjTcYNmwYvXv35s477+T999+v9XU4OjrSrVs3jh49avSze++99/jss89QqVRs376dTz75pNb7O3ToQFpaWq3nTIFCashfaCujoKAANzc38vPzTTpVujThGOn6X4LAuFiUDg5G50tKSjh69CiJiYlUVFQYnXN2dpY9MdVFjkAguDaVlZV1ehfMxT333HPNv9X58+dz+vRp1q1bJx976KGHkCSJzz77DIDHHnuMAwcO4OzsjJ+fH999912ta2m1Wnx9fdmxYwfBwcGEhoZy66238uKLL8rX5Ofn4+bmxsSJExkyZEi9SZFdunThqaee4sEHH5TXd3FxYePGjYwaNUq+Licnh4MHD9a5jrOzMyNGjKj1XF5eHh4eHuzbt48BAwYYnRs4cCDPP/88N998c43XOWbMGNzd3fnll19qrPnAAw9w9uzZOu0xiIirGThwIMHBwWzdupUXX3xR9mpdTb9+/Zg9ezaPP/64fMzf35/Jkyfz0UcfAbqfc8eOHVm7di0dOnQgJCSE9PR0/Pz8AJ1XqW/fvhQUFBAbG8vIkSMpKiqiqKiItm3b8uuvvzJ27FgA9uzZw4gRI8jKyqKkpISAgAD27t3LwIEDAfj4449ZvXo1u3bt4p133mH16tVG/x4jRozA3d2ddevWMWfOHPz9/Xn99dcBnagvLCys9Tk3ZMgQbrvtNubNmwdAWlqaLEaCg4MBWLFiBffddx95eXm1/q7feuuttG/fXv65TJw4EXd3d3788UdAJ6Y8PT1r/Vn/8ssvzJs3r9GJ0A19fouwUSNQtfeR/788Oxu7gACj846OjgwYMICwsDCOHj3K8ePHZRFTVFTE7t27iYuLIzw8nMDAQCFiBIJWTlxcHNnZ2Uau9zNnzuDp6Sl//+6779KtWzcqKytZv3690f3r169n06ZNXLhwgYqKCrRaLWlpaQQGBpKYmMjIkSONrndzc2uQXYcOHSIzM9MoNGRra0tVVVWNPAVvb+86QwdN5eLFiyQkJNSwX6PRMH36dCorK/nhhx9qvXfGjBkUFRXVuXZt3hTQhcJuvfVWFi9eXKdwAfD09Kw1R2fw4MHy/7u5udG7d2+OHz9Obm4uKpWK+++/Xz4vSRIajYbTp08brZGamopWq2XYsGHyMYNISUpKolOnTgCyeADdz99gT2JiYq0i8MSJEwAsWLCAO+64g6NHjzJ06FAmTpxIUFBQna+1OidOnMDd3d1o76FDh1JUVMSZM2fo0qVLjXvCwsLYtm0bAJs2bWLPnj1GydmJiYncfvvtte535coV2rZt2yDbmoIQL43Apk0bFCoVUnk5FbWIFwMODg7079+f0NBQEhISOH78OOXl5YBOxOzZs4e4uDgiIiKEiBEIWjHOzs706dOnRv6Dh4eH/P+rV6+mrKyMyspKfv/9d6ZOnQrAN998w0svvcSCBQsYPXo0dnZ2xMfHU1JSgkqlwtHRkfz8/CbZlZCQgI+PD87OzvKxPXv24ODgUMNrcejQIRYvXlznWh07duTzzz9v1P6bNm3ihhtuwKGad7qkpIQpU6YgSRJbtmwxyiepTl1envrIyMjg6aefpnfv3nVWPhkICwsjKSmpxvHCwsIa37u4uODs7Iy9vb3swTDwwAMP0L59eyMh5OLiAui8B4bXV1JSQmVlpXwOqFG9agiAuLi41BBu1b8PCwvjyJEjpKens2vXLiZMmMCiRYu48847633NhrUNthieOYaKr+q2VScsLIz3338fjUbDY489xmuvvUa7du3k88ePH+fChQsMHz6cbt268b///U/O8Tx+/Djh4eHXtKupCPHSCBRKJbY+PpSfPUv5+exrXu/g4EC/fv0IDQ2VPTEGEVNcXCyLGIMnxpCwJRAIjFEqlXWWbZpzz2sxefJklixZwkcffSR/ykxNTZUrgZKSknjggQdYvXo1paWlzJ07l/DwcAICAti/fz9jx47lscceA+DAgQNGLvYJEybw3nvvMWzYMBwcHCgsLGTfvn2MHj1adqvXhYuLC9nZ2eTl5eHu7o5Go+HZZ5/lP//5Tw3PS+fOnWs8mKtTXQA1lA0bNhiFiwoKCpgwYQIeHh6sWrUKOzu7Rq9ZF9nZ2YwcOZJHHnmE6dOn07NnTw4ePEi/fv1qvX7s2LHMmTPHqIIU4Ouvv2bOnDnY29uzc+dOTp48yfDhw/Hw8MDGxoaysjJuueUWAMrLy/nuu+9qhDX8/f3p1q0b//3vf+XQzocffoifnx9BQUGcP3++3tdy0003cc8995CVlYWvry85OTmsWrVK9sb8+OOPTJo0CX9/f/z9/dmxYwc7duyoVbxc/TsSERGBs7MzX375pZzH9NFHHxEeHl4jJ8lAWFgYly9f5vHHH8fNzU0OQQJUVVWRnJxMSkoKb775Jg899BBJSUmyJ2j37t08/PDD9b7eZiFdh+Tn50uAlJ+fb/K102+7XUoM7CnlfPppo+8tLS2VDh48KH3zzTfSF198YfT1448/SgkJCVJ5ebnJbRYIBOahqqpKmj9/vuTu7i4NGzZMCg8Pl4KCgqR9+/ZJxcXFUnBwsPTyyy/L18+fP1/q06ePpNFopJ07d0rOzs5S//79pcGDB0vBwcFSx44dpeXLl0uSJEnZ2dnS8OHDpXbt2knDhw+XOnfuLG3evFmSJElatWqV5OLiIo0aNUq6++67a9hVVFQkBQcHSyEhIdKjjz4q9erVS5o6daqk0WhM8rrPnj0rTZgwQRozZowESIMHD5YmTJgg7du3T9JoNJKLi4t07tw5+fp77rlHUigU0pgxY6QJEybIX3/99Vez7Lh06ZIUHBwsvfjii/Kxp556Sho8eHCd91RVVUk9evQw2rtz587S+PHjpc6dO0tDhw6VHB0dpQ8//FA+v337dqlTp05SSEiINHz4cKl9+/bSSy+9JEmSJO3bt09ycnKSr/37778lX19fKSwsTIqMjJS8vb2l7du3S5IkSadPn5YAqbCwUL5++fLlUu/eveXv77zzTsnd3V0aPny41KlTJ2nYsGFSdHS0JEmS9NVXX0kdOnSQhgwZIkVFRUleXl7S/v37a32dX331leTm5iaNGTNGeuSRRyRJkqT169dLbdq0kfr16ycFBQVJ/v7+UlxcXL0/4zZt2khKpVI6cOCA0fGUlBQpKChI/n7ChAlSZmamJEmSlJqaKrm7u0vFxcX1rl0bDX1+i4TdRpK5YAEF6zfgPn067V995do31EJZWRnHjh2TSwmr4+joSHh4OD179hSeGIGglZCbm8uJEyfkXAkbGxsyMzM5cuQIY8eOlb045eXlbN26lT59+tCuXTsuXbrEkSNHcHR0JCoqiv3799O1a1d8fX3ltU+dOkVWVhbh4eFGOS8ZGRkkJydjY2PDjTfeWMOmwsJCNm3axMWLF4mMjDTK6WguxcXF/PXXXzWO9+3blyNHjvDcc88RExMjH4+LiyMzM7PG9REREXISbFNITk7m7NmzRrk1BQUF7Nq1ixtvvLHO0NRPP/3EihUr2LBhAwB//vknwcHBsjehe/fuRv8GoEs0TkpKIi8vj9DQUNzd3QFd0rLBi2agtLSUY8eOUVlZSUhICE5OTvLx7du3M27cODl0k5WVRUpKilGezMmTJ8nJySEkJISnn34apVLJ//3f/wG6n31cXBy2trZERETU68VKS0sjJSUFBwcHhg4dCujCUAkJCajVaoKDg6/pBTNUMxnuN7BhwwY2b94sJ6ZHREQQFxcHwMMPP4yvry/PP/98vWvXRkOf30K8NJKc9z/g0pdf4jR0KJ3+92Wz1jIMeTx27BhardbonIODA+Hh4fTq1UuIGIFA0Gp45JFHaNu2bZ2lxNaApM+7GTNmjNX14frtt98YP348AOnp6URERPDll19a3fykt956Cx8fH+bOnUt6ejqPP/64XHW3bds2Bg8ebJTz1FCEeDGTeLmyfDnZi1/Brns3uuhVe3PRaDSyJ6Y2ERMWFkZQUJAQMQKBwOrZu3cvXbt2NUrsFDScl156iWXLluHn50dsbCyzZ8/miy++qJHka2mys7NxcXHBycmJkpISCgoK8PHxufaN10CIFzOJl8K//uLcAw+idHYm8HDTppDWRUNETK9evVCpVCbdVyAQCATWQ2ZmJqdOnaJr167X/RiDqxHixUzipSwpidNTdKWOPQ4fwqYJmfjXQqvVyiLm6tkf9vb2sidGiBiBQCAQXE809PltXcG+VoCqmlus4hplb01FrVYTGRnJrFmz6Nu3r1FCVVlZGQcOHGD58uXEx8fXSPgVCAQCgeB6R4iXRqJ0c0OhT0Iqz752r5fmoFariYiIqFPEHDx4UBYxV4eZBAKBQCC4XjF7BuiVK1f44YcfOHPmDN27d+eOO+6os3ztapKSkvjkk0+IiIho8QZVdaFQKFD5+KA9fZpyM3lersYgYoKDgzl+/DhHjx6VJ8QaRMyRI0cIDQ2ld+/eqNXqFrFLIBAIBAJLYFbPy4ULF4iMjGTVqlW4u7vz5ZdfMmDAgHrnVhgoKytjxowZrFy5ks2bN5vTzEZjmHFUYWbPS419VSrCw8OZNWsW/fv3N+qUqdFoOHToEMuXLyc2NlZ4YgQCgUBw3WJW8fLaa69hb2/Ptm3beOmll/jrr7/Izs7m448/vua9TzzxBIMHD66zxbMlsfVpD0B5Vst4Xq5GpVIRFhYmi5jqtfQajYbDhw8LESMQXAd8/vnnXL582dJmWDW5ubl88cUXljZD0MKYVbz8+uuvTJ8+Xc7VcHNz4+abb+bXX3+t9761a9eyY8cOlixZYk7zmoyqvV68tFDYqE47qomYAQMG1Cpili1bRkxMTI2qJYFAYN38/vvvLF26lDZt2ljaFKvG09OTL774otaOv4LrF7PlvGg0Gs6ePUvAVZOXu3Tpwi+//FLnfRkZGTz44INs3ry5wd35NBqN0cPZMCnTXKj07azLa2l3bQlsbW0JDQ0lKCiIEydOEB8fT2lpKaAru46JiSEhIYHg4GBCQkJMOhRNIBCYhzfffLPGYLsVK1aQnp4O6EaJ9O/fn/79+zd67R9++EFu1+/s7MzAgQOJiopqts3m5O+//yY+Ph4HBwduvPFG/P395XMPPfQQr7/+OjfccIPlDBS0KGbzvBgenleP2nZ1daWkpKTWeyoqKpg1axZPPPEEERERDd7rzTffxM3NTf4yd1MflX7mRXl2NlJlpVn3agy2traEhIQwa9YsBg4caJQYrdVqiY2NZdmyZRw+fFhO+BUIBNbH6dOn2bdvH1OnTjU6/vTTT7Nv3z7y8vKIjY1l6NChfPTRR41e/7HHHuPw4cPyXJ7+/fvz9ddfm8p8k3PPPfcwdepUEhMT2bp1K7169WLNmjXy+ejoaHbu3ElGRoYFrRS0JGbzvDg5OaFQKMjLyzM6fuXKlTobz/z+++8cPnyYsLAw+RPH8ePHsbGx4eGHH+aVV16p1YX63HPP8cQTT8jfFxQUmFXAqDroB4lVVFBx8aJR7xdrwCBievXqRVJSEvHx8bJgLC8vJzY21sgTUz3xVyAQNJ6dO3eye/duLly4gKHvp0ql4oMPPmjSetu3byckJATnak0ws7OzOXfuHCtXrmTgwIEAeHh48OWXX/LYY481eO2UlBSuXLnCwoULCQkJAXQdvL/88stmVXX+9NNPcrXj9u3b8fLy4pZbbmn2WJOcnBy++eYb9uzZIw+XfOaZZ3jnnXeYNm0aoPs59OzZk7/++os777yzWfsJWgdmEy8qlYru3btz4sQJo+MnTpwgKCio1nuCgoJq5Lns2LEDlUpFz5496+woa2dn16KhEFW7dqBUQlUV5ZmZVideDNja2hIcHEzPnj1JTk4mPj6e4uJiQCdi4uLiOHbsGL179yY0NFSIGIGgCcyfP5+1a9cyf/58bG1teeWVV+jevTuTJk2isrKSd955h+eee65Rax4/fpyuXbsaHTt0SDeOpHfv3vIxLy8vqqqqGrX2oUOHsLW1JTAwsFnrXM2TTz5JaGgoSqWSHj168Nprr7Fv374anqFvv/2WCxcu1LnO1KlTjWyztbVFpVIZtYCws7Or0XKje/fuHD9+vFmvQdB6MGuflxkzZvDNN9/wwgsv0KZNG86cOcPGjRt555135Gt+++03duzYwTvvvENAQECNGO+WLVuwt7evcdySKFQqbNu1o+L8ecqzssDKY8W2trb07t2bnj17yp6Y6iImPj6e48ePk5OTw48//khKSgo9evRg4cKFREdHW9h6gQCkigoqcnNbbD9bT08UDfAYbN++nW+//ZaEhATZ25ubm0tcXByvvvoqycnJ7Nmzp9H75+Xl1Qi5Hzp0iE6dOsme66qqKjZs2MCQIUOQJImffvqJo0ePMn78eEaMGFHn2ocOHaJHjx6yGKioqGDTpk3ceOONAGzatIndu3czbtw4hg8f3iB7MzMzuXDhAmPGjOHJJ58EdB9Gv/322xrXFhYW1vDIV+fqruFt2rThs88+Y968eYwbN47Lly9z4MCBGmu7urrWu67g+sKs4uWZZ55h27ZtREVFMWjQIP766y+GDx/OvffeK19z8OBBvvzySyNB0xpQ+frqxIuVJO02BBsbG1nEJCcnExcXJ4uYAwcOyJNLJUkiISGBadOmsWbNGiFgBBanIjeXlBEtl4zZbcdfDfKofvfdd8yaNcsoTO3p6Sk/gBMTE40SSxuKq6trDe+EwWPy1ltvUVpayu+//05RURFvvPEGzzzzDNnZ2fTq1Ytbb72VjRs30q9fPyMvR1hYGOPGjePQoUNUVVXx1ltvUVxczG+//YaNjQ2LFi3i999/56mnnuKOO+5g7ty5rFmzhsjIyGvaGxMTg4eHB48++qh8LCcnBz99cUN1ql/TUJKTk9FoNOTn51NUVERRURGnT58mNDRUvqagoIAuXbo0em1B68Ss4sXJyYndu3ezbds2MjIyuO+++xg+fLjRaO/x48fX+gtu4IEHHsDGxsacZjYJla8vpTExlGdmWdqURmNjY0NQUBCBgYGcPHmSuLg4Nm7cKAsXAEmSUCgULFq0SIgXgaAO9u7dyzPPPGN0LCYmRq7cOX78OB07duSjjz7iypUrLFiwoEEdxnv16sX+/fuNjh0+fJjw8HDy8vJwc3Pj2WefZeLEidja2vLYY4/J76MFBQUkJibSr18/Iy9HaWkplZWVxMXFMXz4cPLy8nB3d2fRokVMmDABpVLJ0qVL+fjjjxk1ahQdOnTgxx9/bJB4iY2NZfDgwUah/cOHD9daCdXYsNGWLVv4+OOPycrKom3btgAsW7aM2267jcuXL8t7nj59mvHjx1/TVsH1gdnHA9jY2DBmzJg6z/fr16/eRnQTJkwwh1nNRuWnrzjKan3ixYCNjQ29evWiR48e3HXXXVw9YFySJE6cOMH+/fsJDQ1t8FgHgcDU2Hp60m1Hy/XxsPX0bNB1VVVVRqGKpKQkNm/ezO7duwGd5yUzM5Nx48bx559/0q5dOx544IFrrnvDDTcwf/58SktLcXBw4PTp0+Tm5vLSSy8xbNiwGtcbhEtRURE7duyQE3iv9nIkJCRQXFzMq6++WqsoycjIoFevXoBOQK1bt04+9/777zNq1Cg5ybc6sbGx9O3b1+jY4cOHaw33NzZsVFpailKpNMp5sbe3R6vVyu9ZhYWFHDt2TJRK/4swu3i5XpHLpVtR2KgubGxs6NmzJwlHj1JdvigUCnx8fDh69CjHjx8nKCiIsLAwIWIELY7C1tYqE+Ojo6N59913adu2LVqtltdee42XX35ZbvVw/PhxvvrqK/r27YuXl1eDWxQEBgYSHh7Ohg0buPXWWzl06BB2dnb19nQpKytjzpw5vP322/jU8bM6dOgQrq6uhIWF1XrexsaGSn37h4qKCrlSqKqqiueff57OnTvXKl5iYmKYN2+e/H1mZibnz5+nT58+Na5tbNho9OjRdO/enUGDBjFlyhQKCwv5/vvvefjhh2VBs27dOgYMGFAjyVlw/SKmSjcRuVHd+fM1PBatkYULFyIBhoCeIYQ0ceJEACorK0lISGD58uXs3bu3zl49AsG/iXfeeYeXXnqJ33//nYSEBJYtW8aCBQsA3d/MlStXZI9EbGws4eHhDV77mWee4dNPPwV0eTTvvPNOnVWVJSUlTJs2jQcffLBWz4yB9u3b8/bbb9cZig8MDJTDVfv376dHjx6ATpz4+fkxadKkGveUl5dzxx13MGDAAPlYWVkZr776qkm6Azs5OXH48GFeeukl1Go1fn5+rF+/3qgy9fPPP290RZegdaOQrocn71UUFBTg5uZGfn5+nT1lmovm9GnSxuniq91378LWy8ss+7Qkn/fvz0dHjnBaq6WHvz+L3n2X4OBg4uLianQtNoScwsLCcHJyspDFAoH1cvLkSZ544gk2btwIwKBBg9i+fXuDO4eDrgHn/ffff00RcPfdd3Po0CG5D8rUqVPrDdfXRUxMDJMmTWLo0KHs2bOHXbt20aVLF/bv309paalVhmVyc3P5+uuva+QeCVonDX1+C/HSRKo0GpLDwgHw/3kFDnW4YVsLUmUlyZFRSPoxC56PPoLXgw8COpdxSkoKcXFx5OfnG91nCDmFh4cLESMQVCMzM5OMjAy5odzq1au55ZZbzLLX5s2bOXPmjPz9gAEDGuXlqc6pU6c4fPgwgwcPplOnTiayUCBoGEK8mFm8AJwcOpTKi7n4vb8E11ae5a49e5bUUaPl710nTsTvvXeNrqmqqiI1NZXY2NgaIkapVMoipnpXUIFAIBAIGkpDn98iYbcZqHx9qbyY26orjgxoUlONvtempdW4RqlU0r17d7p27VpDxFRVVZGYmEhSUpIQMQKBQCAwK0K8NAO1nx9lR45eF+JFm3ba6HtNerrc6+VqqouYtLQ0YmNj5dLH6iImMDCQiIgIIWIEAoFAYFKEeGkGhnJp7XVQLq09q5vGat+7N2XHjyOVlFBx4UK95alKpZJu3brRpUsXTp8+TWxsLFeuXAF0IubEiRMkJyfTo0cPIiIiarQ7FwgEAoGgKQjx0gzkcunrQLyUn9O9BscB/dGcPIlUXo42La1BvTWUSiVdu3alS5cusiemuohJSkoiOTlZ7l1hzjwkgUAgEFz/iD4vzeAf8ZLV6nu9lJ87B4C6U2fU+lksmqtCSddCoVDQtWtXbrnlFkaOHGlU3ilJEklJSfz888/s3LmzRum1QCAQCAQNRXhemoGqg24Ym1RaSuWlSw1uKW5tSFVVsvdI1cEPdUAAmlOnak3abQgKhYIuXboQEBBAeno6MTExXL58WbeXJJGcnMzJkyfp3r07kZGRwhMjEAgEgkYhxEszUHXwA4UCJAltxtlWK14qLl5E0moBUHfogLpLAACa000TLwYUCgUBAQH4+/uTnp5ObGwsly5dAnQi5uTJk5w6dYru3bsTERGBm5tb816IQCBoEb777jtyc3N58sknLW2K4F+KEC/NQKlWY+vjQ8X585SfOwuREZY2qUkYQkYoFKjat8dOP1ZeezrdJOtXFzFnzpwhJiamVhHTrVs3IiMjhYgRCKycEydOcM7wviEQWACR89JM1B06ALomb60Vg3ixbe+DQq1GHaATLxXZ2VQWFZtsH4VCgb+/P9HR0YwZMwbPap4qSZI4deoUK1eu5M8//6x36qxAcD3z999/N2jytMDyPPPMM2zevNnSZvwrEZ6XZqLq1BEOHaI8o/WKF60hWddPJ8TUAQH/nDudhkMtU2Sbg0KhoHPnznTq1ImMjAxiY2O5ePEioBMxKSkppKam0rVrVyIjI3F3dzfp/gKBNXP+/Hn++usvS5shaAB///03AdXeLwUth/C8NBN1R13SrrYVu1DLz+psV+m9SDbOTtj6tgdAcyrFbPsaRMyUKVMYO3YsXtWGWxpEzMqVK9m+fbtcei0QNJS1a9cSFhaGg4MDYWFhrF271qTra7VaoqOjSUpKMjr+ySef8N5778nfX758mVdffZVZs2bxxBNPkJLyz9/U7t27eeihh/j777955JFHmDt3Lm+99Rbnzp1j4sSJTJw4kYSEhHrXqays5D//+Q/Lli2T101ISGDKlCmcP3++ht1JSUlER0ej1ee5GXjkkUfkIZK//PILEydO5Oabb+aee+5h5cqV9f4slixZwo8//mh07D//+Y88oRogOTmZBQsWMHPmTBYtWmQ0YuTChQu8/PLL3HrrrSxYsICMjAyj+yZOnCiHmmujrrUvXrzI1KlTiY2Nla9dvXo1d999N+Xl5QC89dZbTJw4kUmTJvHwww+zb9++Guv/9ddfPProo9x+++3y79GHH35IYmIin332GRMnTuTRRx+t92ckMC1CvDQTlV68lFf7Y2ttGEJe6k4d5WN23boBoEkxn3gxoFAo6NSpE1OmTGHcuHF4e3sbnU9NTWXVqlVs27ZNrloSCOpj7dq1TJs2jYSEBMrKykhISGDatGkmFTBqtRpJkvj666/lYxUVFbz++ut00eeN5eTkEBwcTHp6OlOmTMHJyYm+ffuSqh/HkZmZydKlS3nmmWfo378/M2fOZNSoUXh4eDBv3jzmzZuHr69vvevY2Nhw2223MW/ePBISEigqKmL69On069eP9u3b17C7R48eHDx4kA0bNsjHTp06xeeff06YfsBsWFgY8+bN4/7776dv374sWLCAjz76qM6fRVxcXA0R9+eff5KdnQ3oRNqAAQOws7MjOjqatLQ0+vfvj0Y/CHby5MkkJCQwffp0fH19mTx5srzOlStX2LRpE6WlpbXuXd/aXl5eDBkyhOnTp5Ofn09SUhL33HMPs2bNQqVSATBq1CjmzZvHvffei6+vL+PGjWPv3r3y+h988AG33HILvr6+jB07lp9++ol169Zxww030L59e4YOHcq8efOYPn16nT8fgRmQrkPy8/MlQMrPzzf7XiVHjkiJgT2lxMCeUmVJidn3MwfJQ4ZIiYE9pbyNG+Vj2W+/IyUG9pTO3Htfi9tTVVUlZWRkSOvWrZO++OKLGl9//PGHdOnSpRa3S9B6CA0NlRQKhQTIXwqFQgoLCzPpPr/++qvUrl07qby8XJIkSVq/fr3k5eUlabVaSZIk6fHHH5cmT55sdM+DDz4oPfDAA5IkSdLy5cslOzs76eLFi/L5VatWSYGBgUb3XGsdSZKkN954Q+rRo4d0yy23SGPHjpWqqqrqtPu5556TJkyYIH///PPPS2PGjKnz+m3btkldu3aVv3/mmWekOXPmyN/PmTNHeuGFF4zu6dq1q/TLL79IkiRJAwYMkN58802j85GRkdJ3330nSZIkubq6Snv37pXP5eXlyf9/5coVacOGDVJpaWmttl1rbUmSpMmTJ0uTJ0+WevfuLS1cuLDO1ylJkvTaa6/Jr62oqEiyt7eXNmzYYHSNwb7BgwdLn3/+eb3rCRpHQ5/fIuelmRg8L6DrtGvwWLQWqoqLqbyYC+ga1BloSc/L1SgUCjp27EiHDh3IzMwkJiaGCxcuyOfT0tJIS0sjICCAqKgoo2Z4AgHAyZMnazSOlPQ9hkzJeP00+c2bN3PzzTfz7bffctttt8mf6uPi4sjOzmbixInyPWfOnDFKVu/cubPR97XRkHWeffZZVq1axebNm0lPT691LpmBu+66i6CgIM6fP0+7du34/vvvWbJkiXy+qKiIb775Ru7RVFJSwpkzZ6iqqkKpbJzDXpIk4uPjAdizZ498/MKFC5w4cQKA999/nzvuuIOhQ4cyePBgI8+Lu7u70etu7NoA33zzDR07dqRnz568/PLLRmucO3eOr776iuTkZIqKisjMzMTe3h7QhaPKysq46aabjO4RFZGWR4iXZmLj7o7S2ZmqoiK0GWdbnXipnqtjFDbq3h2AivPnqSwqwsYCwxUVCgUdOnTAz8+PrKwsYmJiZDc0wOnTpzl9+jT+/v5ERUXRtm3bFrdRYJ306NGDhIQEIwGjUCgIDAw06T62trbcdtttfPvttwwcOJCNGzcSExMjn3d2dqZPnz7MnDnT6D4PDw/5/w1Cpz4ass727dtJTU3F3d2dVatW1Vux1L17d/r3788PP/xAWFgYJSUlRoJh9uzZFBUVcfvtt9OmTRuysrL4888/0Wq18oO9OjY2NlRWVhodKykpAXQ/d2dnZ8aPH09EhHE7CUN47Z577mHu3LkkJCSwadMmevXqxcGDB6+ZDNuQtQE+//xzvL29OXnyJIcPH6Zfv36ALm9pwIABjBo1iptvvhkXFxe2b9/Ojh07gH9ESkFBAQ4ODvXaImhZRM5LM1EoFP/kvZxrfRVH2jNnALBxc8Om2qcJu67//OFrLeB9qY5CocDPz4+bb76ZCRMm1Ijjp6ens2bNGrZu3Upubq6FrBRYEwsXLjSaiq5QKJAkiYULF5p8r7vuuotNmzbx4YcfEhYWRki16rzJkydz+PBhBg4cKCfg9urVi6KiojrXc3NzM0pmbcg6WVlZzJkzh88++4yVK1eyYMEC4uLirmn30qVLWbp0KbNnz8bOzk4+t3//fubPn89dd93F5MmTa038rU7nzp05dOiQLBa3bdtmdM+kSZOIj49n9OjRsv2urq5UVVVRVlbG0qVLUSqVhIeHs2DBAtRqtZxke62E3frWBti5cydvv/0269ev5/XXX+fWW2+VCwDOnTtHZmYmb7/9NrNnz2bs2LGcOnVKXjsgIIDevXvz2muvyetlZGRw6NChOv+tBC2DEC8m4J9eL62v4qhcn6yr6tTJ6LjS0VGuPrJE6Kg2qouYiRMn1ipi1q5dy++//y5EzL+c6Oho1qxZQ2hoKPb29oSGhrJ27VqmTp1q8r169+5NeHg4b775JnfffbfRuXvvvZfx48fTrVs3hg8fTkREBJMmTao37NC/f39sbW3p06ePXG1U3zqVlZXMmjWLKVOmMGfOHAYNGsRLL73ErbfeWu8MsRkzZpCRkcHKlStr2H3//fcze/ZsRo4cSVBQkFHVUG385z//ITk5mZCQEIYMGcJrr71Gu3bt5PNLlixBq9XSqVMnbrrpJnr06MGiRYvw8PBApVKxe/duOnbsyMiRI+nevTudOnVi1KhRwLUTdutb+8KFC8yaNYuPPvqIkJAQHn30UaKiopg7dy4A/v7+jBkzhtDQUEaNGkX37t2N9lEqlfz000/88ccfdOvWjaFDhzJ27Fic9Z7o6Oho3nzzTcaOHSuqjVoYhXR1YPg6oKCgQFbELTE358K773L5629wHj6cjl/8n9n3MyXnX15I3sqVuE6YgN+S94zOnZ33AEU7dtDmzjtp99yzFrKwfrKysoiNjSUrK6vGuU6dOhEVFWVUgi0QmIMTJ06QmprKiBEj5AdbdXJzczlx4gRubm707t0bGxsbQPf7m5KSwrBhw4yuLy4uJjY2loKCAgYMGCCHRGtbJzc3l/379zNy5Eg5pCNJEr///jtBQUF0uuqDSXX27dtHfn4+Y8eOrXEuJSWFtLQ0uSfT9u3bGT9+PEqlkqSkJMrKyggPD69hs4eHB0FBQezcuZOgoCAjEZORkUFqaiqdOnWia9euRvtlZ2eTmJiIt7c3wcHB8vG8vDz27Nlj9Ppqo7a1U1JSOHPmjFHOSkFBAbt27WL48OG4uLgAyKNLevfuDegqHIcOHSrfU1FRwYkTJygsLCQyMtLIjrS0NFJSUnBwcDC6R9A0Gvr8FuLFBFxZsYLsRYtRd+1K100bzb6fKTlz112U7NtP2wfm4f3YY0bncpYs4dL/vsJp0CA6ffN1HStYB+fPnyc2NpZM/YDJ6nTq1InIyMgaJdgCgUAgsC4a+vwWYSMTYJguXX72LJI+LtpaKD+j609TvdLIgCUrjhpL+/btmTBhApMmTcLPz8/oXEZGBuvWrePFF1+kd+/eZmtaJhAIBIKWQYgXE2Co0pG0WipycixsTcOp0mop1yfVqTvXdC3LFUc5OVTWEzu3Jnx8fGQR00GfswM6t/Drr7/OiRMnzNa0TCAQCAQtgxAvJkDVvj3Y6qrOtWdaT6fd8nOZoI8aqqv1qzGg7tIF9D0dWoP3pTo+Pj6MHz+eyZMn07FjRzZu3ChXnAByJcpLL71kYUsFAoFA0FiEeDEBCpXqn4qj9HTLGtMItBm6MmmFoyM2tTTJUtrbo+qorzg6earG+dZAu3btGDduHLm5ubU2LTt58iSbNm0y6h8jEAgEAutGiBcToe6syxkx9E1pDRgmYas7dqyzG6d9D11Tr7LkpFrPtxYCAwNrvEaFQoGPjw+ZmZmsX7+eTZs2XbOfhUAgEAgsjxAvJkLt7w+0Ns+LIVm37lJKu5468aJJMm1b9ZZGblqm/16BzvNSve14ZmYmGzZsYOPGjbWWXgsEAoHAOhDixUSoA/yBViZezurEi6pTzXwXA/Y9ewKgSU5udZVU1YmOjuaLUaPoYWeHWqGgh50dq3/4gVdeeYXOnY0rrbKysti4cSMbNmwQIkYgEAisECFeTITseTl7FqmiwrLGNBC5TLpjPZ6XQJ14qSopofxc6+sgbECqqGDYhRx+8Q8gvkcgv/gHMNrXFy8vL8aMGUN0dDT++n9DA+fPnzcSMddhSySBFeLp6XnN1v6tkUGDBvHzzz9b2gzBdYIQLybCIF4oL6e8FXxal8rL5aGM6nqGn6n8fFHqu1CWJbXevBft6dNIGg0A6m667pslBw7K5z09PRk9evQ1RUxmZqYQMYImk5GRgbOzM8XFxXVeU1RUVGPIYVNIS0vD2dkZjf733tKUlJRQXl5uaTME1wlmFy/x8fHccccdDB8+nHvvvZeUa5TcFhUV8cEHHzBlyhQmTJjAokWL5CFa1oyttzcK/dTR1hA60p49C3oPkSHkVRsKhQK7wB4AaFqxeDEILxtPT1zHjQOgRD9crToGETNt2rQaE22zs7PZtGkT69ev59y5c0LECBpNx44dyc7OxsnJyex7BQQEkJ2dbTRwsT5Gjx7N119bdydta+SWW25h2bJlRse+//57nJ2deeGFFyxk1fWPWcXLsWPHGDJkCM7Ozjz77LMUFxczYMAAztUTfhg2bBjnzp1j7ty5zJs3j61btzJo0CAKCwvNaWqzUSiV/1QctQbxcvo0AEpnZ2yvMfvHXh86KmvFSbsG8WIfGIhTv37ysbqa77Vt25ZRo0Zxyy230KVLF6NzFy5c4LfffhMixoopz8qi9PjxGl+W9ooqFIpaZx9Zw17CM1I35eXl5OXlyV+Gad6lpaVs2bKFkSNHyteeO3eOxYsX06VLF6vxel2PmFW8vPrqq4SHh/PZZ58xbtw4fvzxR9zd3VmyZEmd9+zatYslS5YwZcoUbr75ZtavX09ycjKbN282p6kmoTVVHGnS0gBdyKiuMmkD9r30Sbut2POiOaGz3a5nIPahoSjUaqiqoiQmpt772rRpw8iRI+sVMb/++itnz54VIsZKKM/KInXsONKn3VLjK3XsOJMKmDVr1tCnTx9efvllgoKC8PX15cUXXyQ9PZ2JEyfi5eXFgAEDSEhIAGqGjUpKSnjggQfo2LEj4eHhfPfdd0brjxs3jpdeeonJkyfj6+vLwIED+fvvv+Xzly9f5r777iMgIIBu3brx6KOPyg/Wq8NG06dP57nnnmP27Nl06NCB4OBg1q1bB8C8efPYt28fjz/+OM7OzgwaNKjGa128eDHTpk0zOhYTE4Onpyd5eXl89913ODs74+zsjL+/P3fffXedXvOPP/6Y6dOnGx27/fbbeeONN+TvN27cyJAhQ/Dx8WHgwIFs3bpVPvfTTz8RHh5Ou3btGDVqFEeOHJHPde3alS+//LLWfRvLxo0bCQkJoVOnTnh6euLh4YGHhwd33nknANu3byckJMRobto999zD4sWLxSw1M2NW8bJ9+3Zuvvlm+XsbGxsmTJjAtm3b6rzn6k8KDg4O2NjYoNVqzWanqVD7tybPSzoAdl3qzncxYEjaLc/KajVjAq6mLFnnNbLv2QulnR0OYWEAlBysGTqqDYOImT59eo1puDk5OWzevJlff/2VjIwMIWIsTMWVK0h1vF9IWi0VJgxDl5eXExMTQ1lZGdu2beOHH37gnXfeYcyYMTz22GMcO3aMyMhI7r33XgCqqqooLi6Wf0ceffRRYmJi2LJlC6tWreL77783+rReWlrKkiVLuO2224iJiSE6OprRo0fLVXDTp0/nzJkzbN26lfXr13Pw4EHuvvvuWvcqLS3l008/ZebMmcTGxnL//fdz++23U1BQwEcffUT//v15++23yc7OZvv27TVe68yZM1m/fr1RQ8dvvvmGm266CXd3d+bMmUN2djbZ2dls27aNK1eu8Mgjj9T6c9NqtZSWlhodKy0tld/nN2/ezO23387TTz/N0aNHWbBgAbfccguJiYlkZ2dzxx13sHjxYo4fP84LL7zAZ599Jq9TXFxskufFgQMHmDlzJm+//Tbnz58nLS0Nd3d33n//fTm8tmHDBqNn3Oeff45KpeK2225r9v6C+rE118LFxcVcunQJX19fo+O+vr6caUQjt3feeQc7Ozsjt9zVaDQaoz/4Ags9YA2eF02rEC+6sFF9yboG7Lp3040JqKpCk5yMY9++5jbPpJTn5FB56RIA9vq+NY79+1Ny6BDFB/Y3ai0PDw9uuukmIiMjiYuLIzU1VX445OTksGXLFry8vIiKiqJjPc3/BNcPXl5evP322ygUCnx9fYmIiGDUqFGMGjUKgAceeICoqCiqrmo1UF5eznfffcfOnTvp3bs3oHv4BQYGGl132223yV6Kp59+mlWrVrFs2TKmTJnCn3/+SXp6ulzu/8knn9C3b18u6X/fr+auu+5i0qRJADzyyCM8++yzJCUl0a9fP5RKJXZ2dnWGmgIDA+nbty8//vgjTz31FBqNhhUrVvDTTz8BYGtrK9/brVs33n//ffl1NZb33nuPp59+msmTJwMwdepUNm3axPfff8+jjz6KSqWie/fueHp6MmLECEaMGCHfm5aWhkqlqnPt3r171/sM2rhxIyNGjGDRokU89NBDjB8/HtBNp+/fvz8XLlzA3d0dSZLYuHEjW7Zskfd9/fXXOXDgQJNes6BxmM3zYoidXp0s5uDg0OC46urVq3n99df58ssv8fHxqfO6N998Ezc3N/mrYy1zeloCO714qcg6T1VZmUVsaChaQ9jI/9riRWlvL4uc1pj3Unb8OACKaq/DadBAADSJJ6i4fLnRa3p4eHDjjTcyffp0unXrZiRSLl68yJYtW1i3bh1nzpwRnpjrHG9vb6N/fwcHB6P3K8N73tXegKysLCoqKozEytW/S0ANMdOjRw8yMjLIyMjAzs7OqE9RT31fproezu3btzf63sHBgZKSkoa8TEAnfr799lsAfv31VxwdHRk9ejQAp0+fZsaMGfj7++Pq6kpISAilpaV1Cqn6SElJ4ZVXXpHDUM7Oznz//fekp6fj6+vLihUreOihh+jfvz8PPfQQx/V/4wCOjo71ipdDhw7JHqLavoYOHUpxcTHbt28nOjra6N68vDy89DmCsbGx2NraEhISAsDcuXNZvHhxjan2AvNgNvHi4uKCSqXi8lUPhkuXLtG2bdtr3r9+/XrmzJnDJ598wuzZs+u99rnnniM/P1/+Onv2bLNsbyqqam8i1jygseLKFSrz8wFQNyBsBLpEV4CypBNms8tclB1PBHQN9xT6AZoOISEo9RUfJfsb532pjru7uyxiunfvXkPE/P777/zyyy+kp6cLESMwwsfHB6VSaSQ0asudulqInDlzBl9fX3x9fdFoNFy4cEE+d1rvUW3KA7QhXsKZM2eSnp7OwYMH+fbbb7nzzjtR6oe33nXXXbi6urJ9+3bOnj3LIX01X20fVu3t7Wsks1YXOR06dOCNN94wEhWXL19m6dKlAEyZMoW//vqLPXv2MGDAAAYOHEi+/j3tWjg6OhqJoqu/bGxsOH/+POXl5bJQAV0H7sOHD8thoo0bNxqFjHbv3s1jjz0mr7Njxw4+/vjjGmFmgWkwm3ixsbEhNDRU/gU2cODAASIiIuq9d8OGDdx666188MEHzJs375p72dnZ4erqavRlCWw9PLBxcwOsO+/F4HVBoZArpK6FfVAvAMoSW6N40X0qs6/mwlaoVDjqq46K9u5t9h7u7u7ccMMN3HrrrfTo0cPoQZCbm8vWrVtZu3atEDECGTs7O6Kjo3nuuee4dOkSBQUFPPXUUzWu++677/j777+prKxk2bJlHDx4kFtuuYXAwED69OnD448/TkFBAZcuXeKpp55i9OjRtGvXrtH2tGvXjuTk+j2rLi4uTJs2jTfeeINt27Yxd+5c+dylS5fw8/MjICCAyspKXnnllTrXCQkJYf/+/Rw/fpzKykpWrlzJrl275PMPPvggH3zwAUePHsXR0ZGysjK+/fZbNm3axP79+3n55ZfJyspCpVLh7e1NYWEhZXpvtykSdn18fFCr1fz2228AFBYWMnfuXG677TZ69NC1jrg636WwsLCGB+eBBx6Qk7UFpsWsCbv33HMPq1evJjFR98l3z549bN++nXvuuUe+5quvvpJjwwCbNm1i+vTpvP/++zz44IPmNM8sGMIS2tNpFrakbgz5Lio/P5QN7AFhHxwMgObUKasPiV1NbeIFwElfUVG8d6/JBIWbmxsjRoyoVcRcunRJFjGnT58WIsZM2Hp46KrJakGhVmPr4dHCFtXN559/jr29PR06dKBnz55ERUXVCLXfdtttzJ8/HycnJ55//nmWL18uh5dWr15NQUEBvr6+BAQEyOGVpvD444+zdu1a7O3ta602MnDXXXfx66+/MmjQILp16yYf/+CDD/juu++ws7MjMDCQXr161bnG8OHDue++++jfvz+enp78+uuv3HTTTfL5OXPm8Morr3D33Xfj6OhIcHAwaWlp3HTTTURERKBWq+nXrx9qtZp58+bx9ddfy4LNFAm7zs7O/Pe//+XJJ58kJCSEzp07ExAQICcGnz9/npMnTxrl2tTmwVGpVDg6OjbLFkHtKCQzvoNKksTDDz/MV199hb+/P2fOnGHBggVGinzRokV8+OGH5OXlAeDq6kpVVRXB+oelgXvvvVfO2L8WBQUFuLm5kZ+f3+JemKwXXyR/9RpcJ07E7713W3TvhnLh3Xe5/PU3OA0bSqcGfkKpLCjgZL/+APj/vEKu1rF2Ki5e5NTQYQAE/LpODn8BaFJTSZugG8zYdcvmf7okm5CCggLi4uI4efJkDbHSpk0bIiMjCWhAubqgcZRnZdVaVWTr4YHqqiKC5lBRUUF5eTkO+gaVoKuasbW1lfMuJEmiuLgYZ2dno/+vTmVlJTY2NoDu4evg4IBSqWTEiBHccsstPPzww1RVVckhmqupqqpCoVAY/R5dvVdZWZn8QDVQXFyMvb29vDfoCiCqqqqMXtPVFBUVoVarUdciErVarXy8qKgIJycnFAoFpaWlqFQqbG3/qRMx/E0oFArKyspQKpU11qz+s7maiooKo/VAV36uUqnqzXtpKIWFhaSmptKlSxejZ8n//vc/Nm/ezNq1a+u8t67XI6ifhj6/zVZtBLpfyE8//ZRFixZx7tw5/P398bjqU8+9995rNNl327ZtNbLyQRcDbQ3YddV9EtGkplrYkrrRpuk8L3YNqDQyYOPqirpzZ7RnzlB67FirES+lhmRdOzvsroo9q7t0wbZdOyouXKB43z6ziBdXV1eGDx9OREQE8fHxJCcny2/Yly9fZtu2bXh4eBAVFSVEjAlR+fqaVKTUha2tbY2H59UP/erN4upqHFf94VxX9926hEtd567ey97evsY1te3VkI689TW/q/6wrn5dbWKo+u97bfYBdQoXoMbPHjCpp8PFxYXw8PAaxzds2MDUqVPrvbeu1yMwDWYVLwa8vLyMEp+q06FDByNh0k+fh9BasdPPzdGmpSFVVqKo5w/PUvxTJt3lGlcaYx8cjPbMGcqOHb/2xVZCmT5kadczUE7WNaBQKHAaOJD8deso3rsXj1mzzGaHq6srw4YNIyIigri4OCMRc+XKFVnEGDwx9T2oBAKBZVmxYoUQJxZGvEOaGDt9DFjSaKxyCrOk1ermGtGwHi/VMeS9lB07ZnK7zIWh0sihjn4TToP1eS/7DyC1QGt0FxcXhg0bxsyZM+nVq5eRSLly5Qrbt29n9erVpKSk1OqBFPy72Lx5M/fff7+lzRBchaOjo/iAYWHET9/E2Pr4oNS7La0xdKQ9exb0E2sb0l23Og7BOgGgSU2lqhG9ISxJXcm6BpwGDwaFgqrCQkrj41vMLhcXF4YOHcrMmTMJCgoyeiPMy8vjzz//FCJGgIODg0lyNwSC6w0hXkyMQqFArfe+aFKsT7xoTp0CwMbdHRtPz0bda9crCBQKqKpqFc3qKi5dokLfyrwu8WLbpg32obomU0U7d7aYbQacnZ0ZMmRIvSJm1apVnDp1SogYgUAg0CPEixmwk8XLKQtbUhPNyZMA2F3VUK0h2Dg7odYPJ2wNoSO5s65aXSNZtzrOw3TVSEU7d9V5jbmpLmJ69+5tlKSYn5/PX3/9xapVqzh58qQQMQKB4F+PEC9mwPCg1Fqx58Wue/cm3W8IHZUdt37xUnrkKAD2vXqhqMf17jx8BKD72Zhy4nBTcHZ2ZvDgwcycOZPg4OAaImbHjh2sXLlSiBiBQPCvRogXM2CoONKkpSFZ2QNGc1IvXno0TbzY99Yl7ZYmtALxos9hcail1LE69kG9sPHShdCKdlnO+1IdJycnBg0aVKuIKSgokEVMcnKyEDECgeBfhxAvZkCt7/UilZVRnplpYWv+oaqsDG2GbuZSUz0v9iE68aJNS6PSQtO7G4JUVUXpUZ3nxSG8/p40CqUSZ30ju6IdLZ/3Uh8GETNr1ixCQkJqiJidO3fy888/k5SUJESMQCD41yDEixlQ+bZHYag4SkmxsDX/oElNBX1vEbtqbb0bg33v3nIIpvTIEZPZZmq0aWlUFRYC1/a8ADgPHw5A8f79Vjn+wNHRkYEDBzJr1ixCQ0ONRExhYSG7du3iySefpGfPntjb2xMWFlZv90+BQCBozQjxYgYUSiV2+sRWrRWVSxvyXWx9fOQBko1FqVbLlTulcfGmMs3kGISVrbc3tj4+17zeafAgUKmQysoo3rfP3OY1GUdHRwYMGCCLGEOH0djYWD788ENOnjyJRqMhISGBadOmCQEjEAiuS4R4MROGpF3NKSvyvDQzWdeAg34qeGl8XLNtMhfV810aUlVl4+yMU3/d7KbCbdvMaZpJqC5iwsLC2LRpEwqFQu7aK0kSCoWC559/nkp9Xx+BQHB9UZ6VRenx4zW+LF140BK0yHiAfyN2+rHpZfrSZGvAdOIlHL6F0vgjVjsCoTRe53lpzAwml1GjKN6zh6LtfyItrqgxTsAacXBwoH///ly8eLHG4EdJkkhNTWXFihVEREQQGBhY75wYgUDQeijPyiJ17DikWiZoK9Rqum7Z3CLzvSyF8LyYCbueuunFmpSUWn+5LIHBC9Rs8aLPIakqKZEFkTVRWVgo5xo5RIQ3+D6Xm24EhYLKvDxKYmLNZJ15CAwMrOFhUgA+Pj4UFxezZ88eVqxYwfHjx6moqLCMkQKBwGRUXL5S57NF0mprnap+PSHEi5mw79lT9z/l5WjS0ixrDLoHesX580DzxYvK2xuVnx9Ai7bUbyhlCQm6xGRbW+yDghp8n62nJw6RkQAU/vGHucwzCwsXLtSFivTfKwAJmNP3n0GnxcXF/P333/z8888cO3ZMiBiBoJVSvP8AmU89ZWkzLIoQL2bCtm1bbL29AShLSrKwNdVybxQK7Lo2bpp0bch5L3HWl/dSohdU9r16oWzk5FeXkSMBXd7L1WEYa2by6DF81KkzPezssFOp6Nm2LR/7+nH3lctEBgUZzccpLi5m7969rFixQogYgaCVcfn7H8i46y7K09MtbYpFEeLFjNj10nlfNCesQLyc1M0iUnXqiNLBodnrGcIxJVZYcVSqD/k0Jt/FgMsonXipyM5uFSMQDBTv2skoBwfW9exFSV4eR44cYZSXF5Xns+l04CCzZ88mMjIStVot31NSUsLevXtZvnw5CQkJQsQIBFZO3rp1XHjjDZAk7AIDLW2ORRHixYzY9+wFWIfnpex4IkCjwij1Ych7Kc/IoCI31yRrmgKpvJwSvTfIsV/fRt+v7tABuyDdv1vB5i0mtc2cGGx1HjoUpaMjKj8/PB98EIBLS5dCRgZ9+vRh1qxZREVFGYmY0tJS9u3bx/Llyzl69KgQMQKBFaJJSyN74SIAnIYOpf3ixfVeL5VaX78qUyLEixmx1yftliUlWTwEUZZoWvFiHxiIUt+Ir+TwYZOsaQrKjh9HKikBwLFv48ULgNv48QAU/Pab1Y13qI3KggJ5IrbrhAny8bZz70TdtSuUl5P13PNIFRXY2dkRFRXF7Nmz6dOnTw0Rs3//flnElJeXt/hrEQhMyfVSSixVVZH17HNIGg2qzp3o8OEH2Hp7oaj293s1eevWtZyBFsD6a0FbMXb6pN2q/Hwqzp+3WNmapNXK06RNJV4UtrY49ImieNduig8cwHXsWJOs21yKDx0CwK57N2w9PJq0huv48eS8t4SK7GxKDh/GqV+/a99kQQr/2Iak1aJ0csJ5xHD5uEKtxveN10mfNZuyY8e49NXXeM67HwC1Wk1kZCTBwcEcO3aMhIQENBoN8I+IiY+PJywsjKCrcmYEgtbA9VRKnL9+PWVHj4JCgd/bb6N0ckLp5ETXLZtrVBXlr9/Ale++I3/dOtreNVfuOXa9ITwvZkTdqZM8JqAsKdlidmhSUpD0n6JNJV4AualbyYGDJluzuZQc1IkXx75NFxwqX18c+kQBULBxk0nsMicFmzYCuj41VycoO4SF0faeuwG4+OmnNfoOGUTMrFmz6Nu3L3Z2dvK5srIyDhw4wPLly4mPjxeeGEGrofzCBfI3/XZdlBJXaTRc/PAjANymTDEad6Ly9cWhd2+jr3ZPPanzuFZUyPddjwjxYkYUNjbYG5rVJZ2wmB2GkJGtb/smeyNqw7GfTrxo09Ioz8kx2bpNRaqooDQmBmhavkt13CZOBKDg99+tpk9PbZTn5FC8/wAArnqbr8bz4YdRd9OFj84/97wsZKujVquJiIhg1qxZ9OvXr4aIOXjwIMuWLSM+Ph6tFf88BP9uyi9c4Nz8+aTccCMXlyyxtDkmIf+XdVRkZ6NQq/F67NFrXq9QqfCe/ziga/lgTY1STYkQL2bGGiqODOLFQT+TyFTYB/VC6eIC/OPxsCRliYlUNTPfxYDLmDFga0tVfj5Fe/aYwjyzULh5M1RVYePpidOA/rVeo7Szw/eNN0CppOz4cXK//LLO9dRqNeHh4cyePZt+/fphX82To9FoOHjwIMuXLycuLk6IGIFVUbx/P2kTb6Zw8xaoqkLh5GRpk5qNVFnJpW++AcBtWjSqBsxpA3C+8Ua5GunS/74ym32WRIgXM2MfqBMvZScs53kpPX5cZ4sJQ0ag8yw59ukDQMmBAyZduymU6PNd1F27Ytu2bbPWsvXwwHnwYEAXQ7ZW8jfoQkau48bVO87AITSUtvfeC0Dup59dM8lapVIRHh7OrFmz6N+/fw0Rc+jQIZYvX05sbKwQMQKLU7R7D2fv+w9VhYXYennh9/4SOn/7jaXNajZFO3ZQnpEBSiVt7767wfcplEr5771gyxarqgg1FUK8mBn73jrBUH7unEVirFJFBRp9vo2pxQuAY39dbknxQSsQL4Z8l2aGjAy4TZ4EQOH27VRcvmySNU1JWVKS3IvGbdLN17ze6+GHsA8LhaoqMp98qkG/jyqVirCwMGbNmsWAAQNwqNYjSKPRcPjwYSFiBBal7ORJMh9/HKm8HLvAQALWrsF1/Hi4xhwvqcz6S4mv/PwzAM433IC6Y8dG3es6ZjQ2np5QXk7e6tXmMM+iCPFiZuwDA+VytrKEhBbfX5OWhqSvIjGHeDEk7ZafyaBcP37AEkharex5cWpmyMiA88iR2Li7Q3k5+b+uN8mapiRv5SpAV9VmHxx8zesVajV+S95H6eJCxYULZD37bINLwVUqFaGhofWKmGXLlhETEyNXLQmsm+uhjLiqtJTM+U9QVVyMqkMHOn3zNbZeXoDOe1pfKfHl775vKTObRHlmJsW7dSFrjxm3Nvp+hVqN+/RbAN17RWto+9AYhHgxMwq1GvteuqZnpUeOtvj+crKul5f8R21K7AIDsXFzA5ATRy1BSVy8Lt9FqcRx4ECTrKlUq3GbPBmAvFWrLN6rpzpVpaXkb9CFs9yn31JjKGNdqDv40f711wAo3rmLy98ubdS+tra2sogZOHCgkYjRarXExMSwfPlyDh8+LESMFWMoI06fdkuNr9Sx41qNgMn54AO0qakoVCo6fPqJUbhY5etL1y2b8V+z2uir7X/+A0Dh1q0U/vmXpUy/Jvnr14MkYevbHid9CLuxuEdHA3qhqi9muF4Q4qUFsA8LBaD0qAXEi4k7616NoppYKLZgYmvx7l0AOISEmLSiyv3W6YCuoqo01nomTRf8/jtVhYUo7O1xu/naIaPquI4ejcecOQDkvP8+xfv2NXp/W1tbQkJCmDVrFoMGDcJR3xIAdCImNjaWZcuWcfjwYcpagXv+ekeqqqJ43z4uvPU2GXffQ8Z/7m/1ZcRlJ05w5cefAPCaPx/7Wtrl11ZK7DX/cZxHjAAg+5VXqCoubkmzG4QkSXI+m9vNk1BcIwRWF+qOHeW2D3m//moy+6wBIV5aAIdQ3YydsqNHW/zTe+nRIwANCis0FedhwwAo2rMHyUKt5Yv07lWnYUNNuq5d167ypOm8VdYTNzaEjFzHjMHG1bXR93sveBr70FCorOTc4/PRnD7dJDtsbW0JDg5m5syZNURMeXk5sbGxLF++nEOHDgkRYwEkSaLg962kTZhIxl13c3npUor37kWbkmJp05qFJElkv/Y6VFVh16sXbe64vcH3KhQKfBa+jMLRkYrsbC59860ZLW0amhMn0KalAeA2ccI1rq4ft0n63L0tv1N1HeWlCfHSAjjoPS+V+fmUnznTYvtWlZVRlqircnKIjDDbPs5Dh+j2y8+3iHep/MIFNMnJeltMK14A3KfrvC8Fm2t2s7QEZcnJshfI4BlqLEo7Ozp88l9s27WjKj+fcw88SGV+fpNtqi5iBg8ejFO1MtXy8nLi4uJYvnw5Bw8eFCKmhagsKCDz0cfIfOwxtHpx6tinD23/8x887rzTwtY1j+Jdu+QwiM+LL9RbaVcbqvbtaXvvPQBc+uYbKi5eNLmNzcEwq8wuMBC77t2btZbr6NG6tg9FRZQ0wctqrQjx0gKoOnTARh/KaMmHe9mxY1BeDgpFkyYsNxRbLy/Zs1O0c5fZ9qkLQ7jKxsPDLB4m13FjsfHwQNJoyFuxwuTrN5bLS78DdG9sBq9QU1B5e9Px889QODigTU8nc/78WhvYNQZbW1t69+7NzJkzGTJkSA0REx8fL0RMC6A9l0n6jJkU/vEHAC6jR9Nl8290/vEHvJ+Yf83qtIrs7JYws0lIksTFjz4GwHnECByjopq0Ttu5c7H18kIqLbU670vRjh2Armt2c7Fxd5dHnBRs3drs9awFIV5aAIVCgUOoPu+lBZN2DdOV7Xr0wMbZ2ax7yaEj/YDAlqRo124AnAYPRqE0/a+00t4ej1mzALj80zKLul7Lc3LI36iLhbeZO7fBibp1YR8UhO/bbwFQvHcf51962SRVCTY2NgQFBdUrYpYtW8aBAwcoLS1t9n6Cf9CePcuZ229He/o0Cnt7fN99lw4ff4RdQECD1zj/0sto09PNZ2QzKDlwQC5EaEjH2bpQOjrK3pcrK1ZYhVcVdMJTc+oUAM43jDDJmi6jRwNQtP1Pi4X2TY0YzNhC2IeFUrRzZ4t6Xkrj4gFwiAg3+17OI4aT+9lnaJKSKL9wAVW7dmbfE3R9bIr37tXZYOJ8l+p4zJ7Fpa++ojI3l4ING3GfFm22verjyrJlUF6OrZcXbhPGm2RN19Gj0T75BBeXvE/+unUoXV1o99xzzRZG8I+ICQwM5OTJk8TFxVFUVARARUUFR44c4fjx4wQFBREWFmZUvWQqKnJzKT2aQPnZDMqzsqgsKNR5mBQKbNzdsW3jgTqgC3Y9eqD27ywL4PKsrFofaLYeHlY70K+ysJCz98+j4vx5lK6udPrfl7V6XQ1lxHUl7VZevkzGvffhv3yZWaoUm4OhxNlp0EC5krOpuE+fTu4XX1J5+TJXflqG18MPmcLEZlG0cwcAtt7eJiu0cBl5E9mLFlGZl0fp0QQczZhG0FII8dJCyEm7SUlUaTQoq82OMQeSJFGq97w4Rpj/F9U+OBibNm2ovHyZol278JjetFyMxlJyOIaqwkJQKJpcTtgQbD09cZ10M/mr13B56VLcoqea5OHeGKpKS8lbrgtbecyZU28Pi8bS9t57qbx8hcvffsuV739A6eSE16OPmuw12tjY0KtXL3r06MHJkyeJj4+nsLAQ0ImYo0ePkpiYSFBQEKGhoUaJv41Fqqig5HAMhVu3UrR7N+VnzzbcTnd3nAYNwj4slItL3m9VE4mligoyH5+PNi0NhZ0dnb76n+zxvRpDGXFt4qzy0iUy5z9B+blznHt8Pp2/W9ronBJzoU1Pl0MqHnfc0ez1lI6OeMyZTe5/P+HKihW0/c99KE34d9UUiv7aAYDz8OEm+/uz9fTELqgXmsQTFO/Zc12IFxE2aiEcQkN0/1NeTpm+Xb85KT9zhkr9G5NDC4gXhVIpJ8sWtWDvhILfdYltjlFRzR4JcC3a6pMcNadOyW+gLUneqtVU5uejcHDAY+YMk66tUCjwXvA0bnqP0qXP/4+L779v8uo4g4iZMWMGw4YNw0U/Gwv+ETHLly9n3759lOjnVDWUikuXyP2//yNl5Cgy5s7lyrJlsnBRODhgF9QLl1EjcZ9+Cx6zZ+E+YwYuY8fiEBWFUt+rqDIvj4LffiPnzbdaXSlxzpL3Kf77bwB833qzTuFioLYyYofevXEeNgy/jz4ChYLSmBgu/veTljC/QVz+8SeQJNSdO8uh6ubiMWMGCpWKytxcCrdsMcmaTaWquFgetWKqkJEB58G6wgrD70hrx+ziZdmyZQQFBeHk5ER4eDi//fabWe6xdmxcXeVBWS0xxLBEHzKy8fRE1ci20k3F+aYbAV0CbWVBgdn3kyorKfxjG6AfpGhm7Lp3x/lG3Wu8+N//tmjHyqrSUnmgoset03Wdf02MQqGg/Suv4DZ1KqAb6HbhzTfN8jqVSiU9e/ZkxowZDB8+3EjEVFZWkpCQwPLly9m7d+81RUzFlStceOddUm68iYsffqRLNlUocOgThfezzxCwdg2Bhw7SZe1aOvz3v7R/9VV8Xn6Z9osX0eHDD/D/6Ud67N9Ht+3baP/667iMGwsW/vTdWIoPHuTyt7qkU8+HH8Z13Lhmrec8ZDCeD8wD4NKXX1K0x/IPvMrCQvLXrgXA447bTZbfZuvpiet43c/L0u0QivftQyovR6FW4zRggEnXNnimSxMSmlVZaC2YVbxs3bqVO++8k6eeeor09HRmz57NlClTiNOHM0x1T2vBUZ/xXdICc4D+CRmFt1h4w3nYMJSOjkjl5RRu/9Ps+5XGxlKpHzjmMrr5WfkNwevRRwDQJJ6QhVNLcGXZMipzc1E4OMgdQs2BwsaG9q+/hvutunbkV77/Qdd+3UyVQUqlksDAQFnEuFbrWVNZWcmxY8dkEVN8VTMxqaqKK8uXkzpyFJe/+QZJo8GmbVvaPjCPbn/9if+PP9J27lzsg4KuGfZQKBSo/PxwnxZNhw8+oNM3X5vl9ZqDqpISzr/wIgAOkZGy6Ggung8+qBu8KklkPfcslfown6Uo+G0zVSUlKJ2dcZ8yxaRru03VeRxLDh1Ce+6cSdduDIV6j67jgP4omxE6rQ3HyAgUjo6gb1jY2jGreHnvvfeYMGECd999N15eXixYsICwsDA++OADk97TWnDSDzEsiY2r0yVtKkrjdH1AHMJbLraptLfH+aabACjYbH5vWcGW3wHdG3ZLJQjb9+yp+2QO5H7yX6TKSrPvWVlUJI+1b3P77WYPjymUSnwWL5IrMQp//52MuXdRnpNjtj0NIubWW29lxIgRuOnDOPCPiFmxYgVvvPEGISEh2NvbE+TpyQ9PPU1VcTE2Hh54P/sM3f7cjvdjj6Hy8WmePddIHC7eaz1v/jlL3qf87FldZdEbrze5G+vVKGxt8V3yHkpnZyov5nLxv/81ybpNpUA/DsN13FiU1arXTIFjv75yDlO+BTvRluh/r5yHmiYkVh2FWi2XTBdZsBu6qTCbeJEkib1793LDDTcYHb/pppvYq68OMcU9rQnHPn1AoUAqK6PUjEMaK3Jz0ZzSddB0jGp6H5CmYHBXF+/dZ9a8AKmqikJ9zwLXMaPNtk9teD38MCiVaE6lyPOFzMnlb5dSmZeH0tmZtnffZfb9QJ8D89RT+CxaCEolpfHxnJ42jZLDh826r1KppEePHkyfPp0bbrjBSMQcOnSIF154gePHjqHRaEi+coXHsjLZGxJC1z+20nbuXLMnwhu4uGQJmU8+RaW+cspSlCUm6irQAO8n5qP29zfp+qp27fB67DEArvz4E2UnTph0/YZSnpkp/+65NnIcRkNQKJW4TdHNMctf96tF5phpz2XKM6Uc9R90TY3TEEPey16rmtXWFMwmXgoLCykuLsbb29vouJeXF9l1NEBqyj2gm2pbUFBg9GWN2Li7Y9ezJwAlBw+abR+DS1Dp4mLWsQC14TRkMEoXF6iooHCb+cIqpXFxcldMQw+DlsKua1d5YGPOe0vMmt+jzcjg0v/+B0Cbu+8yS65LfXjMnEnHL77Axs2Nyou5nLlzLhc/+RSpvNysU4mVSiXdu3c3EjGb1q9HARjeciV0IuvDpBMm72NU70Rifa5FwaZNpN8ynbKTJ026d0ORJIkL77wLkoRdUC88brvNLPt4zJqJXa9eUFVF9qLFFplOnL9xEwC27dvrPgSaAcPfdPnZsxYZYlhySJcLaePujl23bmbZw3mILu+l4vx5efxAa6XFq42USmWjFd+17nnzzTdxc3OTvzq2UIJqUzC47YoPmFG8/K3zUjkN6N/iJY5KtRqXkSMBKNy82Wz75P3yCwAO4eGo2rc32z514f3EfJ07PTeXix+bx50uSRLZr7+OpNWi6tiRtnffbZZ9roXz0CEErF2jE8KVleR+8glp0dGkjhlr9qnESqWSbl26MFKrJSc7m6vfBSRJIikpid27d8v9Y0xBXROJ/despusfW2n/5pso7O3RpqeTPmMmBfpOti1J0Y4dlOzfD0C7Bc+YpUEj6MJH7Re+DEDpkSMU/v67WfapC0mSdBOW0c35MdfrVHfu/M8cs3XrzLJHfRjEi2PfPmZ7jarOnVH5+QFQrP/daa2YTby4uLjg6OjIxatmRuTk5NCujvyEptwD8Nxzz5Gfny9/nW1EX4eWxuAOLI2LM0unVkmS5KZtToMGmXz9huA6Xtc8rXjffsozM02+fmVRMQW/6YSRm4Waxdl6eeH1+OOALpm21Azl74XbtlGsH7fg89KLKO3tTb5HQ1H5+eG/7CfazrsfbGzQnkqpc5SAKUuJiw8eJH3mLHJefwN/lZqrU88VCgU+Pj6cOHGCFStWsGvXLrl/THOpq5RY7eeH+9Qp+P+8AlXnTkilpWQ++hiXv//eJPs2BKm8nJx33wPA+YYbcBrQ36z7OYSH4zpBNyDw4sf/bdEurZoTJ9CmpgLmCRlVx+B9KfxjW7NHZTSWf8RLX7PtoVAocNRPmS6Nbd1FMGYTLwqFggEDBrDzqnbxf/75J4PqeKg25R4AOzs7XF1djb6sFcc+fUCpRNJoKDtyxOTra1NSqNAnVlpKvDgNGoitb3uQJK6sWmXy9Qt+24RUUoLS0RG38abpMtsUPGbNxC5I504///wLJq3Iqbh8mQuvvAqAy6iRJutp0RwUajXejz+O/4rlqEycW1EdSZIo3r+fjLvvIeOOOylLSACFggVTp8qhIvT/lSSJiRMnAlBVVUVSUpIsYswdPrYPDCRg5Urdw0aSuPDGm1x4990WySXI37BR5/a3tcX76afNvh+A58MPgVKJ9vRpeURFS1Cgb5VhFxiIfY8eZt3LZdRIUCqpys83e35Xdcqzs+WeROYULwAOETrvUom+qKO1Ytaw0RNPPMH69etZtmwZhYWF/Pe//yU2NpbH9AlgAK+88gqenp6Nuqc1Y+PqKre0LjJDsyCD10Xl54eqUyeTr98QFDY2eOhLbfNWrzF5ZZWhF4PrhAkmrzpoDAobG9ovfgVUKjTJyVx4+22TrCtVVZG14BkqLl5E6epKu+efN8m6psIhJAS/d9+p9xqpCUJOe/Ysl77+hrSJN5Mx9y75d9lxwAD8V/7MPSuWs2bNGkJDQ7G3tyc0NJQ1a9bw9NNP06ZNm3/21oeSfv75Z3bu3GlWEWPj5kbHr7/CVS+gLn/9DRfeeNOsAkaqrOSSvueP2+RJ2HVp+Myi5mAXEICbvkQ599PPWswzUahveuk61vy9nGzbtMFB3322cNt2s+9nwOB1Ubq6YmdmgWZ4fRVZ5yk/f96se5kTsyZETJgwgS+++IIXXniB22+/ne7du7Ny5Ur6VlOWVVVVVFRzQTbkntaO8/BhlB0/TtFfO/DWhx5MRVG1kFFLt6+vjvu0aVz85FNd18rt25vdNMtAWVKS7pM44H5ry4wgqA+HkGC8n3yCnLfeJm/5Cpz698d17NhmrXnp66/lSdm+b7xukZyea3KNmPyZO+fiGBWFQ1iYbmZQp47YeHjo+gBVVFBVUqJLGjx7jrJjCZTExaFNSTVaw3HgADz/8x+cBg6Uj0VHRxMdXTNUGBAQQHp6OjExMVy+fBnQiZjk5GROnjxJ9+7diYyMNItXVqlW4/vO2yhdnMlbvoIrP/yAVK7FZ+FCs/wNFv7xh25oolKJ5333mXz9+vB88AHy16+n/OxZ8tevx33aNLPupzl9Wk4sNTSINDcuN42k9HAMhdu30+7FF1rkfdTQuNSxTx+TlbrXhV23bihdXakqKKAkNhY3fTiwtaGQLFwvVVVVhSRJ2JjwH6ygoAA3Nzfy8/OtMoRUmnCMdP3sn67b/kDdoYNJ1pW0WpIHDEQqKcHvww+a/RBtLufmz6dw8xYc+/Wj8/ffmWTN84sXk7d8BXY9exLwy1qLCjQDkiRx7oEHKdqxA6WjI51//KHJA9UKd+zg3EMPQ2UlHnfcjo+VeV0MlB4/Tvq0W0y+rm379riOHoX7rbdi17Vro++XJIn09HRiY2O5dOmS0TmFQkH37t2JiIgwKsE2FZIkceHNN7ny/Q8AtL3vXryffNLke5yeGo0mKQnX8ePxe3+JSddvCOdfeom8VatRd+tKlw0bzPo3eOmbb8l55x1dAvX2bS3y9649e5bUUboKRv/Vq3EI7m32PVPHjkObno73M8/Q9q65Zt8v4/77Kd65C485c/B56UWz79cYGvr8tvhsI6VSaVLh0hqwD+6Nrb4cvOhP03WiLTl8GKmkBJRKHPubN4GvIXjMnAXoysJNUU5acfky+b+s0689wyqEC+jb6r/5BqqOHakqKSHjP/dTvHdvo8uIiw8eJHP+E1BZiUNYGN5PPdWCr6Jx1FdKrFCraf/mG7SZOxfHgQOwqRbSMbpOpULVuRMuo0bhvWAB/qtW0e3P7bR77rkmCRfQ/VsEBAQQHR3N6NGjaVutoZ8kSZw8eZKVK1fy119/kW/iFukKhYJ2zz2Hx5w5gG68Qq6+zN1UFO/ahSYpCYC295uv03J9tJk7FwBtSirFZh4bYHh/dL7xxhb7e1d37CiPcincZv4qsoorV3SeNDBbGfjVOBryXmJbb96LdYwK/ZehUChwvvEG8lb8TOGff9HGBNNR4Z8YrWNkJLYeHiZZszk49uuLXffuaE6d4tL/vrpmnsS1uPLjT0hlZdi0aSPH3q0FWw8POv3vS9Jnz6EyN5eMu++p9bq6JhIXbttG5lNPI5WVoQ4IoMP/fW7x6bb1Ud9UYlsPjxqvT9JqqcjLQyotRaFSobCzw8bDw3zlvQoF/v7+dO7cmYyMDGJiYsjVj5KQJIlTp06RkpJC165diYyMxN1E/XMUCgXtXnieysICCtZv4OKS91H5tMft5okmWf/ydzoPpvOIEdjrH7AtjV3XrjgNH0bxzl1cXroU56FDzLJPxZUrlOjHnLjceMM1rjYtLjfdhCY5maLt200e2r8aQxhcoVZj37Nl/k0NzUs1yclUFhWZvE9SS2Bxz8u/FRd9G/2SQ4dMMiRLkiQKt+vEi/PIm5q9nilQKBS60lp0Db0Mny6aQmV+Ppd/0LnjPW6bY9Gy4bpQ+/vT+bulKOt5EF5dRlxVWsqFt9/h3MOP6IRLly50+m6pVYjPa1FXKfHVwgV0b8wqb2/UnTuj8vXFtm1bswkXo30VCjp37szUqVMZM2YMXl5e8jlJkkhJSWHVqlX8+eef5OXlmWZPpRLf11/HabiuQuz8Cy9QaoLKQk3aaXksQZs7bm/2es2hrd77Uvz332Zr0le8ezdUVqJ0dm4xj4QBl1G6XlWaUylmn3VUelQnXux79UKhUpl1LwP2ISGgUkFVFaXxpq96bQmEeLEQjv31g7cqKynatbvZ65UlJFBx4QKA3CTOGnAdO1bXsryqipyPPmryOpe+/oaqwkKUbm60ud2yb9z1Yde9O75vvln/RZIuETH3f/8jdcxYeRqw44ABdP7pR1RXdZgWNB+DiJkyZQpjx46tVcSsXLmS7du3c8UEPWoUKhV+S5Zg170bklbL2YcebnZlx5XlywFQBwTgWC2J2RI4Dhggh1YM3iBTY6gych42tO5ux2bCrmdPbLx0VbCGpp/mojThKAD2oaFm3ac6Snt7HPR5eaWtNHQkxIuFUKrVOA0dCmCSNvoFm/S9EIJ6mSwB2BQobGzkZm6Fm7dQGh/f6DW0587Jb5Ce992LjYuLCS00PbbeXvWeT58xg7Rx47m45H0qcnJQ2Nvj/dSTdPrqf63C49KaUSgUdOrUiSlTpjBu3Lgao0hSU1NZtWoV27Ztk6uWmoqNszMdPv8cGw8PKnNzOfvAg1RdNRm7oVQVF5Ov7yrtMXu2xfO9FAqFHO4u+G2zyWc8SZWV8pgT5+HDTbp2Q1AoFDjr+2QVm6GlhQFJkijTe14cQkPMtk9tOISHA5h1zp45EeLFgriMGgVA0V9/NWvcvFRZKTdycpto3g6UTcFlzGj5D+X8osWN6g8h6Zt/SRoNKj8/s81vaVH0k6htvb1pe+89dN36O23vvbfFRzn8m1EoFHTs2JHJkyczbty4Gh2809LSWL16dbNFjLpDBzp88l9dL6CkJM6/9FKTesDkb9hAVVERCkdHeYCgpXEdNxalszNSaSkFJm5aV5Z4gip9ON1SXiZDk8/i/fvNNj2+PDOTSr2nzyGkZcWLvb6Kquz48VY5pFGIFwvictONKJ2ckLTaZs0LKTl4UDekUKHAdbxp+qmYEoVCgc/LL4GNDZqkJC599VWD7y3YsEGuOGj3wvNWmevSWHwWL6brtj/otnMH3k89JcJEFsQgYiZNmsT48ePrFDF//PFHk0WMY1QU7RctBHReirwmdJ2+8vNKANwm3Ww1nkeloyOuE3U9QgyNI02Fweui7tbVYn8fBtFUVVBA2bFjZtnDkKyrdHVF1bmzWfaoC/veOvFSefkyFfUMPrZWhHixIEoHB1zG6LpGNmcQWN7qNQA49uuHysfHFKaZHPugINrefRcAF//7Cfm/bb5mKbEmNZVsfYt81/HjcL6hZSsOmsq1yoidhw5B3aGDxV3/gn9QKBR06NCBSZMmMWHCBHyu+js6ffo0q1evZuvWrTX6xzQE92nT5Lk5F15/g7Lkhie5lp04gebECQC5c7W14K7vV1V2/DhliYkmW7dkv068OA20zIgTAJW3t9zt1tDt2dQYknUdQkJa/P1A7e+vy7tE9+/X2hB+agvjPnUK+WvXUno4Bs2pU9h1796o+yuuXKFw61bdWtMt33G2PrweeYTigwcpO3KUrCeeqPUaQykxNjacvX8eVUVF2Pq2N1u3UnPQ2DJigfWgUCjw8/PD19eXrKwsYmNjOV8t0TY9PZ309HT8/f2JjIw0Gm1yLXxefonSo0fRnj5N5vz5BKxeJT886iNf/8HGrkcP7PSjRawFh969sQvqhSbxBFdWraL9woXNXrOqrIySwzEARt2VLYHT4MFoTp6k+O+9eD7wgMnXl5N1Q4JNvva1UCiV2AX1ovRwDKXHj1tVoUdDEJ4XC+PQpw923bsBcHnZskbfn792LVJ5OTbu7riMHmVq80yKQq2mw4cfYlOtcdjVSFotxYdjODNrNuXnzqF0cqLj559jY4aOqOakMWXEAuvDIGJuvvlmJk6cSPurRjSkp6ezdu1afv/9d7l/zLVQOjnh9+EHKNRqtGlpZL/62jXvkcrLyd+gyydxmzrVKgW8h/5DU8GGjVRpNM1erzQuTjcPzcYGx36WHQtjyHspiY+nsqhpydZ1IVVUUHZc561yaMFKo+o49P4n76W1IcSLhVEoFHjMng1A/q/ra/20XhdVWi2Xv/seAPdbpll1UzMDqvbt8Vn4cr3XnH/2WcqzslC6uNDxi/+zWDMugQDA19dXFjG+V4nPM2fOsHbtWrZs2cLFixevuZZ9YKA8aDP/l1/k3kx1UbR7D5WXL4ONjcka3Zka1wkTUKhUVBUVUbRjZ7PXK963H9CFUizdPM2xT5Su90pFBSWHDpp0bU1qKlJpKQD2wS3veam+b9nxxFaXtCvEixXgNmkSNu7uSCUlXP7++wbfV7Bhg67UVqXC43bTdOltCVR+fvVfUFWFXa9e+C/7qcWbUwkEdeHr68vEiRO5+eab8bvqdzgjI4NffvmFLVu2kJOTU+867jNuxXnECADOL1xU7wcWQ8jIecgQbBsRompJbFxdcR6hK2c2RdWRIVnXceCAZq/VXJQODjhE6KYwlxw+bNK1yxJ1eUy23t4WS0qWk3YvXZL7hLUWhHixApROTvK8kCs//Ngg70tVWRkXP/0UALcpk1G1u34qVjwfe5SAlT83Ov9HIGgJ2rdvz4QJE5g0aVKtImbdunVs3ry5ThGjUCjwWbwYpasrlbm5XHj9jVqvq8zPp+gvXaM2t6lTTPoaTI3rBJ1XqGjnTioLCpq8TmVhoRzCcBpg2XwXA459ogAo1efhmArDjCq7Xj1Num5jMEraNVNFlbkQ4sVK8LhtDjbu7lQVFXHxgw+vef3l73+gIus8Cjs7sySSWRLnYcNarE22QNBUfHx8ZBHT4arGkGfPnpVFzIVaPtGq2nnj8+ILgM5bUVujysLtfyKVl6N0crL6SjvnEcP/afvwR9ObbpbGx0NVFQqVCofwMNMZ2AwcovTi5fhxqvRhHlNQlpwMgH2g5cSLIWkXdK+vNSHEi5Vg4+yM1/z5AOStXMmVVavrLCPWpKWRq/e6tLnzzlaXBHqtUmLRZVbQmvDx8WH8+PFMnjyZjh07Gp07e/Ysv/76K7/99hvZV/XScL35ZpxvvBGoPXxU8PsWAJxvuhGlnZ0ZX0HzUdrby003CzY1PXRUEqPzbtgHB1vNa3YMDwcbG6ioMMmMKtA13zR4XlpqGGNdtNakXVEqbUW43zKNKyuWozmRRPZLL9U4r1Cr8V+7lqynn9Z1nO3cCU/94MPWhCglFlyPtGvXjnHjxpGTk0NMTAxnz56Vz507d45z587h5+dHVFQUPj4+uvDRooWkxcRQeekSOUuW4PuargKpMj9fnqnjOtb6Gk/WhuvEieSvW0fx/gOU5+Q0KY+jNFY3Rdow9dgaUDo5Yd+rF2XHjlFyOAanAc3PxanIyaFSPwjUrqflPC/wT95Lmb6XUGtBeF6sCIWNDV6PPlrneUmr5fyCBbqGVba2+L39doP6RFgjopRYcL3i7e3NuHHjmDJlCp06dTI6l5mZyfr169m0aRPnz59H5e1NuwVPA5C/eg0l+od34bbtUFGB0tkZpyGDW/w1NAWnAf2xadMGqqooukYVVW1IWi2lR3V9Txwio0xtXrNw1IeOSmJMk7Rr8Loo7O1Rt3Bn3asxNOKrvJjbqGpXSyPEi5Vhe41PK2WJiaBQ4PvWW/K8IIFAYH14e3szduxYpk6dWquI2bBhAxs3bqS4f3+5oiV78WKkigo5ZORy042togUCgMLWFpebbgKg8I8/Gn1/2YkTSGVlADhEhJvStGbj2FdX9Vgaf6RRs9nqoixJl+9i1707ChubZq/XHNRduujCYoDm5CmL2tIYhHhpZSgcHfH7+CPc9DNFBAKBdePl5cXYsWOJjo6m81WfsrOystj0228kjxgOSiWa5GRyP/uc4r26cmGXsWMtYXKTcRml69JafPCQHBZpKCUxsYBunpG15b0Zknal0lKThFc0ydaR7wKgVKtRB/gDoDnZ8LEVlkaIl1ZGh/9+jOso6+6kKxAIauLp6cmYMWOIjo7G39/f6NwZhYLMsDD+KCxkyNNPE554nKkZZ/i9AY3vrAnHAQNQOjtDRQWFO3Y06t7SOJ14cYywnnwXA7YeHqi7dgWg5FDzQ0ey58WClUbVsdeHjoR4EZgNG3d3S5sgEAiagaenJ6NHj64hYtY6OfJYViantBq0ksTJ0jKmz5jB2rVrLWdsI1Gq1TgP1zWsq638uy4kSZI9Lw5WlKxbHUPeS2l8XLPWqSorQ5ueDliH5wX+yXsR4kXQZEQZsUDw78AgYqZNm0ZAQADrt2xBARiatEtIKBQKXnrppVbVut1QMl28ew9VJSUNukebnq4bg8A/IsHacAjTzR8qPXK0WetoTp2CqioA7Kxk9IksXk6dQtLbZu2IUmkrQ5QRCwT/Ltq2bcuoUaPIzc3laokiSRInT55k/fr1REVF4efnZ5XDGavjPHQICrUaSaOhaPceXMeMvuY9pfG6/ik2Xp6ormr4Zy0YhidW5ORQfuECqnbtmrROmb7SSOXnh42Li8nsaw4G8VJVUkJ5VhZqK/03qI7wvFghooxYIPj3ERgYWEOYKBQKfHx8uHDhAr/99hu//vorZ8+etWpPjNLJCafBuvJuw3iDa1GWkACAQ0io1YozdZcuKJ2cAJrVrM5Q0WMtXhfQPXMMbTdaS+hIiBeBQCCwAhYuXIgkSfLDW6FQIEkSEyf+M006JyeHzZs3s27dOjIyMqxWxBgGTxbt3t2gMESpQbyEhpjTrGahsLH5Zwrz0aaHjrSpKQDYdetmErtMgUKplGfJCfEiEAgEggYTHR3NmjVrCA0Nxd7entDQUNauXcsbb7xBt27djDwSFy9eZMuWLVYrYpyHDQV004qv1Xa+SquVQyn2wdYrXuCf0FHp0YQmr6FJTQPArltXk9hkKux6tC7xInJeBAKBwEqIjo4mOjq6xvEbb7yRyMhIYmNjSU1NlcWKQcR4eXkRGRlJp06drCLsomrfHrvAQDTJyRTt3IVDSN2iRJOcDPrGbw7BvVvKxCZhr/cMlR07hlRZ2egGc5WFhVToB3UaSq+tBbvuuryXslYiXoTnRSAQCFoB7u7u3HjjjUyfPp3u3bvX8MT8/vvv/PLLL6Snp1uFJ8Z52DAAinburPc6w0gAVedOVt8KwiFUN+m6qqQETUpqo+/Xpv5zj11AgMnsMgWGpF3t6XSqtFoLW3NthHgRCASCVoS7uzs33HADt956Kz169DASMbm5uWzdupW1a9daXMQ4j9D1eylLSKAiN7fO68oSjgG6ZF1rR9XOG1sfHwDKEhqf92IIGan8/KxuLp0hbERlJdrT6Ra1pSEI8SIQCAStEDc3N0aMGFGriLl06ZIsYk6fPm0REeMQFobS1RWAot176rxOTtYNCW4Ru5qLIQTWlH4vGr3nRd21i0ltMgW2bdrIni/t6dOWNaYBCPEiEAgErRiDiJkxY0aNcutLly7xxx9/sGbNGtLS0lpUxChsbXHWT8SuK3RUWVSENk3njbBvBZ4XqNasrgkVRxpDpVFX66k0qo66i05UaU+nWdiSayPEi0AgEFwHuLq6Mnz4cGbMmEHPnj2NRMzly5fZtm0bq1evblER4zRUl/dSsm8fUmVljfNlx46DJIGNDfZBvVrEpuZiqIjSpKRQpdE06l6todLICj0vAOouujwcTZrwvAgEAoGgBXF1dWXYsGHMnDmTnj17olT+8zZ/5coVWcSkpqZSZeZW8E6DBgJQmZ9PWWLNacyl+rwRux49UNrbm9UWU2HfSz9MsbJSbjjXEKpKSijPzASsr9LIgF2A3vOSJjwvABQVFXHq1ClKS0sbfM/Fixc5f/68Ga0SCASC6xcXFxeGDRvGjBkz6NWrVw0Rs337dlavXk1KSorZRIyqXTvU+n4mxfv21ThflpgIWH+JdHVs3NxQ+fkB/9jfEDSnT+u8TICdlYoX2fNioTypxmB28fLMM8/g6enJiBEjaNu2Le+++2691y9btozevXsTFBREREQEHTt2bFVTVQUCgcCacHFxYejQocycOZOgoCAjEZOXl8eff/5pVhHjNHAQAMV799Y4pzmha05n16t1hIwM2AcFAVB2ouHixeDNsPXywkafyGxt2OlzXqSSErkfjbViVvHyzTff8Omnn/L333+TmZnJunXreO6559i8eXOd98TGxrJ69WouXrxIdnY2jz32GDNnziRJ34FRIBAIBI3H2dmZIUOG1CtiVq1axalTp0wqYgyho9KYGKqqed+riovRnjkDgH3P1iZedPbWFgqrC0NfGLWVddatjsrPD4VKBVh/6Mis4uWLL77glltuIUo/4nz06NGMGDGCL774os573nvvPXpVU+Hz588HYM+eukvtBAKBQNAwqouY3r17Y1OtS2x+fj5//fUXq1at4uTJkyYRMY59+4GtLVJ5OSUxsfLxsuSTujCKQoF9YI9m79OSGDwvmuRkpIqKBt0jVxp1sV7xorC1Re3fGdCHuawYs4mXqqoq4uPj6d+/v9HxQYMGERMT0+B1Tp48SXl5OR07djS1iQKBQPCvxdnZmcGDBzNz5kyCg4NriJgdO3awcuXKZosYG2cnHMJ0nWmL9/0TOjKEXNSdO8vTmlsLhjCXpNE0uCeKVl/BY409XqqjlpN2rVu8NGq2UU5ODjk5OfVeExAQgJOTE4WFhWi1Wtq2bWt0vm3btuTW022xOlqtlnvvvZeoqChGjhxZ53UajQZNtZK1goKCBq0vEAgE/3acnJwYNGgQ4eHhxMfHc+LECSr1Zc0FBQXs2LGD2NhYIiIi6N69u1G4qcF7DBxIaUwMxXv/SdrVGIYxtpIS6eqovL2x8fKk8mIuZYmJ8kTmupAqK9GePQuAnb9/C1jYdAxJu9be66VR4uXnn3+uN+QD8O2339K3b19sbXVLa6+akaDRaFDpY2r1UVFRwezZszl37hy7d+82+lRwNW+++SaLFy9uwCsQCAQCQW04OjrKIubIkSMkJiYaiZidO3fKIqZHjx6NEjFOgwaR+8knaE6coOLKFWw9POR8EbtWlu9iwL5XL4ov7qYs8QRukyfXe215VpY8fFLduXNLmNdkDEm71t7rpVHi5ZFHHuGRRx5p0LVOTk54eHjUKHc+f/78NUNAlZWV3HbbbRw8eJAdO3bQqVOneq9/7rnneOKJJ+TvCwoKRJhJIBAImoCjoyMDBw4kLCyMo0ePcvz4cVnEFBYWsmvXLuLi4mRPTH0fLA04hASjcHBAKi2l5PBhXEaMQHNK1yPFvpVVGhmw7xVE8a7dlJ24dtKuNl2XmKxQq7Ft397cpjULQ9ioIjubyqJibJytM6Rn1oTdG264gd9++03+XpIkfvvtN2644Qb5WE5ODieq/eMbhMvevXvZsWMHXbpcOz5oZ2eHq6ur0ZdAIBAImo6joyMDBgxg9uzZhIaGyt50+EfE/Pzzz0ZhprpQqFQ4RkQAUHLoEJq000h6r7zc9K2V8U+59Ilr9kQxVFWpO3dC0YSwW0uirjbtWpuebjlDroFZf4ovvvgi+/btY8GCBezevZv77ruPCxcu8NRTT8nXfPbZZwwcOFD+/q677mLjxo188sknlJSUcOzYMY4dO3bNXBuBQCAQmB4HBwcGDBjArFmzCAsLMxIxRUVF7N69m59//tkozFQbjv36AlBy8JCcrGvr5YWtp6d5X4CZMOTqVBUWUn7uXL3XGkSA2srzXUCXYG3brh1g3XkvZhUvERER/PXXX5w6dYrHH39c/kX3r/YP6O3tTZBewQIcP36czp078/zzzzNz5kz5a/Xq1eY0VSAQCAT14ODgQP/+/Zk9ezbh4eFGuYtFRUXs2bOHFStW1CliHPv1A3TlxSX7DwBg10q9LqDriWKoktKcPFnvtf94Xqw738WAQWQZwl3WSKNyXprCwIED+eWXX+o8/+CDD/Lggw/K3zemjFogEAgELYu9vT39+vUjNDRUzokp1yejFhcXs2fPHuLi4ggPDycwMFD21DgEB6Owt0cqKyN/3TrdWr2C6trG6lEoldh1705pfDyakydxuemmOq81eF5UrUW8dOpEyYEDaDMyLG1KnVh38E0gEAgEVolBxMyaNYuIiAgjT0xxcTF///03K1as4NixY1RUVKBQq3EIDzdawyHC+PvWhl0PXXO9sno8L5JWKw9ktPYyaQPqzroiGW2G9XpehHgRCAQCQZOxt7enb9++zJo1i8jISCMRU1JSwt69e2URY98nSj6ndHTEqVq+Y2vEIF7qmy6tPZcJ+iZ/rcXzotJX+JafEZ4XgUAgEFzH2Nvb06dPH2bPnk1kZCRqtVo+ZxAx7+/dy5TTpwk/mcyU9NOs27TJghY3H7seuuZ02vR0qq7qaWbAEDJSOjpi6+XVUqY1C0NuTmVeHpX5+Ra2pnaEeBEIBAKBybCzs6NPnz7MmjWLqKgoWcTExsbyzvLlnNJq0EoSyVeuMG3aNNauXWthi5uOvd7zQmUl2tTUWq8xJOuq/DujUChayrRmoa7WJ02bcdaCltSNEC8CgUAgMDl2dnZERUUxe/Zs+vTpw6ZNm1AoFBg6okiAQqHgueeekxN+Wxs27u7YensDUJacXOs12jPpQOupNAJjL5G15r0I8SIQCAQCs6FWq4mMjOTixYs1mrlJkkRaWhrLly/nyJEjrVLEXCvvxVBu3Bp6vFRHpU/aLbfSiiMhXgQCgUBgdgIDA2uETRQKBT4+PpSVlXHgwAGWL19OfHx8qxIx/4iX2iuOWluPFwPqTjp7tVaatCvEi0AgEAjMzsKFC5EkSRYwCoUCSZKYOnWqfE1ZWRkHDx5k2bJlxMfH1xjsa40YknZrEy9VZWVU6Of7tT7xoi+XPiPCRgKBQCD4lxIdHc2aNWsIDQ3F3t6e0NBQ1q5dy5IlS+jXrx/29vbytRqNhoMHD7J8+XLi4uKsWsQYknYrcnKozMszOld9bECrEy9yrxfr9LyYvcOuQCAQCASgEzDR0dE1joeHh9O7d28SExM5cuQIZWVlgE7EHDp0iKNHjxISEkJwcLBRCbY1oO7aFWxsoLKSspMncdKPQQDQ6sWL0tERGw8PS5nYJAy9XiovXaKyqAgbZ2cLW2SM8LwIBAKBwOKoVCrCwsKYNWsWAwYMwMHBQT6n0Wg4fPgwy5cvJzY21qo8MUo7O9mrcnXSbvlZnXhRdezYasqkDRjCRmCdSbtCvAgEAoHAalCpVISGhtYrYpYtW0ZMTAwajcaClv6DXbduAGjTjHu9lJ/T9UhRdezQ4jY1FxsXF2zatAGsM3QkwkYCgUAgsDpsbW0JDQ0lKCiIEydOEB8fT2lpKQBarZaYmBgSEhIIDg4mJCQEOzs7i9mq7toFAE1qmtFxrd7zou7QscY9rQF1p06UXr5slRVHwvMiEAgEAqvF1taWkJAQZs2axaBBg3B0dJTPabVaYmNjWbZsGYcPH5ZzZVoau646z4vmas/L2dbreYFqSbtWWHEkPC8CgUAgsHpsbW0JDg6mZ8+eJCUlER8fT0lJCQDl5eXExsYaeWKqVy+ZGzu956XyYi6V+fnYuLkhSZKcsFu93X5rQqX3GFWvmrIWhOdFIBAIBK0Gg4iZOXMmgwcPxsnJST5XXl5OXFwcy5cv5+DBgy3miVEHBIA+IdcQOqq8dAlJH+ZStdKwkaqDzmMkxItAIBAIBCbA1taW3r17M3PmTIYMGVJDxMTHx7eYiFHa28sPek1qCgBafcgIhQKVn69Z9zcX6g5+AJRfuIBUUWFha4wR4kUgEAgErRYbGxuCgoLqFTHLli3jwIEDcsKvObDrogsdafWeF4O3wtbbG6UFk4mbg0GQUVlJeXa2ZY25CiFeBAKBQNDqqS5ihg4dinO1pmoVFRUcOXKE5cuXs3//frOIGHW3rgBoUnVJu9pWnqwLOuGFSgVYX+hIJOwKBAKB4LrBxsaGXr160aNHD06dOkVcXByFhYWATsQcPXqUxMREgoKCCA0NNapeag52XXTiRasXL+WtvEwaQGFjg6p9e8ozMoR4EQgEAoHA3NjY2NCzZ0969OjByZMnaxUxx48fJygoiLCwsGaLGDu956U8K4uq4uJWXyZtQN3Bj/KMDLlyyloQ4kUgEAgE1y1KpVIWMQZPTEFBAQCVlZUkJCSQmJhIr169CA8Pb7KIUetzXgA0p9PRZmbqjrfSMmkDKj9DxVGmhS0xRogXgUAgEFz3KJVKAgMD6d69OykpKcTGxhqJmGPHjnHixAl69epFWFiYUeJvQ7BxccHW25uKnBw0SSeo0Ce4ttYyaQNyuXSmEC8CgUAgEFgEpVJJjx496NatGykpKcTFxZGfnw8Yi5iePXsSHh7eKBFj160rFTk5FO3aDZIEgLqVh41UhnJpETYSCAQCgcCyVBcxqampxMbGGomY48ePG4mY6tVLdaHu0pXivfso2rkTAIW9PTaenmZ9HeZG7acTLxUXL1JVVoayBTsX14cQLwKBQCD416JUKunevTtdu3YlLS2NmJgYWcRUVVWRmJhIUlJSg0SMOsAfAEk/7Vrt749C33m3tSL3ekGXjGxXLbfHkgjxIhAIBIJ/PUqlkm7dutGlSxfS0tKIjY0lLy8PMBYxgYGBRERE1Cpi7AICjL+3kgd9c7Bp2xaFgwNSaSnl585ZzWsS4kUgEAgEAj3VRczp06eJjY3lypUrgE7EnDhxguTkZHr06EFERAQuLi7yveqrxIu6q3U86JuDQj/eQJuSalXl0kK8CAQCgUBwFUqlkq5du8oiJiYmxkjEJCUlkZycTGBgIOHh4bi6umLbrh0KR0ck/bRru67dLPkSTIbarwPalFSrKpcW4kUgEAgEgjpQKBR06dKFgIAA2RNz+fJlACRJkkWMwRNj4+pKhSxeWr/nBayzXFqIF4FAIBAIrkF1EZOenk5sbCyXLl0CdCImOTmZkydPUpyZyecZGaSXa+k5bRoLFy0iOjrawtY3D5Wf9ZVLi8GMAoFAIBA0EIVCQUBAANHR0YwePZq2bdvK52JiYpifcopTWg1aSSLh2DGmTZvG2rVrLWhx87HGXi9CvAgEAoFA0EgUCgX+/v5GImbjxo0oAEl/jSRJKBQKFi5caElTm43KVydeKvPzqdKHxCyN2cWLJEns3buX5cuXc/jw4Ubde+nSJZYuXcqOHTvMY5xAIBAIBM2guojJzc2VhYsBQ17Mn3/+KZdetzZUvu3l/y/Xjz2wNGYVL2VlZYwePZpp06axYsUKxowZw4wZM6isrGzQ/XPnzuWBBx7gk08+MaeZAoFAIBA0C4VCQWBgYI2mdAqFAh8fH1JSUli1alWrFDE2Hh4o7OwAKM86b2FrdJhVvLz33nscO3aM+Ph4fv31V/bv38/GjRtZunTpNe/98MMP0Wg03HjjjeY0USAQCAQCk7Bw4UI5VAQ64SJJEhMnTgR0XpiUlBRWrlzJ9u3b5dJra0ehUKBqr/O+lGdZR8WRWcXLsmXLmDFjBu3atQOge/fujB8/nmXLltV7X2xsLO+++y7fffddq2+tLBAIBIJ/B9HR0axZs4bQ0FDs7e0JDQ1lzZo1PP/883h7extdm5qayqpVq9i2bZtcem3NGEJH5eetw/NitlLp8vJykpOTmT9/vtHx4OBgPvvsszrvKyoqYubMmXzyySe0b9++zuuqo9Fo0OhnSQDymHOBQCAQCFqS6OjoWkujO3ToQGZmJjExMVy4cEE+npaWRlpaGl26dCEyMpI2bdq0pLkNxlb/PK6wkrBRo8RLXFwcR44cqfeasWPH4uPjQ1FREVVVVbi7uxudb9OmjTz0qjYeeOABhg8fztSpUxts15tvvsnixYsbfL1AIBAIBC2JQqGgQ4cO+Pn51StiAgICiIqKsjoRo2rvC7RSz8vp06evWfkzYMAAfHx8cHBwAKC4uNjofGFhoXzuanbu3MnKlSt5//335byYc+fOoVKpWLp0Kbfcckutw7Cee+45nnjiCfn7goICOnbs2IhXJhAIBAKB+akuYrKysoiJiSG7WgXP6dOnOX36NP7+/kRFRRn1kbEkKt9WLF7qcofVhr29Pb6+vqSnpxsdT09Pp2vXrrXe4+LiwqxZszh06JB87NKlS9jY2LBjxw4mTpxY+yRPOzvs9JnQAoFAIBBYOwqFAj8/P3x9fTl//jwxMTGcryYM0tPTSU9Px9/fn8jISDw9PS1obbWcl+xspKoqFErLtolTSJJ0dVm6ybj//vvZu3cvcXFx2NraUlxcTLdu3bj//vtZtGgRAPHx8SQmJjJ79uxa15g4cSL29vasXr26wfsWFBTg5uZGfn4+rq6upngpAoFAIBCYFYMn5nwt3o3OnTsTFRVlMRGjPXOG1DFjAei2ayeqqxKQTUVDn99mlU4vvfQSly5dYsKECXzwwQeMGjUKFxcXHn/8cfmadevW8eCDD5rTDIFAIBAIrB5fX19uvvlmJk6ciK8+TGPgzJkzrF27li1btnDx4sUWt83Wx0f+/4qsrBbf/2rMOpixQ4cOxMXF8dVXX5GUlER0dDT33Xcfbm5u8jXh4eF1el0ARo4ciUqlMqeZAoFAIBBYDb6+vnI4KTY2lsxq05wzMjLIyMigU6dOREZG1ijBNhdKOztsPD2pzM2l/Px5HMLDW2TfujBr2MhSiLCRQCAQCK4XsrOziYmJMRIxBjp27EhUVFSLiJjT02+lLCEB76efpu09d5tlj4Y+v83qeREIBAKBQNA8fHx8mDBhAtnZ2cTGxnKu2nTns2fPcvbsWTp27EhkZKTcFNYcqNq3pywhwSoqjoR4EQgEAoGgFeDj48P48eO5cOECsbGxnD17Vj5nEDEdOnQgMjISn2o5KqZCHhEgxItAIBAIBILG0K5dO8aNG0dOTg4xMTFGIubcuXOcO3cOPz8/oqKiTCpiVH6GXi/XecKuQCAQCAQC8+Dt7S2LmNjYWDIyMuRzmZmZZGZm4ufnR2RkZIPH7dSHNY0IEOJFIBAIBIJWjLe3N2PHjuXixYvExMTUKmJ8fX2JjIysUYLdGAwjAirz8qgqKUHp6Nhs25uKEC8CgUAgEFwHeHl5MXbsWHJzc4mJieHMmTPyuaysLLKysmjfvj1RUVFNEjGGLrugy3uxq6Nbfktg2f6+AoFAIBAITIqnpydjxowhOjoaf39/o3Pnz59n48aNbNiwgczMTBrTLcXGwwOFWg3oxgRYEiFeBAKBQCC4DvH09GT06NFMmzatVhGzadOmRokYhUIhd9qtyL5wjavNixAvAoFAIBBcx7Rt21YWMQEBAUbnsrOz2bRpE+vXr+fcuXPXFDEqfR+Z8mzLJu2KnBeBQCAQCP4FtG3bllGjRnH58mViY2NJS0uTz124cIHffvuNdu3aERkZSYcOHVAoFDXWsG1vHZ4XIV4EAoFAIPgX0aZNG0aOHFmniNm8eTPe3t5ERUXVEDGqdjrxUn7BsjkvQrwIBAKBQPAvxCBirly5QmxsLKmpqfK5nJwcNm/ejJeXF1FRUXTs2FGX82LwvJwX4kUgEAgEAoGF8PDw4KabbiIyMpK4uDhSU1Pl3JeLFy+yZcsWvLy8iIyMxMOQ83LBsmEjMVVaIBAIBAKBTF5enuyJuVoinDl4kNU//Eh6uZaevXuzcPFioqOjTbZ3Q5/fQrwIBAKBQCCoQV5eHnFxcaSkpCBJErGxsXzxxRcoAAld6bQkSaxZs8ZkAkaIFyFeBAKBQCBoNvn5+cTFxXH77bfX6AmjUCgIDQ0lPj7eJHs19Pkt+rwIBAKBQCCoEzc3N0aMGMHFixdrhJEkSSI5ObnFbRLiRSAQCAQCwTUJDAys0ftFoVAQGBjY4rYI8SIQCAQCgeCaLFy4EEmSZAFjyHlZuHBhi9sixItAIBAIBIJrEh0dzZo1awgNDcXe3p7Q0FDW/n979x9TVf3/AfyJF0EQUFFUfgl6NwhENgi6oCE/YmEwhy5B08IcaaamTFCbk1CrlW7+2Ny0Yk5j64fTaeUiK6mUwpkwoWmmiwi4k18Cdm8gcIHX9w/z7HM/cBE/X+7FA8/H5h/ndc5LXr7u9d4X5573PadPY/HixTavhRfsEhER0WOBF+wSERHRiMThhYiIiFSFwwsRERGpCocXIiIiUhUOL0RERKQqI/Ku0g8WUBkMhmGuhIiIiAbrwfv2wxZCj8jhxWg0AgB8fX2HuRIiIiJ6VEajERMmTLC4f0R+z0tvby9u374NV1fXPl9l/P9hMBjg6+uL2tpafn+MFbHPtsNe2wb7bBvss21Ys88iAqPRCC8vL4wZY/nKlhF55mXMmDHw8fGx2t/v5ubG/xg2wD7bDnttG+yzbbDPtmGtPg90xuUBXrBLREREqsLhhYiIiFSFw8sjcHR0RF5eHhwdHYe7lBGNfbYd9to22GfbYJ9t43Ho84i8YJeIiIhGLp55ISIiIlXh8EJERESqwuGFiIiIVIXDy3+4c+cOsrOzERcXhyVLluDbb7+1Ss5o19nZiffeew+JiYlISUnB8ePHH5pz8eJFrFmzBs888wxWrlyJ4uJi6xeqciKCo0ePIjk5GYmJidi7dy+6uroGnZ+VlYWoqCj8/PPPVqxyZDh37hyef/55xMfHY8uWLWhpaXlozt27d7F79248++yzWLp0KZ/Tg1BeXo6MjAzExsbilVdewR9//PHQnBMnTiA9PR1xcXFYvnw5vvrqKxtUqm41NTXYsWMH5s6di8OHDw8qR6/XY926dYiLi8MLL7yAkpISq9bI4eVfnZ2diI2NRVlZGXJychAWFobk5GR8+eWXQ5pDwPLly3H06FGsXbsWixcvxsaNG/H2229bPP6DDz5Abm4uIiMjsX37dmi1WsTHx+PTTz+1YdXqs3PnTmRnZ2PJkiV49dVXceTIEWRkZAwq99ixYygqKsLly5fR2tpq5UrV7fTp01i4cCEiIyOxefNmXLp0CXFxcQMOio2NjYiMjMSPP/6I119/HRkZGdi1a9eg3oxHq2vXruHpp5+Gi4sL3njjDbS1tSEqKgp6vd5izv79+7Fq1SrExsZi586dCA0NRWpqKj755BMbVq4uFy5cQGxsLBwcHFBXV4eampqH5rS2tiI6Ohp6vR5bt26Fr68v4uLirDvACImIyIcffiiOjo7S2tqqxDIzMyUkJGRIc0a7K1euCAC5dOmSEjt06JA4OTmJwWDoN8doNPaJrVq1SiIjI61Wp9q1traKo6Oj5OfnK7EffvhBAEh5efmAub///rt4enoqj9XZs2etXa6qBQYGyvr165XtxsZGGTt2rBw/ftxizsqVKyUwMFA6OjqUWG9vr9k2mUtPT5d58+Yp293d3aLVaiUrK8tiTkxMjLz88stmsaSkJElLS7NanWpnNBqlu7tbRERmz54t27Zte2jO7t27xcPDQzo7O5VYSkqKJCYmWq1Onnn5V1FREWJiYjBx4kQllpqaimvXrqGhoWHIcka78+fPY8qUKYiKilJiqampuHfvnsUp3cXFpd/Yo3wEMtr89NNP6OzsxMKFC5XY/PnzMXHiRJw/f95iXmdnJ5YuXYq9e/fC39/fBpWqm16vx82bN8367OHhgejoaIt97urqwokTJ5CZmWn2PRl2dnb8fpIBFBUVmfVZo9EgJSVlwOdzREQEysvL0d7eDgBoaWnBjRs38NRTT1m9XrVycXGBRqN5pJyioiIkJSXBwcFBiaWmpuLChQswmUxDXSIAfmykqK6uhpeXl1nswXZ1dfWQ5Yx2/fXM09MTdnZ2g+5ZbW0tCgoKsGjRIitUODJUV1dDo9Fg6tSpSmzMmDGYPn36gH3Ozs5GUFAQXnzxRVuUqXoPetnf64ClPldVVaGjowOzZs3Chg0bEB8fj4yMDF5bNIC2tjY0Nzc/Up8BYM+ePYiJiYGPjw/CwsIwc+ZMZGZmIjs729oljyqW3gtNJhPq6uqs8jM5vPzLZDL1+a3HyclJ2TdUOaNdfz2zt7eHvb39oHr2999/IzU1FcHBwdi+fbu1ylQ9k8kEBweHPndVd3Jystjnzz//HGfPnsWRI0dsUeKI8KCX/b0OWOpzR0cHAGD9+vXw9/dHbm4uZsyYgZiYGHz99dfWLVil/pc+A8Bnn32GgoIC5OXlYd++fdi2bRv27dvHhRVDbDjeC0fkXaX/F+7u7n1WCDQ3NwMAJk+ePGQ5o11/PTMajTCZTA/tmdFoxIIFC6DRaFBYWGh2ipLMubu74969e+jo6MC4ceOUeHNzs8U+f/TRR+jp6cGCBQsAAN3d3QCAnJwcFBYWDnrVwWji7u4OAP2+Dgz0ugEAaWlpyMnJAQAkJCSgoqIChw4dwnPPPWfFitXJ1dUVY8eOfaQ+A8DmzZuRlZWFTZs2Abjf55qaGmzZsgVJSUlWrXk0Gei98MHzfajxzMu/wsPDUVpaaha7fPkyXF1dodVqhyxntAsPD8dff/2FO3fuKLHLly8DAMLCwizmPRhcTCYTvvvuO7PrjKiv8PBwAMCVK1eUWF1dHfR6vcU+79mzB6dOncLBgwdx8OBBvPPOOwCA1atXY8OGDdYvWoUCAwMxfvx4sz6LCEpLSy322dfXF9OmTev341Ou7OqfRqNBaGioWZ+B+68dlvrc09MDg8EAb29vs7iXl9eglrLT4IWHh/f72Pj7+2PSpEnW+aFWuxRYZX777TfRaDRy7NgxEbm/YmDmzJlmqwgqKipEp9PJjRs3Bp1D5oxGo0ydOlU2bdokIiJdXV2SkJBgtoqgvb1ddDqdfPHFF0rOvHnz5Mknn5SWlpbhKFuVdDqdJCUliclkEhGRtWvXiqenp7S1tSnHLFq0SA4cONBvflNTE1cbDcKaNWskICBAmpubRUTk/fffF3t7e7l165ZyzNatW+W1115TtnNzcyUkJER5PtfU1Mi0adNkx44dti1eRQ4fPixubm5y/fp1EREpLi4We3t7OXPmjHJMfn6+2QqXuLg4iY6OVlYyNjU1SWBgoCxbtsymtauVpdVG33//veh0OmlsbBQRkYsXL4qdnZ0UFhaKiEhVVZV4eHjIW2+9ZbXaOLz8h4KCAnFxcRGtVitOTk6SnJws//zzj7K/uLhYAMiVK1cGnUN9FRcXi5eXl3h7e8ukSZMkNDRUqqqqlP1Go1EAKMt8d+3aJQAkKChIdDqd8ic+Pn6Y/gXqUFlZKSEhIeLu7i5eXl7i4+MjJSUlZsf4+fkpg+R/4/AyOAaDQZKSksTZ2Vm0Wq24urrKxx9/bHZMamqqxMbGKtudnZ2yYsUKcXNzkzlz5oizs7O89NJLXCo9gN7eXlm3bp04ODhIQECAODo6Sm5urtkxeXl5MmHCBGW7srJSoqKixM3NTUJDQ2X8+PGSmJgoDQ0NNq5ePZqbm5XXWCcnJ/Hy8hKdTieZmZnKMSdPnhQAUltbq8T2798vTk5OymOzYsUK6erqslqdvKv0f2lvb8fNmzcxefJkzJgxw2yf0WjE9evXMWfOHIwfP35QOdS/7u5u3LhxA46OjggICDDb19vbi19++QVarRYeHh7Q6/X9fhGVRqNBZGSkrUpWrZs3b6KrqwtBQUGwtze/zO3q1atwd3eHn59fn7zu7m6UlpbiiSee4Md0g1BTU4OWlhYEBATA2dnZbN+tW7fQ3d2N4OBgs3hDQwPq6+vh5+fHHg9SU1MT9Hp9vx9J6PV61NfXIyIiwixeX1+PhoYGeHt7Y8qUKbYsV3VMJhPKysr6xF1dXTF79mwA96/xunXrFsLDw82uPTQYDKisrMT06dPh6elp1To5vBAREZGq8IJdIiIiUhUOL0RERKQqHF6IiIhIVTi8EBERkapweCEiIiJV4fBCREREqsLhhYiIiFSFwwsRERGpCocXIiIiUhUOL0T02Pvzzz8RHx+P8vJyJVZRUYGEhATlruRENHpweCGix96sWbMwbtw4vPnmmwDu30coOTkZ8+fPh06nG+bqiMjWeG8jIlKFq1evIiIiAufOnUNWVhbmzp2L/Pz84S6LiIYBhxciUo1ly5bh1KlTSE5OxpkzZ6DRaIa7JCIaBvzYiIhUoaOjA7dv30ZPTw9Wr17NwYVoFOOZFyJ67PX29iItLQ16vR7BwcGoqKhAWVkZ7Ozshrs0IhoG9sNdABHRw2zcuBG//vorSkpK0NXVBa1Wi5MnTyI9PX24SyOiYcAzL0T0WHv33Xdx4MABXLp0CVqtFsD9Yeabb77B9evXYW/P38GIRhsOL0T02Ort7UVxcTH8/Pzg7++vxO/evYvy8nJERETAxcVl+AokomHB4YWIiIhUhauNiIiISFU4vBAREZGqcHghIiIiVeHwQkRERKrC4YWIiIhUhcMLERERqQqHFyIiIlIVDi9ERESkKhxeiIiISFU4vBAREZGqcHghIiIiVeHwQkRERKryf6QMDZIrDUriAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x320 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "midpoint / exact ratio per element: [-0.25 -0.25 -0.25 -0.25 -0.25 -0.25 -0.25 -0.25]\n"
     ]
    }
   ],
   "source": [
    "# Cell 10 [LECTURE] -- P2-P0: the flux oscillates (n = 8, f = 1)\n",
    "xs = np.linspace(0, 1, 400)\n",
    "nodes = np.linspace(0, 1, 9)\n",
    "mids = (np.arange(8) + 0.5) / 8\n",
    "plt.figure(figsize=(6.4, 3.2))\n",
    "plt.plot(xs, 0.5 - xs, lw=2.2, color=\"0.6\",\n",
    "         label=\"exact $\\\\sigma = 1/2 - x$ (belongs to $V_h$)\")\n",
    "plt.plot(xs, eval_at(sig_h, xs), lw=1.6, color=\"tab:red\",\n",
    "         label=\"$\\\\sigma_h$  ($P_2$-$P_0$, $n=8$)\")\n",
    "plt.plot(nodes, eval_at(sig_h, nodes), \"o\", color=\"k\", ms=4,\n",
    "         label=\"vertex values: exact\")\n",
    "plt.plot(mids, eval_at(sig_h, mids), \"s\", color=\"tab:red\", ms=4,\n",
    "         label=\"midpoint values: $-\\\\sigma/4$\")\n",
    "plt.legend(frameon=False, fontsize=9); plt.xlabel(\"$x$\"); plt.show()\n",
    "\n",
    "ratio = eval_at(sig_h, mids) / (0.5 - mids)\n",
    "print(\"midpoint / exact ratio per element:\", np.round(ratio, 4))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-11",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 11 [LECTURE] -- the errors converge to nonzero limits\n",
    "case = \"f=pi^2 sin(pi x)\"\n",
    "print(f\"{'n':>5} {'P1-P0: rel sig':>15} {'P2-P0: rel sig':>15} \"\n",
    "      f\"{'P2-P0: rel u':>13}\")\n",
    "# Each solve_pair call constructs its own mesh object. A form must not mix\n",
    "# Functions defined on two different mesh objects, even when the meshes\n",
    "# coincide geometrically (UFL: \"Multiple domains found\"). Each error is\n",
    "# therefore assembled on the mesh of its own solve.\n",
    "for n in (8, 16, 32, 64, 128):\n",
    "    mesh1, w10 = solve_pair(n, (\"CG\", 1), (\"DG\", 0), case)\n",
    "    x1, = SpatialCoordinate(mesh1)\n",
    "    e10 = rel_L2(w10.subfunctions[0] - CASES[case][\"sig\"](x1),\n",
    "                 CASES[case][\"sig\"](x1))\n",
    "    mesh2, w20 = solve_pair(n, (\"CG\", 2), (\"DG\", 0), case)\n",
    "    x2, = SpatialCoordinate(mesh2)\n",
    "    e20 = rel_L2(w20.subfunctions[0] - CASES[case][\"sig\"](x2),\n",
    "                 CASES[case][\"sig\"](x2))\n",
    "    eu20 = rel_L2(w20.subfunctions[1] - CASES[case][\"u\"](x2),\n",
    "                  CASES[case][\"u\"](x2))\n",
    "    print(f\"{n:>5} {e10:15.3e} {e20:15.4f} {eu20:13.4f}\")\n",
    "print(\"\\nlimits: sqrt(5/6) =\", float(np.sqrt(5 / 6)), \"  5/6 =\", 5 / 6)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-12",
   "metadata": {},
   "source": [
    "## Self-study cells\n",
    "\n",
    "The following cells contain the material omitted from the lecture; they are\n",
    "required by Exercise Sheet 1. The add-back numbers refer to the timing plan of\n",
    "the course materials.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-13",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 13 [SELF-STUDY S1 / ADD-BACK 5] -- solvability and the parity of n\n",
    "# The singular P1-P1 system is solvable iff (f, chi_h) = 0, with chi_h the\n",
    "# alternating nodal mode. The value of (f, chi_h) is computed here.\n",
    "def chi_pairing(n, case):\n",
    "    mesh = UnitIntervalMesh(n)\n",
    "    Q = FunctionSpace(mesh, \"CG\", 1)\n",
    "    x, = SpatialCoordinate(mesh)\n",
    "    xd = Function(Q).interpolate(x)\n",
    "    chi = Function(Q)\n",
    "    chi.dat.data[:] = (-1.0) ** np.rint(xd.dat.data * n)\n",
    "    return assemble(CASES[case][\"f\"](x) * chi * dx)\n",
    "\n",
    "cases3 = [\"f=1\", \"f=pi^2 sin(pi x)\"]\n",
    "print(f\"{'n':>4} \" + \"\".join(f\"{c:>22}\" for c in cases3))\n",
    "for n in (4, 5, 8, 64, 65):\n",
    "    vals = [chi_pairing(n, c) for c in cases3]\n",
    "    print(f\"{n:>4} \" + \"\".join(f\"{v:22.3e}\" for v in vals))\n",
    "# Compatibility depends on the parity of n (Exercise 1(c)); note that f = 1\n",
    "# is compatible for every n, while f = pi^2 sin(pi x) is compatible only for\n",
    "# odd n.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 14 [SELF-STUDY S2 / ADD-BACK 1] -- regularized solves, two data\n",
    "# Zero block replaced by -1e-8 * M_Q; identical code and mesh, two data.\n",
    "def eps_solve(n, case):\n",
    "    mesh = UnitIntervalMesh(n)\n",
    "    V = FunctionSpace(mesh, \"CG\", 1)\n",
    "    Q = FunctionSpace(mesh, \"CG\", 1)\n",
    "    W = V * Q\n",
    "    sigma, u = TrialFunctions(W)\n",
    "    tau, v = TestFunctions(W)\n",
    "    a_eps = (sigma * tau + tau.dx(0) * u + sigma.dx(0) * v\n",
    "             - Constant(1e-8) * u * v) * dx\n",
    "    x, = SpatialCoordinate(mesh)\n",
    "    L = -CASES[case][\"f\"](x) * v * dx\n",
    "    w = Function(W)\n",
    "    solve(a_eps == L, w, solver_parameters={\n",
    "        \"mat_type\": \"aij\", \"ksp_type\": \"preonly\",\n",
    "        \"pc_type\": \"lu\", \"pc_factor_mat_solver_type\": \"mumps\"})\n",
    "    return mesh, w.subfunctions[1]\n",
    "\n",
    "n = 64\n",
    "nodes = np.linspace(0, 1, n + 1)\n",
    "fig, ax = plt.subplots(1, 2, figsize=(9, 3))\n",
    "for a_, case, ttl in ((ax[0], \"f=1\", \"f = 1: apparently correct\"),\n",
    "                      (ax[1], \"f=pi^2 sin(pi x)\", \"f = pi^2 sin(pi x)\")):\n",
    "    _, u_h = eps_solve(n, case)\n",
    "    vals = eval_at(u_h, nodes)\n",
    "    a_.plot(nodes, vals, lw=1.0)\n",
    "    a_.set_title(ttl + f\"   max|u_h| = {np.abs(vals).max():.3g}\", fontsize=10)\n",
    "    a_.set_xlabel(\"$x$\")\n",
    "plt.show()\n",
    "# At n = 65 the two cases exchange roles (parity; Exercise 1(d)).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-15",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 15 [SELF-STUDY S3 / ADD-BACK 2] -- P2-P0: the scalar converges to u/6\n",
    "case = \"f=pi^2 sin(pi x)\"\n",
    "n = 32\n",
    "mesh, w = solve_pair(n, (\"CG\", 2), (\"DG\", 0), case)\n",
    "u_h = w.subfunctions[1]\n",
    "x, = SpatialCoordinate(mesh)\n",
    "xs = np.linspace(0, 1, 400)\n",
    "mids = (np.arange(n) + 0.5) / n\n",
    "plt.figure(figsize=(6.2, 3.1))\n",
    "plt.plot(xs, CASES[case][\"u_np\"](xs), lw=2.2, color=\"0.6\", label=\"exact $u$\")\n",
    "plt.plot(xs, CASES[case][\"u_np\"](xs) / 6, \"--\", color=\"k\", lw=1.4,\n",
    "         label=\"$u/6$\")\n",
    "plt.step(mids, eval_at(u_h, mids), where=\"mid\",\n",
    "         color=\"tab:red\", lw=1.4, label=\"$u_h$  ($n=32$)\")\n",
    "plt.legend(frameon=False); plt.xlabel(\"$x$\"); plt.show()\n",
    "\n",
    "rel = rel_L2(u_h - CASES[case][\"u\"](x) / 6, CASES[case][\"u\"](x))\n",
    "print(\"|| u_h - u/6 || / || u || =\", float(rel))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Cell 16 [SELF-STUDY S4 / ADD-BACK 4] -- the kernel-coercivity ratio alpha_h\n",
    "# alpha_h = min over K_h of ||tau||^2 / (||tau||^2 + ||tau'||^2), computed as a\n",
    "# generalized eigenvalue problem on an explicit basis of K_h\n",
    "# (the global constant and the elementwise quadratic bubbles).\n",
    "def alpha_h(n):\n",
    "    mesh = UnitIntervalMesh(n)\n",
    "    V = FunctionSpace(mesh, \"CG\", 2)\n",
    "    uu, vv = TrialFunction(V), TestFunction(V)\n",
    "\n",
    "    def dense(form):\n",
    "        Am = assemble(form, mat_type=\"aij\")\n",
    "        ai, aj, av = Am.petscmat.getValuesCSR()\n",
    "        return sp.csr_matrix((av, aj, ai)).toarray()\n",
    "\n",
    "    Mm = dense(uu * vv * dx)\n",
    "    Kk = dense(uu.dx(0) * vv.dx(0) * dx)\n",
    "\n",
    "    x, = SpatialCoordinate(mesh)\n",
    "    xd = Function(V).interpolate(x).dat.data\n",
    "    mid = np.isclose((xd * n) % 1.0, 0.5, atol=1e-8)  # midpoint dofs = bubbles\n",
    "    Z = np.zeros((V.dim(), 1 + int(mid.sum())))\n",
    "    Z[:, 0] = 1.0\n",
    "    Z[np.where(mid)[0], 1 + np.arange(int(mid.sum()))] = 1.0\n",
    "    lam = sla.eigh(Z.T @ Mm @ Z, Z.T @ (Mm + Kk) @ Z, eigvals_only=True)\n",
    "    return lam.min()\n",
    "\n",
    "print(f\"{'n':>5} {'alpha_h':>12} {'alpha_h / h^2':>14}\")\n",
    "for n in (8, 16, 32, 64):\n",
    "    a_ = alpha_h(n)\n",
    "    print(f\"{n:>5} {a_:12.4e} {a_ * n**2:14.4f}\")\n",
    "print(\"\\nreference: 1/60 =\", 1 / 60)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-17",
   "metadata": {},
   "source": [
    "## Pointers\n",
    "\n",
    "Exercise Sheet 1 relies on cells `S1`–`S4`. The two-dimensional experiment of\n",
    "the lecture is contained in the companion notebook `L1B_q1p0_cavity.ipynb`.\n",
    "The theory answering the questions raised by these experiments is the subject\n",
    "of Lecture 2.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Firedrake",
   "language": "python",
   "name": "firedrake"
  },
  "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}