maedoc · June 9, 2022 11:57 · maedoc · Jun 9, 2022
diff --git a/Implicit methods by hand.ipynb b/Implicit methods by hand.ipynb
 {
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "# Implicit methods by hand\n",
    "\n",
    "We usually use explicit methods for solving ODEs, DDEs, SDEs, etc because they are easy to code. Unfortunately, the tradeoff is poor stability. Implicit methods handle stability well, but we need to adapt them to our problems, sometimes autodiff'ing all the way through the method.\n",
    "\n",
    "The trusty Fortran codes aren't built for this, so we need to build up our algorithms by hand. A few elements are required:\n",
    "\n",
    "0. the model itself\n",
    "1. model jacobian\n",
    "2. implicit scheme such as backwards Euler or theta method\n",
    "3. linear system solver, such as Jacobi iteration\n",
    "\n",
    "In our case, we want each part to propagate derivative information usefully, so we will go directly for an iterative solver because direct solvers are more challenging to differentiate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## Linear system solver\n",
    "\n",
    "Here, we'll set up a weighted variant of the iterative Jacobi scheme for solving $A x = b$\n",
    "\n",
    "$ x_{i+1} = w D^{-1} (b - LU x_i) + (1 - w) x_i $\n",
    "\n",
    "where $D$ is the diagonal of $A$, and $LU = A - D$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def jacobi_np(A, b, w=2.0/3.0, tol=1e-9, max_iters=100):\n",
    "    w_m_invD = w * diag(1.0 / diag(A))\n",
    "    LU = A - diag(diag(A))\n",
    "    x = b.copy()\n",
    "    dx = ones(x.shape)\n",
    "    n_iter = 0\n",
    "    while n_iter < max_iters and linalg.norm(dx) > tol:\n",
    "        xn = w_m_invD.dot(b - LU.dot(x)) + (1.0 - w)*x\n",
    "        dx = x - xn\n",
    "        x = xn\n",
    "        n_iter += 1\n",
    "    return x, n_iter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "We set up a test linear system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[10., -1., 2., 0.],\n",
    "              [-1., 11., -1., 3.],\n",
    "              [2., -1., 10., -1.],\n",
    "              [0.0, 3., -1., 8.]])\n",
    "\n",
    "b = np.array([6., 25., -11., 15.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Finally, we compute the solution via the Jacobi method and compare with the result of a direct solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 1.,  2., -1.,  1.]), 40)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_j = jacobi_np(A, b)\n",
    "x_j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.,  2., -1.,  1.])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_d = linalg.solve(A, b)\n",
    "x_d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## Implicit scheme\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def _theta_dy(f, j, h, th, y0, y1, f_y0, pars, tol=1e-4):\n",
    "    return jacobi_tt(\n",
    "        np.identity(y1.size) - h*th*j(y1, *pars), # J\n",
    "        y0 + h * (th * f(y1, *pars) + f_y0) - y1, # -F\n",
    "        2.0/3.0, tol, 10\n",
    "    )\n",
    "\n",
    "def theta(f, j, y0, h, tf, *pars, th=0.5, sigma=0.0, tol=1e-4):\n",
    "    ys = [y0.copy()]\n",
    "    n = int(tf / h)\n",
    "    add_noise = sigma or any(sigma > 0.0)\n",
    "    for i in range(n):\n",
    "        y0 = ys[-1]\n",
    "        # compute standard Euler as initial guess\n",
    "        f_y0 = (1 - th) * f(y0, *pars)\n",
    "        y1 = y0 + h * f_y0\n",
    "        # refine following theta method for theta > 0.0\n",
    "        if th > 0.0:\n",
    "            dy = _theta_dy(f, j, h, th, y0, y0, f_y0, pars)\n",
    "            y1 = y0 + dy\n",
    "            n_iter = 1\n",
    "            while np.linalg.norm(dy) > tol:\n",
    "                dy = _theta_dy(f, j, h, th, y0, y1, f_y0, pars, tol=tol)\n",
    "                y1 += dy\n",
    "                n_iter += 1\n",
    "        # first-order Wiener increment\n",
    "        if add_noise:\n",
    "            y1 += sqrt(sigma) * np.random.randn(y0.size)\n",
    "        ys.append(y1)\n",
    "    return np.r_[:n+1]*h, np.array(ys)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def osc(xy, tau):\n",
    "    x, y = xy\n",
    "    return np.array([tau*(x - x**3/3 + y), (1.1 - x)/tau])\n",
    "\n",
    "def osc_jac(xy, tau):\n",
    "    x, y = xy\n",
    "    return np.array([[tau*(1 - x**2), tau], [-1/tau, 0]])\n",
    "\n",
    "\n",
    "x0 = np.r_[1.0, 0.0]\n",
    "tf = 100.0\n",
    "dt = 0.01\n",
    "\n",
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "t, y = theta(osc, osc_jac, x0, dt, tf, 10.0, th=0.5, sigma=1e-3, tol=1e-9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x11c3ba7f0>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEACAYAAABBDJb9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXlcE3f+/18fQgCPAS88QEFBQELwQCuKYrXaerfb1m3d\nHtt27b22Xbvbb/uz3a11t3drW7u1q629t9fay7XarrZYUCyteBEQQRBQ4gUejCJX8vn9MSSZM5lJ\nQog4z8cjDzKf+cxnPpmE93zmfRJKKXR0dHR0Lm5COnsCOjo6Ojq+owtzHR0dnS6ALsx1dHR0ugC6\nMNfR0dHpAujCXEdHR6cLoAtzHR0dnS6Az8KcEBJOCCkghOwmhBQRQp70x8R0dDoKQshgQsiPhJDi\n9t/sgwr9VhJCygkhewghowM9Tx0dLYT6OgCltJkQMo1S2kgIMQDYTgjZRCn9xQ/z09HpCNoAPEwp\n3UMI6QmgkBDyP0ppqaMDIWQ2gERKaRIhJBPAvwBM6KT56uh4xC9qFkppY/vbcHA3CD0SSSdooZQe\no5TuaX9/DsB+ALGibtcA+KC9TwGAKELIgIBOVEdHA34R5oSQEELIbgDHAGymlP7qj3F1dDoaQshQ\nAKMBFIh2xQI4zNuuhVTg6+gEDf5amdsppWMADAaQSQgx+WNcHZ2OpF3Fsg7AQ+0rdB2dixafdeZ8\nKKUNhJAcALMAlPD3EUJ01YtOh0IpJWr7EkJCwQnyDyml38h0qQUwhLc9uL1Nbiz9t63Toaj5bfvD\nm6UfISSq/X03AFcCKJXrSykN+OvJJ5/slPN25rkvxc/sBe8AKKGUvqawfz2A37f/ricAOEMpPa40\nWFe4roH67rrKZwnU9VKLP1bmgwC8TwgJAXdz+IxSutEP4+rodAiEkEkAbgZQ1G7roQCWAogHQCml\nayilGwkhcwghBwGcB3BH581YR8cz/nBNLAKQ4Ye56OgEBErpdgAGFf0WB2A6Ojp+octHgE6dOvWS\nO/el+Jm7OoG4roH67rrKZwm23zrRopPx6USE0ECdS+fSgxACqsEA6udz679tnQ5D7W+7y6/MdXR0\ndC4FdGGuo6Oj0wXQhbmOjo5OF0AX5jo6OjpdAF2Y6+jo6HQBdGGuo6Oj0wXQhbmOjo5OF0AX5jo6\nOjpdAF2Y6+jo6HQB/JoCV0dHR+dSh2VZFBRwtU4yMzPBMExAzqsLcx0dHR0/wbIsUlNTUVvLpb5P\nSEhAXl4eYmJiOvzcuppFR0dHx0/k5OQ4BTkAVFZW4vLLLwfLsh1+bl2Y6+jo6PiJqqoqSVtlZSWK\ni4s7/Ny6MNfR0dHxE5MnT5a0hYaGIi4ursPPrQtzHR0dHT+xc+dOSZvdbkdNTU2Hn1sX5jo6Ojp+\nwmg0StoMBoO+Mu9qsCywYwf3VwewWoH77gOSk4Gvvurs2ejo+IbVasWiRYsk7W1tbQFZmeuuiR0M\nywI5OUBpKbBqFVBdDSQkAHv2AAFyPw1KrFYgNta1fd11wJdfAtde67pmVVXAggVAALy6dHR8Zt26\ndZCrOJWYmIi0tLQOP78uzDsQlgUuuww4cEDYXlkJ/POfQI8ewNChwLRpl5Zgt1qBO++Utv/xj8CM\nGcDo0dw1AoBHHgEOHdIFuk7w079/f9n2tra2gJz/olOzOFQVZWXAmjWcYAhWLBapIHewdCnw0EPA\nNdcAo0ZdOqoXqxUYNgzYtEm678wZoKDAJcgBoKUF+Phj17auqtK52KipqQmIa+JFtTK3WoHx4wGe\nTz5CQznVxblzwKuvAoMGAYsWBcdKrm9fdf0OHeJUDLfd1rHzUcJqBTZsAObN6/jr9vzznICWo6UF\n+PxzafuzzwJhYcCsWdwcDx4EkpKAnTsvrScaneDmxIkTsu0MwwREzQJKaUBe3Km8p7aW0h49KAWk\nr1tvFW4bjVz/zqShgdJeveTnK/eaNKlz5llby10vgNKIiI69bgcOqL8eal6ffsrNd/VqStt/X1p+\nj2sBHAewT2H/5QDOANjV/nrCzVgdd9F0LhoOHDhAAUhe3bt3pw0NDV6Pq/a3HbRqFquVU6Ps2gWs\nXAmkpwPnz8v3/eIL4XZrK7B2bcfP0R0WC6c2UMu+fcCnnwJLlnAqpI5ATkXx0Ufc9QKApibuCaGj\neOYZ/4732WecEfWee7w6/F0AMz30yaWUZrS//uHVWXQuGXbt2iXb3tjYiF9++aXjJ6BG4vvjBQ2r\nF/5q0dvX8OHc6rizqK31bf6rVvl3/g0NlKakUBoaSumoUdx2YaH0vH36dNx1mzjRvyvzyEj+traV\nOeV+k/FwvzL/r8pxOuaC6VxU/P73v5ddmQOgW7Zs8Xpctb/toFyZr13rWi16y8GDnP+yVkOZvwxs\n77/v2/H33w+MHes/Q5/DGNvWxr1fvpwbX8ypU4Aviwh31+/sWe/HlaOhwb/jyTCRELKHEPItIcTU\n4WfTuWhhWRbr1q2T3de7d2+MHz++w+cQlAbQ/fv9M86xY8CKFcCTT6rrz7KcZ8mhQ0BiIrB7t/cG\ntn/9y7vj+JSXAxs3Ajfe6PtYfGOszQa89JJy3y1bgOnTtZ+DZYEhQzihHR8PFBW5rt/HHwMlJdrH\n7EQKAcRRShsJIbMBfA0gWanzsmXLnO+nTp2KqVOndvT8dIKIgoICNDY2yu6Tiwp1x9atW7F161bt\nk1CzfHf3AjAYwI8AigEUAXhQoZ/qx4oZM/z3KB4Xp/5x5o03pKoObxk3zj/zj4ryj9rjwQe1nTc3\nV/s5nnlGOMbixZTecIO0vWNe6h5F+S+4UbPI9D0EoI/CPi+/FZ2uwt/+9jdFFQsAumPHDq/HVvvb\n9kqAU+EPeSCA0e3vewI4AGCETD9VE29o4Lwq/PVPfvvt6i9aYqLw2KQk9ceK+d3v/PcZVqzwfh4O\n4uK0nTMiQvtNZMiQjhDSal/qfvD8F4ChAIoU9g3gvR8PoMrNOL58NToXObW1tW4FOQBaWFjo9fhq\nf9s+68wppccopXva358DsB9ArPujhDj0rFYrYDZzXhX+4uOP1eudKyqE2+Xl3p/3yiu9P1bMI4/4\nrjvXqq9uagK0POmVlQGHD2s7R2dCCPkYQD6AZEJIDSHkDkLIPYSQu9u7LCCEWAghuwG8CsAPyi6d\nrsiGDRs89nnnnXc6fB5+1ZkTQoYCGA2gQO0xLAukpnL67YEDhQFB/qClhRNK8+f7d1xPDBniv7Fs\nNk7vv2KF92MwjHaBfuutnJ5bTSDRiy96N6/OglJ6k4f9bwB4I0DT0bmIGTdunMc+vXr16vB5+M2b\nhRDSE8A6AA+1r9BVUVDACXCbzf+C3AE/HFwr3qYLkMt4+d57wExPns0KvPIK8P333h0LCJNaqeXs\nWWDECHXXoLRU+/g6Ol0BNSvzzMzMDp+HX4Q5ISQUnCD/kFL6jVK/ZcuWOV8Oa21RkfbzPfss5y2h\nlk8/VSeQsrKkbXfdpf48fMQqiiee4ML1V670bjwAmDvX+5uLt8KWZbnQeU+BTAWqn8WkKC1aRo50\nd9RWAMt4Lx2dzmHz5s1u9yclJQXGu0mNYt3TC8AHAFZ46ENrayl9/nnuVVvrXXh3796cYU6rQe/x\nxz0bGr75Rv5Yb+B/NkK4bf6+xx6jNDpa++f3NqjHYPDNwBgaKg3153+fvox74ACl69dL92nzwFFn\nJOqIF3QD6CVLQ0MDJYS4NX6++OKLPp1D7W/bHz/kSQBsAPYA2A0uj8UsmX40NNT1zxceTunkydr+\n8QcMcAkUhtF27NVXq/li5I/1hvx8SkNCXAJLzjMpIcE7Afjss9rn07u3d+fiv5Yvd43na4RrUhKl\nb70lvEE8/riwzwsvcFGq/fqpGVPdD74jXrowv3TJz893K8gB0PXr1/t0joAJc7Uv7oN594+/ZAml\nn30mXJGmp8v3JUS+/YEHPF80OQEVEeHdF1BbS2lYmGsMuQRWCxZ4dz2io7XPx9PNLzzc83kNBm4V\n/d13vt8c5Hz4+dcsLMx1zdTdONT94DvipQvzS5eGhgYaFRXlVpjX+pi9Tu1vOyjD+fkYjZwXxw03\nCKMxzWZp35QULi2qHP/+t+dzydkxevdWN08x1dWcURfgQujlqkY9/bRwOzQUMBg8j31OtXnZxfDh\nyvtMJnURqzYbd41nzQJOn9Y+BweRkcAtt0jbY2K46Nu33hIWpIiJ4TyedHSCjaNHj+KsBzexwsLC\ngMwl6IX5nDny7dnZwu1HHgF+/RXIyAB69pT2P3XK87nmzZO21dV55+NtNnNFGAAgLY17iUlO5vKl\nLFnCZQCsrga2bQOmTHE/tjc3GKXrGB4O/PwzMHiw9jG95ddfldMkxMRwVYjE7pCPPNLx89LR0cqL\nKnxyq6qqOn4iQPCrWZTUTa+9Juz3+uuufUp5zz0ZDpV05t4mPPvf/yg1m7UbLJOS3F8To5Ebc/Nm\n7qVmfCVj81//yu1X+uz+foWHa7+Ojvk5rovYmMu1q3sU7YgXdDXLJUlDQwNlGMatigUAPcD3fvAC\ntb/toF6ZJycDSh49CxZwq0qA+3vdda593bvLH+NJ1WKxyLd7m0WQUoAQ7cd5WiUzDFcn88oruddl\nl7l/emBZTk0lhyNHPMMAKmIffMbxnWmFYYDCQi5SuKYGWL8eeOEF7m+AnmJ1dASsW7cOrIrH9gNK\ntSP9jRqJ748XNK7MTSbPK87aWqlHBKWUzp2rPKan8ZTmozW1QkODK9eLI3+4WnJz3V8bOS+gu+7i\nXCvlzpOfTwWeRPzXp5+6+v32t/5bgSu9+vblfc5DufSqD66iuYe8yOolAvrKXCeANDQ00G7dunlc\nlQOgK1eu9Olcan/bQbsy//BDz+lnlfSrcsY1QN4Iyae6Wnnfgw+6P1aMxeIar6QE0FLP1ZMhVC4H\nyltvccWh09Kkq/T4eG5MMX36CHXpUVHq5+gtd7dnPsmrysOU96fgf5X/w5T3pyCvKq/jT66j4ycs\nFgsuXLigqu+kSZM6eDYcQSnMu3fnDJneMneufFShJy8Qs1le6AHaVS1mMydECeG8RbTUczWbuZfR\nyOVXTxZl0XZ30zl8WFoyr7pavtjHI48Ib5j33ad+ju5wdxN25El/IucJQfvSH5f65+Q6OgHAbDaj\nr8qK7Url5PxNUArz9HTfjmcY4Oqr5fd5CofnnpqltLZy9Ui1zOG55ziXwI0btRW5YBggLw/IzeX+\nas0iuWqVcNtxcxCzd69wu7lZ23mUGDnSdWNMTeVeBgPX5ii4MmrgKMExJSdLwDb7qaySjk4AMKjx\nIwYwR8mVzM8EpTC/917fx5gwQb7d3erTYnH5hsuxfLn687Ms8NhjXBrdOXO0uzcyDPcZGAYIC9N2\nrLg/w3A3FDFidz+x0NdYIMXJypVAfj53Iyoo4F7btnFtjpta1hBhIpxTTaewtWqrdyfU0QkwFosF\nJ06c8NivW7duYLwtV6aRoBPmw4cD11/v+zh1dfLt7jIPulOzAFw6XbVYLIDDvVSrzlzMb3+rrb84\n/S7LunTjiYnAwoWcB4hYlcUwnMDdsoV7vfmm9rk+/DA3Lv9mxH/voKlV+rixsUzmjqOjE4SoVbFc\nuHABv/hSVFcDQSPMDQZOgOza5X3dTT5K6hK51LQOGIYzAiqtSJVuEHL4ojMXo1BaUJHNm4XqJIvF\ndTOpqQEeekjZJsEwnF57+nTvUhJHR6vrJye4N1dw2eesDVasKVwDa4OXKSJ1dDoYLTU6lWqD+pug\nEeaPPMIJEH89kYweLd/ev7/yMSwLnDkjbywEuNzeamEYLg1B+kQrbnptJXJq18PaYMWOwzs064Z/\n/3tN3WGzCVfVZrPrZpKaqv7G4k164vfeU9dv7/G9krZmWzOsDVYMfXUo7tlwD4a9NkwX6DpBBcuy\n2LJlC9ra2lQf010p8MXP+LXSkC/4UqJNjmnTOC8QcR5ud+53nnTmn37K5SUXe5coUd9sxb7p8di3\nlfviQxACO+wY1msY9t67F0y4ujuXN4bJzZuBv/+de+8wqEZGckZVtTdMb86rduz6C/WSthZbC9bu\nXotWyt1NW+zc9l8v/6v2ieh4BcuysFgsMJvNAdP1XiywLIuMjAwcPHhQ9TEDBw7EeIfVv4MJ+Mo8\nROGMcnlRfIFhuKRbYs8Oh4CTw2x279/d2sqpTNQUiGBZ4M8frgVCXHdwO+wAgENnDuHZ3GfBNrOq\nVupmMxAR4fmcfPbvF247/i8jI9WPsWyZ/BjuUFvabvbw2ZK2JnsTjp07JmgrrdNLGAUKq9WKjIwM\nTJkyBdnZ2aqiGy8VKKV48sknNQlyAHj11Ve7rgFUybjpD6OnGIbhAmkcIfVGI1dn1F3/qChgzRrl\nPjYboKY2q8UCnOmm7Mv4bP6zSHo9CVnvZCFjdYZbgc4w2lfJDQ3y21oyLmZkcIbSW2/l/n77rXy/\nq67iMinm5koToCkxPWG6pG1O4hwsTFsoaPvY8jHK6jyUOdLxGZZlkZ6ejoMHD6KtrQ0lJSUo9sVq\n30VgWRa5ubm444478N133yEpKUn1sQaDIWBuiUAnCHO5cm/p6f7TlYtZt85lDG1tBb780n1/Srnc\nKIWFyrlKjh2Tb+djNgMh/Svd9jl+/jgA4ODpg/ii5AuP8/IWlnUJ2exsbW6SGRnABx9wf0ND5Z+s\nZswANm1SL8gBoKGlQdIW3T0aW6u3StoXrV+kfmAdzbAsizvvvBOneKlFhwwZgjRfrPYXGSzLYseO\nHYKnEZZlMXHiRFx++eX46quvsGXLFhQWFmL9+vVBqYIKuDCPjJQmvHqjA2ugi9UK7r4DluWKGF99\nNfCHPyinmv3DHzyfl2GAYbHdVM/z06JP3e7X6vPNd7G0WDj3SMA3N0mzmbvx8lVR4eHAzTdrG6es\nrgw7rdLE82sL12KXVfo0s+3wNq1T1VGJQw/8+eefC9oDqR4Qz0csVANxTrPZjMmTJyMzMxMsy4Jl\nWbzxxhvOp5MLFy6gpqYGDMNg6tSpsNvtHse12WwBfboJuAF09Ghg/nzOF/qZZ4ClS7Wt6LQivpbu\nrq3FAtjt3KukhMtdIkdhIVfk2NNvvYchElCZNfFM0xm3+wcMAI4cUTcWwBWRcGA2c7r+fft8c5N0\nGFKLi7lrk5vL+a+Lc+O4o6yuDClvpMjua6SN2HtC6uWi03FYLBaJHjg6Ohpjx45VPMZqtWLDhg2Y\nN28eYtx8+SzLoqC90ndmZqbHm8OZM2eQmJiIM2fOwGQyIT8/3683FLFxl2VZFBYW4oUXXkBNe+Km\n/fv3Y/78+aiqqkJ1dTVCQ0NBCIHJZHI+qVgsFpx3pBv1QJw7X2g/E9CVudEIOH4j2dnaH829QYsw\nN5s5NYLRyAk9R3EJMXff7VldwbLA4c2/4fKmqeDUBffVM9xXqpfC97hhGCAnh/tseXm+qbQcAUDJ\nyfJJzjzxTN4zbvcfP3dc0hbaAWsOQshaQshxQsg+N31WEkLKCSF7CCEKzq4XN2lpaQgJCYHRaERq\nair69OmD+vp6zJkzB1arVbJKtlqtGDp0KO655x4kJibCarUqqiiysrJw5ZVX4sorr0RmZibWr1+P\nLVu2SPo5jl21ahVOnToFu90Oi8WC2267DXv37hX08XblfvbsWYwbNw6TJ09Geno6Pv/8c4wZMwbT\npk3D9u3bYTKZYDQaYTabkZqaimpeAqQVK1bgpZdecm6rDRgihDhvEoEgoCvz1lbOy0KrAPCFjAzu\npuFgzBjlvgwD9OjB6dUzM5Wr8wDcTaG4WDltgMUCnD6m7ksHXPpzJcaMkQ/JV0KsY2cYTph3tqrP\nyrp3BYpj4lB6WujBEhoSii2VW5AZm6nanVMF7wJ4HcAHcjsJIbMBJFJKkwghmQD+BUDh2754OXv2\nLPr164dvvvkG586dw6xZs2C321FUVITs7GzU1NQgLS0NeXl5YBgGGzZsQGt7IEZTUxO++OILrFq1\nCuXl5YiPj0deXh5iYmJgsVhg4RUI2L9/P6655hoAQK9evZCZmYlu3bph69atYFkWI0aMQHNzM+Li\n4nD06FEkJSXBZDJh7ty5qK+vR1NTE4xGI+x2O2w2G0wmE37++WfFlXtDQwO+/fZbWK1W5OfnY/Pm\nzc4bQE1NDf7617+ioqICAKdCWblyJXr06OFcfe/YsQMlJSVISUnBihUrcPjwYed1eP/991Vd2wED\nBgTU7hA0QUNqcLjyWRusWH9gPVYWrHQGlZTVleHRLY9KPB/uv9/13mAQbstBKZcMimGA427kq8Hg\nPprUbAZ6D3YvoPmEG9xXbejRQ/VQAIBvvhG6UDY0cJ+ts73N2uzugy3EghzgXBav/PBKjFszzm/J\nuCil2wC4q2R6DdoFPaW0AEAUIWSAX04eJLAsiy+++AIpKSmYMGECMjMzMaQ9F4TdbkdlZaXEs2Xq\n1KkIbTfIEEKwc+dOlJaWwmazobKyEpMmTYLVasW2bcp2jjNnzuD777/Hjh070NDQ4NQtHzx4EAzD\nYNOmTfj555/xj3/8A5988okzQIevpy4pKUF+fr7s+A0NDYiNjcVNN92Ep59+GrNnz8bOnTsxatQo\nGI1GjBw5Ejk5Oc7tlJQUUEqRlpYGhmHAMAzy8vKQm5uLpUuX4tChQ2hra0NRURE++ugjfPzxx6qu\n79VXXx1Qu0NAhTk/a55W2GYWya8nI+udLMS+EotrPr0GD333EOJeiUNeVR5S3kjBC9tfQMobKYqu\nbGo8QiivOpC78P3mZqkvNx+GARYvcPMYIGJy/GS3+93llJGDUpePPctyQVQ2m3ZvFn/CNrPIqc7x\n+viyU2WBTMYVC4CfOb62va1LYLVaYTKZsGTJEuzZswdWqxUMw+BuR8L5dkJCQpz6YpZlce211zqF\nK6UUH3wgfLCpqqrCmDFj8H//938ghGD9+vVYv349UmUqctfV1SEhIUGQfbCsrAw9evRwCsHRo0cj\nLS3NqQZKTU2F0WhEdHQ0HnvsMRw+fFiidvn0009xrt0Ht6GhAcOGDUNycrJTQDueHvLy8rBp0ybY\n7XbMnj1b4FvPMAzS0tKwdu1aREdHO8/54IMPqq7pOXHiRFX9/EVA1Sz8rHlaeTb3WRw7L/UJtMGG\nO76+Q9C2/Kfl+Oj6jwAAr7ziarfbgVdf5cqNKWG3c7nLx48HPOWeLy935ecWwzazePnINaoNoEsm\nLHG7/+qrgZ9+UjeWA0c4vsXiuvE4vFmU1EMdybdlCo7qGsirysP8lPl+mI1/WcaLsJo6dSqmKtU7\nDAJYlsWoUaNQ175aYVkWkyZNwvbt27GGF2TRt29fjBo1Ch9++CEA4O2330aJwy3KDY5sgrR99TR/\n/nxMnToVv/zyCxobG7F06VKUlnJPYKtWrUJpaSn+9a9/oby8XGBoBOBcJRcXFzvbi4uLYTKZ8Oij\njyIlJQUtLS0YPnw4tm/fjrCwMHz88cfo168f6urqYLPZcN999yEvLw8DBgzAhAkTwLIscnJyUFlZ\nidWrVzs/U1FRET777DOMGDECR44cwX333YczZ84gNTUVn3zyCSilmD1bGuymhFrdupitW7dqyv3i\nIGjC+T3xav6rivsOnT0k2D54ymWd37xZ2Fe8zYdludzhs2ZxHh/JyZwHiBJff62crvejXevQSE+r\nFubfV3yP7KHK1mA5A+iIEdxTxIEDwKBB0sRY06Zxf81mLidLUZHvSb984YfKH3we49ejv4JtZv2p\nO1eiFgA//+Tg9jZZlonDZYMYi8XiFOQOampq8O233+JIu8tUSEgICCHIyclBdnY2DAYDysvLERIS\nosotz8GDDz6Ibt26ITMzE9PbVz5Tp05FcXExfvzxR8yfPx9tbW1ISkrCF198gW7dpO68DMNgAm/1\n4Xh/6623Ys2aNbDb7Thw4AAGDhwIm83mvIk4KC8vx5AhQ9CrVy/ExsaioqICLMuiZ8+eeO6559DU\n1IT9+/cjMjIS9913H9ra2hAWFuZ8Ajl48KBTn56WloaioiJV18DbnCzixcBTTz2l6riAqll8ecS/\nAOVlsiNM3sGAHi7Vpnh17e78Fgunnmhr41aw8z0sAN2t3D/d9V/3B4vYfXS32/2ZmdKAneZmLlf4\n9u3StAUAUN+e/oRhuJtYWJjv3iy+EBri+9ohtyYX2e9m+0t3TqB8u10P4PcAQAiZAOAMpVS9ESSI\niY+Ph1EUuBAXF4e5c+c6VRoJCQmoq6sDpRSVlZUob0+eRAjBCy+8gK+//hpRUVHoo+S/205VVRVm\nzZqFrKwspyeLQzhPmzYNra2tsNlsKC0txXXXXYerrroK48aNc6o73HmvpKenIz09HUajEaNGjcL6\n9eud+0JDQzF8+HDnvpMnT6KoqAiLFy92ZjFsbm7G2LFjsX37duTl5eGzzz5zHm+325GQkACj0eh8\nWmAYBhs3bkSkypwYcqqljiSgwtyXgJUwqK/QsL58vVNvLs5BXlGhnFvFbOZWug7XxPvv5/J/K6GU\nXREAJiVq8yXs2839IxnDcLne+XTr5nIVlItWPcR7YOnZkzPaigV5INPNVpyp8HkMO7Vj7/G9HiNm\nPUEI+RhAPoBkQkgNIeQOQsg9hJC7AYBSuhHAIULIQQCrAXgwnV88VFdXw8bLKDdw4EBs377dqUfO\nzc3FOwo5K2w2G95++23ce++9OHv2rCBqlBDhfdGxbbPZYLFYMHPmTIFe2mw2O4VxYvs/GqUUZWVl\nuP3222GxWDBmzBhkZ2cr5op56aWXsGnTJuTl5SE5ORmEEISGhiItLQ0//fSTU0ceFRWFAQMG4IYb\nboDZbJYIaYcB2HEzEx/v0OGXlJTgzBn3MSEO9rszqnUAARXmnjxA5HB4sPSOUAjHVODP3/8ZgHwt\nUKWiCwzDRU5u3sytYGNigN27gbvuku+/Y4fySt9yapdqFQsAHDxz0GPiLfGNnh8YJFcXlK924Rt2\nHVgbrIh/LR73bLgHca/GdXgOlKhw/1WMvmP9HbLRomqhlN5EKY2hlIZTSuMope9SSldTStfw+iym\nlA6nlI6ilAamkGMA4K/Mw8LC8NNPPyEmJsYZVJOWloa3335b8fiysjIcE+W0CAsLwyeffOIUlGaz\nGd988w1obYU8AAAgAElEQVTMZrPTwGm32wWeMXyvkY0bNzq9ZEJDQ/Hll18iPT0dFRUVsNls2Ldv\nH958803k5eVh//79uO2229CnTx9ceeWVmDlzJuLj45GcnAy73Y6hQ4di48aNiImJwYQJEwQeJfxz\n8oW03D6544OZgArztjauOIJa2GYW49aMw+R3JuN4k7Yn3O2HtwPgws3FuCv8IeefLV4R81Hy/ba5\ny6Urw+j+o5G1NgvZ72Yja22WrEAXB+Xxfeblanzu5mluGho44y7/5rP8p+VOV0EbtfnV9U+OGMa/\nAQZP/fSU8waoo57q6mqnr7jNZsPhw4dhtVoxevRoZGdnY+TIkU6jp1psNhvi4+ORn5+P3Nxc5Ofn\nY/78+cjPz8f3338Ps9kMQggGDhzoNGQ6bh4mkwkbN25EU3uxW5vNhgULFiCEp1eklOLRRx/FlClT\nYDKZ8MEHHzh12iEhIXj88ced/aurq90G6zhW4nJC2t0+ADCZTIJ5KTFkyJCApb51EFBhPmKENuNb\nTlUOyk6VSXTiajjXzLkmydkOlKrhsCx3w5kxQ6jf/81vlM/zXwXVeP4ReR9YJXqE9YDlpAU2aoPl\npAW/1ErvOIsWuep7hoVx2w4YRnrjcqiYWBaYMZdFU78dyJrGOj+X2LuEbWWxdtdaTfPWwsnzJ/06\nXvHxYsS+HIusd7I8d9ZxEh8f7zQ02mw23HnnnZg0aRIqKyths9lQVVUlMSJ6IjU1VaCycAhDhmEw\nffp05Ofn4/nnn0e/fv1gsVhQVlbmVKEMGjQIb7zxBuLi4px+4K+99ppABfPQQw9J1DgAt4o3mUz4\n3e9+51SRpKSk4Ny5cxK1jD/yvlRXV6syfi5fvjzgK3q/CHM1odFa2WXdhb98/xdVfftESI0wreBW\nHnK2ikmT5MdxBKw5DKAO/X6pm5TaSi6nQ3oOkd+hwEs/vyTYLq+XVuuIieH04G+9xf0VR9KKUyNM\nbnddL9jDYn/2GOAPWbBMScXWQk4/LrdSfjn/ZU3z7kwqzlaAbdVzbmuBZVnMmDFDkFukqqpK0Xfa\nYDAgXi7VaTshISFYuXKlxzwqDMPgrrvuQlFREbKyspCSkuJUoTQ3N+O9996DxWKR+IHn5uZi9+7d\n+Pvf/46RI0fCaDQivH3VkpSUhO+++07iN97c3IzZs2dj8uTJzpQEFosF6enpbvXvaq6d2vD8zsgF\n76+V+bsAZnrqdOCAOgPoLusujH1rLMpPqys/dLpJOZDPbJbmMFfK6e1QVTgMoI6nCHcxAkpG0Jhe\nvqkU/ntAfskfE6OcE+X554Xbzz3H/S3CR0CfCu7bjqzF3TsvA9vMIq2/9DHp5Dn/rp75xDJdJubm\nooK/IrVYLDhw4IDHY0JCQtC7d2/YbDbUuikGazKZcPvtt6tahe7fv192xT9s2DCMHDlSdlXv2Obr\nsysrK/HAAw9gyJAhuOKKK5z9KaV4/fXXUV5ejra2Nuzbtw8JCQnIysrCyJEjUVNTA5vNhqKiIuze\n7d57TExtbS3S09OxcOFCz50BTFJaMXYgfhHmKkKjAXAGOzVqlue3P++5E//8MtmsHKt1huEChfgo\n1vIMY4HBO/DtZlbgwjdrlvK5b79dvr30hG/GRDvVrlo6JcrVdbr9G/mo+G2XMZYAx85b8UvtL5Kq\nPh1NZIS8S5cBbso76fgEy7KYMGECsrOzkZWVhb59+6pyrSOE4HT7D0ip3qXBYNCUKtdsNsNsNiM0\nNBQRERFO98GffvpJ1RgO4R4TE4OXX34Zx48fx0cffYQdO3Zg/fr1GDlyJCIjI51G2GHDhjnnbjAY\nkJiYCKPRiJ49e2LJkiVuV9mOG6DVasW///1vxMfHC5JveeJLT4UTOoCgDBpK6JXg8xg9jK5kJnz1\nQ1iY1JAIOIyt44E7DmDxvhHYOaEAAPcDU1KzEMJVMhJTVleG6nM8v0AKgADPXvEs/t+P/0/V/KvP\nVmsOjhEXAXdsn26U3mfrGuvQr3s/SXtYqHoXUK3wvxM+w3sPx4HTnleLOtpgWRb//Oc/nRGOFosF\nl112GRrEZahksNls6NmzJ5qamtDW1gZCiHNV7XD/M5lMmox8/GjOuLg4ZxIvb3TLRqMRL774Iq6+\n+mrY7XaEhITg888/x7XXXguWZZ3nmDNnDkpKSpxG1pqaGphMJqxevRrjx4/HmjVrEB0dDXP7Y7nF\nYkF8fDyuuOIKHDx4EKGhoWj2ohiu1pW/PwioMC8pWYYlS7hKPu5Cnv9X8T+fz9W/R3/eeV3tLS3y\nmRu/LclBWX0pYADK6vdj4/6tuHE0FzWUmyt/Dkrlx3pzp8j3sX1VPHXYVGdRZ0+U1pci+91s5N2R\n53O04+lmqTAvqy9D8QmpzottY2FtsCIm0r+eJ2V1Zfjz5j/L7pOrOuSRQwCqfJpSl4ZlWUyePBn7\nRCHMagS5g27duuHjjz/GwoULnYE2AKeCWbVqFW688UbNgpgfzekuF7oaevXqBbvdDrvdDoPBgEGD\nBknOwU8FwDCM85yPPPIITCYTfvOb36CtrQ1RUVGIiIjAyZMnYTAYBN4+Dvg3NE9kZGT49Nm8IaDe\nLOnpy/DKK8uwbNkyt7krbHZtbn1yPH3F05r6r9+bJ9h+t+A/zvcKydkAAIcPS9syB2dK2sIN4UiL\nVu/KQ0Gx7/g+FJ9UH2Uljh52bPeO6CfJq3604ahiOtpVv8qEk/oA28wi8y3pNXHQancTfaXEMADT\neC8dARaLRVUeFXfU1dWhtLTU6TIIcN4jZrPZK0Hub/hBR+KcLg7cuRryo1dZlsWJEyecN4fhw4c7\n1UEOv/nNmzerLjZx6623ev/BvMSfwtxdaDQA9aHkI/opKbU5bk3zfKFK6lw/ZJPJ1W40SoNvAKC6\nVegK+P2xD51BKe4q/OyUVj/DlLgpwgYK9A7rAyacwa3p6r9kCoqzF86q7m8yucrFhYe7PmcPSH0x\n39z1Jo41yuvMdx/z7yPi24Vv40yLctScls+oo474+HhV/tDuoJTiiSeewIgRI5wCzeE90tmCHHAf\nAKQGs9nsdGd03Awcn/Onn35CXl4eKioqnH7z48ePd+au8cQpsQErAPjLNVESGi3XT+21rjztvhBy\nz249PY7x8P8edgpjvsBtbeXKvomxnqsRGAkB4IFNDwCQpgTgM2WKtK36rNRQ0tzCPW28Pvd1j6H7\nfP703Z9U962u5gKDAO6vw76TGK3Ni2RQj0Ga+ntiTeEat/sdbqRiRvfvksV9AkJ1dTVa3P1wRQwb\nNszpx82/CbS0tOC+++5zCrTp06cHhSB34CnIx9OxjptBfn6+M+CJH/3JjwLNyclR5WNOCOmUYtj+\n8maRhEb7Mp6nSEG1equp708F28xKDJhyBs2689Lk5Q6dsoIxHyEhLl9uPk5hTeFUb7x4JZeLlwln\ncOihQ5g82H3+cgeHTh3y3Kkds5nT3xMi9Bzq21N9LVIAql1CPeGIzlQyfDqIMsqH+R9p0FD0VEeA\nO99wORoaGpz/V5RSZ2h9REQErrvuuosqrF0LYvdHpc/JsizuvPNOVWNedtllnXKtgrLSUGS4e9ep\nA3XqPB/YFhZflnwpyQcjp/bqFiZNvdnYwhl9Fi+WH58Q+fQEb/z6RnsH7s8VcVdhUeZNzv1MOIMx\nA9UVrjCGGj134mG3uzI/Oiiu05bdLLcmF3lVeZ47ysCvBnXZmssw+Z3JKDpR5PaYm9Jvkm1njF1P\neAQKvhtdSEgIfvvb37rtX+9IsQluZfn0008jOTkZFRUVPhsquwIWiwUnT6qLwSguLr6og4b8SkFt\ngeK+UIR6LD3G553d76CfyANPvL3LugsnL0i/KBs41cif5Z0wYLfL3xg2lgsTthw+J11dR/dUyCkg\notnWrDpfSk6OKyNkaSngyG+f1k/7I9/SH5dqPoZtZmF6w4TJ70zG+LfG48CpA7DDjhaq/LgfFhKG\nJy5/Qnaf3Heiow6Hoc5RKWiYUnXydvih8na7HUuXLkVYWFiXXI17g5ZCE+fPn3cmEwskQSnM3WVI\nvCblGiT0Vu+Hfuj0IYhTLvfmDV9WV4axb8k4nsMVzMIw8ilmKQX+J+dFKdYCyWiFRg0YpTxpHq20\nVXJzUGLvXvnttrYQ2Tm4o2eYZ7uEmJyqHBxhj8AOO2rPKUcNOggPCcehhw4hJjIGRKQHIiB+zbJ4\nKcGyLK688koAcGZDfMFdeS1IVZeOtLXehr53NbRU/klISLh4deb+hG1mFQs1EBCsnLMSS7PVrxpP\nXzgNcQK4f//b9f7h7x5WPHZo76HO96tXy/eRq80ZYYgQbBtJhKTPtGHq/emUQvs94fj/jG7j+byq\nFOreRKD+csRNOkoZhkUNc/qzh4h+igYYkBWnJ9DyhoKCAmcu7SNHjqiOXJQrNMFPWXspM2/ePNlE\nX3LMnj1b15kDnIqlBdLH8sxBmTiy5AhiImOQ3C8ZI6PVFX84ZzuHJLNwZZGe7nq/44hy+tQWm2se\nGRmcF8wVVwj73C9TsiB9ULpgu3bXSEnecyackQgwJapOV6nql5Qk3E5O5v6SyFqJp44nUvtqr5Ky\ns1bGT9MND0580Ple/DTWp1sf/OOKf2iegw5UF08Q061bNyxfvlwgtOLj4ztllRlsxMTE4Ikn5NWB\nYi677LIOno08QSfMy+vkPSk+uO4DQVRik61Jtp8cB1jhipFf7q2bUWr4dBBmEIa2Z2QADz0k7CP3\nf/OXie3ZHttXwQ2bH5ZNMDawx0Bpowz1F+o9d4Jy0NDuEz+rOp6P0aDN8Mo2s/j+kMxjigI9DT1x\ny8hbnNtLsoQFrf808U9I7peMCTGdUHn6IqahoQG33Xabc3vo0KGqj62trcXo0aMxcuRIzXlTujpW\nqxV///vfVfUdPHhwB89GnqAS5tYGK/646Y+S9j7hfZDcL1nQJpdvRIkfLUJvih94dYXHDZJRhrcz\nM1GaCFKcQVEuo+KppvaAgfYFzrDU07IJxrob1RV8VWsAlcvNUlZXhh+qtRdSPtuiLZBHnBvdE6Nj\nRgvSFDww/gHERXJGu7jIOCwez7kQHW6QCbHVUSQ/P98Zem8wGLBs2TL0kiu3pUD37t2Rl5eHvLw8\n7Nq1S/dkaWfdunWq+g0YMCDgRSkcBJUwX7d/nWwGxL9e/ldJW1Ob+pU5O+zfgu1581zvF44UpbTk\nnf69ve9JBOmCBa73oaHA9ddLzyf2Wf/bs3WyAVOmaJOkTa7W6bHGY6oE+okT0u3Hf3jc43Fy9IuQ\nJuFyx2eWzzx34vFEtvCRlQlnYLnfgh2LdsByv8Up6If2Gqpp3Esdu92OsLAwhIaGIikpCXfccYdq\ntUtCQgLGjx/vUyBOV0VtROef//znTrtuQSXM6xul6oSosCgsylgkac+IlSayURJAPZg2Z2X7sDDg\nqqt44wwUjcPTKZ9rPYetVVsFu/m50G02+dzoYt3x7uMyIacAdh2TlpWUU/vYqA2byjfJjsGHn6qX\nEGDmTGDfCe/qhbz484uaSsjtO67tPKEGaY43JpzBhMETBCv2RWOk372OPCzL4t5770VLSwva2tpw\n/Phx1QF2cXFxQROmH2ywLIsXX3xRVd+bb765g2ejTFAJ869Kv5K0XW+6XjZr4OQh0gjK5D7JkjYA\nSO4+0fmeUmGgz2fF7leUYiG6lldVjVLgXZlYV3HJtzXfFsgWfmabpI1/miAfvr+uxPNjXn09nDct\ng4HLbx7dg+fPLvN/PTJ6JKbGTZW0t9nbJDcyd2gtM7atepuqfgvSFqBfN21PCZcqFotFUEjCkY9c\nDW1tbbogV8BisQiyRioRERHRqdcwqIT51PipkrafD8sb7/p0l7pRTUuQd/cb1M8Ih9pw+HBXmDvb\nzOLZn54VdhbJpJxDOYLtG28U7pcLrBMbZy80N8kaQG8ZdYtg+/aRt+PPk/4sW5HnYP1B6QAi4uM5\nIQ5wf+PiRJWDZDxZ0vunK4bNFxxRDt4So1XH/vZu5ervfJhwBmuv7ri6pF0Js9msyeDJ5/jx47oL\nogJmsxlRUZ5jHpqamrBRqcJ7AOh0Ye4I/y6rK8PgKKkV+Ch7VPY4saCZkzgH9192P8IN4ZK+JSdL\ncbZd1lRWugo1F9QW4AK9IOnPZ3jf4YJtpWo+fBamC/XwA+sXyhpAxUbdsbFjwYQzePkqaR1OuXQD\nYkpKXGXsmpu5XOup/dy7GFadrhLkfuez9/he2XY52mzqo3IBoHuoOuMvwPnkR4XJ/zOJg40uZRyJ\noxIS1AfVMQwDg8HgzCCoI8+FC+7lhIOCAvULIH/TqcKcbWYxce1EZL2ThZQ3UvDolkelfVrk9baP\nThL2/fsVf0dMZAwqH6zE+EFCazJ7IhqOHPPNzcAXX3Dvf6z4UTqwKJV6ZqwwD3edKB+XeBsARg4Q\n+sCveXqkrAF0QeoC580n3BCO60ZcB0A+6+LeY3s16bAdXDOCVwpJRhNyuvk0xsXIe/Tk1+TD2iCf\n81xMU6t6gzQAPDfjOdV9mXAG//2dfOAUY/DusZYQMosQUkoIKSOESH54hJDLCSFnCCG72l/qnIw7\nmZiYGOzZswe33HKLx75GoxE7d+7Etm3bdH25GywWi+oMlPfee28Hz0aZThXmlhMWj8UXDES+PmRG\nTAYK7yrErSNvReFdhciI4QyZMZExuO+y+wR9x/e9UrDdv30h+uavoopAFIDleoHQuyHtBkGXMlFp\nz3IZt3hTP5NzxRiCMIyNk18dO24+b81/C5UPVjr96KO7S/O2nG87jy9L3NcVNJk4Ay/gymfuyR+/\nu7E7YiPlU+Seaj6Foa8OVSXQtebOHhSpLc2uXLUkAG7zvihBCAkB8E9wRcjTAPyOECKXRD+XUprR\n/rpoIpgYhsHdd9+tuP+RRx7B66+/jqqqKiQnJ+ueKx6Ij4+XjY4Vc+211yI5Wd5uFwg6VZjLqUTE\ntFLlKjQZMRn44NoPnILcwekm4T9+S4hw2xFMc8Em8+jUVyidxX7OEaLI/HCZj7DTutPpYmlHC7ZV\nynuzAJxAvzPjTkFAlJJv9Yr8FYrjAPL5zD8p+sTtMTeNvEkQvCOmlbZ6rDzENrNosqtfmSf2TtRU\ndQkASuvkC7G22LULcwDjAZRTSqsppa0APgUgU8314tThWK1WTJ8+XXF/QkICFi9erPuQq4BlWWRl\nZalyTexsNVWnCvMP9n3gsY/aKEk+Hxd9LNjOaxAKNIcwDw+RkcR99ws2D58VCtbf/EbYXa6gc2m9\nUPA89PdSWW8WJVra5AXUvrp9blPTms2AIzmeI595bYP7hFcxTAxiImOwJHOJYh9PnieWExa3+8XM\nGj5Lc11TR0CRmGFR7rMBKhALgP/FHmlvEzORELKHEPItIUQaFBBksCyLLVu24Mknn3TWsJRj0ybP\nbq46HBaLRXVum8TExA6ejXsCWtBZTEqfFI99eoWrj15zEBYiDLzpy4ThVF/OeDliBOAI0OLnXgHA\nrcMMrYL12KayTbh9zO3O7fp6zoebUi5oSO6GLRY8x8viUFwMTPBDZPqi9YtQ9mCZ547txEbGovKs\ncuUmy3ELbjTfiKemPYU1hWtwvu28pM8RVt7bxdpgxYbyDW6jaOXIjsvW1B8A+vWQd0/s37M/Ks5W\naB5PBYUA4iiljYSQ2QC+BqD4DL1s2TLne3fFyjsKlmWRmZnpTLDljvz8fLAsq6tWVBAfH6+6kLPa\nfOee2Lp1q6YsjQ46VZgrPTrzGRI1RPO4vSJ6ibajYDBwQtjAU8H36dYH1kaRPlj0YC2uh+lw/2tr\n44S5XD7z7mFCT434mO6y3ixKyBlAHRw8reyiaLEAh9pTpx84ABQXA78z/w55h5VX847Px4QzuG/c\nfXjp55ckfWrOSCtwWBusiH81Hm20TXXCMIDLR++NMM+MzYQ52gzLSddTAAHB8zOex4wPZ8gmZ3ND\nLQD+Nze4vc0JpfQc7/0mQsgqQkgfSqns8zZfmHcGFosFZWKDjgJ1dXUoLi52VrDXUaa6ulpT4JU/\nEC8GnnrqKVXHdaqaxZ2QcfDMjGc0jyv2gDl2+hzq6jg9skPIAUr6VqHB9fj544Lt6mo4PWPa2uQr\nDYlzrqxc0V11/VMA2FKxRXGfXLoDB3Jqltgo9/U/S064Cl9HRshXeGpFK74qEQZ0rfp1Fdoo545o\nh/p0uTbYUNMgc9E8wIQzyF+Ujx2LdqDwrkI8NvkxlP6xFNlDs3HoIfWl9dr5FcBwQkg8ISQMwEIA\n6/kdCCEDeO/HAyBKgjwYMJvNGDBggOeOAIYPH97p+t2LBbPZjNhYdTV0u4sz3QWYgApza4MVawrX\nuLwjbO77Tx48WWLcVENEqNBKGRrehOhoLjrSZHIFDclmBmwS1qsUh9fzA3OUVuZ9Itot3+1yt1e4\ncrENORxCUgl3Lor8lMtsC4sHNz2o2BfgUgU4cFer8+7/Cr0jNldsdjuuO7QaPx04wv0zYjLw7PRn\nnX76fOOxGiilNgCLAfwPQDGATyml+wkh9xBCHB90ASHEQgjZDeBVADcqDBcUsCyLOjk/WRHdu3fH\nt99+q6tYVMIwDNauvTiC1gIqzONeicM9G+5B3CtxKKsrw7HGY277p/T1rFOXY/bw2YLtnScKMPP2\nXZg5E9i4Ec5V8vWpwixZ5EJfwCj0yjhzQZikiO8xYrPJr8xf3t4e9NMuWOc9/bImA2hiL/eGFKUw\ne4uFC4oCuCeQVT9+4VZlAwC9wlwqqe8qvlPsV9dUJ3BRFD+xqIWCajZ+dgSU0u8opSmU0iRK6XPt\nbasppWva379BKTVTSsdQSrMopZ0XDeIBlmUxceJEVb7QjY2NmD17tl49SCVlZWWYxU96FMQEVJg7\namraYEPWO1k4e8F9CHjhMWWXPnfIVSr6tPYpfP89MGcOFwFqbbDinzv/KehD8/8E2IQeLuGhwm2z\nGYiO5lbADlWGmJK6EsF2Q2iJbDi/EjMSZ7jdr1TRR6xmyT3r2Vto0tBJzvee0gqnv5nufCqYlzTP\nbV+dwGGxWHDkiLyRWo6qqio9dF8lr732muq+l5SahU/9hXqEhUrTvfKxtXnQwyjwXaV0hdnS2gq7\nnQt5Ly6GrO90NNMTOC6Uzmn9pdLaU/WobqFC1UyEsYcmA+j9l92PUDe26cozyt4p/LkdVfBC4TOI\ncQXvnGBPuOnJ5Wl/NpfLZfPr0V89jq0TGMxmM4YPH+65YztDhw7VdeYq6dbNcxoNABgxYkSn5TF3\n0KkGUHdVfgDtulAHcnk8SGSdYDW98YA0IU76vG0w9hAZGEWC22Lh8oRTKjSm8hkcKcwxM2tyjCYD\naExkDKqXVOPFK+XTbirpti0WoKLdS6+0FKg/L5OfV8TogaOd71uh7Jvs4JX8V7DLugs7rdpKxDmI\nCtWLNPsblmVRUaHOPTMkJAT/+c9/dJ25StSG8S9btqzTr2mnCnNjiPvSZOFhniNE5RgzaIykzXAu\nXrAtDuwBgKyYGehdK7Rz3WAShvObzcDAgVJjKp+kfsJinEm9tev+YyJj8Jesv+B3pt9J9in558fH\nc0ZZgPvbr6f7EGTGyGDq0KmuBhUeWE1owti3xmryYOETE6VHHfqbDRs2wGZT9xRrt9tx44036jpz\nlahNnPXYY491+jXtVGGulK3PwZNTnvRq3OtSr5O0tfX9GZRyK1YldWF/QxJC6oXFmPlqCIAzni5e\nDEyeLDSm8kkfIBzjs5VmTQZQPsumLZO0rfhZPqy/uppzlwS4v2f4NgkZQZ3YK1FgjJQrGOFvXpn5\nSoef41Jj3rx5MBrV12zVdebqyMvLwy+/yNunxBw+fLjTr2mnCnOlJFoA8OVvv/TKLREA5ibPlTZ2\n56KzHO6EYvdFAoKVj47HsZ2ZMJ51Bfrd9vVtAi8OlgVWrQK2bXMZUyWnEqV3rT3YT5MBlE9yv2SY\n+gojyZUSY4lX5q12XjSnjJ5f/ATx3HT1mQzdIdb3X598PSbETsB3N32HmUnSuqo6vhETE4Pnn39e\ndl94eDhCQkKQnJyMhIQEhIaGIi0tTdeZe4BlWcydKyNHFIiOju70a9qpwrzitLyer3tId1xrutbr\ncZlwRmpAJJxawOFOKM4lntJzHKrKGKCFQduvrlJlzbZmfLH/C+e2xQJYrRAYU/mwzSwe2/IYt0EB\nnDQhNWq8JgOomOnDhEmTsgZnyfbjBzTZbEBrm3tVyBXDrhBsD+vjVZ4TCQQERsKtFI3EiJVzV2LH\nnTt0Qd5BlJWV4eGHH5bdZ7fbsXr1auzcuRN79uxxFmvubP1usJOTk6NJbfLaa691+jX1izD3lBta\nibom+SAHpUhELYQTkb7dbhQYQMUh5dOGXe4suUYuCFPQMmGuL8lsBmJipDrzr0q+gukNE1bkr8D+\neld+jPFnn8P36xlNBlAx4uvRM7ynbD+zGXBk4Bw2rgzn4d71U1xhqG+3vt5PkocddlT9qQpvzX8L\nVX+q8tqQreMZlmUxbpxybpyUlBTceOONYBhGL9SsgaqqKtV9hw0bhtmzZ3vu2MH4LMw15IaW0CtU\nPomW1nqSsmN3F40d3gga5rrTJvQRVmPp0ZLoXNXaewiFHD8lLcMAf/kLcNllLp35VyVf4br/XIf9\ndfuxLHeZ4Nhf9jQqqmPUIq4VqhQByjDAv//N/U2653GP45bXC9P9ugsa0oIBBtnUvjr+x2KxKK4g\nQ0JC8Mwzz+jC2wsWLFgA4skHuZ1gWJUD/lmZq80NLeFM2xnZdkp8F+bi/CgwAJj6F6cBVJy29WyE\nBfHtDi8944QJi/hCj2WBl14CCgqAmTO57T99L1+EGQAwoFBWHaOF8tNCoauUbItlgZtv5v5utHgu\nmHyuVei6ODRqqNdz5CO+Uep0DCzL4sCBA4r77XY7li5d2uleFhcjDMMgNFSdQ8Bf/vKXoLjG/hDm\nanNDq+aZadqTa4k51SyTEynjPYSGAn0Gsgglwi9qWJQrhP78gB8E+3445NouKAAcwXYWC/DLL1wl\ne3Wkl9EAACAASURBVEWiixVdGNUizmWS1k9+MIvFVQmJGt2nSgCA+cnzBdvThk2DEeq9IpRYM3+N\nz2PouMeR8vaOO+5w26+0tLTTvSwuRr799lu3OeH5lJWVqfZ66UgCmwKXX+h+KAAFe9vQPkN9PtXz\nVzyPO7+9U9jYZESbgcW0z0fDekEYRfnVvo2oruYKNIjzXNmpe0PiyXPKeYyTojIUXRjVIo5ClYtK\nBTideWIiF8wEJUchCqdnS1JfoTcLE87gs99+huv+I3XtVEvmoExkD9We4lYr3uZ87ipYLBa3q3IH\nlFK/pWa9VGBZFmvWXHwLEn8Ic4+5oZ1MUzdgY2ujr3PConGLsPyn5ag5x8uEdToVsZNyUH1BGg5f\n3rgLBkO7N0ir0MBY31gPtpkFE84gMxMYPJhbnaemcoUuWrcp38HLv7gVcwqAvDzvBbq4jFxNQw3K\n6sqwaucqTBw8EXOS5jj9xRXVfBRAYx+gcQAM0WVIjU7F+Fhp+PFARntlJz4FRwOTj8rbnM9dBbPZ\njN69e6O+vt5tP5vNhv379+sl4lTCsixGjx6NykrllBlikpOTOz2UH/CPmsVjbmitSPTdXiJIWEWA\nfmQUZt+dK9s3PsLsNIDiB6HPbittFWQqdHi9OFLhRkBUGJQPc9hnnbn1rLCARsnxEqS8kYLXCl7D\nwi8WYuyasWCbWVgswEGHOv28jHH5/c0w/1yA7xduQ/6ifNnsheb+ZowaMArGECN6hCqnxNXpXBiG\nwfXXX++5o44mCgoKNAlyQghWrVrVNQygSrmhvZtMCMzRZtkVozfsO7bPtUGBOuzCez/I67ZO0xpn\nbvKBw49K9udWcTcBi8WlMy8tBTYWlKEJbooZuwn7V8vmQ8Lc4d8c+EawXX6qHL/U/gKzGXDkW4o8\nIbRBX93/jzD1y0B+DoPpKRMU09Ay4Qzy7shD7h25WD5tufeT1ulwGhs9P8GGh4cjNTU1ALO5NBk4\ncGBQrMoBP/mZy+WG9oZXZ72quGL0Bkn+lcG70RQlrwYY2X8MHB6R9QmrJfu31XDeIeLiFBvrX3U7\nh6syh/isM689K9RaXbBfkPQpq+csn4QAiM1Dw/D3BfuTe5sRFaVuHo4iEMP7qM/E56B7SOemAb1U\nYFkWn3/+ucd+drsdNXJJ93VkyczMREKCem+sV155JShW5UAnR4CKuX307X4tXHCujed6R8B92nD5\nLGg3D3gOh9tV060ywrL8FOceKC5OceG8e++P/1V/47OfeYTRjRqnnRX5K2CxAOXlAObeKwnf//nk\nFpw7p20etay86cMdA3qqK12m4xsFBQUeM/oZDAaYTKZODzO/mGAYBvn5+ar7JyUlee4UIIJKmPu7\nAk2fcJmsgTIGwqjQKMydkIwh7bWjmbPSQrcXmjkBbzZzuV0I4dQnreFSlYyAHrU+68xHDhjpsc/B\nMwfRFr0LI0YAiKqSJNbav34eioqA7Gz1Av3L/V9qnuvVI67WfIyOdjypWHJzc7Ft2zY9dN8LtBSZ\nuOaaa4LCxxwIImHeLURdEngtTImfoqqfoV1v4lCz8CNFHdAQbifDAM8/D0yYwHmoNFMPOcNjCnzW\nmf/jin+o6jf7syl47wsrQI2Sm9bpfM5YpuXGwq8PqgYjMeL/Jv2fpmN0vMNT/vI9e/boofteojbt\nLQAcO3YsaPz4g0aYj40Z6/cxH54on3xITP+e/QWGzXNG6T9K/+7CdL0OwT96wGhJXz5kUInPOnO1\nqWnPt53HnK/HAdZMQfsVQ65C/EBuAkql7uQYO1Ddd9LL2EvPwxJAdu3apZhYy0H//u7TS+vIw7Ks\npviFYMpAGXBhHttDPjg0tb//Le451TmeOwGIZWJhNnP+4wBgCJNGdDa0NADgVBSPPspFgmZnAz8e\nknd1dEBDGnzWmZv7m0Hk9EMynGg8CgzYJWg723IG3qS7EQf/LExdiKGRQ8GECu9MRqNRz8MSQJTS\n3fLp169fAGZyccOyLHbs2OFUk7Asi8jISDz99NOqjr/qqquCSo0VUGEeilDFR/d7x97r9/MdPnvY\ncycApy9whYwdAs/Y2lvSh21hnb7cNTVc35ISgGkb4n7wEG2qDTmYcAaxPTVkSIgUZqMMbenrNO4q\nlbqTY9rQaUjtm4oQhCC1byrWXLMGh5YckhSYeHbas+rnpuMznqrFB0sQSzBjtVphNpuRnZ2NjIwM\nWK1WWCwWzwfyUFuqL1AEVJhXL6lGfaN8xFqLXV2tPS3MSHBf5d5BZEQkLBagtt15o+mMNOCmjbZh\nY/lGiQH0FPHwhZ7v5bPOHACOnvNgaHXgWMDzVuKthtOCQCe10d1MOIOCuwqwfdF2FNxV4DRQmwYI\ni2WkDQyOx8xLAavV6tbVMDIyEjk5OUGzWgwG5FbgY8eORU1NDWw2Gw4ePIjLL78c4eHaylRWVlYG\njb4cCLAwj4mMQZ/uUg+TgT0GSpJJ+YM5SXPQp5v7OpgAYGWtMJuBIUM4Id0zUl4ncd+G+8BSK36/\ndAeSzCynCw93b/nu1ZaBzz/3TWcOADZoM0bytTLHz9U5o1vb2rgnC7U4fM75nkb8KNFRA0Z1yHfX\n0ajJwU8IWUkIKSeE7CGEuDeOBACr1YrExEQsW7ZMsU9DQwP27/cqZq9LwrIssrOzMWXKFGRlZeHd\nd9/F7NmzceyYMBFddXU1PvroI83jB1Pem4DrzCfESt3+/jj+j353SwQ4QbT5lvboSQrFgsXDooaB\nYYDHHwdGjQL+34J5sv1ON5/G2NXj8FRNNsqyMzHzahbH3CTZAoAzn72IUaO46kS+0CPE+9D69AFm\nQaCTr78/fpRo3h15HfLddSRqcvATQmYDSKSUJgG4B8C/Aj5RERs2bEBTk5toYx0JFosFRUVFaGtr\ng8ViwR/+8AccPHgQJpMJBoPBmeY2JSUFISHaxCGlNKhunAEX5gvTF0raRg0Y1WHnO9XUngqXQNbH\nHACa7c1gWeAf/wD27gXWvpAs3xHAscajgMEGRO9HccS/kBIxSfnklZcBJzLQ1MQVsvCFByY84PWx\nY3tc4/XKXAm5FftFhJoc/NcA+AAAKKUFAKIIIZ0aEcVPLOaOIUM82HG6OA61itVqxb59+0AIcbof\nA8CpU6ewcuVKbNu2DQcPHkRUVBSOHj2Kl19+WfO51KRUCBQBF+Zzk+ciqbcrampE3xGYOnRqh51P\nTQZGywkLLBbg8GHOsFmzIxMDI4Z6PI7O+D+sP/qWcoceXPENQoAp6lzeFZmeON1zJz68p5D0oUP8\nujLvAqjJwS/uUyvTJ6CoDcv//vvvO3gmwYujjN7kyZORkJCAe++9FwMHDsRXX30Fs9kMo9EIk8mE\n8ePHY8KECejTpw/sdrvH7JNKaAkw6mgCm88c3Iqu8J5C/FLLJbwaHzu+Q1d3jS2ehXlMzxjEx3PZ\nEO12wNDG4NvrtmPix8PQAjeGWU/egj2PA+BuEIcPu+pzekNcpPcSuLG+t2RlrmdE9S98PbY4Pa+/\nULsKvJR9zC0WCyoqKmC329Hc3AwAOHHiBKKjo5Gfn4/i4mKkpaU5DcRaCzfziYuL6xCvIW9z9Qdc\nmAOcQJ+eoHGl6SU1DZ5XM/+c80+UlHCCDgCam4HTNTHIHJKJvMN5PpzdZbSsk69drRpf6nMWtn6E\n2NgMHDkCjBjhu2dNF0BNDv5aAEM89HHizijpL9SuAi9lH3Oz2Qyz2Yzi4mKEhobCZrM589M4Clrz\n0VK4WczixYs7xGvI21z9QRMB2lEorWhDEYoZw2Yg97ZcZA/NhnjR09gItNrUlY1SxOZydTrp3k7q\nkf49NK62eE8N0d36g1J4FTjURVGTg389gN8DACFkAoAzlNLjgZ2mkMzMTI9FJmJjYy9pH3OGYZCX\nl4e8vDxUVFQgNzfXbWCPJ599JcLDw3HzzTf7MlW/0+WFeb8e8quUcYPHYfPvN7stcSbIh+4NB651\nvr3qKt+GqjpdJWlbNHIRPr3+U4/H/mff105vGi1BQ10VpRz8hJB7CCF3t/fZCOAQIeQggNUA7u+0\nCfPwpNv98ccfL3kfc8cKPCYmxmN+mvr6eoFxVA0zZsxAZWVl0FVv6hQ1SyAx9TPJti+bskywLX6C\n7d6d83JRxGN0PQFyXcUdfNWZf2b5TNJWUl+Ct699G6caT+H+TTxZI5rb0eYKpz0gJEQ3gAJcDn4A\nKaK21aLtxQGdlAdycnKcemAlTp2SKWSuo4jZbEZycrImF8O//e1vQSfIgUtgZf7+3vdl28URpyaT\nq+hEWBhX39MtblQWA3sMxDtjSoFz/vvCJ8VJXSDb7JyS/5ZRt6BPhHJwVGL4BIE9oLDQb9PSCSCe\nKsATQoIm6dPFxEmNOtCffvqpg2biG11emL/565uuDd6K1VEGzkFJCZweHy0twP79QO8IaY4WT9yQ\negPKHijDzHHJzuLK4eEqbg4eKDpRJGlzPB4y4QxuSLtB8djLm4X+s6WlCh11ghar1eoxwdbYsWMv\neRWLViwWi2a3xOPHO9V0okiXF+YGyOvDfjj0g2Bb7G1SVwekD0hXHlhBzXLLqFvAhDOornYZHG02\n3wN1FpqlwVYvTH/B+f5A3QHZ424feTuy04T6nREjZLvqBDEbNmxAW5s0myeflJQUt/t1pJjNZphM\n8qpYJR54wPsAvo6kywvzCzZpCTgAKKsrE2yLhW1NDfDSVS9pOlf30O7OACiH3zqgLbmVEreMvAWJ\nvRMBAN0M3fDdTd8JjLdNrTJh3hT4tXYXpk0DBrTHLiYnAx3gAq3TwajxW7/6ar3Kk1YYhsGXX6qv\nqMUwDAYNGtSBM/KeLi/MlSJAxRZs8Wp1xAggIyYDhXcVKuZgF7PQvNAZAFVS4qoV2tzMqW18gQln\nsPue3dixaAeOP3IcM5NmCvb36i7N9AgA1jPcI6HBwEWihoX5Ng+dzkGNKmD27NkBmEnXgmVZzJkz\nR1P/YMqUyKfLC/MxA8fIto+PFfriTpsGOGItkpJcq9eMmAx0M6oraWendm+nqQp3+VD+MU2+tNzw\nfkNhsQDHjnFqH9018eIkPj7eY5+jR1WmSdZxUlBQoClwKDExMWiNzF1emGfHy/uRi8u9sSxwmqtR\ngZoaYWUgY4hR1bn4qg6+Gs4fBlBPZMRk4O7RdwsbCTCif7JA5aO7Jl6clJSUeOwzderUoCkufDHA\nsiyWLFkCm019eul77703aI3MXV6YKxkxxWXYNmxwebM0NwuzHCb2TVR1rm2HtznfV1e72v1hAFVD\nZHiksIECPx/+RZKqIIiydur4kePHjwetCiAYsVgsqm6SfILZyNzlhXnf7n1l2080nhBsz5vHZRQE\nuJU0X402sOdAVecKD3GF7/flnbatDejjuUaGzxTXS/+Rjzf4mEdAJyjIzMz02MdgMARVsYRgx2w2\no1cveVuTHIQQjB3r/8Lz/qLLC/PMWPl/gl4Rwi+RYYDI9oXt4MHCykCPZD2i6lzNNld03vuiWKUP\nP1Q1hE/ERUn/kQdHDdIeEKUTdKhRn9hsNtVpcnW4a6olYpZSGtTXt8sLcyacQYjMx7zQJnRZ3LgR\ncHyvFRXApk2ufcn9knHgjwewZMISRCBC8Vxnm84634tVGYFQbVxoFblhEgCEYudOYUCUHgF68bF8\n+XKPfRzZAXXU8corr3juxKN79+5BfX27vDAH5CM5exp7Cra3bIHb7eR+yVgxcwXM/ZWrIt019i7n\n+yVLhPseekjdXH0hKjxK0hYWGgaxsZ6vz9cJfqxWK1avXu22z6xZs5Cfnx+0xrlgRGvO8C+//DKo\nr69PwpwQsoAQYiGE2AghGf6alL8ZFjVM0vZV6VeCbXGlLSXV48lz8o+7g3oOwrJpy5zb4nxILW5q\nXPiLilMVkrb0AemYPFnYlpXV8XPR8R8bNmzw2MdT3hYdIVarFXv27FHdPzo6GjNnzvTcsRPxdWVe\nBOBaAMGZeaYdU39puG58lNBvd56ohrNSHAExSCMtk3on4cDiAwL/b7n0AB1NaZ0o6QoFjrP1kqeM\nH4SZDHSCGJZlERkZ6bHfqVOndE8WDaxbt85jegQ+06cHppiOL/gkzCmlByil5VCRELYzeWjs4/+/\nvXMPk6I68//37csMDDSXAebCDM4wzAzDTCNgNg4gsPxYJASzWVR2XdFE+WliEmKia1QSsuqP9VHU\nTfJgwAtGMU8CotFdGbna86g/gSHihRF7gAFEhPQIGwhYBRoj8O4f1dVT3V1VXTNTVdOX83meerqr\n6vR5z+k5863T73nPOUmrHHopfgaoVcGT/nY66VpfX9+kiTx6ywM4Tb43P+naoAvVkKT4a59+mpRM\nkIbIsowpU6bg2muvTZm2rKwsrf256YQsy7jtttu69Jmurt/SG+SEz/yLjlpg3xVx14oD8Rut/+lP\n8Z+JGGwQxjpr344tSo5lT3TTuBExNqhvcphVsKI8aQp/H+MxXEEa0ZU46MWLF6e1PzedeP3118Fd\n3HYrE0I+U25OQUQhAFrlIyj93MXM/EpXjLmx6a0ewSBQfOJ6HOcNsd8Q3xwdvyhR4iBhorirzBsz\nD0+1PhV37Z7p9ySlS9yG0Y1tGT/+NGFkk4Cb/u46vBkfUo+aGufL4jTd3fQ2kwgGg2hoaMD777+f\nMu2uXbtcKFHmI8syFixY0K3PpTspxZyZL7fLmBub3uoRCABP/uQKXLWxFhhyENWF1ZhT0+kUl2Ug\nHI7/jFH0ybEzJ5RHWfSRNqNiFmqHJm8h1NioLGzl8SgrFbqxLeMjMx/BdS9fF3Mplby7HAEarruL\nUqbT3U1vMwl1P8uf/exnWL58uWnaEYkj+DmILMsIh8MIBoOGv1LC4XC3dmOa1dN9H13ATjdL2vrN\nZRm468cBXHjyHVS9sR2vz38nzscdDgMffRT/GaPok7F58T36yf3nG9p1exPl+ePm456G1cDJkcDz\nq3Fyy0K0tSW7jIxcSIL0ZMOGDab3PR4PbrrpJpdKk550dHRgwoQJmDZtGqZOnWrYk87PTx5XssLR\no0d7UjxX6Glo4lwiOgpgIoD1RLQp1Wd6g3AYOHQIwN8C+Hj7RBw5EP/UThwgBIADB/TzWjjjavil\nGuCCB3lnavD9GVfppnv9deX1wgVlwpBbkWM/+dp8eB8/BP+H81FfDzQ0AIneiDTd9UqgQzgcTjnr\n8I477kjLPSndQpZlTJ48GR9++CHOnTuHcDisG6opyzIuv9w2R0P6wcyuHIqp3kGSmEePVvrJ48Yp\n51pqatQ+dOcxZ45xfs1bJS796g6OnJB070sSc2lpfH5r19pYoRQUFjJv2NBZz9Wr48uyerV7ZXGL\naPtyrT1rDyfbtiRJPHz4cIbiPNM9mpqaHLOfCbS0tLDH44n7ToLBIEuaf3RJknjJkiWm36PRUVtb\nG5eX21ht2zkRzRIIKGujVFQAW7fGr7uyf79+L7xfP+P8igcFMPjsRAwfou+X27ABSFxauot7xvYI\nnw/4ylc665k4eXDlSvfKIugZgUAADz/8sOH9fv36uRZIkK4Eg0HU1saPW7W3t8fi7mVZxsUXX4x7\n7kkOVEjFlVdeiXfeeScjIoVyQswBRdj69IkXcgAwGjczm6fx+efA2bPxa55r2ZTgbCICrtL3xjgC\ns+LWUcv3YcLE0EOH3CuLoGfIsozrr7/e8P7AgQMzQmicJBAIYPPmzXHXtJtIbNiwoUsbUGiZNGlS\nxny/OSPmRJ3buGkxWgDr7bf1r8sy8K1vKeubTJ2qL+iJE8vmzgXccmmqm2xcdZVSvo6O5DGBK67Q\n/2wuQESDiehVImonoi1ElLygjZLuMBG9T0S7iKjX5sq/9dZbpvez2gdsAVmWsWPHDhxI+Hk9b948\nBAIBdHR0YP584yCFVFx33XU9LaJr5IyYezz6kSXFxcnXAKXnrUc43NnT3bNHfwu2ROGurrZezp4S\nDisPk3PnlPKtXp38wElcBCzHWASgmZlHA3gNwE8N0l0AMJ2ZJzCzC4Gl+qTa+/OFF17IiBhoJ1Bn\nyE6bNg233HILfD4ffD4fSkpK0NraClmW8eKLL3Z5gpDKgw8+mFEDyzkj5mfPAp99lixsZBBQabQY\nVTAIjIpuPKRGiySS6KJx81daMKj4zH0+pXwffJCcJsfnl/wTAHW1+d8CmGuQjpAG/x8f6P0BNXz+\n+ec5uyZLOBxGOBzGuXPncOjQIZw7dw4TJ06Ez+fD+vXr0djYiBUrVnQr77KyMixcuNDmEjtLrzdW\nN5Bl4NprFZdDomtknMGKtkY+80AAeO45ZQOLxMFUlcQZlrXJc4ocIxAAioqAF19UyqeGSGp5/nn3\nypOGFDHzcQBg5mMAigzSMYAQEb1NRN8xSOM4hywMcGTCVHMnCAaDKCqK//Nt27YNf4pO3967dy/2\n79/frbzvvvvujPGVq6ScAZoNhMPAwYPKe9U1MnGiEsmydKn+Z370I+P8AgFlxx6jv/Xnn5ufOw1R\n52YUxcXJSxMkDohmGyZLUPxcJ7nRb/DLmPkTIhoGRdT3MvM2g7SOLFUhyzKOHTuWMt2RI0cyyh1g\nF4FAANdee22XN5mwwtVXX217nlbp9lIVVuIX7TjQy3HmY8Ykx5n/8z8nx5cDzPffb57fwYPMI0ca\n37/55vj8vvtd++qSCkli9vuZfT6lrvPmJddv82b3yuMWsBiLC2AvgOLo+xIAey185l4A/2Zy3/b6\nSJLE1dXVKWOgKysrezUGurdZvHhxt2LHzY4ZM2b0drXisNq2c8LNEggobodhw+JdI++9p59+8WLz\n/M6eVXrbRuNOaq/Y6NxJwmHgyy87B0Avvjj+/qJFQJqvse80TQBujL6/AcC6xAREVEBE/aPv+wGY\nBSCcmM5J3nrrLRxUf04aUFpaiu3bt2ecO8BO8hKXBLWBadOm2Z6nG+SEmAOK6yFxUFsvkqVvX/N8\nZBm45hrg2DHj0MTERe4sLHpnG8Eg4Pd3DoAmhiF2Y8G4bOMhAJcTUTuAfwCwFACIqJSI1C19igFs\nI6JdAP4I4BVmfrVXSmvAlVdeifb29px0r6jIsozf//73tue7b9++1InSkJwQc1kG5s1TdvvRCrBR\nWKIZev73RBLj2d3smQcCQEkJcN99yibV2xK8vFu2uFeWdISZ/8LMM5l5NDPPYubT0eufMPM3ou8/\nYubxrIQljmVmg5EV52hsbMRQk3WTr7nmmpzukQPKr5fE+HI7+MEPfmB7nm6QE2IeDiuDnUC8AL/z\nTnJavYlFWoLBzmgVo9DEU6fiz7ux4ma3kWXg+HFFzOfMAc6cib/fhZ2yBL2ILMs4b9IL+Otfk7cv\nFPQcn88Hv9/f28XoFjkh5sFgZ3igKsAdHYDeqpYzZ5rnFQgAL72U7H/X4o3fkQ4+F2OGwmFl+V7V\nZ/7yy/H3c71nngnIsozx48fjVGKvQEO2b8xhhcbGRlRVVdmaZ11dXcZuv5cTYh4IAGvWAP37K66H\nQAAwcrVZ2erPaDapSmJUmpvrIAWDStik6jMflLCTXOKDRpB+hMPhlDM/N27cmLMzP7V4ow168ODB\ntuT3wAMPZKz7KifEXJaB+fMVl8OcOcq50QYNn32WOq+rrkr2v2spLIw/t6mdWSIQUCY0/fu/Kw+u\nG2+Mv59BS03kLBUVFfB4zP81T548mbMzP1XC4TA+iu4qc/p08kbr3aEgg7fhygkx1/OZJ+4spGI2\nWcgor0QSHxSJy+E6ifqg+o//UB5czz4bfz/HZ39mBHv27MG5FIMbpaWlGesOsItgMIj66E9p7ub6\nK1oqKytxqRv7OzpEToh5MAiMHq28V33mFRXJ6fr1Sz31Xs//nkji+Imb4ynhMPDFF4rPvK0teVch\n4WbJDtatW5ex7gC7CAQC2LRpU8pfMVZZunRpRn+nOSHmgQDwwgvKeuaqz3z9+uR0VkQ3EABeeUXx\nRRsNgCaG/paWdq/c3UH1mXu9ShkT9zK1oQMjcJjGxkZUVlaapol0YyNXdbnYbPK179mzBxdShaBZ\nxCwUNBPICTFXJ/r89a+dPvPy8uR0Vvd6DQQUsTR6iDc1mZ87jbo2y4kTyfdMAiQEaUIgEMD27dvh\nMwmDam1tBdAp0B0dHTGh1hNtWZZx6aWXptzwONM4odfIu0FdXV1Gu1iAHBHzcBhQJ3Wpfm69tjxp\nkrX8zp5VerzaPPbvB+6+W3kdMSI+/ciR3St3d1DdLEakGuAVpAfDhw/HaNU3qMMzzzyDjo4OjB8/\nHlOmTEF5eTkmT56M4cOHY+TIkZg6dSq++tWv4ujRo3j22Wcxe/Zs7Nu3D+fOnUNbW1tWDJ7Ksoy7\n7rrLlryyIm7fygIudhzo5YW2gsHOhbYiEeby8uQFqNrbreXV0BC/aFd7u/6CXerxk584X0eVSMS8\nLBUV7pXFTZBFGzpLksTr1q0zXQzK4/HwLbfconudiJKul5SUcH19PXu9XvZ6vfzUU0/ZWubeIBQK\nJW3k3N3D5/Pxjh07ertKulht2xnb4LvKgQPKaoLvvss8YECyyBUWWsunpUVZkRBQ8tuxg3nKlPQR\n85YWIeZuH3a2bUmSuKGhIaX4eDweHjhwINfU1LDf7+c+ffqwz+fjYDDIwWCQ/X4/V1dXs9frZQDs\n9/u5ubmZd+zYwS0tLTxq1Cj+4Q9/yCdOnOCWlpaMW3lRkiQOBoO2rZQYDAbT9jsQYq5BkpjHjjUX\nOavF0+vlp8p34UJn66clVXkmT3avLG6SLWIeCoUsC1BRURFHIhHesWNH7FWSJJYkKXZt3Lhx7Pf7\nedy4cXFiderUKZ4zZw4XFBSwz+dLum+GJEm9+gCIRCJ822232dYrr6qq4kgk0it1sYIQcw3a3rTR\n4fFYz+/oUea+fZlDIeZnnkkt5lbcN3aRqmd+553ulcVNclHMrbgGVGHXE95t27bFBNHj8fDLL7+c\nJNKJwn3q1Cmuq6uLPQAikQiHQiEOhUKxB4mdQq/mF4lEuKWlhdvb2zk/P9+2Hrn6qyVdXSzM+UXz\nYwAAEe5JREFUQszjkCTm+npzkTPbbCKRffuUz3i9zAMHmuc7YIBz9dJDkszL09TkbnncIlvEXJIk\nDgQCKQXI6/V2qTdtZGvcuHHs8/l40KBB7PF42OPxcElJCS9fvpyff/75mKtm6NChPGnSJO7bt29c\nOfLy8mLvhwwZwiUlJV3u6ZuVb+zYsez1ejk/P5+9Xi+XlpbaKuR2fI9OI8Rcg3bQ0uiYNMl6XiNG\nmOelPYqLna2bXvmMylJQ0LnLUraRLWLOzDx48GBTAVq0aJFhb7urqD33UCjEPp8v1kufM2cOT5o0\nKWaTiPj222/ntra22ANg2LBhpiLZ3Nzco7K1tLTEfP5OHJWVldzc3JzWQs4sxDyOlhalF20mum++\naS2vUMi6kKuHm5i5Wd59192yuEm2iHkkEkkpQoWFhbYLkNpL1/rX1UFGbc94xIgRfP/998eEXCu2\nFRUV3KdPn9i5dlCxO+6XSCTCBQUF3RLqwsLCpGuPP/54rD7p7ifXIsRcQySiRJ7YIbjpLuZGPfNU\n+5pmOtki5k8++WRKoSIiXrlypSOCru3xS5LE9fX1caGOHo+HL7vsspiv3e/386OPPhrr4YZCoVh6\nn8/HK1eu5EgkwhdffLGp+0Ur9mo+VqJ6zH4ZVFZWxs5Hjx4dNzic7r1xLa6IOYCHoWyQ2wrgJQAD\nTNK6UG19nnvOPsGVJOb8fOtC7vM5Vy+j8umVI43Hd2whW8Q8EolYGuDT+nqdii5pbm7WtWkWJaP2\n8LXir3XH+Hy+JPeL1ndfU1PDtbW1PY5UCQaD3N7ezoMHD2YiSuvQw1S4JeYzAXii75cCeNAkrfO1\nNuC228wFNy/Pel6SxNynj3Uxr6lxrl566P1y8Hqz11euki1iLkmSZdeC1+vlRx99lGtqamwbdFT5\n8ssv+Wtf+xoXFhay3+/nYDAY51/W68Vre9YrVqxIKqvWHbNu3bpYhIrWX69Xx/z8fCai2GBrVVUV\njxw5Uje93+/ntWvXJk28SveIFTNcd7MAmAvgdyb3na6zIW++aS6406dbzytV6F/isXmzc/XS49vf\nTi7D+PHulqE3yBYxTxWa+MADD8Sda10gqmD1tKd+4cIF/v73v8+zZs3ikydP6roltDYOHz7Mo0aN\ninugSJLEAwcO1C2ntpfu9Xo5GAxyQUFB7L066Ul9gEQiEW5ubuZf/vKXPGTIEI5EIvzcc88Zfke/\n/vWvORQK8YABA+J66qJnbr1BNwGYb3Lf6TobsmyZueAuWWI9LyuThLSH250Bjye5DEOGuFuG3iAb\nxFySJB4zZoyhSA0cONB0CntJSQk/8cQTPGTIEPZ4PFxXV8ehUCjWA7YiZpIk8cKFC7mhoYE//fRT\n3TThcJhLSkqYiNjn88WFKxIR/+pXv+L58+dzWVmZJXeJ6orZsmVLkl87EonwsmXLuKysLJa+rq4u\n7kGR2JMfNmxY3C8BPddOJmGbmAMIAditOT6Ivv6jJs1iAC+lyMelqiezYoW54HZlUk+qnrnHo8SW\n+3yda7e4iV6ZRo92twy9QTaIeUtLi6G7AQDfe++9LEkSV1dXM4DYFP76+vpY/LXqktDrAadyw0iS\nxBdddBED4DFjxsS5UJqbm/mxxx7j6dOn84ABA2Ii7fP5OBQKxXzoZWVl7PP5mIh45MiRSaGFiefa\nSUuJ5TMaP9Dr5RsdmRBHngqrbTvlVsPMfLnZfSK6EcAcADNS5XXffffF3k+fPh3TXdoc88AB43s+\nH/CXv1jPK9Xu9tdfDyxfrqzM2NBgvEyum8xI+ZfJPN54442s29Q4GAzioosuwqFDh3Tvh8NhBAIB\nvPrqqxg/fjz27t2LI0eO4MyZM/j6178OALhw4QJGjRqFw4cPx3YrUl93796N3/zmN7jpppvQ1taG\nYDCIQCAAWZbR2tqK5cuX48iRIwCAgwcPYs2aNRg0aBC+973v4fTp0xgwYABWrFiBP/zhD5g5cyb2\n7NmD+vp6NDY2YuvWrWhra4uVhZlx9OhRVFRUxNWHmTF06FCcPHkSRIQ777wTjzzyCC5cuIA9e/Zg\n586dKCgoQDAYxPr16/GF2RKgCXg8nqS1zSsqKrBx48aM3nTCMlYU3+gAMBtAG4AhFtI6+/gywWhV\nQ7+/673nf/kX8555b8+w1Iun74obKVOB1Z+iwDwAYQDnAVxikm42gH0A9gO4O0WettXDbLXExx57\njJkVn3b//v158+bNMVdEdXV1LLpE9TGrvme1Bz9ixAiuq6vj/Px89ng8XF5ezg899FDMhdGvXz+u\nq6tjv98f+4WQl5cX6wlrBxGNQvwS49UjkQg3NTVxVVVV3KSk6upqvuOOO+JmkA4bNozz8vLY4/Fw\nfX09r1q1KslNc/PNN+t+N5WVlbH6anvzmTzwqWK5bVtJZPhh4ACAjwG8Fz0eM0nrQrWNSRQ4IsWf\n3dVfX6tWGQs5Ue9HjfznfyaXa+3a3i2TG3RBzEcDqAHwmpGYQ1nn/yCACgB+KKG3dSZ52lYPSZIM\nZ1aqfl9JkrhPnz7s8Xi4qKiIi4qK2Ov1cnV1ddxEGO2CW6rwbt++Pebq8Hg8PGXKlLiY8ZUrV3JJ\nSUmcGGofFFb97olCL0kSr1y5MmabiHQn9qQ6vvOd78QmMOXn57PP54vVW3UHaR8QmTzwqeKKmHfl\n6G0xT1y/vLy8e/mYTZf/xjfsLXN3+e5348uVwWM/lrHa4LmzPb5uIuYTAWzSnC8y653b3bb1IjVG\njBgREyXtNHdtz9Vo4a3EsMFgMBhbLjcSicTOS0tLubCwkJ944omk3rUdE220vfbi4uIuC7n2ICKu\nr69Pmo7f0tIS58/P5IFPFSHmCdTWxgtcTwYFS0r0xdzN1RHNUJfp9fmU1wzvmFjCZjG/GsBKzfn1\nAB41ycvWumgHItWjSeO/04qiGsqnbjqxYcOGpLzGjBnDPp+Pq6qq+Omnn+aysjImIi4rK+Nly5Zx\neXk5A8qA6rvRNR+cmimZuDyv0dorfr8/paDruVD0liXIdKy27ZQDoNnCDTcAixfHn3eXYBA4dqzz\n3OsFdu4Eamu7n6edBAJAS0t6DcK6CRGFABRrL0ERgMXM/IoTNu0c3A8EAtixYwemTp2Kjz/+GGPG\njInLLxAIxAYcGxoaAABtbW3485//jAULFmDJkiWor6/HkSNH8Itf/AJ79+4FAHz00UdYtmwZOjo6\nwMz45JNPsGbNmtjm0OfPn8ffojuABwIBTJw4sdt1MKubmu/WrVuxc+dO3HrrrbEy1tbW4vz58zhz\n5gwmTJiALVu2qA/MGPn5+Th//jzq6+tj9Tf6bjJx4LPbg/tWFN+OA73cM0+cGdmTX1833xyf1403\n2ldOQfeA/W6WzZpzV90sKt3pHe/cuTPWqy0oKOAVK1bw2LFj41wmiS6U3u7Jqr5udYKQOruzrq4u\nNqhZVVXFa9eujaXJtPVVeoLVtp0zYm6n6yExOiZd3Cu5TDfF/CsG97zoHADNgzIAOsYkLxdrao42\nVl07IzRx6r3ZeW+SWH51q7t0KFtvYbVtk5LWeYiI3bJlhCzb53rYvx9YtQpYsCB93Cu5DBGBmclC\nurkAfg1gKIDTAFqZ+etEVArgKWb+RjTdbADLoES2PM3MS03y7PW2rSLLMqZOnRqLAd+6dWtGuRoy\nvfxOYLlt55KYC7IXqw3eIdtp1bZlWc5on3Gml99uhJgLcgoh5oJsxWrb9rhRGIFAIBA4ixBzgUAg\nyAKEmAsEAkEWIMRcIBAIsgAh5gKBQJAFCDEXCASCLECIuUAgEGQBQswFAoEgCxBiLhAIBFmAEHOB\nQCDIAoSYCwQCQRYgxFwgEAiyACHmAoFAkAUIMRcIBIIsQIi5QCAQZAFCzAUCgSALEGIuEAgEWYAQ\nc4FAIMgChJgLBAJBFiDEXCAQCLIAIeYCgUCQBfRIzIloCRG9T0S7iGgzEZXYVTCBwAmIaB4RhYno\nPBFdYpLusKZt73SzjAJBd+hpz/xhZh7HzBMAbABwrw1lspU33ngj52znYp27wAcArgTw/1OkuwBg\nOjNPYOZLnS+WOW58r2797bKlLunW1nsk5sx8RnPaD8o/QFqRi8KWi3W2CjO3M/MBAJQiKSGN3JDZ\nJE7ZUpd0a+u+nmZARPcD+DaA0wD+T49LJBCkBwwgRETnAaxk5qd6u0ACgRkpxZyIQgCKtZegNPTF\nzPwKM/8cwM+J6G4AtwK4z4mCCgRWSdVmLWZzGTN/QkTDoIj6XmbeZndZBQK7IGa2JyOiEQA2MvNY\ng/v2GBIIDGDmVK6TGET0OoA7mPk9C2nvBSAz8y8N7ou2LXAUK227R24WIqpm5oPR07kA9vakMAKB\ny+i2SSIqAOBh5jNE1A/ALAD/zygT0bYF6UBPB3iWEtFuImoFMBPAj20ok0DgGEQ0l4iOApgIYD0R\nbYpeLyWi9dFkxQC2EdEuAH8E8Aozv9o7JRYIrGGbm0UgEAgEvYfjoVdENJuI9hHR/uggqZO2yono\nNSJqI6IPiOhH0euDiehVImonoi1ENNAh+x4ieo+Imly2O5CI/kBEe6N1b3TDNhHdHp2As5uIVhNR\nnlN2iehpIjpORLs11wxtEdFPiehA9DuZZUcZTMrmyuQ5Ino4Wp9WInqJiAY4YMPSpKpu5u24Fui1\nEwds6OqMzTbyieitaJv6IDpuYw4zO3ZAeVgcBFABwA+gFUCdg/ZKAIyPvu8PoB1AHYCHANwVvX43\ngKUO2b8dwO8BNEXP3bL7LIAF0fc+AAOdtg1gOIBDAPKi588DuMEpuwCmABgPYLfmmq4tAPUAdkW/\ni8poGyQH211/zftbATzukJ2ZUHz5ALAUwIMO2BgNoAbAawAusTFfV7RAr504YENXZxywUxB99UJx\n911qlt7pnvmlAA4w88fM/CWAtQD+ySljzHyMmVuj789AGZAtj9r8bTTZb6EM1toKEZUDmAPgN5rL\nbtgdAGAqM68CAGY+x8yfumEbSiPrR0Q+AH0BRJyyy0pY4KmEy0a2vglgbfS7OAzgAJS26Ajs0uQ5\nZm5mZjXvP0Jp23bbsDqpqqu4ogUG7cRuG3o6U+aAnc+ib/OhdExMfeJOi3kZgKOa8z/BgUrrQUSV\nUJ7QfwRQzMzHAeUPAaDIAZO/AnAn4r9wN+yOBHCCiFZFXTwro9EYjtpm5g4AvwBwBIqIf8rMzU7b\nTaDIwFZiu4vA4XZHRPcT0REA8wHc46StKP8XwCYX7NhFr2mBk2h05i0H8vZEB+GPAQgx89tm6dNm\nurKdEFF/AC8C+HH0yZn4RLN11JeIrgBwPPq0NuvRODHa7ANwCYAVzHwJgLMAFunYsrvOg6D0rCqg\nuFz6EdF1TttNgWO2iCgUHRtQjw+ir/8IAMz8c2a+CMBqKK4WR+xE0ywG8CUzr3HKhiA1OjpjK8x8\ngZV1r8oBNBJRvVn6Hk/nT0EEwEWa8/LoNceI/uR/EcDvmHld9PJxIipm5uPRwan/sdnsZQC+SURz\noLgbAkT0OwDHHLYLKD2co8z8TvT8JShi7nSdZwI4xMx/AQAi+m8Ak12wq8XIVgTACE26Hrc7Zr7c\nYtI1ADaimzOhU9khohuhuPNmdCd/KzYcwnUtcBIDnXEEZpZImeQ2G8Aeo3RO98zfBlBNRBVElAfg\nXwE0OWzzGQB7mHmZ5loTgBuj728AYOuXz8w/Y+aLmLkKSh1fY+ZvAXjFSbtR28cBHCWi2uilfwDQ\nBofrDMW9MpGI+hARRe3ucdguIf6Xj5GtJgD/Go2uGQmgGoBjy9gSUbXm1HTyXA/tzIbiyvsmM3/h\nhI1Ekzbm5aYWJLYTJ9DTGdsgoqFqdBYR9QVwOYB9ph9yasRXMyI7G8po7wEAixy2dRmA81BGyncB\neC9qvxBAc7QcrwIY5GAZ/h6d0Syu2AUwDso/SyuA/4ISzeK4bShLHu8FsBvKAKTfKbtQerwdAL6A\n8iBZAGCwkS0AP4USPbEXwCyH292L0e+gFcoDpdQhOwcAfBxt1+8BeMwBG3Oh+LY/B/AJgE025u24\nFui1Ewds6OqMzTbGRvNtjbatxak+IyYNCQQCQRaQlQOgAoFAkGsIMRcIBIIsQIi5QCAQZAFCzAUC\ngSALEGIuEAgEWYAQc4FAIMgChJgLBAJBFiDEXCAQCLKA/wVp2YVH10YmoAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11c148198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure()\n",
    "subplot(121)\n",
    "plot(t, y, '.-')\n",
    "subplot(122)\n",
    "plot(y[:, 0], y[:, 1], 'k.-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Function profiling\n",
      "==================\n",
      "  Message: <ipython-input-13-7d8070171b87>:35\n",
      "  Time in 30145 calls to Function.__call__: 9.484848e+00s\n",
      "  Time in Function.fn.__call__: 6.588749e+00s (69.466%)\n",
      "  Time in thunks: 6.155458e+00s (64.898%)\n",
      "  Total compile time: 3.325101e+00s\n",
      "    Number of Apply nodes: 23\n",
      "    Theano Optimizer time: 5.740211e-01s\n",
      "       Theano validate time: 5.067110e-03s\n",
      "    Theano Linker time (includes C, CUDA code generation/compiling): 5.047009e-01s\n",
      "       Import time 6.441069e-02s\n",
      "\n",
      "Time in all call to theano.grad() 0.000000e+00s\n",
      "Time since theano import 25.343s\n",
      "Class\n",
      "---\n",
      "<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>\n",
      "  82.9%    82.9%       5.102s       8.46e-05s     Py   60290       2   theano.scan_module.scan_op.Scan\n",
      "   7.6%    90.5%       0.471s       7.80e-06s     Py   60290       2   theano.tensor.basic.Diag\n",
      "   4.6%    95.1%       0.282s       9.37e-06s     Py   30145       1   theano.tensor.nlinalg.ExtractDiag\n",
      "   1.8%    96.9%       0.110s       5.19e-07s     C   211015       7   theano.tensor.elemwise.Elemwise\n",
      "   1.0%    97.9%       0.062s       2.06e-06s     C    30145       1   theano.tensor.subtensor.IncSubtensor\n",
      "   0.7%    98.6%       0.045s       4.94e-07s     C    90435       3   theano.tensor.elemwise.DimShuffle\n",
      "   0.5%    99.1%       0.030s       3.31e-07s     C    90435       3   theano.compile.ops.Shape_i\n",
      "   0.3%    99.5%       0.020s       6.69e-07s     C    30145       1   theano.tensor.subtensor.Subtensor\n",
      "   0.3%    99.7%       0.016s       5.30e-07s     C    30145       1   theano.tensor.basic.AllocEmpty\n",
      "   0.2%    99.9%       0.011s       3.68e-07s     C    30145       1   theano.tensor.basic.ScalarFromTensor\n",
      "   0.1%   100.0%       0.007s       2.22e-07s     C    30145       1   theano.compile.ops.Rebroadcast\n",
      "   ... (remaining 0 Classes account for   0.00%(0.00s) of the runtime)\n",
      "\n",
      "Ops\n",
      "---\n",
      "<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>\n",
      "  45.0%    45.0%       2.770s       9.19e-05s     Py    30145        1   do_while{cpu,scan_fn}\n",
      "  37.9%    82.9%       2.332s       7.74e-05s     Py    30145        1   do_whileall_inplace,cpu,scan_fn}\n",
      "   7.6%    90.5%       0.471s       7.80e-06s     Py    60290        2   Diag\n",
      "   4.6%    95.1%       0.282s       9.37e-06s     Py    30145        1   ExtractDiag{view=False}\n",
      "   1.0%    96.1%       0.062s       2.06e-06s     C     30145        1   IncSubtensor{InplaceSet;:int64:}\n",
      "   0.5%    96.6%       0.030s       3.31e-07s     C     90435        3   Shape_i{0}\n",
      "   0.3%    96.9%       0.020s       6.69e-07s     C     30145        1   Subtensor{int64}\n",
      "   0.3%    97.3%       0.020s       6.60e-07s     C     30145        1   Elemwise{Composite{Switch(LT((i0 + (i1 - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1))), i3), (i0 + (-i4)), Switch(GE((i0 + (i1 - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1))), (i1 - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1))), (i2 + i4), Switch(LE((i1 - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1)), i3), (i2 + i4), ((i0 + Switch(LT(i2, i4), i2, i4) + i1) - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1)))))}}[(0, 1)]\n",
      "   0.3%    97.6%       0.020s       6.53e-07s     C     30145        1   Elemwise{sub,no_inplace}\n",
      "   0.3%    97.9%       0.017s       5.57e-07s     C     30145        1   InplaceDimShuffle{x,0}\n",
      "   0.3%    98.1%       0.016s       5.31e-07s     C     30145        1   Elemwise{add,no_inplace}\n",
      "   0.3%    98.4%       0.016s       5.30e-07s     C     30145        1   AllocEmpty{dtype='float64'}\n",
      "   0.2%    98.6%       0.015s       5.03e-07s     C     30145        1   Elemwise{Mul}[(0, 1)]\n",
      "   0.2%    98.9%       0.015s       4.89e-07s     C     30145        1   Elemwise{Sub}[(0, 1)]\n",
      "   0.2%    99.1%       0.014s       4.72e-07s     C     30145        1   Elemwise{Cast{int64}}\n",
      "   0.2%    99.3%       0.014s       4.63e-07s     C     30145        1   InplaceDimShuffle{x,x}\n",
      "   0.2%    99.6%       0.014s       4.61e-07s     C     30145        1   InplaceDimShuffle{x}\n",
      "   0.2%    99.7%       0.011s       3.68e-07s     C     30145        1   ScalarFromTensor\n",
      "   0.2%    99.9%       0.010s       3.27e-07s     C     30145        1   Elemwise{Inv}[(0, 0)]\n",
      "   0.1%   100.0%       0.007s       2.22e-07s     C     30145        1   Rebroadcast{0}\n",
      "   ... (remaining 0 Ops account for   0.00%(0.00s) of the runtime)\n",
      "\n",
      "Apply\n",
      "------\n",
      "<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>\n",
      "  45.0%    45.0%       2.770s       9.19e-05s   30145    16   do_while{cpu,scan_fn}(Elemwise{Cast{int64}}.0, IncSubtensor{InplaceSet;:int64:}.0, Elemwise{Mul}[(0, 1)].0, b, Elemwise{Sub}[(0, 1)].0, tol, Elemwise{sub,no_inplace}.0)\n",
      "  37.9%    82.9%       2.332s       7.74e-05s   30145    17   do_whileall_inplace,cpu,scan_fn}(Elemwise{Cast{int64}}.0, IncSubtensor{InplaceSet;:int64:}.0, Elemwise{Mul}[(0, 1)].0, b, Elemwise{Sub}[(0, 1)].0, w, tol)\n",
      "   4.7%    87.6%       0.288s       9.54e-06s   30145     7   Diag(ExtractDiag{view=False}.0)\n",
      "   4.6%    92.2%       0.282s       9.37e-06s   30145     2   ExtractDiag{view=False}(A)\n",
      "   3.0%    95.1%       0.183s       6.07e-06s   30145    13   Diag(Elemwise{Inv}[(0, 0)].0)\n",
      "   1.0%    96.1%       0.062s       2.06e-06s   30145    14   IncSubtensor{InplaceSet;:int64:}(AllocEmpty{dtype='float64'}.0, Rebroadcast{0}.0, Constant{1})\n",
      "   0.3%    96.5%       0.020s       6.69e-07s   30145    22   Subtensor{int64}(do_while{cpu,scan_fn}.0, ScalarFromTensor.0)\n",
      "   0.3%    96.8%       0.020s       6.60e-07s   30145    20   Elemwise{Composite{Switch(LT((i0 + (i1 - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1))), i3), (i0 + (-i4)), Switch(GE((i0 + (i1 - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1))), (i1 - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1))), (i2 + i4), Switch(LE((i1 - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1)), i3), (i2 + i4), ((i0 + Switch(LT(i2, i4), i2, i4) + i1) - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1)))))}}[(0, 1)](TensorConstant{-1}, Shape\n",
      "   0.3%    97.1%       0.020s       6.53e-07s   30145     6   Elemwise{sub,no_inplace}(TensorConstant{(1,) of 1.0}, InplaceDimShuffle{x}.0)\n",
      "   0.3%    97.4%       0.017s       5.57e-07s   30145     3   InplaceDimShuffle{x,0}(b)\n",
      "   0.3%    97.6%       0.016s       5.31e-07s   30145     9   Elemwise{add,no_inplace}(TensorConstant{1}, Elemwise{Cast{int64}}.0)\n",
      "   0.3%    97.9%       0.016s       5.30e-07s   30145    12   AllocEmpty{dtype='float64'}(Elemwise{add,no_inplace}.0, Shape_i{0}.0)\n",
      "   0.2%    98.1%       0.015s       5.03e-07s   30145    15   Elemwise{Mul}[(0, 1)](InplaceDimShuffle{x,x}.0, Diag.0)\n",
      "   0.2%    98.4%       0.015s       4.89e-07s   30145    11   Elemwise{Sub}[(0, 1)](A, Diag.0)\n",
      "   0.2%    98.6%       0.014s       4.72e-07s   30145     5   Elemwise{Cast{int64}}(max_iters)\n",
      "   0.2%    98.8%       0.014s       4.63e-07s   30145     1   InplaceDimShuffle{x,x}(w)\n",
      "   0.2%    99.1%       0.014s       4.61e-07s   30145     0   InplaceDimShuffle{x}(w)\n",
      "   0.2%    99.3%       0.013s       4.30e-07s   30145    19   Shape_i{0}(do_whileall_inplace,cpu,scan_fn}.0)\n",
      "   0.2%    99.5%       0.011s       3.68e-07s   30145    21   ScalarFromTensor(Elemwise{Composite{Switch(LT((i0 + (i1 - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1))), i3), (i0 + (-i4)), Switch(GE((i0 + (i1 - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1))), (i1 - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1))), (i2 + i4), Switch(LE((i1 - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1)), i3), (i2 + i4), ((i0 + Switch(LT(i2, i4), i2, i4) + i1) - Composite{Switch(LT(i0, i1), i0, i1)}(i2, i1)))))}}[(0, 1)].0)\n",
      "   0.2%    99.6%       0.010s       3.34e-07s   30145     4   Shape_i{0}(b)\n",
      "   ... (remaining 3 Apply instances account for 0.38%(0.02s) of the runtime)\n",
      "\n",
      "Here are tips to potentially make your code run faster\n",
      "                 (if you think of new ones, suggest them on the mailing list).\n",
      "                 Test them first, as they are not guaranteed to always provide a speedup.\n",
      "  Sorry, no tip for today.\n"
     ]
    }
   ],
   "source": [
    "_jacobi_tt.profile.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "## Theano implementation\n",
    "\n",
    "Implementing with NumPy is fine, but \n",
    "\n",
    "- we need to derive Jacobian of system by hand\n",
    "- if we wanted an SA of the method, need to handle ourselves\n",
    "\n",
    "But, with Theano, if we can express the whole thing in terms of its tensor ops, we get these things for free."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "button": false,
    "collapsed": true,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import theano.tensor as T\n",
    "import theano\n",
    "import numpy as np\n",
    "\n",
    "theano.config.floatX = 'float32'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def _make_jacobi_result_for_vars(A, b, w, tol, max_iters):\n",
    "    dA = T.diag(A)\n",
    "    LU = A - T.diag(dA)\n",
    "    w_m_invD = w * T.diag(1.0 / dA)\n",
    "\n",
    "    # seqs, outs, args\n",
    "    def scan_fn(i, x, w_m_invD, b, LU, w, tol):\n",
    "        xn = w_m_invD.dot(b - LU.dot(x)) + (1.0 - w)*x\n",
    "        norm_dx = T.sqrt(T.sum(T.sqr(x - xn)))\n",
    "        return xn, {}, theano.scan_module.until(norm_dx < tol)\n",
    "\n",
    "    results, updates = theano.scan(\n",
    "        fn=scan_fn,\n",
    "        sequences=[T.arange(max_iters)],\n",
    "        outputs_info=b,\n",
    "        non_sequences=[w_m_invD, b, LU, w, tol]\n",
    "    )\n",
    "    \n",
    "    return results[-1]\n",
    "\n",
    "def _make_jacobi_theano():\n",
    "    A, b = T.dmatrix('A'), T.dvector('b')\n",
    "    A.tag.test_value = np.zeros((5, 5))\n",
    "    b.tag.test_value = np.zeros((5, ))\n",
    "\n",
    "    w = T.dscalar('w')\n",
    "    w.tag.test_value = 2.0/3.0\n",
    "    tol = T.dscalar('tol')\n",
    "    tol.tag.test_value = 1e-9\n",
    "    max_iters = T.iscalar('max_iters')\n",
    "    max_iters.tag.test_value = 100\n",
    "\n",
    "    result = _make_jacobi_result_for_vars(A, b, w, tol, max_iters)\n",
    "\n",
    "    return theano.function([A, b, w, tol, max_iters], result, profile=True)\n",
    "\n",
    "_jacobi_tt = _make_jacobi_theano()\n",
    "\n",
    "def jacobi_tt(A, b, w=2.0/3.0, tol=1e-6, max_iters=100):\n",
    "    return _jacobi_tt(A, b, w, tol, max_iters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "Maybe more optimization occurs when we build larger and larger graphs? The main interest "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.00000034,  2.00000069, -1.00000024,  0.99999916])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jacobi_tt(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.,  2., -1.,  1.])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "linalg.solve(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "for{cpu,scan_fn}.0"
      ]
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def _():\n",
    "    A, b = T.dmatrix('A'), T.dvector('b')\n",
    "    A.tag.test_value = np.zeros((5, 5))\n",
    "    b.tag.test_value = np.zeros((5, ))\n",
    "\n",
    "    w = T.dscalar('w')\n",
    "    tol = T.dscalar('tol')\n",
    "    max_iters = T.iscalar('max_iters')\n",
    "    \n",
    "    result = _make_jacobi_result_for_vars(A, b**2, w, tol, max_iters)\n",
    "    \n",
    "    gls = T.jacobian(result, b)\n",
    "    \n",
    "    return gls\n",
    "    \n",
    "_()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "button": false,
    "new_sheet": false,
    "run_control": {
     "read_only": false
    }
   },
   "source": [
    "What about using scan to compute ODE solution? SA for free.."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
 }
No results found