Bayes, Kalman, and Sobol

HW7.ipynb

{
  "cells": [
    {
      "metadata": {
        "collapsed": true,
        "trusted": true
      },
      "cell_type": "code",
      "source": "import numpy as np\nimport matplotlib.pyplot as plt\nimport pandas as pd\nimport seaborn as sns\n%matplotlib inline",
      "execution_count": 1,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "# Q1\nThe Boy-Boy paradox (as made famous by Martin Gardner)\n\n## Part 1"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "- child 1 = boy\n    - child 2 = boy\n    - child 2 = girl\n- child 1 = girl\n    - child 2 = boy\n    - child 2 = girl\n    \nlet \"A\" be \"Both are boys\" and \"B\" be \"one is a boy\". Then let our prior assume 4 possible unique ways to have children, and all be equally likely: {BB, BG, GB, GG}. Thus, P(A)= P(BB)= 0.25. Additionally, P(B) = P(BB or BG or GB) = 0.75. \n\nAs a consequence, P(B|A) is obviously 1.0 (one is definitely a boy if both are). We can now update our posterior via bayes' rule. "
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "def bayes(a,b,b_a):\n    a_b = b_a*a/b\n    return a_b\n\nA = .25\nB = .75\nB_A = 1.\nprint(f'P(A|B) = {bayes(A, B, B_A):.2f}')",
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "text": "P(A|B) = 0.33\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## Part 2\nNow The issue becomes modifying one of the terms to include being on a tuesday: B|A = P(B,T|BB). If all days are equal, probablility of something on tuesday is $t=1/7$. The chance that neither boy (given that their are two of them)is born on tuesday is $(1-t)^2$ (not tuesday, twice). Then \"at least one on tuesday\" is \"not both\", so P(B,T) = $1-(1-t)^2$\n\nP(B,T) is the chance of *having a boy on a tuesday*, so we again compare that over all possible ways to have children {BB,BG,GB,GG}. Having a boy on a tuesday is impossible with two girls, and $t$ with only one boy to pick from (i.e. two of the cases). With two, we can't double count on the union (one boy plus the other boy, not both together since that's already counted), so it's $t+t-t^2$. Each of these four options is equally likely (1/4), so P(B,T) = $0.25t+0.25t+0.25(t+t-t^2)$"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "A = .25\nt=1/7\nB = .25*t+.25*t+.25*(t+t-t**2)\nB_A = 1-(1-t)**2\nprint(f'P(A|B) = {bayes(A, B, B_A):.2f}')",
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": "P(A|B) = 0.48\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "collapsed": true
      },
      "cell_type": "markdown",
      "source": "# Q2 \n## Find $P(XY-2000<0)$\n\nLet $$X = \\mathcal{N}(38, 3.8)$$ $$ Y = \\mathcal{N}(54, 2.7) $$\n\n$$G(X,Y) = XY - 2000$$\n\nThen, we can find the gradient: \n$\\nabla G = \\{\\frac{\\partial X}{\\partial X}Y + \\frac{\\partial Y}{\\partial X}X, \n              \\frac{\\partial X}{\\partial Y}Y + \\frac{\\partial Y}{\\partial Y}X \\}\n          = \\{Y, X\\}$\n\nIn Rackwitz-Fiessler, we don't want to use the gradient directly, but one calculated in the standard normal space. If some $z_i$ is the standard normal representation of a variable $x_i$, the gradient $$\\nabla_z G = \\sigma_x\\odot \\nabla_x G$$ where $\\sigma_x$ is the standard deviation vector of the $x$ variables and $\\odot$ is the Hadamard (elementwise) product.\n\nOur gradients now become $$\\alpha = \\sigma_x\\odot \\nabla G = \\{\\sigma_X Y, \\sigma_Y X\\} = \\{3.8Y, 2.7X\\}$$\n"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "from sympy.stats import P, E, variance, Normal, std, Expectation, cdf, Uniform\nfrom sympy import *",
      "execution_count": 4,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "x, y = symbols('x, y')\ng = x*y - 2000\nprint(f'Part A:\\n\\tfunction: {g} \\n\\tgradient: [{g.diff(x)}, {g.diff(y)}]')\n",
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Part A:\n\tfunction: x*y - 2000 \n\tgradient: [y, x]\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "from scipy.stats import norm\n\nmu = np.array([38., 54.]).reshape(-1,1)\nsd = np.array([3.8, 2.7]).reshape(-1,1)\n\n# xi = mu.copy()\nxi = np.array([1,1]).reshape(-1,1)  # initialize\ng = lambda x: x[0]*x[1] - 2000\ngrad = lambda x: np.flipud(x)\n\ndef r_f(xi, mu, sd, g, grad):\n    ui = (xi - mu)/sd\n    beta = np.sqrt(ui.T.dot(ui))\n    diff = 1.\n    print('Optimizing...')\n    while diff >= 1e-6:\n        ui = (xi - mu)/sd\n        alpha = grad(xi)*sd\n        ui = (1/np.linalg.norm(alpha))*(alpha.T.dot(ui) - g(xi))*alpha/np.linalg.norm(alpha)\n\n        beta_i = np.sqrt(ui.T.dot(ui))\n        diff = np.abs(beta - beta_i)\n        beta = beta_i\n        \n        xi = ui*sd + mu\n        print(f'beta: {beta[0,0]:.2f}, \\tx:{xi[0,0]:.2f}, \\ty:{xi[1,0]:.2f}, \\t diff={diff[0,0]:.2f}')\n    print(f'done!\\nProbability CDF(-beta): {norm.cdf(-beta[0,0]):.2f}')\nr_f(xi, mu, sd, g, grad)",
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Optimizing...\nbeta: 409.52, \tx:1306.57, \ty:694.43, \t diff=387.61\nbeta: 184.40, \tx:457.73, \ty:452.68, \t diff=225.12\nbeta: 78.98, \tx:281.74, \ty:178.42, \t diff=105.42\nbeta: 29.71, \tx:113.12, \ty:113.88, \t diff=49.27\nbeta: 8.39, \tx:64.06, \ty:67.07, \t diff=21.32\nbeta: 0.94, \tx:40.95, \ty:55.42, \t diff=7.46\nbeta: 0.20, \tx:37.32, \ty:53.75, \t diff=0.74\nbeta: 0.23, \tx:37.22, \ty:53.73, \t diff=0.03\nbeta: 0.23, \tx:37.22, \ty:53.73, \t diff=0.00\nbeta: 0.23, \tx:37.22, \ty:53.73, \t diff=0.00\ndone!\nProbability CDF(-beta): 0.41\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "This makes a lot of sense... just the max. likelihood value alone is 38\\*54 = 2052, which is fairly close to 2000. \n\n## Find $P(X-Y/2000 <0)$"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "g = x/y - 1/2000  # rearranging a bit to match previous problem\nprint(f'Part B:\\n\\tfunction: {g} \\n\\tgradient: [{g.diff(x)}, {g.diff(y)}]')",
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Part B:\n\tfunction: x/y - 0.0005 \n\tgradient: [1/y, -x/y**2]\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "xi = np.array([1,1]).reshape(-1,1)  # initialize\ng = lambda x: x[0]/x[1] - 1/2000.\ngrad = lambda x: np.array([1/x[1], -x[0]/x[1]**2]).reshape(-1,1)\nr_f(xi, mu, sd, g, grad)",
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Optimizing...\nbeta: 3.22, \tx:47.97, \ty:48.97, \t diff=18.69\nbeta: 7.14, \tx:15.74, \ty:65.01, \t diff=3.92\nbeta: 10.54, \tx:-1.47, \ty:58.82, \t diff=3.40\nbeta: 9.96, \tx:0.16, \ty:53.52, \t diff=0.58\nbeta: 9.99, \tx:0.03, \ty:54.06, \t diff=0.03\nbeta: 9.99, \tx:0.03, \ty:54.01, \t diff=0.00\nbeta: 9.99, \tx:0.03, \ty:54.01, \t diff=0.00\ndone!\nProbability CDF(-beta): 0.00\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "This is a little less obvious, given that we are shrinking closer to $38/54\\approx0.7$ for the max. likelihood estimate for G, and the 1/2000 doesn't really change anything. Let's check our approximations with a naive monte-carlo sampling over the variable space (a million should do...)\n\n## Check with MC"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "x_mc = norm(mu[0], sd[0]).rvs(int(1e6))\ny_mc = norm(mu[1], sd[1]).rvs(int(1e6))\n\n\ng1 = x_mc*y_mc - 2000.\ng2 = x_mc/y_mc - 1/2000.\n# g2 = (x_mc-y_mc)/2000.\n\nprint(f'Part A:\\n\\tNaive Monte Carlo, {x_mc.shape[0]:.0e} samples: P = {(g1<0).sum()/g1.shape[0]:.2f}')\nprint(f'Part B:\\n\\tNaive Monte Carlo, {x_mc.shape[0]:.0e} samples: P = {(g2<0).sum()/g2.shape[0]:.2f}')",
      "execution_count": 18,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Part A:\n\tNaive Monte Carlo, 1e+06 samples: P = 0.42\nPart B:\n\tNaive Monte Carlo, 1e+06 samples: P = 0.00\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "It's a match!"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "# Q3\n\nHere we implement a kalman filter in 1D with 5 observations. We will assume no inputs (i.e. linear quadratic estimation without controls)\n"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "a = 1.05\nb = 0.1\nc = 10.\n\nQ = 1e-6\nR = 1\nv = norm(0,Q)  # process noise\nn = norm(0,R)  # meas. noise\nP = 1.  # 1d prob => 1d cov (init)\n\n#observations and init\nxk = 0.5  # x0\ny_meas = np.array([\n    7.479527812,\n    9.015371756,\n    13.76472827,\n    14.10629548,\n    17.28189501,\n])\n\nx_est = []\nerr = []\nfor y in y_meas:\n    #predict:\n    xp = xk   # no input/controls\n    Pp = a*P*a + Q  # Predict the cov\n    \n    #correct:\n    K = Pp/(Pp+R)  # kalman gain\n    xk = xp + K*(y - xp)  # Correct x from observation\n    P = (1-K)*Pp  # correct cov from gain\n    x_est += [xk]  # store state\n    err += [P]  # sore uncertainty\nx_est, err = np.array(x_est), np.array(err)  # make array\n\nplt.figure(figsize=(12,7))\nplt.plot(y_meas, 'k.:', label='measured output')\nplt.plot(x_est, 'o-', color='dodgerblue', label='estimated state')\nplt.fill_between(range(5), x_est-err, x_est+err, alpha=.2, \n                 color='xkcd:rust', label='estimated uncertainty')\nplt.legend(loc=0)",
      "execution_count": 19,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 19,
          "data": {
            "text/plain": "<matplotlib.legend.Legend at 0x1230d4adbe0>"
          },
          "metadata": {}
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "<matplotlib.figure.Figure at 0x1230d765240>",
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAGfCAYAAAC6KN9TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4lNX9/vH3MzNJJnsgEBYjBCwggYQQCUtVCEQBC1oR\ncKMKUhewWJWfVK21tmqLVSx1oSpfF5TibkWq1I3KomyyVwQBISB7NpJJMllm5vz+GJgKBEhCyCTk\nfl2XF3lmnuWewCUfTj7nHMsYg4iIiIhIU2QLdgARERERkWBRMSwiIiIiTZaKYRERERFpslQMi4iI\niEiTpWJYRERERJosFcMiIiIi0mSpGBYRERGRJkvFsIiIiIg0WSqGRURERKTJctTnw1q0aGGSkpLq\n85EiIiIi0gStXr061xjT8lTn1WsxnJSUxKpVq+rzkSIiIiLSBFmWtbM656lNQkRERESaLBXDIiIi\nItJkqRgWERERkSarXnuGq1JZWcnu3bspKysLdhRpxJxOJ4mJiYSEhAQ7ioiIiDQiQS+Gd+/eTXR0\nNElJSViWFew40ggZY8jLy2P37t106NAh2HFERESkEQl6m0RZWRnx8fEqhKXWLMsiPj5eP10QERGR\nGgt6MQyoEJbTpj9DIiIiUhsNohgWEREREQkGFcONXGZmZp1sZLJw4UKWLl16Wvf485//fNo5RERE\nROqTiuEGyOPx1PszVQyLiIhIU9Qoi+Fly5YxdepUli1bdtr3ys7O5vzzz+fmm2+me/fujBkzhs8/\n/5wLL7yQTp06sXLlSgBKSkoYP348GRkZ9OzZkw8++CBw/cUXX0x6ejrp6emBgnLfvn3079+ftLQ0\nunfvzpIlSwCIiooKPPvdd99l3LhxAIwbN47JkyczcOBA7r333hM+z+12c+2115Kamso111yD2+2u\n8nMtWLCAnj17kpKSwvjx4ykvLwf8W2Ln5uYCsGrVKjIzM8nOzub5559n+vTppKWlsWTJEsaNG8eE\nCRO4+OKL6dy5Mx9++CEAs2bNYtKkSYHnDB8+nIULF3LffffhdrtJS0tjzJgxp/37IiIiIlIfGlwx\nnJmZyaxZswD/GsSZmZn84x//AKC0tJSePXuSmZnJgw8+SFZWFj179uSf//wnALm5uWRmZvKvf/0L\ngP3791frmdu2bePOO+9kw4YNbN68mddff50vv/ySadOmBUY7//SnPzFo0CC+/vprvvjiC6ZMmUJJ\nSQkJCQl89tlnrFmzhrfeeotf//rXALz++usMGTKEdevWsX79etLS0k6ZY8uWLXz++ec8+eSTJ3ze\nc889R0REBBs2bOCBBx5g9erVx92nrKyMcePG8dZbb/Hf//4Xj8fDc889d8LnJiUlMWHCBO6++27W\nrVvHxRdfDPgL/UWLFvHRRx8xYcKEk67W8NhjjxEeHs66deuYM2fOKT+riIiISEMQ9HWGa6qwsBCP\nx4PP56OiooLCwsLTvmeHDh1ISUkBoFu3bmRlZWFZFikpKWRnZwPw6aefMm/ePKZNmwb4C85du3bR\ntm1bJk2axLp167Db7WzZsgWAjIwMxo8fT2VlJVdeeWW1iuHRo0djt9tP+rzFixcHCu7U1FRSU1OP\nu893331Hhw4d6Ny5MwBjx45lxowZ3HXXXTX6vlx99dXYbDY6depEx44d2bx5c42uFxEREWnoGlwx\nvHDhwsDXISEhRx1HREQwZ84csrKyqKioIDQ0lDlz5tCvXz8AWrRocdT5rVu3rtYzw8LCAl/bbLbA\nsc1mC/TvGmN477336NKly1HX/uEPf6BVq1asX78en8+H0+kEoH///ixevJiPPvqIG264gSlTpnDj\njTcetQTYsSOtkZGRga9P9Dw49TJixpgTvudwOPD5fFU+/1TPsSzrqOurcw8RERGRhqzBtUmcSr9+\n/ViwYAGPPPIICxYsCBTCZ9qQIUN45plnAoXm2rVrAf9IdZs2bbDZbMyePRuv1wvAzp07SUhI4JZb\nbuGXv/wla9asAaBVq1Zs2rQJn8/H+++/X+Pn9e/fP9CG8M0337Bhw4bjrj3//PPJzs5m27ZtAMye\nPZsBAwYA/paII60V7733XuCa6OhoXC7XUfd555138Pl8fP/992zfvp0uXbqQlJTEunXr8Pl8/PDD\nD4GeavD/46WysrJa308RERE5u9XlHK8zqcGNDFdHv3796q0IPuLBBx/krrvuIjU1FWMMSUlJfPjh\nh9x+++2MHDmSd955h4EDBwZGdxcuXMgTTzxBSEgIUVFRvPbaa4C/t3b48OGce+65dO/eneLi4ho9\nb+LEidx0002kpqaSlpZG7969j7vW6XTyyiuvMHr0aDweDxkZGUyYMAGAhx56iF/+8pf8+c9/pk+f\nPoFrLr/8ckaNGsUHH3zAM888A0CXLl0YMGAABw4c4Pnnn8fpdHLhhRcG2kq6d+9Oenp64B633nor\nqamppKenq29YRESkCVu2bNlRP8mvzwHMmrJO9iP1utarVy9z7Jq4mzZtomvXrvWWQapn3LhxDB8+\nnFGjRgU7SrXpz5KIiEhw/etf/+KFF16gX79+PPTQQ3i9Xux2O4888gj3339/vWaxLGu1MabXqc5r\ndG0SIiIiItJwGGOOmo+0d+9eUlJSCA0NxW63ExoaSmZmZnBDnkSjbJOQM+/I8nYiIiIiJ1JQUMCw\nYcMYN24ct956K6NGjWLkyJHYbDYWLFjAwoULyczMbLAtEqBiWERERERq6NChQ8TFxREXF0diYiIx\nMTGAf+WpI6tRBWOOV22oTUJEREREqm3q1Kmcf/75uFwuLMvi7bff5tprrw12rFrTyLCIiIiInNTu\n3buJjo4mNjaWrKws3G73Kfc9aCw0MiwiIiIiJ3Tw4EE6d+7M448/DkDv3r15+OGHiYqKCnKyuqFi\nuIZmzZrF3r17A8c333wz33777WnfNzs7m9dff73G140bN4533323WucuXLiQpUuX1tl5IiIicnYq\nLi7ms88+AyAhIYHp06dzyy23BDnVmdHoiuG5m+GnL0PSU/5f526u3+cfWwy/+OKLJCcnn/Z9a1sM\n14SKYREREamO+++/n8svv5zc3FwAbrvtNpKSkoIb6gxpVMXw3M1w3wLY4wKD/9f7Fpx+QfyPf/yD\n3r17k5aWxm233YbX68Xr9TJu3Di6d+9OSkoK06dP591332XVqlWMGTOGtLQ03G43mZmZHNlIJCoq\ninvvvZcLLriASy65hJUrV5KZmUnHjh2ZN28e4C96L774YtLT00lPTw8Unffddx9LliwhLS2N6dOn\n4/V6mTJlChkZGaSmpvLCCy8A/rX8Jk2aRHJyMsOGDePgwYNVfqann36a5ORkUlNTufbaa8nOzub5\n559n+vTppKWlsWTJEv71r3/Rp08fevbsySWXXMKBAweqPC8nJ4eRI0eSkZFBRkYGX3311el9w0VE\nRKRB8Xq9vPrqq3z//fcA3HvvvSxcuJAWLVoEOdmZ16Am0P1xEXybc+L31+yHCu/Rr7k9MOVzeOOb\nqq9JbgkPDTjxPTdt2sRbb73FV199RUhICLfffjtz5syhW7du7Nmzh2++8d/4yBIizz77LNOmTaNX\nr+M3NCkpKSEzM5O//OUvjBgxgt/97nd89tlnfPvtt4wdO5YrrriChIQEPvvsM5xOJ1u3buW6665j\n1apVPPbYY0ybNo0PP/wQgJkzZxIbG8vXX39NeXk5F154IYMHD2bt2rV89913/Pe//+XAgQMkJycz\nfvz447I89thj7Nixg7CwsED2CRMmEBUVxT333AP41wZcvnw5lmXx4osv8vjjj/Pkk08ed97111/P\n3XffzUUXXcSuXbsYMmQImzZtOvE3VURERBqVnJwcJk6cyOTJk3n00UdJTEwkMTEx2LHqRYMqhk/l\n2EL4VK9Xx4IFC1i9ejUZGRkAuN1uEhISuPzyy9m+fTt33HEHw4YNY/Dgwae8V2hoKEOHDgUgJSWF\nsLAwQkJCSElJITs7G4DKykomTZrEunXrsNvtbNmypcp7ffrpp2zYsCHQD1xYWMjWrVtZvHgx1113\nHXa7nbZt2zJo0KAqr09NTWXMmDFceeWVXHnllVWes3v3bq655hr27dtHRUUFHTp0qPK8zz///Ki+\n6KKiIlwuF9HR0af8noiIiEjDtGjRIhYsWMDDDz9M69at+frrr+uk9bOxaVDF8MlGcMHfI7zHdfzr\n50TDW6Nq90xjDGPHjmXq1KnHvbd+/Xo++eQTZsyYwdtvv83LL7980nuFhIQElhmx2WyEhYUFvvZ4\nPABMnz6dVq1asX79enw+H06n84S5nnnmGYYMGXLU6/Pnz6/WUiYfffQRixcvZt68eTzyyCNs3Ljx\nuHPuuOMOJk+ezBVXXMHChQv5wx/+UOW9fD4fy5YtIzw8/JTPFRERkcZh0aJFzJo1i8mTJxMXF0e3\nbt2CHSkoGlXP8G9+CuHHlO/hDv/rtZWVlcW7774b6L3Nz89n586d5Obm4vP5GDlyJI888ghr1qwB\nIDo6Gperioq8mgoLC2nTpg02m43Zs2fj9XqrvO+QIUN47rnnqKysBGDLli2UlJTQv39/3nzzTbxe\nL/v27eOLL7447hk+n48ffviBgQMH8vjjj3Po0CGKi4uPe0ZhYSHnnHMOAK+++mrg9WPPGzx4MM8+\n+2zgeN26dbX+/CIiIhIcBw4cYPTo0SxZsgSAe+65hy1bthAXFxfkZMHVqIrhK8+Hx7L8I8EW/l8f\ny/K/XlvJyck8+uijDB48mNTUVC699FL27dvHnj17yMzMJC0tjXHjxgVGjseNG8eECRMCE+hq6vbb\nb+fVV1+lb9++bNmyhcjISMDf1uBwOOjRowfTp0/n5ptvJjk5mfT0dLp3785tt92Gx+NhxIgRdOrU\niZSUFCZOnMiAAccPp3u9Xn7xi1+QkpJCz549ufvuu4mLi+Pyyy/n/fffD0yM+8Mf/sDo0aO5+OKL\nj2qQP/a8p59+mlWrVpGamkpycjLPP/98Lb/bIiIiUt98Ph/gH+zasGED27dvByAiIuKEP6FuSixj\nTL09rFevXubIygtHbNq0ia5du9ZbBjl76c+SiIjI0R5//HHmzZvHkiVLsCwLr9eL3W4Pdqx6YVnW\namPM8SseHKNRjQyLiIiIyMm5XK5AG2abNm3o3LkzpaWlAE2mEK4JFcMiIiIiZ4mtW7fSsWNH3njj\nDQBuuOEGXn755UBbphxPxbCIiIhII+bxeNi6dSsA5513HqNHj26SS6TVVoNaWk1EREREauamm25i\n4cKFbN26FafTyd///vdgR2pUNDIsIiIi0sj85z//obi4GIBf/epXPPXUU4H9DaRmVAyLiIiINCLf\nfPMNWVlZzJw5E4C+ffty1VVXVWtTLjleg2uTyN24gsriwjq7X0hULC269amz+82aNYvBgwfTtm1b\nAG6++WYmT5582r052dnZLF26lOuvv75G140bN47hw4czalQtt+CrI9XNv3fvXn79618Htpk+kT//\n+c/89re/rcuIIiIijdamTZvYsGED11xzDd27d2fu3LkMHTo02LHOCg1uZLiyuJCwuBZ19l9dFtbg\nL4b37t0bOH7xxRfrpEk9Ozub119//bTvEwwej6fa+du2bXvKQhj8xbCIiIj4Pfzww9x9992Ul5cD\n8POf/1xtEXWkwRXDwfCPf/yD3r17k5aWxm233YbX68Xr9TJu3Di6d+9OSkoK06dP591332XVqlWM\nGTMmsANdZmYmRzYSiYqK4t577+WCCy7gkksuYeXKlWRmZtKxY0fmzZsH+Iveiy++mPT0dNLT01m6\ndCkA9913H0uWLCEtLY3p06fj9XqZMmUKGRkZpKam8sILLwBgjGHSpEkkJyczbNiwwDbSx/pxrtzc\nXJKSkgB/MX/VVVcxdOhQOnXqxG9+85vANR9//DHp6en06NGDrKwsAEpKShg/fjwZGRn07NmTDz74\nIHCf0aNHc/nllzN48ODj8p/oc2ZnZ9O9e/eTZrnvvvtwu92kpaUxZswYHnzwQZ566qlAzgceeICn\nn376NH/XRUREGq7c3FwmT57M7t27AZg2bRobNmxQAXwmGGPq7b8LLrjAHOvbb7896njv8k9M3ubV\ndfbf3uWfHPfMY58/fPhwU1FRYYwxZuLEiebVV181q1atMpdcckngvIKCAmOMMQMGDDBff/114PUf\nHwNm/vz5xhhjrrzySnPppZeaiooKs27dOtOjRw9jjDElJSXG7XYbY4zZsmWLOfI9+eKLL8ywYcMC\n933hhRfMI488YowxpqyszFxwwQVm+/bt5r333jOXXHKJ8Xg8Zs+ePSY2Nta88847x32uH+fKyckx\n7du3N8YY88orr5gOHTqYQ4cOGbfbbdq1a2d27dplDh48aBITE8327duNMcbk5eUZY4y5//77zezZ\nswPfg06dOpni4mLzyiuvmHPOOSdw3rH5T/Q5d+zYYbp163bSLMYYExkZGbjXjh07TM+ePY0xxni9\nXtOxY0eTm5tb5e+liIjI2WDHjh0mPDzcvPrqq8GO0mgBq0w16tMG1zNc3xYsWMDq1avJyMgAwO12\nk5CQwOWXX8727du54447GDZsGIMHDz7lvUJDQwP9OykpKYSFhRESEkJKSgrZ2dkAVFZWMmnSJNat\nW4fdbmfLli1V3uvTTz9lw4YNgZaCwsJCtm7dyuLFi7nuuuuw2+20bduWQYMG1fgzZ2VlERsbC0By\ncjI7d+6koKCA/v3706FDBwCaN28eyDFv3jymTZsGQFlZGbt27QLg0ksvDZx3rOp+zqqynHvuuUed\nk5SURHx8PGvXruXAgQP07NmT+Pj4Gn9uERGRhuzvf/87W7Zs4W9/+xtJSUn88MMP+vuuHjT5YtgY\nw9ixY5k6depx761fv55PPvmEGTNm8Pbbb/Pyyy+f9F4hISGBmZw2my3wowybzYbH4wFg+vTptGrV\nivXr1+Pz+XA6nSfM9cwzzzBkyJCjXp8/f361Zos6HA58Ph/gL2B/7Mc/YrHb7Xg8HowxVd7XGMN7\n771Hly5djnp9xYoVJ93Nprqfs6osVbn55puZNWsW+/fvZ/z48Sd8roiISGPi8/mw2fxdqzt37mTz\n5s14PB4cDocK4Xpyyp5hy7JetizroGVZ3xzz+h2WZX1nWdZGy7IeP3MRz6ysrCzefffdQO9tfn4+\nO3fuJDc3F5/Px8iRI3nkkUdYs2YNANHR0bhcrlo/r7CwkDZt2mCz2Zg9e3Zg7/Bj7ztkyBCee+45\nKisrAdiyZQslJSX079+fN998E6/Xy759+/jiiy+qfE5SUhKrV68GqNaEtX79+rFo0SJ27NgR+D4c\nyfHMM8/g/2kDrF27tsrrj81/os9ZXSEhIYHPDjBixAg+/vhjvv766+P+gSAiItIYffPNN3Tr1i3w\n9/Wf/vQnPv74YxyOJj9WWa+q892eBTwLvHbkBcuyBgI/B1KNMeWWZSXUVaCQqFjKD+XW1e0IiYo9\n6fvJyck8+uijDB48GJ/PR0hICDNmzCA8PJybbropMLp6ZOR43LhxTJgwgfDwcJYtW1bjPLfffjsj\nR47knXfeYeDAgYHR1dTUVBwOBz169GDcuHHceeedZGdnk56ejjGGli1bMnfuXEaMGMF//vMfUlJS\n6Ny5MwMGDKjyOffccw9XX301s2fPrlYrRcuWLZk5cyZXXXUVPp+PhIQEPvvsMx588EHuuusuUlNT\nMcaQlJTEhx9+eNz1x+Y/0eesrltvvZXU1FTS09OZM2cOoaGhDBw4kLi4OOx2e43uJSIi0lAYY3C5\nXMTExNCuXTtatWqF2+0GUBEcJNaREb+TnmRZScCHxpjuh4/fBmYaYz6vycN69epljqxwcMSmTZvo\n2rVrTW4jTZDP5yM9PZ133nmHTp06VXmO/iyJiEhDN2bMGH744QcWLVqkTTLOMMuyVhtjep3qvNou\nrdYZuNiyrBWWZS2yLCujlvcROaVvv/2Wn/zkJ2RlZZ2wEBYREWmotm7dGmg3HDp0KKNGjQr85FmC\nr7bj8Q6gGdAXyADetiyro6limNmyrFuBWwHatWtX25zShCUnJ7N9+/ZgxxAREamxxYsXk5mZyXvv\nvceIESO44YYbgh1JjlHbkeHdwD8PL+O2EvABLao60Rgz0xjTyxjTq2XLllXerDqtGiInoz9DIiLS\nUBw4cICvv/4agJ/+9Kf86U9/on///kFOJSdS25HhucAgYKFlWZ2BUKBWs96cTid5eXnEx8erd0Zq\nxRhDXl7eCZdvExERqU8jR44kNzeXb7/9FofDwf333x/sSHISpyyGLct6A8gEWliWtRt4CHgZePnw\ncmsVwNiqWiSqIzExkd27d5OTk1Oby0UA/z+qEhMTgx1DRESaILfbzUsvvcTNN9+M0+nkb3/7GzEx\nMYH1g6VhO2UxbIy57gRv/aIuAoSEhAR2PRMRERFpbJYvX84dd9xB69atGTVqFL16nXIBA2lAtKCd\niIiISA298847uFwuxo8fz8CBA1m9ejXp6enBjiW1oPF7ERERkRp67bXXeOWVVwITuFUIN14qhkVE\nREROYePGjQwbNoyDBw8C8Oqrr7Jw4UJN/j8LqBgWEREROYEjI78Oh4P169fz3XffAdC8eXPsdnsw\no0kdUc+wiIiIyDGMMdxyyy2EhYUxY8YMunTpQnZ2Ng6HSqezjUaGRURERA4rLi4GwLIsmjdvTlxc\n3FGjw3L20e+qiIiICPDxxx9z9dVXs2zZMrp168bjjz8e7EhSDzQyLCIiIk1WaWkpe/bsASAjI4OR\nI0cSGRkZ5FRSnzQyLCIiIk2SMYa+ffvStm1bPv74Y+Lj43nllVeCHUvqmYphERERaTJ8Ph+ff/45\nl156KZZl8cADD9CmTZtgx5IgUpuEiIiINBlz5sxhyJAhLFq0CIBrrrmG/v37BzmVBJNGhkVEROSs\n9tVXX+H1eunfvz/XXHMNYWFhKoAlwDqyXEh96NWrl1m1alW9PU9ERESaNp/PR3JyMueccw4LFiwI\ndhypR5ZlrTbG9DrVeWqTEBERkbPKjh07mDx5MpWVldhsNt5//33mzZsX7FjSQKkYFhERkbPKN998\nw/PPP8/atWsB6Nq1q5ZLkxNSMSwiIiKNmsfj4be//S0vvPACAMOHDyc7O5vevXsHOZk0BppAJyIi\nIo2SMQbLsrDb7Xz99deUlJQA/q2UExISgpxOGguNDIuIiEij8+9//5uUlBQOHTqEZVnMnz+fp556\nKtixpBFSMSwiIiKNgs/nw+12A9C2bVvi4+PJy8sDICQkJJjRpBFTMSwiIiINXkVFBRkZGfzud78D\noEePHixatIjzzjsvyMmksVMxLCIiIg3W9u3bAQgNDWX48OH07ds3yInkbKNiWERERBqkGTNm0KVL\nF77//nsA/vjHPzJ69Oggp5KzjVaTEBERkQZj27Zt2O12OnTowFVXXUVJSQlt2rQJdiw5i2k7ZhER\nEWkQ3G43iYmJDB48mDfeeCPYcaSR03bMIiIi0uAVFBTw0ksvARAeHs4//vEPpk+fHuRU0pSoGBYR\nEZGgeemll7jlllv47rvvALjsssto3bp1kFNJU6JiWEREROpNZWUlM2bMYNGiRQDcfvvtrF+/ni5d\nugQ5mTRVKoZFRESk3vh8Pp544gnefvttACIiIkhJSQlyKmnKVAyLiIjIGfWf//yHa6+9Fq/XS1hY\nGMuXL+fZZ58NdiwRQMWwiIiInCFHVqzKyclhzZo17N69G4DWrVtjWVYwo4kEqBgWERGROlVYWMhl\nl10WWCVi9OjRbNy4kfbt2wc5mcjxVAyLiIhInSgtLQUgJiYGu90eeN1msxESEhKsWCInpWJYRERE\nTtvTTz9N586dKS4uxrIsPvzwQ26++eZgxxI5JRXDIiIiUit5eXm4XC4AevfuzYgRI/B4PEFOJVIz\nKoZFRESkxvLy8vjJT37C448/DkDfvn155plniIuLC3IykZpRMSwiIiLVUl5ezldffQVAfHw8Dz74\nINdcc02QU4mcHhXDIiIiUi33338/WVlZHDhwAIDJkyfTvXv3IKcSOT0qhkVERKRKxhg++OADduzY\nAcCdd97JvHnzSEhICHIykbqjYlhERESqdPDgQa677jpmzJgBQPv27Rk8eLA2zJCziiPYAURERKTh\nWLNmDZ9++in33XcfrVq1YtGiRfTs2TPYsUTOGI0Mi4iISMAHH3zAk08+SV5eHgAZGRk4HBo7k7OX\nimEREZEmatmyZTz44INcccUVgVUipkyZwrZt24iPjw9yOpH6ccp/6lmW9TIwHDhojOl+zHv3AE8A\nLY0xuWcmooiIiNREfn4+ZWVltG3bFoBFixbh8XjIysoC4G9/+xt79uxhxowZVFRU4PP56Ny5Mxde\neCFRUVHBjC5S76ozMjwLGHrsi5ZlnQtcCuyq40wiIiJNjs/nC3y9a9cuVq9eHTj+/PPPeemllwLH\nTz/9NHfddVfg+NZbb2Xw4MGB4+uuu46rrroqcPzwww/z0EMPBY7//e9/M3/+fCoqKvB6vdhsNo0E\nS5N1ymLYGLMYyK/irenAbwBT16FEREQak4qKCvbv3x8oaLdt28Y///lPjPH/FblgwQJ+85vfBM6f\nOXMmgwYNChzfc889tGzZMnA8depULrvsssDxm2++ye9///vA8Q8//MDmzZsDxz169KBfv36B47vv\nvpvf/va3geMXX3yRN954I3D8ySef8OKLLxIaGordbic0NJTMzMzT+RaINFq16hm2LOsKYI8xZn0d\n5xEREal3LpeL9evX43a7Afj+++957rnnOHToEABffvklN9xwQ2BS2ZtvvknHjh3JzfV3CM6YMYM2\nbdpQVFQEwD//+U9GjhxJaWkpACtWrGDGjBmUl5cDYLfbcTgcgWJ50KBBR4303n777UcVr9OnT+e7\n774LHD/xxBN8/PHHgeNf/epX/PGPfwwcDx06lCuuuCJw3KFDB84999yjPnO/fv1YsGABjzzyCAsW\nLDiqmBZpSmpcDFuWFQE8APz+VOcePv9Wy7JWWZa1Kicnp6aPExEROY4xhuLi4kBx6XK5+PTTTwM7\no+3evZtPWbpbAAAgAElEQVQ//OEPbN26FYANGzYwePBg1q/3j+F8/vnnREZGsmLFCsA/cpuWlsam\nTZsAWLt2Lbfffju7dvk7AXNycvjqq68CxW6rVq248MILA3mysrL4+9//TmhoKAA33HAD69atIyws\nDPDv3FZSUhI4/uUvf8mnn34aWK/3Zz/7GQ8++GDgfikpKYH+XoDo6Ogz0svbr18/7r//fhXC0qTV\nZmT4PKADsN6yrGwgEVhjWVbrqk42xsw0xvQyxvT68Y+ARESk6fL5fGzfvj0wslpeXs7s2bPZuHEj\nAIWFhUycOJGFCxcCsGfPHnr06MH7778PwHfffUd0dHTgODs7myFDhrBkyRLAv1nEH//4x6NGU10u\nV6B4bt++PRMnTgy0JvTp04d3332XpKQkwF+c7tu3j+TkZABGjBjB9u3b6dChAwADBw5k9uzZtGjR\nAoDU1FQmTpxIREQEAG3atKFHjx6BJcm0SYVIw1XjYtgY819jTIIxJskYkwTsBtKNMfvrPJ2IiDQY\nPp+PioqKwPHKlSuPKjafeuopPvvss8Dx1VdfzaxZswLXtm7dmqlTpwJQWVnJeeedx8yZMwPHN954\nI/Pnzwf8I7/vvfce27dvByAiIoIOHToERkfbtGnDX/7yF1JTUwE477zzWLJkSaAPNy0tDY/Hw/Dh\nwwF/sbps2TJ69+4NQKdOnZg2bRodO3YM3G/kyJE0b9488LzWrVtrfV2RJqA6S6u9AWQCLSzL2g08\nZIx56eRXiYhIQ5Ofn095eTlt2rQB4IsvvgD8o5wA06ZNIyYmhltvvRWAa665hg4dOvDYY48B/oJz\nwIABgQL3yiuvZNiwYfzf//0fAH/6058YOXIkl156KeBvVSgoKADAZrMxatQounf3r9AZFhbGq6++\nygUXXABAZGQkW7dupVWrVgDExcVx8ODBQPZmzZoxd+7cwHFsbOxRE9IiIiK46KKLAsc2m5bRF5Hq\nsY4079eHXr16mVWrVtXb80REGrulS5eyaNEiMjMzadu2LXl5eaSnpwPw8ccfs2/fPm666SYAnnzy\nSfbu3cuTTz4JwNixY8nNzeWjjz4CIDMzE5/Px+LFiwG46KKLCAsLY8GCBYH3W7duzZtvvgnAhAkT\naNeuXWBVgmeffZZzzz2Xn//85wAsWbKEhIQEunTpAkBxcTEREREqREWkQbAsa7UxptepztPPf0RE\nGqgZM2Zwxx13YLPZCA0NZdCgQaxbt47du3cDMGfOHL788stAMbxr1y527NgRuD49PT0w4QvgN7/5\nDT8eAHnttdcCE76AQH/uEc8///xRx5MmTTrq+OKLLz7qWJs1iEhjpJFhEZEG6tprr+Wtt94C/Etx\n3X777YwYMSLQ1lBUVITdbicyMjKYMUVEGqTqjgzrZ1kiIg1ISUkJK1euBODOO+8kPDw8sCnCdddd\nFyiEAWJiYlQIi4icJrVJiIg0IBMnTuTDDz8kOzs7sCnCwoULyczM1FqwIiJngNokREQaAGMMlmWx\nc+dOtmzZEliRQUREakcT6EREGolHH32UPXv28Nxzz9G+fXvat28f7EgiIk2GeoZFRIKstLSUkpIS\nPB5PsKOIiDQ5GhkWEQmCjRs3YrPZ6Nq1K48++iiWZWnLXhGRIFAxLCJSzzweD1deeSVt27Zl0aJF\n2qRCRCSIVAyLiNQTr9eLzWbD4XDw+uuv07Zt22BHEhFp8jQcISJSD1wuF0OHDuXvf/87ABkZGZxz\nzjlBTiUiIiqGRUTqQWRkJM2aNdMmGSIiDYzaJEREzqC33nqLSy+9lObNm/PWW29pkpyISAOjkWER\nkTMkOzubG2+8kb/+9a8AKoRFRBogjQyLiNSxsrIynE4nSUlJLFy4kIyMjGBHEhGRE9DIsIhIHfrm\nm2/o3LkzX3zxBQD9+vXD4dC4g4hIQ6ViWESkDrVv3560tDRatGgR7CgiIlINKoZFRE6T2+1m2rRp\neDweoqOjmTdvHikpKcGOJSIi1aBiWETkNH388cdMmTKF//znP8GOIiIiNaRGNhGRWioqKiImJoYR\nI0awYcMGjQaLiDRCGhkWEamF1157jZ/85Cfs2LEDQIWwiEgjpWJYRKQWLrzwQi6//HJatmwZ7Cgi\nInIaVAyLiFTT1q1beeKJJwA477zzeOmll4iKigpyKhEROR0qhkVEqumVV17hL3/5CwcOHAh2FBER\nqSMqhkVETsLn83Hw4EEAHn74YdatW0erVq2CnEpEROqKimERkZO46aabGDhwIG63G4fDQWJiYrAj\niYhIHdLSaiIiJ3HjjTfSt29fnE5nsKOIiMgZoGJYROQYb775JmVlZYwbN46srCyysrKCHUlERM4Q\ntUmIiPyIMYZXXnmF2bNnY4wJdhwRETnDNDIsIgLk5OQQGhpKbGwsb7/9NuHh4ViWFexYIiJyhmlk\nWESavLKyMvr27cutt94KQGxsLKGhoUFOJSIi9UEjwyLS5DmdTn7/+99rS2URkSZII8Mi0iSVl5cz\nadIkvvrqKwDGjh1Lenp6kFOJiEh9UzEsIk2S2+3m008/ZcmSJcGOIiIiQaQ2CRFpUtatW0dKSgpx\ncXGsWbOGqKioYEcSEZEg0siwiDQZGzduJCMjg7/97W8AKoRFREQjwyLSdCQnJ/P0009z/fXXBzuK\niIg0EBoZFpGz2o4dOxgyZAg//PADlmUxceJEYmNjgx1LREQaCBXDInJWKy8vZ/PmzWzfvj3YUURE\npAFSMSwiZx1jDF988QUA559/Plu3bmXAgAFBTiUiIg2RimEROeu8+OKLDBo0iKVLlwJoNzkRETkh\nTaATkbOGMQbLshg7diwRERH069cv2JFERKSBO+XIsGVZL1uWddCyrG9+9NoTlmVttixrg2VZ71uW\nFXdmY4qInNy//vUvMjMzKS0tJTQ0lDFjxmBZVrBjiYhIA1edNolZwNBjXvsM6G6MSQW2APfXcS4R\nkRqx2+2Ul5dTVFQU7CgiItKInLIYNsYsBvKPee1TY4zn8OFyIPEMZBMROan8/Hw+//xzAH72s5+x\ndOlSWrduHeRUIiLSmNTFBLrxwL9P9KZlWbdalrXKsqxVOTk5dfA4ERG/u+++m1GjRlFYWAiAzaY5\nwSIiUjOn9TeHZVkPAB5gzonOMcbMNMb0Msb0atmy5ek8TkQEAK/XC8Djjz/OJ598ok00RESk1mq9\nmoRlWWOB4UCWMcbUXSQRkRO75557+OGHH3jzzTdp1aoVrVq1CnYkERFpxGpVDFuWNRS4FxhgjCmt\n20giIieWkJCA1+vF5/Nht9uDHUdERBq5UxbDlmW9AWQCLSzL2g08hH/1iDDgs8NLFy03xkw4gzlF\npAlbsWIFdrudXr16MWXKFC2ZJiIideaUxbAx5roqXn7pDGQRETmO1+vlxhtvpHXr1ixatEiFsIiI\n1CntQCciDVJ5eTkhISHY7Xbmzp2r3mARETkjtA6RiDQ4LpeL/v37M3XqVAC6du1K8+bNg5xKRETO\nRiqGRaTBiYqKIj09neTk5GBHERGRs5zaJESkQTDG8NxzzzFixAjatGnDc889F+xIIiLSBGhkWEQa\nhF27djFlyhSef/75YEcREZEmRCPDIhJUhYWFxMbG0r59e1auXKnWCBERqVcaGRaRoFm3bh0dO3Zk\n3rx5AHTr1k1Lp4mISL1SMSwiQXP++edzxRVX0L1792BHERGRJkrFsIjUq8LCQu6//37KyspwOp28\n8sordOzYMdixRETkDPB5PMGOcEoqhkWkXi1dupQnn3ySr776KthRRESkjvk8lZQX5ePavY2D678k\nZ/2X+DyVwY51UppAJyL1Yt++fbRp04bLLruMbdu20a5du2BHEhGR0+TzevGUuqgoLqQsbz+VJUVg\nDJbDjt0ZgbfMjTEm2DFPSsWwiJxxM2fO5P/9v//HqlWr6NKliwphEZFGyvh8eMpKqHAVUl5wgPLC\nPDCAzcLhjCQ0tjmWZTF/bzOeXd2W/WWhtPkG7r0Qrjw/2OmrpmJYRM644cOHs23bNpKSkoIdRURE\nasAYg7eslIqSIsoP5VBekIPxejEWOJwRhMbGH7cK0Py9zXjkm3aU+ewA7C2G+xb432uIBbFVn0PX\nvXr1MqtWraq354lI8Kxdu5a33nqLqVOnark0EZFGxFteRmWpi/JDuZTl7z88Cc5gC3XicIZj2ewn\nvf6yhd3YXxZ23OvnRMPS8WcodBUsy1ptjOl1qvM0MiwiZ8T8+fOZM2cOd911F61btw52HBEROQGf\np5LKkiLKi/IpyzuAp6wUC4MtJAy7M5IQx6nLxf3uEJbnxbA8N5r9ZaFVnrPXVdfJ64ZGhkWkzlRU\nVLB3716SkpLw+XwcOnSI5s2bBzuWiIj8SGDSW1EBZQUHqCxxYYwPm8OBwxmJLaTqYvbHXJU2VuVH\nsyIvmuW5MewsdQLQIqySEo8Nt/f40WONDIvIWe/GG29k5cqVbNy4kfDwcBXCIiINgPH58LiLqSgu\n8k96K8oHnwGbDYczgrC4+FPeo9Jn8d9DEYHR342FkfiwCLd7uaB5MaPa5dI3vojzosr4976je4YB\nwh3wm5+eyU9ZeyqGRaTO3H333ezYsYPw8PBgRxERabKOmvRWkEN5wUGMz2As45/0FtP8lHM5jIHt\nxU6W50WzIi+GVflRuL12bBi6xZYy/rz99I13kRpXQojt6C6Dn7UtAODZLYdXk4iCey+0GuTkOVCb\nhIicBmMMzz//PGVlZdx9993BjiMi0mT5J70VUXYol/L8A3grK7Esqj3pDeBgWQgr86IPF8DR5Jb7\n2yXaRZTRt4WLPvFFZDQvJjrEW+1c5YfySEgfgL0arRd1TW0SIlIvFixYQEVFBXfeeSc2mza1FBGp\nD97KCjylLsoL8/yT3srdWIAtJNQ/6S3q1CVeqcfG6vwof99vXgzfF/t/qhcXUknveNfhAthF2/CK\nGuczxuAtd2N8vhpfW99UDItIje3atYvw8HBatmzJ7NmzCQsLUyEsInIG+TwePG4XFUWHjpn0FoLD\nGYEzIvKU9/D44NuiCJbnxrA8L5r/HorCYyzCbD56Nivm8nPy6BPvonO0G1stVsT0VpbjdZfi83qx\nLAiNaU5Umw7YHCG1+MT1R8WwiNRIWVkZF110Eenp6cydO1f9wSIiZ8CPJ72V5e+noijf/4ZlwxEe\nWa1Jb8bArtIwluf6+36/zo+i2OPAwnB+TCm/SDpA3xYu0uKKCbPXvG3W5/XgLSvFW1EOWDicEUS2\naU9oTHNCImKwVWNJtoagcaQUkQbD6XTy9NNP07Vr12BHERE5awQmvRUXUlaQQ8WhHHw+H5xkp7eq\n5Fc4/H2/uf6+3yObX7QNL2dw60P+vt94F81Cq9/3e1TGcrd/HWIDlt2Os3krnM0TCImIwR7mrPE9\nGwIVwyJySiUlJdx8882MHTuWoUOHcuWVVwY7kohIo+cpd+MpdVFWkEN5/kG8nkosy8IWGkZIdBxW\nNdrPyrwWawuiWJ4bw4q8aL5zRQAQ7fDQO97F+I4H6NuiiMTwCmqzGaivsgJPWQnG44XDrQ+RrZMI\njY7FER51VuwwqmJYRE7Jsiy2bNnCtm3bgh1FRKTR8lZW4PnRTm/e8jIMPuyhYdgjIgmxn7os8xrY\nXBTBirxoVuRGs+5QFBU+Gw7LR1qzEiZ12kOfeBddY0ux16JONT4vHnfJ0a0PrRtf60NNnH2fSETq\nzKJFi+jXrx8REREsX76ckJCGPQlCRKQhOWrSW/5+/6Q3C2x2Bw5nBI5qTHoD2FMaGljvd2VeNIWV\n/vKtU3QpV7fLoW+8i/RmxYQ7ar5yQ1WtD2HNEghvnoAjMgZH2Nk/L0TFsIhUaePGjQwcOJCpU6dy\n7733qhAWETmFwKQ3VyHu/P1UugowBizb4UlvzVpU6z6FFXZWHt7qeEVuDLvd/r7fhLAKBiQU0je+\niN7xLuLDPLXKeVTrAxAS04zYI60PzshqtWecTVQMi8hRjDFYlkW3bt144403uOKKK4IdSUSkQTLG\n4HGXUFlSRFn+QSoKc/D5/P8PtTvDqz3prcJnsb4gkuV5/r7fbwsjMFhE2r30indxfdJB+sYXkRRZ\nXqu+X+Pz4jm86oNl/Nn+1/oQ3eCXPjvTVAyLSMCmTZu48cYbmTNnDp07d+aaa64JdiQRkQbFU+7G\nU1Lkn/RWkIPP6wFjsDvDCYluVq1RVZ+Bra5w/8hvXjRr8qMp89mwW4aU2BJu+8k++sS76BZbQkgt\nBmnV+lAzKoZFJCAyMpLy8nIKCgqCHUVEpEHwVpQf3uktn7L8/f5JbxbYQ0KrPekN4EBZCMtz/Tu9\nrcyLJr/CPxrbIdLNiHNz6RPv4oLmLqJq0fcLVbQ+RMcR2779/1Z9aGKtDzWhYlikifN4PMydO5dR\no0bRrl071q1bp93kRKTJOjLprbywgLL8/XjcJf/b6S08EkdEVLXu46q0sepI329eDNkl/jV440Mr\n6Rvvok+LIvrEu2jlrKxVzqNbH8AeFk5Eq3aExcar9aGGVAyLNHGzZs3illtuYcmSJVx00UUqhEWk\nSTE+H5WlLipchf4VH1yHMIDNbsfujKjWTm8AlT74pjAysNXxxsJIvMbCafdyQbNiRp6bS5/4In4S\nVVa7vl9j8FWU4XGXggHLblPrQx1RMSzSRHk8HhwOBzfddBOJiYlcdNFFwY4kInLGHT3p7QAVhbmY\nwzsRO8IjCI2r3qQ3Y2B7iZMVh7c6XpUfRanXjg1Dt9hSbuqwnz4tXKTGlRBqq/lWx1B160NM+/PV\n+lDHVAyLNEFvvPEGjz32GIsXLyY2NpahQ4cGO5KIyBnjKSv9305vBTn4PJVgWdjDnNWe9AaQU+Zg\nZX5MYKvjnPJQAM6NKGNY23z6xLvIiHcRE1LzrY7hBK0PCef6Wx8iY9T6cIaoGBZpghITE2nbti0e\nT+3WqBQRaciOTHorO5RHWf5+fBXltZr0VuqxsbogKjD6u63Y34oQF+Khd7y/57dvvIu2ERW1yvnj\n1gdjDDabjbDmCYQ3b6XWh3pkGVO7ofva6NWrl1m1alW9PU9E/mfv3r0sW7aMkSNHBjuKiEid8nkq\n/X2/RQW48/fjLS3BYLCFhOJwRlR7RNXjg01FESzP84/+bjgUicfYCLX56Nms2D/xLb6ILjFubLXo\n+4WqWx/C49uo9eEMsCxrtTGm16nO08iwSBPxu9/9jrlz53LJJZcQGxsb7DgiIrXm83oP7/R2iLL8\nA1S6DgH+9XTt4RHV3unNGNhVGhZY8WFlXhTFHgcWhvNj3Pwi6SB94l2kNSvGaa/d4OHxrQ9OtT40\nMBoZFjmLGWMoKysjPDycQ4cOsXfvXpKTk4MdS0SkRvyT3oqpLHH9b9Kbz4DNhsMZji3UWa1JbwD5\nFQ6+zotmeV40y3Oj2V/m3+q4jbOcvi1c9I0vIiPeRbPQWvb9qvWhwdDIsIhw2223sWvXLj766CPi\n4uKIi4sLdiQRkWo5MunNnX+Q8oIcjNdTq0lvZV6LdQVRga2ONxdFABDt8JAR7+KmjgfoG+/i3Ija\nbXUM/jYNj7sE4/VgDIRGxxHdrgthMXFqfWgEVAyLnMUyMjJITEys9oiJiEiweCvKqSwpovxQHmUF\nB/436S00jJCoaCybvVr38RnYXBTOijz/er/rCqKo8NlwWD56xJXwq0576RtfRNfYUuy1/F9joPWh\nvAwL63DrQ6JaHxopFcMiZ5nPPvuMkJAQMjMzueWWW4IdR0SkSv+b9JaPO28/XncpAFZICA5nBCGR\n0dW+197SUJYH+n6jOVTpL286Rbm5ul2Of6vjZsWE13Kr4ypbH5ol4GyXQEhkDA5nRK3uKw2DimGR\ns4jX62Xy5MkkJCQwYMAAjQiLSIPh83rxlLqoKC48etKb4/BOb9Wc9AZQVGkP9P2uyIvmh1L/Vsct\nwyq4uGUhfVsU0TveRYuw2i8feaT1weepxLJshETFqvXhLHXKYtiyrJeB4cBBY0z3w681B94CkoBs\n4GpjTMGZiykiJ1NcXIzT6cThcPDhhx8SH1+9HZRERM6Uoya95e2nvDAX/z7H/klv1d3pDaDCZ7Gh\nIDLQ9/ttYQQ+LCLsXno1d3Ftuxz6tnDRIbJ2Wx3DkdYHN76KMsDCHhpGRMtEwuLU+nC2q87I8Czg\nWeC1H712H7DAGPOYZVn3HT6+t+7jicipuFwu+vTpw7Bhw3jiiSdo3759sCOJSBNkjMFb7vZvc1xw\nZNKb1z/pzRlOaEzzao+mGgNbi8NZkesf/V1TEEWZ147dMnSPLeGW8/bTt0UR3WJLCKnlAO2JWh/C\nmrUkNCpWrQ9NyCmLYWPMYsuyko55+edA5uGvXwUWomJYJCiio6MZOXIkgwYNCnYUEWlivOVlVJa6\nKD+Ue3jSWwVYYAsNIyQqptqT3gAOlIUcLn79fb95Ff6R2A6RZVx5Th59Wrjo1dxFVC37fqGK1ofI\nGKLbdSY0Oo6QiGi1PjRR1Vpn+HAx/OGP2iQOGWPifvR+gTGm2QmuvRW4FaBdu3YX7Ny5sw5iizRt\nXq+XqVOncv3119OxY8dgxxGRJsLnqaSypIgKV4F/0luZG4wPKyQUhzMSm6P6U5GKPTZW5UcHtjre\nUeLv+20eWnl4m+Mi+rRw0cpZWeu8xufDU1bqb30w/iI9PL4NYXHxOCKisYeE1vre0vA1mHWGjTEz\ngZng33TjTD9PpCnYt28fTz75JJZl8cADDwQ7joicpY6a9Ja3n4riQix+NOktLr7a96r0wcbCSJbn\n+vt+/1sYiddYOG0+0pu7GJGYS58WLjpFuWvd9wv+0WpPWSkYg2VZan2QU6ptMXzAsqw2xph9lmW1\nAQ7WZSgRqdq+ffto06YNiYmJbNiwgcTExGBHEpGziH8ktYQKVyHlBQcoL8w7POnNwuGMJKwGk96M\ngR0lTlYc3ultdX40JV47NgzJsaWM63CAvvFFpDYrIdRW+7GyKlsfzu2k1gepttoWw/OAscBjh3/9\noM4SiUiVVq9eTf/+/XnxxRe57rrrOPfcc4MdSUQaOWMM3rJSKkqKKC/IobzgIMbnw1jgcEbUaNIb\nQG65g5V50YHR34Pl/jaEcyPKuKxtPn3iXfSOdxETUrutjuEErQ8tzsEZF48jMkatD1Jj1Vla7Q38\nk+VaWJa1G3gIfxH8tmVZvwR2AaPPZEgRgR49ejBx4kQGDhwY7Cgi0kj5vF487mI87hLKCnKoKMrz\nr/iA8U96i46t0aQ3t8fG6oIoVuRFsyI3hq3F4QDEhXjoHe+iT/w++sS7OCei4rRyB1offD4sm0Vo\nXALOczup9UHqRLUm0NWVXr16mVWrVtXb80Qau4MHD/K73/2Ov/71r0RFRQU7jog0Mt6KcjzuYn/b\nQ2Guf6MLy8JgsIc6cTjDa1T8eg1sKoxgeV4My3OjWX8oEo+xEWrzkdasmL7xLvrEF3F+jBvbafT9\n+jyV/tFfTyUWFiGRMTjjW6v1QWqkwUygE5Ha+/bbb3njjTcYM2YMAwYMCHYcEWnAjrQ8eNzFlBfm\nU34oF2+5G4PBZndgD6vZRhf+e8IPpWGBnd6+zovG5fGXDufHlDImKYe+8UWkNSvGaa/94FrVqz60\nVeuD1AsVwyIN0JYtW+jcuTOZmZns3LmT5s2bBzuSiDQwPo8HT1kJlSUuyg/lUFGU7+/3xWB3hGJ3\nhuOIiDzh9fP3NuPZLW3ZXxZKa2cFkzrv5WdtCyioOLLVsX/0d19ZGACtneVktT5E33gXGfEumofW\nfqtjOEnrQ2QMdmeEdtGUeqM2CZEGZsaMGUyePJnVq1fTvXv3YMcRkQbCW17mb3koLqT8UA6VxUVg\n+Rd7cISFYw9zVrvlYf7eZjzyTTvKfP873275aBVWyb6yUAwWUQ4PveOL6RNfRJ94F+0iyk9rybMj\nxbuvssLf+hAVgzO+DaHRcTjCo7DZq9+uIVIdapMQaaSuvfZaiouL6dq1a7CjiEiQHGkb8Lc85Plb\nHirK/bu72e3YQ2ve8nCE22Pjyc2JRxXCAF5jI6cihImd9tE3voiuMaU4TqM11/h8eMvdeMvLALCF\nhKr1QRokjQyLNACLFy/mtddeY+bMmdg0MUSkyTmyVu6xLQ9Y/iLSHhZRo93dfiy/wsH6gkjWFESx\ntiCKzUUReE3VRbSFYc3QtbX+HMe1PsS2xNk8gdCoWLU+SL3TyLBII7J+/Xq+/PJLcnJyaNWqVbDj\niMgZ5il34yktpqL4EOUFOXhKi8ECsLA7wwmJjqvVignGwF53KGsPF75rC6IC2xyH2nx0jy3hpg77\neW93CwoqQo67vrWzZkugVdX6cGTDC7U+SGOhkWGRICktLWXHjh1069YNYwxut5uICK2XKXK2ObKr\nW2WJi4qifMoP5eCrrATAcvhXebCHhtXq3j4D24rDWVcQyZp8f/F7ZKOLaIeHtGYl9GxWTM9mxSTH\nlgZ2equqZ9hp8/Jg9138rG3BST/LkdYHYwz2kBCczdvgbNZCrQ/S4GhkWKSBu+GGG1i+fDnbtv3/\n9u48SLKzvPP9983l5L7V0quqte9CqGWFsYCWZLDssUYIkDQMESPrCi8ylwDDncEGObyEPKBgCJvL\njLhcArANMsQ1tuESAgsxY4SQZAtdC7WEERIKEFq6W7VnZVVuJ/Oc894/TlZ19abuqsyqzKr6fSI6\n1N11dPKtN05X/urUc57np6RSKQVhkU0iLHmo0qrOhw+6LZSxQRhCI4kE0WSGeHZ1b7+twPDjSnrp\nru+T5cxSq7NtidZS8N07VOOc7Il7/S4G3uN1kzhaWPpQg8AulT5kd5+t0gfZNBSGRfrkzjvv5MCB\nA6RSqX4vRURWyVqL7zbwGjVa82WalSn8Rh1rw+AYTaSI50qrHhJR9SI81Qm9T5SzPF3J4Abhuc7M\nNLl2xxyXlapcXqqyK9VaUbeH63aVjxt+jyl9yOTI7T4HJ18ils6p9EE2HYVhkXVireWuu+7CdV3+\n7BpWVisAACAASURBVM/+jEsuuUSt00Q2mMD38Zs12vUqbmWa1twMge+BtWHJQzJFoji86vNPu7Ej\n6n2fm08RYIgaywX5Ov9hzxR7S1UuK9W67vO7aKn0odnAAtF4nNRi6UM6t+oSDpGNQmFYZB09//zz\nuK5LEATqGiGyAfjtFl6jSnuhQnNuinZ1DmvBGEPESRBNZ4hHV/dWai28VE8slTs8Uc7ycj182C0Z\n9bm0UON3zh5n71CV1xRqpGNBTz4nGwT4rbDrg8FgDDiFEZU+yJalMCyyxp555hmy2SxjY2N85jOf\nIRaL6Y1GZAAdLnmo4s7N4lam8Zv18IORCLFkCqewut6+AL6F5+ZTR9z5nel0dCjGPS4rVbl5bJrL\nS1XOz9eJ9+j7ZRv4+G4z7PfbKd9w8sNkdpyBk80TTWZU+iBbmsKwyBpyXZdrr72Wyy67jG9+85vE\n48e2MhKR/gh8H69RxWvUaJanaFWmw3HGNiASd4glUsRKI6s+f9M3/KiSWbrz+1Q5S80PQ+fOpMsv\njcx3HnircUamecKH3VbKBj5es0HgNgHTeegtvPMbz+SIJTOrrmEW2YwUhkXWgLUWYwyJRIJ77rmH\n8847r99LEtny/JYbdnlYqOBWpmkvzIExWCxRJ0k8mz/lccbHM9+O8mQ5s3TX9+lKGs+GofOcbIPr\nds0udXvYkWr36tMi8D18t0HQcrE2nFCXKI6SOG2EeCZPLJlW+BV5FQrDIj1WqVR4xzvewW/91m/x\njne8gze96U39XpLIlmOtxV8aZzwbjjN2G1gskWjY23e144wXTTTjS719nyxn+Wk1icUQMwEXFerc\ncsYke0tVXlusUXD8nn1ugefhu3X8dgtjIRKL4xRHSJZGiadzqvkVWSGFYZEeS6VSeJ5HvV7v91JE\ntozA88IH3erVI8YZWyzRmEM0mSKWzqz6/NbCz2tJ9i+783uoEXZZSEd9Xlusce3OMntLVS4p1EhG\nezfQKvDa+M0GftvFmAiRWIxEaRvJYqfbQyKl8CvSBYVhkR75+te/zq/+6q+STqf5p3/6J705iawh\n3212BltUwsEW1fnD44wT3Zc8tAP4yXyaJ8qH7/zOtcO3zCGnzd5Slf90enjn99xcg1gPqxCCdgvP\nbWDbLSyGqJMgObSdRHGYWDpHLKHe5CK9pDAs0qVHH32Uv/u7v+OTn/wkd911F3fccYeCsEgPheOM\nF0seZsKSh5YLWCKxGFGn+5KHhhfhh5UM+2ez7C9n+GElQ7PzsNtYuslV2ypL9b570u6KhlucjN92\n8ZsNAq+NMRGiTpLU8C4ShaGw7CGR7N2LicgxFIZFVungwYN84Qtf4KMf/SitVotEIsEb3/jGfi9L\nZMMLxxnXaNcWjih5wEAk7hBNpIlncl29RrkV5clydunO77PzaXxrMFjOyzV42+4ZLh+qclmxymiy\nN8MtFvktF9+tE3heGH6TKdLbTiORH9KQC5E+UBgWOUXT09N8+9vf5uabbyaRSPClL32JP/qjPyIa\njeL74cMxjzzyCPv27evzSkU2Fs9t4NWrtKpzuOUpvHr1cMlDMkU8V+yqG4K18ErDWQq++8tZfl4L\n77Y6kYBLCjVuO3OCvaUqlxar5OK9GW4RvrYlaLnhw3u+D8YQS2XI7DgdJxeON47GnZ69noisnMKw\nyAksLCzw7W9/mze84Q3s3LmTf/7nf+aWW25hbGyMq666iltvvZXdu3dz++2302q1cByHa665pt/L\nFhloYclD565vZYZWZZqgHbYZM7Gwy0Oii96+AIGFn1WT7C9nl7o9TLph4MzFPC4r1bh+9wyXl6pc\nVKjjRHr3sFsYfpt4zQY2CDBAPJsnu+ss4rkCsVRW4VdkwBhre/dF4GSuuOIK+/jjj6/b64msRKvV\n4v777+fss8/m4osv5tlnn+XCCy/kc5/7HL/9279NtVrlueee47WvfS3RZdOaHn30UR588EGuueYa\nrrzyyj5+BiKDJyx5qNKqzocPui2UsUH4vhNJJIg6KSKx7u7LtALDjyvppbu+T5YzLHjhOUcTLfaW\nqlxeqrJ3qMY52UbPhlvA4al1vtuETilHPFskWdqG0wm/kZiG7Yj0gzHmB9baK056nMKwbFXWWu6/\n/36y2Sz79u2j0WhQLBZ5//vfz8c//nGstTz22GP8wi/8gibHiZyCw+OMa7TmyzQrU/iNelinEDFE\nEymiTrLrARBVL8JTndD7RDnL05UMbhCe88xMk8s6D7pdXqqyK9Xq6cNuNgjwW038ZiP8vAzEcyWS\npe042Xwn/OqHriKDQGFY5DgefvhhZmdneetb3wrAeeedx0UXXcTXv/51AJ588kkuvPBCEgk9wCJy\nMoHv4zdrYW/fyjStuRkC3wNrw5KHZIpovPt/S9NubOmu7/5ylufmUwQYosZyQb7OZcVq+LBbqcaQ\n09uH3WwQhAHfbYCFSMTg5IdIlLaHo41TWSLR1bdwE5G1ozAsAuzfv5+nnnqK2267DYC3vOUt/Oxn\nP+PHP/4xAM8++yx79uwhnU73cZUiG4PfboWDLRYqNOemaFfnwpujxhBxEkQTSSLR7u6KWgsv1RNL\n5Q5PlLO8XA8fdktGfS4t1NhbqnFZqcqlxRrpWO8edgOwgY/XDEcbYy0mEsEpDJMsbVsKvxptLLIx\nKAzLlvTcc8/xj//4j3zgAx/AGMMHP/hBPv3pT1Mul0kkErzwwgsMDw+Ty3XXlklks1saZ9ys4c7N\n4lbCccZhyUOEWDJFxEl23VPbt/DcfOpwve9clmk3LEsqxr2lkoe9pSoX5OvEe5xDA9/DdxsErVYY\nfqNRnEJntHEmRyyZUfgV2aAUhmVLOHjwIF/5ylf4zd/8TYrFIp/73Oe4/fbb+elPf8rZZ5/N+Pg4\njuMwNDTU76WKDLTA95eNM56mVZleGmccicWJJVJEetAFoekbflTJLN35faqcpdYZbrEz6XL5UBh8\nLyvVODPT7OnDbtAJv806QbuFtRCNxXGKwySKo8TTOWKpjIbmiGwSCsOyKZXLZb785S9z7bXXcv75\n5/PQQw9x9dVXc//99/Nrv/ZrVCoVms0m27dv7/dSRQaa33LDLg8LFdzKNO2FOTAGiyWWSBFNJLsa\nZ7xovh3lyXJm6c7v05U0ng3vtJ6TbSzd9d1bqrIj1e769Y4WeGH49b0WxkIkHidRHCVZGg3bnCXT\nCr8im9SphmE98ioDzXVdvvjFL3LhhReyb98+XNflfe97H3fffTfnn38+V155JQcOHGD37t0AFAoF\nCoVCn1ctMliWSh4aVdzKLO7cFL7bDO/6RsPevt2OM1400Ywv9fZ9spzlp9UkFkPMBFxUqHPLGZNc\nVqpxWbFKwfF78NkdKfDaeM06ttO7OOIkSAxtI1kcCcNvIqXwKyJHUBiWgfOlL32JdDrNjTfeSDwe\n58Mf/jC/8Ru/wb59+9ixYwcvvfQSY2NjAMTj8aUgLCKhwPOWlTwcHmdssURjDtFkilg62/XrWAs/\nryXZv+zO76FG2D0iHfV5bbHGtTvL7C1VubhQIxXt/U8ig3YrDL+ehwViiSSpoZ0kisPE0lliiVTP\nX1NENheFYem7b3zjGxw6dIjf/d3fBeDuu+9mZGSEG2+8kUgkwtNPP82OHTuWjl8MwiIS8t1mZ7DF\nHO7cNO3q/OFxxokk8Wy+JyUP7QB+Mp9eGmv8ZDnLXDt8Gxly2uwtVflPp0+yt1Tl3FyD2Bo8d+a3\n3bDm1/cxFqLJNOltp5HID3Xu/CZ7/6IisqmpZljW3UMPPcT3vvc9/viP/xiAW2+9lccee4yf/OQn\nAExNTTE8PExET3CLHCMcZ7xY8jCDOzeN33IBSyQWI5pIE4k7PSkFaHgRfljJsH82y/5yhh9WMjQ7\nD7uNpZvsLdWW6n33pN2eDrdY5LfC8GsDHyxE0xmSpW2Hw6+jnuAicnx6gE4GxlNPPcU999zDxz72\nMeLxOHfddRcf+9jHOHjwILlcjnK5TD6fP2LEsYiEAt/Hqy/Qri0sK3mwYCyRuBOG3x5NPCu3ojxZ\nzi7d+X12Po1vDQbLebnwYbfLh6pcVqwymuztcAsIa5uDlovXrEMQYAEnkyNR2oaTL4XhtwcdLURk\na1AYlr554YUX+PSnP8373vc+xsbG+OpXv8ott9zC448/zsUXX0y1WiWRSGjEscgJWGtp1+Zpzk5Q\nnziA9T0wJpzo1oNxxuFrwCsNZyn47i9n+XktLDFwIgEXF2pc3rnze2mxSi7e2+EW4RosQasZhl9r\nwUI8WyBZ2oaTL3ZGG+vrhIisjrpJyLqZnZ3lE5/4BG95y1t43eteR7Va5ZOf/CRXXXUVY2NjXH/9\n9czNzS2NOM5mu39wR2Qz8pp1muUp6uMv4bsNTCxKLJPreqobQGDhZ9Uk+8vZpW4Pk254lzUb87is\nVOP63TNcXqpyUaGOE+n9jRJrLb7bwG82gPD88WyR3Nh5ONlCJ/zqbUlE1pe+6siKeZ7HXXfdxWte\n8xre/va3k0gk+MQnPsH27dt53etex8UXX8zc3NzSiOPFECwix/LbLdzKDPXxl2hXK2BMOPksnenq\nvK3A8ONK+vBkt3KGBS/8kj+aaIUlD6Uqe4dqnJNt9Hy4BYT1zX7nzq/BABYnP0Rm++nEs/kw/Ko8\nSkT6TGFYTsndd99NNBrlPe95D7FYjHvuuYcbbriBt7/97WQyGWZmZkilwhZGxpilICwixwp8n/ZC\nmfrUIZqzE+Ggi2SaRGlk1eesehF+2Glx9kQ5y9OVDG4QllOckWnyKzvmlgLwrlRrTR52s4GP7zbx\n3SZYMBFw8sNkdpyBk80TTWYUfkVk4CgMy3F9+ctf5plnnuEjH/kIAN/61rdIJBK85z3vAeDpp58+\n4o7vYhAWkeOz1tKuVmjMTNCYOkDgeUSTSZzC0Ak7P9x3qMSnntvFeNNhR7LFe887xHW7ygBMu7Gl\nu777y1mem08RYIgaywX5OjePTYUPu5VqDDm9f9gNwvDrNRsEbhMwmIjBKYyQ3X12eHc7melJfbOI\nyFrSA3QCwP3338/f//3f8/nPfx5jDB/4wAf47ne/y/79+4lEIrRaLRxHT3GLrJTXqNGcnaQ28RJ+\nyyUSixHP5E7a9/e+QyX+64/20AwOHxc3Aa8pVJlqObxcDx92S0YCLi2GoTd82K1GOtb7h90AAt/D\ndxsELZcw/EZIFEdJFEeIZ/LEkmmFXxEZGHqATl7V448/zl/8xV/wmc98hkKhwAsvvMADDzzAzMwM\nIyMj/Pmf/zmxZQ+yKAiLnDq/5XbqgF+kXVtYqgOOZ3KnfI67n9t1RBAGaNsI++dyXL2tws1j0+wt\nVbkgXye+Rvkz8Dx8t47fbmEwRKIxnOIIydIo8XSOaDKt0cYisuEpDG8Rzz//PHfeeSe///u/zyWX\nXMLCwgLf+973eP7559m7dy+/8zu/w7vf/e6l42N6oltkRQLPo7VQpj55ELc8CQZiqcyK64B/Vk3y\njYPDjDdP/A3o/3n5890u97gCr43fbOC3XYyJEInFSAxtJ1kYJpbOEU2kFH5FZNPpKvEYY/4P4LcJ\ne+T8G/Aua22zFwuT7szNzfHBD36QG2+8keuuu45kMsl9993HTTfdxCWXXMLVV1/NwYMHl97YNPBC\nZOXCOuA5GtPjNKYPEvg+0UQKpzi8otBYaUX59niJew8O83QlQ8xYEpEANzj23+WOZKtn6w/aLTy3\ngW23gAgRxyE5tJ1EMQy/sYSeBRCRzW/VYdgYsxv4PeAia23DGPN3wDuBL/RobbIC1lp+7/d+j4sv\nvph3v/vd5HI5vvOd77B3714Adu3axeTk5NIbtEYdi6xeu16lOTtJffIlglYLE48TzxZOWge8nBfA\nozN5vnFwmAcnCrRthHNzdT54wcv8+s4y35/JHVMznIz4vPe8Q6tet9928ZsNAq+NMRGiTpL0yG6c\nfCkse0gkV31uEZGNqtufhceAlDGmDaSB1X+VlhX7yEc+QqPR4KMf/SjGGH70ox+RyYS9SaPRKM8/\n//wRd6f0402R1fPdJm5lhtr4i7TrVUwk0qkDzq/oPItlEP94aIhpN04x3uY/7Jnmht0znJ9vLB23\n2DXiRN0kTmnNLRffrRN4Xhh+kynS28dI5Eph2YOjHuAiIqsOw9bag8aYPwdeAhrA/7TW/s+erUyO\n8bnPfY6HH36Ye+65BwjHHjcah988H3jgAYVfkR4K64BnqU8cwJ2bBmOIpTIkV1gHfLwyiDeOVrhh\n9wxvHJ0nfoJpb9ftKp9y+A1HG7v4bgPr+0trzew4AydXDMNvXA/CiogcrZsyiRLwVuBMYA74e2PM\nLdbaLx113O3A7QB79uzpYqlbz7333svdd9/Nt771LWKxGLOzs7z00ku0223i8Tif//znjzhe4Vek\nezYIaFXnaEy/QmP6EDYIiCXTK64DPm4ZRLbBf7ngANftnGUo0V3v3zD8NvGaDWwQYIB4Nk9211nE\ncwXi6RyRWLyr1xAR2Qq6KZP4FeDn1topAGPM14DXA0eEYWvtZ4HPQthnuIvX2/Qee+wxPvzhD3PP\nPfcwNjZGu91mYWGBiYkJdu/ezYc+9CE+9KEP9XuZIptSu75Ac3aC+vjLYU1tPI6TK624b+6xZRAe\nNy+WQeQaq578Zq3FdxvhdLcgAAPxbJHc2Gk4i6ONFX5FRFasmzD8EvBLxpg0YZnEmwFN1FiBn//8\n57z73e/mD//wD7n66qtJpVJUKhUmJiYYGxvjpptu4qabbur3MkU2Lc9t4M7NUJt4Eb9eg0g0rANe\nYWvB1ZZBnIwNfNq1eaznh+E3VyK9bc+y8KsWiCIi3eqmZvgxY8w/AE8AHrCfzh1gOb5arcY73/lO\nbrzxRt71rncxMjLCxMQECwsLAFx66aU88cQTfV6lyOYWeO1wIMbkAdzKLMYYYunsivsBr2UZhNes\n4zdqmGiU1OgYqeHtxNI5ImqBKCLSc13dVrDW/inwpz1ay6Z06623cuaZZ3LnnXeSTqep1Wq0220A\ncrkcTz75ZJ9XKLL5LdYB1ycP0px5BawlmsqQWGEdMKxdGUTgebRr8xD4xHMlcueeQ6Iworu/IiJr\nTF9le+yOO+7gwIED/M3f/A0QtjhbfLM1xvDAAw/0c3kiW4a1Fq++QGNmnPrEAQLfIxp3cPJDK64D\nnm9Huf+VEt84OMyPKhmiy8og9nVTBmEtXqNG4DaJxOJkd50V3gVOZVZ1PhERWTmF4S596lOf4m//\n9m955JFHAEgmk6TT6aWP//Vf/3W/liayJXnNOs3yNPWJl/CbdUw0SiyTIxJd2Ze7VyuD+PWdswx3\nUQYRtFu0awtgLYnSNjJnXrSqkC4iIt1TGF6hr33ta9xxxx08/vjj5HI58vk8u3btol6vk06n+dM/\nVdWIyHrz262wDnjiZVoLZYyJEMusvA4Y1rAbRBDg1asE7RbRRJLcnvNJDW3T1DcRkT5TGD6Jxx57\njNtuu42vfOUrXHrppYyMjHDhhRdSLpfJ5XLceuut3Hrrrf1epsiWE/g+7cU64NkJrA06AzFGV3yu\ntSqDgHByXbtRxWBIje4iPbqLeLaovuAiIgNCYbjDWosxhpdffpnrr7+eP/mTP+Gmm25i165dnHnm\nmUsPvV111VVcddVVfV6tyNZkraVdm6cxM0Fj8mWs7xNxHJz8yvsBewF8v1MG8d0el0HYwKddncf6\nPrF0luJZl5AojmgCnIjIANqyYXgx/LZaLfbt28fb3vY27rjjDnbu3MnY2BiZTPgAy9jYGPfdd1+f\nVyuytXmNGs25aerjL+K5DSKx2KrqgCEsg/jmwSG+eWiIadfpWRkEhPXKXqNGJBolvX0PqZGdxNO5\n1Z9QRETW3KYPw48++igPPvggV199Na9//esBePvb306pVOKv/uqvcByHSy65hN27dwMQi8X45je/\n2c8liwidOuC5aeoTL9OuzoGJEM/kSKazKz7XicsgXu66DCLwPNrVCliLky+R33MeTm5ILdFERDaI\nTf3V+tFHH+XNb34zzWYTYwyPPPIIV155JZdeein5fH7puL/8y7/s4ypFZFHg+7QXytQmDuDOTWGx\nxFMZEquoA17TMghr8RpVfNclGouTO+0cksPbiSXTJ/+fRURkoGzqMPzggw/SarWwNrzr893vfpcr\nr7ySO++8s88rE5FF1lra1QqNmXEaUwcIfJ9oIolTGFrVQ2ZrWQbht1286gIYQ7K0jfRZp+HkVl6v\nLCIig2NTh+FrrrkGx3FotVo4jsMv//Iv93tJItLRrldxy1PUJl7Cd5tE4nHi2QImsvKRw2tZBmGD\ngHZ9gaDdJpZIkT/zIpLFEbVEExHZJDZ1GL7yyiv5zne+w4MPPsg111zDlVde2e8liWxpfssN+wGP\nvxgOnYiEdcDxzMofMlvLMggA323QrlcxxpAePY3U6K4wrKslmojIprKpwzCEgVghWKR/As+jtVCm\n3qkDxkAstbqBGLC2ZRCB7+HVFgh8HyeTo3jOpSSLI0Ri8dWfVEREBtqmD8Misv5sENCuVWhMvUJ9\n+iA2CIgmUjjF4VXdWV3LMggIW7d5jTqRWIz09jG1RBMR2UIUhkWkZ9r1BZqzU9QnXiRotTCOg5Mr\nrqoOeK3LIAKvHQ7GsJZEYYj86Rfg5IeIRFe+VhER2bgUhkWkK77bpDk3TW3iRbx6DbNYB5wtrOp8\na1kGYa3Fqy+ELdESCbVEExERhWERWbnAa9NaKFObeJnW3AwYQyydJbnKOuC1LoPwWy5evQpAcng7\n6W2n4WSLaokmIiIKwyJyamwQ0KrO0Zh+hcbUQbCWaDK96jrgtS6DsIFPu7ZA4LWJJTPkz7iQZGmU\nqJPo6rwiIrK5KAyLyAktlhU0y5PUx1/Gb7eIOgmc/NCq76quZRkEhC3RvHoNIstaomXyaokmIiLH\npTAsIsfw3AZueZra+Iv4zXqnH3CeeG51dcBrXQYR+B7t6gLW93ByRQrnvEYt0URE5JQoDIsIENYB\nu5UZ6hMv486XMZ064NX2A17zMghr8Ro1fLdBJBonu+t0kkM7iKezXZ1XRES2FoVhkS0s8H3a1Tnq\nU4doTr8CQDSVXvWDcLD2ZRDLW6IliyMUzrhQLdFERGTVFIZFthhrLe3aPM3ZCeoTBwj8dlgHXFh9\nHfBal0HYIMCrVwnaLSJOgtye80gObSOWSHV1XhEREYVhkS3Ca9Zplqeoj7+E7zYw0SixTI5IdHVf\nBpaXQTw4WaAV9LYMAhZboi0AhuTwDtLbdqslmoiI9JTCsMgm5rdbS3XA7YW5sB9wJksi3dsyiJvG\nelcGEbZEm8d6PtF0hsKZF5MojqglmoiIrAmFYZFNJvB92gvlsA54dgKLJZZMr/pBOFj7MggI71z7\njRpEIqS37SE9upNYOqeWaCIisqYUhkU2gcU64Mb0OI2pA1jfJ5Lo1AGvMkyuRxlE4Hm0a/MQBMRz\nRXLnnkOiMEIkpi9NIiKyPvSOI7KBeY0azfJU2A+45RKJdVcHDPB8Nck31rIMotMSLXCbmGiM7K4z\nSQ3vIJbKdHdiERGRVVAYFtlg/JYb1gGPv0S7VgETIZ7JEc/kVn3O9SiDCNot2rUFsJZEcZTMmRcR\nz5XUEk1ERPpKYVhkAwg8j9ZCmfrkQdy5KcASS2VIlEZXfc71KIM4tiXa+SSHRtUSTUREBobCsMiA\nstbSrlZoTL9CY/ogge8TTaS6qgOGtS+DAPDdJu1GFYMhObKTzLbdxLNFPQwnIiIDR2FYZMC061Wa\ns5PUJ18iaLUwsRjxbAETWX05wXqUQSxviRZLZymedUnYEi3udH1uERGRtaIwLDIAfLeJW5mhNv4i\n7doCJhrt1AHnV33O9SiDgLAlmteoEYlGSW3bQ3pkR1frFhERWU8KwyJ9EtYBz1KfOIBbmSGsA86S\nHFp9HTAcWwZRiHvceNo0bz2td2UQy1uiOfkS+T3n4eSG1BJNREQ2HL1ziawjGwS0axUaU69Qnz6I\nDQJiyXTXdcDHK4N4w0iFG3YfYN+2Ck4vyiCsxWtU8d0m0ZhDdtdZpIa3qyWaiIhsaArDIuugXV+g\nOTtBffxlAq+NicdxcsWu6oB9C9+fznPvwSEenCzSCiKck23wn88/wHW7elcG4bddvFoVgGRpG+mz\nTsPJlTCRSE/OLyIi0k8KwyJrxHMbuHMz1CZexKtVMdFYWAfcZSnBicogbtg9wwX53pRB2CCgXV8g\naLWIJdPkT7+AZGmUaCLZ/clFREQGiMKwSA8FXpvW/Cy1iZdpVWYBiGVyXdcBr0cZBCxviQbp0dNI\nje4KO1moJZqIiGxSCsMiXQp8n3atQn3qEM3pQ2At0VQGpzjcVYhcrzKIwPfwagsEnoeTzVM86xKS\npVEisXhPzi8iIjLIFIZFVih8kKxGuzZPc3aSVmWKILBE4w5OfqjrWtrFMoh/PDTE1BqVQQB4jRpe\ns04kGiO9fYzUyE7i6dWPdBYREdmIFIZFToHfcmnX5nHnpmnOjuO32xggmkwR78HDZPPtKN/ulEH8\n27IyiA/1uAwi8Nq0q/NYG5AoDJM//QKc/BCR6Oof5BMREdnIFIZFjiPwPLzGAm6lTGP2Ffx6DbBE\n4gmiyQzxbPf/dNarDMJai1dfwG+1iDoOudPOITm8nVgy3ZPzi4iIbGQKwyIc7qHbWqjQnB2nNT+L\ntRCJRomm0iRKIz17rfUqg1hqiWYtyeEdpLefhpMtqiWaiIjIMl2FYWNMEfg8cAlggd+01j7ai4WJ\nrDXPbeDV5mnMTuKWp7C+DyYsfXAK3T38drT1KoMIh3rME3htYskM+TMuDFuiOYmenF9ERGSz6fbO\n8H8H7rfW3myMcQD93FUGVuC1w7rfygzNmQl8t4k14YNv8WyuqwEYx7NeZRAAvtvAq9fAQHrbWNgS\nLZNXSzQREZGTWHUYNsbkgauA2wCstS2g1ZtliXRvcXDEYulDe2EOABOLEktmiKW7HyN836ESPyUx\nMQAAEh5JREFUn3puF+NNhx3JFu897xAX5BvrUgYR+B7t6gLW93CyBQrnvIZkcUQt0URERFagmzvD\nZwFTwF8bY14L/AB4v7W21pOViayQtRa/WadVm8ctT+GWJwmCAGMMsVS6676/R7vvUIn/+qM9NIPw\njvIrzQR/9MMzsJi1K4PofI5hS7Q42V2nkxzaQTyd7cn5RUREtppuwnAMuBx4n7X2MWPMfwc+DPzx\n8oOMMbcDtwPs2bOni5cTOZbfbuHV5mnOzdCcHSdouWAg4iSJ5wo9L31Y7n88t3spCC+yGHIxj/93\n3497WgZxuCWaJVFUSzQREZFe6SYMHwAOWGsf6/z5HwjD8BGstZ8FPgtwxRVX9Ob2mGxZge/j1Rdo\nzZfD0ofaAhZLJO4QS6aJZ9ZuaIS18PNakoen8jw0WWCiefxyhKoX7UkQXmyJFrRaRJwEuT3nkRza\nRiyR6vrcIiIiElp1GLbWjhtjXjbGnG+t/QnwZuDHvVuayFHT3mbGcSvTYd+SSIRYKtPTlmfH0woM\nP5jN8tBkgYenChxshF0Zzs/VycZ8qt6x/4R2JLsrnfdbLl59ATAkh7eT3qaWaCIiImul224S7wO+\n3Okk8Tzwru6XJFud7zZp1+dplqdwZycJvDYYQzSR7Mm445OZdmM8MhWG30enczT8KIlIwOuGF7jt\nzAneOFphR6p9TM0wQDLi897zDq34NW3gh3e5PY9oKk3hzItJFEfUEk1ERGSNdRWGrbVPAlf0aC2y\nRQWeR7s+T2t+lsbMOH6zATYsfYimM8SjazsbJrDw7HyKh6cKPDRZ4MfzYZeJHckW1++aZd9ohSuG\nF0hFj6zyuW5XGeCYbhKLf38qvGYdv1GDSCRsiTayUy3RRERE1pEm0Mm6s0EQTnurztOYeYX2QvnI\naW/F4TVfQ92L8NhMjoenCjw8lWfadTBYXlOs8d5zD7Jv2zznZk/eBu26XeUVhV/ohP/aAgQ+8VyR\n3LnnkCgMqyWaiIhIHygMy7rwmnW8+kJn2tsk1vexBmLJdM+nvZ3IwbqzFH7/dSZH20bIxnxePzLP\nG0cP8YbReYac3nWAWG6x9tlvNojE1BJNRERkUCgMy5o4ctrbeGfaG0SdBPFsfk1bni3yAvjhXKYT\ngAv8rBp2YTg93eQ/nj7FVaMVLitVia9hCXLQboV3ga0lURyleOZFxHMltUQTEREZEArD0hOB74el\nDwtzNGfGadfmsTYgEosTS2WIrdMd0Eoryr9M53l4qsA/T+eZb8eIGcvlQwu8/bRp3jg6z+kZd03X\nYIMAr14laLlEEkm1RBMRERlgCsOyKkvT3qoVmrOTtCpT2MCGLc+SaZzC0LqUPlgLz9eSPDwZlj88\nNZfFt4aS0+bq0Qr7tlX4peF5cvFgzdfiu028RjiAMTmyk8y23cSzRT0MJyIiMsAUhuWU+S03LH2Y\nm6Y5O47fbmOMIeIkiOdK69YH99V6/77rrHH2jVa4uFAnug4ZNPDaePVq2BItk6VwVqclWtxZ+xcX\nERGRrikMywkFnofXWMCthNPevEYNrMXE48SSGeLZ9bt8ppoxHpku8MhJev+uh6Ddol2vQhAQcRJk\ndp5BsjRKLJ3TXWAREZENRmFYloQdD6q0q/M0Z8dx52fDOgTTmfa2Di3PFq229+9a8Vsu7foC2LAD\nRm732SRKI8RSWQVgERGRDUxheIvz3AZefeHYaW/JVDjtbR2D3mLv34emCjzSRe/fXvHdBl6jDtYS\nTWfI77mARGFI7dBEREQ2EYXhLSbw2rTrC+G0t+lxvGYdsESdxLpMeztaP3v/Hs1ai+828Jt1wBDP\n5MmfcSGJwhCxZHpd1iAiIiLrS2F4k1ua9rZQ6Ux7mwPAxKLEkhmSpZF1Xc8g9P5dbrErhtesYzA4\nhSGyu84ikR8imkiuzyJERESkbxSGNyGv0/LMLU/hlicJggBjDLFUGqe4PtPelhuE3r/L2SDAa9YI\nXBeMIVEYJjd2Lk6uRNRJrNs6REREpP8UhjcBv93Cq83TnJuhOTtO0HLBQMRJEs8V1mXa23KD1Pt3\naU2BH45DbrkYY0gO7SB1xg6cXIlILL5u6xAREZHBojC8AQW+j1dfCKe9zY7Trs5jDUSiMWKpDPFM\nbt3XNEi9fxfZwKddrxK020SiEZLDO0kNbSeeLRKJ6dIXERERheENIWx5VqNdm6c5O0GrMo21Fowh\nlsyQWOe630WD1Pt30WJvZOv5mGiM1MhOksPbiWcKRKLre4dcREREBp/C8IDy3Sbt+gLNuSnc2Ylw\n2hsQTabWddrbcoPW+3dpXZ0pcIHvE43FSW8bI1kaJZ4p9GWfREREZONQGB4Qh6e9zdKYHcev1wBL\nJJ4gus7T3pYbtN6/iw5PgbNEHIfMztNJFEeJZ/IagiEiIiKnTGG4TxanvbUWKjRnx2nNz2ItRKJR\noql030ofYLB6/y6nKXAiIiLSawrD68hzG3i1eRqzk7jlKazvYU0Y7JzC+rc8W1rXst6/D00VeL7P\nvX+X0xQ4ERERWUsKw2so8Nq0a/O4lRmaMxP4bgNrIOokiGdz697ybLlX6/17Yx96/y5aPgXOWnCy\nBU2BExERkTWjMNxDNgjCUccLc2HXh4U5jAETDae9xdKZ/q1tWe/fh6byPFXOEtDf3r+H16YpcCIi\nItIfCsNdWAxxrdo8zdlJWnNTYcszIJZKk+jDtLflXq3372+e3Z/ev4sOT4FrgoloCpyIiIj0hcLw\nCvktF6++QLM8TbM8cdS0t2LfW3kt9v59eLLA92cGo/fvohNNgYtni0TjTl/WJCIiIlubwvBJLE17\nmy+H097qVawNiMQdYsn+THs7Yn0WnplP88hUfqB6/y5anAJn221MNEJiaAfp4R2aAiciIiIDQWnk\nKEdMe5sZx61MgwUiEWKpDInicL+XeFTv3wLTbnwgev8u0hQ4ERER2SgUhlmc9jZPszyFOztJ4LXB\nGKLJFE5+qO+lD3Cy3r+VvvX+XaQpcCIiIrIRbckwHHge7fo8rflZGjPjeI06Bog4DtF0hni0/9sy\nyL1/Fx2eAhcQcRKHp8ClcwrAIiIisiH0P/WtA2ttWPdbnacx8wrt+TIAJhYlmkyT7OO0t+UGtffv\ncpoCJyIiIpvJlgjDtVdeZP6FZyEaCae99bnl2aJB7v27nKbAiYiIyGa1JcJw4HtEHAcnV+z3Uga6\n9+8iTYETERGRrWJLhOF+G+Tev4uOmQKXL2kKnIiIiGx6CsNrYNB7/y7SFDgRERHZ6hSGe2TQe/8u\n0hQ4ERERkcMUhrtwsO7w0FRY/vD4bHbgev8uOmYKXGk76ZGdmgInIiIiW56S0ApshN6/izQFTkRE\nROTkFIZPYrH370NTBf5lQHv/LtIUOBEREZGVURg+ykbp/btIU+BEREREVk9hGHD9Tu/fzsNvg9j7\ndzlNgRMRERHpjU0dhr/+LHz8X+DQwtlsT+zmfeePc92ucBTzRuj9u5ymwImIiIj03qYNw19/Fj78\nHWh4AIZxN8WdPzqd74wXeaXp8MyA9v5dpClwIiIiImtv04bhj//LYhA+rBVEeGCyyKUD1vt30dFT\n4OK5ItmzzsLJl4glUv1enoiIiMims2nD8KGF4/+9Ab74S8+t61pejabAiYiIiPRP12HYGBMFHgcO\nWmuv735JvbErBwePE4h3JFvrv5ijaAqciIiIyGDoxZ3h9wPPAPkenKtn/uD1y2uGQ8mIz3vPO9SX\n9WgKnIiIiMjg6SqFGWNOA/498FHgP/dkRT3ytgvC/4bdJCzbE80jukmsB02BExERERls3d6S/CTw\nB0DuRAcYY24HbgfYs2dPly+3Mm+7IPw1//LPqE+8hJMrrvlragqciIiIyMax6jBsjLkemLTW/sAY\nc82JjrPWfhb4LMAVV1wxGH3LekxT4EREREQ2pm7uDL8BuMEYcx2QBPLGmC9Za2/pzdIG2xFT4BKp\ncApccZhYOqcpcCIiIiIbxKrDsLX2DuAOgM6d4Q9u9iB8oilwsVRGAVhERERkA1Ibg1ehKXAiIiIi\nm1tPwrC19kHgwV6cq9+WpsC5jXAKXLagKXAiIiIim5TuDKMpcCIiIiJb1ZYNw5oCJyIiIiJbKgxr\nCpyIiIiILLdlEqDf6QKRGtlFcmhbGIA1BU5ERERkS9sSYTg9spNkcVhT4ERERETkCFsiDMdSGSDT\n72WIiIiIyIDRbVIRERER2bIUhkVERERky1IYFhEREZEtS2FYRERERLYshWERERER2bIUhkVERERk\ny1IYFhEREZEtS2FYRERERLYshWERERER2bIUhkVERERky1IYFhEREZEtS2FYRERERLYshWERERER\n2bIUhkVERERky1IYFhEREZEty1hr1+/FjJkCXly3FzxsBJjuw+tuVNqvldOerYz2a2W0Xyuj/VoZ\n7dfKaL9Wpp/7dbq1dvRkB61rGO4XY8zj1tor+r2OjUL7tXLas5XRfq2M9mtltF8ro/1aGe3XymyE\n/VKZhIiIiIhsWQrDIiIiIrJlbZUw/Nl+L2CD0X6tnPZsZbRfK6P9Whnt18pov1ZG+7UyA79fW6Jm\nWERERETkeLbKnWERERERkWNsqjBsjPl3xpifGGN+aoz58HE+njDGfKXz8ceMMWes/yoHxyns123G\nmCljzJOdX7/dj3UOCmPMXxljJo0xPzrBx40x5n909vOHxpjL13uNg+QU9usaY0xl2fX1J+u9xkFi\njBkzxnzXGPOMMeZpY8z7j3OMrrGOU9wvXWMdxpikMeb/M8Y81dmvO49zjN4jO05xv/QeeRRjTNQY\ns98Y883jfGxgr69YvxfQK8aYKPB/AdcCB4B/Ncbca6398bLDfgsoW2vPMca8E/hvwH9c/9X23ynu\nF8BXrLXvXfcFDqYvAJ8C7jnBx38dOLfz63XA/93571b1BV59vwAettZevz7LGXge8F+stU8YY3LA\nD4wx/+uof5O6xg47lf0CXWOLXOBN1tqqMSYOPGKM+Za19vvLjtF75GGnsl+g98ijvR94Bsgf52MD\ne31tpjvDvwj81Fr7vLW2Bfwt8Najjnkr8MXO7/8BeLMxxqzjGgfJqeyXLGOtfQiYfZVD3grcY0Pf\nB4rGmJ3rs7rBcwr7JctYa1+x1j7R+f0C4RvK7qMO0zXWcYr7JR2da6ba+WO88+voh4b0Htlxivsl\nyxhjTgP+PfD5ExwysNfXZgrDu4GXl/35AMd+YVw6xlrrARVgeF1WN3hOZb8Abur8OPYfjDFj67O0\nDetU91QOu7LzY8hvGWMu7vdiBkXnx4d7gceO+pCuseN4lf0CXWNLOj/CfhKYBP6XtfaE15feI09p\nv0Dvkct9EvgDIDjBxwf2+tpMYfh4310c/V3cqRyzVZzKXnwDOMNaeynwTxz+jk6OT9fXyjxBOCrz\ntcDdwNf7vJ6BYIzJAl8FPmCtnT/6w8f5X7b0NXaS/dI1toy11rfWXgacBvyiMeaSow7R9bXMKeyX\n3iM7jDHXA5PW2h+82mHH+buBuL42Uxg+ACz/ruw04NCJjjHGxIACW/fHuCfdL2vtjLXW7fzxc8Av\nrNPaNqpTuQalw1o7v/hjSGvtfUDcGDPS52X1Vac28avAl621XzvOIbrGljnZfukaOz5r7RzwIPDv\njvqQ3iOP40T7pffII7wBuMEY8wJh2eWbjDFfOuqYgb2+NlMY/lfgXGPMmcYYB3gncO9Rx9wL/G+d\n398MPGC3bqPlk+7XUbWINxDW5MmJ3Qvc2nni/5eAirX2lX4valAZY3Ys1osZY36R8OvRTH9X1T+d\nvfhL4Blr7SdOcJiusY5T2S9dY4cZY0aNMcXO71PArwDPHnWY3iM7TmW/9B55mLX2DmvtadbaMwjz\nxAPW2luOOmxgr69N003CWusZY94LfBuIAn9lrX3aGPNnwOPW2nsJv3D+jTHmp4Tfjbyzfyvur1Pc\nr98zxtxA+NT2LHBb3xY8AIwx/w9wDTBijDkA/CnhQxVYaz8D3AdcB/wUqAPv6s9KB8Mp7NfNwP9u\njPGABvDOQfnC2CdvAH4D+LdOnSLAHwJ7QNfYcZzKfukaO2wn8MVOJ6EI8HfW2m/qPfKETmW/9B55\nEhvl+tIEOhERERHZsjZTmYSIiIiIyIooDIuIiIjIlqUwLCIiIiJblsKwiIiIiGxZCsMiIiIismUp\nDIuIiIjIlqUwLCIiIiJblsKwiIiIiGxZ/z+hh9MhqI5lIAAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "# Q4\n\nGiven the function of random variables $$Y = X_1 + X_2 + X_3$$ for $X_i \\sim \\mathcal{N}(\\mu_i, \\sigma_i)$, \nwe will perform sensitivity analysis for each input variable using Sobol's indices. The Sensitivity index $$S_i = \\frac{V_i}{{\\rm Var}[Y]}$$ where each $V_i$ is the variance of the expected value of $Y$ conditioned on the i'th variable: $$V_i = {\\rm Var}[{\\rm E}[Y|X_i]]$$\n\nThis is really elegantly solved using the Sympy Package, which we used earlier for general gradient calculation. Unfortunately there is a bug in the current release in Uniform Distributions that does not return a uniform random variable after calculating an expected value, so we will first find those expectations and then plug them into the variance calculations. "
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "from sympy.stats import Uniform, E\n\nX1 = Uniform('X_1', 0.5, 1.5)\nX2 = Uniform('X_2', 1.5, 4.5)\nX3 = Uniform('X_3', 4.5, 13.5)\n\nY = X1 + X2 + X3\nV_y = variance(Y)\nprint(f'E(Y|X1)={E(Y, X1)}\\nE(Y|X2)={E(Y, X2)}\\nE(Y|X3)={E(Y, X3)}')",
      "execution_count": 21,
      "outputs": [
        {
          "output_type": "stream",
          "text": "E(Y|X1)=X_1 + 12.0\nE(Y|X2)=X_2 + 10.0\nE(Y|X3)=X_3 + 4.0\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# after the above\nfor n, exp_vals in enumerate([X1+12., X2+10., X3+4.]):\n    # hold one of the rv's at its mean:\n#     V_n = variance(Y.subs({rv:E(rv)}))\n    V_n = variance(exp_vals)\n    # calculate the index: \n    print(f'S_{n+1} = {V_n/V_y:.2f}')",
      "execution_count": 23,
      "outputs": [
        {
          "output_type": "stream",
          "text": "S_1 = 0.01\nS_2 = 0.10\nS_3 = 0.89\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "We would expect these values to sum to one (which they do) since there are no interaction terms in the function... being linear, we do not need to calculate variance contributions of non-linear terms. Since the indecies are normalized over all contributing terms, the linear term contributions will sum to 1. "
    }
  ],
  "metadata": {
    "hide_input": false,
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "language_info": {
      "name": "python",
      "version": "3.6.2",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "latex_envs": {
      "eqNumInitial": 1,
      "eqLabelWithNumbers": false,
      "current_citInitial": 1,
      "cite_by": "apalike",
      "bibliofile": "biblio.bib",
      "LaTeX_envs_menu_present": true,
      "labels_anchors": false,
      "latex_user_defs": false,
      "user_envs_cfg": false,
      "report_style_numbering": false,
      "autocomplete": true,
      "hotkeys": {
        "equation": "Ctrl-E",
        "itemize": "Ctrl-I"
      }
    },
    "varInspector": {
      "cols": {
        "lenName": 16,
        "lenType": 16,
        "lenVar": 40
      },
      "kernels_config": {
        "python": {
          "delete_cmd_postfix": "",
          "delete_cmd_prefix": "del ",
          "library": "var_list.py",
          "varRefreshCmd": "print(var_dic_list())"
        },
        "r": {
          "delete_cmd_postfix": ") ",
          "delete_cmd_prefix": "rm(",
          "library": "var_list.r",
          "varRefreshCmd": "cat(var_dic_list()) "
        }
      },
      "types_to_exclude": [
        "module",
        "function",
        "builtin_function_or_method",
        "instance",
        "_Feature"
      ],
      "window_display": false
    },
    "gist": {
      "id": "",
      "data": {
        "description": "Bayes, Kalman, and Sobol",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}