wiso/tangent.ipynb

## tangent.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tangent\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "import scipy\n",
    "from scipy import optimize\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 360,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>production (4mu)</th>\n",
       "      <th>ggH</th>\n",
       "      <th>VBF</th>\n",
       "      <th>VH</th>\n",
       "      <th>top</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>category</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>ggH_0J_CEN</th>\n",
       "      <td>0.084853</td>\n",
       "      <td>0.011622</td>\n",
       "      <td>0.016070</td>\n",
       "      <td>0.000215</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ggH_0J_FWD</th>\n",
       "      <td>0.149223</td>\n",
       "      <td>0.022524</td>\n",
       "      <td>0.033279</td>\n",
       "      <td>0.000509</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ggH_1J_LOW</th>\n",
       "      <td>0.070651</td>\n",
       "      <td>0.054836</td>\n",
       "      <td>0.045405</td>\n",
       "      <td>0.001756</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ggH_1J_MED</th>\n",
       "      <td>0.034186</td>\n",
       "      <td>0.060485</td>\n",
       "      <td>0.039882</td>\n",
       "      <td>0.001872</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ggH_1J_HIGH</th>\n",
       "      <td>0.006470</td>\n",
       "      <td>0.019141</td>\n",
       "      <td>0.012041</td>\n",
       "      <td>0.000801</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "production (4mu)       ggH       VBF        VH       top\n",
       "category                                                \n",
       "ggH_0J_CEN        0.084853  0.011622  0.016070  0.000215\n",
       "ggH_0J_FWD        0.149223  0.022524  0.033279  0.000509\n",
       "ggH_1J_LOW        0.070651  0.054836  0.045405  0.001756\n",
       "ggH_1J_MED        0.034186  0.060485  0.039882  0.001872\n",
       "ggH_1J_HIGH       0.006470  0.019141  0.012041  0.000801"
      ]
     },
     "execution_count": 360,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LUMI = 79808.76\n",
    "BR = 2.270E-03\n",
    "\n",
    "store = pd.HDFStore('data.h5', mode='r')\n",
    "EXPECTED_BKG_CAT = store['bkg']\n",
    "xs_x_BR_P4 = store['P4/xs_x_BR']\n",
    "eff_CP4 = store['/P4/eff']\n",
    "store.close()\n",
    "\n",
    "NCATEGORIES, NPROCESS = eff_CP4.shape\n",
    "\n",
    "eff_CP4.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 375,
   "metadata": {},
   "outputs": [],
   "source": [
    "# these are just random number as example\n",
    "\n",
    "# 5% systematics on all the generated number of events due to luminosity\n",
    "systematics_nsignal_gen_lumi = np.ones(NPROCESS) * 0.05\n",
    "\n",
    "# 6% systematics from btagging on ttH categories (last seven)\n",
    "systematics_efficiency_btagging = np.zeros_like(eff_CP4)\n",
    "systematics_efficiency_btagging[-7:, :] = 0.06\n",
    "\n",
    "# 2% systematics from identification\n",
    "systematics_efficiency_id = np.ones_like(eff_CP4) * 0.02\n",
    "\n",
    "systematics_nsignal_gen = [{'name': 'lumi', 'values': systematics_nsignal_gen_lumi}]\n",
    "systematic_efficiencies = [{'name': 'btag', 'values': systematics_efficiency_btagging},\n",
    "                           {'name': 'id', 'values': systematics_efficiency_id}]\n",
    "                 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_expected(ntrue, efficiencies, bkg):\n",
    "    expected = np.dot(efficiencies, ntrue) + bkg\n",
    "    return expected\n",
    "\n",
    "def chi2(mus, obs, ntrue, efficiencies, bkg):\n",
    "    expected = np.dot(efficiencies, ntrue * mus) + bkg\n",
    "    return np.sum((obs - expected) ** 2 / expected)\n",
    "\n",
    "# asimov dataset\n",
    "obs = compute_expected(xs_x_BR_P4.values * LUMI, eff_CP4, EXPECTED_BKG_CAT.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [],
   "source": [
    "# gradient of the chi2\n",
    "chi2_dmus = tangent.grad(chi2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 143,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute chi2 in the minimum\n",
    "chi2([1, 1, 1, 1], obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 0.])"
      ]
     },
     "execution_count": 144,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute gradient in the minimum\n",
    "chi2_dmus([1, 1, 1, 1], obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hessian(f):\n",
    "    vhp = tangent.grad(tangent.grad(f))\n",
    "    last_arg = vhp.__code__.co_varnames[vhp.__code__.co_argcount - 1]  # bad solution\n",
    "    def hf(x, *args):\n",
    "        H = []\n",
    "        for i in range(x.size):\n",
    "            v = np.eye(1, x.size, i)[0]\n",
    "            H.append(vhp(x, *args, **{last_arg:v}))\n",
    "        return np.array(H)\n",
    "    return hf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[117.28122811,  11.20455909,   5.29496405,   0.89595748],\n",
       "       [ 11.20455909,  14.86204398,   0.75786282,   0.20062095],\n",
       "       [  5.29496405,   0.75786282,   4.4902303 ,   0.80467926],\n",
       "       [  0.89595748,   0.20062095,   0.80467926,  10.86344109]])"
      ]
     },
     "execution_count": 146,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# hessian of the chi2 / 2\n",
    "hessian(chi2)(np.array([1, 1, 1, 1]), obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values) / 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Likelihood with Poisson"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "from scipy.special import gammaln, digamma\n",
    "\n",
    "def poisson(lam, n):\n",
    "    return np.exp(n * np.log(lam) - lam - gammaln(n + 1.0))\n",
    "\n",
    "def logpoisson(lam, n):\n",
    "    return n * np.log(lam) - lam - gammaln(n + 1.0)\n",
    "\n",
    "nmu = 4\n",
    "all_sys_ntrue = [np.array([0.1, 0.1, 0.1, 0.1])]\n",
    "\n",
    "def nll(pars, obs, ntrue, efficiencies, bkg):\n",
    "    mus = pars[:4]\n",
    "    ntrue_sys = ntrue\n",
    "    #for sys_ntrue in all_sys_ntrue:\n",
    "    #    ntrue_sys *= (1 + sys_ntrue * pars[4])\n",
    "    #ntrue_sys = ntrue\n",
    "    expected = np.dot(efficiencies, ntrue_sys * mus) + bkg\n",
    "    #p = logpoisson(expected, obs)\n",
    "    p = obs * np.log(expected) - expected - gammaln(obs + 1.0)\n",
    "    return -np.sum(p) + pars[4] ** 2 / 2.\n",
    "\n",
    "# define the derivative of gammaln\n",
    "from tangent.grads import adjoint\n",
    "@adjoint(gammaln)\n",
    "def dgammaln(result, x):\n",
    "  d[x] = d[result] * digamma(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "101.51786000815788"
      ]
     },
     "execution_count": 151,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# value of the nll in the minimum\n",
    "nll(np.array([1, 1, 1, 1, 0.]), obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {},
   "outputs": [],
   "source": [
    "nll_grad = tangent.grad(nll)\n",
    "nll_hessian = hessian(nll)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 0., 0.])"
      ]
     },
     "execution_count": 153,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nll_grad(np.array([1, 1, 1, 1, 0.]), obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[117.28122811,  11.20455909,   5.29496405,   0.89595748,\n",
       "          0.        ],\n",
       "       [ 11.20455909,  14.86204398,   0.75786282,   0.20062095,\n",
       "          0.        ],\n",
       "       [  5.29496405,   0.75786282,   4.4902303 ,   0.80467926,\n",
       "          0.        ],\n",
       "       [  0.89595748,   0.20062095,   0.80467926,  10.86344109,\n",
       "          0.        ],\n",
       "       [  0.        ,   0.        ,   0.        ,   0.        ,\n",
       "          1.        ]])"
      ]
     },
     "execution_count": 154,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = nll_hessian(np.array([1, 1, 1, 1, 0.]), obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values)\n",
    "#H = H[:4, :4]\n",
    "H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 9.63197130e-03, -6.74079472e-03, -1.02363153e-02,\n",
       "         8.83199446e-05, -0.00000000e+00],\n",
       "       [-6.74079472e-03,  7.25888777e-02, -4.21809477e-03,\n",
       "        -4.72149714e-04, -0.00000000e+00],\n",
       "       [-1.02363153e-02, -4.21809477e-03,  2.38489021e-01,\n",
       "        -1.67432791e-02, -0.00000000e+00],\n",
       "       [ 8.83199446e-05, -4.72149714e-04, -1.67432791e-02,\n",
       "         9.32935111e-02, -0.00000000e+00],\n",
       "       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
       "         0.00000000e+00,  1.00000000e+00]])"
      ]
     },
     "execution_count": 155,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cov = scipy.linalg.inv(H)\n",
    "cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.14011414, 0.28738281, 0.49848673, 0.32139308, 1.        ])"
      ]
     },
     "execution_count": 258,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# errors on mus\n",
    "np.sqrt(np.diag(cov))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "metadata": {},
   "outputs": [],
   "source": [
    "def correlation_from_covariance(covariance):\n",
    "    v = np.sqrt(np.diag(covariance))\n",
    "    outer_v = np.outer(v, v)\n",
    "    correlation = covariance / outer_v\n",
    "    correlation[covariance == 0] = 0\n",
    "    return correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 222,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  1.,  1.,  1., -1.],\n",
       "       [ 1.,  1.,  1.,  1., -1.],\n",
       "       [ 1.,  1.,  1.,  1., -1.],\n",
       "       [ 1.,  1.,  1.,  1., -1.],\n",
       "       [-1., -1., -1., -1.,  1.]])"
      ]
     },
     "execution_count": 222,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# correlations between mus\n",
    "correlation_from_covariance(cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 101.5178600081266\n",
       " hess_inv: array([[ 0.00943962, -0.00722565, -0.01093327,  0.00118301],\n",
       "       [-0.00722565,  0.07109579, -0.00669575,  0.00278579],\n",
       "       [-0.01093327, -0.00669575,  0.23395622, -0.0114784 ],\n",
       "       [ 0.00118301,  0.00278579, -0.0114784 ,  0.08613256]])\n",
       "      jac: array([-2.61610462e-06, -8.84158003e-07, -3.65705652e-07,  1.00346926e-06])\n",
       "  message: 'Optimization terminated successfully.'\n",
       "     nfev: 14\n",
       "      nit: 12\n",
       "     njev: 14\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([0.99999998, 0.99999995, 0.99999993, 1.0000001 ])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# minimize nll from a random point\n",
    "optimize.minimize(lambda mus: nll(mus, obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values),\n",
    "                  jac=lambda mus: nll_grad(mus, obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values),\n",
    "                  x0=np.array([0.,0.,0.,1.]),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def profile_likelihood_ratio_gridscan(nll, index_profile, xrange, nll_grad, nll_hessian):\n",
    "\n",
    "    start_point = np.array([1.1, 1.1, 1.1, 1.1])\n",
    "\n",
    "    global_minimum = optimize.minimize(lambda mus: nll(mus, obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values),\n",
    "                          jac=lambda mus: nll_grad(mus, obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values),\n",
    "                          hess=lambda mus: nll_hessian(mus, obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values),\n",
    "                          x0=start_point,\n",
    "                          method='Newton-CG')\n",
    "\n",
    "    result = []    \n",
    "    start_point_profiles = np.delete(start_point, index_profile)\n",
    "    for value in xrange:\n",
    "        res = optimize.minimize(lambda x: nll(np.insert(x, index_profile, value), obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values),\n",
    "                                jac=lambda x: np.delete(\n",
    "                                    nll_grad(np.insert(x, index_profile, value),\n",
    "                                             obs, xs_x_BR_P4.values * LUMI, eff_CP4.values, EXPECTED_BKG_CAT.values), index_profile),\n",
    "                                x0=start_point_profiles)\n",
    "        result.append(2 * (res.fun - global_minimum.fun))\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_intersections(xs, ys, val=1):\n",
    "    xs, ys = xrange, yscan\n",
    "    xs, ys = np.asarray(xs), np.asarray(ys)\n",
    "    argsort = np.argsort(xs)\n",
    "    xs = xs[argsort]\n",
    "    ys = ys[argsort]\n",
    "    argmin = np.argmin(ys)\n",
    "    xs_left = xs[:argmin]\n",
    "    ys_left = ys[:argmin]\n",
    "    xs_right = xs[argmin:]\n",
    "    ys_right = ys[argmin:]\n",
    "\n",
    "    def solve(xs, ys, val=1):\n",
    "        closest_to_one = np.abs(ys - val).argmin()\n",
    "        x_closest = xs[closest_to_one]\n",
    "        range_closest_to_one = max(0, closest_to_one - 3), min(closest_to_one + 3, len(ys))\n",
    "        ys_closest_to_one = ys[range_closest_to_one[0]:range_closest_to_one[1]]\n",
    "        xs_closest_to_one = xs[range_closest_to_one[0]:range_closest_to_one[1]]\n",
    "        roots = np.roots(np.polyfit(xs_closest_to_one, ys_closest_to_one - val, 2))\n",
    "        x_solution = roots[np.argmin(np.abs(roots - x_closest))]\n",
    "        return x_solution\n",
    "\n",
    "    return solve(xs_left, ys_left, val), solve(xs_right, ys_right, val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3419bd9c7f4b47d2bc3c2936f0c29e48",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=4), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ggH  : 1 ± 0.098\n",
      "VBF  : 1 ± 0.269\n",
      "VH   : 1 ± 0.488\n",
      "top  : 1 ± 0.305\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAx8AAALaCAYAAABQ9u4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzs3Xl4nFd59/HvmRmNNmuzNlveN9mOnTjOQjYCYQuBsLe0QKA0lMJVlhYKLeUtpQttaSktLaGUtlBCIaQQCISGkkA2yEI2O3a8yrsty5YlS7J2abbz/vHMyIrxomVmnu33uS5dkaXR6M4j6Zy5n3PucxtrLSIiIiIiIoUWcTsAEREREREJByUfIiIiIiJSFEo+RERERESkKJR8iIiIiIhIUSj5EBERERGRolDyISIiIiIiRaHkQ8RDjDGPGGMeOMfnvmqMOVTkkEREJA+MMT8yxhwzxpz1tZcx5lpjjDXG3JqdC+yktxFjzHZjzEeMMeaMr7PnePtUcf7PRKYn5nYAIiIiIiHwTeD1wMuBs91kugUYBb4PvBt4Bvj97OfmADcBXwBKgH8442v/Hbj9jI+15yNokXxT8iEiIiJSeP8L9APv5IzkwxhTAvwG8CNr7UB2cWPAWvvkpIc9YIy5FPh1fjX5OHrGY0U8S9uuRPLAGPM+Y8xBY8yoMeZxY8wlxpiUMeYvsp83xphPZZfch40xP5m0xP7b7kYvIiKFZq0dA74HvMUYU37Gp18NNOCsjpzPAM7Kh4hvKfkQmSVjzBtwlrwfA94MfAe4C5i8L/fDwGeyn3sz8ATwrXM/pYmd+XbG84mIiP98E6gC3nDGx98JdAP3T/rY5LmgxhjzNuBmnHnkTJGzzBkinqRfTpHZ+xTwlLX2Xdl/32eMSQD/BmCMiQKfBL5jrf1o9jE/NcZUAp84y/O9HEie43sdzl/YIiJSZL/AGcdvIZtEGGNyychXrbWpSY8921zwDeBzZ3nev8y+TaYbVuJJWvkQmYVsYrER+OEZn7pn0vsLgXkXeMxkTwNXnuXtXI8XEREfsNZa4NvATcaY+uyH3wyU86tbribPBS8BPpZ97FfO8tT/xq/OGSKepJUPkdlpxPk76j7j45P/Pf8sHwPoOsdzDlprnz3zg8aYkzOKUEREvOSbOKvhv4GTNLwTaLPWPnPG486cCx7NHrP7eWPMF621OyZ97tjZ5g0RL9LKh8jsnATSOEnIZJP/3XmWjwE0FSooERHxJmvtLmATcIsxZh7O9qpz1QCeaWf2v+sLEZtIMSj5EJmF7P7czcCbzvjU5H+34yQg53uMiIiEx7eAa4H/h/NabKrJx8XZ/565ki7iG9p2JTJ7fw3cY4z5Js4Eshr4EJABMtbatDHms8C/GGOOA/8HXI2z5E72cSIiEh53Ap/HmSses9YeOstjqo0xV2ffLweuwDngZBtO4bqIL2nlQ2SWrLU/At4PXI9TFP524D04f1/92YfdBnw6+7l7so/9UPZz/YiISGhYa08AP8U5kepcqx5XAr/Mvt0HvA/4T+BlZ5yKJeIrxjl4QUTyyRjzeuBHwA3W2p+f4zEfBr4ILLXW6ghdERERCTxtuxKZJWNMDc7WqwdxVjEuAf4Up6DwF9nHtOKshjwOjADXAH8C/ECJh4iIiISFkg+R2UsCS3C6nM8FenHqOv7Inl5aHAUuB94L1AAngK/hFBuKiIiIhIK2XYmIiIiISFGo4FxERERERIpCyYeIiIiIiBRFaGo+Ghoa7NKlS90OQwqgZyjBsf5RSqIR1syrcjucvLDWMjo6CkB5eTnGGMWSZxlr2XFsAIDlDZVUloZmODynTZs2nbTWNrodh9tmO18kUhnaTgwCsKJxDhXxaJ4iE68L6njpB7r2xTWb+SI0s+3SpUt59tln3Q5DCuCBnSd4738/SzRiePIzNxGL+n9Bb3x8nO3btwOwfv16SktLFUue7TkxyI1fcPp0/fSPX8aiuRUuR+Q+Y4xOXmP284W1luv+7iGO9Y/xoZtW84EbVuYxOvGyoI6XfqBrX1yzmS/8/ypNQm/h3HIA0hnLicFxl6MRv+joc+6QRQzMqylzORoJEmMMV6+oB+CX+3tcjkZExFuUfIjvLagtn3j/aO+Ii5GInxw95SQfzdVllARgtUy85doVDQA8e6iPRCrjcjQiIt6hGVd8r6qshJryEgA6si8oRS4kt/IxOXkVyZdrsisfo8k0W4+ecjkaERHvUPIhgZB7AXm0T8mHTE0uUV1Qp+RD8m9BbTmLs3VE2nolInKakg8JhIXZF5AdAUk+IpEI9fX11NfXE4m4+2fqpVjyqaPP2aKnlQ8plGuzqx9P7D/pciRSLEEdL/1A194/QnPalQRb7u710VPBqPkoKSnBK0dDeymWfNLKhxTaNSvq+Z9n2tl85BRjyTRlJTpyN+iCOl76ga69fyg1lEBYWOdsbwjKyocUViKVoSt7MppWPqRQrlnurHwkUhk2H+lzORoREW9Q8iGBkHsBeezUGJmMdTka8brj/aPY7K/JQq18SIE0VZexorESUN2HiEiOtl1JIOReQCbSGbqHxmmu9nffhkQiwaFDhwCn4Vk8HlcseTR5haxFKx9SQNesqGd/97CSj5AI4njpF7r2/qGVDwmEyXevj/b5v+7DWsvg4CCDg4NY6+5KjpdiyZdcj4+5lXEq4roHI4VzzXKn38eW9lMMj6dcjkYKLYjjpV/o2vuHkg8JhJryEuaUOi8iddyuXIh6fEixXL18LgCpjOWZQ70uRyMi4j4lHxIIxhj1+pApmzjpSsmHFFj9nFLWzKsC4PF9OnJXRMTV5MMY8zZjzKPGmAFjzK+sRxtjmowx3zDG9GQfs8UY0+JGrOJ9E70+1OVcLmBi5UPF5r7g97ni+lXO1qvH9qnuQ0TE7c3OfcCXgXLgPyZ/whhTBjwIPAmsBnqBtcDQTL6RtZbx8fFZBSveNq/aKS470jPs+5/15Pjd/n/xUiz5kqsLap5TEpj/p4Ar2lwB+Z8vrlpay38+CruOD9DRM0DDnNK8Pbd4SxDHS7/QtfcPV5MPa+39AMaYG87y6XcDtcAHrLXJ7Md2TOf5jTH1QD1Aa2sr27dvn3mw4nnRsWEADp44Faif9d69e90OYYKXYpmpjLUcz66OpQa62L59wOWI5EIKPVdkn7tg80VFyhKLQCoDd/3iea5frBW3MAjCeOlXuvbe5uWaj5cBe4Hbs0vpu40xH53mc3wYaAPaentV6Bd0jRVO9+DukbROupBz6hvLkMr+euR+Z8TX8jFXQAHni9KYYW2DszK79UQir88tIuI3bm+7Op8GnEnlI8CtwCXAfcaYLmvtHVN8jtuAbwO0tLS0rV+/viCBijdk6vrhyWdIpKFl2Wrq5/j3jO9kMsmxY8cAaGlpoaSkRLHkyeYjp4BuAF5y+XpqK/z9/yN5mSugwPPFjX2H2PbAPnb2Zli3bh3GmLw+v3hD0MZLP9G19w8vJx+DQIe19l+y/37WGPMt4I3AlCYUa20P0ANwxRVXUFqqfbZBtqypZuL97pE0LfX+/XmXlpbS2trqdhiAt2LJh65hp165Mh6lqbZSLwL9b9ZzBRR+vnjpmmb+8YF9nBgY5+hAkpVNVXl9fvGGoI2XfqJr7x9e3na1BTjb3hntp5Gzqq+MU1bi/ErruF05l4ljduvKlXgEgy/minUtNROrbI/u1ZG7IhJebh+1G82eVBLP/rss+2aA24F6Y8wHs4/bANwC3O1exOJlk3t9dJzyf5dzKQw1GPSfIMwV0YjhuhXZI3eVfIhIiLm98vEuYBS4H4hm3x8FllhrDwOvBd4LDADfA/7CWvsdl2IVH1hYVwH4f+UjmUyyf/9+9u/fTzKZvPAXhCSWfJi88iG+EYi54sXZfh9PHughmc64HI0UQtDGSz/RtfcPt4/avR3nrtW5Pv8IsLFI4UgA5F5Qdvg8+chkMpw6dQqAhQsXKpY8Or3yUeFyJDJVQZkrXrzSST6GE2meO3KKFy2b63JEkm9BGy/9RNfeP9xe+RDJq1yXc7+vfEhhWGu18iGuWTS3gqX1TtL72D5tvRKRcFLyIYFyuuZjVL0+5Ff0jSQZSaSB04mqSDHltl49trfb5UhERNyh5EMCJVfzMTSeon9Uez7lhY70nj6IYPFcbbuS4nvxykYAth7tZ2BMY5SIhI+SDwmUyXeztfVKzpRLPiriUeor/duEUvzrmhX1RAykM5Zf7u9xOxwRkaJT8iGB0jinlHhUvT7k7NqzycfiuRXq8SGuqCkv4ZKFtYCO3BWRcFLyIYESiRhaasuA00eqiuQc6XGSj0XaciUuuj5X96GicxEJISUfEjine334t9GgMYaqqiqqqqpcv0PvpVhm68iklQ8Rt+SO3D14ctjX45T8qiCNl36ja+8frvb5ECmEiROvfLztKh6P09ra6nYYgLdimS0lH+IFGxfXURGPMpJI89jek7ztRYvdDknyJEjjpd/o2vuHVj4kcNTrQ84mmc5wvN/5nVDyIW6KxyJcvbwegEe19UpEQkbJhwTORJdz1XzIJMdOjZLJtn5RzYe4Lbf16ol9J8lk1JNIRMJD264kcHI1H/2jSQbHklSVlbgc0fQlk0k6OjoAWLBgASUl7v0/eCmW2Zjc40MNBsVtuaLzvpEkO48PsH5BjcsRST4EZbz0I117/9DKhwTOgkkvLP26+pHJZOjp6aGnp4dMJqNY8iCXfDRXl1JWEnU5Ggm7lU1zmFftnMz3C3U7D4ygjJd+pGvvH0o+JHCaq0qJRZyTLvxcdC75pWJz8RJjzMTqxyNtSj5EJDyUfEjgxKIR5tU4dxRVdC45uQaDqvcQr7hhdRMAmw/3MTCWdDkaEZHiUPIhgbRQRedyBq18iNe8eGUDEQOpjOUJnXolIiGh5EMCaUGt8wKzvVcNvMSR626u5EO8oqaihMsW1wHaeiUi4aHkQwIp9wLziJIPAfpHkgyMpQAlH+ItN6xuBJzkw1oduSsiwafkQwJpSX02+egZ0YQuL0hClXyIl+TqPjoHxmg7MehyNCIihafkQwIpV1Q8OJ7i1Ij/CjmNMVRUVFBRUYExRrHMUi75KI1FaKwqdTkakdMuml9Nw5w4AD/X1ivfC8J46Ve69v6hJoMSSLmVD4DDvSPUVcZdjGb64vE4a9eudTsMwFuxzFR73+l6D01K4iWRiOElrY3cvbmDR9q6ef9LV7gdksxCEMZLv9K19w+tfEgg1VfGqYg7jeRU9yE66Uq8LLf16tnDvQyNp1yORkSksJR8SCAZY04XnfcMuxyNuE09PsTLrs8euZtMWx7XkbsiEnBKPiSwJorOfbjykUqlaG9vp729nVTK3TuhXoplpo4o+RAPq6uMs2FRLQA/36O6Dz8LwnjpV7r2/qHkQwIrt/JxuMd/yUc6naarq4uuri7S6bRimYVUOkNHttO9tl2JV93Q6my9+rmO3PU1v4+XfqZr7x9KPiSwFtdXAmo0GHbH+8dIZZwXc0o+xKty/T46To2yr2vI5WhERArH1eTDGPM2Y8yjxpgBY8w518iMMb9njLHGmE8VMz7xt9wLzeMDY4yndBckrCYnn4vmlrsYicxUGOaKixfUMDd7Kp+6nYtIkLm98tEHfBn4yLkeYIxZAnwM2FasoCQYlmSTD2uhvXfU5WjELbl6j4Y5pVTEdbq4TwV+rohEDC9Z1QCo7kNEgs3V5MNae7+19k7gwHke9jXgT4He6T6/MabeGNNqjGlV8VH4tNSWE8m2dNDWq/A6fcyuVj38qtBzBXhjvsgdufv0wV6GdeSuiASU2ysf52WMeT8wbK39zgyf4sNAG9DW1dWVv8DEF+KxCC21zgvOwzpuN7TU4yP48jBXgAfmi5e0NmIMJNIZfrm/x5UYREQKzbPJhzFmMfAp4AOzeJrbgNXA6qamprzEJf4y0etD265Cq10nXQVanuYK8MB8MbcyziULnSN3H9mjG2YiEkyeTT6ArwJ/ba3tmOkTWGt7rLV7rLV7YjHt9Q6j070+/LXyYYyhtLSU0tJSjDGKZRbUYDDwZj1XgHfmi5e2OqdePaIjd33J7+Oln+na+4eXX5G/CrjcGPM32X/XAFcaY15trb3exbjERxbN9WejwXg8zvr1690OA/BWLNM1OJakdzgBaOUjwAI1V9ywupEvPriXo32jHDg5zIrGOW6HJNPg5/HS73Tt/cPV5MMYEwVKgHj232XZT40Di854+F3Ao8A/Fi1A8b0lc51eH0d6R7DW6m5IyEw+5UwrH/4Vprliw8Ja6ipK6BtJ8khbt5IPEQkct7ddvQsYBe4Hotn3R4El1tqjk99wJpkBa+0J98IVv8ltuxpLZugaHHc5Gim23IpXPBqhubrsAo8WDwvNXBGNGF6S3Xr10G5f/i+IiJyXqysf1trbgdun+NgbChmLBNPku91Hekd88wI0lUpx4oTzwqO5uRk396B7KZbpytV7LKwrJxrRqpdfhW2uePmaJu7ZcoynDvQyMJakuqzE7ZBkivw8Xvqdrr1/uL3yIVJQNeUl1FY4E/fhHv/UfaTTaTo7O+ns7CSddrc7u5dima4jKjYXH7qhtYloxJDKWB7dc9LtcGQa/Dxe+p2uvX8o+ZDAW+zTonOZPfX4ED+qqSjhyqV1ADy4S1uvRCRYlHxI4E0kH2o0GDrtSj7Ep165thmAh9u6SGd05K6IBIeSDwk8rXyEUzpjOZptMKhtV+I3r8gmH30jSTYf6XM5GhGR/FHyIYF3utGgko8wOTEwRiKdAbTyIf6zrKGS5Q3OUeEPaOuViASIkg8JvNxd75NDCYbGUy5HI8XSPinZXDS33MVIRGbmFWubAHhwV5fLkYiI5I+SDwm8JfWVE++3a/UjNHIrXXMr41TpqFLxodzWq31dQxxWzZqIBISSDwm8edVllESdHg9+Om43Fot55pxyL8UyVe06Zld87oolddSUO4nzA1r98A0/jpdBoWvvD/oJSeBFI4ZFdRUcODnsm5WP0tJSNmzY4HYYgLdimY6JHh912nIl/hSLRrhhdSP3bDnGg7tO8DsvXuZ2SHIBfh0vg0DX3j+08iGhsDhbdH64V1sXwkI9PiQIcluvnj7odDsXEfE7JR8SCqeP2x11ORIpltzPWsmH+NlLWxuJZbud/2JPt9vhiIjMmpIPCQW/NRpMpVJ0dnbS2dlJKuXuCV1eimWqRhIpTg6NA0o+xN9qyku4culcQKde+YEfx8ug0LX3DyUfEgq5F6BH+0ZJZXs/eFk6naajo4OOjg7S6bRimaYjLzhmV8mH+FvuyN2H27p8MX6FmR/Hy6DQtfcPJR8SCrnjdlMZy/H+MZejkUI7dNJZ4YpHI7TUquBc/O2V2bqPUyNJNh855XI0IiKzo+RDQmFykzl1Og++A9nkY3F9BdGIcTkakdlZ2lDJikbnBsqD6nYuIj6n5ENCoSIeo7GqFFDyEQa5lY9lDZUXeKSIP+RWPx5Q8iEiPqfkQ0IjV/fhp0aDMjOHTjo/YyUfEhS5I3f3dw9PJNciIn6k5ENCY0k2+fBLo0GZudy2q6X1Sj4kGC5bXEttRa7buVY/RMS/lHxIaKjRYDgMjiUnjtnVyocERSwa4WWrnVOvlHyIiJ8p+ZDQmLztylrrcjQXZozBGG8US3splguZvK1OyYcEyY0Xne523juccDkaORc/jZdBo2vvDzG3AxApliXZlY/BsRT9o0lqK+IuR3RupaWlXHbZZW6HAXgrlqnIbbkqL4nSXF3qcjQi+fPS1Y2UxiKMpzI8sPMEv3HlIrdDkjP4bbwMEl17/9DKh4TG5GZzKjoPrlwx7tKGSt0Bk0CpiMd4SWsjAPfv6HQ5GhGRmVHyIaHROKeU8pIooON2g+zgxDG76mwuwfPqdfMAeHTfSYbGUy5HIyIyfUo+JDSMMRN1H15PPtLpNN3d3XR3d5NOpxXLNBxUjw8JsFeubSIaMSRSGR5p63I7HDmD38bLING19w9Xkw9jzNuMMY8aYwaMMakzPvdaY8xDxpiTxpi+7OOudytWCYbciVdHPL7tKpVKceTIEY4cOUIq5e7dTS/FMhWHenTMbtBorjittiLO1cvnAnD/Dp165TV+Gy+DRNfeP9xe+egDvgx85CyfqwNuA1YCjcC3gZ8YY1RhJzM2ceKVjtsNpL7hBKdGkgAsb1TyESCaKybJbb16eHcX4ynd4RURf3E1+bDW3m+tvRM4cJbP3WGt/YG19pS1NmWt/TdgCLiy6IFKYCzNbsXJdcCWYDnYczqp1MpHcGiueKEbL3KSj6HxFE/s63E5GhGR6XF75WPKjDEXAw3Atml8Tb0xptUY06olOAFYnk0+OgfGGFaxZuDkTrqqKosxt9K7RylL4cxkrsh+nW/mi3k1ZVy6qBaA+7br1CsR8RdfJB/GmCbg+8DnrbV7p/GlHwbagLauLhXmyQu34uQKkyU4cj/T5TpmN5RmMVeAz+aL3NarB3adIJ3xftNUEZEczycfxpgW4GHgp8Anp/nltwGrgdVNTU35Dk18qLmqbOK4XSUfwXNwUo8PCZdZzhXgs/ni1eucbuc9wwmePdTrcjQiIlPn6eTDGLMUeBT4ibX2Q9baad3esdb2WGv3WGv3xGJq5i4QiZiJF6YHupV8BI1Ougqn2c4V4L/5YnnjHFqb5wBwnxoOioiPuH3UbtQYUwbEs/8uy74ZY8wa4DHgTmvtx92MU4IlV/dx8OSQy5FIPllrOZhNKHXSVbBorji73Narn+44wQzyLRERV7i98vEuYBS4H4hm3x8FlgCfABYAHzHGDE16u8W1aCUQci9MvbztKh6Ps3HjRjZu3Eg87m7htJdiOZ/uoXGGE86xo1r5CBzNFWeRSz46To2yvWPA5WgE/DNeBpGuvX+4fdTu7dZac5a3Q9baW7Pvzznj7Q43Yxb/y3W+PnBy2LN3C40xRCIRIpGI64XTXorlfA5O2kanmo9g0VxxdutaqllQWw7A/dp65Ql+GS+DSNfeP9xe+RApulzyMTiW4uRQwuVoJF9y9R71lXFqyktcjkak8IwxE6sfSj5ExC+UfEjoLG+YM/G+V7depdNpent76e3tJZ12t4Oxl2I5n4PZxpFa9ZAwuWm9k3zs7Rpif7fq2Nzml/EyiHTt/UPJh4ROTUUJ9dkGdF4tOk+lUhw8eJCDBw/idsMzL8VyPrmf5TIlHxIily+pmxjPtPrhPr+Ml0Gka+8fSj4klJbpuN3AOZRd+VDyIWESjRhuzPb8ULdzEfEDJR8SSrkTrw54dNuVTE8mYydqPpR8SNjctH4+AM8f7ae9d8TlaEREzk/Jh4TSsmzdh1drPmR6jg+MMZ7KADpmV8Ln2hX11FU4hyzc+/xxl6MRETk/JR8SSrm744d7hklnvHncrkzdoZOTj9mtcDESkeIriUYmCs9/vO2Yy9GIiJyfkg8Jpdy2q2TacrRP2xT8LreCNa+6jIp4zOVoRIrv5otbANjeMfCCZFxExGuUfEgoLamvINeDSHUf/pdLPrTqIWF19fK5E6de/Xibtl6JiHcp+ZBQKo1FWVjndAY+qBOvfC93p1fF5hJWsUlbr1T3ISJepv0JElrLGubQ3jvKAQ/2+ojH42zYsAGAaDSqWC7goJIPEW6+ZD53PHWEXccH2N89xIrGORf+IskrP4yXQaVr7x9a+ZDQWp59oerFE6+MMcRiMWKxGCa3P0yxnFUqneFI9nhRnXQlYXbVsnoa5jhbr/5Pqx+u8Pp4GWS69v6h5ENCK1d0rm1X/tZxapRU9sSy3M9UJIyiEcNrsj0/VPchIl6l5ENCK7dF51j/GKOJtMvRvFAmk6G/v5/+/n4ymYxiOY/cgQERA4vmquBcwu11lzjJx+7OQfZ1DbocTfh4fbwMMl17/1DyIaG1fNJ+aK9tvUomk+zbt499+/aRTCYVy3nkis1basspjWmfr4TbFUvn0lRVCqjw3A1eHy+DTNfeP5R8SGjNry6jNOb8CXgt+ZCp00lXIqdFI4bXXpzdeqXkQ0Q8SMmHhFYkYiZesB704IlXMjUHlHyIvEBu69XeriH2nNDWKxHxFiUfEmq5F6wHVHTuW4d6sg0GddKVCACXLa5jXnUZAPduPeZyNCIiL6TkQ0ItdzqSupz703gqTUffKADLdNKVCOCs6ua2Xt277TjWWpcjEhE5TcmHhNqyBqfo/ED3kCZoH2rvHSF7yi7LtPIhMuHm7NarA93D7O7U1isR8Q4lHxJquW1XA2Mp+kZ0Oobf7M9ulyuJGhbWlbscjYh3XLa4lgW1zt+ECs9FxEuUfEiorZi0VedAt4rO/WZfl/MzW1pfSSyq4UwkxxjDay+eB8C9zx/Tyq6IeEbM7QBE3FRbEaeuooS+kSQHTg5zxdK5bocEQDweZ/369RPvK5azyyUfq5rnXOCRIuHzukta+M9HD3KoZ4StR/u5dFGt2yEFnpfHy6DTtfcP3SqU0Dt93K53is6NMZSWllJaWooxRrGcw95sB+eVTVUuRyLiPZcsrGF5dnz74XMdLkcTDl4eL4NO194/lHxI6OU6nWvblb9kMpb9XU7CuLJJKx8iZzLG8KaNCwD4363HSKYzLkckIuJy8mGMeZsx5lFjzIAxJnWWz99kjNlhjBk1xmw3xtzoRpwSbF5c+chkMgwNDTE0NEQm4+4LBi/FMlnHqVFGk2kAVin5CDTNFTP3xktbAOgZTvDYvpMuRxN8Xh0vw0DX3j/cXvnoA74MfOTMTxhjlgN3A58FarL//YExZmkR45MQyG1LONQzQjrjjaLMZDJJW1sbbW1tJJPunsLlpVgmy9V7RIy6m4eA5ooZWlJfyWWLnVqPe7T1quC8Ol6Gga69f7iafFhr77fW3gkcOMun3w1sstZ+y1qbsNbeAWzOflwkb3LN6RKpDMdOjbocjUxVLvlYPLeCspKoy9FIIWmumJ3c1qv7d5xgePxXFo5ERIrK7ZWP89kAbDrjY5uzH58SY0y9MabVGNOaSmnAlbNbWl9JrjZNnc79Q8XmkjXruQKCPV/cfPF8YhHDaDLNz3aecDscEQk5LycfVUD/GR87BVRP4zk+DLQBbV1dXfmKSwKmrCRKS43TjOugis59Q8fsSlY+5gq5CWJcAAAgAElEQVQI8HxRP6eUl7Q2AvADbb0SEZd5OfkYxNm/O1ktMDCN57gNWA2sbmpqyldcEkDLG71XdC7nZq1lbzb5WNmo5CPk8jFXQMDni9zWq0f3dtM9OO5yNCISZnlJPowx1caYD+bjuSbZClx2xsc2Zj8+JdbaHmvtHmvtnlhM/RTl3HJF59p25Q9dg+MMjjlbY7Ty4S8FmC9mPVdA8OeLV61tpjIeJWOdjuciIm6ZVfJhjLnOGHM7cBz4+xl8fdQYUwbEs/8uy74Z4L+BK4wxbzfGlBhj3g5cDnxjNjGLnE3utKQD3Uo+/CC35QpghVY+fGE284Xmitkrj0d59fp5gBoOioi7pp18GGPmGmM+aozZATwKtOIcfzh/Bt//XcAocD8Qzb4/Ciyx1u4H3gJ8Cmf5/FPAm621h2bwfUTOK1e03HFqVKfB+MDeE06x+YLacipLg3eXOijyOF9orsiDN2e3Xm092q+mqiLiminP2saYlwG/C7wZ58z1bwK/Zq3dPdNvbq29Hbj9PJ+/D7hvps8vMlWtk7bu7O0a4tJFtS5GAyUlJaxdu3bifcXyQhP1Hmou6En5ni80V+THtSsaaKwqpXtwnB9uOcYfvqrV7ZACx4vjZVjo2vvHlFY+jDF7gZ/gLHm/FVhkrf3EbBIPES9prCqlptwZrPZk76q7KRKJUFFRQUVFBZGIu+dCeCmWnImTrpR8eI7mC++KRgxv2OB0PL9nSwfWeqOpapB4cbwMC117/5jqT2cFcBB4BHjcWpsuWEQiLjDGTLyQ3euB5EPOb59WPrxM84WH5bZeHe4Z4bn2Uy5HIyJhNNXkYxnwfeCTwDFjzHeNMTdli/1EAmFVs1P3seeE+3uhM5kMo6OjjI6OkslkFMskvcMJeoYTgE668ijNFx62rqWaFdmjxVV4nn9eGy/DRNfeP6aUfFhrD1trPwUsBt4OVAD3Au3GmL81xmjjqPheru5j8klKbkkmk+zcuZOdO3eSTCYVyySTfz4rG9Xd3Gs0X3ibMWZi9ePe54+TSOlFWj55bbwME117/5jWpjhrbdpa+0Nr7euApcDXgHcCu4wxjxcgPpGiaW0+feLVkE688qy9Xc62uMaqUmoqVFToVZovvOuNlzrJR+9wgofbgtXNXUS8b8YVOdbao9baP8eZVN4IdOcrKBE3TN7Co7oP79p7QsXmfqP5wlsWza3gmuX1ANz17FGXoxGRsJn1cQDW2oy19l5r7ZvyEZCIWxrnlFKbvZO+1wN1H3J2+7uVfPiV5gvveOsVCwF4uK2LrsExl6MRkTCZVncuY8zDwNnO5rPAGLAXuN1auyUPsYkUlTGG1qYqnj7U64njduXscomhTrryNs0X3vaa9fP59D07GBpP8cPnOnjfS1a4HZKIhMR0Vz72AJfhFBJ2Zt8WZT/WA7wceNoY89J8BilSLLmtV3s8UHQuv2pgLEnngHOXNteVXjxL84WHlcejvH6D02j+rmePqueHiBTNtFY+gFPAncAHrbUZgOzxiV8CBqy1v2WM+QrwN8CL8xrpLFlrGR8fdzsM8bjl9eUA7OkccPX3ZfL3dvv31kux7O7on3h/cW2J6/HIeWm+8Lg3XtLMnU+3s7driGcOdLNhYY3bIfmel8bLsNG194/pJh+3AtfnJhIAa601xnwReAznXPcvAz/PX4j5MTo6yvbt290OQzwuOuwMWJ0D4zz13PNUlrjfJXXv3r1uhzDB7Vh+fnAEgDlxw/GDe+hU6wgv03zhcTFrWVAVpWMwzVcf3M77L1fykU9uj5dhpmvvbdN9ZVWKc1rJmZYC8ez7I4BeEYgvLa4+nY8fHdBxu15zdNBplr2oOoZ61nme5guPM8bw8qXOau9jR8YYT2nrlYgU3nRXPr4PfNUY8yfAU9mPXQV8Frgr++9rgN35CS9/ysvLWb9+vdthiA/UPfhz+kaSZOY0s379ArfDkUn6t24Bhrl4SRPr1691Oxw5P80XPtC0ZJxv73iMkZSlI9LAG9bPdzskEQm46SYfHwQ+j9MsKtfdKwl8Ffh49t/PAb+Tl+jyyBhDaWmp22GID6xqruLpg70c6BnT74zH7D85DMDq+TX62Xif5gsfWFRayg2tjTy4u4sfbu3krVcudTskEQm46XY4H7XWfhCoBzZm3+qttR+y1o5lH7PdWrsj/6GKFEdr9sSrXCdtN1hrSSQSJBIJ10+h8Uoso4k0R/tGAfX48APNF/6R6/nx+L4e2ntHXI7G37wyXoaRrr1/zKia1lo7bK19Pvs2nO+gRNzU2uwc4epmo8FEIsG2bdvYtm0biUTCtTi8FMv+7iFy84l6fPiH5gvve/maZuZWOmU439+sjuez4ZXxMox07f1j2smHMWadMeYbxphnsm+3G2MuKkRwIm5Yle0f0TkwRv9o0uVoJGdftvdKZTzK/Joyl6ORqdB84Q/xWIQ3XerUt31v01EyGd01FpHCmVbyYYy5GdgCrAQeyL6tArYYY16T//BEii+37Qpgn4tbr+SFcsnHyuYqnXTlA5ov/CW39epo3yhPHuhxORoRCbLpFpz/NfBP1tpPTP6gMeZzOI2ifpKvwETcUj+nlPrKOD3DCfacGOLyJXPdDkk4XYOzslFbrnxC84WPrJ1fzfoF1WzvGOCuTUe5dmWD2yGJSEBNd9vVWpyTSs70VUBL6RIYq7KrH3tOaOXDK/ZmVz5WNSv58AnNFz7zG1csAuAn248zMKYtpyJSGNNNPnpxJpQzrcl+TiQQvFB0LqclUhkO9zin8OikK9/QfOEzb9jQQjwWYSyZ4YfPdbgdjogE1HSTj28D/2GMudUYsyr79h7gP4Bv5T88EXesyiYfWvnwhkM9w6SzRbA66co3NF/4TG1FnNdd7DQZvOPJIzquVEQKYro1H5/ESVj+DYhnP5YAvgT8aR7jEnFVa/YFbtfgOP0jSWoqSi7wFVJIuRWo0liEhXUVLkcjU6T5woduuXoxdz/XQduJQTYf6VPNm4jk3bSSD2ttEvhDY8yngBXZD++z1o7mPTIRF+W2XQHs6RrkyqXFnYBjsRgrVqyYeN9NXohld+cA4Kx6RCM66coPNF/402WL61jdXEXbiUHuePKIko9p8sJ4GVa69v5xwZ+OMeb/pvAYAKy1r81DTCKuq6uM0zCnlJND4+w5UfzkIxqNUltbW9TveS5eiGXXcSf5WDu/2tU45Pw0X/ifMYZbrl7Mp+/Zwb3bjvNnr7uIusr4hb9QAG+Ml2Gla+8fU6n5ODGNt7wyxswzxnzHGNNtjOkzxjxkjNmQ7+8jcja5wmYVnbtv13Gn9kbJh+dpvgiAN21cQHlJlEQqo47nIpJ3F1z5sNbeWoxAzuHLQBXQCgzjnBt/rzFmsVUlnBRYa/Mcfnmgx5Wic2stmUwGgEgk4mpTPbdj6R9J0nHK2alzkZIPT9N8EQzVZSW88dIW/ueZdr791BF+58XL1NhzitweL8NM194/pnvaVbGtBO6y1vZZaxPA14CFQL27YUkY5E68yvWXKKZEIsGWLVvYsmULiUSi6N/fS7HsytZ7gJIPOS/NF3n0jqsWA3Dg5DC/VMfzKXN7vAwzXXv/8Hry8Q/ArxljGo0xZcD7gMestSen8sXGmHpjTKsxpjWVShU0UAmeXNF59+A4p0Y0kLklV+/RUlOmU8fkfDRf5NElC2u5eEENAHc8dcTlaEQkSKacfBhjlhtjPmiMeYcxpvKMz1UbY/4r/+HxOBAFuoAh4C3A707j6z8MtAFtXV1d+Y9OAq11UiftPar7cI2Kzf1H80Uw3JJd/fjpjk66B8ddjkZEgmJKyYcx5mpgK/AXwFeAncaYSyc9pBx4dz4DM8ZEgAeAPUANUAH8DfCoMaZ5ik9zG7AaWN3U1JTP8CQEaiviNFaVAmo26CYVm/uL5ovgeP2GFuaUxkimLXdtanc7HBEJiKmufHwG+C7QBMwD7gUeNMZsLFRgwFxgGXCbtXbAWpuw1n4VJ+ZrpvIE1toea+0ea+0enfksM5Fb/dir5MMVqXSGthNKPnxG80VAVJbGePPGBQDc+fQRMhnV7YvI7E01+bgM+Jx1jFhrPwh8AXjgjDtaeZPdp7sH+IAxptIYEzPGvAfnNJPnC/E9Rc60qsmp+9C2K3ccPDlMIuWcXrJ2ftUFHi0eofkiQHKF5+29ozy6b0rlMyIi5zXV5CMKvKDS01r718CXcJa6C3WW+ptw7mYdBnqADwJvtdYeKND3E3mB1okTr7Ty4Yad2XqP8pIoS+orL/Bo8QjNFwGydn41ly+pA+COJw+7HI2IBMFU15Z3A1cC2yd/0Fr759lTRb6f78Cyz78LeF0hnltkKlbPc7ZdnRxK0D04PlEDIsWRq/dYPa+KaERntvuE5ouAeceLFrPpcB8P7u7ieP8o82vK3Q5JRHxsqisfdwO/drZPWGs/AdwO6JWBBM6aedXk+hTtONZftO8bi8VYunQpS5cuxe39527GopOufEnzRcDcfMl8aitKSGcsdzypY3fPx0tjd9jo2vvHlJIPa+3nrLXnvKNkrf2wtdbrPUNEpq2yNMay7Haf3BagYohGo9TX11NfX080Gi3a9/VaLLnk4yLVe/iG5ovgKSuJ8rYrndqPbz99hLFk2uWIvMtLY3fY6Nr7x4wnAGPMB4wxuh0pgXdRi/NrvuNY8ZIPgZ6hcbqyvQW08uFvmi/8713XLCEaMfQOJ/jR1mNuhyMiPjabu0+34RyjKBJo61qcLr87lXwUVa7eA2CNkg+/03zhcwtqy3n1OqdlytcfP4S1OnZXRGZmNsmH9uxKKORWPg6eHGZoPFWU7zk+Ps6mTZvYtGkT4+PudhZ2K5adx50am8VzK5hTqv27Pqf5IgBuvW4Z4GyHfPpgr8vReJOXxu6w0bX3D+27FbmAdS2n77rvKmLdR9id7myueg8RL7hiSd3EePj1xw+5G4yI+JaSD5ELaJhTSnO1c8Tujo7inXgVdjrpSsRbjDETqx8/3dlJe++IyxGJiB8p+RCZgom6D618FMV4Ks2+LqervJIPEe94/Yb5NMyJk7HwLTUdFJEZUPIhMgXrdOJVUe3rGiKVcQpaL1LyIeIZpbEo73iRc+zunU8fYSRRnDo4EQkOJR8iU5BLPvacGCSRyrgcTfDl6j2qSmMsrFM3ZREveefVSyiJGgbGUty9ucPtcETEZ5R8iEzBRfOdbVfJtGVv1+AFHi2zlav3WDO/CmN0UJKIlzRVl3HzxfMBuP0JHbsrItMzm+TjL4GT+QpExMsWzS2nqsw57lVbrwpPxeaBo/kiYH47W3i+r2uIx/bpRysiUzfjw/OttX+Zz0BEvMwYw0Xzq3nqYG9Rmg3GYjEWLVo08b6bih2LtVbJR8BovgieSxfVsnFxLc8dOcXXHz/E9asa3Q7JE7w0doeNrr1/6KcjMkXrWmqKlnxEo1GampoK/n2motixnBgYp28kCSj5EPGyW69bxnNHnuOh3V0cPDnMsoZKt0NynZfG7rDRtfcP1XyITFGu6Hzn8QEyGe1xLpTcqkfEwOpmNRgU8arXrJ830QPpa48dcDkaEfELJR8iU7RugZN8DI2nOKLmWgWT66WytKGS8njU5WhE5FxKohHek639uOvZo5wcGnc5IhHxAyUfIlO0onEO8ZjzJ1PoovPx8XG2bNnCli1bGB93d0Ivdiyq9xDxj7dftZiq0hjjqQzfeOKQ2+G4zktjd9jo2vuHkg+RKSqJRia2Ae041l/w75dOp0mn0wX/PlNRzFhyyYeaC4p4X3VZCe+42mk6+N+/PMzwuJoOemnsDhtde39Q8iEyDZPrPiT/xpJpDp4cBpR8iPjFe65bRjwaoX80yXeeaXc7HBHxOCUfItOQSz7U66Mw2joHydXya9uViD80V5fx5o0LAPjaYwdJpjMuRyQiXqbkQ2QaLmpxOp13D47TNTjmcjTBk9tyVVdRMnGKjoh43/teuhxjoOPUKPc+f8ztcETEw5R8iEzD2vlVGOO8r9WP/NvW4dTSXNRSjcldaBHxvBWNc3jV2mYA/v3nB7BWx5GLyNkp+RCZhop4bKKRVjGaDYbN1qOnANiwsNblSERkut7/0hUA7O4c5JE93S5HIyJepeRDZJrWZbdeKfnIr7Fkmt3HBwHYsEjJh4jfXL6kjhctnQvAv/98v8vRiIhXKfkQmabTReeFO243Go3S0tJCS0sL0ai7jfaKFcuOYwOkstXmlyr5EPGl9790OQBPHuhlS/spl6MpPi+N3WGja+8fnk8+jDGvNMY8aYwZMsacNMZ82e2YJNxyycehnhEGx5IF+R6xWIz58+czf/58YrFYQb6H12LZmn2hMq+6jObqsoJ9HwkuzRfue9nqJlY1zQHCufrhpbE7bHTt/cPTyYcx5gbge8DngXpgIfBVN2MSyW27AtiV3SYkszdR77Go5gKPFPlVmi+8IRIxE7Uf9+3o5ED3kMsRiYjXeDr5AD4LfMVa+z1r7bi1dsxau3mqX2yMqTfGtBpjWlMpdV2V/JhbGWd+jXNnvhidzsMit/Kheg+ZIc0XHvGGDS3MrynDWvhKCFc/RPzg4d1djCbc6Qbv2eTDGFMJvAiIGWM2Z5fQHzHGXDGNp/kw0Aa0dXV1FSROCadc9+1CHbc7Pj7Otm3b2LZtG+Pj4wX5Hl6K5dRIgkM9IwBcqpOuZJo0X3hLPBbhd693aj/u3txBe++IyxEVj5fG7rDRtZ+67R39/M43nuHV//wL9ruwOunZ5AOow4nv7cBvAy3AT4H/M8ZM9dXJbcBqYHVTU1MhYpSQytV9FPLEq0QiQSKRKNjzT0ehY9l61FlBMgbWL9S2K5k2zRce846rFtNYVUoqY/nXh/e5HU5ReWnsDhtd+wtLpjP88feeJ2MhYmBBbXnRY/By8pHbTP91a+3z1toEzrJ6CXDtVJ7AWttjrd1jrd2j4iPJp1yn871dgyRSGZej8b/clqsVjXOoLitxORrxIc0XHlNWEuX9L3FWP7636WioVj9EvOxrjx1k53HnxunfvuViykqKfzKYZ5MPa20/cAg4s02qPcvHRIpq/QJn5SOZtuw6rn4fszVR76EtVzIDmi+86ZarltAwx1n9+PIjqv0Qcduhk8N84Wd7AHjblYu4dkWDK3F4NvnI+jJwqzHmImNMDPgjYBx4wt2wJOwW1JbTMKcUgOeO9Lkcjb9ZaydOurpUJ13JzGm+8Jjy+OTVj3Y6To26HJFIeFlr+eTd2xhPZWisKuWTr1nrWixeTz4+D/wX8BBwEngN8JrsXS4R1xhjuGyxc5f+uRA20sqnjlOjnBxy9ujqpCuZBc0XHnTL1Yupr4yTTFu+HLLaDxEvuevZo/zyQA8Af/WGddRUuLfF2dPJh3V82lo7z1pba619mbV2i9txiQBctqQOgM1a+ZiVre3Oa8N4LMKaedUuRyN+pfnCmyriMd6XXf347rPtHNPqh0jRdQ2O8dc/3gnAq9c185qL57saj6eTDxEv25i9S9/eO0r3oI71m6nclqt1LdXEYxqSRILmXdcsYW529ePfVPshUnR/+aOdDIylqCqN8VdvXO92OEo+RGbqkoW1RCMGyH/dRzQapbm5mebmZqLR4p9EUcxYtqjYXCTQKuKxib4f33mmneP9wV398NLYHTa69mf30x2d/HjbcQA++dq1NFeXuRyRkg+RGSuPR1k7vwqAzUfyW/cRi8VYuHAhCxcuxO1jPwsZSyqdYVu2x8elqvcQCazfumYJdRUlJNIZvhLg1Q8vjd1ho2v/qwbGkvzZPdsBeNGyubztykUuR+RQ8iEyC5ctduo+dOLVzOzrHmI0mQZUbC4SZJWlMd6bXf2485l2TgyMuRyRSPD97Y93cWJgnHgswmffcjGR7G4Ntyn5EJmFjdkTr54/2k8qrWaD05Xr71FdFmNpfYXL0YhIIb372qXUVpSQSGV08pVIgf18Tzf/80w7AB955SpWNM5xOaLTlHyIzEJu5WM0mWZ35+AFHj11iUSCnTt3snPnThKJRN6e12uxTNR7LKrFGG/ckRGRwphTerr249tPH+FIT/C6nntp7A4bXfvT+keTfOJ7zwPO/Pq+7N+dVyj5EJmFxXMrmFsZB/K79cpay+joKKOjo1jrboPmQsaypV31HiJhcut1S2msKiWZtnzhgT1uh5N3Xhq7w0bX/rS/vncnnQNjxGMR/vGtlxCLeuvlvreiEfGZyc0G8110HnQjiRR7TjirRTrpSiQcKuIx/uAVqwD44ZYOdh0fcDkikWB5eHcXd206CsDHb2xlZVOVyxH9KiUfIrO0UUXnM7Lj2ADpjHN36pJFNS5HIyLF8ptXLmJpfQXWwufu2+12OCKB0T+S5E/udrZbXba4lt95sbe2W+Uo+RCZpVzR+aGeEXqG1GxwqnLF5gtqy2mqcv/ccREpjpJohI/duBqAh9u6eepAj8sRiQTDX927kxMD45TGIvzDWzdM9CLzGiUfIrO0YWEtub/vXAG1XNjpYnOteoiEzc0Xz2f9gmoA/v6+3aHfoy8yWw/sPMH3Nzvbrf7o1as9dbrVmZR8iMxSZWmM1fOcSXSztl5N2daj6mwuElaRiOGPX70GcOrlfrbzhMsRifjXqZEEn/zBNgCuWFLHrdctczmi81PyIZIHE0Xnh7XyMRU9Q+O0944Cai4oElbXr2rg2hX1APzD/W0TNWAiMnXWWv7snh10D45TVuLt7VY5Sj5E8iBXdL716Km8TKCRSISGhgYaGhqIRNz9My1ELM8fdY7YjRi4eIG2XYmEkTGGT9zkrH7s7Rri7uyWET/z0tgdNmG99vdsOcb/bj0GwJ/ctIZlDZUuR3RhMbcDEAmC3MrHSCJNW+cgF7VUz+r5SkpKWLJkST5Cm7VCxJKr91jVVEVlqYYhkbDasKiW16yfx0+2d/KFn+3h9RtaKCuJuh3WjHlp7A6bMF77o30j/Nk92wF4SWsjv3XNUncDmqLwpIYiBbSsoZLaihIAnmtX3ceFPHu4F1BzQRGBj796NdGI4Vj/GN968rDb4Yj4Qjpj+dh3tzI4lqKuooR/+PVLiHh8u1WOkg+RPDDGsHGR6j6mIpHKsOmwk6BdvWKuy9GIiNtWNM7hN65YCMAXH9xL73DC5YhEvO8/Hz3AUwedG3mffcvFNFf758h6JR8ieXJZrtlgHlY+EokEbW1ttLW1kUi4OxHnO5bnj55iLJkB4Kpl9bN+PhHxv4++spXKeJSBsRRf+Nket8OZMS+N3WETpmu/vaOff/xpGwBvvXwhN62f73JE06PkQyRPckXnB7qHOTUyu4HPWsvQ0BBDQ0Oun3+f71hyd2oWz62gpbZ81s8nIv7XVF3Gh16+CoA7njrM7s4BlyOaGS+N3WETlms/lkzzke9sIZm2LJ5bwZ+/YZ3bIU2bkg+RPNmwqAaT3W753BFtvTqXJ7PdjK9api1XInLae168lMVzK8hY+Kv/3RnoF5AiM/V3P9nNvq4hIga+8JsbmOPDQ1uUfIjkSVVZCa1NVQA8p2aDZ5VMZ3j2ULbeY7m2XInIaaWxKH9681oAntjfw0/VeFDkBR5p6+L2Jw4B8MGXreTyJf68iafkQySPLluSLTrXysdZPX+0n9FkGoCrlvtz0BSRwrnxouaJxoN/8+NdjKfSLkck4g1dg2N8/K6tAFyysIbff8UqlyOaOSUfInm0cZFT97GlPT/NBoMmt+VqYV05C+sqXI5GRLzGGMOnX38REQNHekf4r8cOuR2SiOvSGctHv7OFk0MJKuNR/uVtGymJ+vclvH8jF/Gg3MrH0HiKvV2DLkfjPblic51yJSLnsmZeNe+4ajEAX3poL10DYy5HJOKuf3tkH4/vc27e/e1bLvZFF/PzUfIhkkfLG+ZQl202+OT+Hpej8Ran3sNJPq7WlisROY8/fNVqqstiDCfS/MP9bW6HI+KaZw718k/Z46d/44qFvPHSBS5HNHu+SD6MMRFjzBPGGGuMWeh2PCLnEokYrsnuV358FslHJBKhrq6Ouro6IhF3/0zzFcv2jn5GEs7+bRWbS6FovgiGuZVxPvqqVgDu2nSUre3+qKPz0tgdNkG89n3DCX7/zufIWFjZNIe/8OGxumfjl/O5PgqMzOYJrLWMj4/nKRyRc7tqSS3/t62TJw/0MDwySmyG+zIXLHDubmQyGdd/d/MRy2N7ugBoqSmjsSLi+v+TBJbmi4B468Z5fOvJw+zvHubP79nOnb9zBZGIcTusC/LS2B02Qbr21lo+9t2tHO8fozQW4Qu/vp6oTTM+7v9DGDyffBhjWoEPAL8GPDfNr60H6gFaW1vZvn17/gMUOUN9OgXA4FiKHz22hdb6uMsRecND250tV6tqYceOHS5HI0Gk+SJ43rEmzme6h9lytJ9//t9nuHGFDqqQcLh37zAPtTm1o7+9YQ6J7kNs73Y5qDzx9LqUMSYC/BfwcWAma64fBtqAtt7e3nyGJnJO8yqjNFQ4f1rbumbX6Two0hnLrpNJANY1KhmT/NN8EUyXzivlukVlAHxz2yB9Y/6/6ytyIfv7knxzq5N4XLuwjFctK3c5ovzy+srHHwCd1tofGGOWzuDrbwO+DdDS0tK2fv36PIYmcm4v2b+Du587zoHhEmbye5dMJuno6ACcZeSSkpJ8h1jUWJ4/2s9YymkY9ubr1rN4ru5eSt5pvgiov1syzmu/9EsGx1L84GCEf3qrd382Xhq7wyYo175/NMkf/PvTpCwsqivni++6iqoyr79cnx7P/t8YY1YCHwOumOlzWGt7gB6AK664gtLS0jxFJ3J+L13dzN3PHWdzez82EqOsJDrt5xgaGgKgpKTE9d/d2cay+ahzB2d+TRkr59VijPf3bYt/aL4ItkWlpXzipjV86ofb+fH2E7z1ysXcsLrJ7bDOyUtjd9j4/dpnMpZPfHsrR/tGiUcj/Ostl9FQ4+9jdc/Gy9uuXgw0AkOVXsQAACAASURBVNuNMSeBzdmPP2+M+YB7YYlcWO7Eq0Qqw6bDfS5H475cc8Grl9cr8ZBC0HwRcO940WI2Lnb6KP3ZPdsZTWj7lQTPbQ/t4+E2p7DjL9+4jksW1rocUWF4Ofn4LrACuDT79trsx28E/tutoESmoqmqjNbmOQA8tu+ky9G4K5XO8OwhJwG7apn6e0hBaL4IuEjE8Nm3XEwsYmjvHeWLD+11OySRvHqkrYt/ftDp5/HWyxfytisXuRxR4Xg2+bDWjlhrj+begM7spzqttUNuxiYyFdeuaADgiZAnHzuPDzA47pwApv4eUgiaL8Jhzbxq3nv9cgD+8xcH2N054HJEIvnR3jvCH/zPFqyFdS3VfOZN6wO9S8CzyceZrLWHrLUmO7GIeN51K53kY1tHP/2jSZejcU9uy1VzdSlL6lVoLoWn+SK4/uAVq1hYV04qY/l/d28jk7FuhyQyK2PJNL93xyb6R5PUlJfwlXdePqM6UT/xTfIh4jdXLZ9LNGLI2NMvwMPoqQPOsaWq9xCR2SqPR/nMm5zTrjYfOcWdzxxxOSKR2fnze3awvWMAY+Cf33Ypi0JwGqSSD5ECqS4r4ZKFNQA8HtKtV+mM5emDTvJx1TJtuRKR2XvZ6iZed8l8AP7u/3Zz7NSoyxGJzMx3njnCd55tB5xVvZd5+BS3fFLyIVJA12XrPqabfEQiEaqrq6muriYScffPdDax7HpBvYeKzUUkPz79+ouoKS9hcDzFH3/veaz1xvYrL43dYeO3a7/pcC9/9sMdANywupHff/kqlyMqHu//dER87NqVzt3+/d3DdPaPTfnrSkpKWLVqFatWrXK9UdJsYsltN2usKmVZQ/DOKhcRdzRVlU1sv3ps30m+9ZQ3tl95aewOGz9d+45To7z/m5tIpDMsqa/gn3/zUiKR8GxLVvIhUkCXLa6jrMT5M3tif/i2Xv1yv/p7iEhhvP6S+dx8sbP96m9/vIvDPcMuRyRyYSOJFO/9xrOcHEpQVRrja+++gtqKuNthFZWSD5ECKiuJcuVSZ7vR4/vCVXQ+lkzzeDbhevFK1XuISH4ZY/jMm9bTMKeU0WSaj313K2mdfiUelslYPvbdrew67hSYf/HtG1nZVOV2WEWn5EOkwK6dVPcx1X3JyWSSAwcOcODAAZJJd4/pnWksj+87yVgyA8DL1oSjiE5EimtuZZzPvuViAJ493MfXHjvgajxeGrvDxg/X/l8e3MtPtjttiD75mjWhnRuVfIgU2HXZu/6dA2McODm1bQGZTIa+vj76+vrIZDKFDK9gsTywqwuADYtqaaoqK1R4IhJyr7qomV+/fCEAn79/D3tODLoWi5fG7rDx+rX/8fPH+ZcH9wLwa5ct5HezDTPDSMmHSIGta6mhuiwGhKfbubWWh3afAOCVIb2zIyLF8+nXX0RLTRmJdIY//O4WkmnvvfiU8Nre0c/H7toCwGWLa/nbtwS7g/mFKPkQKbBoxEzaehWOuo/tHQOcGBgH4BVrm12ORkSCrrqshM/9+gbAGX/+9eF9Lkck4jjeP8p7v/EsY8kM82vK+Mq7Lqc0FuwO5hei5EOkCHJbr355oCcUBZEP7HJWPVpqylg7P3zFdCJSfC9e1cC7rl4CwG0P7WPT4V6XI5KwGxhLcuvXn6FzYIzykij/+VtXaBsySj5EiuLalc7KR/9okh3H+l2OpvAezG65esXa5lAvLYtIcX3ytWtY3lBJOmP5/Tu3cGok4XZIElKJVIbf+9YmdncOEjHwpXdsZP2CGrfD8gQlHyJFsLyhkpYa527Hg9lC7KDq7B9je8cAAK+8SFuuRKR4KuIxvvSOy4jHInScGuXjd3mn+7mEh7WWP7n7+Ymt1p9503ptQZ5EyYdIERhjuHHdPADu39HpcjSFlVv1qIxHuXr5XJejEZGwuailmk+/7iLA2QL6tccOuhyRhM0XfraHuzd3APCBG1Zwy1VLXI7IW5R8iBTJTeud5GN35yAHL3DkrjGGyspKKisrXd+2NN1Ycis7169qDH1RnYi445arFnPzJU7387+/bzdb2k8V5ft6aewOG69c+/95+ghffMg58OCNl7bwR69e7VosXqXkQ6RIrlw6l/rKOAA/2X78vI+Nx+OsWbOGNWvWEI/HixFeXmIZTaR5PHuc8CvW6ohdEXGHMYa/e8vFLKmvIJm2fOjbm+kfLXzjOS+N3WHjhWv/cFsXf/rD7QBcvXwun/v1S5SEnoWSD5EiiUYMN65z9nzetz2YW68e23eS8VQGY9TVXETcVVVWwr++4zLi0QhH+0b5xPdU/yGF89yRPj54x2bSGUtr8xz+/V1XaPX/HJR8iBTRTeudbQDPH+3naN+Iy9Hk34PZI3Y3LqqlYU6py9GISNitX1DDn968FoD7dnTy37887HJEEkS7jg/w219/hpFEmubqUr5+64uoKS9xOyzPUvIhUkTXLK///+zdd3xc1Zn/8c+jLtuyZWTL3ZabDLgbgykJMQQIkEAIhEDoJQmbDtlssmm/sJtsdtM2IWTTKDGhhUASSgoQCMU0GxvbuOHem2TZlixZbWbO7497ZYSQZZWZufdqvu/Xa16W7tyZeXTLPH7uOfecw7Odd9T60dzczJYtW9iyZQvNzanvKtCRzsaSSDiefcu730OjeohIWFxzyhjO9Qf8+K+/ruaNrftT9llh+u7ONEFt+01767j6roVU1zczsE8u9904hxHFhWn7/ChS8SGSRnk5WYeHn+2o+EgkEuzdu5e9e/eSSCTSFV6PYlm+o5rKg96s5mep+BCRkDAzvv/RaYw+pg9N8QQ33buY3dUNKfmsMH13Z5ogtv3OA/VcdecC9tY2UpSfw+9umMPEIZpY92hUfIik2Xl+16vFW/dTUZOaBBiEli5XIwcWUj6kX8DRiIi8bUBhLndcM5u+edlUHmzkU/cuoqE5HnRYEmF7axu56q4F7DhQT0FuFndddyJTR2oSwc5Q8SGSZu+dOIi+edk4B0+t2hN0OEnzjD/E7lma1VxEQmjS0CJ+evlMzLz77r6iG9Clm6rrm7nmroVsrKwjN9v41VUncNJYzWvVWSo+RNKsIDf78EhQTx5lyN2o2HmgnlW7vFnNNcSuiITV2ccP4cvnePMuPL5sJ798YUPAEUnU1DXGuHHe66zaVUOWwU8vm8ncScp7XaHiQyQALV2vXtu4j/11TQFH03MtN5r3y89hztiSgKMRETmyz8wdzwXThwPww6fW8EwvaoGW1KptjHHdbxeyaIs3aMH/XDzt8GSW0nkqPkQCMHfSYPJzsognHP/oBYmv5X6P08sHkZejrxURCS8z4weXTGPqiAE4B1/8/RLW7jkYdFgScgcbmrn27oW8vtkrPP7jwsl87MRRAUcVTaH+X4KZfd/MVppZjZntNLM7zEyd6iTy+ubncHr5YODos52HXcXBBuav82Y1P/t4jXIl6adcIV1VmJfNb645gcFF+dQ1xfnEPYuoqm0MOiwJqZqGZq65eyGL/RaP71w0hWtPLQs2qAgLdfEBxIGrgBJgOjASmBdkQCLJct4Ub9z5l9bvpabhnWOSmxkFBQUUFBQEfvP20WJ5dMkO4glHv/wcPuCPpS+SZsoV0mXDBhTy66tPIC8ni637DnHDvNepa4z16D3D9N2daVK17avrm7n6roUs2XoAgO99ZCpXnzwmae+fiUJdfDjnvu6cW+Kca3bOVQK3AXM7+3ozKzGzcjMrj8V69oUikmzvP24IudlGc9zxT3+kqBZ5eXlMnjyZyZMnk5eXF1CER4/FOccfFm0H4ILpw+iTlxNEiJLheporQPkiU80aPZAfXzodM1i2vZp/uW8xTbHuzxERpu/uTJOKbV99qJmr71rAsm0HMIP/uXgqV8wZnZT3zmShLj7a8X5gWRfW/zywBlhTUVFxtHVF0mpAYS6njh8ERLfr1ZJtB1hfUQvApbPV91VCo6u5ApQvMtYF04dz6wWTAZi/bi//+vAyEgkNwZvp9tc1cdVdC3hzezVm8P1LpnH5SSo8kiEyxYeZXQL8C/DFLrzsdmASMKm0VMOgSfi0dL16YW0lh5qid7X14UXbAJhQ2o+Zo4oDjkak27kClC8y2rWnlvGFMycA8MSynfzHEys1B0gG23mgnkt//SrLd3iFxw8/Op2P6QJb0kSi+DCzS4E7gAudc2909nXOuSrn3Frn3NqcHHUHkfA5+/ghZBk0NCd4fk3l4eWxWIzt27ezfft2gu4CcqRYDjXFeGKZ12Lzsdkj1b9ZAtfdXAHKFwK3nF1+uEvNPa9u4ef/XN/l9wjTd3emSda2X7fnIJf88hXWV9SSk2X89LIZfPSEkUmMVEJffJjZ9cCvgQucc88FHY9IMpX0y+eU8d68GL9/fdvh5fF4nD179rBnzx7i8XhQ4XUYy9+X76a2MUZ2lvGRmfpilmApV0hPmRnf+fCUwy3SP/7HWu5fsKVL7xGm7+5Mk4xt/8bW/Vz661fZVd1AYW42d147mw/PGJHkSCXUxYeZfQH4EfAB59zLQccjkgpXzvFGzXhxbSUbK2sDjqbz/uB3uTrz2FIGF+UHHI1kMuUKSZbsLOMnl83g5HHeSM3fenQFjy3dEXBUkg7PrangyjsWcOBQM8V9cnngk3M0c3mKhLr4wBuxpD/wnJnVtjyCDkokmc4+fghD+xcAcO9rXbvKFpTNe+tYsGkfgPrBShgoV0jSFORmc8c1s5k8vD8JB7c8tJQ/Lt4edFiSQo8u2cEn71lEfXOc4QMKeORfTmHm6IFBh9Vrhbr4cM6Zcy7XOdev9SPouESSKTc7iyv9fsaPLNre43Hm0+ERPxEP6pfP3EmDA45GMp1yhSRbUUEu990453AB8uVHlvHQ61uDDkuSzDnH/z23npsfWkos4ZhQ2o9HPn0qE0qLgg6tVwt18SGSKS4/aTS52cbBxhiPhryJP55wh4uPS2aNIDdbXyMi0vsM7JvHA584mekjB+AcfPWPy7kvIq3TcnQNzXG+9Idl/PCpNQDMGFXMwzedwvDiwoAj6/30vwaREBhclM8Hpw4D4HevbAn1EI/z11Wyu6YBgEtn60ZzEem9BvTJ5d5PzGHWaG8o8W8+uoJ5L28KOCrpqcqDjXz8jtf48xLvYt+Hpg3jwU+ezMC+mhgyHVR8iITE1aeUAbBmz0Fe33Ig2GA68LA/o/ms0cVqmhaRXq9/QS6/u3EOJ5Z59wDc+sQq7py/MeCopLtW7azhwz9/iSVbvTx7y1nl3P7xmRTmZQccWeZQ8SESErNGFzNlRH8A7luwjdzcXHJzcwOOytMSy766Jp5etRvQjeYikjn65ecw7/qTDo+C9d2/ruZ/n17Tbit1mL67M83Rtv1TK3fz0V+9ws7qBgpys/i/K2bxxbMmap6qNFPxIRISZsY1fuvHM29VMnjMRKZNm0Z+frDD2Obn5zNt2jSmTZvGk6v30hx3FOZm88FpwwKNS0Qknfrm5/Db607ivRMHAfCzf67nS39YRmPs7TklWn9fBv3dnWk62vbxhOOnz6zlX+5bzKGmOEP65/PwTacqjwVExYdIiFw4fTjFfXKJJxwPLAjXyCrOOR7yJ0I8f+owigp0ZU9EMkthnjfx3EUzhgPw5yU7uPquhRw41BRwZHIklQcbufbuhfz0mXU4B9NGDuDxz72HqSMHBB1axlLxIRIiBbnZXHai153pwYVb33FFLWhPrdzDW7sPAnD5SepyJSKZKT8nm59cNoMvvH8iAAs37ePiX7zClqq6gCOTtl7bWMUHfzafl9bvBeDjJ43iDzedwhB/bi0JhooPkZC5as4YzGBvbRMPzn+LWCzYeT9isRjbd+zk+39bCcB7JgzixLJjAo1JRCRIZsaXzi7nR5dOJyfL2Li3jo/84hUWbtzLrl272LVrV+Df3ZkmFosd3vZNTc3833PrueKO16g42EhhbjY/uWw6/33xNApydWN50FR8iITMqGP6MLfc61P84KIdxOPBtn7E43Eeem0jm/Z5w+t++QOTAo1HRCQsPnrCSH53w0kUFeSwr66Jq+5+nd+/uoGdO3cG/t2daeLxODt37mT1pu3c+LvF/PCpNSQclA/pxxOfP42PzNTQ8GGRE3QA6eKco7GxMegwRDrlsplDeW7NXtZUNbNk815mlg0KLJbaQ/U8tLIWgPeXl3BcaaHOJenVlC+kK04YVcRDnziRT963hB0HGrhtYTUrK5v4XtkhdFdB+jQ2NrJgRwO/XlxDdWMCgItnDONbHzyWPnnZOqdDJGOKj/r6elasWBF0GCKdUuIcw/tls7M2zk+eXs1XTh0YWCxPrj9ExaE4BnyoDJ1H0uspX0h3fOe9RfxkQYLlFU08s6metb9ZyJdPKWZEUcb8VyswB5sS3L2khhe3ei30BdnGjTOLOHOsY+Pa1QFHJ22p25VICGWZcfFx/QBYsKORJbuDuWLTEEvw8Gqv1eP0MQWMHqARrkRE2jOgIJtvnT6Qyyb3w4Ct1TG+8o8q5m+tDzq0Xm3RzgZufmrv4cJj8uBcfnxOCWeO7RNwZHIkGVOOFxYWMmXKlKDDEOmUxsZGEm4t/9h4iDVVzfxuZQMfe98M8tN8o9xv5m/mQEMFOQaXTe7HxIkTNXa99HrKF9JdjY2NfMzWcdygXG5fVEtVXTM/XVDNzlg/vnFeuW52TqKa+ma+9+Ra/rzUm6m8ICeLK6b05bwJfZhUXq5cFWIZU3yYmQ5EiZQsMz45qz9feaaKrfvqmbdgx+GhHdOhur6ZO1/eAsBZ4/owpG8O+fn5Oo+k11O+kJ6aWprPY5+ewlcfXcXL66v4w+IdLNp6gP+5eBonjdVogT2RSDj++MZ2vv/kW+yt9eZXOWHMQL734eOo3bURQLkq5NTtSiTExhbncqU/p8b/PbeebfsOpe2z73hxI9X1zRTkZvHR4/qm7XNFRHqDwUX5/O6GOdxyVjlZBhsr6/jYr1/lG39eTk1Dc9DhRdLy7dVc8qtX+LdH3mRvbRP5OVl84/zj+MNNp1BWom5WUaHiQySksrKyyMrK4otnjmdQv3waYwn+44lVafnsyoON3P3yJgCunjOKkr65ZGXp60JE5GhavrsBsrOML541kT9/5jSOHVoEwP0LtnL2/77A0yt3BxlmpOyra+Jrf1rOhf/3Eku2et2sPjB5CM986X188vRxZGcZ8M5tL+FlzrmgY0iL2bNnu0WLFgUdhki3/HnJdm55aBkAd107m/cfNySln3fr4yuZ98pmigpymP+VMyjuk5fSz5NwMLPFzrnZQccRNOULSYXmeILfvLiR255dR1PMGwr2/KlDufWCyZRqxu12NTTHeWDBVm57dh3V9V5r0bjBfbn1gsmcXj444OgyW0/yhcpDkQi4aMYITvJnFb/1iZU0NKdu8qrXN+/j/gXevR43nT5OhYeISBLkZmfx2TMm8OQX38sc/76Pvy3fzdwfPc+Pn16jrlitNMbi3PvqZub+8Hn+8y+rqK5vpm9eNl8//1ie/OLpKjwiTsWHSASYGf950WSys4xt++r5xfMbUvI52/Yd4qZ7F9Mcd5SV9OH608am5HNERDLVuMH9ePCTJ/M/F0+lf0EOh5ri3P7P9bzvB89x5/yNKb24FHZNsQQPLNjKGT98nm89tpLdNQ1kZxmXnjCSf355Lp86fTx5Ofqva9RlzGhXIlESj8epqqoCoKSkhOzsbI4d2p/rTy3jzpc28asXNnDxzBGUDUrejeAHG5q58Z7X2VfXRFFBDndddyJ983PajUVERN6ts9+XWVnG5SeN5twpQ/nlCxuY9/Jm9h9q5rt/Xc3dL23i5rPLuWTWyMP3MvR29U1x/rxkB794fj3b93vzomQZXDRzBF84c2Kncp1yVXSo+BAJoVgsxrZt2wAYMGDA4S/RL541kceX7aTiYCOfuf8N7r3xJEr69Xw4wXjC8fkHl7B2Ty3ZWcYvrzyB8YP7dRiLiIi8U1e/L4v75PG1847j+lPHctuza/nDou3srG7gK4+8yS+f38C1p4zhkhNGUlTQOyd43VJVx32vbeGh17dR0xADwAwunD6cL7x/4uE81BnKVdGhtiuRCCkqyOU7F03BDFbtquFjv36VXdU9nz33v/66mufXVAJw64WTec/EQT1+TxER6ZyhAwr474un8fQtp3P+1KEAbNpbx61PrOLk7z3Ltx9bwYbK2oCjTI5EwvH8mgpumPc6c3/0PHfM30RNQ4ycLOPC6cN5+ubTue3ymV0qPCRa1PIhEjEfmDyUn10+k1seWsqGyjou/dWr3P+JOYwp6V4XrAcWbD08rO61p4zh6pPHJDNcERHppPGD+/GLK09g5c5q5r28mceW7aSuKc49r27hnle3cHr5YK6aM5r3TRpMfk50ruw751i5s4a/vLmLvy7fybZ9b180G1yUzxUnjeaKOaMZolG/MoKKD5EIumD6cPrmZ/Pp+95g+/56Lv3Vq9z3iTmUDynq0vu8smEv/++xFQC8d+IgvvWh41MRroiIdMHk4QP44aXT+ffzjuX3r2/j3le3sLumgRfXVvLi2kr65efw/uNKOW/KUN5XXkphXvgKkZaC42/Ld/HX5bvYUvXOSXJnjxnINaeWce7kobqJPMOo+BCJqDOPHcK860/iE/e8TsXBRj7261f53Q0nMW1k8VFfe+BQE799eTN3vbSJWMIxfnBffn7FLHKylQBERMKipF8+nz1jAp86fRxPr9zDPa9uZuGmfdQ2xnhs6U4eW7qTwtxszjh2MGcfP4TZY45h5MBCzIK5UX37/kMs2LiPBZuqeHVj1TtaOABGDizkg1OHceGM4UwePiCQGCV4oS4+zCwb+B/gOqAAeBq4yTm3N8i4RMLilPEl3P/Jk7n27oUcONTMFXcs4OvnH8cp40soK+nzrgRUcbCBu+Zv4r7XtlDX5A3neEzfPO6+7kQGFPbOGxolMyhfSG+Wm53FB6cN44PThrG7uoEnV+zi7yt28/rmfdQ3x/nb8t38bbk3Y3ppUT6zRg/khDEDmTWmmMnDB1CQm/yWkQOHmthQWcvaPbW8vnkfCzbuY8eBd9+DOKK4kA9NG8b5U4cxbeSAwAojCY9QFx/AvwMfBuYAVcDdwL3AeUEGJRImM0YV89BNJ3PVnQvZW9vI1/+8HIBB/fKYNXogs8sGMmX4AJ5auZvfv76NRn9m3T552Vx18hg+8d6xlBapn61EnvKFZIShAwq47rSxXHfaWCoPNvL0qt08uWI3CzftozGWoOJgI0+u3M2TK71iJMtg2IBCRg4sZNQxfbx/B/ZhWHEBffJyyMvOIj83i/ycLPJzssnJMmobY1TXN1PT0ExNfYya+mb2H2pic1UdGyrq2FBZS1VdU7vxFeRmccKYgcwZW8Lp5YOZroJD2gh78fEp4D+dcxsBzOwrwHozG+Oc23K0F5tZCVACMH369JQGKhKkY4f255F/OYX//MsqFmysoq4pzt7aJp5etYenV+15x7pFBTlcf2oZ1582loF9NXu59BrKF5JxBhflc+WcMVw5ZwxNsQSrdtXwxpb9LN66nyVb9rOzuoGEgx0H6tlxoJ4Fm/YlPYbiPrlMG1nMnLHHcPK4Y5g6olj3cEiHzDkXdAztMrNiYD8w0zm3tNXyauBq59zjnXiPW4Fv+78eAlZ3I5RsYAiwB4jCtKNRixeiF7PiTS3Fm1odxTvGOTc4/SH1jPJFtyne1FK8qRW1eCF6MackX4S5+BgFbAXGOec2tVq+BfiGc+6+TrzH4StZQJVzrqobcZQDa4BJzrm1XX19ukUtXohezIo3tRRvakUt3s5QvugexZtaije1ohYvRC/mVMUb5m5XB/1/2w6HUAzUdOYN/OTR5QQiIiKRonwhIhIRoe2U55w7gHcla1bLMjMbB/QH3gwqLhERCRflCxGR6Aht8eH7DfBVMxtrZv2B7wNPOec2pzGGKuA/iM4VsajFC9GLWfGmluJNrajF21nKF12neFNL8aZW1OKF6MWcknhDe88HHB63/ft447bnA/8APqVx20VEpDXlCxGRaAh18SEiIiIiIr1H2LtdiYiIiIhIL6HiQ0RERERE0kLFh4iIiIiIpIWKDxERERERSQsVHyIiIiIikhYqPkREREREJC1UfIiIiIiISFpkZPFhZtlm9kMzqzSzg2b2RzMb1MH655rZSjOrN7MVZnZOm+cnmNkzZlZnZtvN7F+DitfMzjezf5rZXjPbb2bzzey9bdZxZnbIzGpbPQYEFO9cP57WsbzSZp0wbd+vt4m11o//Z63W2WxmDW3WmZrEeC/392uNmcU6sf5sM1vo7/MNZnZVm+dLzexP/t9eaWbfN7OkfTd0JV4zO9nM/mpme8ys2swWm9lFbdYJzfY1szJ//9e1imV7m3XCtH2vbOf4jZvZ463Wed7MGtus86FkxRs1Xfl+8NdXvkhdvMoXXYtVuSKFuaIbMStftHDOZdwD+AawFhgHDAD+CPz9COuOAw4BVwF5wJVAHVDmP58NrAZuB/oAs4AK4LKA4r0S+AhQDOQAnwZqgVGt1nHAe0KyfecCsQ7eK1Tbt53XlgMJ4KRWyzYDV6Vw+34A+DhwQ0fbzl93AFAJfBVv1uez/ePhlFbr/AP4k7/uOH9bfDWgeM8HrgEG4V0cuQioB04M6fYt88+nkR2sE5rte4Tjow74WKtlzwPfTNX2jdqji99nyhepjXduR8d42LZvO69Na77o4neZckXqYy5D+cJ7XSp3SlgfwBbgxla/j/cPiDHtrPsfwPw2y+YD3/Z/PgMv2fRr9fx3gOeCiPcIr98NXNzq91Qnk65s37kdnQBh377Aj4DFbZal/AuvM9vOX+d6/++zVsvuBX7r/zzW/1vHt3r+RmBTEPEe4XWvAV8K6fbtMJmEffsCn/O/H3JbLetWMumtD+UL5Ytkbd+g8oVyReq2bRe3sfKF/8i4bldmVgyMBha3LHPOirMHngAAIABJREFUbQBqgOntvGR663V9b7Radzqw1jlXe4Tn0x1v29dPxbsysLzNUw/7Te0LzOziZMTag3izzWybme32m1Fbrxfa7Wtm+cB1wK/befp/zWyfmS01s5uSEWs3TQeWOP9bwtf2+K32/+bWz5eZWf80xXhEZjYUmAwsa/NUWLZviwV+E/nzZja31fJQb1/gJuBu51xzm+U3+9t3pZl9zcxygwguaMoXhylfdD/elteGPV8oV6RPxueLjCs+gCL/3+o2yw8A7e3coqOse7Tne6qr8R5mZqV4TcI/cs6ta/XUWXgV9kjgf4H7zezc5ITb5XjfAmb48RwLvAn808yGt3q/UG5f4KN4XSseaLP8Wrzm0iHAvwHfC/BLr7vHLyRvG3eLmfXFO37/6px7ttVTYdq+e4FT8I7fMvwuGGY2zX8+zNv3NOB44I42T30NmAgMxrvq9gngP9MbXWgoXyhftNab84VyReopX/gysfg46P/b9oa5YryrF+2t39G6R3u+p7oaLwD+l/FzwNN4B8dhzrlnnXMN/uMh4D68vr9pj9c5t9s5t8w5F3POHXDOfQ3YB5zX6v1Ct319NwH3t7nKhnPuBedcrXOu2Tn3D7yEfVW775B63T1+W54LhJkVAX/H6699TevnwrR9/Thec841OefqnHO3Ay8Bl/qrhHL7+m4CnnbObWq90Dn3qnNuv3Mu7px7Dfh/BHf8Bk35Qvmi2/G2EfZ8oVyRYsoXb8u44sM5dwDYincjGgBmNg6vqnyznZcsa72ubyZvN+0tA8r9yru959MdL2ZWhtfP+O/Ouc+1aUZtTwKwoOI9Sjyh277+OscD7wV+1YmPSdr27YZleFcKW2t7/A7w/+bWz292zrW9ApMWZlYCPAvsBC51zjUd5SVBbt/2tD1+Q7V9AczsGLyEF/bjN1DKF+1Svuid+UK5IhiZmS+6epNIb3jgjVaxBq/pqz/wMPDkEdYdj3cD28eBXP/f9kYvuQ0oxDt59wCXBxTvscB24LtHeH4KcBJe828u3ggRh4ALA4r3TGACXiHcD7gVr5lxVBi3b6vX3Aa82s7yMXg3PRb4sb8P74rM55MYb7b//ucAMf/nAlrdKNhq3WK8EUz+zd/n76f9EUwe8f/2sf62+PeA4h0KrADuAbIjsH1P9s+pHH+dTwENwAlh3L6tXnMLsK3tNvaPlw/556LhJb41wI+TFW/UHl35fkD5ItXxKl90LVblihTmim7ErHzR8tpk7oSoPPyN/yO8/ncH8YY1G+Q/dyVQ22b9c4GVeMO4rQTOafP8BLzq+xBeBf7loOIFfos3WkJtm8eV/vNn+H9DHbAfWEQSv5i7Ee8teCNs1PlfDE/Saqi8sG1ff1mhv+2ubee9jgeW+O9Tg/fl+Lkkx3udv4/bPsrwrq7VAqNbrX8isNA/fjfSZvQPoNT/mw/62+AHQFYQ8QLf9p+ra3P8fj2M2xfvP5fr/Xir8K4gnx3W7dvqNavxR2Bqs3ww3ogx1X68a/19kpfMbRylRze+H5QvUhev8kXXYu3SdwPKFamOWfnCf5j/BiIiIiIiIimVcfd8iIiIiIhIMFR8iIiIiIhIWqj4EBERERGRtFDxISIiIiIiaaHiQ0RERERE0kLFh4iIiIiIpIWKDxERERERSQsVHyIiIiIikhYqPkSSwMxmmJkzsyGtlpmZHTCzy4KMTUREwkP5QjKdig+R5JgJ7HTO7Wm1bDwwAHgjmJBERCSElC8ko6n4EEmOmbw7acwCDgLru/OGZnaTmf2yzbIVZnZc90IUEZEQ6DBfmNl4M1ve+kkzyzezTWY22c8Nu8xsqZmtN7NHzSyvzfKlZnZfuv4gka5Q8SGSHDOBJW2WzQKWOudcN99zKq0SlJkVAGXA2m6+n4iIBO9o+WITMNLMWv8f7VPAi865lXi54evOuRlAOTAFmOYv/6Zzbob/uCrVf4hId6j4EOkhMzO8L/62V7JOpFWCMbPjzOxFM3vTzP7NzNZ3tLyd95wKrHXOxVP314iISKp0Jl845xLAVryLTZhZIfCvwLf9dafxdm6ZABjeRalpwNLURS+SHCo+RHpuPNAfL1kAYGZDgffgJxgzywHuB77onJsGjANWHGm5/zaTgT+Z2WYz2wz8HXgzLX+RiIikwlHzhW81cKz/82eBJ5xzm/3fJwO/M7PVwCLgeudcjb/8t36Xq2dS+leI9EBO0AGI9AIz/X8/b2b/DQwDvg/kAVlmlgdcBCxzzrVcrVoFVAAXt7fczEYBlc65luSDmf0crzleRESi6aj5wjnXhFd8TDKzF4HPAXMA/NxQ4V+swsyuAb5lZjcAu1uWi4SZWj5Eem4mMB/oi9cy8Rvgx8AuvKbyZt7dHD7F//1Iy6cCK9t8zvGo5UNEJMo6ky/g7ZaPLwL3txoZayreRaoWy4BS2s8ZIqGklg+RnpsBvOGcu7nN8odbfjCzKrwbAzGzGcBVeFe7yo6w/GO8M8GA16S+HBERiaqj5gvfauBrwFnACa2WT/Ofa7l/5FrgGX+5ig+JBLV8iPTcTI5+k9+9wGx/+MQbgc3OuY0dLH/H1S0zOwYw59zuVPwBIiKSFp3JF+DdQD4V+I1z7kCr5VOB68xsCd79HgXAt3h3i4hIaFn3RwEVEf9GwV3ATOfcUjN7HJgNjPRHLGlZr59zrtbMTgVeBp7Au+nQOefO8Nf5N2CAc+6bZnYncJZzrizNf5KIiKRAO/niImC0c+5nAYcmklZq+RDpAefcbuecOedarmTdi3cD4ZltVr3FzFYCjwEx4AZ/+WgzW2lmS/G6YH0nDWGLiEiatZMvLgK+EGRMIkFQ8SGSXE8A1Xj3bhzmnPsOXl9fgD865/b6P29yzk32J4T6rHOuMY2xioiIiKSVig+RJHLONQCPABf7E0O19gFgEF7riIiIZCgzm4d3s/h4M3P+43n/uVPN7HkzO2RmNWb2uJlNavP6zWZ2p5ndbGZbzazezJ41s4np/2tEukbFh0jy3QsUARe2WX4VUAk81WqZmVlO2wfejLUiItI7fQf4G7AdOMV/fMbMZgL/BHKBK4BPAZOAl8xsWJv3OBe4Bm843hvwJql92szy0/IXiHSThtoVSb4XgS3AlcBDAGbWUozc6ZyLtVr3TN4e172tLakMUkREguGc22BmlUCjc+61luVm9kegBjjHOVfnL3sVWIdXZPx7q7cpAWY55yr89VbizftxLd78ISKhpJYPkSRz3hByDwDnmlmJv/gjQCHv7nK1EDixncdj6YlWRERC5HTgsZbCA8A5twV4xX+utRdbCg9/vTfxipST0xGoSHep+BBJjXvxms0/5v9+FbDGOfd6m/UOOucWtX0AexERkUwzEGhvPqfd/nOtVbSz3h68ERdFQkvFh0gKOOdWA4uBK/2x3c8E7gs2KhERCbn9wNB2lg8F9rVZVtrOekPw5hIRCS0VHyKpcx9wKvB1vHNNxYeIiLRoxJuhvLUXgQtbj5ZoZqPwcsn8NuuebmalrdabBkwEXkMkxFR8iKTOg0AC+BzwknNuc7DhiIhIiKwGhpvZdWZ2oj+c7neBAXijVl1kZpcBT+PNH3Vbm9dXAU+a2UfM7HLgUWAzMC9df4BId2i0K5EUcc7tMbOngfNQq4eIiLzTncAs4Ad4c0C96Jyba2ZnAv+NN3BJDHgeuMg517Y71ZPACuBnwGDgZeDTzrmm9IQv0j3mDcwjIiIiIlFgZpuBZ5xznwg6FpGuUrcrERERERFJCxUfIiIiIiKSFup2JSIiIiIiaaGWDxERERERSYuMGe1q0KBBrqysLOgwpBOcc9TX1wNQWFiImQUckciR9abjdfHixXudc4ODjiNoqcoXGytrqWuKM6hfPsMGtJ3eQaKoN53/vZ32VXL1JF9kTPFRVlbGokWLgg5DOqGxsZEVK1YAMGXKFPLz8wOOSOTIetPxamZbgo4hDFKVL259fCXzXtnMnLHH8NBNpyT9/SX9etP539tpXyVXT/KFul2JiIikwZQRAwBYtbOGREL3W4pIZlLxISIikgZTRvQH4GBjjK37DgUcjYhIMFR8iIiIpMGEwf3Iz/HS7oqd1QFHIyISDBUfIiIiaZCTncWxw7zWjxU7agKORkQkGCo+RERE0mTKcK/4WKmWDxHJUBkz2pVER1ZWFiUlJYd/FgkzHa/SFS03na/YUY1zTsN9RpzO/+jQvgoPFR8SOrm5uWhOFokKHa/SFVOGe8XH/kPN7KxuYERxYcARSU/o/I8O7avwUOknIiKSJuVD+5GT5bV2rNihrlciknlUfIiIiKRJfk42E4cUAbBSxYeIZCB1u5LQaWpqYvPmzYA303BeXl6wAYl0QMerdNWU4f1ZvauGlTs14lXU6fyPDu2r8FDLh4SOc46DBw9y8OBBnNMswBJuOl6lqw7fdK4RryJP5390aF+Fh4oPERGRNGqZ6XxPTSMVBxsCjkZEJL0CLT7M7HIzm29mNWYWa+f5UjO7x8yq/HWWmtnwIGIVEZFg9LZccdyw/rSMsKuuVyKSaYJu+dgP/AK4ue0TZlYAPAs0AZOAYuBKoDadAYqISOB6Va7ok5fD+MH9AN10LiKZJ9Abzp1zTwGY2dx2nr4WL4l8xjnX7C9b2YPPorGxsbsvlzRqvZ+0zyTsdLymXjpzhf95Kd+Xxw3tx/qKWt7ctl/HTYTp/I8O7avwCPNoV2cA64B5ZnYuUAn82jn3k86+gZmVACUA5eXlrFixIiWBSuqsW7cu6BBEOk3HayB6nCsg/fniGKsDYOmWfcpNvYTO/+jQvgpW0N2uOjIIL6ksBIYBVwHfMLMru/AenwfWAGv27duX/AhFRCRoycgVkOZ8Ma7Yu/ZXcSjOwaZEyj9PRCQswtzycRDY4Zy7zf99kZndB3wYuL+T73E78ADA8OHD10yZMiX5UUrSNTc3s3PnTgCGDx9Obm5uwBGJHJmO18AlI1dAmvPF6Ppmvv3CC94vxSOZMu6YlH6epIbO/+jQvgqPMBcfS4HZ7Szv9ODMzrkqoApg9uzZ5OfnJyk0SaX8/HzKy8uDDkOkU3S8Bq7HuQLSny8G5+czpqQPW6oO8VbFIeYeNyylnyepofM/OrSvwiPooXaz/ZFK8vzfC/yHAfOAEjP7rL/edLwRTP4UXMQiIpJuvTVXTBtZDMDSrQcCjkREJH2CvufjaqAeeArI9n+uB8Y457YA5wOfAGqAR4BbnXMPBRSriIgEo1fmihmj/OJjm4oPEckcQQ+1Ow/vqtWRnn8emJmmcCQkmpub2bp1KwCjR49Wv0wJNR2vqddbc0VL8bG7poFd1fUMG1AYcETSVTr/o0P7KjyCbvkQeZdEIsGBAwc4cOAAiYRGgZFw0/Eq3TV5eH9ys72pztX1Kpp0/keH9lV4qPgQEREJQEFuNscP6w+o65WIZA4VHyIiIgFp6Xq1RMWHiGQIFR8iIiIBmTHaKz6Wb68mFldXEBHp/VR8iIiIBGTGqIEA1DfHWbPnYMDRiIiknooPERGRgJSV9KG4jzfqju77EJFMoOJDREQkIGb29nwfGvFKRDJAoPN8iLTHzCgqKjr8s0iY6XiVnpoxqpjn11Sq5SOCdP5Hh/ZVeKj4kNDJy8ujvLw86DBEOkXHq/RUS8vH+spaahqa6V+gyc+iQud/dGhfhYe6XYmIiASopfhwDt7cVh1wNCIiqaXiQ0REJEDFffIYN6gvAEu37Q84GhGR1FK3Kwmd5uZmduzYAcCIESPIzVUXBAkvHa+SDDNGFbNxb53u+4gYnf/RoX0VHmr5kNBJJBJUVVVRVVVFIqFJtyTcdLxKMrRMNrhk6wGccwFHI52l8z86tK/CQ8WHiIhIwFru+6iqa2L7/vqAoxERSR0VHyIiIgE7dmh/8nO8lLxEXa9EpBdT8SEiIhKwvJwspowYAGiyQRHp3VR8iIiIhEBL16slGvFKRHoxFR8iIiIh0FJ8rNxZQ1NMN8SKSO+k4kNERCQEZvojXjXFEqzeVRNwNCIiqaHiQ0LHzOjTpw99+vTBzIIOR6RDOl4lWUYUFzKoXz6A5vuICJ3/0aF9FR6aZFBCJy8vj+OOOy7oMEQ6RcerJIuZMWNUMc+s3sOSrfu59tSyoEOSo9D5Hx3aV+Ghlg8REZGQaOl6pZYPEemtVHyIiIiExEz/pvPNVYfYX9cUcDQiIsmn4kNCJxaLsW3bNrZt20YsFgs6HJEO6XiVZJo6cgAt3dGXblfrR9jp/I8O7avwUPEhoROPx6moqKCiooJ4PB50OCId0vEqyVRUkMvE0n4ALNFkg6Gn8z86tK/CI9Diw8wuN7P5ZlZjZkcsQ83s02bmzOyb6YxPRESCl2m5YtbogQAs2rwv4EhERJIv6JaP/cAvgJuPtIKZjQH+FVierqBERCRUMipXnFh2DOC1fDTHNdmgiPQugRYfzrmnnHMPAhs7WO0u4BuALgGJiGSgTMsVJ431io/65jgrd2qyQRHpXYJu+eiQmd0E1DnnHurm60vMrNzMynVzkYhI79TTXOG/R2jyxciBhQztXwDA65siX0uJiLxDaIsPMxsNfBP4TA/e5vPAGmBNRUVFUuISEZHwSFKugBDlCzPjRL/1Y6Hu+xCRXia0xQdwJ/Bd59yOHrzH7cAkYFJpaWlyohIRkTBJRq6AkOWLk8revuk8kXABRyMikjxhLj7OBr5nZnvNbC9wGvA1M5vf2TdwzlU559Y659bm5OSkLFBJLjMjPz+f/Px8rGXAe5GQ0vEauB7nCghfvmhp+dh/qJkNlbUBRyNHovM/OrSvwiPQb1gzywZygTz/9wL/qUZgVJvVHwbmAz9OW4ASiLy8PKZMmRJ0GCKdouM19TIxV5SXFjGgMJfq+mYWbt7HxCFFQYck7dD5Hx3aV+ERdMvH1UA98BSQ7f9cD4xxzm1v/cBLMjXOuT3BhSsiIgHIuFyRlWXMHuN1vdJN5yLSmwTa8uGcmwfM6+S6c1MZi4iIhFOm5ooTxx7Ds29V8Prm/UGHIiKSNMF3bBVpIxaLsWePd9FyyJAhhKH/tciR6HiVVDnRv+l8x4F6dhyoZ0RxYcARSVs6/6ND+yo8gu52JfIu8Xic3bt3s3v3buLxeNDhiHRIx6ukytQRxeTneGlaXa/CSed/dGhfhYeKDxERkRDKy8lixqhiQPN9iEjvoeJDREQkpE7yh9xVy4eI9BYqPkRERELqxDKv+FhXUcv+uqaAoxER6TkVHyIiIiE1a8xAsvz50F5X1ysR6QVUfIiIiIRUv/wcJg8fAKj4EJHeQcWHiIhIiLV0vVqo+T5EpBdQ8SGhlJOTozG4JTJ0vEoqnTTWm+9j5Y5qDjXFAo5G2tL5Hx3aV+GgPSChk5+fz/Tp04MOQ6RTdLxKqs32Wz5iCceSrQc4bcKggCOSFjr/o0P7KjzU8iEiIhJig/rlM25wXwAWashdEYk4FR8iIiIhd5Lf+qGbzkUk6lR8SOjEYjF2797N7t27icXUv1nCTcerpEPLTedLth6gOZ4IOBppofM/OrSvwkP3fEjoxONxduzYAcDAgQN1c5iEmo5XSYeWmc7rm+Os2FHNzNEDA45IQOd/lGhfhYdaPkREREJu5MBChvYvANT1SkSiTcWHiIhIyJkZJ/qtH7rpXESiTMWHiIhIBJw8zis+FmzcR0z3fYhIRKn4EBERiYDTxnvzexxsjLF8R3XA0YiIdI+KDxERkQgYU9KHEcWFALyyoSrgaEREukfFh4iISASYGaeOLwHg5fV7A45GRKR7VHxIKJkZZhZ0GCKdouNV0uW0CV7Xq0Vb9tPQHA84GgGd/1GifRUOGuRYQic/P59Zs2YFHYZIp+h4lXRqafloiiVYtHk/75k4KOCIMpvO/+jQvgoPtXyIiIhERGn/AiaW9gPg5Q3qeiUi0aPiQ0REJEJaul69ovs+RCSCVHxI6MTjcSorK6msrCQeV59mCTcdr5JuLV2vlu+oprq+OeBoMpvO/+jQvgqPQIsPM7vczOabWY2Zxdo8d76Z/dPM9prZfn+99wYVq6RPLBZj69atbN26lVgsdvQXiARIx2vqKVe805xxJWQZJBy8tlFD7gZJ5390aF+FR9AtH/uBXwA3t/PcQOB2YAIwGHgA+LuZjUpfeCIiEgLKFa0MKMxl6shiQF2vRCR6Ai0+nHNPOeceBDa289z9zrk/O+cOOOdizrlfArXAiWkPVEREAqNc8W6ntcz3ockGRSRigm756DQzmwoMApZ34TUlZlZuZuVqYhMR6f26kyv810UqX7TcdL6+opY9NQ0BRyMi0nmRKD7MrBT4I/Aj59y6Lrz088AaYE1FRUVKYhMRkXDoQa6AiOWLE8YMJC/HS+GvaMhdEYmQ0BcfZjYceA54GvhaF19+OzAJmFRaWprs0EREJCR6mCsgYvmiIDeb2WMGAvDyenW9EpHoCHXxYWZlwHzg7865zznnXFde75yrcs6tdc6tzcnRZO4iIr1RT3MFRDNftJ7voxt/sohIIIIeajfbzAqAPP/3Av9hZnYs8BLwoHPuy0HGKSIiwVGuaF/LfB87qxvYXHUo4GhERDon6Ms7VwO/bfV7vf/vWOCrwAjgZjNrPbziTc65+9MUnwQgLy+PmTNnAmBmAUcj0jEdr2mhXNGOqSMGUJSfw8HGGC+v38vYQX2DDinj6PyPDu2r8Ah6qN15zjlr57HZOXe9/3O/No9enUzE+1LIysoiKytLXxASejpeU0+5on052VnMGee1fuim82Do/I8O7avwCPU9HyIiInJk75ngFR+vbqgikdB9HyISfio+JHTi8Tj79u1j3759xOPxoMMR6ZCOVwlSy03n+w81s2pXTcDRZB6d/9GhfRUeKj4kdGKxGJs2bWLTpk1EYbIvyWw6XiVIE0r7UVqUD6jrVRB0/keH9lV4qPgQERGJKDM7POrVS5rvQ0QiQMWHiIhIhLV0vVqwsYqGZnUnEZFwU/EhIiISYe+bNBiAxliCVzeq9UNEwk3Fh4iISISVFhUwZUR/AF5YUxlwNCIiHVPxISIiEnFnTCoF4J9vVeCchtwVkfBS8SEiIhJxc/3iY+u+Q2zaWxdwNCIiR6biQ0REJOJmjCqmuE8uAM+r65WIhFhO0AGItJWXl8f06dMByM7ODjgakY7peJUwyM4yTp84mMeX7eS5NRXc8J6xQYeUEXT+R4f2VXio5UNCx8zIyckhJycHMws6HJEO6XiVsJjrj3q1YNM+DjVpErV00PkfHdpX4aHiQ0REpBc4vXwwZtAUS/DqBg25KyLhpOJDQieRSFBdXU11dTWJRCLocEQ6pONVwmJQv3ymjSwG4Lk1FQFHkxl0/keH9lV4qPiQ0Glubmb9+vWsX7+e5ubmoMMR6ZCOVwmTueVe16vn11RqyN000PkfHdpX4aHiQ0REpJc441hvyN3t++vZUFkbcDQiIu+m4kNERKSXmDZiACV98wANuSsi4aTiQ0REpJfIyjJO97te6b4PEQkjFR8iIiK9SMuQuws37aO2UUPuiki4qPgQERHpRU6fOJgsg+a445X1e4MOR0TkHVR8iIiI9CID++YxY1TLkLu670NEwkXFh4iISC9zxiRv1KsX1lRoyF0RCZWcoAMQaSsvL48pU6Yc/lkkzHS8ShjNnVTKj/+xlp3VDazdU8ukoUVBh9Qr6fyPDu2r8FDLh4SOmZGfn09+fj5mFnQ4Ih3S8SphNHl4fwb1ywc06lUq6fyPDu2r8FDxISIi0stkZdnhUa+ee0vFh4iER6DFh5ldbmbzzazGzN41HqCZnWtmK82s3sxWmNk5QcQp6ZVIJKitraW2tpZEIhF0OCId0vGaesoV3dNy38eiLfs5cKgp4Gh6J53/0aF9FR5Bt3zsB34B3Nz2CTMbB/wJ+G9ggP/vn82sLI3xSQCam5tZs2YNa9asobm5OehwRDqk4zUtlCu64fTyQeRlZxFPOJ5drdaPVND5Hx3aV+ERaPHhnHvKOfcgsLGdp68FFjvn7nPONTnn7gfe8Jd3ipmVmFm5mZXHYppoSUQkilKdK6B35ouiglxOm1ACwNOrdgccjYiIJ+iWj45MBxa3WfaGv7yzPg+sAdZUVOiqj4hIL5SMXAG9NF98YPJQAF5YW0l9UzzgaEREwl18FAHVbZYdAPp34T1uByYBk0pLS5MVl4iIhEcycgX00nzx/uOGYAYNzQleXKcJB0UkeEkpPsysv5l9Nhnv1cpBvP67rRUDNZ19A+dclXNurXNubU6OpjQREQlaCvJFj3MF9N58Mbgon9ljBgLw1Ep1vRKR4PWo+DCz08xsHrAL+H5SInrbMmBWm2Uz/eUiIhIhKcwXyhVH0dL16tnVFcTiGuVHRILV5eLDzI4xs1vMbCUwHyjHG4FkWDfeK9vMCoA8//cC/2HA74DZZvZxM8s1s48DJwD3dPVzREQk/ZKVL5Qreuac473io7q+mYWb9gUcjYhkuk63LZvZGcAngY/gDXt4L3CJc+6tHnz+1cBvW/1e7/871jm3wcwuBn4M3I03yslHnHObe/B5IiKSYinIF8oVPTC6pA/HDi3ird0HeWrlbk6dMCjokEQkg5lz7ugrma0DRgF/AeYBf3fORWrYjNmzZ7tFixYFHYZ0QiKRoKGhAYCCggKyssI8LoJkut50vJrZYufc7B6+h/JFCP3kH2u57dl1DBtQwCv/fiZeo5H0VG86/3s77avk6km+6OyWHw9sAp4HXo5aIpFoycrKok+fPvTp00dfDhJ6Ol7fRfkihFru+9hV3cDyHW0HB5Pu0vkfHdpX4dHZrT8W+CPwNWCnmf3BzM41XToREZF3Ur4IoeOGFTHqmEJAo16JSLA6VXw457Y4574JjAY+DvTBa1LfZmbfM7PyFMYoGSaRSFBfX099fT2JhEZmkXDT8fpOyhfhZGaHbzx/auWegKPpPXT+R4f2VXh0qd3JORd3zj3qnPsQUAY6TXYPAAAgAElEQVTcBVwFrDazl1MQn2Sg5uZmVq1axapVq2hubg46HJEO6Xhtn/JF+LR0vVpfUcuGytqAo+kddP5Hh/ZVeHS705tzbrtz7tt4SeXDgKZOFRGRd1G+CIcTxgykpG8eAE+r9UNEAtLjO26ccwnn3F+ccxclIyAREemdlC+ClZ1lnH38EED3fYhIcDo9zweAmT0HtDc2rwMagHXAPOfc0iTEJiIiEaV8EU7nTB7C71/fxtJtB9hd3cDQAQVBhyQiGaarLR9rgVl4NxLu9h+j/GVVwJnAQjN7XzKDFBGRyFG+CKFTxw+ib142AP9Yra5XIpJ+XS0+DgAPAuXOuSucc1cA5cADwA7n3DS8GWb/K7lhiohIxChfhFBBbjZzjy0F4Gl1vRKRAHSp2xVwPfBe59zhMcqcc87Mfga8hDeu+y+AF5IXYnI452hsbAw6DOmE1vtJ+0zCTsfrESlfhNT7y0v465u7eHVDFXv211LcJzfokCJL5390aF+FR1eLj3y80UrWtFleBuT5Px8CQjeZVH19PStWrAg6DOmidevWBR2CSKfpeH0H5YuQKo0lyMuCpoRj3jNLOGtcn6BD6hV0/keH9lWwutrt6o/AnWZ2pZlN8B9XAncAD/vrnAK8lcwgRUQkcpQvQqowN4tZw/IBeGlbQ8DRiEim6WrLx2eBH+FNFtXSTtsM3Al82f99CXBjUqJLosLCQqZMmRJ0GCIimUL5IsSuyNrDaw8tZ0VlE6VjJlJalB90SCKSIbpUfDjn6oHPmtlXgPH+4g3OubpW64SyrdrMyM/Xl6uISDooX4TbOVNG0O/R1dQ2xnhmTRXXnzY26JBEJEN0a5JB51ydc+5N/1F39FeIdJ5zjqamJpqamnCuvWkCRMJDx2vHlC/CqSA3m3P8CQefWLYz4GiiS+d/dGhfhUeXiw8zm2xm95jZ6/5jnpkdn4rgJDM1NTWxfPlyli9fTlNTU9DhiHRIx+uRKV+E2wXThwPwxtYDbNt3KOBooknnf3RoX4VHl4oPM/sgsBSYADzjPyYCS83svOSHJyIiUaR8EX6nTRh0eJjdv7y5K+BoRCRTdLXl47vA/zrnTnPOfc1/nAb8FE0UJSIib1O+CLm8nCzOmzIMUNcrEUmfrhYfx+GNVNLWnYCa0kVEpIXyRQRcMN0rPlbtqmF9RW3A0YhIJuhq8bEPL6G0daz/nIiICChfRMKcsSWHh9lV64eIpENXi48HgN+Y2fVmNtF/3AD8Brgv+eGJiEhEKV9EQHaW8cFpb3e90ihAIpJqXZ1k8Gt4BcsvgTx/WRPwc+AbSYxLRESiTfkiIi6YPpzfvryZjXvrWLmzhikjBgQdkoj0Yl1q+XDONTvnvgQcA0z3HwOdc192zjWnIkAREYke5YvomDmqmJEDCwF44k11vRKR1Dpqy4eZ/a0T6wDgnDs/CTFJhsvJyWH8+PGHfxYJMx2vb1O+iCYz44Lpw/nl8xv4y7JdfPUDx5KVZUGHFQk6/6ND+yo8OrP196Q8CpFWsrOzKS4uDjoMkU7R8foOyhcRdaFffOw4UM+Sbfs5YcwxQYcUCTr/o0P7KjyOWnw4565PRyDtMbOhwG3AmXixLgFucc4tCyomERFpn/JFdB07tIgJpf1YX1HL40t3qvgQkZTp6mhX6fYLvP7C5cAQYBHwF2tpt5deyTlHPB4nHo9r5BUJPR2voaF80QNmxgXThgPw1+W7iMUTAUcUDTr/o0P7KjzCXnxMAB52zu13zjUBdwEjgZLOvNjMSsys3MzKY7FYKuOUJGpqamLp0qUsXbqUpqamoMMR6ZCO19BQvuihC2d4xcfe2iZeWr834GiiQed/dGhfhUeniw8zG2dmnzWzK8ysb5vn+pvZ3ckPjx8Cl5jZYDMrAD4FvOSc6+y34ueBNcCaioqKFIQnIiJtKV9E09hBfZk52usT//Di7QFHIyK9VaeKDzM7GVgG3Ar8ClhlZjNarVIIXJv06OBlIBuoAGqBi4FPduH1twOTgEmlpaXJj05ERN5B+SLaLj1hFAD/WLmHA4d0dVhEkq+zLR/fAf4AlAJDgb8Az5rZzFQFZmZZwDPAWmAA0Af4L2C+mQ3pzHs456qcc2udc2s1rJqISFooX0TYh6YPoyA3i6Z4gseXac4PEUm+zhYfs4AfOM8h59xngZ8Az7S5opVMxwBjgdudczXOuSbn3J1+zKek6DNFRKRnlC8irH9BLudOHgrAI+p6JSIp0NniIxvIbb3AOfdd4Od4V5umJzku/H66a4HPmFlfM8sxsxuAIuDNZH+eiIgkhfJFxF062+t69eb2atbsPhhwNCLS23S2+HgLOLHtQufct/FGFPljMoNq5SK8q1lbgCrgs8ClzrmNKfo8ERHpGeWLiDtlXAkjigsBeHjRtoCjEZHeprPFx5+AS9p7wjn3VWAekPSx1J1zq51zH3LODXLODXDOneCceyzZnyMiIkmjfBFxWVnGJbNGAPDo0h00a84PEUmiThUfzrkfOOc+1MHzn3fOhX3OEImInJwcysrKKCsrI5Nv/JRo0PH6TsoXvcNH/VGv9tY28dxbmTn0cGfo/I8O7avw6HYCMLPPmFn/ZAYjApCdnU1JSQklJSVkZ2cHHY5Ih3S8Hp3yRfSMLunDnLHHAJrzoyM6/6ND+yo8enL16Xa8YRRFREQ6onwRQS03nj/3VgV7axsDjkZEeoueFB9J77MrIiK9kvJFBJ0/dSh987KJJRyPLtkRdDgi0kuo362ETmNjI4sXL2bx4sU0Nv7/9u47PIrr7Pv49+yuKkhCgCQQEghEbwIM2DTbce/GjnuvuDy2kzhxEqc5yZvHSZ5Ux4njjm3cYmzc4t6pNr13hIQESEJCqKGu8/6xC5ZljCXQama1v8917QWamV3dmnbvPefMGV1tE3fT/iqdVWykj7NH9wb8z/yw1jockfvo+A8d2lbuoeJDREREDulA16uNBRWs3VnucDQi0hmo+BAREZFDGt8vkf49uwAwe5me+SEiR0/Fh4iIiBySMYaLjkkD4PWVu6ipb3Q4IhEJdSo+RERE5BtdOK4PxkBZdT3vrStwOhwRCXEqPkREROQb9U6I4aQhyQA8+3muw9GISKg7muLjN0BxewUiIiKdlvJFiLtqUj8AluSUsmG3bjwXkSN3xMWHtfY31tq97RmMiIh0PsoXoe+EQUmkd48B1PohIkdH3a7EdXw+H+np6aSnp+Pz+ZwOR+SwtL9KOPB4DFcd62/9eHXFTipq6h2OyB10/IcObSv3UPEhruP1eklOTiY5ORmv1+t0OCKHpf1VwsXF49OJ9HnYX9fIq3riOaDjP5RoW7mHig8RERH5Vt27RHJO4Innsxbl6onnInJEVHyIiIhIq1x9nL/r1ZaiSr7Yrtt4RKTtVHyI69TW1rJy5UpWrlxJbW2t0+GIHJb2VwknY9K7MbJPPOBv/Qh3Ov5Dh7aVe6j4EFdqbGyksVFP0pXQoP1VwoUx5mDrx3vrCigsr3E4Iufp+A8d2lbuoOJDREREWu28rD7ERftoaLK8uDjP6XBEJMSo+BAREZFWi4n0cvEx6QA8vziX+sYmhyMSkVCi4kNERETa5Mrj+gJQWF7LRxsKHY5GREKJig8RERFpk8ykrkwd2BOAWXriuYi0gYoPERERabOrAjeeL9hawtaiSoejEZFQoeJDRERE2uyUYcn0TogGYOaC7Q5HIyKhQsWHuI7X6yU1NZXU1FS8Xq/T4YgclvZXCVc+r4drJ2cA8PKyfEoqw+/ZCTr+Q4e2lXu4vvgwxpxijPncGFNpjCk2xjzkdEwSXD6fj969e9O7d298Pp/T4YgclvZX91C+6HhXHNuXrlE+ahuawvLeDx3/oUPbyj1cXXwYY04EXgb+DPQA0oDHnYxJRETcR/nCGfHREVw2wT/s7jOLcqmp1wPcROTwXF18AL8HHrbWvmytrbXW1lhrl7f2zcaYHsaYwcaYwQ0NDUEMU0REHKZ84ZDrp/bH6zHsrarjleX5TocjIi7n2uLDGNMFmAj4jDHLA03onxpjxrfhY+4ENgGbioqKghKntL/a2lrWrFnDmjVrqK0Nvz7EElq0vzpP+cJZfbrFcM7o3gA8Pm87jU3W4Yg6jo7/0KFt5R6uLT6ARPzxXQ5cB6QC7wNvG2O6tfIzHgSGAEOSk5ODEaMESV1dHXV1dU6HIdIq2l8dp3zhsJunDQBge3EVH4bZQwd1/IcObSt3cHPxURH4d6a1drW1tg5/s3oEMLk1H2CtLbHWbrbWbtbNRSIinZbyhcNG9klgysAeADw2N9vhaETEzVxbfFhry4AcoGX7rT3ENBERCVPKF+5woPVjaW4py3JLHY5GRNzKtcVHwEPA9caY4cYYH3APUAssdDYsERFxGeULh50wOIkhKXEAPD5PrR8icmhuLz7+DDwJfAwUA2cCZwaucomIiBygfOEwYww3H+9v/Xh3XQE5xVUORyQibuTq4sP6/cpa28ta281a+x1r7Uqn4xIREXdRvnCH87JSSYmPwlp4Yv52p8MRERdydfEhIiIioSPS5+G6yf0BmL0sj71VGllIRL5KxYe4jtfrJSUlhZSUFLxer9PhiByW9leRr7ri2L50ifRSU9/EUwtznA4nqHT8hw5tK/dQ8SGu4/P5SEtLIy0tDQ15KW6n/VXkqxJiIrjyuH4AzFywnbLqeocjCh4d/6FD28o9VHyIiIhIu7p52gCiIzxU1DQwc4Hu/RCRL6n4EBERkXaVFBfFVcf6Wz+emN+5Wz9EpG1UfIjr1NXVsX79etavX09dnW5WFHfT/ipyaDNO+LL146kFOU6HExQ6/kOHtpV7qPgQ17HWUl1dTXV1Ndbq4cTibtpfRQ4tOS66WetHdqds/dDxHzq0rdxDxYeIiIgExYwTBhDl81DeiVs/RKRtVHyIiIhIUCTHRXPVcV+2fpTXdL7WDxFpGxUfIiIiEjS3qPVDRJpR8fEtiipq+Nmra6iua3Q6FBERkZCTHBfNlYF7Px6fp9YPkXCn4uMwqusaueyRz3n+ix3c+uwyahtUgIiIiLTVrc1aP55W64dIWFPxcRgxkV4unZAOwGeb9/C9F1bS0NjkcFQiIiKhJTk+miuO7QvA4/O3q/VDJIyp+PgWt5yQyV0nDwLg3XUF/Pjl1TQ1aYi2YPJ4PPTs2ZOePXvi8WgXFXfT/irSOredkEmUz0NZdT0z5+c4HU670PEfOrSt3MPndACh4AenDKKqtoEn5m9nzoqdxER6+d30kRhjnA6tU4qIiKBfv35OhyHSKtpfRVrnQOvHzAU5PDYvmyuP60vPrlFOh3VUdPyHDm0r91Dp1wrGGH5x9jAun+hvMn7uix3c//YGPaRGRESkDe74zkDionxU1jbw4EdbnA5HRByg4qOVjDH8bvpIpo9JBeCxedt5QCdOERGRVuvRNYpbT8wE/BfysvdUOhyRiHQ0FR9t4PUY/nxxFqePSAHg7x9u4ZHPtjkcVedTV1fHpk2b2LRpE3V1dU6HI3JY2l9F2uaGKf3pFR9NQ5PlT+9tcjqco6LjP3RoW7mHio828nk9/OPysRw/OAmA37+zkUfnqgBpT9ZaKisrqaysVNc2cT3tryJtExPp5e7TBgPwztoCluWWOhzRkdPxHzq0rdxDxccRiPJ5eeSqY5gysAcA97+9kcfmZjsclYiISGj47rg0hvaKA+D3uodSJKyo+DhCMZFeHr9mwsEC5H/f3qACREREpBW8HsNPzxwKwNLcUt5fX+hwRCLSUVR8HIVDFSCPz1MBIiIi8m1OGJx0MH/+8Z2N1OshviJhQcXHUTpQgEzO9J9Af/eWChAREZFvY4zh3jOHAZBdXMWLS/IcjkgkfDj5wGwVH+0gJtLLE9d+tQBRFywREZHDG9kn4eAQ9g98uJnK2gaHIxLp/GrqG7l25mJmLcpx5Per+GgnLQuQ/317Aw98uEU30YmIiBzGD08bQqTXQ3FlnYavFwmy+sYm7nxhBfO2FPPL19exJr+sw2MIieLDGOMxxiw0xlhjTJrT8XyTAwXItEE9Afjbh5v5w7sbVYC0kcfjITExkcTERDyekNhFJYxpf3WXUMkX8qX07rFcNyUDgEfmZrOjZL+zAbWBjv/QoW3l72r1o9mr+CAwwMMPTx3MqLSEDo/D1+G/8cj8ADiqs5G1ltra2nYK55t5gIcuG833Z6/ho417eOSzbCr21/HLs4bg8Zig//7Ook+fPgA0NTV1yHYTORraX10lZPKFfGnGlL68ujyfPZV1/Or1NTx8RRbGhEbO1PEfOsJ5W1lrue+/G3l95S4AbprSj5unpDuyHlxffBhjBgO3A98FVrTxvT2AHgCDBw9m7dq17R/gN5gxwkNtVTTz82p4fkk+u4qKuX18Al4VICIiQRGq+UL8rhgewwOL6/h0czFPfbCMCanRTock0ilYa3lmdQVvbPZflzltQAxn9K5h3bp1jsTj6nYnY4wHeBL4EbDvCD7iTmATsGnv3r3tGdq38nkMdx2bwMn9YwD4NLeGv31RRr2DowuIiHRWoZwvxG9a32hGJEUA8MSKcmoblC9F2sMrG6oOFh7H943m5nHxjrYsur3l43tAgbX2VWNMxhG8/0HgeYDU1NRNI0eObMfQWuefIyz3v7uZWV/ksSi/hsiYLjxwyWhiIr0dHkuoqK+vZ+fOnYC/iTQiIsLhiES+mfZX1wj5fCHwx5RKpv/7C/bsb2JeSSzfPznT6ZAOS8d/6AjXbfX0oh28sK4AgFOGJvHAJaPweZ1te3Bt8WGMGQj8EBh/pJ9hrS0BSgDGjx9PVFRUO0XXNr+dPor42Ej+9ck2PttSwg2zVvDkdRPoFhvpSDyhoLKyEoCIiAjHtptIa2l/dVZnyhfhbmR6FDdM7c+jc7N5YkEul0zsR/+eXZwO67B0/IeOcNtWLyzewf3vbgZg2qCe/OuqY4jyOX/x283drqYCScBaY0wxsDwwfbUx5nbnwmo7Ywz3nD6Un501FIDlO/Zx0cOL2LWv2uHIREQ6hU6TLwTuOnkQKfFR1DU28avX12rESJEj8NKSPO6dswaA8f0SeeRqdxQe4O7i4yUgExgTeJ0VmH4a8IxTQR2NGcdn8tdLsvB5DFuLKvnuvxeytajC6bBEREJdp8sX4axrlI9fnjMcgHlbinl3bYHDEYmElleW5fOTOasBGJPejZnXTyA20j2dnVxbfFhr91tr8w+8gANnnwJrbaWTsR2NC8el8di144mJ8LK7rIaLHl7EstxSp8MSEQlZnTVfhLOzR/Vm6kD/M7N++9/17K/Tk89FWuO1FTv50cursBZGpyXw9A0TiYt21/0tri0+WrLW5lhrTSCxhLTvDEnmuZuPpVtsBPv213Pl45/zycYip8MSEekUOlO+CFfGGH593ggivIbdZTU88NEWp0MScb03V+3i7pdWYi2MSI1n1g3HkhDjrsIDQqj46GzG9U3k5VsnkZoQTU19Ezc9s5Tnvsh1OiwRERFXGJjclZumDQDg8XnbWZNf5nBEIu71zprdfP8/K2myMKx3PM/eeCwJse4rPEDFh6MGJsfxyu2TGZISR2OT5eevruX3b2+gSc8CERER4a6TBtG/Zxcamyw/mr2K2oZGp0MScZ131xZw5wsraGyyDEmJ49kbJ5LYxb0jqqr4cFjvhBhm3zaJaYP8fVsfmZvNHS8sp6Y+fE+wHo+H+Ph44uPj8Xi0i4q7aX8VCZ6YSC//d9FojIFNhRX88+OtTof0FTr+Q0dn3VZvrtrF/zy/nIYmy6Dkrjx387H06OruYYRNuAxhN378eLt06VKnw/hG9Y1N/PK1tby4JA+AsX278dg14+np8h1IRDoPY8wya+0RPyujs3B7vghHv31zPU8u2I7XY3jt9imMSktwOiQRx726Ip8fvrSKJguDU7ry7E3HkhwX3SG/+2jyRecp/UJchNfD7y8cxY/PGALAih37uOChBWwt0kAtIiIS3u45fQgZPWIPdr+qa2hyOiQRR720JI+7A4XHsN7xvHDzcR1WeBwtFR8uYozh9hMH8uDlY4n0ecjbW82FDy1g7uY9TocmIiLimJhIL3+6OOtg96sHP9boVxK+nv08lx+/shprYVSfBF4Iga5Wzan4cKFzs1J5/qZjSYyNoLymgetmLuaJ+dvD5imv9fX1ZGdnk52dTX19vdPhiByW9leRjjEhozvXT+4PwEOfbnPF6Fc6/kNHZ9lWMxds5xevrQX8DxB89qZj6Rbr3pvLD0XFh0uNz+jOG3dMZUhKHE0W/t9/13PPy6vDYqSPpqYmSktLKS0tpalJTevibtpfRTqO27pf6fgPHZ1hWz382TZ+8+Z6ACZkJDLrxomufI7Ht1Hx4WLp3WN55fbJnDY8BYCXl+Vz+aOfU1RR43BkIiIiHU/dryQcWWv547sb+cM7GwE4bkB3nrrefU8uby0VHy7XNcrHw1cdw10nDwJg+Y59nP/PBa5obhYREelozbtf/euTrSzevtfhiESCp7HJ8ovX1vLvT7cB8J0hScy8biJdonwOR3bkVHyEAI/HcPepg3noynHERHjZXVbDRQ8vZPbSPKdDExER6XD3nD6EQcldabLwvRdXUFpV53RIIu2uvrGJ7/9nJc99sQOA87JSefSa8cREeh2O7Oio+AghZ43qzcu3TaJPtxhqG5q45+XV3DtnTVg/kFBERMJPTKSXf14xjiifh91lNfxo9qqwGZRFwkN1XSMznlnKm6t2AXDlsX3526VjiPCG/lf30P8LwsyI1ATevHPqwSeiv7B4B5c8soi8vfsdjkxERKTjDOkVx6/PGwHARxuLeHJBjrMBibST8pp6rn1yMZ9s8j9q4fYTM/nd9JF4PcbhyNqHio8Q1L1LJE9dP/HgfSCr88s495/z+XRTkcORiYiIdJzLJqRzzujeAPzhnQ2szt/ncEQiR6ewvIbLHvmcxTn+e5nuPXMoPz5jKMZ0jsIDVHyELG/gPpCZ100gISaCffvruf6pJfz9w800NoV207Mxhi5dutClS5dOdbBJ56T9VcQ5xhh+f+Eo+naPpb7RcucLK6io6bhnOOj4Dx2hsK22FlVw4UMLWb+7HI+BP1w4iltOyHQ6rHZnwqWP5Pjx4+3SpUudDiMo8vbu57bnlrF2ZzkAkzN78PdLx5AcH+1wZCISSowxy6y1452Ow2mdOV90Vqvz9/Hdfy+kvtFyblYq/7hsjGu/YIocytKcvdz49FLKquuJjvDw4OXjODXwqAU3Opp8oZaPTiC9eywv3zqZyyemA7BwWwlnPjCPT9QNS0REwsDotG789MxhALy5ahf/WaLRICV0vLu2gCsf/4Ky6noSYyN47qbjXF14HC0VH51EdISX3184mn9cPpauUT5Kquq4fuYS7n97g+NPgBUREQm2G6ZkcMqwZADue2Od7v+QkPDMohxue24ZtQ1NpCXG8PJtkzmmX6LTYQWVio9O5rysVN66aypZaQkAPDo3m4sfXsiOktAZDau+vp7c3Fxyc3Opr++4vrsiR0L7q4g7GGP400VZB4ejn/HMMorKa4L6O3X8hw63baumJssf3tnIr15fh7UwIjWeObdPJjOpq9OhBZ2Kj06oX48uzL51MjOOHwDAqvwyzvrHPOYszw+JcdCbmpooLi6muLiYpia12oi7aX8VcY/ELpE8es0xxER4KSivYcasZUF9FpaO/9Dhpm21v66B259bzsOf+Z9aPm1QT/5zyySS48LjXl0VH51UpM/Dz84axszrJ9C9SySVtQ3c/dIqbn9uOXv1JFgREemkRqQm8NdLsgBYmbePn726JiQuvEl4KCir4ZJHFvHuugIALh2fzhPXTqBrlM/hyDqOio9O7jtDknn3e9M4cUgSAO+sLeD0v8/lk426GV1ERDqnM0f15vun+J+FNWf5Th6bl+1wRCKwdmcZ5/9rPmt3lmMM/PysYfzhu6OI9IXX1/Hw+mvDVHJ8NDOvm8Dvpo8kJsLLnoparn9qCT97dQ1VtQ1OhyciItLu7jppEGeN6gXA79/ZqItu4qj31hVw8cOLKCyvJSbCyyNXHcPNxw8IyyGhVXyECWMMVx3Xj7e/N42xfbsB8PwXOzjrH/NYEniKpoiISGfh8Rj+fHEWw3vHYy3c9cIKthZVOB2WhBlrLf/+dBu3PruM6vpGeidEM/vWSZw2opfToTlGxUeY6d+zC7NvmcSPThuMz2PILdnPJY8s4r7X11KpVhAREelEYiN9PHbteHp2jaSitoGbnl5KSWWt02FJmNhf18AdL6zgj+9uxFoYnZbA6/8zhZF9EpwOzVGuLj6MMX80xqwzxpQbY3YZYx4zxnR3Oq5Q5/N6uOOkQbz2P1MY2isOa+HpRbmc/re5zN28x+nwRETaRLlCDqdPtxgevuoYIryGnJL9XP/UEl1sk6DbUbKfCx9ayFurdwNwzuje/GfGJJLjw2NEq8NxdfEBNAJXAT2ALCANeMrJgDqTkX0SeOOOqdx96mAivIad+6q55snF/Gj2Ksr2OzcGtjGG6OhooqOjw7IvpIQW7a+uoFwhhzU+ozt/uWQMxsDq/DJmPLOU2oajH4JXx3/o6Mht9dnmPZz7z/lsLKjAY+DeM4fy4OVjiYn0BvX3hgoTSsPPGWPOAF6y1sa3cvke+JMRWVlZm1auXBnM8ELa5sIK7nl5Navy/E+ETYqL4tfnjuCsUb10QhUJE8aYZdba8U7HcbTamisC71G+CAOzPs/ll6+tBeDMkb345xXj8HqU46R9WGt5+LNs/vTeRposdIuN4MHLxzJtUJLTobW7o8kXbm/5aOlkYFUblr8T2ARsKirSKBeHMzgljjm3TeYXZw8jOsLDnopa/uf55Vw7cwk5xVVOhycirVDX0MRn6joJbc8VoHwRFq4+rh93nzoY8A89/4vX9AwQaR8VNfXc8bz//o4mC8N6x/PmHVM7ZbRcYVUAACAASURBVOFxtEKm+DDGfBe4FfheG972IDAEGJKcnByUuDoTr8dw07QBvPu945k6sCcAczfv4bS/z+WBD7cE9SmxInLk6hqaeO6LXE780ydc++RiNhWE74g+R5grQPkibNx50kCum5wBwAuL8/jTe5ucDUhC3rpdZZz74HzeWuO/v+O8rFTm3DaZ9O6xDkfmTiFRfBhjLgYeA86z1i5v7fustSXW2s3W2s0+X/g8OfJoZfTswqwbJ/Lg5WNJjouirqGJv324mTMfmMe8LcG/qtrQ0EB+fj75+fk0NOimQHE3J/fX5kXHz19dy66yGjwGPs8u6dA43OJIcwUoX4QTYwy/Omc408ekAvDQp9t4bO6RPYRQ+Sp0BGNbWWt57otcLnhoITkl+4nwGn55znAeuGyM7u84DNefYY0x1wN/Ac611i5wOp5wYYzh3KxUThiSxF/f38wzi3LYXlzF1U8s5uxRvfnpmUODVtE3NjZSWFgIQFJSEvoiIG7mxP5a19DE7GV5/OvjrewqqwHAY/xX2+44aRADk7sGPQa3Ua6QtvB4DH+6OIvymgY+3ljE/769gagID9dMymjT5yhfhY723lYVNfX87NW1vLlqF+AfVe1fV45jTHq3o461s3P1UWKMuQu4DzjdWrvE6XjCUXx0BL8+bwQXHZPGz19by6q8fby1ZjcfbChkxrQB3HZiJl2iXL0biXQa1XWNvLhkB4/OzWa3io6DlCvkSER4PfzrinFc8+QXLMkp5Vevr6OmvpEZx2c6HZq43LpdZdzx/Aq2B+6JPXV4Cn++KIuE2AiHIwsNbv/W+ADQAHzSfMQla214ZlgHjeyTwKu3TWb2Mn//2OLKOv75yVZeWprHj88YyoVj++DRiCEiQVFeU8+sRbk8OX87JVV1gIqOFpQr5IjERHp56vqJ3Pj0Ej7P3sv9b2+kuq6Ju04eqJEe5WuamiwzF+bwx3c3UtfQhM9j+OmZQ7lxan/tL23g6uLDWqst6SIej+HSCX05a1Rv/vnJVmbOz6GoopYfzV7FrEU5/PKc4YzP0HO9RNpLSWUtMxfk8PSiHCpq/H2UfR7DBWP7cOuJmWQm6bs1KFfI0ekS5eOp6ydyy6xlfLZ5D3/7cDM1DY38+PQh+kIpBxWU1fCj2auYv7UY8Hez+ucVYxnbN9HhyEKPq4sPcae46AjuPXMYV0zsy/1vb+C9dYWsyi/joocXccqwFH58xhAGp8Q5HaZIyMopruKJ+duZvSyPmvomAKJ8Hi6bkM6MEzLp0y3G4QhFOpfoCC+PXnMMdzy/gg/WF/LvT7dRXdfIfecOVwEivL1mN/fOWUNZtf8BzNPHpPKb80eSEKNuVkdCxYccsX49uvDI1eNZuLWY//fWBjbsLufDDYV8vLGQC8el8YNTB+tLkkgbLMst5bG52by3voADjx7oGuXj6kn9uGFKf5LiopwNUKQTi/J5eejKcfzgPyv57+rdPLUwh9qGRn43fZQeRBimKmrque+NdcxZvhOAuGgfv5s+kvPH9HE4stCm4kOO2uSBPXnrzqm8sWoXf/lgE3l7q3l5WT5vrNrFNcf14/bvDKR7l0inwxRxpcYmy4cbCnlsbjZLc0sPTk+Jj+L6Kf25fGJfXV0T6SARXg8PXDaW6AgvLy/L54XFeZRU1vH3y8YQG6mvTOFk4dZifvzKavJLqwGYNKAHf7kki1RdVD1qOpKkXXg8hulj+3DWqN48/0UuD368lZKqOh6fv50XFu/guikZ3DR1AImtLEIiIvRlS0LHkeyvZdX1zF6ax9OLcsjbW31w+pCUOG4+fgDnZaUS6QuJRzGJdCpej+H/vjua2EgvzyzK5f31hVzyyCIev2YCvRKiv7a88lXoaM22Kq+p5/63NvDikjz/e7yGe04fwk1TB2hgnXZi7IG2/U5u/PjxdunSpU6HETYqaxt4fF42j83NpqrO/2T0LpFerp2cwU3TBqglRMLW1qIKnlqYwyvLdlJd33hw+uTMHsw4fgAnDE5yrI+5MWaZtXa8I7/cRZQvBPwPkHti/nb+9+0NWOtvjXzi2gmM7JPgdGgSJB+uL+Tnr62hsLwWgFF9Evi/i0YzrHe8w5G5z9HkCxUfElR7q+p4bF42zyzMOViExAaKkJtVhEiYaGhs4uONRcz6PJd5W4oPTo/yebhwXB+unZzB0F7OJzcVH37KF9LcB+sL+d6LK9hf10hMhJe/XzaG00f0cjosaUcllbX85s31vBF4YGCkz8Pdpw7mpqn98XnVAn0oKj5aQcnEWXur6nh8XjZPNytCYiK8XDohnZum9SctMThPSxdxUmF5DS8uzuPFJTsOPhQQIDUhmqsnZXDZhPRWd0XsCCo+/JQvpKV1u8q48amlFJTXYAz89IyhzDh+gEbCCnFNTZY5K3Zy/9sb2Bt4htKEjET++N3RDNBQ5oel4qMVlEzcobSqjsfnZ/PUgi+LEK/HcO7o3sw4PpPhqfE0NDSwZ88eAJKSkvD5dGuSuFfL/dXj8bJwWwnPfp7LBxsKaWz68hw7aUAPrp7Uj9OGp7jyapqKDz/lCzmUwvIabnp6KWt2lgFwblYqvz13GDUV/oEilK/creW5enPRfn71+tqDA310ifTykzOHctWx/XRvRyuo+GiFY445xi5cuNDpMCSgrLqe5xfnM+uLvINPbAaYOrAH10zoTff6PRhjGDRoEFFRGl5U3Ku2tpYtW7awZ38jq8tjeH1NETv3fdnKER/t44Ixvbl0fBqZSV0cjPTbRUdHq/hA+UK+WXVdIz95dR3vrS8CID0xmjvGdWFg9wjlK5c7cK6uqmvi3Z0+Xly2iwPXhk4dlsS9Z+jxAG1xNPkibEr06upq1q5d63QY0sy0HjDx9EQ+y6nm9c1VFFQ2Mn9rCfO3ltA33seZA2OpadhEtEb8EZeqa7Qs2VXDR9urWV1YR/NLOYO6R3BaZgxT0mOI8tZTXbidtYWOhSptoHwhh3PzcEPviDieXV1BXmkNP/+4hqtHx3G23axuWC7WZC2f5dYwa3UFZbX+h7f26urlxjHxjOvtpTR/G6X5DgcZJsKm+BB3ivIaTsuM5eQBMSzeWctrG6vYWlrPjvIGHllezqw1FZyUEcOZA2Pp1VW7qzjPWsumknrm7qhmwY4aKuu/LDniIg3H94vhpIwYMrpp+E2RzsgYw7mDuzCsZyR//XwfhVWNzFxVweqiOu6YkEB8lC6Yuc3aolpmra5ka6n/CeWRXvjusK6cN7gLkV4VjB1N3a7EVay1LMku5pFPNrEov4bGwO5pDJwwqCeXju/D8QN7uLK/vHRu24ureHN1AW+sLiCv9MvnchhgTK9ITu4fy2XTRhDXJXSb7dXtyk/5QlqrpLyKe15awYI8f1fLlPgo/jB9OJMzezgcmQBsKarkLx9s5ZPNX44yOCE1it9OH82AFA2ZfDTU7aoVjDHqixkiJmYmEVtdyHXVjays7MqLS3dRXFnLp5uL+XRzMclxUVx0TBqXTkinXw9396GX0FZQVsPba3bz+qpdrMrb95V5A5O7csHYPpw9IomSvK0AxHWJ0XmmE1C+kNbqEQ8/ODaB0SmRzFxVSWF5Ldc/s4ILx/XhF2cP13DyDiksr+FvH2zmpaV5B+/rGJEax8WDfIxKjmJASoKOcQeFTfEhoScxxsudEzK565ShvLN2N899sYPF2/dSVFHLQ59u46FPtzFpQA8um5jOacN7ERPpdTpk6QQKy/0Fx1urdx8cBeWAnl2jOH9MKheM7cOI1HiMMdTW1lKS51CwIuI4Ywyn9I/lnGOH87PXN7Aqv4w5y3fyycYifnnOcC4Y20f3gnSQ4spaHgsM619T77+vIy0xhh+fMZRTB3dn/fp1DkcooOJDQkCkz8P5Y/pw/pg+ZO+p5KWl+by8LJ/iyloWZZewKLuELpFeTh/Zi+lj+jA5U92ypG3y9u7ng/WFvLu2gCW5e2neGzUuysepI1K0b4nIYQ1M7sqc26fw9MIc/vz+Jkr313P3S6uYs3wn/3vBSLXUB1FReQ2Pzs3m2S9yDxYdCTER3HnSQK6e1I8on5fa2lqHo5QDVHyIK3k8h/6CNyCpKz89cyg/PG0wH28s4j9L8vh0UxFVdY3MWb6TOct30rNrFOdm9Wb6mD6MTkvQFSf5Gmst63eX8/66Qj5YX8j63eVfmd81ysdpw1M4e3Rvpg7qSZTv8K1q37S/ikjn1/z493oMN0ztz+kje/Gr19by0cYi5m8t5rS/zeW2EzO5adoAukbpq1d72V1WzSOfZfPC4h3UNviLjrgoH9dNyeCmqQNIiP3qwB86V7tD2NxwrodGdV5FFTX8d9VuXl+5k1X5ZV+Zl5YYwxkjenHGyF6M65uoBweFsZr6RhZll/DZpj18sL6QnfuqvzI/PtrHSUOTOXt0KtMG9SQ6Ivy68ekhg37KF9IerLW8vaaA+95YR3Gl/6p7jy6R3HnSQK44th+RGkb+iG0qqGDmgu3MWb6TukZ/0REf7ePGqQO4bkoGCTEabTDY9JDBVlAyCQ/Zeyp5feUuXl+5k5yS/V+ZlxQXxekjUjhjRG+OHdCdCHWf6fRyiqv4ZFMRn27aw+fZJQevjB3QOyGa04ancNqIXkzsr31CxYef8oW0p7Lqeh74cAvPfp578ItyWmIMPzxtMOdn9dFFsVZqarJ8vLGIJxdsZ+G2koPTu8VGcPO0AVwzqR9x0So6OoqKj1ZQMgkv1lrW7izn3XW7eWdtAdl7qr4yPy7Kx9RBPfnOkGROGJJESny0Q5FKeyqprGXhtpLAq5jcFgUowKg+CXxnSBKnjeh18KZx8VPx4ad8IcGQX7qfv32whTkr8g/eVza0VxzfP2Uwpw5Pwasi5JAqauqZvTSfpxflfOWcnpYYw3WTM7hsYl91ZXOAio9WUDIJHY2NjZSU+K9q9OjRA6/36Lu/bC2q4J01Bby7roB1u8q/Nn9473hOHJLEtEFJjO3bLSy73ISiffvrWJZbysJtJSzYWszGgoqvLZMQE8G0QT05cUgyJwxOIimufYdXDMb+6hQVH37KF9JaR3L8byqo4E/vbeLDDYUHp/XtHst1kzO4eHyart7jb+VYlF3Cy8vyeWft7oM3kQMc2787N0ztzynD2lawdaZztRuo+GgFJZPQUVtby9q1awEYOXJku4/FnV+6n0837eHTTUUs2FpCdX3jV+ZH+Twc0y+RSQN6MCmzB6PTuqlvrkvkl+5naU4pS3L2sjSnlE2FXy82IryGsX0TmZzZg2mDepKV1i2oI1QFe3/tSCo+/JQvpLWO5vhflruXv36wmQVbv+xC1DXKxyXj07l+Sgbp3WPbPV63yy2p4pVl+byyfOdX7suL9Hk4PyuV66ZkMCL1yB4O2JnO1W5wNPlC7VQSdtISY7nquH5cdVw/auobWZKzl0827uHTzUVk76mitqHpYNcdPoCYCC/H9EtkXN9ujO2XyLj0xK+NoCHtr6KmnjU7y1iVV8aqvH2syt/H7rKary3nMTA8NZ4pmT2ZPLAnEzISiY3UqU1E3O2Yft157qbj2LC7nJkLtvPayl1U1jbw5ILtPLVwOycOSeb8MamcOjylU5/TckuqeG9dAe+tK2RZi2crjeqTwEXHpHFeViqJemBjp9F592aRVoiO8DJtkL+71a8YTkFZDZ9n++8XWJRdQt7eaqrrG5m/tZj5W4sPvi8zqQvj+iaSld6NEanxDO0Vr4ccHoXSqjo27C5nfeC1Or+MbXsqOVTDbJTPw5j0bkzs353xGd0Z27cb8eqmICIhaljveP7voix+fMZQnvt8B7M+z6W4spaPNxbx8cYiYiK8nDI8hfOyUjlhcFLIt8Rba1m3q5z31xfy/rqCr3WX7dk1kulj+nDR+DSG9op3KEoJJhUfIs30Sohm+tg+TB/bB/B381m0rYRluaUs31HKliL/F+Jte6rYtqeK2cvyAf/V98ykroxIjWdEagLDesczMLkrKfFRuqG5maraBrL3VLFtTyVbiirYsLuCDbvLD9micUDf7rFkpXcjKy2BsX0TGdUnIeSTr4hISz27RvG9UwZx64kDeHvNbl5bsYv5W4uprm/kzVW7eHPVLuKjfZw6vBdTB/VgSmZPkkNksJS8vftZlF3C59v8DwZuec7vFhvBKcNSOGNEL04YkhT2Iw92dio+RA4jLTGWi8fHcvH4dADKa+pZuWMfy3eUsnzHPtbuLGNvVR1NFrYUVbKlqJLXVu46+P6uUT4yk7qQmdyVgcldGdCzK327x5LWPabTXq2vqm0gr3Q/eXur2bF3PztKqgLFWuVhiwyAPt1iGNY7jlF9upGVnsDotG50V1O7iISRKJ+XC8amccHYNEoqa3l7zW7eWLWLJTmllNc08MryfF5Z7r/wNSi5K1MG9mRyZg+O7d/DFV2C6xqa2FJUwbpd5SzevpdF20q+9lwl8J/vTx2ewukjejEhIzGo9+aJu7i6+DDGeIE/ANcB0cD7wC3W2uLDvU8kWOKjIzh+cBLHD04C/M3HheW1rNtVxrpd5Qf/zS/1n2graxtYlV/2tYcfgn8UpvTuMaQnxpKWGENKfDRJcVEkx0WTHB9FSny0q4YPbGqylNfUU1xZR1F5DQXlNRSW11JYXkNRRQ279tWQt3c/JVV13/pZB4qyIb3iGNY73v/qFe+KxCmhSflCOqMeXaO4elIGV0/KYOe+at5avYvPNu9hSU5p4Eu+/6LXUwtzAP8X+qG94hjaO46hveIZ1juOjB5dgvLFvrahkd37asgvrWZjQaDb7K5ytu2ppL7x631mYyO9TMjoznED/IOBaKjz8OWebzaH9lPgfOBYoAR4EpgFnOlkUCIHGGPolRBNr4RoTh6WcnB6VW0D2/ZUsm1PJVuLvnzt2Lv/4Em5rLqesp31rN359aF/D4iJ8JIYG0F8TAQJMRF0i/X/mxATQWykj5hILzERXqIjPERH+P/v8xqMMXiNwWMMHgMej6GpyVLfZGlobKK+0dLQ1ERDo6WmvpGqukb21zZQVddIVW0DVXUNlFc3sG9/HXv317Fvfz379vtbeFrL5zH0SfQXVwdafzKT/C1AyXHqjibtTvlCOrU+3WKYcXwmM47PpKa+kWW5pSzYWsyCbSWsyd9Hk4Wd+6rZua+ajzYWHXyf12NI6hpFSnwUyfHRJMf5L2717BpFlM9DhM9DpNcQ4fUQ4fXg8xiq6xuprG1g/4GcUNtIRU09u8tr2Flaza591RRV1B423rgoH2P6duO4AT04bkAPRqclqDuVAO4vPmYAv7XWZgMYY34MbDXG9LPW5n7bm40xPYAeAFlZWUENVKS5LlE+Rqd1Y3Rat69Mb2yyFJb7WwjySqvJD3RP2rlvP0UVtewpr6WituHg8tX1jVSXNbLrW7orOaFbbAQpcdGkJESTEhdFr4Ro0hNjSe8eS3r3GHrFR6sZXTqS8oWEjegIL1MG9mTKwJ6A/2LWup1lbCyoYGNBORsLKthUUEFtQxONTZaCQGs1fL0Vvj306RbD8NR4hveOP/hvWmKMLjLJIbn2OR/GmG5AKTDWWruy2fQy4Gpr7Rut+IxfA/cFftwPbDiCULxAClAINH7Lsm4QavFC6MWseINL8QbX4eLtZ61N6viQjo7yxRFTvMGleIMr1OKF0Is5KPnCzcVHOrADGGCt3d5sei7wc2vts634jINXsoASa23J4Zb/hs8YDGwChlhrN7f1/R0t1OKF0ItZ8QaX4g2uUIu3NZQvjoziDS7FG1yhFi+EXszBitfN3a4ODPzc8lGW3YBv7iTfTCB5tDmBiIhISFG+EBEJEa7tkG2t3Yf/Sta4A9OMMQOAeGC1U3GJiIi7KF+IiIQO1xYfAY8CPzHG9DfGxAN/BN6z1uZ0YAwlwG8InStioRYvhF7Mije4FG9whVq8raV80XaKN7gUb3CFWrwQejEHJV7X3vMBB8dt/yP+cdujgA+AGRq3XUREmlO+EBEJDa4uPkREREREpPNwe7crERERERHpJFR8iIiIiIhIh1DxISIiIiIiHULFh4iIiIiIdAgVHyIiIiIi0iFUfIiIiIiISIdQ8SEiIiIiIh0iLIsPY4zXGPMnY8weY0yFMeYVY0zPwyx/hjFmnTGm2hiz1hhzWov5A40xHxpjqowx+caYHzoVrzHmLGPMx8aYYmNMqTFmnjFmWotlrDFmvzGmstkrwaF4TwzE0zyWhS2WcdP6/VmLWCsD8f+j2TI5xpiaFsuMasd4Lwts13JjTEMrlh9vjFkc2ObbjDFXtZifbIyZE/jb9xhj/miMabdzQ1viNcYcZ4x5yxhTaIwpM8YsM8ZMb7GMa9avMSYjsP2rmsWS32IZN63fKw+x/zYaY95otsynxpjaFsuc017xhpq2nB8CyytfBC9e5Yu2xapcEcRccQQxK18cYK0Nuxfwc2AzMABIAF4B3vmGZQcA+4GrgEjgSqAKyAjM9wIbgAeBWGAcUARc6lC8VwIXAN0AH3AbUAmkN1vGAlNdsn5PBBoO81muWr+HeO9goAmY2GxaDnBVENfv6cDlwA2HW3eBZROAPcBP8D/1+dTA/jCp2TIfAHMCyw4IrIufOBTvWcA1QE/8F0emA9XABJeu34zA8ZR2mGVcs36/Yf+oAi5pNu1T4BfBWr+h9mrj+Uz5Irjxnni4fdxt6/cQ7+3QfNHGc5lyRfBjzkD5wv++YG4Ut76AXODGZj9nBnaIfodY9jfAvBbT5gH3Bf7/HfzJpmuz+f8P+MSJeL/h/QXAhc1+DnYyacv6PfFwB4Db1y/wZ2BZi2lBP+G1Zt0Flrk+8PeZZtNmATMD/+8f+Fszm82/EdjuRLzf8L7Pgbtdun4Pm0zcvn6BOwLnh4hm044omXTWl/KF8kV7rV+n8oVyRfDWbRvXsfJF4BV23a6MMd2AvsCyA9OstduAciDrEG/Jar5swPJmy2YBm621ld8wv6Pjbfn+UfivDKxpMWt2oKn9C2PMhe0R61HE6zXG5BljCgLNqM2Xc+36NcZEAdcBjxxi9l+NMXuNMSuNMbe0R6xHKAtYYQNniYCW+29Z4G9uPj/DGBPfQTF+I2NML2AEsKrFLLes3wO+CDSRf2qMObHZdFevX+AW4ElrbX2L6d8PrN91xph7jTERTgTnNOWLg5QvjjzeA+91e75Qrug4YZ8vwq74AOIC/5a1mL4PONTGjfuWZb9t/tFqa7wHGWOS8TcJ/9lau6XZrFPwV9hpwF+B54wxZ7RPuG2OdyMwJhDPUGA18LExJrXZ57ly/QIX4e9a8XyL6dfiby5NAe4B7nfwpHek+y+03zo+IsaYLvj337estR81m+Wm9VsMTMK//2YQ6IJhjBkdmO/m9TsFGA481mLWvcAgIAn/VbebgN92bHSuoXyhfNFcZ84XyhXBp3wREI7FR0Xg35Y3zHXDf/XiUMsfbtlvm3+02hovAIGT8SfA+/h3joOstR9Za2sCr/8Az+Lv+9vh8VprC6y1q6y1Ddbafdbae4G9wJnNPs916zfgFuC5FlfZsNZ+Zq2ttNbWW2s/wJ+wrzrkJwTfke6/B+Y5whgTB7yDv7/2Nc3nuWn9BuL43FpbZ62tstY+CMwHLg4s4sr1G3AL8L61dnvzidbaRdbaUmtto7X2c+BXOLf/Ok35QvniiONtwe35QrkiyJQvvhR2xYe1dh+wA/+NaAAYYwbgrypXH+Itq5ovGzCWL5v2VgGDA5X3oeZ3dLwYYzLw9zN+x1p7R4tm1ENpAoxT8X5LPK5bv4FlhgPTgIdb8Wvabf0egVX4rxQ213L/TQj8zc3n51hrW16B6RDGmB7AR8Au4GJrbd23vMXJ9XsoLfdfV61fAGNMd/wJz+37r6OULw5J+aJz5gvlCmeEZ75o600ineGFf7SKTfibvuKB2cC737BsJv4b2C4HIgL/Hmr0kgeAGPwHbyFwmUPxDgXygd99w/yRwET8zb8R+EeI2A+c51C8JwED8RfCXYFf429mTHfj+m32ngeARYeY3g//TY/RgdhPwH9F5s52jNcb+PzTgIbA/6NpdqNgs2W74R/B5J7ANj+ZQ49g8nLgb+8fWBc/dSjeXsBa4GnAGwLr97jAMeULLDMDqAGOceP6bfaeHwB5LddxYH85J3AsGvyJbxPwl/aKN9RebTk/oHwR7HiVL9oWq3JFEHPFEcSsfHHgve25EULlFVj5f8bf/64C/7BmPQPzrgQqWyx/BrAO/zBu64DTWswfiL/63o+/Av+RU/ECM/GPllDZ4nVlYP53An9DFVAKLKUdT8xHEO8P8I+wURU4MbxLs6Hy3LZ+A9NiAuvu2kN81nBgReBzyvGfHO9o53ivC2zjlq8M/FfXKoG+zZafACwO7L/ZtBj9A0gO/M0VgXXwf4DHiXiB+wLzqlrsvz9z4/rF/+VyayDeEvxXkE916/pt9p4NBEZgajE9Cf+IMWWBeDcHtklke67jUHodwflB+SJ48SpftC3WNp0bUK4IdszKF4GXCXyAiIiIiIhIUIXdPR8iIiIiIuIMFR8iIiIiItIhVHyIiIiIiEiHUPEhIiIiIiIdQsWHiIiIiIh0CBUfIiIiIiLSIVR8iIiIiIhIh1DxISIiItKBjDFjjDHWGJPSbJoxxuwzxlzqZGwiwabiQ6SdKamIiMi3GAvsstYWNpuWCSQAy50JSaRjqPgQaX9KKiIicjhj+Xo+GAdUAFsBjDFpumAlnZGKD5H2d9ikYozJNMasaT7TGBNljNlujBlhjLnFGPPvFvPXGmOGBTdsERHpIGOBFS2mjQNWWmtt4OeTA9NEOhUVHyLt79uSynYgzRjT/PibAcy11q4DRtGseDHGRAMZwOZgBi0iIsFnjDHAaL5+kWoCgdxhjJkK/BW4yBiz0hgzwBgz1BjzceDnD40xPQPLvmCM+Y8xZrExJtcYc3ZH/j0ibaXiQ6QdtSapWGubgB34CwqMMTHAp01ewQAAAoBJREFUD4H7Asu2fP8oYLO1tjFogYuISEfJBOLx5wEAjDG9gKkEzv3W2vnAEuB8a+0YYCfwCnB34OcPgB8E3p4FZFtrJwJX8mUuEXElFR8i7etbk0rABmBo4P//A7xprc0J/DwCmGOMyTHG5ADvAKuDG7aIiHSQsYF/7zTGDDbGnAC8BkQCHmNMZGD+EGBj4P/TgfnW2pWBn9cDyYGW8STgN82mJwb7DxA5Gio+RNpXa5PKBmCIMaYrcAfwOwBjTDqwx1rbz1qbYa3NAF4EvnKPiIiIhKyxwDygC/4LS48CfwF2428Frw90qSqz1jYE3jOcr+aBUfgLjZHAFmttTWD6OGBV0P8CkaOg4kOkfX1rUgksd6Dl43vAc81GxhoFrGvxmcNRy4eISGcxBlhurb3EWhttrR1irZ1trU211o4M3BuYAexq9p6d+HMBxpgBwNXAM/i7XPU1xkQbY7rgbwH5W0f+MSJt5XM6AJFO5kBS+X6L6bNb/LwBuBc4BTim2fTR+K9mNTcCtXyIiHQWY4GXvmWZjUBPY8xa/AOSzALOCoyUWA3cYK0tMcZkAXOAL4AI4H5r7YLghS5y9FR8iLSv1iQV8I9cNQr4ubV2X7Ppo4D/HvjBGNMdMNbagnaNUkREOlzgHsBewMrDLWetrQQmtpg8/RCLZgEzrLV3tU+EIsFnvhxOWkSORiCp7AbGNrspUEREJCiMMflA38AoiiIhQcWHiIiIiIh0CN1wLiIiIiIiHULFh4iIiIiIdAgVHyIiIiIi0iFUfIiIiIiISIdQ8SEiIiIiIh1CxYeIiIiIiHQIFR8iIiIiItIhVHyIiIiIiEiHUPEhIiIiIiId4v8DjUedMf1OX0sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 936x864 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xrange = np.linspace(0., 1.7, 50)\n",
    "with plt.style.context('seaborn-talk'):\n",
    "    fig, axs = plt.subplots(int(np.ceil(len(eff_CP4.columns) / 2.)), 2, figsize=(13, 12))\n",
    "    for ipoi, ax in enumerate(tqdm(axs.flat)):\n",
    "        yscan = profile_likelihood_ratio_gridscan(nll, ipoi, xrange, nll_grad, nll_hessian)\n",
    "        ax.plot(xrange, yscan)\n",
    "        ax.set_ylim(0, 17)\n",
    "        \n",
    "        l, r = find_intersections(xrange, yscan)\n",
    "        error = (r - l) * 0.5\n",
    "        \n",
    "        print (\"{:5}: 1 ± {:.3f}\".format(eff_CP4.columns[ipoi], error))\n",
    "        \n",
    "        ax.axvline(l, ls='--', color='0.8', zorder=-100)\n",
    "        ax.axvline(r, ls='--', color='0.8', zorder=-100)\n",
    "        ax.axhline(1, ls='-', color='0.8', zorder=-100)\n",
    "        ax.axhline(4, ls='-', color='0.8', zorder=-100)\n",
    "        ax.axhline(9, ls='-', color='0.8', zorder=-100)\n",
    "        ax.axhline(16, ls='-', color='0.8', zorder=-100)\n",
    "\n",
    "\n",
    "        ax.set_title(eff_CP4.columns[ipoi])\n",
    "        ax.set_xlabel(r'$\\mu_{%s}$' % eff_CP4.columns[ipoi])\n",
    "        ax.set_ylabel(r'$-2\\log\\Lambda$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# tensorflow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 379,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5,)\n",
      "(29,)\n",
      "[1. 1. 1. 1. 0.]\n",
      "2.658964e-10\n",
      "[[234.56238     22.40912     10.589928     1.7919149   13.467677  ]\n",
      " [ 22.409115    29.724094     1.5157256    0.4012419    2.702509  ]\n",
      " [ 10.589924     1.5157255    8.980461     1.6093585    1.1347737 ]\n",
      " [  1.7919148    0.40124193   1.6093584   21.726883     1.27647   ]\n",
      " [ 13.467681     2.7025092    1.1347739    1.2764698    2.9290712 ]]\n",
      "[[ 0.01213198 -0.00424079 -0.00773632  0.00258832 -0.05000007]\n",
      " [-0.00424079  0.07508886 -0.0017181   0.00202785 -0.05      ]\n",
      " [-0.00773631 -0.00171809  0.240989   -0.01424328 -0.05000001]\n",
      " [ 0.00258832  0.00202785 -0.01424328  0.09579351 -0.05000002]\n",
      " [-0.05000009 -0.04999998 -0.04999997 -0.05000002  1.0000005 ]]\n",
      "[0.11014528 0.27402347 0.4909063  0.30950525 1.0000002 ]\n"
     ]
    }
   ],
   "source": [
    "NSYS_GEN = len(systematics_nsignal_gen)\n",
    "NSYS = NSYS_GEN\n",
    "\n",
    "all_vars = tf.Variable([1.] * NPROCESS + [0.] * NSYS, dtype=tf.float32)\n",
    "\n",
    "sigma_lumi = systematics_nsignal_gen[0]['values']\n",
    "\n",
    "eff = tf.constant(eff_CP4.values, dtype=tf.float32, name='eff')\n",
    "ntrue = tf.constant(xs_x_BR_P4.values * LUMI, dtype=tf.float32, name='ntrue')\n",
    "observed = tf.placeholder(shape=(29, ), dtype=tf.float32)\n",
    "\n",
    "print(all_vars.shape)\n",
    "\n",
    "ntrue_sys = ntrue * (1 + sigma_lumi * all_vars[4])\n",
    "expected = tf.matmul(eff, tf.expand_dims(ntrue_sys * all_vars[:4], 1))[:,0] + EXPECTED_BKG_CAT.values\n",
    "print(expected.shape)\n",
    "\n",
    "chi2 = tf.reduce_sum((expected - observed) ** 2 / expected) + all_vars[NPROCESS] ** 2\n",
    "\n",
    "hess = tf.hessians(chi2, all_vars)[0]\n",
    "\n",
    "cov = tf.linalg.inv(hess / 2.)\n",
    "\n",
    "errs = tf.sqrt(tf.linalg.tensor_diag_part(cov))\n",
    "\n",
    "with tf.Session() as sess:\n",
    "    writer = tf.summary.FileWriter('logs', sess.graph)\n",
    "    sess.run(tf.initialize_all_variables())\n",
    "    print (sess.run(all_vars, feed_dict={observed: obs}))\n",
    "    print (sess.run(chi2, feed_dict={observed: obs}))\n",
    "    H = sess.run(hess, feed_dict={observed: obs})\n",
    "    print(H)\n",
    "    print(sess.run(cov, feed_dict={observed: obs}))\n",
    "    print(sess.run(errs, feed_dict={observed: obs}))   \n",
    "    writer.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 344,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 9.6319737e-03, -6.7407964e-03, -1.0236318e-02,  8.8319968e-05],\n",
       "       [-6.7407945e-03,  7.2588868e-02, -4.2180931e-03, -4.7214961e-04],\n",
       "       [-1.0236314e-02, -4.2180917e-03,  2.3848900e-01, -1.6743278e-02],\n",
       "       [ 8.8319648e-05, -4.7214958e-04, -1.6743276e-02,  9.3293503e-02]],\n",
       "      dtype=float32)"
      ]
     },
     "execution_count": 344,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.linalg.inv(H[0:4, 0:4] / 2.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 306,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = tf.Variable([1., 1.], dtype=tf.float32, name=\"x\")\n",
    "f = (x[0] + x[1] ** 2 + x[0] * x[1]) ** 2\n",
    "hessian = tf.hessians(f, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 307,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9.0\n",
      "[[ 8. 18.]\n",
      " [18. 30.]]\n"
     ]
    }
   ],
   "source": [
    "with tf.Session() as sess:\n",
    "    sess.run(tf.initialize_all_variables())\n",
    "    print(sess.run(f))\n",
    "    print(sess.run(hessian)[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "metadata": {},
   "outputs": [],
   "source": [
    "def createCons(x):\n",
    "    return tf.constant(x, dtype=tf.float32)\n",
    "\n",
    "def hessian_tf(func, variables):\n",
    "    matrix = []\n",
    "    for v1 in variables:\n",
    "        tmp = []\n",
    "        for v2 in variables:\n",
    "            tmp.append(tf.gradients(tf.gradients(func, v2)[0], v1)[0])\n",
    "        tmp = [createCons(0) if t is None else t[0] for t in tmp]\n",
    "        tmp = tf.stack(tmp)\n",
    "        matrix.append(tmp)\n",
    "    matrix = tf.stack(matrix)\n",
    "    return matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 310,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<tf.Tensor 'stack_15:0' shape=(1, 1) dtype=float32>"
      ]
     },
     "execution_count": 310,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hessian = hessian_tf(f, [x])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 311,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1. 1.]\n",
      "[1. 2.]\n",
      "9.0\n",
      "[array([[ 8., 18.],\n",
      "       [18., 30.]], dtype=float32)]\n"
     ]
    }
   ],
   "source": [
    "with tf.Session() as sess:\n",
    "    sess.run(tf.initialize_all_variables())\n",
    "    print(sess.run(x))\n",
    "    print(sess.run(y))\n",
    "    print(sess.run(f))\n",
    "    print(sess.run(hessian))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}