notes/adjoint/Adjoint-method.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recurrence relation\n",
    "\n",
    "This is an example from my [notes on adjoint methods for recurrence relations](http://math.mit.edu/~stevenj/18.336/recurrence2.pdf).\n",
    "\n",
    "In particular, we have a simple linear recurrence:\n",
    "$$\n",
    "x^n = Ax^{n-1} + \\begin{pmatrix} 0 \\\\ p_n \\end{pmatrix}\n",
    "$$\n",
    "in terms of parameters $p$, with initial condition:\n",
    "$$\n",
    "x^0 = b = \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix}\n",
    "$$\n",
    "and\n",
    "$$\n",
    "A = \\begin{pmatrix} \\cos\\theta & \\sin\\theta \\\\ -\\sin\\theta & \\cos\\theta \\end{pmatrix}\n",
    "$$\n",
    "We execute this recurrence for $N$ steps, and at the end we compute $g(p) = (x^N_2)^2$.  The goal is to compute the partial derivatives $d g / d p_n$ for $n = 1,\\ldots,N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Int64,1}:\n",
       " 1\n",
       " 0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "θ = 0.1\n",
    "A = [ cos(θ) sin(θ); -sin(θ) cos(θ)]\n",
    "N = 100\n",
    "p = rand(N)\n",
    "b = [1,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "recurrence (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute xᴺ\n",
    "function recurrence(A, b, p)\n",
    "    x = b\n",
    "    for i = 1:length(p)\n",
    "        x = A * x + [0,p[i]]\n",
    "    end\n",
    "    return x\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       "  9.29829\n",
       " -3.75176"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "recurrence(A, b, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the implementation of $g$ and $dg/dp$ by the adjoint method, as explained in the notes.   This involves computing the recurrence \"forward\" to obtain $x^n$, *saving* the results for each $n$ (using a $2\\times N+1$ matrix `X` here), then computing the adjoint recurrence \"backward\" to obtain $\\lambda$ and accumulate $dg/dp$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "recurrence_adjoint (generic function with 1 method)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute g(x) = x[2]^2 from last iteration, and dg/dp, by adjoint method\n",
    "function recurrence_adjoint(A, b, p)\n",
    "    X = Array(Float64, 2, length(p)+1)\n",
    "    X[:,1] = b\n",
    "    for i = 1:length(p)\n",
    "        X[:,i+1] = A * X[:,i] + [0,p[i]]\n",
    "    end\n",
    "    λ = [0, 2*X[2,end]]\n",
    "    dgdp = Array(Float64, length(p))\n",
    "    for i = length(p):-1:1\n",
    "        dgdp[i] = λ[2]\n",
    "        λ = A'λ\n",
    "    end\n",
    "    return X[2,end]^2, dgdp\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(14.07573606941016,[6.67207,6.98148,7.22113,7.38863,7.48231,7.50123,7.44519,7.31477,7.11126,6.83669  …  -4.66427,-5.22776,-5.73902,-6.19293,-6.58497,-6.91121,-7.16839,-7.35396,-7.46604,-7.50353])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "recurrence_adjoint(A, b, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For comparison, we will compute $dg/dp$ by a center difference approximation: $\\partial g/\\partial p_n \\approx [g(p + \\delta e_n) - g(p - \\delta e_n)] / 2\\delta$ for some small $\\delta$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "recurrence_FD (generic function with 2 methods)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function recurrence_FD(A, b, p, δ=1e-5)\n",
    "    g = recurrence(A, b, p)[2]^2\n",
    "    dgdp = Array(Float64, length(p))\n",
    "    for i = 1:length(p)\n",
    "        pold = p[i]\n",
    "        p[i] = pold + δ\n",
    "        g₊ = recurrence(A, b, p)[2]^2\n",
    "        p[i] = pold - δ\n",
    "        g₋ = recurrence(A, b, p)[2]^2\n",
    "        dgdp[i] = (g₊ - g₋) / (2δ)\n",
    "        p[i] = pold\n",
    "    end\n",
    "    return g, dgdp\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we might have hoped, the answers are pretty close (limited by roundoff error and/or finite-difference error):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.5216061289064273e-6"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "norm(recurrence_FD(A, b, p)[2] - recurrence_adjoint(A, b, p)[2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In practice, it is really easy to make mistakes in computing derivatives of complicated functions by hand, so it is always important to check your results against a finite-difference approximation.\n",
    "\n",
    "It is interesting to ask what the influence of the finite-difference step $\\delta$ is in this problem, by plotting the error vs. $\\delta$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtgAAAIzCAYAAAA6UG5EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XtY1HX6//HncBpAQPPABriamdhBV8WvxgZZfTcKD2gHrCXCTFu37WCb22paZn6pa4vV9adr7WYFi1qkmQc8lKtpJzIt1LZ2S6wQClFQQZFg5PT747MMoGA4zDDM8HpcF5cz7/kMc89eu9urt/fnfpvq6urqEBERERERu/BwdgEiIiIiIu5EAVtERERExI4UsEVERERE7EgBW0RERETEjhSwRURERETsSAFbRERERMSOFLBFREREROxIAVtERERExI4UsEVERERE7EgBuwXTpk0jJCSEoKAgLr/8clJTU51dkoiIiIi4AJOOSm/eV199Rb9+/fD19eWzzz7j2muv5auvvuKSSy5xdmkiIiIi0oF5ObuAjuqKK65o8jwoKIiAgAAnVSMiIiIirkItIufxwAMP4O/vz7XXXsvLL79Mz549nV2SiIiIiHRwbhGwy8vLmTdvHrGxsXTv3h0PDw/S09ObvdZisTBr1ixCQ0Px9/cnMjKS7du3N3vtiy++SHl5OatXr2bKlCl8//33jvwaIiIiIuIG3CJgFxcXk5yczIEDBxg6dCgAJpOp2WsnT57MokWLSEpKYsmSJXh6ejJmzBiysrKavd5kMhEXF8c111zDhg0bHPYdRERERMQ9uEUPdmhoKEeOHCE4OJjs7GxGjBjR7HV79uxh1apVLFiwgBkzZgCQlJTEoEGDmDlzZoshG6CqqoouXbo4pH4RERERcR9usYPt4+NDcHAwAOcbirJmzRq8vLyYNm2adc1sNjN16lR27dpFQUEBAKdOneL111+nvLyc6upq3nzzTT755BNiYmIc+0VERERExOW5RcBurX379hEeHn7ONJD6He/9+/cDRlvIK6+8Qu/evQkODmbJkiVs3LiR3r17t3vNIiIiIuJa3KJFpLUKCwsJCQk5Z71+7fDhwwAEBgayY8eOdq1NRERERNxDpwrYFRUVmM3mc9Z9fX2tr9vi2LFjbN26lUsuuQQ/P7821SgiIiIi9ldRUcGhQ4e4+eabHT56uVMFbD8/PywWyznrlZWV1tdtsXXrVu6+++421SYiIiIijrdy5UoSExMd+hmdKmCHhIRY20AaKywsBIxpJLaoPz595cqV55wA2VHFxMSwbds2Z5fRaq5WL7hezarX8VytZtXrWK5WL7hezarX8Vyp5q+++oq7777bmtscqVMF7GHDhvHee+9RVlZGYGCgdX337t0A1hnaF6p+5/uKK64gIiKi7YW2A29vb5epFVyvXnC9mlWv47lazarXsVytXnC9mlWv47lize3RztuppojEx8dTU1PDsmXLrGsWi4W0tDQiIyMJCwtzYnXty9W+q6vVC65Xs+p1PFerWfU6lqvVC65Xs+p1PFesuT24zQ720qVLKS0ttbaAZGZmkp+fD8D06dMJCgpi5MiRTJw4kdmzZ1NUVET//v1JT08nPz+ftLQ0Z5bf7lztfxCuVi+4Xs2q1/FcrWbV61iuVi+4Xs2q1/Fcseb24DYBe+HCheTl5QHGHOt169axdu1aTCYTkyZNIigoCIDly5czd+5cVqxYQUlJCUOGDGHTpk1ER0c7s3wRERERcRNuE7Bzc3NbdZ3ZbCYlJYWUlBQHV9SxJSQkOLuEC+Jq9YLr1ax6Hc/Vala9juVq9YLr1ax6Hc8Va24PprrznS0urbJ3716GDx9Odna2yzX6i4iIiHQG7ZnXOtVNjiIiIiIijqaALSIiIiJiRwrYIiIiIiJ2pIAtIiIiImJHCtgiIiIiInakgC0iIiIiYkcK2CIiIiIidqSALSIiIiJiRwrYIiIiIiJ2pIAtIiIiImJHCtgiIiIiInakgC0iIiIiYkcK2CIiIiIidqSA3YwzZ84wZcoU+vbtS9euXfnlL3/JJ5984uyyRERERMQFKGA3o7q6mn79+pGVlcXJkyf5/e9/T1xcHOXl5c4uTUREREQ6OAXsZvj7+zN37lx69+4NwJ133omPjw85OTlOrkxEREREOjoF7FY4ePAgJ06c4LLLLnN2KSIiIiLSwblFwC4vL2fevHnExsbSvXt3PDw8SE9Pb/Zai8XCrFmzCA0Nxd/fn8jISLZv397i766oqODuu+9mzpw5BAYGOuoriIiIiIibcIuAXVxcTHJyMgcOHGDo0KEAmEymZq+dPHkyixYtIikpiSVLluDp6cmYMWPIyso659qqqiomTpxIeHg4c+fOdeh3EBERERH34BYBOzQ0lCNHjpCbm8uf//znFq/bs2cPq1at4rnnnuP555/nvvvuY8eOHfTt25eZM2c2uba2tpakpCQ8PT1b3A0/29NPt+VbiIiIiIg7cIuA7ePjQ3BwMAB1dXUtXrdmzRq8vLyYNm2adc1sNjN16lR27drFDz/8YF3/7W9/y5EjR1i9ejUeHq37j+nDD2HKFBu/hIiIiIi4BbcI2K21b98+wsPDCQgIaLI+YsQIAD7//HMA8vLyePXVV/n000/p2bMngYGBBAYGNttG0lhpKbz/vmNqFxERERHX4OXsAtpTYWEhISEh56zXrx0+fBiAvn37Ultba9NnFBXBF1/A4MG21ykiIiIirqtT7WBXVFRgNpvPWff19bW+3hYxMdC/P1xzDWRmtulXiYiIiIiL6lQ72H5+flgslnPWKysrra+3RXZ2DF5e3nh6hjFhQhhdu4KfXwKjRyeQmtqmXy0iIiIirZSRkUFGRgYABQUFFBQUUFVV1W6f36kCdkhIiLUNpLHCwkLAmEbSFtu2bSMiIoLaWoiIgM8/h5MnYeNG4+ZHhWwRERERx0tISCAhIaHJ2t69exk+fHi7fH6nahEZNmwYOTk5lJWVNVnfvXs3gHWGdlt5eEDjbpNjx+Cjj+zyq0VERESkg+tUATs+Pp6amhqWLVtmXbNYLKSlpREZGUlYWJjdPisqCnr1anheXW38iIiIiIh7c5sWkaVLl1JaWmptAcnMzCQ/Px+A6dOnExQUxMiRI5k4cSKzZ8+mqKiI/v37k56eTn5+PmlpaXatJzXVaAvJyoKwMGNGdlISrFgBXm7zn7qIiIiInM1tot7ChQvJy8sDjGPS161bx9q1azGZTEyaNImgoCAAli9fzty5c1mxYgUlJSUMGTKETZs2ER0dbfeaGvdcv/UW/PrXxmOFbBERERH35TYxLzc3t1XXmc1mUlJSSElJcXBFTd1+O6xaBXfeCbt3GwE7Olo3PoqIiIi4m07Vg+1st90Go0ZBbi4cPAgbNuhodRERERF3o4Ddzn74oeHxiROwejXs3Al1dc6rSURERETsRwG7nTWeLhIQYLSK/O//wqBB8MILxo2QAwdqZ1tERETEVSlgt7PUVBg3DsLDYeJEKCmBHTvgiivgoYdg5UrIyYFNmxSyRURERFyR29zk6ErOvrHxhhuMn0svNfqzAYqLjRF/IiIiIuJatIPdgVx/PfTs2fD80kudVoqIiIiI2EgBuwNJTYW4OLjsMggJMY5X/+8p7iIiIiLiIhSwO5jUVGOE38GDMGQIxMbC5587uyoRERERaS0F7A6qSxfYvNloE4mJMQ6q0XQRERERkY5PAbsD69oVtm6FM2dg/Xpjusi6dXDrrXDoEBw7BvfcY0wkUfAWERER6Rg0RaSD69kTevSAkyeN56WlRthev77pdYWFUFsL//hHu5coIiIiIo1oB9sFXHddw3SRiy6CsWPhn/80boSsd/o0vP46vPGGEbRFRERExDkUsF1A/XSR8HC45RbjEJqYGOMGyPpTIS+6CH72M0hIgIgI43W1joiIiIi0P7WIuIizD6epX5syxTiQJirKeJ6VZdwQWT955Nixlt8vIiIiIvanHewW/O1vfyMiIgIfHx/mz5/v7HJalJoKBw40BOioKAgKani9pAS2bIG6OufUJyIiItLZKGC3IDQ0lPnz53P77bdjMpmcXc4FiY5uaB3x9YWjR43WkdJS59YlIiIi0hkoYLdgwoQJxMXF0a1bN+pcbPs3NRXGjTN6sBMSjBsf334bhg6Fjz92dnUiIiIi7k092G7q7J7rq6+Gu+6CUaOMEyLLyoydbvVmi4iIiNiXW+xgl5eXM2/ePGJjY+nevTseHh6kp6c3e63FYmHWrFmEhobi7+9PZGQk27dvb+eK298ll8AHH8CgQbB3r3EU+6ZNmjIiIiIiYm9uEbCLi4tJTk7mwIEDDB06FKDFvunJkyezaNEikpKSWLJkCZ6enowZM4asrKz2LNkpvLygoqLheXGxMXVEREREROzHLVpEQkNDOXLkCMHBwWRnZzNixIhmr9uzZw+rVq1iwYIFzJgxA4CkpCQGDRrEzJkzWwzZrnaT4/lERcGJE8b4Pk9P+OUvnV2RiIiIiHtxix1sHx8fgoODAc57Q+KaNWvw8vJi2rRp1jWz2czUqVPZtWsXBQUF1vWamhoqKyuprq6mqqqKyspKat3giMT6Q2v69DFOfLz0UmdXJCIiIuJe3CJgt9a+ffsIDw8nICCgyXr9jvf+/futa8nJyfj7+/Pqq6/y7LPP4u/vz8qVK9u1XkdJTYW8PHjiCXjmGfjiC2dXJCIiIuI+OlXALiwsJCQk5Jz1+rXDhw9b155++mlqa2ub/EyaNKndam0PTz4JAwbAvfdCdbWzqxERERFxD50qYFdUVGA2m89Z9/X1tb7emZjNkJYG+/bBggXOrkZERETEPbjFTY6t5efnh8ViOWe9srLS+npbxMTE4O3tTVhYGGFhYQAkJCSQkJDQpt/rSCNHwh/+AE8/DRMmwBVXOLsiERERkbbJyMggIyMDgIKCAgoKCqiqqmq3z+9UATskJKRJG0i9wsJCwJhG0hbbtm0jIiKiTb/DGebPh/Xr4brroFs3HUAjIiIirq25Dc69e/cyfPjwdvn8TtUiMmzYMHJycigrK2uyvnv3bgDrDO3Oxs/P6MUuLjYOoNmwQQfQiIiIiNiqUwXs+Ph4ampqWLZsmXXNYrGQlpZGZGSkta2jM/rmm4bHJ07AG2/AypXw449G2B44UKFbREREpDXcpkVk6dKllJaWWltAMjMzyc/PB2D69OkEBQUxcuRIJk6cyOzZsykqKqJ///6kp6eTn59PWlqaM8t3uqgoKCkxdrEDA8HfH5KSjFDt6QmVlcbrU6aofURERETkfNwmYC9cuJC8vDzAOHlx3bp1rF27FpPJxKRJkwgKCgJg+fLlzJ07lxUrVlBSUsKQIUPYtGkT0dHRzizf6VJTjfCclWWE7dRU+PZb4ybIEyeMa4qLYccO59YpIiIi0tG5TcDOzc1t1XVms5mUlBRSUlIcXJHrOXtnun9/Y7LIxo3G0eoeHsYBNZMnG1NHLrnECUWKiIiIdHCdqgdbLlz90erh4UbLyNKl8M47xvOrrjJCuHqzRURERBq4zQ62OM7ZO9v33GOM9Nu3D+rq4OhR489O3sYuIiIiAmgHW2wQEACnTxuhGqC8HDIyYP9+59YlIiIi0hEoYItNoqKgVy/jcbduxizt//kfmDXLGO0nIiIi0lkpYItNUlNh3DijF/vWW402kf/7P1i8GAYPhthYzc4WERGRzkk92GKzs3uz58yBiRNh1CjYutVY0+xsERER6Wy0gy12NWCAcVBNveJiY7a2iIiISGehgC12Fx3d0J8NTQO3iIiIiLtTwBa7a9yfPXQoZGfDwoXOrkpERESkfagHWxyivue6rg6efBIeewx8fODhh51bl4iIiIijKWCLQ5lM8MwzcOYMTJ9uzMs+ftwY86cbH0VERMQdqUVEHM5kgpQUuPJK2LULcnJgwwaN8BMRERH3pIAt7cJkgqqqhucnThi72S+8YMzQnjJFc7NFRETEPShgS7tpPF0kMBB69IDf/x4uvhhWrjR2tjdtUsgWERER16aA3YLi4mLGjh1LQEAAl19+OTt27HB2SS6v8XSR+Hj44Qc4cgSCgxt2t4uL4cMPnVuniIiISFvoJscWPPjgg4SGhnLs2DG2bdvGHXfcwcGDB7noooucXZpLO/vGxh49YOxYY+e6uNhYO3UKioqM4C0iIiLiarSD3YzTp0+zYcMG5s+fj6+vL3FxcQwePJgNGzY4uzS31Hhne9w4o1975Ej417+cXZmIiIjIhdMOdjMOHjxIQEAAoaGh1rXBgwfz73//24lVubfGO9v5+TBhAlxzDbz2mjFxJCtLo/1ERETENbjFDnZ5eTnz5s0jNjaW7t274+HhQXp6erPXWiwWZs2aRWhoKP7+/kRGRrJ9+/Ym15w+fZqgoKAma0FBQZw+fdph30Ea9OkDH30EsbFwyy2werVugBQRERHX4RYBu7i4mOTkZA4cOMDQoUMBMJlMzV47efJkFi1aRFJSEkuWLMHT05MxY8aQlZVlvSYgIIBTp041ed/JkycJDAx03JeQJrp0MYJ19+5QXm6sFRfD5s1w6JDxXKP9REREpCNyi4AdGhrKkSNHyM3N5c9//nOL1+3Zs4dVq1bx3HPP8fzzz3PfffexY8cO+vbty8yZM63XDRgwgNOnT3P48GHr2hdffMFVV13l0O8hTXl4GK0i9feVensbIbtfP+jaFd54QzvbIiIi0vG4RcD28fEh+L8jJ+rq6lq8bs2aNXh5eTFt2jTrmtlsZurUqezatYsffvgBMHawJ0yYwLx586isrGTjxo18+eWXTJgwwbFfRM6Rmmq0iYSHw913Q2kpvPWW8VpFhfFncbHRoy0iIiLSEXSqmxz37dtHeHg4AQEBTdZHjBgBwOeff07v3r0BePHFF7nnnnvo0aMHP//5z1m9ejXdunVr95rl3Bsbb7vN2LXOzITjx421nj3bvy4RERGR5nSqgF1YWEhISMg56/VrjVtCevbsyebNm9utNrkwqalGW8hHH4GXF3z8MUyfDn/5i/FcRERExFk6VRSpqKjAbDafs+7r62t9XVxH453tv/8dHnoIDhwwjmP/9FON9RMRERHn6FQB28/PD4vFcs56ZWWl9fW2iImJwdvbm7CwMMLCwgBISEggISGhTb9Xftr998OAATBmDNTUGD8lJcYut0K2iIhI55KRkUFGRgYABQUFFBQUUFVV1W6f36kCdkhISJM2kHqFhYUATQ6WscW2bduIiIho0+8Q2/3qVxAa2jDGTzc/ioiIdE7NbXDu3buX4cOHt8vnu8UUkdYaNmwYOTk5lJWVNVnfvXs3gHWGtriuG26AHj0anptMDdNGRERERNpDpwrY8fHx1NTUsGzZMuuaxWIhLS2NyMhIa1uHuK7UVBg/3mgXueYayMuDq6+Gr792dmUiIiLSWbhNi8jSpUspLS21toBkZmaSn58PwPTp0wkKCmLkyJFMnDiR2bNnU1RURP/+/UlPTyc/P5+0tDRnli921Ljn+osv4I47YPhwiIiAo0chOlp92SIiIuI4bhOwFy5cSF5eHmAck75u3TrWrl2LyWRi0qRJBAUFAbB8+XLmzp3LihUrKCkpYciQIWzatIno6Ghnli8OMngwfPYZDB1qjPQD3fwoIiIijuU2ATs3N7dV15nNZlJSUkhJSXFwRdJRdOliHLte79ixhrAtIiIiYm+dqgdbOq+oKGM+dmN1dc6pRURERNybArZ0CqmpMG4chIfDtdfCwYPwhz8oZIuIiIj9uU2LiMhPadxz/cILxsmPQUHw9NNOK0lERETckAK2dEoPPghlZTB7NgQGGrvZIiIiIvaggC2d1uOPGyH7scdg7Vrj5seoKE0XERERkbZRD7Z0as88A1dcAR9/DDk5kJlpjPATERERsZUCtnRqJhNUVzc8P37cCNmnTzuvJhEREXFtCtjS6UVHN4zw8/ODEyfgssvg73+HyZNh4EDtaouIiEjrKWBLp9d4hN+vfw25uXDTTfC738HKlUbryKZNCtkiIiLSOrrJUYRzb2xcvhzefx/y843nxcWwY0f71yUiIiKuRzvYIi341a8aWkc8PSEvDx54wJg2IiIiItISBWyRFjRuHbn7bli4EF5/HQYMgKuvNtbVNiIiIiJnU4uIyHmc3TqSlGQctb5nj/G8pMQI2ZqdLSIiIvW0gy1yAXr1grq6hufHjkFWlvPqERERkY5HAbsFf/vb34iIiMDHx4f58+c7uxzpQKKiGnqzAX72M+fVIiIiIh2PAnYLQkNDmT9/Prfffjsmk8nZ5UgHUt+bPWCAMS87Oxu+/NLZVYmIiEhHoR7sFkyYMAGALVu2UNe4J0CEhp7rH3+EyEiIj4dPP4XAQOfWJSIiIs6nHWyRNvD3hzVr4PBhuO++pv3ZIiIi0jm5ZMAuLy9n3rx5xMbG0r17dzw8PEhPT2/2WovFwqxZswgNDcXf35/IyEi2b9/ezhWLOwsPN3a0V6+GpUudXY2IiIg4m0sG7OLiYpKTkzlw4ABDhw4FaLFPevLkySxatIikpCSWLFmCp6cnY8aMIavR6IfXXnuNwMBAAgMDeeCBB9rlO4h7iY+H3/8eHnkE+vTRfGwREZHOzCV7sENDQzly5AjBwcFkZ2czYsSIZq/bs2cPq1atYsGCBcyYMQOApKQkBg0axMyZM60hOzExkcTExBY/Tzc5SmuUlBgnPn7/PaxaBUePwvr14O3t7MpERESkPbnkDraPjw/BwcEA570Bcc2aNXh5eTFt2jTrmtlsZurUqezatYuCgoIW31tTU0NlZSXV1dVUVVVRWVlJbW2t/b6EuJ1du6C62nj844+wZQv07Am33w4vvwx33gkDB2p3W0RExN25ZMBurX379hEeHk5AQECT9fod7/3797f43uTkZPz9/Xn11Vd59tln8ff3Z+XKlQ6tV1xb4/nYvXpBXBz88Y/GTva0aUaPdk4ObNqkkC0iIuLO3DpgFxYWEhIScs56/drhw4dbfO/TTz9NbW1tk59JkyY5rFZxffXzscPDjT8zM+HJJ+Gjj6B//4briot1+qOIiIg7c8ke7NaqqKjAbDafs+7r62t9XcSe6udjn23UKDh1ygjXAJde2n41iYiISPty64Dt5+eHxWI5Z72ystL6uj3FxMTg7e1NWFgYYWFhACQkJJCQkGDXzxHXk5pqtIV89JHRn52VBfv2wbBhzq5MRETE/WRkZJCRkQFAQUEBBQUFVFVVtdvnu3XADgkJabYNpLCwEDCmkdjTtm3biIiIsOvvFPdRv7t9+jTccAOMHm3cGNmvn3PrEhERcTfNbXDu3buX4cOHt8vnu3UP9rBhw8jJyaGsrKzJ+u7duwGsM7RF2lNAAGzebByrfvPNDW0jIiIi4h7cOmDHx8dTU1PDsmXLrGsWi4W0tDQiIyOtbRwi7S04GLZuhZMn4corYcAATRYRERFxFy7bIrJ06VJKS0utLSCZmZnk5+cDMH36dIKCghg5ciQTJ05k9uzZFBUV0b9/f9LT08nPzyctLc2Z5Ytw6aUQGWlMGzl2zBjnV1cH+q+miIiIa3PZgL1w4ULy8vIA46TFdevWsXbtWkwmE5MmTSIoKAiA5cuXM3fuXFasWEFJSQlDhgxh06ZNREdHO7N8EQC+/rrhcVkZvP46jBljHE7j4dZ/vyQiIuK+XDZg5+bmtuo6s9lMSkoKKSkpDq5I5MJFRRlHrBcXQ7du0KUL3HEHDB5snAL5ww8QHd3y+D8RERHpeLRHJuJEjQ+nufVWI1BnZRktIzt3wsGDsGYNTJ7s7EpFRESktVx2B1vEXZy9O33NNcaEkf9Ok6SsDF57Da6+2gjadh7fLiIiInamHWyRDigqCnr1Mh5fdBH8/Ofw0ENwySUwfDhcdpmmjoiIiHRUCtgiHVDj1pFbboHvvoOcHCNs790L334L69YpZIuIiHREahER6aDObh3p398Y41evtBR27GjfmkREROSnaQdbxIU0bh3x9ISiIjh0yKkliYiIyFkUsEVcSOPWkfh4CAmBm24yDqkRERGRjkEtIiIupnHryHffGXOyR482xvp17eq8ukRERMSgHWwRF3bppbB1K+TmwvjxUFHh7IpEREREAVvExQ0eDJs3w6efwsCBMGCApouIiIg4kwK2iBu45hq49lr4/nv45hvjYJpRo4wTIcEI3AMHKniLiIi0B/Vgi7iJxtNEzpyBDz+En/0MevaE8nLjp6TECNlnjwAUERER+9EOtoibaDzCr1cv+PWv4eWXG8I1QHExZGU5r0YREZHOQDvYIm4iNdXYnc7KMsJ2/S71Rx9BZiYcPw4mk3FgjYiIiDhOmwL2xo0befvttzn037+bvuSSSxg9ejRxcXH2qM1pzpw5w/3338+7775LaWkpV155JYsWLSIyMtLZpYmcV3OtH/XB+4MP4McfjdMf1683jmAXERER+7MpYJeWlnLLLbfwwQcf4OXlRUhICHV1dWzbto2///3vXHvttWzYsIFu3brZu952UV1dTb9+/cjKyqJ3796sWrWKuLg4Dh06RJcuXZxdnsgFqw/eFgskJcHtt8NLL8F99zm3LhEREXdkUw/2I488wkcffURKSgolJSXk5eWRn59PSUkJzz//PFlZWUyfPt3etbYbf39/5s6dS+/evQG488478fHxIScnx8mVibSN2QwZGfDb38JvfgPDhxunQmq6iIiIiP3YtIO9fv16fve73/HYY481WQ8ICOCPf/wj+fn5LF++3C4FdgQHDx7kxIkTXHbZZc4uRaTNPD3hhRdg1y7Yu9dY03QRERER+7FpB9vLy4vLL7+8xdcHDhyIl5fj7p8sLy9n3rx5xMbG0r17dzw8PEhPT2/2WovFwqxZswgNDcXf35/IyEi2b9/e6s+qqKjg7rvvZs6cOQQGBtrrK4g4lclk9GPXO3YM3n/fefWIiIi4E5sCdnx8PG+++SY1NTXnvFZdXc2bb77JxIkT21xcS4qLi0lOTubAgQMMHToUAJPJ1Oy1kydPZtGiRSQlJbFkyRI8PT0ZM2YMWY1mlb322msEBgYSGBjIAw88YF2vqqpi4sSJhIeHM3fuXId9HxFnaDzWz8MDDh+IMiGUAAAgAElEQVSGnTudW5OIiIg7MNXV1dVd6Js+/PBDHnroIcxmM9OmTWPAgAEA5OTksGzZMqqqqli6dCn+/v5N3hcREWGXos+cOUNpaSnBwcFkZ2czYsQI/vGPfzBp0qQm1+3Zs4fIyEgWLFjAjBkzAGNHe9CgQQQHBzcJ2Werra3lrrvuoqKignXr1uHh0fK/i+zdu5fhw4eTnZ1tt+8o0h7qx/pFREBREbz3HjzzDMyaZYRuERERd9Geec2mPo7rrrvO+vizzz5r9ppRo0Y1eW4ymZrd8baFj48PwcHBAJzv3w/WrFmDl5cX06ZNs66ZzWamTp3KnDlzKCgoICwsrNn3/va3v+XIkSNs3br1vOFaxJU17rmuqYGnn4Y5c+DFF40bIkeNUl+2iIjIhbIpYKe6yD9x9+3bR3h4OAEBAU3WR4wYAcD+/fubDdh5eXm8+uqr+Pn50bNnT+v6O++8Q1RUlGOLFnEST09ITobdu2H7dqirM3a16+ogLc3Z1YmIiLgOmwL25MmT7VyGYxQWFhISEnLOev3a4cOHm31f3759qa2tdWhtIh1VXp4RqgHKyuC11yA2FiZOVNuIiIhIa7j1Py4rKiowm83nrPv6+lpfF5GmGt/82K0bBAfDr38Nv/gFvPkm3HsvDByo2dkiIiItadUO9r333tvilI7m1NXVYTKZnN5K4ufnh8ViOWe9srLS+ro9xcTE4O3tTVhYmLX1JCEhgYSEBLt+jogj1R+tnpVlhO3UVGNm9vz5cMcdRitJTY1mZ4uISMeVkZFBRkYGAAUFBRQUFFBVVdVun9+qgL1z585zAnZ5eTnHjh0D4KKLLqKuro7S0lIAevbs2SGOFA8JCWm2DaSwsBCA0NBQu37etm3bNEVE3MLZofmXv4R33oE+feD774214mIjhIuIiHQ0zW1w1k8RaQ+tahE5dOgQubm51p9Nmzbh4+PDnDlzKCoq4vjx45w4cYKjR48ye/ZsfHx82Lx5s6Nr/0nDhg0jJyeHsrKyJuu7d+8GsM7QFpHWufFGqL/v12SC/07oFBERkUZs6sF++OGHiY2N5ZlnnmkyZaNXr148++yzxMbG8vDDD9utSFvFx8dTU1PDsmXLrGsWi4W0tDQiIyNbHNEnIs1LTYW4OLjsMqM3+/334aOPnF2ViIhIx2LTFJHdu3ef96TGYcOGWfteHGXp0qWUlpZaW0AyMzPJz88HYPr06QQFBTFy5EgmTpzI7NmzKSoqon///qSnp5Ofn0+a5o6J2KS+faS83AjbsbGwZYsxM1tERERsDNgXXXQRW7Zs4Xe/+12zr7/99tt069atTYX9lIULF5KXlwcYh9isW7eOtWvXYjKZmDRpEkFBQQAsX76cuXPnsmLFCkpKShgyZAibNm0iOjraofWJuLsuXWDTJiNkjx5tBOzvvmu4MVJERKSzsqlF5P7772fz5s2MHz+ebdu2cejQIQ4dOsQ///lP4uLi2LJlC/fff7+9a20iNzeX2tpaamtrqampoaamxvq4T58+1uvMZjMpKSkcPnyYiooKPvnkE2JiYhxam0hn4e8PGzdC167GTZA5OUbo1gg/ERHpzGzawX7yySexWCykpKSwadOmJq95e3sze/ZsnnzySbsUKCIdm7+/sZtdr7hYfdkiItK52RSwAZKTk5k+fTrbt2+3tmr07duXmJiYJjc+ioj7u/ZaKC2F/07u5NQpY5zfz3/u3LpEREScodUBe8GCBYwdO5YrrrjCutarVy8doiIiTQ6nuewy+OILGDLEWL/lFmdXJyIi0r5a3YP9/PPPc9VVV3HppZfy8MMP8/bbbzd7SqKIdE6pqXDgAGzeDPv3w/XXw623whVXGPOy1ZctIiKdRasD9tGjR/n444+5++67+fjjjxk3bhzdu3cnLi6Ov//979YReSIi3bvDW28ZJ0B+/TV88w1s2KCQLSIinUOrA7aHhweRkZH83//9H9nZ2RQUFPDXv/4Vs9nM448/ziWXXMLgwYN5/PHH+fDDD6mtrXVk3SLSwZlMcPx4w/MTJ2DrVufVIyIi0l5sGtMHcPHFFzNlyhTWrFlDcXEx7777LjfffDOZmZlcd9119OjRg1//+tfWY8lFpPOJioJevYzHPj5w+DA8/jhUVzu3LhEREUeyOWA35u3tzQ033MCCBQv4z3/+w7fffsszzzzDqVOn+OCDD+zxESLiglJTYdw4CA+HxET4859hwQK4+WYoKnJ2dSIiIo5h85i+8+nXrx8PPvggDz74oCN+vYi4kLNPdfyf/4E774SICBg2zDicRqc/ioiIO7EpYM+fPx+TydTi6yaTCV9fX3r37s2oUaMICwuzuUARcS/XXw979xrhuv6cqhMnjBsgFbJFRMQd2BywW8vT05P77ruPF154AQ8Pu3SkiIiLCwuDbt2MUx/BOKBmyxaorQX934SIiLg6m/5R9v333/OLX/yCyZMnk52dTWlpKaWlpXz22Wfcc889DB06lK+//pq9e/eSmJjISy+9xLPPPmvv2kXEhUVHN9wAaTbD0aNGq8innxq72QMHaqyfiIi4JlNdXV3dhb5pwoQJ+Pn58cYbbzT7+p133onFYmH9+vUAjBkzhm+++YacnJy2VdtB7d27l+HDh5OdnU1ERISzyxFxGfWnP0ZFwaRJ8Mgj8K9/ga8vVFYaAXzcOLWOiIhI27VnXrNpB3vnzp1cf/31Lb5+3XXX8e6771qfjx49mry8PFs+SkTcWP3pj6mpRm92djYEBxvhGowWko8+cmqJIiIiF8ymgO3j48Mnn3zS4uu7d+/GbDZbn1dXVxMQEGDLRznFtGnTCAkJISgoiMsvv5xUbZ+JtAsvLxg7Fnr0aFgrLYXPP3deTSIiIhfKpoB91113sXz5cv7whz/w7bffUltbS21tLd988w0zZsxgxYoVJCQkWK/fuXMnV155pd2KdrRHH32U3NxcTp06xcqVK3nwwQc5dOiQs8sS6RRSU2H8eGN29ujRRpvI8OHw2GNGG4l6s0VEpKOzaYrI888/z9GjR1m0aBGLFi2yTgepPx799ttvJyUlBYDKykqGDx9OVFSUnUp2vCuuuKLJ86CgIJfagRdxdY3/0ujMGfjLX+DJJ6Guzpg0UlKisX4iItJx2RSw/fz8WLVqFbNmzeKdd96x9lf37duX2NjYJo3jvr6+zJs3zz7VtqMHHniAf/zjH9TV1bFq1Sp69uzp7JJEOiUfH+N49Zdegvq/SCouhvffd2pZIiIiLWrTxNmIiAjmzJnDSy+9xEsvvcScOXMcfldmeXk58+bNIzY2lu7du+Ph4UF6enqz11osFmbNmkVoaCj+/v5ERkayffv2Vn3Oiy++SHl5OatXr2bKlCl8//339vwaInKBbrihYayfyQTffw+vvGLsaouIiHQkLnekQ3FxMcnJyRw4cIChQ4cCtHiq5OTJk1m0aBFJSUksWbIET09PxowZQ1ZWlvWa1157jcDAQAIDA3nggQeavN9kMhEXF8c111zDhg0bHPelROQnpaYaI/vCwyEhARIT4Te/MYL37berN1tERDqOVrWIeHh4YDKZaDwyu3GorV+vX6urq8NkMlFTU2PPWgEIDQ3lyJEjBAcHk52dzYgRI5q9bs+ePaxatYoFCxYwY8YMAJKSkhg0aBAzZ860huzExEQSExPP+5lVVVV06dLFvl9ERC7Y2T3XiYlw661w+rTxXL3ZIiLSEbQqYD/11FPnrK1bt45///vf3HzzzQwcOBCAAwcOsHXrVgYPHsytt95q30r/y8fHh+DgYADOd0bOmjVr8PLyYtq0adY1s9nM1KlTmTNnDgUFBYSFhZ3zvlOnTrFp0yYmTJiA2Wxm3bp1fPLJJ7z88sv2/zIi0iY33ggXXwzffGM8Ly42Dq4RERFxplYF7KeffrrJ82XLllFUVMSXX37J5Zdf3uS1r776iv/93/8lNDTUbkXaYt++fYSHh58z/aN+x3v//v3NBmyTycQrr7zCgw8+iMlk4qqrrmLjxo307t27XeoWkQtz7bVw8qQRrsE4qEZERMSZbJoikpKSwkMPPXROuAZjxN1DDz1ESkoKv/nNb9pcoK0KCwsJCQk5Z71+7fDhw82+LzAwkB07dji0NhGxn9RUoy2k/sTHTz6Bf/4TbrrJuXWJiEjnZVPALigowNvbu8XXvb29nT51o6KioslpkvV8fX2tr4uIe6jvua6uhgkTjJseP/gAhg1zbl0iItI52RSwBw0axN/+9jfuuuuuc1onvv/+e1588UUGDx5slwJt5efnh8ViOWe9srLS+rq9xcTE4O3tTVhYmLX9JCEhocmpliLiOF5esGqVMVlkzBjYtQsuucTZVYmISHvLyMggIyMDMDaGCwoKqKqqarfPtylgL1q0iJtuuomBAwdyyy23MGDAAABycnJYv349ACtWrLBflTYICQlptg2ksLAQwCE94tu2bXP4HHAROb+AANi0Ca65xjhqPSICPvsMoqI0XUREpLNoboNz7969DB8+vF0+36aAHR0dze7du3nqqadYt25dk13h2NhY5s+f7/Qd7GHDhvHee+9RVlZGYGCgdX337t0A1hnaIuJ+fvYzePtt+MUvjAkj1dUa4SciIu3H5oNmBg8ezLp16ygrK+Pw4cMcPnyYU6dOsXbtWqeHa4D4+HhqampYtmyZdc1isZCWlkZkZGSzE0RExH2EhxsTRaqrjefFxfDee04tSUREOgmbdrAb8/T05OKLL7ZHLa22dOlSSktLrS0gmZmZ5OfnAzB9+nSCgoIYOXIkEydOZPbs2RQVFdG/f3/S09PJz88nLS2tXesVEee48UbIzITjx43j1XNzjcNpZs6ExYuNmdlqHREREXsz1Z3vtJYOql+/fuTl5QHnnh6Zm5tLnz59AGPHeu7cuaxcuZKSkhKGDBlCcnIyMTExdq2nvqcnOztbPdgiHcyUKUaQvvpqGDkSFi6EQ4fA2xuqqqBXL+MIdoVsERH31p55zSUDdkejgC3iOqqr4ec/hyNHGtbCw+HAAefVJCIijteeec3mHmwREVfk5WVMF+nZs2Ht0kudV4+IiLgfBWwR6XRSUyEuDi67DMLCYOdO2LjR2VWJiIi7UMAWkU4pNRUOHoRvvzUOpbntNuOQGhERkbZq0xSR//znP3z33XeUlJTQXCv3pEmT2vLrRUQczmyG1auNmyETEuCvfzVG+mm6iIiI2MqmgP3tt9+SmJjInj17znudAraIuAIvL/jHP2DPHmPiCOhgGhERsZ1NAfu3v/0tX375JYsXLyY6OpqLLrrI3nWJiLQrDw+orW14XlwMH33kvHpERMR12RSws7KymD17Ng8//LC96xERcZroaCgtNcI1wKlTxuNevZxbl4iIuBabbnLs0aMH3bp1s3ctIiJOlZpqHDoTHm7c+FhXByNGwL/+5ezKRETEldgUsH/3u9+xcuVKqqur7V2PiIhTpaYah85s3gyffgoXXQTXXAO/+hUMHGj0ZYuIiJyPTS0i4eHh1NTUMHToUO6991769OmDp6fnOdfddtttbS5QRMRZ+vQx+rAHDYIdO4y1oiLjNMjlyxuuqz+OXZNHREQEbAzYd955p/XxH//4x2avMZlM1NTU2FaViEgH0aULeHs3PC8thRUr4OOPYfBgyM01fk6d0uQREREx2BSwd9Rv5YiIdAKNb37s1s3Y0R45Er74Ar78Eur3EoqLG8b8iYhI52VTwL7++uvtXIaISMeVmtpyG8iUKfDWW8YOttls9GuLiEjn1qaTHEVEOouW2j7q1995BwoLjce1tcZcbRER6ZxaFbBvuOEGTCYT//znP/Hy8rI+b0ldXR0mk8nlW0l27dpFVFQUycnJPPHEE84uR0Q6qPqQ/frrcPfd4OcHL7wA5/m/SRERcWOtCth1dXXNPj973Z3U1tby6KOPcvXVV5/3XyZEROrddRdUVsLUqeDrCwsXKmSLiHRGrQrY77333nmfu6Nly5YRGRnJyZMn3fpfJETEvqZMgYoKeOghePddI3BrfJ+ISOeiLsFmHD9+nMWLFzN//nxnlyIiLujBBxtOgMzJgY0bdUCNiEhn4nIBu7y8nHnz5hEbG0v37t3x8PAgPT292WstFguzZs0iNDQUf39/IiMj2b59+09+xhNPPMGjjz5K165dAdQiIiIX7OTJhsfHjsHWrc6rRURE2pfLBezi4mKSk5M5cOAAQ4cOBVoOwJMnT2bRokUkJSWxZMkSPD09GTNmDFmNBtW+9tprBAYGEhgYyAMPPMC+ffv47LPPuO+++wCjz1wtIiJyoaKioFcv47GPDxw+DBMmQF6ec+sSERHHc7kxfaGhoRw5coTg4GCys7MZMWJEs9ft2bOHVatWsWDBAmbMmAFAUlISgwYNYubMmdaQnZiYSGJiovV9ixcv5sCBA4SFhQFw8uRJvLy8+O6773j11Vcd/O1ExF2cPTt77FiYPh2uvNL4OXnSOMBGvdkiIu7H5QK2j48PwcHBwPmnmKxZswYvLy+mTZtmXTObzUydOpU5c+ZQUFBgDdGNTZs2jYSEBOvvf+SRR7j00kt5/PHH7fxNRMTdnR2eb7oJIiPhs8+M58XFzV8nIiKu7YJbRCwWC5mZmXz++eeOqMdu9u3bR3h4OAEBAU3W63e89+/f3+z7/Pz8CA4OJjg4mJ/97Gf4+fkREBBAUFCQw2sWEfcWGAjV1Q3PS0thzZqGA2pERMQ9XHDA9vb2Jj4+nl27djmiHrspLCwkJCTknPX6tcOHD7fq96SlpTFnzhy71iYinVfj3uzAQDhzBgYOhEWLoKrKubWJiIh9XHDA9vDwYMCAARw7dswR9dhNRUUFZrP5nHVfX1/r6yIi7S01FcaNg/BwiI83bn5MSoI//AGCg+HnP9dIPxERV2dTD/acOXOYMWMG8fHxXH755fauyS78/PywWCznrFdWVlpft7eYmBi8vb0JCwuz9ncnJCRYe7pFRODcnusXXoAffoAtW4y2kdWrm79ORERaJyMjg4yMDAAKCgooKCigqh3/mtCmgP3JJ5/Qs2dPBg8ezHXXXUe/fv2aDaxLlixpc4G2CgkJabYNpPC/zY6hoaF2/8xt27YRERFh998rIu7v668b+rPLy+Gtt+Bvf4Nm/iJORER+QnMbnHv37mX48OHt8vk2BewXXnjB+njHjh0tXufMgD1s2DDee+89ysrKCAwMtK7v3r0bwDpDW0SkI4iKgpISY7JIYCCUlUFMDKxdCz17Ors6ERG5EDYdNFNbW9uqH2eKj4+npqaGZcuWWdcsFgtpaWlERkY2O6JPRMRZzu7N/ugjY1c7MtL4U0REXIfLzcEGWLp0KaWlpdYWkMzMTPLz8wGYPn06QUFBjBw5kokTJzJ79myKioro378/6enp5Ofnk5aW5szyRUSadXbP9e7dRuiOjDR+cnONnW71ZouIdGymujacA/7dd9/x9ttvW8Nt3759GT16NP369bNbgc3p168fef89b7j+mPS6ujpMJhO5ubn06dMHMHas586dy8qVKykpKWHIkCEkJycTExNj13rqe3qys7PVgy0idnXyJFx1FRQUGM979oS4OIVsEZEL1Z55zeYd7BkzZrB48eJzTlP08PDgkUceYeHChW0uriW5ubmtus5sNpOSkkJKSorDahERcaSuXcHfv+H5sWOwc6fz6hERkZ9mUw/2woUL+X//7/9x++23s2vXLkpKSigpKWHXrl3Ex8ezaNEi/vKXv9i7VhGRTik6uuFwGg8PY6Tfq6+C7X//KCIijmTTDvbLL79MXFwcq+uHtf7X1VdfzRtvvEFlZSXLli1jxowZdilSRKQzS001Dp/JyoKRI43RfffdBxs2QJcusHeverNFRDoSmwL2oUOHeOSRR1p8/aabbuLtt9+2uSgREWnq7PA8fjzceSdYLMZOdkmJEcIVskVEnM+mFpFevXqxf//+Fl//17/+Ra/6v88UERG7Gz8eQkMb2kSKi9WbLSLSUdgUsO+44w5eeeUV/vSnP1FeXm5dP336NM899xwvv/wyd955p92KFBGRc113XdPe7Lw8eOopqKhwbl0iIp2dTWP6ysvLGT9+PDt37sTLy4vQ0FDq6uo4fPgwNTU13HDDDWRmZtKlSxdH1NzhaEyfiDhLfW92ZCT07QvPPw9hYdC/P+TnqzdbRKRehx/T16VLF9599102bNjAli1brDOpY2NjGTt2LHFxcdb51CIi4jhnh+ekJLj+eti+3Xh+4oR6s0VE2tsFB+wff/yRxMRE4uPjSUxMZMKECY6oS0REbDBggDFZpN6xY8ax6yIi0n4uuAfb39+fd999lx9//NER9YiISBs1nptdTzOzRUTaj003OUZHR7Nr1y571yIiInaQmgrjxkF4OIwaBQcPwmOPKWSLiLQXmwL20qVL+fDDD3niiSf44Ycf7F2TiIi0UWoqHDgA778Pf/0r/OUv8Oyzzq5KRKRzsOkmxyFDhlBdXc2f/vQn/vSnP+Ht7Y3ZbKaurg6TyWT989SpU/auV0RELtBDD0FpKcydC127wsMPO7siERH3ZlPAvv3223/yGk0RERHpOJ54wgjZ06cbO9ljxmiyiIiIo1xwwK6rq2Px4sX4+Pjg5+fniJpERMTOTCZjZJ+vLxw9ChkZcPo0rF7t7MpERNzPBfdgWywWunfvzl//+ldH1NNhXH/99fj5+REYGEhgYCBjx451dkkiIm2SlQWVlcbjykp4801jR/vIEefWJSLibi44YPv6+nLxxRdjNpsdUU+HYTKZePXVVykrK6OsrIzNmzc7uyQRkTaJimoY39ezJ0REwPLlxqmPjz8OiYkwcKBxMI2IiNjOpiki9957L8uXL8disdi7ng7FhlPkRUQ6rMbj++LiIDsbcnPh97+HhQuNtpGcHNi4USFbRKQtbLrJcfDgwaxfv55BgwZxzz330K9fv2b7sW+77bY2F+hMjz76KI8++ihDhw5l4cKFDB482NkliYi0ydk3Nl50kXHT4xtvwHffGWvHjjUctS4iIhfOpoCdkJBgffzUU081e43JZKKmpsa2qn5CeXk5KSkp7N69mz179lBaWkpaWhr33HPPOddaLBaeeuopVqxYQWlpKb/4xS945plnuPHGG8/7GSkpKVx11VV4eHiwdOlSRo8ezddff01AQIBDvpOIiDNddx2UlUFxMXh5wfffw29+A889Bz16OLs6ERHXYlPA3rFjh73ruCDFxcUkJyfTt29fhg4dynvvvdfiWMDJkyfz1ltv8eijjzJgwADS0tIYM2YMO3fuJCoqCoDXXnuN+++/H4CkpCRefPFFRowYYf0djz32GKmpqXzyySc/GcxFRFxRaqrRFpKVBddcAyNGwOzZsH49XHGFcSNkdLRG+4mItIZNAfv666+3cxkXJjQ0lCNHjhAcHEx2dnaTMNzYnj17WLVqFQsWLGDGjBmAEaAHDRrEzJkzycrKAiAxMZHExMR2q19EpCM6OzzfdptxY+SHHxrPS0qMEK6QLSJyfjbd5FjPYrGwa9cuNmzYQHFxsb1q+kk+Pj4EBwcD578Rcc2aNXh5eTFt2jTrmtlsZurUqezatYuCgoJm33fy5Em2bduGxWLhzJkzLFq0iNLSUq6++mr7fhERkQ7s4ouNdpF6x441hG0REWmZzQF78eLFXHzxxURFRXHbbbfxxRdfAEb7Ro8ePXj11VftVqSt9u3bR3h4+Dl90/U73vv372/2fVVVVcyZM4devXoREhLC5s2b2bJlC4GBgQ6vWUSkI2k82g+goqJhlraIiDTPpoCdlpbGo48+yujRo0lNTW2yi9yrVy9+9atfsWrVKrsVaavCwkJCQkLOWa9fO3z4cLPv69mzJ59++imnTp3i+PHjbN++naFDhzq0VhGRjqjxaL/YWOM0yAkTjKAtIiLNs6kHe+HChYwfP57XX3+dY8eOnfN6REQES5YsaXNxbVVRUdHsgTi+vr7W10VE5Pwa91zv3Aljx8L48bBhA/j7O68uEZGOyqaA/c033zB9+vQWX+/evTvHjx+3uSh78fPza/YwnMr//v1mc7O72yImJgZvb2/CwsIICwsDjJGGjccaioi4shtugC1bjJB92WXQpQtce61ufBSRjiUjI4OMjAwACgoKKCgooKqqqt0+36aA3bVr1/Pe1PjVV19x8cUX21yUvYSEhDTbBlJYWAgY00jsadu2bURERNj1d4qIdDTXX2/MzX77beP5iROaLiIiHUtzG5x79+5l+PDh7fL5NvVgjx07lpdffpmSkpJzXvv3v//Nyy+/zPjx49tcXFsNGzaMnJwcysrKmqzv3r0bQH3VIiI2+vbbhscnThi72ucZ6iQi0qnYFLCTk5Opqalh8ODBzJ07F4D09HQSExMZPnw4vXr1avGEx/YUHx9PTU0Ny5Yts65ZLBbS0tKIjIy0tnGIiMiFaTxdxGyGo0dh9GjjBEgRkc7OphaRsLAwPvvsM5544gneeOMNAFasWEFgYCB33XUXzz33HL0az3VygKVLl1JaWmptAcnMzCQ/Px+A6dOnExQUxMiRI5k4cSKzZ8+mqKiI/v37k56eTn5+PmlpaQ6tT0TEnTU++TEqCuLjYdo0uOoq+Mtf4OOPG15T64iIdDamuvOd1NIKdXV1FBcXU1tbS69evfD09LRXbefVr18/8vLyAKzHpNfV1WEymcjNzaVPnz6AsWM9d+5cVq5cSUlJCUOGDCE5OZmYmBi71VLf05Odna0ebBHptEpL4Q9/MAK1tzdUVRm73OPGKWSLiPO1Z15rc8AWBWwRkcZ694bGB+WGh8OBA86rR0QE2jevtemodBERkbPddBN07248Npng8sudW4+ISHtTwBYREbtKTTVOe7z0UujRA3bsgG3bnF2ViEj7UcAWERG7S001RvkdOgSjRsGYMbBypbOrEhFpH60K2KdOnaK6utrRtYiIiJvp0sU4Uvhu2YUAACAASURBVH3SJEhKghEjjJ7sKVOcXZmIiOO0KmB369aN1atXW5/fe++91sNaREREzsfLC155BYYMgc8+g4MHYf16hWwRcV+tCthmsxmLxWJ9np6ezreNj/ESERE5D5MJKioanpeUwJo1kJvrvJpERBylVQfNDBw4kFdeeYW+ffvStWtXAHJzc9m7d+9536eRdSIiUi8qygjWxcUQEGDMyR44EH7zGzh+HPbt08E0IuIeWjUH+5133uGOO+7g9OnTrf/FJhM1NTVtKs5VaA62iEjrND79celS42fuXDhzxni9Z0+Ii1PIFhH7a8+81qod7NjYWHJzc/n0008pKipi8uTJTJs2jcjISIcWJyIi7uXs4DxzJixbZkwcATh2DLZsgdpa8NCcKxFxUa0K2AA9evQgNjYWgNTUVOLj47nxxhsdVpiIiHQOo0bBqVNG64jZDEePwi9/CYsXg/ZxRMQV2bQ/8N577ylci/x/9u49Lso6////YwbkJKCRkoB5yNQOmqKLkVqWpRlmVoLGmofN1i0rS9sy7Wfqx+y7kYbrmttqgohGmXm2LbWyg5q2HsqsRbdFScTEFEXlzPz+uJYRFBBHhmsGnvfbrRvM+7qGeQ6h18s3r+v9FpEakZAA999vLN/3+9/Dl18a/dm33QbXX29sWKMVR0TEnVR7BvtCp06dIj4+nvXr13Po0CEsFgstW7akf//+jBs3jsDAwJrMKSIiddiFrSPffgu9esHWrWCzGbPbNhskJpqTT0Tkcjg0g33kyBHCw8P5v//7P86ePUuPHj3o3r07Z86cYdq0aYSHh5OZmVnTWUVEpJ7w8DhfVAOcOQPvvQc//mhuLhGR6nCowJ4wYQK//vor69at48cff2TlypWsXLmSH3/8kfXr13P06FEmTJhQ01lFRKQe6dEDmjY1Pm/cGLy8IDwcpk8/v+qIiIgrcqjA/vjjj3n22WeJioq66Nh9993Hs88+yz//+c8rDmemuLg4WrRoQWBgIF26dLmsJQpFROTKle3Nfugh4+bHP/8Zpk2D3/3OWM6vfXv1Z4uI63GoB/vs2bM0a9as0uPXXHONWxekb731Fhs2bGDr1q00b96cH374AS8vL7NjiYjUOxf2Zs+YATExcPfdsHevMXbypFFka+1sEXEVDs1g33jjjbz77rsUVPA7uoKCAt577z1uuummKw5nhuLiYl577TUWLFhA8+bNAejQoYMKbBERF9G5M1x99fnHWVnG5jUiIq7CoQL7pZdeYseOHURERPCPf/yDzZs3s3nzZt5++226devG9u3beemll2o6a604fPgw586d44MPPqBZs2bccMMNvPPOO2bHEhGRMnr2PN+fDeDjc/6GSBERszlUYMfExJCQkMDRo0d58skn6d27N71792bMmDEcPXqUxMREYmJiajorYLSnTJkyhX79+hEUFITVaiUpKanCc/Pz85kwYQKhoaH4+fkRGRnJpk2bqvz6GRkZnDp1igMHDnDo0CE++OADJk2axNdff+2MtyMiIg4o258dEQHffw9jxxo7QIqImM3hjWhHjhzJ4cOH2bJlC++++y7vvvsuW7Zs4fDhw4wYMaImM5aTlZXF9OnTSU1NpXPnzgBYLJZKM8bHxzNs2DDmzJmDh4cHUVFRbCnzu8SlS5cSEBBAQEAAY8aMwc/PD4BXXnkFb29vOnbsyCOPPMJHH33ktPckIiKXLyEBUlNhxw74xz/grbfgD3+AoiKzk4lIfefwRjMADRo04LbbbuO2226rqTyXFBoaytGjRwkODmbnzp1ERERUeN6OHTt4//33mTlzJuPHjwdg2LBhdOjQgRdffNFeZA8dOpShQ4fan3f27NkK+60rK+JFRMR8o0dDQAAMHw6bNxtL+t1+u258FBFzODyDbRYvLy+Cg4MBsFXRcLd8+XI8PT0ZPXq0fczb25tRo0axbds2MjIyKnxew4YNiY6OZsaMGRQUFPDTTz+xbNmyCpckFBER1xEbC3feCenp8J//wNq1WsJPRMzhdgV2de3evZt27drh7+9fbrx0xnvPnj2VPvett97i+PHjNGnShP79+/Pqq6/So0cPp+YVEZErl55+/vPjx0G3z4iIGa6oRcSVZWZmEhISctF46diRI0cqfW6jRo1Yvny507KJiIhz9OhhrIudlWU8Liw0erI96+zVTkRcUZ2dwc7NzcXb2/uicR8fH/txERGpW8quLtK7N/zyi3HjY3Gx2clEpD6ps/+m9/X1JT8//6LxvLw8+/Ga1qdPHxo0aEBYWBhhYWEAxMbGEhsbW+OvJSIiFSt7Y+OyZUZvtpcXLFgA1jo7rSQiZaWkpJCSkgIYSzBnZGRQWFhYa69/xQV2SUkJp06dqvCGw6CgoCv98g4LCQmpsA0kMzMTMFYjqWkbN26kS5cuNf51RUTEMYMHQ0GBsbrI1q3GTHbPnlpdRKSuq2iCc9euXXTt2rVWXt+hArugoIDXX3+dhIQEfvnlF0oqWNnfYrFQbOLv5MLDw9m8eTM5OTkEBATYx7dv3w5gX0NbRETqtkcfhbffPr+d+q+/Gn3Zixebm0tE6i6HCuwnn3ySxMREbrvtNh588EEaNWp00TlmrxsdHR3NzJkzmT9/Ps8//zxg7OyYmJhIZGSkvYVDRETqvtKbHgFOn4bkZKPQvu8++Oor2LtXM9siUnMcKrCXLVvGsGHDKt2i3Nnmzp1Ldna2vQVkzZo1pP9vbaaxY8cSGBhIt27diImJYeLEiRw7dow2bdqQlJREeno6iYmJpuQWERFzlF1dJCgI2rY1xsePh9IOx99+M9bNVpEtIlfKoQLb19e3VndvvNCsWbM4dOgQYMyUr1y5khUrVmCxWBg+fDiBgYEALF68mMmTJ5OcnMzJkyfp1KkT69ato2fPnqZlFxGR2peQYBTPW7YYxXZpEd22rbEpDcCJE/Dll+ZlFJG6w6EC+/e//z3r1q3jiSeeqOk81ZKWllat87y9vYmLiyMuLs7JiURExNVVNDN9++1w6pQxs22xGO0jZ89Cw4a1n09E6g6HCuy//OUvDBs2jPvvv5/HHnuMa6+9Fg8Pj4vO04oaIiLiysrObLdvD59/DjExsHo1NGhgdjoRcVcOFdj5+fk0aNCADz/8kI8++qjCc8xeRURERKQ6ys5sb9wI/fvDqFGwaJHWzRYRxzhUYD/++OOsWLGC2NhYunXrVuEqIiIiIu6mTx9j+b7YWLjmGnjjDbMTiYg7cqjA/vjjj3nmmWeYPXt2TecREREx1SOPwLFj8OyzsHmz0Zdd9sZIEZFLcajADggIoG3pGkciIiJ1zNixsHAh/OtfxuPMTCgsNNbPFhG5FIe6yx5//HFSUlLUYy0iInVWbu75z3NyYMkSiI6GDz+EESOMmyIfe8y8fCLiuhyawe7QoQNr166lS5cuDB8+nBYtWlS4isjDDz98xQFFRETM0LMnZGef35zmuuvgv/81imyLxdig5uRJbU4jIhdzqMB+5JFH7J+/8MILFZ6jVURERMSdVbY5TevWcPCg8XlWlnFcRKQshwrszz77DIvFgq10f1kREZE6qKKZ6bvugrVr4fhxYxm/iIjazyUiru2yC+y8vDy+++47wsPD6dWrlzMyiYiIuKzSme3PP4eMDOPmR5vNaBsREQEHbnL09vbmpZdeIjU11Rl5REREXF5CAqSlwdKlsGyZseKIiEipyy6wLRYLN998MwdLG9BERETqqZgYGD3aWNZv3z6z04iIq3Bomb4ZM2bw9ttvs3HjxprOIyIi4lbi440VRoYMKb+0n4jUXw7d5Dh37lyuvvpq7r33Xq677jpat26Nr6/vReetWbPmigOKiIi4Mj8/eO8942bH8HCjH1s7P4rUbw4V2Hv37sVisdCiRQuKioo4cODARedY3PhuD39//3L5z507x8yZMxk3bpyJqURExFV16ABdusDWrcZjrY8tUr85VGDX9f7rM2fO2D/PzMykRYsW2jRHRESqlJVV/vNPPoGSEmMpPxGpX/TH/hKWLl1K9+7dadmypdlRRETEhfXsCU2aGJ97esKRI3DTTfD22zB8uLZWF6lPHJrBLrV582bWr19Peno6AC1btqR///51an3s5ORkxo4da3YMERFxcRfu/DhqFLz5Jjz55Pmt1U+cUOuISH3g0Ax2QUEBDz/8ML1792bWrFls3LiRDRs2MHPmTO666y4GDRpEYWFhTWcF4OzZs0yZMoV+/foRFBSE1WolKSmpwnPz8/OZMGECoaGh+Pn5ERkZyaZNm6r9Wt9//z0HDhwgJiampuKLiEgdlpAAqanGxx494MMPoVUro7gGY/dHLcAlUvc5VGBPmzaNVatW8ec//5nMzExOnDjByZMnyczM5IUXXmDlypVMmzatprMCkJWVxfTp00lNTaVz585A5TdUjhw5kvj4eIYNG8acOXPw8PAgKiqKLVu22M9ZunQpAQEBBAQEMGbMmHLPT05O5oEHHiAwMNAp70VEROq+u+6Cpk2Nzz094fBhGDcOzp0zN5eIOI/FZiv9d3X1tW7dml69erFo0aIKj48cOZLNmzc75WbIgoICsrOzCQ4OZufOnURERLBo0SKGDx9e7rwdO3YQGRnJzJkzGT9+PGDMaHfo0IHg4OByRXZFSkpKaNGiBf/4xz/o379/lefu2rWLrl27snPnTrp06XJlb1BEROqc0taR7t3hlltg0iS49lpYtAjeeed8W4laR0ScpzbrNYdmsDMzM4mMjKz0eLdu3cjMzHQ4VFW8vLwIDg4GoKp/GyxfvhxPT09Gjx5tH/P29mbUqFFs27aNjIyMKl/n008/pbCwkPvuu69mgouISL1V2jqSmGjMXu/ZA1dfbRTV770H+/fDunW6CVKkrnCowA4LC+Pzzz+v9PiXX35J8+bNHQ5VE3bv3k27du3w9/cvNx4REQHAnj17qnz+kiVLiI2Nxar1lUREpIa1bw9ff22sOlK6+2NWFqxZA3v3Gj3bjz2mlUdE3JVDq4iMHDmSKVOm0LhxY8aPH8/111+PxWJh//79zJ49m2XLljmtB7u6MjMzCQkJuWi8dOzIkSNVPr+yGydFRERqgocHDBhgFNW//QZeXnDqlNFC0rAhFBdDXp42rRFxRw4V2BMnTuTnn39mwYIFLFiwwD7LW1JSAsCIESOYNGlSzaV0QG5uLt7e3heN+/j42I+LiIiY6cKl/f7+d/jiCxgyBM6eNc7JyjKOi4j7cKjA9vT0ZNGiRYwbN46PPvqIQ4cOAefXwb7llltqNKQjfH19yc/Pv2g8Ly/Pfrym9enThwYNGhAWFkZYWBgAsbGxxMbG1vhriYhI3XDhzHTfvvDQQ7B2rbGsH8BVV9V+LhF3lpKSQkpKCgAZGRlkZGQ4bQnpilzRRjOdOnWiU6dONZWlRoWEhFTYBlJ682VoaGiNv+bGjRu1ioiIiFyx0pntr782Wke2bzdujpw502gtEZGqVTTBWbqKSG24ogIb4MyZM5w8ebLCFT1atGhxpV/eYeHh4WzevJmcnBwCAgLs49u3bwewr6EtIiLiisrObM+bB888A2lpsHSp0aMtIq7LoSUycnNzeemll2jatCmBgYG0bNmSVq1alfuvdevWNZ31skRHR1NcXMz8+fPtY/n5+SQmJhIZGWlv4RAREXF1Y8bA6tWwaZOxM2SbNlpdRMSVOTSD/dRTT7Fo0SIeeughevbsyVW13Bw2d+5csrOz7S0ga9asIT09HYCxY8cSGBhIt27diImJYeLEiRw7dow2bdqQlJREeno6iYmJtZpXRETkSt1/P9x9t7Fe9vHjcOKEMa7VRURcj0M7OTZu3JjBgweXmx2uTa1bt7bfWFm6TbrNZsNisZCWlmZvTcnPz2fy5MksWbKEkydP0qlTJ6ZPn06fPn1qNI92chQRkdrQvr2xKU2ppk3h11/hf5dCEalCbdZrDs1gWyyWWmsSr0haWlq1zvP29iYuLo64uDgnJxIREXG+Hj2MdbGzssDX1/gYGwsLFkCZ241ExGQO9WAPHDiQTZs21XQWERERqUJCgtEq0q4dPPIILFsG69dDRISxtJ92fhRxDQ7NYE+ePJnBgwfzxz/+kSeeeIIWLVrgUcG6QUFBQVccUERERM67sOf6llsgMhJSU43H2vlRxHwOFdht27YFYPfu3SxcuLDCcywWC8XFxY4nExERkUtq3x6aNIHsbONxVhZs2AA2m3qzRcziUIH9yiuvXPIci/5Ui4iI1Irbb4dTp4zi2tMTMjLgzjvh//4PevUyO51I/eNQgT116tQajiEiIiKOKt35ccsW40bIhx6CV14xiuy77wZvb/jPf4xjah0Rcb4r3slRREREzHdh4Xz//bBqlVF4l7aPqD9bpHY4tIqIiIiIuDaLxZjJbtr0/FhWFnz9tXmZROoLFdgiIiJ1WM+e5YtsLy/jBkgRcR4V2CIiInVY2bWzIyNh3z54/nkV2SLOpB5sERGROq5sz/W8efDUU0aB/eabWspPxBlUYIuIiNQjY8YYRfWYMcZ62YWFRhuJbnwUqTlqEREREalnnnwSbrsNfvwRDhyANWu0xbpITVKBLSIiUg/99lv5z9etM2azReTKqcAWERGph3r0OL+6iI+PsYRf587w6afm5hKpC1RgV2Dfvn3ccccdNG7cmDZt2rBw4UKzI4mIiNSosquLxMbCrl0QFAT33AOtWsF116ltRMRRusmxAiNGjODBBx/kyy+/ZPfu3fTq1YsePXpwww03mB1NRESkxlx4Y+OXX8Jdd8FXX0FJCfz6q9E2kpxsTj4Rd6UZ7Ar89NNPxMbGAhAeHs6NN95IamqqyalEREScy2KBzEyjuAY4dw5SUmDhQiguNjebiDtRgV2Be+65h+TkZIqKiti+fTvp6elERkaaHUtERMTpyvZmBwVBy5bw+OPwu9/B558bbSPt26t9RKQqbldgnz17lilTptCvXz+CgoKwWq0kJSVVeG5+fj4TJkwgNDQUPz8/IiMj2bRp0yVfY9asWSQmJuLr60vPnj2Ji4vjmmuuqem3IiIi4nLK9mYPHAg//wzbthk3QvbuDUuXwv79xqojKrJFKuZ2BXZWVhbTp08nNTWVzp07A2CpZBuqkSNHEh8fz7Bhw5gzZw4eHh5ERUWxZcsW+zlLly4lICCAgIAAxowZQ25uLvfccw/x8fEUFBSwa9cuXnrpJXbv3l0r709ERMRsCQmQmnq+RzsyErZuhWbNoKDAGMvKgjKXUxEpw+0K7NDQUI4ePUpaWhpvvPFGpeft2LGD999/n7/85S+8/vrrPP7443z22We0bNmSF1980X7e0KFDycnJIScnh3nz5vHDDz9QUFDAww8/jMVioWPHjnTv3p0vvviiNt6eiIiIS7JY4L774KqrjMceHhAebm4mEVfldgW2l5cXwcHBANhstkrPW758OZ6enowePdo+5u3tzahRo9i2bRsZGRkVPq9NmzacPn2aNWvWYLPZ+PHHH/nqq6+45ZZbavaNiIiIuJmEBHjwQWMZP29vY2m/w4fNTiXietyuwK6u3bt3065dO/z9/cuNR0REALBnz54KnxcUFMTSpUt5+eWXadSoEVFRUTz//PP07t3b6ZlFRERcXUICpKXB3r1Gu8gdd8DBg2anEnEtdXYd7MzMTEJCQi4aLx07cuRIpc8dOHAgAwcOdFo2ERERd3fddfDFF8aNj716QUSEUXT36HHx+toi9U2dncHOzc3F29v7onEfHx/7cREREXFcy5bG5jTZ2bBypVYXESlVZ2ewfX19yc/Pv2g8Ly/Pfrym9enThwYNGhAWFkZYWBgAsbGx9k1rRERE6pqwMGjSBE6fNh5rdRFxBSkpKaSkpACQkZFBRkYGhYWFtfb6dbbADgkJqbANJDMzEzBWI6lpGzdupEuXLjX+dUVERFxZr15w6hT89pvxuEULc/OIVDTBuWvXLrp27Vorr19nW0TCw8PZv38/OTk55ca3b98OYF9DW0RERK5MQgI88ABcfz1ce62x42NiotmpRMxTZwvs6OhoiouLmT9/vn0sPz+fxMREIiMj7S0cIiIicuUSEuDAAfjvf2HUKKMP+//9P6hiRV2ROsstW0Tmzp1Ldna2vQVkzZo1pKenAzB27FgCAwPp1q0bMTExTJw4kWPHjtGmTRuSkpJIT08nUf+sFhERcQpPT3j7bQgJgUmTIDkZioqgZ0+tLiL1h1sW2LNmzeLQoUOAsU36ypUrWbFiBRaLheHDhxMYGAjA4sWLmTx5MsnJyZw8eZJOnTqxbt06evbsaWZ8ERGROs1igalTYcMG2LbNGDt61Ci0Fy82NZpIrXDLAjstLa1a53l7exMXF0dcXJyTE4mIiMiFSm96BMjJgSVLoEEDGD4cFi2CrVu1brbUTW5ZYIuIiIjr69EDTp40lu4LCjJWF/n8c6OgtlqhpAROnDD6tVVkS11SZ29yFBEREXMlJMD990O7djBwIOzeDT//bKw0UlJinHP8uLEjpEhdogJbREREnCYhAVJTz89QWyxwzz3QtKnx2GqFI0fgq6/MyyhS01Rgi4iISK0qO7M9eDDceiv07g1//auW9ZO6QT3YIiIiUuvK9lwXFsJLL8Fzz8GOHeDhAdu36wZIcV8qsEVERMRUDRrArFnQrRs8+qgxi11cbNwgqRsgxR2pRURERERcwpAhEBZmFNdgrD6yZYu5mUQcoQJbREREXEbv3saSfmC0inTpYm4eEUeowBYRERGXkZBgLOnXqhV4ecHevcZMtog7UYEtIiIiLiUhAdLSYNcuY53se+4pvyukiKtTgS0iIiIu6YYb4NNPjXWy+/Y1boBs39648VHElWkVEREREXFZN98MmzYZK4x8/z0UFWl1EXF9msEWERERl9apEzRrZhTXoNVFxPWpwBYRERGXd/fd51cXAfD0hHPnzMsjUhUV2BXYt28fd9xxB40aNeLmm2/miy++MDuSiIhIvVa6ukjbtka7yH//C507ayZbXJMK7AsUFhYycOBABg8eTHZ2NnPmzCE6OpoTJ06YHU1ERKReS0iA/fuNbdT37IGrr4bbb4cOHYzCWzc/iqtQgX2B1NRUsrOzefrpp7FYLNx9992Eh4ezcuVKs6OJiIjI/7RvD19/DV27wr598J//wIoVKrLFNajAroDNZiv3uKSkhH379pmURkRERCri4QGnT59/fOoULFsG6enmZRIBNyywz549y5QpU+jXrx9BQUFYrVaSkpIqPDc/P58JEyYQGhqKn58fkZGRbNq0qcqv3759exo3bkx8fDyFhYV8/PHHfPnll5zTnRQiIiIup0cPaNrU+DwgAEpKjPWzX3sNRozQutliDrcrsLOyspg+fTqpqal07twZAIvFUuG5I0eOJD4+nmHDhjFnzhw8PDyIiopiS5k7IpYuXUpAQAABAQGMGTOGBg0asGrVKlavXk1ISAizZ89myJAhNG/evFben4iIiFRfQgLcfz+0awfR0fDrr/DUU/D//X+wdKnRs71unYpsqV0W24X9EC6uoKCA7OxsgoOD2blzJxERESxatIjhw4eXO2/Hjh1ERkYyc+ZMxo8fDxgz2h06dCA4OLhckX0pPXr0YMqUKfTt27fC47t27aJr167s3LmTLl26OP7mREREpEa0agWHDp1/3LatUWxL/VWb9ZrbzWB7eXkRHBwMXNwrXdby5cvx9PRk9OjR9jFvb29GjRrFtm3byMjIqPS5e/fuJS8vj3PnzvHGG29gs9kqLa5FRETE9fTufb51BIz+7Cou/SI1yu0K7OravXs37dq1w9/fv9x4REQEAHv27Kn0uYmJiYSEhBASEsKOHTtYtWqVU7OKiIhIzSrbOtK/P3h5QXi4se26iLN5mh3AWTIzMwkJCblovHTsyJEjlT73zTff5M0333RaNhEREXG+hITznx8/DkOHQt++MG2asVHN1q3GTZJlzxOpCXW2wM7NzcXb2/uicR8fH/txERERqR+aNIGPPoJXX4VXXoEGDaCwEE6eNG6AVJEtNanOFti+vr7k5+dfNJ6Xl2c/XtP69OlDgwYNCAsLIywsDIDY2FhiY2Nr/LVERETk8nh4wJQpsGDB+X7srCxYtcpYeaRLFxg1yth+XTPb7i0lJYWUlBQAMjIyyMjIoLCwsNZev84W2CEhIRW2gWRmZgIQGhpa46+5ceNGrSIiIiLi4vr2hbVrjbYRHx84dw5+9zvw94fiYsjN1cy2u6togrN0FZHaUGdvcgwPD2f//v3k5OSUG9++fTuAfQ1tERERqV8SEmDAAOMGyNhYOHMGNm4Eq9UorsGY2b6MFX1FyqmzBXZ0dDTFxcXMnz/fPpafn09iYiKRkZH2Fg4RERGpfxISIDXV+OjpCffcA4MGwVVXGce9vY02ERFHuGWLyNy5c8nOzra3gKxZs4b09HQAxo4dS2BgIN26dSMmJoaJEydy7Ngx2rRpQ1JSEunp6SQmJpoZX0RERFxQaTvIP/8JR49CZKS5ecR9uWWBPWvWLA79b3smi8XCypUrWbFiBRaLheHDhxMYGAjA4sWLmTx5MsnJyZw8eZJOnTqxbt06evbsaWZ8ERERcVGlRfaTT8LYsRARYayfLXI53LLATktLq9Z53t7exMXFERcX5+REIiIiUpfEx8OOHRATAzt3QqNGZicSd1Jne7BFREREHOXjAx98YKw08thjYLOZnUjciQpsERERkQpcdx0kJsKKFTBnjtlpxJ2owBYRERGpxEMPwfjxMG4ctGhhzGaLXIoKbBEREZEq/PabsQvkL7/Ahx+qyJZLU4EtIiIiUoVt26CoyPj89GlYvtzY6VGkMiqwRURERKrQowc0bWp8HhAAeXnQsSNs2GBuLnFdKrBFREREqpCQAPffb2ytHh0NP/8MN94I995rfGzbVm0jUp5broMtIiIiUptKN6Ap9cknxsz2N98Yj48dq/g8qZ80sIMASgAAIABJREFUgy0iIiJymaxWOHHi/OPTp+Hdd+Hzz83LJK5DBbaIiIiIA8r2ZjdqZPRn9+5ttI7s2mW0jbRvr/aR+kgFtoiIiIgDyvZmP/yw0Sby4Ydw6BB07WrMaO/fD+vWqciub9SDLSIiIuKgC3uuH34YHngAmjeHX381xrKyYMuW2s8m5tEMtoiIiEgN8vSEqCgICjo/1rKleXmk9qnAFhEREalhCQkwcCBcfz1cey18+im89ZbZqaS21NsC++9//ztdunTBy8uLadOmlTuWlZVF//798ff354YbbuCzzz4zKaWIiIi4q4QEOHAA0tLg2Wfh6aeNj8XFZicTZ6u3BXZoaCjTpk1j0KBBWCyWcseeeuopQkNDOX78OG+88QaDBw/mpPZEFREREQd4eMCbb8K8ecYsduvW2pymrqu3NzkOHDgQgI8++gibzWYfP3PmDKtXryYtLQ0fHx8GDBhAx44dWb16NSNHjjQprYiIiLi7J5+EVatg40aw2eDoUSgpgUWLzE4mNa3ezmBX5sCBA/j7+xMaGmof69ixI/v27TMxlYiIiNQFBw8axTXAmTPGUn4ffnh+TOoGlyywz549y5QpU+jXrx9BQUFYrVaSkpIqPDc/P58JEyYQGhqKn58fkZGRbNq0yeHXPnPmDIGBgeXGAgMDOXPmjMNfU0RERATKb05z1VUQHAzR0dCtG/Trp41p6gqXLLCzsrKYPn06qampdO7cGeCiPulSI0eOJD4+nmHDhjFnzhw8PDyIiopiS5kFJ5cuXUpAQAABAQGMGTOmytf29/fn9OnT5cZOnTpFQEDAFb4rERERqe/Kbk7z4INw+LCxvfqhQ/DJJ8bGNKtXq8h2dy7Zgx0aGsrRo0cJDg5m586dREREVHjejh07eP/995k5cybjx48HYNiwYXTo0IEXX3zRXmQPHTqUoUOHVvp6ZYv3tm3bcubMGY4cOWJvE9m7dy9/+MMfaurtiYiISD124eY0d94JjRsbG9IAnDgBH39c67GkBrnkDLaXlxfBwcEA5W5AvNDy5cvx9PRk9OjR9jFvb29GjRrFtm3byMjIqPS5xcXF5OXlUVRURGFhIXl5eZSUlODv78/AgQOZMmUKeXl5rF27lh9++MF+U6SIiIhITevZ83zriJcXZGbC448bfdriflyywK6u3bt3065dO/z9/cuNl85479mzp9LnTp8+HT8/PxYuXMiMGTPw8/NjyZIlAMybN48jR45w9dVX88ILL7Bs2TIaN27svDciIiIi9VrZ1pGhQ+GddyAlBbp0MbZeV2+2e3HJFpHqyszMJCQk5KLx0rEjR45U+typU6cyderUCo81adKE9evX10hGERERkeq4sHXk9tuhe3dYu9Z4fOKEUWRfeJ64Hreewc7NzcXb2/uicR8fH/txEREREXfUrh0EBZ1/fPw4bN5sWhy5DG49g+3r60t+fv5F43l5efbjtalPnz40aNCAsLAwwsLCAIiNjSU2NrZWc4iIiEjd0LMnZGcbN0BarfDLL5CUBMOHQyULrAmQkpJCSkoKABkZGWRkZFBYWFhrr+/WBXZISEiFbSCZmZkA5TaLqQ0bN26kS5cutfqaIiIiUnclJBhtIVu2QESEse36yJFG28jbb0OTJmYndE0VTXDu2rWLrl271srru3WBHR4ezubNm8nJySm3TvX27dsB7Gtoi4iIiLirC3uuBwyAP/0JOnY0/jt0yNjARr3ZrsOte7Cjo6MpLi5m/vz59rH8/HwSExOJjIy0t2mIiIiI1BXR0bB3rzGbvXGjsTnN2rVaZcSVuOwM9ty5c8nOzra3gKxZs4b09HQAxo4dS2BgIN26dSMmJoaJEydy7Ngx2rRpQ1JSEunp6SQmJpoZX0RERMRpQkPBz+/84+PH4bPPzMsj5blsgT1r1iwOHToEGDstrly5khUrVmCxWBg+fDiBgYEALF68mMmTJ5OcnMzJkyfp1KkT69ato2fPnmbGFxEREXGqsjdAengYN0D+7W/w9NO6AdJsLltgp6WlVes8b29v4uLiiIuLc3IiEREREddR9gbIyEho1AjGjjW2WQ8MhF271JttFpctsEVERESkahcWz/36wcMPQ0EB2Gxw8qQ2pzGDW9/kKCIiIiLnRUVBWJhRXIPRPrJli7mZ6iMV2CIiIiJ1SK9e5dfHbtfOvCz1lQpsERERkTokIcFYK/v666FpU2MGe+9es1PVLyqwRUREROqYhAQ4cMBYI7tlS+jTx3gstUMFtoiIiEgd1bgxfPKJ8fGee4yl/MT5VGCLiIiI1GHBwbBpk7E2docO0KaNdn10Ni3TJyIiIlLHNW8Ov/sdrFgBp08bOz/abKCNr51DM9giIiIi9cDeveeX7zt9Gt59F95/H4qLzc1VF6nAFhEREakHevQwVhUBoye7aVN45BGjbWTJEhg5Etq3V/tITVCBLSIiIlIPJCTA/fcb62I/9BAcPgzbtxvL+Q0bBsnJxqoja9eqyL5S6sEWERERqScu3DK9WzejoG7ZEtLTjbHjx2H9emO7dS+v2s9YF2gGW0RERKSeu/vu8+0j3t5w7Jgxs/3WWzBihFpHLpcKbBEREZF6rmz7yO9/D/v2we23w9NPG/3Z+/fDunUqsqur3hbYf//73+nSpQteXl5Mmzat2sdERERE6qKEBEhNNT7edBMsXQqtWkFJiXE8K8vYdl0urd4W2KGhoUybNo1BgwZhsViqfUxERESkvrjrrvOtIwA+PueX+pPK1dsCe+DAgQwYMIDGjRtju+AnpapjIiIiIvVFaetI27YQEQHff28s51dYaHYy16ZVRERERESkUmVXHklJMW56PHoUli+HgADzcrkyl5zBPnv2LFOmTKFfv34EBQVhtVpJSkqq8Nz8/HwmTJhAaGgofn5+REZGsmnTplpOLCIiIlL3xcbCxx/DN98Y/dlt2ujGx4q4ZIGdlZXF9OnTSU1NpXPnzgCV9kKPHDmS+Ph4hg0bxpw5c/Dw8CAqKootZbrwly5dSkBAAAEBAYwZM6ZW3oOIiIhIXdS7t9GbnZ0N//2vseX6/ferN7ssl2wRCQ0N5ejRowQHB7Nz504iIiIqPG/Hjh28//77zJw5k/HjxwMwbNgwOnTowIsvvmgvsocOHcrQoUMrfb2qbmTUTY4iIiIi5f300/nVRfLzjY1pOnSAP/4Rvv0W/vUvY2v2Cze2qS9ccgbby8uL4OBggCpvMly+fDmenp6MHj3aPubt7c2oUaPYtm0bGRkZlT63uLiYvLw8ioqKKCwsJC8vj5L//aRUdUxERESkvuvR4/zqIk2bwr33ws03w/jxxox2fd9y3SUL7OravXs37dq1w9/fv9x46Yz3nj17Kn3u9OnT8fPzY+HChcyYMQM/Pz+WLFlyyWMiIiIi9V3ZjWnuv9/oy162DFq3Pn/O8ePw9dfmZTSTS7aIVFdmZiYhISEXjZeOHTlypNLnTp06lalTp172MRERERGpuP2jVy/IyTE2pQFj3ez6yK1nsHNzc/H29r5o3Od//zdzc3NrO5KIiIhIvVV2ZrtbN9i7F+bNMztV7XPrGWxfX1/y8/MvGs/Ly7Mfr019+vShQYMGhIWFERYWBkBsbCyxsbG1mkNERETELGVntseNg6efhtBQePDB2suQkpJCSkoKABkZGWRkZFBYi7vjuHWBHRISUmEbSGZmJmCsRlKbNm7cSJcuXWr1NUVERERc1axZkJFhrJ/92Wdw222187oVTXDu2rWLrl271srru3WLSHh4OPv37ycnJ6fc+Pbt2wHsa2iLiIiISO2zWmHxYmOb9QEDYNAgaN++7q8u4tYFdnR0NMXFxcyfP98+lp+fT2JiIpGRkfY2DRERERExh48PrFoFRUXGx/37Yd26ul1ku2yLyNy5c8nOzra3gKxZs4b09HQAxo4dS2BgIN26dSMmJoaJEydy7Ngx2rRpQ1JSEunp6SQmJpoZX0RERET+JyjI+O/UKeNxVhaU2XS7znHZAnvWrFkcOnQIMHZTXLlyJStWrMBisTB8+HACAwMBWLx4MZMnTyY5OZmTJ0/SqVMn1q1bR8+ePc2MLyIiIiJl3HknnD4Nv/1mPG7RwtQ4TuWyBXZaWlq1zvP29iYuLo64uDgnJxIRERERRyUkGG0hX30FBQXGTY9z5xqrjNQ1Lltgi4iIiEjdUrqEX3ExvPgiPPMM/PwzZGfD1q3GFuwVbWDjblRgi4iIiEit8vAwlvC77jpjBtvLy5jVPnnSmOV29yLbrVcRERERERH39dRTxiY0BQXG47py86MKbBERERExzb33wlVXGZ9bLHDttebmqQkqsEVERETENAkJxjbqbdoYxfWnn8LLLxt92u5KPdgiIiIiYqrSnmubDeLiYNIk2LEDmjaFnTvd7+ZHFdgiIiIi4hIsFpgwwdhaPSoKCguhpMT9bn5Ui4iIiIiIuJTevY2bH0tKjMdZWcb26idOmJurulRgi4iIiIjLufNOaNLE+NzLC44fN4ru3/8e7rsP2rUzZrVdkQpsEREREXE5CQkwYIBRSA8dCpmZ8OqrsH49fPwxHDgAq1e7ZpGtHmwRERERcUkX9lz/+c8wfz6cPm08PnHC2HLd1WgGW0RERETcRs+exuoiYOwIefQo7NplbqYLqcAWEREREbeRkAD332+0jgwZArfcAnfdBVu3mp3sPLWIiIiIiIhbKds6cvq00avdty+sWWOsQGK2ejuD/fe//50uXbrg5eXFtGnT7OMFBQU89thjtGzZkkaNGnHbbbfxzTffmJhURERERCoTGAj//KfROhIVBX36QPv25t78WG8L7NDQUKZNm8agQYOwWCz28aKiIlq3bs2WLVs4deoUzz33HAMGDODs2bMmphURERGRyvj5GSuKXHMNbNoE+/ebu8JIvS2wBw4cyIABA2jcuDE2m80+7ufnx+TJk2nevDkAQ4YMwcvLi/3795sVVUREREQuwdvb+K/UiRPw4Yfwn//Ufhb1YF/CgQMHOHHiBNdff73ZUURERESkCj17Qna2sfNjw4ZQUGC0iwwebKyjXVtccgb77NmzTJkyhX79+hEUFITVaiUpKanCc/Pz85kwYQKhoaH4+fkRGRnJpk2baiRHbm4ujz76KJMmTSIgIKBGvqaIiIiIOEfZFUYGD4aTJ2HuXGOb9S++qL0cLllgZ2VlMX36dFJTU+ncuTNAuT7pskaOHEl8fDzDhg1jzpw5eHh4EBUVxZYtW+znLF26lICAAAICAhgzZky1MhQWFhITE0O7du2YPHnylb8pF5OSkmJ2hMvibnnB/TIrr/O5W2bldS53ywvul1l5nc8VMyckQGqq8dHHB558EkJCajeDSxbYoaGhHD16lLS0NN54441Kz9uxYwfvv/8+f/nLX3j99dd5/PHH+eyzz2jZsiUvvvii/byhQ4eSk5NDTk4O8+bNu+jrXFi8l5SUMGzYMDw8PCqdOXd3rvgHoirulhfcL7PyOp+7ZVZe53K3vOB+mZXX+dwlc8+e0Lhx7b2eSxbYXl5eBAcHA5S7AfFCy5cvx9PTk9GjR9vHvL29GTVqFNu2bSMjI6PS5xYXF5OXl0dRURGFhYXk5eVRUlICwJ/+9CeOHj3KsmXLsFpd8lskIiIiItWUkAC33157r+fW1ePu3btp164d/v7+5cYjIiIA2LNnT6XPnT59On5+fixcuJAZM2bg5+fHkiVLSE9PZ+HChXz77bc0adLE3lpStuWkLqjqHx+uyN3ygvtlVl7nc7fMyutc7pYX3C+z8jqfO2WeOrX2XsutVxHJzMwkpIKmmtKxI0eOVPrcqVOnMrWS73TpTHZd5k5/IMD98oL7ZVZe53O3zMrrXO6WF9wvs/I6nztmrg1uXWDn5ubiXXbBw//x8fGxH6+tHAA//fRTrbxeTSgsLGTXrl1mx6g2d8sL7pdZeZ3P3TIrr3O5W15wv8zK63zulLm0TquN+tCtC2xfX1/y8/MvGs/Ly7Mfrw0HDx4E4NFHH62V16spXbt2NTvCZXG3vOB+mZXX+dwts/I6l7vlBffLrLzO526ZDx48SI8ePZz6Gm5dYIeEhFTYBpL5v5XEQ0NDayXHvffey5IlS2jVqlWtFfUiIiIiUn25ubkcPHiQe++91+mv5dYFdnh4OJs3byYnJ6fcRjDbt28HsK+h7WxNmjRh6NChtfJaIiIiIuIYZ89cl3LrVUSio6MpLi5m/vz59rH8/HwSExOJjIwkLCzMxHQiIiIiUh+57Az23Llzyc7OtreArFmzhvT0dADGjh1LYGAg3bp1IyYmhokTJ3Ls2DHatGlDUlIS6enpJCYmmhlfREREROorm4tq1aqVzWKx2CwWi81qtdqsVqv980OHDtnPy8vLs73wwgu2kJAQm4+Pj+3WW2+1bdiwwcTk5c2bN88WHh5ua9CggW3q1KnljvXq1cvm4+Nj8/f3t/n7+9uioqJMSnleVXlLbd261WaxWGyvvvpqLaerWFWZ//jHP9qaNWtmCwgIsLVv3962cOFCk1KeV1ne/Px82x/+8AdbixYtbIGBgbbIyEjbtm3bTExqqOr7W52fFzNUlevYsWO2qKgoW8OGDW3t27e3ffrppyalrNwPP/xgu/32222BgYG2m266ybZ582azI1WpNG+jRo1s1113ne2dd94xO1KVGjZsaP9719/f32a1Wm1vvvmm2bGq9Prrr9uuvfZaW0BAgC08PNyWk5NjdqQqueL1rTpc7fpWGVe8tlXFVa9vVbnS65vLtoikpaVRUlJCSUkJxcXFFBcX2z9v0aKF/Txvb2/i4uI4cuQIubm5fPPNN/Tp08fE5OWFhoYybdo0Bg0adNGW7BaLhYULF9q3cV+/fr1JKc+rKi8Ya4SPGzeOW2+9tcLjZqgq87hx40hLS+P06dMsWbKEp556yr7qi1kqy1tUVETr1q3ZsmULp06d4rnnnmPAgAGcPXvWxLRVf38v9fNilqpyPfXUU4SGhnL8+HHeeOMNBg8ezMmTJ01KerHCwkIGDhzI4MGDyc7OZs6cOURHR3PixAmzo1VqxIgR9O3bl+zsbJYvX864ceP497//bXasSp05c8b+9+7+/fuxWq08/PDDZseq1FtvvcWGDRvYunUrp0+fZvHixXh5eZkdq0queH27FFe8vlXGFa9tVXHV61tVrvT65rIFdl0xcOBABgwYQOPGjSvc9r2iMTNdKu/8+fOJjIzkhhtucJnsVWW+8cYb7euiAwQGBl6082dtqyyvn58fkydPpnnz5gAMGTIELy8v9u/fb1ZUoOrv76V+XsxSWa4zZ86wevVqpk2bho+PDwMGDKBjx46sXr3axLTlpaamkp2dzdNPP43FYuHuu+8mPDyclStXmh2tUj/99BOxsbGAcfP5jTfeSGpqqsmpqmfp0qV0796dli1bmh2lQsXFxbz22mssWLDA/ndDhw4dXL7ABte7vl2KK17fKuOK17aquOr1rSpXen1TgW2ycePGERwcTN++fdm7d6/Zcar022+/8de//pVp06aZHeWyjBkzBj8/P26//XYWLFhAkyZNzI5ULQcOHODEiRNcf/31ZkepMw4cOIC/v3+5JTw7duzIvn37TEx1sQv/Mi8pKXG5jGXdc889JCcnU1RUxPbt20lPTycyMtLsWNWSnJzM8OHDzY5RqcOHD3Pu3Dk++OADmjVrxg033MA777xjdqxq0fXNudz12gb14/pWrwvss2fPMmXKFPr160dQUBBWq5WkpKQKz83Pz2fChAmEhobi5+dHZGQkmzZtuqLXj4uL4+DBg6Snp9O3b1/uu+8+zpw547J5X375ZcaNG0ejRo0AqvUrE7MzA8ybN4+zZ8+ybNkyHnvsMX755ReXzgvGWp2PPvookyZNKrcEpavmvRxmZj5z5gyBgYHlxgIDA6v8c1fb+du3b0/jxo2Jj4+nsLCQjz/+mC+//JJz585VO2Nt5gWYNWsWiYmJ+Pr60rNnT+Li4rjmmmtcNm+p77//ngMHDhATE+NQ1trIm5GRwalTpzhw4ACHDh3igw8+YNKkSXz99dcumxku//pmdl5Hrm9m5oXLu7a5Smao/vXNVfI6ql4X2FlZWUyfPp3U1FT7mtmV/aEaOXIk8fHxDBs2jDlz5uDh4UFUVBRbtmyxn7N06VICAgIICAhgzJgxl3z9iIgI/Pz88PHx4c9//jOBgYF88803Lpl39+7d/Otf/+Lxxx8HjBm26vzKxOzvcSmLxcKAAQPo3r17le0ArpC3sLCQmJgY2rVrx+TJk6s81xXyXi4zM/v7+3P69OlyY6dOnbqsv+Sdnb9BgwasWrWK1atXExISwuzZsxkyZIj9V6uXy9l5c3Nzueeee4iPj6egoIBdu3bx0ksvsXv3bpfMW1ZycjIPPPDARf/ocqW8fn5+ALzyyit4e3vTsWNHHnnkET766COXzQyXf30zM6+j1zez8pZV3Wubq2S+nOubK+S9Ild+n6X7ys/Pt/366682m81m+9e//mWzWCy2pKSki87bvn27zWKx2GbNmmUfy8vLs11//fW27t27V+u1nnjiCdu0adOqPOfGG2+0bdy40SXzzp492+bv729r1qyZrVmzZjZfX19bQECA7bHHHqvy67ja97hfv362hIQEl81bXFxsGzJkiO2BBx6wFRcXX/JrmJ23OsdcKXNOTo7Ny8vLlpGRYR/r1auXbdGiRdX6erWdv1T37t1tn3zyyWU9p7by7tixwxYSElJuLDo62hYfH++SeUsVFxfbwsLCbOvWrXMoZ23lPXPmjM3b29uWnp5uH3vmmWdskyZNctnMFbnU9c3MvI5e38zKW5FLXdtcIfPlXt/Mzlvqcq5vZdXrGWwvLy+Cg4OBqm/GWL58OZ6enowePdo+5u3tzahRo9i2bRsZGRmVPre4uJi8vDyKioooLCwkLy+PkpISTp06xcaNG8nPz6egoID4+Hiys7O59dZbXTLv6NGj+fnnn/nuu+/Ys2cPDzzwAE8//TTx8fGVfi2zM58+fZp3332Xs2fPUlRUxAcffHDJVWbMzAvwpz/9iaNHj7Js2TKs1kv/8TQ7b1XHXDGzv78/AwcOZMqUKeTl5bF27Vp++OEHBg4cWGXm2s6/d+9e8vLyOHfuHG+88QY2m42+fftWO2Nt5m3Tpg2nT59mzZo12Gw2fvzxR7766ituueUWl8xb6tNPP6WwsJD77rvPoZy1lbdhw4ZER0czY8YMCgoK+Omnn1i2bBlRUVEum9mR65uZeR29vpmV15Frm9mZ4fKvb2bndeT6Vla9LrCra/fu3bRr1+6iO3QjIiIA2LNnT6XPnT59On5+fixcuJAZM2bg5+fHkiVLKCwsZNKkSTRt2pSQkBDWr1/PRx995HA/krPz+vr6EhwcTHBwMNdccw2+vr74+/tf0a9WnZ3ZYrHwzjvv0Lx5c4KDg5kzZw5r1651+Fftzs6bnp7OwoUL+fbbb2nSpIn911Vlf43lSnkvdcxVM8+bN48jR45w9dVX88ILL7Bs2TIaN25cI5lrKn9iYiIhISGEhISwY8cOVq1aVeP5LuRo3qCgIJYuXcrLL79Mo0aNiIqK4vnnn6d3794umbfUkiVLiI2NveILfXVdSd633nqL48eP06RJE/r378+rr75aK9s9O5rZmdc3Z+R19vWtpvM689rmrMyHDh1y2vXNGXnhyq9vLruToyvJzMwkJCTkovHSsdLdJisydepUpk6dWuGxb7/9tkbyXchZecuq6Z0ynZX5s88+q5F8F3JW3sv51/HlcFbe6v68OMJZmZs0aVIra/JeSf4333yTN99802nZKnIleQcOHHhZvwWoCVeSF6j0hihnuZK8jRo1Yvny5U7LVhlHMzdp0sRp17eqXOnPRKna2gna0bwBAQFOu7ZdiqOZW7Zs6bTrW1Vqox6qjGawqyE3Nxdvb++LxkvXoMzNza3tSFVyt7zgfpmV1/ncMXNZ7pZfeZ3L3fKC+2VWXudzt8xm5lWBXQ2+vr7k5+dfNJ6Xl2c/7krcLS+4X2bldT53zFyWu+VXXudyt7zgfpmV1/ncLbOZeVVgV0NISEiFv0bIzMwEKLdphStwt7zgfpmV1/ncMXNZ7pZfeZ3L3fKC+2VWXudzt8xm5lWBXQ3h4eHs37+fnJyccuPbt28HsK/B6CrcLS+4X2bldT53zFyWu+VXXudyt7zgfpmV1/ncLbOZeVVgV0N0dDTFxcXMnz/fPpafn09iYiKRkZGEhYWZmO5i7pYX3C+z8jqfO2Yuy93yK69zuVtecL/Myut87pbZzLz1fhWRuXPnkp2dbf8Vwpo1a0hPTwdg7NixBAYG0q1bN2JiYpg4cSLHjh2jTZs2JCUlkZ6eXmt3G7trXnfMrLzKfCnull95ldfdMyuvMrtd3svemqaOadWqlc1isdgsFovNarXarFar/fNDhw7Zz8vLy7O98MILtpCQEJuPj4/t1ltvtW3YsEF562Bm5VXmS3G3/MqrvO6eWXmV2d3yWmy2KrbAERERERGRy6IebBERERGRGqQCW0RERESkBqnAFhERERGpQSqwRURERERqkApsEREREZEapAJbRERERKQGqcAWEREREalBKrBFRERERGqQCmwRERERkRqkAltEREREpAapwBYRERERqUEqsEXEbXz77bd0796dhg0bYrVa+e6775g6dSpWq2N/lS1atAir1Up6enoNJ62+zZs3Y7Va+fLLL+1jI0eOpHXr1uXOO3PmDI8//jjNmjXDarUybtw4AH799Veio6O5+uqrsVqtzJkzp1bzi4jIxTzNDiAiUh2FhYXExMTg5+fHX//6V/z8/GjZsiUWiwWLxeLw173wufPmzaNhw4aMGDHiSiM7rKL39Nprr5GUlMQrr7xCmzZtuPHGGwEYN24cGzZsYOrUqTRr1ozf/e53ZkR2G1u3bmXjxo0899xzNGrUyOw4IlJHWWz8HqrbAAAMoklEQVQ2m83sECIil/Lvf/+bm266iXfeeYfHHnvMPl5cXExxcTFeXl6X/TVLSkooKioq99wOHTrQtGlTPv/88xrJfSmbN2+md+/ebN68mTvuuAOAoqIibDYbDRo0sJ8XGRmJl5dXuZlugGbNmtG3b18WL15cK3nd3cyZM3nxxRc5ePAgLVq0qPXXHzJkCOvWravwty4lJSXExMSwaNGiWs8lIjVLLSIi4haOHTsGcNGso4eHh0PFNYDVanX4uc7k6elZrrgG4/1XNOOalZVVozOxRUVFFBYW1tjXc1VmzS2VlJSwdu1acnJyLvpvxYoVFBcXm5JLRGqWCmwRcXkjR47kzjvvBCAmJgar1Urv3r0BKuzBtlqtPPPMM6xatYoOHTrg4+NDhw4d+OSTT8qdd2EPdqtWrfjxxx/54osvsFqtWK1W7rrrLvv52dnZPPfcc1x77bX4+PjQtm1b4uLiql2sHT58mAcffJCGDRtyzTXXMH78ePLz8yt8v6U92KU92gcPHmT9+vX2XElJSVitVmw2G2+99ZZ9/HKyHjx4EKvVyqxZs5g9ezZt2rTBx8eHn376CTB+a1Da3+3r60tERARr166t8Hu4detWxo8fT9OmTfH39+fhhx/m+PHjF723f/7zn/Tq1YvAwEAaNWpEt27dSElJKXfO9u3b6devH40bN6Zhw4bceeedbN26tVrf47/97W/cfPPNNGzYkKCgICIiIuxff+rUqbz44osAtG7d2v49K9uDv2TJErp27Yqfnx9XX301sbGxHD58uNxr3HnnnXTs2JGdO3fSvXt3/Pz8uO666/jHP/5RrYyV/bzoF8oidYd6sEXE5T3xxBM0b96c1157jWeffZaIiAiuueYa+/GKerC//vprVqxYwVNPPYW/vz9z5sxh0KBBpKenExQUVOHr/PWvf+WZZ54hICCAl19+GcD+OufOnaNXr15kZmbypz/9iRYtWrBlyxYmTpxIZmYm8fHxVb6H3Nxc7r77bg4fPszYsWMJCQkhOTmZTz/9tMLzS9/TTTfdRHJyMuPGjePaa6/l+eefB6Bz584kJyczbNgw+vbty/Dhw+3PvdysiYmJ5Ofn88QTT+Dt7c1VV13Fvn376NGjB9deey0TJ06kYcOGvP/++zz44IN8+OGHPPjgg+W+xjPPPENQUBDTpk0jLS2N2bNn8/TTT/Pee+/Zz1m0aBGPPfYYHTt2ZNKkSTRu3Jhdu3bxySefEBsbC8Bnn33GfffdR0REBFOnTsVisZCYmEjv3r356quviIiIqPR7vGDBAp599lliYmIYN24ceXl5fPfdd+zYsYPY2FgGDRrEgQMHSElJYfbs2TRp0gTA/nHGjBm88sorDBkyhNGjR3Ps2DH+9re/cccdd7B79277bwosFgsnT56kf//+DBkyhKFDh/L+++/z5JNP4uXlxR/+8Icqfxau5J4BEXETNhERN/D555/bLBaL7cMPPyw3PmXKFJvFYik3ZrFYbD4+Prb//ve/9rHvv//eZrFYbHPnzrWPJSYm2iwWi+3QoUP2sZtvvtl21113XfT606dPt/n7+9v+85//lBufOHGizdPT0/bLL79UmX/27Nk2i8ViW758uX3s3LlztrZt29osFovtiy++sI+PGDHC1qpVq3LPb9mypW3AgAEXfV2LxWJ75plnHMqalpZms1gstsaNG9uOHz9e7ty7777b1qlTJ1tBQUG58R49etjatWv3/7d3ZyFRfm8cwL+DTj/GaawpclRoc0YrQsjSJJFcijExstzCimjTCoooENs0lzQrYS7SyuzCsoIwIyTKFpGuvFHQLMI0dAhbLE1z1HI7v4twcLbchv//l30/4MV73vc955kzRI+H5xyN1yNzqNVqTZ47duyYcHR0FN+/fxdCCNHZ2SkUCoVYs2aN+Pnzp9U5Gh4eFp6eniI8PNykva+vT3h4eFiMYS4yMlJ4e3v/9pmLFy9afOdCCNHS0iIcHBxETk6OSfurV6+EVCoV2dnZxragoCAhkUiETqcztvX39wsfHx+hUqnEwMCAzfFjYmJERUWF1XuPHz8WO3bs+G38RPRnYIkIEU1L69evNznqztvbG87Ozmhubp5UfyUlJVi7di1mz56Nr1+/Gn/WrVuHoaEhi82H5h49egR3d3dER0cb22QyGRITEycVjz1jjY6Oxty5c43XHR0dqKysRGxsLLq6ukz60Gq1aGxsxMePH036MP8cgYGBGBoagl6vBwA8e/YMBoMBx48ft1n3Xltbi6amJsTHx5uMaTAYEBoaOuYcK5VKvH//HtXV1eOeqxH379+HEAIxMTEmY6tUKmg0GotNr1KpFPv377e4bmtrQ01NzYTHJ6LphSUiRDQtWTshQqlU4tu3b5Pqr7GxEfX19Zg3b57FPYlEgi9fvgD4telw9EY1hUIBuVwOvV4PjUZj8a6Xl9ek4rFHrCPMz9xuamqCEAIpKSlISUmx2kdbWxvc3NyMbebzrVQqAcA43+/evQPw65SW38UNwOYRiRKJBF1dXTY3dSYnJ+P58+dYvXo1NBoNtFottm3bhoCAAJtjjh5bCAFPT0+r9//55x+Ta3d3d8hkMpO2kXf1ej38/f3HHJOIpi8m2EQ0LTk4OFhtF5PcSCaEgFarNW6SMzeSKPv5+ZlsmktLS0NqauqkxpyssWI1TyLNE8Xh4WEAQFJSEsLCwqz2oVarTa7tMd8j4+bm5mLFihVWn5HL5TbfX7p0KRoaGvDw4UOUl5ejtLQUly9fRmpqKtLS0sYcWyKRoLy83OpnmTlz5rg/BxERE2wiolFsbUBTq9Xo7u42nl5iy507d/Djxw/jtYeHBwBg4cKFeP36tcXzDQ0NU4jWuvHGastIzI6OjpPuw1pMAFBfX2/s39YzCoVi0uM6OTkhLi4OcXFxGBgYQFRUFLKysnDy5EnMmDHD5ver0WgghMCiRYtsrmKP1trait7eXjg5ORnb3r59C+DXaTRE9HdjDTYR0ShyudxqGUlcXByqqqrw9OlTi3udnZ3GspCAgACEhoYaf0aSrYiICHz48AH37t0zvtfb24tr165ZjWMqJ02MN1ZbXFxcEBwcjIKCAnz69MnivnmJyXiEhYVBoVDg3LlzVo8mBABfX1+o1Wrk5uaip6dnwuO2t7ebXEulUuNfvBw523tkBdz8O46KioKDgwPS09Mt+hVCoKOjw6RtcHDQ5Fi+/v5+FBQUwMXFBatWrfptnEQ0/XEFm4j+GuMpV/D19cWVK1eQlZUFtVoNlUqFkJAQJCUloaysDBs3bsSuXbuwcuVK9PT0oL6+HqWlpdDr9TaP/wOAhIQE5OXlYefOnaipqYGrqyuKi4ttljxMtpQFwJRjBYD8/HwEBgbC29sbCQkJWLx4MT5//oyqqiq0traitrZ2QjEpFArodDrs27cPfn5+iI+Ph1KpRF1dHfr6+lBUVASJRILr168jPDwcy5cvx+7du+Hu7o7W1lZUVlZi1qxZKCsrszmGVquFm5sbAgICoFKp8ObNG+Tn5yMiIsI4zyN/Sv7UqVPYunUrpFIpNm3aBA8PD5w9exYnTpxAS0sLIiMjoVAo0NzcjAcPHiAxMdF4RCLwqwb7/PnzaGlpgaenJ+7evYu6ujoUFhbaLJchor8HE2wi+mNYW9WVSCTjXu219f5oqamp0Ov1uHDhArq7uxEcHIyQkBDIZDK8ePEC2dnZKCkpwc2bN+Hs7IwlS5YgIyMDzs7Ovx1bJpOhoqIChw8fxqVLlyCXy7F9+3Zs2LAB4eHhY36miaxoTzVWAFi2bBmqq6uRnp6OoqIitLe3Q6VSwcfHB2fOnBlXbObte/bsgYuLC3JycnD27FnjCvPRo0eNzwQFBaGqqgqZmZnIy8uDwWCAm5sb/P39TU7tsObAgQO4ffs2dDodDAYD5s+fjyNHjuD06dPGZ3x9fZGZmYmrV6+ivLwcQgg0NzdjwYIFSE5OhpeXF3Q6HTIyMgD82rwZFhaGyMhIk7HmzJmDGzdu4NChQygsLISrqyvy8/Oxd+/eMed2Kr88EdGfQSL4L52IiGjcgoOD0dHRgZcvX0743djYWBw8eNBqjfmTJ09w69YtFBcX2yNMIvo/Yg02ERHR/5CtdS2udxFNHywRISIimqCpJMObN2+Go6Plf7+Dg4PYsmXLVMIiov8IJthEREQTMJG6f3MlJSV2joaI/otYg01EREREZEeswSYiIiIisiMm2EREREREdsQEm4iIiIjIjphgExERERHZERNsIiIiIiI7YoJNRERERGRHTLCJiIiIiOyICTYRERERkR0xwSYiIiIisqN/AeRKZEaYzzdPAAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x31670ed10>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "PyObject <matplotlib.text.Text object at 0x31730e050>"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using PyPlot\n",
    "δ = logspace(-15, -1, 100)\n",
    "loglog(δ, [norm(recurrence_adjoint(A, b, p)[2] - recurrence_FD(A, b, p, d)[2]) for d in δ], \".-\")\n",
    "xlabel(\"finite-difference step δ\")\n",
    "ylabel(\"norm of error in dg/dp\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At first glance, this seems ridiculous: the errors get *monotonically worse* as $\\delta$ decreases!  How can this be right?\n",
    "\n",
    "The key point is to realize that our $g(p)$ depends *quadratically* on $p$, and a center-difference approximation is *exact* for a quadratic function *in exact arithmetic*.   So, *all* of the difference arises due to roundoff errors.  In the finite-difference approximation, the difference $g(p + \\delta e_n) - g(p - \\delta e_n)$ is the difference of two nearly equal quantities as $\\delta \\to 0$, so there is a large cancellation error in the result (i.e. the result arises from the *least* significant digits of $g$), and hence the error worsens with decreasing $\\delta$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.0-dev",
   "language": "julia",
   "name": "julia 0.4"
  },
  "language_info": {
   "name": "julia",
   "version": "0.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}