nuts6.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def Leapfrog(theta, r, grad, eps, L):\n",
    "    \"\"\"\n",
    "    SUMMARY: \n",
    "    From Algorithm 1 in NUTS, Leapfrog step in Hamiltonian Monte Carlo/NUTS.\n",
    "    Helper function for Algorithm 6. \n",
    "    \n",
    "    INPUTS:\n",
    "    theta = vector of parameters, len(theta) = number of parameters in the model, \n",
    "    np.array or list\n",
    "    \n",
    "    r = input from other function, momentum in Hamiltonian dynamics\n",
    "    \n",
    "    eps = step size\n",
    "    \n",
    "    L = function from model that returns the log likelihood and the gradient like:\n",
    "    logp, grad = L(theta)\n",
    "       \n",
    "    OUTPUTS:\n",
    "    theta_tilde = update theta vector\n",
    "    \n",
    "    r_tilde = updated r momenta vector\n",
    "    \"\"\"\n",
    "    \n",
    "    r_tilde = r + (eps/2) * grad\n",
    "    theta_tilde = theta + eps * r_tilde\n",
    "    logp_tilde, grad_tilde = L(theta_tilde)\n",
    "    r_tilde = r_tilde + (eps/2) * grad_tilde\n",
    "    return theta_tilde, r_tilde, grad_tilde, logp_tilde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def FindReasonableEpsilon(theta, grad, logp, L):\n",
    "    \"\"\"\n",
    "    SUMMARY:\n",
    "    From Algorithm 4 in NUTS.\n",
    "    Heuristic for choosing an initial value of epsilon.\n",
    "    Helper function for Algorithm 6.\n",
    "    \n",
    "    INPUTS:\n",
    "    theta = vector of parameters, len(theta) = number of parameters in the model, \n",
    "    np.array or list\n",
    "    \n",
    "    L = function from model that returns the log likelihood and the gradient like:\n",
    "    logp, grad = L(theta)\n",
    "    \n",
    "    \n",
    "    OUTPUTS:\n",
    "    eps = value for initial value of epsilon, step size\n",
    "    \n",
    "    \"\"\"\n",
    "    #initialize\n",
    "    parems = len(theta)\n",
    "    eps = 1\n",
    "    r = np.random.multivariate_normal(np.zeros(parems), np.identity(parems), 1)\n",
    "    r = r.ravel()\n",
    "    theta_prime, r_prime, _, logp_prime = Leapfrog(theta, r, grad, eps, L)\n",
    "    \n",
    "    prob = np.exp(logp_prime - logp - 0.5 * (np.dot(r_prime, r_prime) - np.dot(r, r)))\n",
    "    a = 2 * int(prob > 0.5) - 1\n",
    "\n",
    "    while prob**a > 2**(-a):\n",
    "        eps = 2**a * eps\n",
    "        theta_prime, r_prime, _, logp_prime = Leapfrog(theta, r, grad, eps, L);\n",
    "        prob = np.exp(logp_prime - logp - 0.5 * (np.dot(r_prime, r_prime) - np.dot(r,r)))\n",
    "    return eps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def BuildTree(theta, r, grad, u, v, j, eps, L, esto0):\n",
    "    \"\"\"\n",
    "    SUMMARY:\n",
    "    builds tree in NUTS sampler, helper function for Algorithm 6\n",
    "    \n",
    "    INPUTS:\n",
    "    \n",
    "    theta = vector of parameters, len(theta) = number of parameters in the model, \n",
    "    np.array or list\n",
    "    \n",
    "    r = input from other function, momentum in Hamiltonian dynamics\n",
    "    \n",
    "    v = direction in creating tree, value from -1 to 1\n",
    "    \n",
    "    j = height of tree, starts at 0\n",
    "    \n",
    "    eps = step size\n",
    "    \n",
    "    r_theta0 = joint probability of theta0 and r\n",
    "    \n",
    "    L = function from model that returns the log likelihood and the gradient like:\n",
    "    logp, grad = L(theta)\n",
    "    \n",
    "    \n",
    "    OUTPUTS:\n",
    "    theta_minus, r_minus, theta_plus, r_plus, theta_prime, n_prime, s_prime, alpha_prime, n_alpha_prime\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    if j == 0: \n",
    "        # base case, take one leapfrog step in the direction v\n",
    "        theta_prime, r_prime, grad_prime, logp_prime = Leapfrog(theta, r, grad, v*eps, L)\n",
    "        \n",
    "        esto = logp_prime - 0.5 * np.dot(r_prime, r_prime)\n",
    "\n",
    "        n_prime = int(u <= np.exp(esto))\n",
    "        \n",
    "        s_prime = int(esto > np.log(u) - 1000)\n",
    "\n",
    "        alpha_prime = min(1, np.exp(esto - esto0))\n",
    "                                                   \n",
    "        return theta_prime, r_prime, grad_prime, theta_prime, r_prime, grad_prime, theta_prime, grad_prime, logp_prime, n_prime, s_prime, alpha_prime, 1\n",
    "    else:\n",
    "        # recursion, build left and right subtrees\n",
    "        theta_minus, r_minus, grad_minus, theta_plus, r_plus, grad_plus, theta_prime, grad_prime, logp_prime, n_prime, s_prime, alpha_prime, n_alpha_prime = BuildTree(theta, r, grad, u, v, j-1, eps, L, esto0)\n",
    "        \n",
    "        if s_prime == 1:\n",
    "            if v == -1:\n",
    "                theta_minus, r_minus, grad_minus, _,_,_, theta_doub_prime, grad_doub_prime, logp_doub_prime, n_doub_prime, s_doub_prime, alpha_doub_prime, n_alpha_doub_prime = BuildTree(theta_minus, r_minus, grad_minus, u, v, j-1, eps, L, esto0)\n",
    "            else:\n",
    "                _, _, _, theta_plus, r_plus, grad_plus, theta_doub_prime, grad_doub_prime, logp_doub_prime, n_doub_prime, s_doub_prime, alpha_doub_prime, n_alpha_doub_prime = BuildTree(theta_plus, r_plus, grad_plus, u, v, j-1, eps, L, esto0)\n",
    "\n",
    "            # Use Metropolis-Hastings\n",
    "            prob = n_doub_prime / max(n_prime + n_doub_prime,1)\n",
    "            if (np.random.uniform(0, 1, 1) < prob):\n",
    "                theta_prime = theta_doub_prime\n",
    "                grad_prime = grad_doub_prime\n",
    "                logp_prime = logp_doub_prime\n",
    "            \n",
    "            ind_1 = int(np.dot(theta_plus-theta_minus, r_minus) >= 0)\n",
    "            ind_2 = int(np.dot(theta_plus-theta_minus, r_plus) >= 0)\n",
    "            s_prime = s_prime * s_doub_prime * ind_1 * ind_2\n",
    "            n_prime = n_prime + n_doub_prime\n",
    "            alpha_prime = alpha_prime + alpha_doub_prime\n",
    "            n_alpha_prime = n_alpha_prime + n_alpha_doub_prime\n",
    "        \n",
    "        return theta_minus, r_minus, grad_minus, theta_plus, r_plus, grad_plus, theta_prime, grad_prime, logp_prime, n_prime, s_prime, alpha_prime, n_alpha_prime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def nuts6_dual_averaging(theta0, M, M_adapt, L, delta = 0.6):\n",
    "    \"\"\"\n",
    "    SUMMARY: \n",
    "    Implemented from Algorithm 6: NUTS with Dual Averaging.\n",
    "    Needs helper functions BuildTree, FindReasonableEpsilon, and Leapfrog.\n",
    "    \n",
    "\n",
    "    INPUTS:\n",
    "    \n",
    "    theta0 = initial values for values of parameters in model. len(theta0) = number of parameters\n",
    "    \n",
    "    M = number of samples desired\n",
    "    \n",
    "    M_adapt = the number of steps for the burn-in,\n",
    "    also how long to run the dual averaging algorithm to find the appropriate epsilon\n",
    "    \n",
    "    L = function from model that returns the log likelihood and the gradient like:\n",
    "    logp, grad = L(theta)\n",
    "    \n",
    "    delta = target HMC acceptance probability.\n",
    "    default value of 0.6\n",
    "    is a value between 0 and 1\n",
    "    \n",
    "    \n",
    "    OUTPUTS:\n",
    "    \n",
    "    samples = np.array matrix of samples of theta from algorithm\n",
    "    dimensions of matrix are M x len(theta0)\n",
    "    \n",
    "    burned_in = same as samples matrix with burn_in samples removed\n",
    "    dimensions of matrix are are M-M_adapt x len(theta0)\n",
    "    \n",
    "    \"\"\"\n",
    "    # initialize samples matrix\n",
    "    # put initial theta0 in first row of matrix\n",
    "    parems = len(theta0)\n",
    "    samples = np.empty((M+1, parems))\n",
    "    samples[0, :] = theta0\n",
    "    logp, grad = L(theta0)\n",
    "    eps = FindReasonableEpsilon(theta0, grad, logp, L)\n",
    "    mu = np.log(10*eps)\n",
    "    eps_bar = 1\n",
    "    H_bar = 0\n",
    "    gamma = 0.05\n",
    "    t0 = 10\n",
    "    k = 0.75\n",
    "    \n",
    "    for m in range(1, M+1):\n",
    "        # resample\n",
    "        norm_samp = np.random.multivariate_normal(np.zeros(parems), np.identity(parems), 1)\n",
    "        r0 = norm_samp.ravel()\n",
    "\n",
    "        esto = logp - 0.5 * np.dot(r0,r0)\n",
    "        # resample u ~ uniform([0, exp(inside)])\n",
    "        u = np.random.uniform(0, np.exp(esto), 1)\n",
    "\n",
    "        # initialize minus's and plus's\n",
    "        theta_minus = samples[m-1, :]\n",
    "        theta_plus = samples[m-1, :]\n",
    "        r_minus = r0\n",
    "        r_plus = r0\n",
    "        j = 0\n",
    "        grad_minus = grad\n",
    "        grad_plus = grad\n",
    "        \n",
    "        j = 0\n",
    "        samples[m, :] = samples[m-1, :]\n",
    "        n = 1\n",
    "        s = 1\n",
    "        \n",
    "        while s == 1:\n",
    "            v_j = np.random.uniform(-1,1,1)\n",
    "            if v_j == -1:\n",
    "                theta_minus, r_minus, grad_minus, _, _, _, theta_prime, grad_prime, logp_prime, n_prime, s_prime, alpha, n_alpha = BuildTree(theta_minus, r_minus, grad_minus, u, v_j, j, eps, L, esto)\n",
    "            else:\n",
    "                _, _, _, theta_plus, r_plus, grad_plus, theta_prime, grad_prime, logp_prime, n_prime, s_prime, alpha, n_alpha = BuildTree(theta_plus, r_plus, grad_plus, u, v_j, j, eps, L, esto)\n",
    "            \n",
    "            if s_prime == 1:\n",
    "                # Use Metropolis-Hastings\n",
    "                prob = min(1, n_prime/n)\n",
    "                if (np.random.uniform(0,1,1) < prob):\n",
    "                    samples[m, :] = theta_prime\n",
    "                    logp = logp_prime\n",
    "                    grad = grad_prime\n",
    "                    \n",
    "            n = n + n_prime\n",
    "\n",
    "            boolean_1 = int(np.dot(theta_plus-theta_minus, r_minus) >= 0)\n",
    "            boolean_2 = int(np.dot(theta_plus-theta_minus, r_plus) >= 0)\n",
    "            s = s_prime * boolean_1 * boolean_2\n",
    "            j = j + 1\n",
    "        \n",
    "        eta = 1 / (m + t0)\n",
    "        H_bar = (1 - eta)*H_bar + eta * (delta - alpha/n_alpha) \n",
    "        if m <= M_adapt:  \n",
    "            eps = np.exp(mu - np.sqrt(m)/gamma * H_bar)\n",
    "            eta = m**(-k)\n",
    "            eps_bar = np.exp((1-eta) * np.log(eps_bar) + eta * np.log(eps))\n",
    "        else:\n",
    "            eps = eps_bar\n",
    "            \n",
    "    burned_in = samples[M_adapt+1:M+1, :]\n",
    "    return samples, burned_in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.7300000000000004"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# generate random poisson data\n",
    "X_pois = np.random.poisson(5, 100)\n",
    "def L_pois(theta):\n",
    "    grad = sum(X_pois) / theta - 100\n",
    "    logp = sum(X_pois)*np.log(theta) - theta*100\n",
    "    return logp, grad\n",
    "np.mean(X_pois)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.1600574868954583"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "# run on poisson data\n",
    "samples_pois, burned_in_pois = nuts6_dual_averaging(np.array([2]), 5000, 2000, L_pois)\n",
    "np.mean(burned_in_pois)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# multivariate normal data\n",
    "cov = np.identity(10)\n",
    "X_MVN = np.random.multivariate_normal(np.ones(10), cov, 100)\n",
    "\n",
    "def L_MVN(theta):\n",
    "    cumsum = 0\n",
    "    grad = 0\n",
    "    for x in X_MVN:\n",
    "        cumsum += ((x - theta).T @ np.linalg.inv(cov)) @ (x - theta)\n",
    "        grad += (x - theta).T @ np.linalg.inv(cov) \n",
    "    logp = -0.5 * cumsum\n",
    "    return logp, grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 39.7 s, sys: 238 ms, total: 40 s\n",
      "Wall time: 40.2 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# run for MVN\n",
    "samples_MVN, burned_in = nuts6_dual_averaging(np.zeros(10), 5000, 1000, L_MVN)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# subset just accepted theta's\n",
    "unique_samples = np.array(samples_MVN[0,:])\n",
    "for i in range(1, M):\n",
    "    if all(samples_MVN[i-1,:] != samples_MVN[i,:]):\n",
    "        unique_samples = np.vstack([unique_samples, samples_MVN[i,:]])\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEJCAYAAAB2T0usAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXdgVFXah5/pKTOphCSUSL9AAAFpshZEBQEFV1Tsn9hQ\nBLvrKq66dkVQQAVUELGsKFhQqUqxB+n9QiAQSgghdSbJ9Pn+mMzNzGSSTCqo5/nHuTe3nInhvOe8\n5feqPB4PAoFAIPh7oz7dAxAIBALB6UcYA4FAIBAIYyAQCAQCYQwEAoFAgDAGAoFAIEAYA4FAIBAA\n2rreIEnSHEAty/JdNVzTGpgBDAPKgcXAw7IsW+s7UIFAIBA0HXXaGUiS9CxQrRGouEYPfA/EAecC\n1wKXA6/Wc4wCgUAgaGLC2hlIktQemAekA4drufxGIBkYKMtyScX9TwH3NGCcAoFAIGhCwt0ZDAay\ngZ7AoVquHQas9hkCAFmWP5BleVC9RigQCASCJkdVVzkKSZLWAvurixlIkrQJ+AGwAjcBHuAL4ElZ\nlm0NG65AIBAImoI6B5DDIAa4A1gGXA20Bt4CkoBbm+B9AoFAIGggTZFa6gDygZtlWd4sy/I3wIPA\nzZIkxTfB+wQCgUDQQJpiZ3AMKJdl2d//tBtQAe2Awupu9Hg8HpVK1QRDEggEgr80DZ44m8IY/ATc\nIUmSRpZlV8W5noCTWoLPKpWKvDxzEwypeUhKMonxn0b+zOP/M48dxPibiuxcM7OW7ABg8tiepCWb\nQl6XlBT6fF1osDGQJEkHJAAFsiw7gDnAJGBhRV1CW7w1Bh/IslztrkAgEAgEgcz5dCND1i8hN7Yl\ns4CpEwc32bvqEzMITj8aDBzHW2CGLMsngQvwGohNwEfA58DE+g9TIBAI/l5oN2Tw7JyJ3PrzQi7c\n+2PTv6+uN8iyPDToeD2gCTq3FxjRsKEJBALB35CyMqJfepbId2aj8niwaQ18fdmdTB7bs0lf2xQx\nA4FAIBDUA91vv2C6fyKaQ1nKOftLr/DQ/41r8ncL1VKBQCA43VgsGB9/hLgxI1BnH8Yd783Ct102\nCust45tlCMIYCAQCwWlE99N6EoYMJnLeOzi7SFheeBWV2YwrOQXz629CM6XbCzeRQCAQnAZUFjPR\n/32KyA/m4dFoKLv/YcrumUTcyEtQOZ2YZ83Bk5jYbOMRxkAgEAiaGd3aHzA9fB+ao0dwduuOecbb\nOHv3xfjQZLQHD1B29yQcQ4bW/qBGRBgDgUAgaCZUJcVEPz2FyI8X4tFqKX34McoefBT0evTfLiXy\now9wpvekdMrTzT42YQwEAoGgGdD/sArjw/ejOX4MZ3pPSmbOxtWzFwDqnOOYHp6MJyKCkjnzwGBo\n9vEJYyAQCARNiKqoEON/Hidi0Sd4dDpK//UEZfc/DDqd9wK3G9OkCagLCzG/Mh2X1PW0jFMYA4FA\nIGgi9CuWYXz0ATS5J3Cc3QfzjLdxdU8PuCby7Vnof1qPbfgIrLfefppGKoyBQCAQNDqqgnyMUx4j\nYslnePR6LFOepvze+0EbOOVqt28l+qVncbVMxvz6W82WRhoKYQwEAoGgEdF/uxTTYw+hzjuJo+85\nmGfMDu36KS3FdPftqBwOzDNn42nRovkH64cwBgKBQNAIqE6dwvj4I0R8/QUegwHLU89Rfve9VXYD\nPoxPT0GbuZ+yCRNxDL2kmUdbFWEMBAKBoCF4PBiWfonx3w+jzs/H0X+gNzbQqXO1t+iXfUvkwvk4\nu/egdMozzTfWGhDGQCAQCOqJ6uRJTI89hOG7pXgiI7E89xLld9wNGk2196hP5GB6aFJlGmlERDOO\nuHqEMRAIBIK64vFg+OJzjE88irqwEPugwZjfeAt3h4413+d2Y5p0N+qCAswvvYara7fmGW8YCGMg\nEAgEdUCdewLjow9gWLEMT1QU5pemYh1/J6hr1/2MnPMW+h/XYrt0ONbb7myG0YZPnY2BJElzALUs\ny3eFef23QFRwUxyBQCD4U+HxYFj0Ccb/PI66uAj7eRdgnj4Ld7v2Yd2u2bGd6Beewd0iCfMbb5/W\nNNJQ1EnCuqKncVhGoOL6CcDIug5KIBAIziTUx48Rc+M1xNx3DzidmF99neLFS8M2BJSVEXP3bd40\n0lmz8SQlNe2A60FYOwNJktoD84B04HCY93QCXgB+rffoBAKB4HTi8RDxyYdEP/UEanMJ9gsv8u4G\n2qbV6THGZ6ag3b+Psjvvxn7xsCYabMMI1000GMgGrgMW1XaxJElq4APgZUACaomqCASCvwtFFhuL\n1mQCMG5oJ+KMzS/KFg7qo0cwPTQZ/bo1uE0xmKfPwnrjLXV27+hXLCNywTyc3bpT+p9nm2awjUBY\nxkCW5Y+BjwEkSQrnlicAtyzLr0mS9G79hycQCP5qLFqTScbuXOV4wuj0M8tAeDxELHyf6P/+B7XF\njO3iS7G8NgN36zZ1fpQ69wSmB+/FYzBQMmf+GZNGGopGzyaSJOkc4EGgX2M/WyAQ/DUJZSBOB+rD\nh7y7gZ/W446No2TmbGzjbqhfsNftxjT5btT5+ZhffBVXt+6NP+BGpFGNgSRJBmAh8KQsy1mN+WyB\nQPDXYNzQTiE/n1bcbiLmvYPxuadRlZViGz4Cy9Q3cKek1vuRke/ORr9uDbaLL8V6+4TGG2sTofJ4\nPHW6QZKktcD+UKmlkiRdAKwFSgGfKTXgzVoqA7rLsny0hsfXbTACgeC0c/BYMc+/nwHAk+MH0qF1\nbJ2fUVBiZd7SnQDcProHCTHN6E45cABuvx3Wr4eEBJg5E26o527Ax7ZtMGAAxMbCjh2QnNx44w1N\ng/NUG9tNlAEEC3K8BKQBNwDHa3tAXp65kYfUfCQlmcT4TyN/5vH/mcf+7Hu/k19iVT5PnTg47Huz\nc83MWrIDgMlje5KWbMJlc5CX52iSsQbgchH53hyiX3wWVXk5tpFXYH5lOp7kZDhlqf9zy8uJH3cd\nWrud4jfewq6Ogib+f5uUZGrwMxpsDCRJ0gEJQIEsyzbgYNDPS4By4TYSCATBzFqyQzEks5bsqNaQ\nNHaAWZO5H9P9E9H9kYE7MRHV++9TctGIRikEM/73SbTyXspvvwv7pZc1+HnNRZ2KzioIduUMxrvi\nP7fhwxEIBH82Jo/tSVJ8JIkxEUwe27NJ3uELMGfszlWMQjBFFhtzl+5i7tJdFFlsoR/kchH51kzi\nh/4D3R8ZWMdcRcFPf8C4cY1iCPSrlhM5/12cXbtheeq5Bj+vOanzziBYVkKW5fVAtRJ9siyfWQIc\nAoGgUUlLNjH/yWH1cnNNHtuTGZ9vp9TqIDUxiiKLrdZVv93hYu7SXUDgLqG2jCTNPhnT/feg27QR\nd4skSt56F/sVY+o85upQ5eZieuBePHo9JbPnQWRkoz27ORBCdQKBoF74u24mXtO7Xs9ISzbRJS2O\njN257MwqYNGazJBppf5ZR3aHq25pqE4nkW/PJHrqS6hsNqxXXYPlhVfxJCbWa8whcbuJuf8e1KdO\nYXnuJVzpPRrv2c2EMAYCgaBe+K/EDUt3cuvwsApS60Wc0aBM+r5dQTChUlY1e3Z7dwNbt+BqmYxl\n6hvYR4xq9PFFzpuLfs332C+6mPI772n05zcHwhgIBILTwq6sfGYt2YEHkNrGEWcyVKk7CBU4rq5O\nwd9g4HAQNe0Voqa/isrhwHrt9VieewlPfEKjfw/N7l1EP/sU7sREzDNnhyVlfSYijIFAIKgX/hPx\n7aN74LLVLR101pId2J1uALJySphz45Aq14SKAwRM+iHQ7NzhzRTasQ1Xaissr73RdFk95eXE3HM7\nKpuNknkLcSenNM17mgFhDASCvyH1TdUMvs83KSfERDRPbUBN2O1EvT6VqBnTUDmdlN9wM6X/fQFP\nbFyTvTL6uafQ7tlN+fg7sA8b0WTvaQ6EMRAI/kKEO8nXVwuotvvqYmQmj+0ZUHAWinClK7TbtmC6\nbyLaPbtwtW6DedpMHEMvqf0LNQD99yuJem8uzi4Slqefb9J3NQfCGAgEfyGaUvCtyGIj82hxo70/\nvX0icx4Zohz7YgjgNQ7p7RNrdQlhs3ljA7NeR+VyUX7LbZQ+/SweU0x4X6qeqE6exHTfxMo00qio\nJn1fcyCMgUDQRJxRssxAxu4TvPvNbgBuuLRSNSZcsbhFazKVauHYaD37sot49O1fvUVnYcgh+H4f\npVYHx06WolarFAkKCIwhzFqyI8BQhEK7eSOm+yeilffiSjsL8/RZOC6o+Z5GwePB9MBE1KfysPz3\nRVw9ezX9O5sBYQwEgibidMgy1+RWefeb3bgr9AM+Wb2f9x6rvi15bYas3OYMmLgX9GhFdq6ZfdlF\n6DRqPB43+7KLyM41ExOt58OVMruyCpR7fFQnQeFwuXn07V9JSzZy83Ap8P3l5US/+iKRs2ehcrv5\n46KxxM98jdjkRqwbqIGI+e9g+H4V9gsvonzCxGZ5Z3MgjIFA8BeiVrdKEKGE4iC0IfM3Lvuyi7AH\nST7MWrKDQr9zhRYbs5bsoFObWLbsP1XrWCaP7cnMxTtwuNx4PJBfYiW/xIpep1G+U/7KNcQ+MImo\n/KPkt2jFaxfdw862PRmYcYIJo5veGGj27sH4zJO4ExIwz5rzp00jDYUwBgJBE9Fcuv3ZuWZF0qFL\n2zhuG9VNWUn7r/BvuLQzn6zeD8CdV3QnO9fMswv+UHYLNQnFQaWhKbLYmL9sDxarAzweUltEsUU+\nSaHZGtZ4Y6P1qFUqxU3ko3WSMaQ8UObRYo4cOol9ypMMWO3tuvt1n8v59MJbsKj1Yb2zUbBaiZlw\nmzeN9J0FDep1cCYijIFA0ESEs0r3BU1VKhWTrupBevvwVrf+k/zew4UUl9oBqkg6+K/wM48W89q9\n/yDOaKhiCIKpyZAtWpPJzoMFyvHOgwXsf39DlWepqFjd77YSFaFFp1GjVqlol2qq6vqpeG6wGwkg\nefcmWs28jZSC4xyLa8WM4ZPZ07obOq2agV2SQo6xKYh+/mm0e3ZRfsttTVLFfLoRxkAgOI3UFDSt\nzoUDgZO8uhqxzW9/ywpw9eSXWHnozV/QadVE6DRVJm+ny012rpm0ZBMlpXYlc6ik1F5r8NsRYhL3\nf3yZ1UliTAT5JVYK99vYlVWAKUpf5Xv50GvVpLc00P+jmYzaugyPSsUX51zJx4Ovx67zjkVqG1fF\n2Na0S2oIujWriXpnNs5OnbE8+2KDn3cmIoyBQHCG4PF4AtQ4w9X6N0XqKS7z7gxio/TKKvmL9aFb\niDic7pCTd3GpnRmfb2fapH8w4/Ptiv//hYUbeeWewcqOYl92EWoVijHRa9WKQVMRXrtCu9NNfomV\nGZ9vp12qiexcC6ktoujTuQV6nYZrVEdI+tcEEvJzON6iLdMvmYTcSiLeaCA91YRepwm5G/CPW9Qk\nfFcXVHl5xEy+B49Oh3nu/L9EGmkohDEQCE4j40d25Z1vduPxgEGvCVjJh0LR8/F46JoWR6zRwHk9\nU1iwXMbpdIMKXli4iVtH1Cwap1JBqI635nI7d7+2LsBd43B5lEk1OEgc7nMBTFFa0pJbBGQVlVod\nSnA5v8TK+R2MTM74H5EL5uFSqfm8/1hWDfs/TpR7r1erVYohCLXid1fn96ovHg+mB+9FnXcSy9PP\n4+x5duM+/wyizj2QmxjPn7X1H/y5WxeCGP/pYO7SXSENwMDuyYwYmMasJTuUSV6rUVNsseKbp3Ua\nFT06JIZM2awJnVbNWclGMo+VhH1PvNGAWq2ipNSGw1U5Z/jvEMLB/3sBpCZGsTPLG384+/A2Hlk3\nm7j8E5xs3YGXLpxIZkonenRIIOdUGSWldhwu7/eMjdbz9Pj+VQzC9EVblefFRul5+rb+IYPpwcak\nur+diPnvYvr3w9jPH0Lx51+dsdlDSUmm5u+BLEnSHEAty/JdNVwzDvg33n7Ix4F5wFRZlsP/ixUI\n/obotWr6dElSJqupEwfz6Nu/Ku4if1xuT1gpm/60bhGF1e7mUE74hkCnVQfsBnQatTIp13UhXlru\nIC3ZpLi8snPNFBw5ydUr3uHS7avwaDSUPvQoJXfeT+IvR0jE27+g9c4NjMzezjd9RlEUHU9xqT2k\nCyg6Uqd87touPmDC94+zZOzOZcLo7pRaHXy0ypthddOwzgzt21a5XiPvxfjMFNzx8Zjf/GulkYai\nTsZAkqRngbuA92q4ZgTwEXAfsALoU3G9Fnih3iMVCM5g6lttPG5oJ/ZlFymTbXr7hBp93DFROqx2\nFwAdWsWwN7uoTuPMyS+r0wSeGBNBWrIxwOg4XfVf0x05aWH6oq3IR4rweDz0ytzE86vfIsmST3Zy\nB1ZOeIohtwZm6vT+bTljlryExuPm8i3fsmTAWL7qOzrgmiKLjQ9XymQeK0YFaLVqzuuZEvDzYCkN\n/yI8gI9W7a80BjYbMXffjspqpWT2PNyprer9nf8shGUMJElqj3d1nw4cruXyCcDnsizPrjjOkiSp\nOzAeYQwEf1HqW20cZzTwn1v78dUvh7DZnJzXM4VH3/4VqMwg8hd0u3WExM87TgAwYmAaz7z/R53G\nqarJqR8Cnw8+NlqvpK/WdLdGBa4aLjCX2dmZVUC01cId6+dzya41ONUaPj73OhYPGIvToiO3wqhm\n7DrBVRu/5KqfFlIeZeKn869i4E9fcfMvH3P5rtUUtn6KuV+5Qa3G7nAFGCyH082C5TJTbjmHRWsy\nyTxaXGV3VZNRjH7hv2h37aD85luxj7qi5l/SX4RwdwaDgWzgOmBRLdc+B5QGnfMA8XUbmkDw9yDO\naODRm/qxaefxkEVgackmZVJbsFwOmNQiIzSUW707BbUK4k0RIV1KPlx19OsUWmwU7q8+YOxDp1Wj\nglpjF24P9D+wgXu/n01iaSGZLTswY/h9ZLdsFzg5u93csX4+YzZ/wylTC6bf/Dyebt2Z12kY12xY\nwpVbviH+35O4JqUz8y68jZPd+lR5l7nMzgsLN1X7+9CoVQG/jzYtowGwL19B0pw3yU9Oo/CxZ/Cc\nYRpTTUVYxkCW5Y+BjwEkqeYsBVmWN/kfS5IUA9wNLK/fEAWCM5/GqDaetWRHwIRYaLYqef/BOw8f\nNpsr4HjKLecw/7s97DtSpMg61Ib/qr++hEpVDcZYbuaude9x0Z71ODRaFv7jRr7o90/cWi13XdGd\nDXtOcijHzIGsPO5fMZPem78nO7EtT1/1FKciktAcLsBliOaD82/hpwGjGPv9+1wg/8yrix5ne5+L\nWHDeLRyKSFTiGb701eoITp4pt7pQ5ecTc989ONRanrvkPlr8ngNQ511fdq6Z1z/bhrnMjilSz4Pj\nzg5ZT3Em0aSppZIkRQJfARHA4035LoGgrvj8zNm5ltCCaHWgrppAwRSUWDGXBU7Ibg9K3r8/ahV0\naRuH3eFCo1bjrpj8TFF64owGHhrXm11Z+cz4fBvOMIyBpbx2Q6B1OnBqdbVeByiSEv5z7aD9vzPx\nhznElxWxL7kzM4ZPJrtFWsWF8MOmo2QeKyHSXs7jS1+md/Y2jnTsyWPD/oUl0juJ+ocqDhoSmTrq\nEZb2uZy7f3qfXlvW8sr2n1ja53KW9h1NgTF0e8vYKD1arTcQ3CIuAtkv5pLWMhrTg/diKM7n/fNv\n4UByJ1qE9Y2rMmvJDsXAFpfZa5X6OBNoMmMgSVIi8A3QFbhEluUj4dwXjhTumYwY/+mlLuNfsFIO\nyHE3VbhrTgdTP9oY0sVSZnOSlGRi4jW92fri99gcLtwe2HekGHfFbOsLAxSX2nn4rV948Pq+vP7Z\n9rADxbXFg6/77VNu/O1T7rhtDrlxtbd1jImq3GnElBUzYe27XCD/jF2j4/3zb+Grc8Zwyc4feObL\nZ7EYjBRHxVAcGYsKDxfIPyvPWTtqPHElxdh0BhzaqhpE8SYDh85K57G2rzJw14/c+tNCxm78iis3\nLSWj4wCWn30Z29J60bV9IqeKvTuEJ8cPpEPrWAC2yCd5/v0MnE4P3TvE85htK4YVy3CcfyGn7riX\nC9Rqbh/dAwDD0p2At71nQkxErb8DtUZV5fhM/7dV5zoDSZLWAvtrSS1tB6wCooFhsizvCvPxos7g\nNPJ3G39wjv/A7slNIjMdnGlUUmqvIjOxYKXMj1uOAWCK1GEu97aQ7Ng6hhaxkZRaHQF6QDWh06rD\nctuEy1evX4XG42Zlj0t5c9i9Yd/3j32/cPcP7xBXXsyeVImZwydzNKENWqeD+e/didFqwabTY7SV\n1fic7IQ23Hvrm1XO+1c+A+gdNobuWcdl21bQMc9bfZ3Voh0Ln36fiePOqXK/f8puL0cez793Hx6D\ngcJ1v+Fu1bra8YSTOdbcbqLTUmdQG5IkJQFrATtwrizL2Y39DoGgMRg3tBN2h0txE9Xm669v+miw\nv98/ldTnBrp9dA9sNifgzcX3FU6dKrRyoJrisMSKFWqwX7wxDQHAil7DGbVtOcN3rg7LGMSWFXH3\nD+9w3v5fsWn0vHfheL7pczlutQaAC+SfiC8rYkm/K1lwwa1Ix2Ve+/QxAPamdOG73iOJsZZw/t6f\n6XpiH0cS21Z5h06jwqDTBBgDu87Ail7DWdFzGJ1P7OeBlTNpf+oQJkfNxkbrdHDX5y+jKi+n5M13\nKIhpwTufbGbfkSI0ahX3Xd0rQEAwnMyxtGQTr08+r9bf1ZlEg42BJEk6IAEokGXZAbxdcTwUsEmS\nlFxxqUeW5ZMNfZ9A0FjEGQ1MHht+l6rGaFazL7uIotLK7JxSq3cHkBATochDv7CwMgej3O6s9lnG\nSC1FFnuNEhCNwf7kSiNpKi/BHFlNS0mPhwvkn5iw5l1irGZ2te7GjGGTyYlvFXDNmC3f4FKp+a73\nKM4+vI0nlr6EGxVzh97Jst4jAUiwFHD9b59Sqo/inSF3BLwmNsq70n79s23KuYDMIJWK/aldyDMl\nkVZwlNGXeWWyg425L2V33KoPOevEAU5ddT2PHEnCPOc3xci4XR5mLtnB3Fq6rtWFM60Dno/6GIPg\nP7vBwBrgIkmSNgD/xKtXtcHvGhXgBJpRfFwgaBx8geZdWeG5afzvW7QmE7vDRbzR4E3TDNL1advS\nqHzelZXP9M+2BUzs7hpm+cO5ljqNp77sbNtD+Tx4/2+s7DW8yjXxlgIm/jCHQQc2YNUaeGfIHXzb\nZyQeVWDVbo+ju+hwMoufuwym6/G9PLhiBqjg1csf4ZculYHyu9a+i9FWxluX3EOBKVDWOyUxipz8\n0oAMqOA0UYBIRzkulZqF6w6j1x8LqEXIPFpMpzaxPN/BQptfl0DnzkyRrg1d6R0UVKlv5pjv7yG4\n5qE5OuCFQ52NgSzLQ4OO1wOahjxTIDiTWbQmM6CgKTEmIqxJwH8nodeGljIoKKk0Dl4BusCfO2uq\n4GomcmOT+aN9P/pnbWTS97MDjYHHw0V71nHn2nmYbBa2t+nBrGH3ciIudOOXMZuXAuBWqfnXsmmU\n6qN4YfTj7EirbHKj8rgZlJnBSVMSK3teWuUZ+44WIR8JrLyePLYnM5fsCHCRRdrLKddHsiUzH/Dq\nK/nIL7Fi33SSmE8exKPV8tFNU8gr1QBVXWyd28QFHNc3c6y69OAzhb+22IZA0AQUma28sHAT2bnh\nB6vtTjfxRgN9OrdA52cYSq0O5i7dxY9bjtZJbK65eeuSu5XP0VbvjiTBnM9TX73AQytmoHE7efvi\nCTx5zbPVGoKUohwGHfA6DC6Qf6YgOp7Hr30+wBAAeFRq9rTqRgvzKeJLq8pthNosrdl8jAidJuBc\nlL2ccl1l5o+lwiXne8jkVW8RU5zP2tF38llxXIAh6ZoWx8Duyd6kgjGNv3JPjIlgYPfkZmnKEy5i\nFS8Q1MK4oZ3YsCdXmYRcFf15q8sd9zVYsVgdAf58b5vHXkrTGnOZHbvTTcbuXP7Yc+auGAHyTZUZ\n9w+snElGxwHcvn4+RlsZW9LO5s1LJ3IyNrmGJ8CVm5Yqn4/FteLpsU+TW809W8/qRY9ju2iXd6ja\nmgF/tuw/hU4TuLaNtJdTFBWrHPtP9sN2rObcAxlkde7Nb5fdCHsDBf/2HSlCq1EzeWzPRvPpB7uX\nzpRYgQ9hDASCWogzGtBp1CFX7qGCgdVp/jtdbtZsPsInq/fj8UB0pLYyUOm32g1VtHUm8NyYJ/jP\n1y8y6MAGBh3YQJk+klmXTGRVz0sJ2bzYj3hLAaO2eUUIDiS156mxz1DiN1EHczKmJQAtzXl1GGHg\nLyzCYaVc762LMOg02Bzeau3WBce4c908LIZoTkx7C3ljCQadhvapJiINWrZmnsLt8e7mGjN43NDC\nxKZGuIkEgjCYPLYneq0anUaNKVJHYkwEk8f2VPzAGbtzFaNQHcWldj5atR+3xzttWcpDZwqltTSS\n3q721XCz4vEQX1oYcOre/5vJql7DajUEMeUlLHznNuX43+NerNEQQKUxSCoJLwFRp1UTaahc22pc\nTvQuB+X6SADFEGhdDh5ZNo0Ip405w+9l7mYzhRYbNoeLvCIrk8f2wv/bBAePq6PIYmPu0l3MXbqL\nohqa/5zJiJ2BQBAG6e0TA/oThyLzaDFFFhuTx/YMaBtZVwrNNlISo+t1b1PQsjiXyavfonf29oDz\np4y1izUklZzk+cVPK8c3T3gfa8UEXRMnY7yN7luW1L4z8PVWDg4eA5TrAt914y+f0OnkQVanX0z+\npZdDftUahM5t4pQAdXDwuDoaI+34dCN2BgJBAxg3tJOSpZJfYuW5BRuJidYzbdI/eHjc2V41T5W3\nqjg+TB+xpdzBuKGdlKKy04XK42bk1mW8ufB+emdvZ0OHfrww+t/Kz/sc3lrj/amFOUz9379pVeQV\ne1vR81KKosMTL843JuJSqcMyBqktqvYkjnR4jYFVX/k77JW9nas2fkVOXAoLLrkT8EqCJ8ZEkBQf\nyeSx3kD2hDHpYQWP/XcDdoer2uv+LIidgUDQAOIq2kH6KKyIIUwYnU56+8QAf3N1LS4hsJG8KdJb\njtOQJjINJaUoh8mr3qLX0Z2YDUamXfYA67pdGHDNDb99ypZ2VaWjAQwOK49/8zKJpZW1GV8HNaSp\nCbdaQ75dHBJeAAAgAElEQVQxMSw3USiZDmVnULELMZWX8OCKGXhUKqaOeIgSTQQ7swrYfaiA9q1i\nKLLYmfH5dlonRRMdqasS4A0VG/LfDfTp3IKB3ZOVn/8ZETsDgaCBpCUbA45DrRKLLDbsDpdXNVNT\n1cfeo2MiEXrvP8fiMjtT/7e5wbLS9UHlcXPF5m+YtfABeh3dyW8dB3Lv/81kXfch3tiASsXXfS4H\noGuOTPejIWTHPB7u/X427U8dZm9KFwA2tevL0RCyEjVxMiaJhNJCNK7qq7Crw2cMynSR4PEwafXb\ntLDk87/B17M/tYtyndsDB46VkF9spdBiY2dWARm7c/n3nN8CUodriw3pdRrFCCxak/mnjBsIYyAQ\nNJCbh0u1uoB8hWvFZfaQhWRRkTqs9sqdQE5+eaOPszZSC4/z0mdTuGvdPOxaPa+OfJgXR/+bwqDU\nzt86n6t8vjZjcZXnjNi+gov2rEdO6Uxehe//67517xaWF5OExuOmhSW/zvdGOrwVvuX6SC7d+T2D\nM39nV+vufN7/qrDutzvdiphgdYwb2klxJ40b2qlOyQRnIsJNJBBUQ7gaMnFGA13S4iqrjYOKnyD0\nbsGfjJ0nGjja+hNttTBs5/fc+MsnGFx2fu4ymDlD76I4KnTwdE+rrhRGxRJfVsw5h7fQ6UQmmSne\nVXGXnH3cuXYexZExzL3oTl5Z9ATZCW3YclbvOo/LP6OounqE6vDtDOJLi7hmwxIshiimjXhQEcvr\n0SGBYydLKbU6aJtspMhix+3yYCl3KM1x/BkxME3poTxioLcPQ7ipor66EqhUqT0TEcZAIKiGcDNE\niiw2Sq0O9Fo10RE6endKrNLH2J/GlpiuLxH2cv77xbN0P75HObe7VTcW9x+LOaL6Ccut1pDRcSCX\n7VgFwLUZn/PimMc55+BGHln+Ohq3i6kjH2bAwY3o3E6+6XtFremnoThpqj2jqPOJ/dy59j3mXzie\nva26Kud9xuDybcsAeGX4I8ouRadVc9vIbgHG/WhBOc/Pz8CDh6gILVabE6fTrXSaW56RregJLc/I\nDvm3UJ1m0awlO5R7z+QmN8JNJBA0kEVrMtl5sAC7002hxcb87/aQX2JVqpSLLDay/UTl+nZJQmpb\nddWtbrAiffh0ydnH529eH2AIALof38MbHz/MJ2/fzFNfPsfYDV8gHd+L1uUIuO5XP1fRuQcyeHzp\nyzzz1fPonXbeuOw+9rTqyojtKyiJMLG225B6jTGvmvRSnVZNTJS361r7k1l0y5F5dskzpPvFL3zG\nAGB9r4v5WaqUk3Y43cxftiegJuCF9zdgd7pxujyUWZ24PZUdyiBwZ1fdLs+3U5gwOv2Mqy4OB7Ez\nEPztqc4dVF91Sn/1TLfbE9CU3SdyF8qn7N/CsqlQedzcv3IWF+9eq5x7ZdQj/CydR1LJSXoc3UWP\no7tIP7qb/lmb6J/lldO2ag3sbSWxq3U6O9t0Z2+qhNlgxGTzGrnBmb8DUBwZQ8eTBxmz+Rtiy0v4\nufNgnOqqbrNwqK7wzOF0Kzurk7HeayIdVl7+bAov3PgCG1LSGXDwD+X6n8Y/BkcC6wl8GUj7Ktpe\n2qqZ4N3htourAZ9Utu/zmYowBoK/DdVN+tW5g8LxCfuyhPxdP13axpFX5J38U1tEVUl9XLQmU/FB\n+0sZh/JVNyZJJSeZ/15lg8KC6HjuvWWm0mM4L6Yla7u3ZG33iwCvhET6sd0VxmEXvbO3K4VnDo0W\nXYgsnyRLPmM2f6Mcn7f/V86d+Tt5MUnkxKaQE5fKibgUcuJSyIlN4URcCjZd6HqKvAo9pGolKTwe\ncmMCYwlTPp7CexeOp29FDcSCIePBZAJCN7gJVRjYpmU0R0+WAl4J8SKLLSAOFComVBNpyaYz1jXk\njzAGgr8N/pO+T8++oTnh/vLWvoyi3IJy2qWauHm4VGUHkF9iJX+3lcyjxaQlGykwW5tFg+i63xZx\n42//U47fGD6ZH9IvrvGeQmMCP0vnKS4WU3kJ3Y/tocfRXVy5eWnIe75PH0qb/KN0PbEPgB+6X0Rq\n0QlSik/QJ3sbfbK3VbknPzqeE3EpnIhN4bjPWFQYjsKo2CpuIoPDxoQ179D30GbmDJ1Q5Xl3rH9f\n+bw1rRcHdud6QxaeylqOmhoClVsrdwnFpXYWrcms9y7xz4QwBoK/Jb5J2e5wcfNwSTkfzj90X7Ob\n7FxLgBtBrVYpK/3C/TalnWafzi3YlVUQIHTniymEg16rxuFy18toRFstfPr2TQHnrpv4EaURxmru\nqB5zZAyb2/VRXDDebmUjsWt0XL3xSwAu2bUm4J4yfSRL+17OrtbdseoiSCnOJaXoBKlFOaQWnyCl\n6AStinLoelwm/dieKu8EiC8r5rFvX+V4XCvUHjdX//GF8rMp37ysfB794BcsfT0wddSi9RadBf/u\navpdut0epRmRjzNdZK4xqLMxkCRpDqCWZfmuGq7pB7wB9AGOAs/LsvxhvUcpEDQCvol+y748ZWLO\nzrXU+R96qGY3ndrEYne4AiZ434Tfo30Cpih92JO/D5XKW5k8fmRX5i/bW20Gkk6jRq32YHMEznAX\n7/yBB1bNUo7fu3A8X58zJuz3q90uWhcep8PJg7TPy6JDXhYdcw8SYzVzoGUHXrriMSXl81hCa+5f\nVbVp/RVbl3HFVm9Gz9H4Vuxsk862tLP5qt+YgC5oWpeDpJI8WhXleI1FsddgDDi4EYDz9v1a63jb\nFByjXBeBW6XmeHwrOudmogsKfPsTE6WjpKzqzwstNqX3hH8x2V+dOhkDSZKeBe4C3qvhmhbACuAj\n4DZgGDBPkqQcWZa/b8BYBYIG4Zv0Zy3ZrkzmwdXD9aFTm1ilfzF4DUxBiVVxSezKKiC9QwIlpTYc\nFQVniTERSj+D6vBUuDU+W3OADq1ikLOrNnqBiliDX/xT63TwyeyblcIrgBvv/qBGpdAIezntTh2u\nnPhPHuKsU4cxuAKroHNjWrK224UsPO9m7LrKjJlgRdNQtCk8TpvC41y2YzXPaZ9gQ8cBys+cGh05\n8a0C+yUDt61/n39u+pqDSe3okHeoYgxJ7E/pTErRCTqdPAiATaPn0e9eI9JhZWmfyxm95VsA7v1+\nNo+PexEIlPzQadU8NK43z7z/B6EotNjooov7y+8G/AnLGEiS1B6YB6QDh2u5/E6gSJblByqO90mS\n1Bd4BBDGQHDauXm4pAQB67PqGze0E3aHi0M5XrkCu8NFkcVGnNHA5LG9AJjw2jplJe/Bm73i73pI\nSzYy5rz2zFqyg5JSe43BY3O5ncLs8OQNOp/Yz/RPHlWOvzxnDPMvHF95gcdDQmkh7U8epEPFar/9\nyUOkFuWg9usH4FBryU5sS1bL9hxMak9WUnuyktph1UXQsiSP9GO7SS3KIaX4BG0KjtM/a2ON4zru\nFzzOTmzL5rNCaxoF41Mv7ZB3iCMJbXjpisc44i9r4fEQYzVz/a+fKjUFPkMA8P4Ftyqf/Y1BWksj\nackmJozuztylu5Vr9NrQfSv+DoS7MxgMZAPXAYtqufY84Megc+uAt+o0MoGgiWio/9c36U9ftJWd\nWQUU7rfhcLp5aFxllW3rFlEcOhHYsN5SHuiS8GWZPDjr52p1iHQadVhZRjqnnRt++5R/bvxKObcl\n7WxW97iEIbvX0S9rE92P7SapFmmHnNhkCowJqN1uWhcep2NeVq3vDocZl93H7tbd63zfgZYdAfip\nyz+YNWySIjznQ+1xBxgCgE1n9WHxgLHsbJOuFLvptWrcHg/uip3ZgeMlZOeaGdg9ha2Z+UpiQXr7\nhAYtFP7MhGUMZFn+GPgYQJKkWq6mDbA56NxxIEqSpARZlqtKDAoEf0L2+TVl9/+cnWuuYgggMHXU\nPz2xzFrVb61WQfd2Cew6VPs/F+n4Xh77dqoy0R+PS6VVUQ59srfx9geTw/syFaQW55Ja3DgtOL88\nZzQ2rYHiqFj2J9dvYt3TuhvX3/Ohkv7qQ+t0MHTPOm7++SPiyr0yEXmmFjw/5gkOtuzgFQN0eVAB\nHVrFEBOtZ1tmYGvLFz/cxMt3n8u4oZ0wGLTYbM4zsh1lc9EU2URRQHCkzLfHPb0C7QJBIxIdocNe\n4faJjtAp52sSOIuN0tO1XTwjBqYxd6m3Ylan0+AIztn3eI1KlawXPzdP1xyZ6zI+r/IOX/+A+lKm\nj+RwYhqHW5zF4RZpHGpxFtktzqIkwlRFVqJd3iFmffiAcvx9+lAGHPgDlcfD/wZdR7mhaq+BuuJv\nCCLt5QzfvpIrNy0NkMde2+1Cpo94UDl2ujyoVfDgtWezYLnMgeMlVZ5rr6hEfuja3jx6Uz/y8sxV\nrvk70RTGoBwINq2+49Labk5KOjNFnMJFjP/00pzjf+auc3n+/QwAnhw/UHm3Okii2j+nXatT8+Tt\ng5j60cYqvQ3Ubhfn7fuF7W17UhQdj8VipW1QNk/7k4eUlXB9sGoN5FTk9PuKv07Eev+bZ0rCpanb\nlHD5lu8Cjj8bcDU5sSnc/OsnjNy2giUDwlMJrY2Y8hKu2Pwto7Yuw2SzUKaLoCgylrjyYna36sob\nw++rco/bA29+sbPa6mKAvYeLOFpQzj3T1gMwZfwA+kgtG2XMfzaawhgcAVKDzrUCLLIs1/pX/Ge2\nzklJJjH+00hzj9+kV/PKhEqNHt+7772yR4D8gL9QWWGJlU07j2OzBe4EBu/7lce/fVU53pfcmbPy\nD2NwBsYSCmvpHQxQFBmrFG4FT/hFUXH1Eo0LhbHczJC965XjDR36kRPfiu96j+SqjV8xZvPXfNNn\nVEDWUV1JKsnjyk1fM3zHKgxOO8WRMXw4+AY2t+vLC58/Sak+iukVaqR6rZouaXEBFd/uWooznC43\nz7z7G75ykefnZ/DEzefUqjIarqJtc9EYi6CmMAY/A7cGnRsK/NIE7xIIzhj8J4gpt5yjTBCTx/bk\n2QV/4PZ4V6szPt9Ou1QTiTERuHJO8OHcW6s8q/2pLLIT05RMnoMt22O0Wrh9/Xygck21ouewgMn+\nRGxKo7hmwuHSnd8HGCtfz4LSCCPf9R7JtRsWM2znar6taIZTF9rkH2Hsxi8Zsmc9WreLPFMLvjzn\nSlb1vASnWsvLi54gymHltREPKrUOdqe7ivRHh1YxZB0vUTKE9Fo1Lrcb/5h8sPxQOCqjwRIm/npT\n1RmHM82ABNNgYyBJkg5IAApkWXbgTUF9VJKk2cAM4FK8WUjDG/ougeBMpchiCxCkg0qNo7RkE7HR\nlWmlFquDrftOMmLbCu5Z807Ac2YMm0RmcieOJrTGqfHGISLs5dz600JGbVuOS6Vmcf+rGJT5O4mW\nAt665J5GW+mDV8huyJ71WAzRbOzQL6AwzB+128WobcuV46wWZ7G9bS/l+Ou+VzB68zdc9ceXrOg1\nXPkutdE5Zx9X//EFgzIzUOMhO6ENS/pfxfquFxATG02bGAPnLp5D1xP7WNf1An5JvxCd2jvGUEV5\ncnZRQGe59PYJlNuc7A2q2VCrvGmnbo+HgjoWB0J4cufhSqKfLupjDIL3XYOBNcBFwI+yLJ+UJOky\nYCberKLDwM2yLK9HIPiLsmhNZo0Vxq2TohVj0O3UQSZ8PY20/CPKz2cMm8T6nhdz7c+f4NDoOJTU\nDoBe2du4b9WbJJfkkZ3YljeGTWZ/ahcSLAUM3bOOC/f+yPqg3sT1JcGcz4MrZyhidNkJbfii/z9Z\n3/WCKpN5/4MbSfZTE10a1LOgJCqWlb2GMWbzNwzdvY5VPS+t/sUeD2dnb+eaDYs5+4jXPSOndGbx\ngLFkdByAR6XmrGQjCTER2Nau5+qMxeTGJDH74gl0PSuBnPwy3G4Pbo8nZIqu0+VBr1XTp0tSFcVY\nX/U4ahUZO08oXejUKog3RVSrMhqsVfRn7GwWTJ2NgSzLQ4OO1wOaoHMbgEENG5pA8OfEJ1MNlV2u\nSsrsRFst3PzLx4zYtiKgwOv+m6ZzsGUHDLiV7KDux/ag9ri4bMdqXCo1iwZew6cDr8Wp9U7Knwy+\nnsH7f2PCmnfJ6DgAa1D+fV0ZvO9XJq1+G5PNwh/t+1ESaeLCvT/ywMpZ3PjLJ3zVbwyrelyqvOfy\nrZWB4+LIGNZ3vaDKM78450pGblvO2A1L+D59qNJlzIfa7WJQZgZX/7GEzrkHANhy1tl83v9qdrTt\nEWBcVCrYt+MQM5e/AcC0EQ+hi4/nSK6F4jKvAYiN0tOjfQJ7swurtBY1RemVlbh/17Krh3Rg8bqD\nVWpA4k0RTLnlHEWDKi3ZyM3DJcW1E1yrEo6Q3ZkudieE6gSCRsD3j9vX+MSndOnzP1+w90fuWDeP\n+DLvJGQxRGO0lbI3VWLCmneJtpUSbatMthu5fQUAx+JaMXXUwxxI7hjwvkh7GRq3i3J9ZLWunHCI\ntJdz19p3uWTXGmxaPW9fPIHlvS4DlYqPB9/AmE1fM3zHau5cN59xv3/Ot71Hsi2tl7J7AFje6zIc\nWn2VZxeYEvk+/WJGbF/J+fLPyg5G63IwZM96xv7xJW0Kj+FGxc+dB7Ok/1VK+0x/4o0GDuWYeeSH\nubQ05/G/QePY07obsSoo9wvEl9udPDSuN4++/WvALi0mShewwvfvWvbuN7uVmIFOo0algki9ltQW\nUQFuv/wSK3qdplrXTjiFjGe62J0wBgJBI+D7hz536a4qKaPRVgsPL3s9YDdgrJj4u+bIuFRqSg3R\nIZ9r0+kpjI4POKd32Hhk2evo3E5mDJ+MrZ7ZOtLxvTy8/HVSi3PJbNmBaSMf4mhCG+XneTFJvHfR\nHSwadC2jti7jii3fccPvi7jh90oRAoday7KzL6v2HUv6X8WwHau5dsNiNnQcwLAdq7ly09e0sOTj\nUGtZ1eMSvuj3T44ltK72Ga2Tojl7w0oulH9ib6rEp4OuBaDM5gyo9G7VInTg/M4rulfbd9g/eGyM\n1DFt0j9C/j/8OyCMgUDQRGQeLSY1MYpdZiPPXjmFGKsZW5QJlymWE24dpRHRlBqMlOsi6Je1iUmr\n3wbgQMsOvHXJPYzctpxLdq1h2ieP8uGEl1jj8ea/3/Lzh5yVn823Z49kU/tz6jwutdvFuN8/Y1zG\n56g8Hhb3v4qPB19fbZDXHBnDp+dex5f9rmTYjtX8c9PXJJm91bw/dj2fQmNCte/KjU1mc7u+9M/a\nyGdvXg9AuS6Cr/qO5qtzRpNf0cCmOuKNBop27OU/a+ZSpovgNb+m9h63m6MnKyu9D52w8PCbv+Bw\nBdYVvP7ZNvp3S1YyeMYN7VSlsRBAy/iqrjZTpI5SqwONRs15PVOU82d6ZlB9UHmao7NG+HhEnvvp\nQ4y/4fgmCf/JpkeHBHJOeTtt+dwVs5bswO320CnKxYjPZ9B/02rvSnnYzXzRfywnS53g8fDPjV9x\n608LcWj1TL/sfsr1kTz7xX/JTmjDQzdOq/OuIMpWyjNfPEu3HJk8UwumX3Y/O9vWsRWjx8OzXzxD\nn8PbeODGaVVcWD4Szae4ctPSgEY4H597Hd/1Hok5MiasVyVF63j0vYfoliMz/bL7lS5s9cEXLPa5\n9IKzv0xROmx2Fx6gQ2oMcSYD+7KLlMB/YkyEkmbqv3sY2D35tLt/kpJMDU4pEzsDgaABhFoh+txF\n+bu9E010hK5KrvrUiYPRL/sW478eRHMyF2efvjwz6C62R6ZAaYUfXKXiy/5eF8ojy6bz+LevUqqP\nwqHWMn3Eg/VyD7U7dZhuOTLrpfOZffGEejW5aVtwlD6Ht7GzdXpIQ9C64BhX/fEFF+1Zj87tJM+Y\nqOgm7UvpEpYhUKvg7E4tuOjbecp413YbEvLacJv/+BoagTetc8ot5/Dcgo2VKb/lDuUZWTklzLlx\nCI++XXsfhb8KwhgIBA0guJXmlFu8bhu7w0ViTARpycYqmSOq/HwM/3oY0zdf4NTqOfWv//BO5+Fs\n31MppBZvNFBcasPtgQ0dB/Cv617mP1+9QEtzHgvOuyXkJKzXqrA7a54Rd7fuXu9OZz6uqJCI9hWZ\n+eiYm8k1G5Zw7v7fUePhaHxrFve/ivXdLqBt/hFmfvQQ1274nM3t+9b6DrcHOhzayZDlH1CckMKG\ne/9DD1UEO7OqCvfNeWRIQI8K8PYr6NE+gTHntWd5RnZIt1Cc0cB/bu3HojWZGAxaftt+XOk34aO6\nZvZnemZQfRDGQCBoJPJLrMouwTcxdWoTq/iTiyw2tr36Dpd9/BpGcyF7UyVmDJtE6y79OJhTKaSm\nVkG7VBNl1kjkCjXUQ0ntePCmaUjH97KxQ7+Q7w9Op6yOpHYpRJQ569x5Dbx9kIfuXsuJ2GQ2dOwP\nHg+9juzg6g1LlP7G+5M7srj/WH7vNFDx72e17MAf7fvRP2sj6Ud3satNzW6VKFspo+Y/DcBXdzyN\nLcrE1QPTyMkvCxj3iEHe3gb+KrCJMREBFeC+xkP+Ozgfvp1cUpKJdRsOVZn4q2tmf6ZnBtUHYQwE\nggZQXTAyGFVeHo5b7uDqTWuxafR8OuxO/tf9MtxqDcf3BGauuD1eYxLvF5TsmhZHXlEE+xJbEKMC\ntUpFy/hIMo8V46pIiQmWVQCvYfGd9/nM/f3gdcWnEfRt75EMOPAHV/+xBOnEfsDbfH5x/7FsS+sV\nsip60cCr6Z+1kfHbvuTRtumKS8a3ivdx8FgJty97naSiXL4+/3oWO1OhYvc1eWxPXvxwkyIvUVDi\nrTEIXqkHB3TDmbzT2ycy55EhdfuF/IUQxkAgaABxRgNTbjkn5KoTYNxFHTF8uRjj44/QoqCA3a26\nMWP4JBL69qR/pI4Ne3KVyVqlggRThGJY1OrKCTXWaOBfN/QNcIe0SzUphsCH/+Sv16oZP7Iri9d5\nW0P6RNfq6wfXuJxcvsXbRGbMpqVKHODXToNY3P8q9qd2qfH+U117s79TbyR5I53O2c/+lM4ASG3j\nlA5xRRYb301+iYv2rmdvShc+HnBNgObB8ozsgE5k9iBF0tJyB88t2IharWLy2J6Yy+wBq/309ol1\n+s5/xayh6hDGQCBoIKFWnRNGp6PKzcU06TYMy7/FExnJsmvuY27rC3GrNSRr1UwYnc4fe3KVuU4F\nAYZlxMA0lmdkA5VGxtdqM/izj/7dkpUYRp8uSWzNzFeMy/KMbMYN7URqiygKzdYqOwl/QxKKczN/\nV3oIxJcV8X36UJb0+ydHE9uiAqIiNJRZq5eLjjfp+aD3VTyfuZVrMz7nhTFPVLlm5aIfuXP1HMp0\nEUwb+SCtW8WRGBup/A4+XCmHfHaw7g94M7b8+0zPWrKjziv/M11PqDERxkAgaGw8HgyLF2Gc8i/U\nRUXYB5+H+fU32byjHHfFxOLzcXduE6fEBaIjvXn+/hOOz9/tk0UIlk24aVhnPlq1X/nct0ulFn8o\nzZwPV8qKsme80UBirIEDx0oqRNq811RnFLQuJxZDND90v4iv+o3hlCmp8itDjYYAIPNYCbTtyd5U\niUEHNnBW3iEOJ7XjWF5F5bXTyZh5/yXaXsYbwydzIi6V+GIbT9zcT1mhBxtA/1iBoGEIYyAQNCLq\nnOMYH30Aw6oVeKKiMb88Deutt4NazbikSj+9b6U/YUy6ku9uLnPw4Uo5oAdvnNHAojWZAZkyPtql\nmhjaty1D+7YNOB9KMye/uJwNu3MDVCbVahXZuZYqypPd2yUgZxcpbTp1GjV6vZp13YewrvsQtBpV\n2MHqKqhUfDbwap766gWu2bCE10Y9TKHFxt2vreO1wrW0z9zOL9J5/NDdK4FWanUohiB45e+vATVu\naCe27MtTdgFqFSHdRHXlr5g1VB3CGAgEjYHHg+HTjzH+53HUJcXYz78Q8/RZuM9qp1wSyp0UZzTQ\nqU2skv+enWsJKYPtj16rxhSlZ8x57Wsdlu+dd7yypsqkH2/SVwl8q1QwfEBbbhvVjUVrMim1Ojh2\nspRSvz7NPTsksiurIMB3H4rYKD3xMfoq/aD/aN+Pg0ntOG/fL3wy+DqOx7emffZu0hbNwJ7Siuh5\nc9EtPYDD6cbudAcYSAgsHvMXjjNFVX6feFOEIkHRkKDwXzFrqDrqr3AlEAgAUB87Suz1Y4m5fyK4\n3Zhfm0Hx4qUBhqAmxg3txKAeKSTGROD0m2A37MnllY83Y3e46NEhgcSYCOKNBuxON/klViWeUBtF\nFlsVt0+P9glet00QHg8sWC4rk2DOqTIKLbaA5jBD+7bm8n+cVet7tVo1yQnRmCKDZC5UKj4beA0a\nj5urN3xBpK2MR5a/jsrjYeqw+/l2bwn+uUjZuRbGDe3EwO7JDOyezJRbzmHC6PSAYG6RxUZqYhR6\nrZp4o6Feu4C/O2JnIBDUF4+HiI8+IPrpKagtZuwXXYx52kzcbdrWfq8fcUYDep2myird40GJJwzs\nnszUiYPrJKLmL43hj0atYs+hqsVbPtxuD3OX7gIIME7g7SY2bdG2au/17/fsq/jVaaquOX/rNIij\n8a25aM86UopPkFKcy2cDxrI5pRv2oO+XlmysdYW+aE2mUpDWJy2uWmE6QfUIYyAQ1AN19mFMD92H\n/se1uGNiMb/xFtbrb2rUrmOhCPZhF1lszP9uD/uOFBEdoeP+a3opE2EoPztQJR3Vn+gILYmxBuU+\n/y5hNZEUH0mbFtFKxa+//97h12PyrGQjJaUOCi02Pht4NQ+tmEHPo7vYl9KZlcPH06VljDKpB2sJ\nCZqWsIyBJElq4AXg/wATsAK4V5blk9VcPxR4CUgHcoB3ZFme2igjFghOJ243ER/MJ/rZp1CXWrBd\nOhzLazNwp7Zq0GNvH90Dm82J3eHCXGbn4PES1GoVHVvFotOq2ZddxKNv/8rksT0DVshzl+5SJk+7\nxVZtz97EmIiw0kkdTjfZuZU+/nACxRNGd+fyCzuzP+uUkr2UFBfBsQpxPn8O51p4eNzZ/LzjBBm2\nAcr5N4fdy7/HD+TDlbIi4+HfTAZqzvn/OwV6m4pwdwb/BW4GbgIKgNnAYqBKeyNJkjoC3wAvAuOA\nvhZgTm0AAB3kSURBVMBCSZIssizPboxBCwSnA/WhLEwPTkL/y0+4Y+MomTUH27XXN8puICEmolo3\nyKNv/6pUDFc32YcieII8lmdh1pIdiirnqWJrFdeU3enGFKkLCA7XVH8QbzSwNTOfPUeKMVtsIbOe\nglmwXGbqxMHcveMoW9LOxuC0cySxbUDWlL+Mh4+acv7/ToHepqJWY1DR8P4+YJIsy2sqzl0HZEmS\nNEiW5d+DbrkMKJNl+YWK40OSJI0DhuM1IgLBnwu3m4j572B8/hlUZWXYLhuFZerruJNTar21qRk3\ntBOl5Q7FTeQfOA2eIOOMhoDMGn9lVX/UahWx0XrMZXZMkXquu6QT7y/bWyV7SK9V07pltDJBx9eh\nOrfIYqNtWiJPXf1fVCq464rubM3MD3mdbzdQ6ldjEVx5/HeqFG4qwtkZ9AaMgNLQXpblw5IkHQLO\nB4KNQR6QUGEwFuF1FV0AvNkI4xUImhXNwUyMD0xC//uvuBMSME+fhe2fVzd5bMCf6pQzwTvBPzSu\nd1jP8fVj9j3Hf+ew93Ch0kxerVLRJS2OjN25FJfZ2ZqZT58uScqk7+/Lr64RvEatqhKb0GpUmCL1\n3DpCCugl0LtTC7Zm5mN3uOjTuQV6nUYZm/9uoCZj83eqFG4qwjEGvj54x4LOHwdCpU0sAeYDHwMf\nAhpgkd9OQSA483G5iHxnNtEvP4eqvBzb5WMwvzwNT8uWtd/byFSnnFlXfP2YfZ+nThysTJrTF22l\nuCL20DopOmAVXlru4LZR3ZRj/4naN4GbjAbMFhuF+73urJgoveLa6tO5haI9BBU7Ej/3lH9tRU2N\nYvy1mkTlceMTjjGIAtyyLAfXmtuAiBDXxwHtgJeBz4CewAxJkp6RZfmZ+g9VIGgeNPv3Ybp/IrqN\nG3AnJlIyaw720f883cOqlsZwkUT71QJER+rYl12kHB/LK63icvJPce3TubJ1pW9lH0pXKRS+YHF1\nqq8jBqax93ChV2Ooot4iOkJXVRBQBJAbTDjGoBxQS5KklmXZ32loAEpDXP8q4JBleUrF8baKuMNs\nSZJmyLJcWNPLkpL+3PnBYvynlwaN3+mE6dPhqafAZoNx41DPmkVsUlLt9zYC9R37gpWyMjEbDFoe\nvSl0v4On7hjE8+9nAPDk+IEB75t4TW8MS3cC3symR2b+qPxMq1NXGZvBUDl1HD1VSl5F4PeCPq2V\n95/TI3SGVfC7AOb5HSfEVK4xF6yUFfeVudxBdq4Zm93N1v2nmDJ+AH0k704tKcnEk7fX3E+5Nv7s\nf/sNJRxjcKTiv6kEuopaUdV1BDAQ+CLoXAagB9KAGo3B6e5h2xDOhB68DeHvPH7N3j2Y7r8H3ZbN\nuJNaYn71deyjKjp5NcPvpCFjt9mcAZ+re45Jr+aVCecqx8HX3TpcAsBlc3DvlT2U+MK9V/aocu2V\n/2invNfucJFXWF7r+0O9Kz/fErCrcdkc5OVVuqj8vxtASWnlz56fn9Fo/Qf+Cn/7DSUcY7ANsAAX\nAp8ASJLUDq8r6McQ1x8FegWd6wm4gAP1HKdA0DQ4HES9NYOo115GZbdjvXocludfxpNQN93700lT\nuEhqi1P4u42KLDZMRgM2m7PO768p8FtksWF3uNBp1AGFa4KmoVZjIMuyXZKkt4HXJEnKx5st9Baw\nVpblDRUuoASgQJZlBzAD+EaSpCl4jUc6MA14S5Zly/+3d+dhUlV3Gse/1ftW2KCEiKAkcXKMy4So\nCUpEFA0uEzU6Ki4wZh6JrAIiCkSigCKKIDuohCRikhnBLTCJKAqigkqCiomJBxdQlpnYQkOv9FJd\n88e9rU3Ry63uqrpV3e/neXjavn1u1/uUfetX595zf9X4o4gkXvr7f3OuDbz3LqFuX6ds9nyqL7rE\n71hRi9Ua+9ZeeygsyOaOwWey9W97mbFiK/DVB+m09BjNaXjfQf0Kpt4nHs2v//TBl48hseP1prMp\n7tgngEzgeWC0+7O+wHrgfOBVa+3zxpir3H0mAv8HPIJzR7KI/6qryZs/h7x5swnU1HDouhspm34/\n4cLOfieLq5Ze7Bu+S/9o98HDPkfYi8ZWK0WKnAl4ndWc2OOoLwten5P9v7+jPfJUDNyVRHe4/yJ/\nthFn+WjDbauB1bEIKBJLGX/dRnDMSDLe/yuh7sdRNmc+1RcM9DtWQkSzFn9fySGeXP9R3NfrNzer\n0QqhxFKjOukYqqrImzuLvAVzCdTWUjn4Jsqn3ke401F+J0sagwacyEe7Dza5zLMlzd0c1/AxGvvv\nxqjFRGKpGEi7l/Hu2wTHjCDjg38Q6tGT0ocXUnPeAL9jJVxLL8SFBdmHfQZztO/Gvdwcpxf45KVi\nIO3XoUPkz36A3MXzCYRCVP70Zsrvnk64oGOuJ/fyQqwX645LxUDapYy/bCE4bhQZ2y2h43tROm8R\nNecc0WQ3pUX2GmrrB7qo2VvHpo+9lPalspL8e+6i8McDydhuqRg6jP2vbG53hQC+Wr2zr+TQl0Wh\nLeovML/193822oDuQFkVj65+n0dXv88Bt+9Qc6IdL/7SzEDaj02b6HzTT8n4+CNCvb5B6fwl1Jz9\nQ79TtRvRdgZVJ9HUomIgqa+8nPyZ02HZI6QDFcNGUT75F5CX53eyuPKyeicaWsrZsakYSErL3Pw6\nwXGjSN+5A4zhwJxF1P6gj9+xEiJWra3rtXTxONpioeKSWlQMJDWVlVFw3z3k/moZ4bQ0KkaPI2/W\n/dSW1ba8r7RKtCuNtDIptagYSMrJfPUVguNvJf2zT6k1J1E6bzG1Z3yfvNxcKEvdzpMiftJqIkkZ\ngdISCiaMo/Dqy0nbs5vycRMofuk1as/4vt/RRFKeZgaSEjI3vOzMBvbspvY7p1C6YAm13/2e37FE\n2g3NDCSpBQ4eoGDcKAoHXUnaP/+P8tsnUrxuowqBSIxpZiBJK2vdWgomjCP9f/dSc+q/Ujp/CaHT\nIj83SURiQcVAkk7gQDEFUyaRs/K/CGdmUj7xLirGjIfMzJZ3FpFWUTGQpJK19k8UTBhL+uf/pOa7\n33NmAydreaJIvHkqBsaYNGAGcBMQBNYCo6y1nzcx/jicj78cCFQCTwG3W2tb1yhd2r3A/n0U/PxO\ncp5ZRTgri7IpU6kcOQYy9H5FJBG8HmnTgCHAYGA/sBTnBf6I7l/GmCzgJWAPcDZwDLACCAFj2h5Z\n2pusNX8gOHE8aV8UUXPGmZTOW0LInOR3LJEOpcVi4H7g/RhgtLV2vbvtOmCHMeYsa+2bEbvcCHQD\n+lhrS9zxdwMjYppcUl6gqIiCyRPIWf0s4exsyu65j8rhoyA9veWdRSSmvCwt7Q0UABvrN1hrPwV2\nAv0aGT8QWFdfCNzxj1trz2pbVGk3wmGyn3uaLuf+gJzVz1Lz/T4Ub9hM5agxKgQiPvFymqiH+3VP\nxPa9QM9Gxn8beNkYMx3ntFIYeAaYYq1VU/MOLvD55wQnjif7j6sJ5+ZSdu9MKocOVxEQ8ZmXYpAH\n1FlrQxHbq4CcRsZ3AoYCfwKuBo4DFgNdgZ+2OqmktnCY7KdXUnDXnaQVF1N99g8pnbuIum9+y+9k\nIoK3YlAJpBlj0qy1dQ22ZwPljYyvAfYBQ6y1YeBt96LySmPMbdba4uYerGvX1P58WuVvxN69MHw4\nrFkD+fmwcCFZI0dydFrsb4BP5ec/lbOD8qc6L8Vgl/v1WA4/VdSdI08d4W6rdAtBvb8DAaAX0Gwx\nKCpK3a6TXbsGlb+hcJjsJ39PwS8mk3bwANX9+lP68ELqTugF+xp7H9E2qfz8p3J2UH6/xaKQeXlr\ntg0oA/rXbzDG9MJ5YX+1kfGvAb2NMQ1PAp8G1OJcdJYOIG3vHjrdcDWdxoyA2lpKH5rHwadWO4VA\nRJJOizMDa221MWYJMNsYsw8owrkGsMFau8VdetoF2G+trQEeAUYDK9yLyD2BWcDjLZ0iknYgHCbn\ndyvIv+cu0kpLqO5/vjMb6Hm838lEpBleT9pOAX4HPAG8DOwArnF/1hdnZdHZAO5dyefiFIitwG+B\nVcDImKWWpJS26zOOuvYnBMffCkDp3EUcXPmcCoFICvB0B7K7kugO91/kzzYC6RHbPgAuiUVASQF1\ndeSs+DX5035BWnkZVRf8iLI5C6jrfpzfyUTEIzV+kTZJ+3QnwdtGk/X6q9QdVUjJgqVUDboBAgG/\no4lIFFQMpHXq6sj59TIK7p1KoKKcqosvpWzWXOq+fqzfyUSkFVQMJGppn3zszAbe2ERd586UzplP\n1VXXaDYgksJUDMS7UIjcZUvJn3kvgcpKqv7tckofmEO4Wze/k4lIG6kYiCfpH31IcOxIMv/8FnVH\nH03pgqVUXX6lZgMi7YSKgTQvFCJ36SLyZ80gcOgQh664irKZswkfc4zfyUQkhlQMpEnp9gOCY0eQ\n+fZW6o7pSsmSX1L948v9jiUicRD7TmGS+mpryZ0/h84XnEPm21s5dNU17H99iwqBSDummYEcJv3v\n7zvXBra9Q6jb1yl7aB7VF1/qdywRiTPNDMRRUwP33kvnH51L5rZ3ODToBopfe0uFQKSD0MxASP/r\newTHjoS/vUfdsd0pmzOf6gsv8juWiCSQZgYdWXU1eQ/cR+eLziPzb+/BzTc7swEVApEORzODDipj\n2zsEx4wk4x/vE+rRk9I5Cyi89ieEU/gDPkSk9VQMOpqqKvJnP0DuonkEQiEqb7qZ8runEQ528juZ\niPhIxaADydj6Z4JjR5Kx3RI6/gRK5y6ipl//lncUkXZPxaAjqKwkf9b95C5dSKCujsqbb6HsrqlQ\nUOB3MhFJEp6KgTEmDZgB3AQEgbXAKPdTzVra93+APGvtgLYEldbJ2PIWwbEjyPj4I0K9vkHpvMXU\n9D3H71gikmS8riaaBgwBBgP9gB7AUy3tZIwZBmihuh8qKsj/xWQKLxtI+icfUzFsJPs3bFYhEJFG\ntTgzcD/wfgww2lq73t12HbDDGHOWtfbNJvY7EWc2sTmGecWDzDc2UTBuFBk7PqH2WydSOm8JtX3O\n8juWiCQxLzOD3kABsLF+g7X2U2AnzizhCO5ppceBB4B/tDmleFNWRsHkCRRecQnpn+6kYuQYitdv\nUiEQkRZ5KQY93K97IrbvBXo2sc/PgTpr7ezWBpPoZL62kS7n9SV3+WPUfttw4I/rKJ96H+Tm+h1N\nRFKAlwvIeTgv7KGI7VVATuRgY8wZwG3AmW2PJ8357J+lLPv9W1y7bjkD/vJHwunpVIwZT/mESZBz\nxP8aEZEmeSkGlUCaMSbNWlvXYHs2UN5woDEmG1gBTLHW7mhNoK5dg63ZLWkkMv/yyY8xbdWDdCsp\nYvfXetHjj6vIO/NM8trwO/X8+yeVs4PypzovxWCX+/VYDj9V1J0jTx31AU4CHjTGzHK3ZeMUkxLg\nZGvt7uYerCiF2yF07RpMTP7qavJn3c+dy+cSDgT47z7X8NKPhjDzBANtePyE5Y+TVM6fytlB+f0W\ni0LmpRhsA8qA/sDvAYwxvYBewKsRY98C/iVi20zgeOAGnOsM0gbpH24nOGIome+9S1WPE5h10Vg+\n7nkyt/77aX5HE5EU1mIxsNZWG2OWALONMfuAImAxsMFau8VdetoF2G+trQI+abi/OyOobO1pI3GF\nw+Q8/isK7vk5gcpKKq8fTPmMBxle0LGntiISG17bUUxxxz4BZALPA6Pdn/UF1gPnc+RMQWIgUFRE\n8LZRZL+4lrrCQkoWPUb1ZVf4HUtE2hFPxcBdSXSH+y/yZxuB9Gb2/Vmr0wlZ69YSHDuKtC+KqD73\nfEoXLqXu2O5+xxKRdkYfbpOsKioouPM2jrrxWgIlBymbfj8HVz6rQiAicaGupUko4713CY4YSsaH\n26n9zsmULPkloVNO9TuWiLRjmhkkk1CI3AVzKbzkAjI+3E7FsJEUv/CKCoGIxJ1mBkkibfcugqOH\nkbX5dULdvk7pgqXUnH+B37FEpIPQzCAJZD+zis7n9SVr8+tUXXoZxa+8oUIgIgmlmYGPAgcPUDBp\nAjlPryScl0/pvMUcun4wBAJ+RxORDkbFwCeZb2wiOOoW0nfvouaMMylZvIy6b37L71gi0kHpNFGi\nVVeTP2MaR/3kUtL27qF8wiQOrH5BhUBEfKWZQQKlf7id4MifkbntHUIn9KJkyTJqv9/H71giIpoZ\nJEQ4TM5vltP5wn5kbnuHyusHU7xhkwqBiCQNzQzi7Mi+Qo9SfdlP/I4lInIYFYM4ynrpBYJjRjp9\nhfqdR+miR9ROQkSSkk4TxUNFBQUTx3PUDdd81Vdo1XMqBCKStDQziDH1FRKRVKSZQayEQvDgg1/1\nFbplhPoKiUjK0MwgBur7CrH5deq+1o3ShY+onYSIpBRPxcAYkwbMAG4CgsBaYJS19vMmxg8CJuF8\nHvJeYDnwkLW2Lhahk0n2M6souHM8aSUH4corKb7/YcJHH+13LBGRqHg9TTQNGAIMBvoBPYCnGhto\njLkE+C3wGHAaTlGYCExua9hkEig5SHDEUDoNv5lAbS2lcxfB00+rEIhISmpxZuB+4P0YYLS1dr27\n7TpghzHmLGvtmxG7DANWWWuXut/vMMacDPwnzuwi5TXVVyioBnMikqK8zAx6AwXAxvoN1tpPgZ04\ns4RI9wLTI7aFgc6ti5hEIvsK3T5RfYVEpF3wcs2gh/t1T8T2vUDPyMHW2q0NvzfGdAKGA8+3JmCy\nSP/oQ4Ijhjp9hY53+wr9QO0kRKR98DIzyAPqrLWhiO1VQE5zOxpjcoHn3HGpec2gvq/QBeeQue0d\nDl13I8UbXlchEJF2xcvMoBJIM8akRawGygbKm9rJGHM0sAY4CbjQWrvLS6CuXYNehiXG55/D0KGw\nZg107gwrVpBz9dXNVsCkyt8Kyu+fVM4Oyp/qvBSD+hfxYzn8VFF3jjx1BIAxphfwIpAP9LPWvu81\nUFFRqdehcXVEX6GFS6nrfhw0k69r12DS5G8N5fdPKmcH5fdbLAqZl9NE24AyoH/9BvfFvhfwauRg\nY0xXYAPOReOzoykESSGyr9A0t69Q9+P8TiYiEjctzgystdXGmCXAbGPMPqAIWAxssNZucZeedgH2\nW2trgCXu9wOAKmNMN/dXhZu6SS1ZZPx1m9NXaLul9qTvOH2FTj3N71giInHntR3FFHfsE0Amzsqg\n0e7P+gLrgfONMVuAK4EAsKXB/gGgFsiKQebYC4XIXbKQ/AfuJVBTQ8UtIyi/ayrk5vqdTEQkITwV\nA3cl0R3uv8ifbQTSo/2dySJt9y6Ctw4na9NrhL7WjdIFS6kZcKHfsUREEqpDdy3NfvYpOp/Xl6xN\nr1F1yY8p3vimCoGIdEgp9S4+VgIlBymYeDs5T68knJdP6dxFHLphCKidhIh0UB2uGGS+udnpK7Tr\nM2pOP4PSJcsIffNEv2OJiPiq45wmqqkh7/7pTl+hPbspH38nB9a8qEIgIkIHmRmkf/QhwZFDyXxX\nfYVERBrTvmcG4TA5j/+Kzhf2I/Pddzg06Ab1FRIRaUS7nRkEvviC4G2jyH7heeoKCylZsJTqy6/0\nO5aISFJql8Ug66UXCI4dRVrR51T360/pwkfUTkJEpBnt6zRRRQUFk253+godPEDZ1BkcXPUHFQIR\nkRa0m5nBYX2FzEmULF2uvkIiIh6l/swgFCJ34TwKLx5AxnZLxc+GU/ziRhUCEZEopPTMIG3PboKj\nh6mvkIhIG6XszEB9hUREYiflZgaBkoMUTJpAzlNPEs7Lo/ThhRy68T/UV0hEpA1Sqhgc1lfoe6dT\nuvSXaichIhIDnoqBMSYNmAHcBASBtcCopj65zBhzJjAP+B6wG7jPWvtEq1PW1JD30EzyFjwMQPn4\nO6m4fSJkZrb6V4qIyFe8XjOYBgwBBgP9gB7AU40NNMYcg1Ms/oJTDBYCy40xrTqhn/7xhxT+24Xk\nz5tNXY/jOfCHtVRMmqJCICISQy3ODNzPOB4DjLbWrne3XQfsMMacZa19M2KXnwEHrLXj3O+3G2NO\nByYAL3lOFg6T88RvKLh7MoGKCg5dez1lMx8iHOzk+VeIiIg3XmYGvYECYGP9Bmvtp8BOnFlCpHOA\nVyO2vQL80GuowBdf0Omm6wlOGEs4K4uSZb+hdNGjKgQiInHipRj0cL/uidi+F+jZxPjGxuYZY7q0\n9GBZL79Il/5nkb32T1T360/xK29QdcVVHmKKiEhrebmAnAfUWWtDEdurgJwmxh9qZCxNjP/S2wMH\ncfq6lYSzsiibOoPK4aMgLWVvhRARSRleXmkrgTR3RVFD2UB5E+OzGxlLE+O/dPq6lRR1/wbFazdQ\nOfJWFQIRkQTxMjPY5X49lsNP/3TnyNNB9eOPjdjWHSiz1h5s7oEuG/8cwH+tOb/vDR5yJaWuXYN+\nR2gT5fdPKmcH5U91XorBNqAM6A/8HsAY0wvoxZEXigFeB34asW0AsKmlB1oz5wrdRiwi4oNAOBxu\ncZAxZibODWf/CRQBi4EKa+0F7tLTLsB+a22NMeZrwAfAk8B84EfAQ8BF1tqNjT6AiIj4yutJ+SnA\n74AngJeBHcA17s/64qwWOhvAvSv5Ypwbzt4GRgJDVAhERJKXp5mBiIi0b1quIyIiKgYiIpLAFta+\ndz5to1bkHwRMAv4F55rKcuAha21dYhIfkSeq/BH7/g+QZ60dEN+UzWaI9vk/DmcBw0Cce1+eAm63\n1kbeEJkQrcg/AJgJnAL8L/CYtfahBMVtkjHmESDNWntLM2OS6thtyGP+pDp2G/KSP2K852M3kTMD\n3zqfxkg0+S8Bfgs8BpyG84c1EZickKSN85y/IWPMMODS+EbzJJrnPwunKWIhzsKGa4EfA7MSkrRx\n0eT/FrAGWA2civO3c48xZkRiojbOGDMdaPZFKEmPXcBz/mQ8dgFv+SPGR3XsJmRm4Fvn0xhpRf5h\nwCpr7VL3+x3GmJNxlubOSFTueq3IX7/fiTh5NycsbOM5os1/I9AN6GOtLXHH3w348mLaivwX4yzd\nrv9b2em+W70IWEqCGWO+gfPu+BTg0xaGJ9WxC1HnT6pjF6LOX79P1MduomYGCe98GmPR5r8XmB6x\nLQx0jlO+lkSbv/60xuPAA8A/4h+xWdHmHwisqy8E7vjHrbVnxTlnU6LNXwR0McZcZ4wJGGNOBc4F\n/pyArI3pC3yG8055Zwtjk+3YhejyJ9uxC9Hlb/Wxm6hrBq3pfPp2I2PzjDFdrLX7Y5yvJVHlt9Zu\nbfi9MaYTMBx4Pi7pWhbt8w/wc5wGhbONMcvilsybaPN/G3jZnVYPxjmYnwGmWGurGhkfb9Hmfxr4\nFV/d25MOPNlgppBQ1trfuVkwxrQ0PNmO3ajyJ+GxG+3zD608dhM1M0hY59M4iTb/l4wxucBz7ji/\nzjtGld8YcwZwG/AfCcjmRbTPfydgKPBN4GpgHDAIeDSeIZsRbf5CnHYvDwBn4vx/GGiMmRrHjLGS\nbMduqyXJsRuVthy7iSoGCet8GifR5gfAGHM0zh3bvXHacexqamycec5vjMkGVuC8i96RoHwtifb5\nrwH24dz5/ra1dg3OATLEGOPHdD/a/LOAGmvtXdbabdba3+Kcc5/kU/5oJNux2ypJdOx61tZjN1Gn\niRLW+TROos1f38zvRSAf6GetfT+eAVsQTf4+wEnAg8aY+tU32TgvZiXAydba3fEM24hon/89QKW1\ntuHt9X8HAjjvuIvjkLE50ebvg3Naq6G3gCzgeBKfPxrJduxGLcmO3Wi06dhN1MygYedTwFPn03Mj\ntnnqfBonUeU3xnQFNuCcqz47Cf6Yosn/Fs766t7Ad91/z+JcvPwuzvnfRIv27+c1oLcxJr3BttOA\nWjxcgIuDaPPvBv41YttpQAj4OC4JYyfZjt2oJOGxG402HbsJmRlYa6uNMUuA2caYfXzV+XSDtXZL\nZOdTnGVUdxhjlvJV59PrcJbWJVwr8i9xvx8AVBljurm/KuzlJi+f81cBnzTc331XUenXaaNWPP+P\nAKOBFe5F5J44p14et9Ym/F11K/LPB9YYY+7CaRt/CjAHWGytLUt0/uYk+7HbkmQ/dlsSy2M3kTed\npXrnU0/5jTE5wJU4Swm3uNv34txFmujTKw15fv6TVLR/P+fiHCRbcW4iWoXzd+SXaPI/D1wFXIEz\nq3gYp8DdntjIjYrsbJkKx25DTeZP4mO3oWaf/7ZQ11IREVGjOhERUTEQERFUDEREBBUDERFBxUBE\nRFAxEBERVAxERAQVAxERQcVARESA/wfvADcRcz5SkgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10e254390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(samples_MVN[:,0], samples_MVN[:,2], \".\")\n",
    "plt.plot(unique_samples[0:50,0], unique_samples[0:50,2], 'r')\n",
    "pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pymc3 as pm\n",
    "import numpy.random as rng\n",
    "import scipy.stats as stats\n",
    "from sklearn.preprocessing import StandardScaler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " [-----------------100%-----------------] 5000 of 5000 complete in 0.6 sec"
     ]
    }
   ],
   "source": [
    "niter = 5000\n",
    "with pm.Model() as test_context:\n",
    "    mu = pm.MvNormal('mu', 0, tau=cov, shape= 10)\n",
    "    y = pm.MvNormal('y', mu=mu, tau= cov, observed=X_MVN)\n",
    "    step = pm.Metropolis()\n",
    "    trace = pm.sample(niter,step=step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "metropolis_samples = trace['mu']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "unique_samples_met = np.array([metropolis_samples[0,:]])\n",
    "# subset just accepted theta's\n",
    "for i in range(1, niter):\n",
    "    if all(metropolis_samples[i-1,:] != metropolis_samples[i,:]):\n",
    "        unique_samples_met = np.vstack([unique_samples_met, metropolis_samples[i,:]])   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEJCAYAAAB2T0usAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXdcVeUfx993socgooIswaMgDpyZpZJZVlrZsGFlplma\npZUt29ufLUfDLCvTysyGLUe591YEPeACJyLKvlzu+v1xuZd7GXJBlva8Xy9ewDnPfc73XuV8zvN9\nvkNhsVgQCAQCwX8bZWMbIBAIBILGR4iBQCAQCIQYCAQCgUCIgUAgEAgQYiAQCAQChBgIBAKBAFC7\nMkiSpBbANOBawAPYAjwly3JyFeN/BG4HLICi9PA/siwPumiLBQKBQFDnVLsykCRJAfwKRANDgCuA\nXOBfSZKaVfGyjsAzQCugZenXHXVhsEAgEAjqHldWBp2BXkAHWZZTASRJug84B9wIzHccLEmSFqtw\nbJNl+UzdmisQCASC+sCVPYMM4CabEJRiLv1e2cqgPaAC9l+kbQKBQCBoIKpdGciyfA74u9zhJwB3\nYHklL+kIGIDXJUkaDOiARcCbsizrL85cgUAgENQHLm0gOyJJ0lDgbeB9WZblSobElX5PAWYC8cCH\nQCjwYC3tFAgEAkE9oqhJoTpJkkYCnwPfybI88gLj/GVZznH4/U7ge6C5LMvna22tQCAQCOoFl1cG\nkiRNAd4AZsiyPPFCYx2FoJSk0u9tgCrFwGKxWBQKRVWnBQKBQFA5F33jdDXP4BngdeBFWZbfrmbs\nQkAjy/Iwh8M9AD1w8EKvVSgUZGXlu2JSgxEU5CNscpGmaJewyTWETa7TFO0KCvK56DmqFQNJkjoB\nbwFzgS8lSQp2OJ2PdbM4ADgny7IB+An4XpKkScBvQALWhLVpsiwXXbTFAoFAIKhzXAktHV46bhRw\nstzXRKBP6c9XAMiyvAgYWfqVhFUIPpRl+ZW6NV0gEAgEdYUroaVTgCnVDFOVe818yiWjCQQCgaDp\nIgrVCQQCgUCIgUAgEAiEGAgEAoEAIQYCgUAgQIiBQCAQCBBiIBAIBAKEGAgEAoEAIQYCgUAgQIiB\nQCAQCBBiIBAIBAKEGAgEAoEAIQYCgUAgQIiBQCAQCBBiIBAIBAKEGAgEAoEAIQYCgUAgwMUeyP8F\nko9kM3NxEgATbosnLjKwkS0SCASChkOsDEqZuTgJt/wcrtz7L6s+mI8y8zRYLI1tlkAgEDQIYmXg\nwJBdf3L35oXWXxa+hNnXDzq0xyeiLcYYCVNMO0zt2mEKjwS1+OgEAsHlg7ijlTLhtnhyV+kBOD34\nVpopTajSZJQ7duC+ZYvTWItGgykyClOMhDGmnVUkYtphjG4H3t6NYb5AIBBcFEIMSomLDMSrRxhs\nA7eJj5PXtRsAQf7unNuehCotFVWajLr0uyotDXWqjFu5eUytQ6zCENMOU+lqwhgjYWnRAhSKhn9j\nAoFA4AJCDByp7Gat0WCKjsEUHQODbyw7brGgPJOJKlVGlZaKulQgVGky2jWr0K5Z5TSN2dcPU0yM\ndTUR3Q5TO8n6u3A5CQSCJoC4C1WGKxvHCgXm4JaYg1tiuKqf86mC/NKVRGrpSiIV1cFU1Ht2o9mx\n3flSGg2mqLalLqeYstVE2xjhchIIBA2GEANHbCuDi4wisnj7YOzaDWPXbugdTxgMqNKPWlcTB1NR\nl35Xpaailg9UdDmFhJa5nEpXE8bodtBciIRAIKhbhBg4UkdiUCWOLidHLBaUmacrdzmtXol29Urn\n8f7++LeNwdhOsouEKSYGU1iEcDkJBIJaIe4cjtS3GFzguuaWrTC3bIXh6v7Op/LzUB1MQ5Uqoy79\n7nbkIOo9u9Ds2OY01qLVVu1y8vJqwDckEAguNYQYNHEsPr4VXE5BQT6cPXkO1dEjFV1OaWmoD+yv\n6HIKbYMp2nk1YYxuhyUoSEQ5CQQCIQZO2FcGjWuGS2g09vwGJywWlKdPlYmEbQM7LZXKXE5mf39M\n0e0cRMKaL2EOjwCVqqHejUAgaGSEGDhgaSw3UV2iUGBu1Rpzq9YY+g1wPpWXW8HlpDqYinrXDjTb\ntzqNtWi1mNpGW4XClljXTsIYFS1cTgLBZYgQA0ds3pJLWQwugMXXD2NCd4wJ3Z2jnEpKrC4n++a1\nQ2Ld/pSKLqc2YZW6nESUk0Bw6SLEwJHSlYHikvAT1SFarTUiqZ1ECUPKjlssKE+dLJd9XepyWvUv\n2lX/Os8TEOAc5RQTgzFGwhwWLlxOAkETR4iBoGoUCsytQzC3DqnocsrNsbqcHBLr3A6nod65Hc22\ncrWc3NwwRUU713GKkTC1jQZPz4Z8RwKBoAqEGDhxGewZNBAWP3+M3Xpg7NbDOcrpRDaqI4fLuZys\ngqHen+w8h0KB2eZyslWFtdVyCgwUUU4CQQPikhhIktQCmAZcC3gAW4CnZFlOrmJ8d+AjoCtwHHhT\nluVv68Ti+uRy2EBubLRaTFJ7TFJ7ShyPWywoT55wFonSTWztyn/QrvzHaRpzs2YOVWGlMpdTmzDh\nchII6oFqxUCSJAXwK9aAyyFAIfAa8K8kSR1kWT5fbnxzYCkwHxgFDAK+lCTplCzLzn/xTQ0hBvWH\nQoE5JBRzSCiG/onOp3Jz7OJQVqJDRr19K5qtm53GWtzdK7qcottZs7o9PBryHQkElxWurAw6A72A\nDrIspwJIknQfcA64EetN35ExQI4syxNLf0+VJCkBeBpo2mIgaBQsfv4Yu/fE2L2nc5STXl+5y+lg\nKuqUfc5zlLqcHOs4mWLawRXdoEI8lEAgKI8rYpAB3GQTglLMpd+bVTK+L7C23LHVwMc1tq6hESuD\npoWbG6b2HTC17+DscjKbrVFOqTLq0kJ/tixst39XwL8rnKYJDAio3OUU2ka4nASCUqoVA1mWzwF/\nlzv8BOAOLK/kJaHAznLHTgKekiQFlM7XNBFicGmgVJa5nAZc43RKkXPeyeXkmXEY875k1Nu2oNmy\nyWms3eXUzrkqrKlttHA5Cf5z1DiaSJKkocDbwPuyLMuVDPEEissds63+3Wt6vQZFiMElj8W/GcYe\nvTD26IUe8Azy4XxWvtXldPhQWYmOVNkqGFW6nMIrikS7dlgCAhvnjQkE9UyNxECSpJHA58B3siw/\nW8UwHRWdtLbfC6u7RlCQT01Mqlu8rVrl7+cBDnY0qk1V0BRtgqZpl9UmHwhtDlf3cj5pNsPx47B/\nPxw4APv3ozhwANWBA6j+WQ7/lFv8Nm8OHTpA+/Zl39u3h/BwUCpraFPTQtjkOk3VrovBZTGQJGkK\n8AYww2FzuDKOAa3KHWsNFMiynFvddbKy8l01qc7xLNTjBeTkFGEotSMoyKdRbaqMpmgTNE27XLLJ\noxkk9LF+OaA4f866erBnXpduYm/YgGLdOqexFg+Pql1O7s4L4kv2c2pgmqJN0DTtqgtxcjXP4Bng\ndeBFWZbfrmb4emBkuWOJwIYaW9fQCDeRwAFLswC7y8mJ4uLSKCdbiY7SOk4HU1EnJznPoVBgDgt3\nquNEjy4ogkKxNAtowHcjEFwYV/IMOgFvAXOx5gsEO5zOBwxAAHBOlmUD8CUwWZKkT4HpWBPV7gKu\nq2Pb6x4hBgJXcHfH1CEWU4fYilFOJ45X6FanTkvFbcUyWLHMPrQ5YG7e3Lp6sEU4lQqGObRNjVxO\nAkFd4MrKYDigxJpANqrcuZewPvGvBAYAa2VZPiNJ0vXADKxRRenAfbIsr6kzq+uJy6KEtaDxUCox\ntwnD3CYMQ+JAp1OK8+esq4c0GZ8TR9HvSUKdlopm62a0mzc6jbV4eGBsG4Op1OVkX1VEta3gchII\n6gpXQkunAFOqGeYUrC3L8lag90XY1TjYa+EIMRDULZZmARh79sLYsxc+QT7k2XzOxcWoDh9yKNFh\nzZtQH0pDs2+v8xxKpbPLyVbwLyZGuJwEF40oVCcQNCbu7phi4zDFxjkfN5tRHj9W1tLU5nI6mIrb\n8qWwfKnz8OZBFZLqTDHtMIeECpeTwCWEGDhR2s9AuIkEjU3pKsAcFo4h8VqnU4pz2XaXk201oU6V\n0WzeiHaTc5yGxdPT6nKKiXHOwo5qC26iTIegDCEGjog9A8ElgCUgEGOvQIy9ynlidTqry+lgallL\n01TrZrYmaY/zHEolpvCI0mJ/EiR0Qt0yDFNMOyz+lVWZEVzuCDFwRIiB4FLGwwNTXEdMcR2dj5vN\nKI9lONVxskU8qR1cTjYJMAe1qNzl1DpEuJwuY4QYCASXO0ol5vAISsIj4JpBTqcU2dmo0lJplplB\n0Y499tWEZtMGtBvXO421eHphjI7BFB1jTaqzCUZklHA5XQYIMXBErAwE/zEsgYEYA6+AoEEU3uyQ\nVavToTp00NnllJaGOvUAmr27nedQqZxcTo69Jix+/g38jgS1RYiBI3YxaFwzBIJGx8MDU8d4TB3j\nnY+bTCiPH6uQVKdKk1Ev+xuWORc4NrUItgpDaaE/Y2kWtrl1iGhr2sQQYuCIWBkIBBdGpSpzOQ10\nLiqgyM4ui3CyryZS0Wxcj3ZDuVpOnl6ljYhiHKrClrqctNoGfEMCG0IMHLHnnAkxEAhqiiUwEENg\nHwy9nQv+UVTknFiXZhUJ9YEUNHt2Oc+hUmGKiLR2qescj1tIhF0wLL5+Dfhu/nsIMRAIBPWLp2fV\nLqdjGZW7nJb+BUv/wtdxeItg6+ohOsap8J+5VWvhcqoDhBg4ItxEAkHDoVJhjoikJCISrr3e6ZTi\n7Fmanz1O/tZdZYX/DqahXb8W1jt31TV7eVuT6pwaEUmYIiKFy6kGCDFwRIiBQNAksDRvDh0iKW7f\nxflEUZE1yqm8y2l/CprdVbmcpNI6Tg5RTsLlVAEhBo4IMRAImjaenpjiO2GK7+R83GRCmZHulFhn\nrekkoz50EJb+6Tw8uGXlLqeWrf6zLichBg6IEtYCwSWKSoU5MoqSyChnl5PFguLs2XIlOkpdTuvW\nwDrnyvp2l5N9NWH9boqMAo2mgd9UwyLEQCAQXL4oFFiCgjAEBWG44krnc4WFqA8fdAiFtRb/U6ck\nV3Q5qdV2lxOdO1qjnGwuJx9fLgeEGDgi+hkIBP8dvLwwxnfGGN/Z+bjJhDL9aCUup1TUB9Pg7z+c\no5xatnKo49TOvqq41FxOQgycEG4igeA/j0qFOaotJVFtYdDgsuMWC4qsLGuU07bdDqGwqWjXrYZ1\nq52mMXv7VCwdHtPOGuXUBF1OQgwcUYh+BgKBoAoUCiwtWkBcW4pjE5zPFRaiPpRWGuEko7blTSTv\nQ7Nrp9NQi1qNKTKqYpRTdEyjupyEGAgEAsHF4uWFsVMXjJ3KhcIajagyjpYm1aU6rSbUaakVpjG1\nam2v4zTvCGQEhHIkKIJCD1++eC6xXt+CEANHRDSRQCCoS9RqTFHRmKKi4bqKLqeyfImKLqexpUML\ntZ48NPrz+je13q9wKSHEQCAQNASlLidDixYYrrzK+VxBAepDaSx78WPu2vIjeR4+FGvc690k0bbI\nESEGAoGgsfH2xtipC9Jpqxvpi/4PYVHV/3O7WBk4IsRAIBA0AbR//0nX9N2UJA7kkc+fbZAQVbEy\nEAgEgqaETof3y89j0WgoeHNqg+UqCDFwRKwMBAJBI+P5yQxUGenoHh6HKTqmwa4rxMARIQYCgaAR\nUR4/hueMDzC1CKboyckNem2xZ+CIEAOBQNCIeL36IgqdjsJpHzV4AppYGTgixEAgEDQSmnVrcF/y\nC4buPdHfcVeDX1+IgUAgEDQ2RiPeU57BolBQ8O57jVLgToiBA6KfgUAgaAw8vpqD+sB+ikeMrFjS\nooEQYuDIJVRuViAQXB4ozp7Fc+rbmP38KXz+pUazQ2wgV4ZYGQgEggbC6+3XUOblkv/ONGvv50ai\nxmIgSdJngFKW5YcvMOZH4HasXWJsj9v/yLI8qFZWCgQCwWWIevdO3BfMw9ghjuIHHmpUW2rkJpIk\n6XWgShFwoCPwDNAKaFn6dUeNrWtoRD8DgUDQUJjNeD8/GYXFQsE700DduI4al64uSVIk8CUQB6RX\nM1YLRAPbZFk+c9EWNiRiA1kgEDQQbj9+j2bHNopvGYahT9/GNsfllUEfIAOIB45WM7Y9oAL2196s\nRkKIgUAgaAAU+Xl4v/EKFk9PCl95s7HNAVxcGciyvABYACBJUnXDOwIG4HVJkgYDOmAR8KYsy/ra\nm9oACDEQCAQNgOd7U1FmnaHw+Zcwh4Q2tjlA/UQTxZV+TwFmYl1NfAiEAg/Ww/UEAoHgkkGVKuMx\n51NM4REUPTqhsc2xU+d5BrIsTwFayrI8Q5blZFmWfwCeAO6XJKlZXV+vThErA4FAUJ9YLHi/+CwK\no5GCN94F9/rvYOYq9bJ9LctyTrlDSaXf2wDnL/TaoCCf+jDJNfw8AfDxdsPHwY5GtakKmqJN0DTt\nEja5hrDJdWpt16+/wuqVcP31+I2486ISXe+d8gfumSc55x2AUaVJ//39m8NrPRn1IAaSJC0ENLIs\nD3M43APQAwere31WVn5dm+QybvnF+AL5eTqKS+0ICvJpVJsqoynaBE3TLmGTawibXKfWdul0BDwx\nEaVGw/mX38J0tqB2BhQU4P7zIl77YjrRZw7zdd/7WdxzWFjtJivjosVAkiQNEACck2XZAPwEfC9J\n0iTgNyABmAZMk2W56GKvV68IN5FAIKgnbE1risY/UaumNap9SXjMm4vbTz+iLMjHBzAplOxt07FO\n7KuNGJS/U/YBVgIDgLWyLC+SJMkNmAy8CZwBPpRl+d2LslQgEAguUWrdtEanw+23n/H4Zi6aHdsA\nMLUOwaxQoMzPY1HP20lr1Q6sof8XRY3FQJblxHK/r8GaV+B4bD4w/+JMawTEykAgENQDNW1ao0qV\ncZ83F/eF36PMzcGiUKAfOMhasqJEj99D92Po3JVrfv6UazQagoJ8Lmq/AEShOidECWuBQFDXODWt\nuX141QP1etz+XIL7vK/QblwPgDmoBYWTnqb43gcwh4WjOHOGgH69sLi7k//x56DR1JmdQgwcse/s\nCzEQCAR1QPmmNcqK0fzKw4fw+PZr3H+YjzI7G4CSqwege+BBSq6/seyGb7Hg89QElNnZFLz5LqZ2\n1SYA1wghBgKBQFBP2JrW6O570LlpjcGAdtnfeHzzJdo1qwAwBwRQNO5xiu8fiSkqusJc7t99i9uy\nvym5qh+60Y/Uua1CDJwQbiKBQFA3VNa0Rnn8GO7zv8Z9/jxUZzIBKOndh+IHRqG/cWiVSWjKo0fw\nevE5zD6+5E//pNIVxsUixMARUcJaILjkycjMZ+Zia57rhNviCQtunMQ1W9OagjffRbNzG+7fzEX7\n7woUZjNmXz+KxjxC8f2jMEntLzyRyYTvhEdQFhaQN2s25tA29WKvEANHxAayQHDJM3NxEtl5xfaf\np43r0+A2qHfvxGP+NwB4zJqO6vQpAAzduqN74CH0Q28FT0+X5vL4ZCaaLZvQ33Qz+jvuqj+b623m\nSxEhBgKB4GIwm9GsXon/XWUFGBT5+ejuH0XxAw9ijO9co+lUyfvwmvom5qAW5E/7qF77tNe940kg\nEAgakQm3xRPo604zbzdaNfdk9pJkcgrqt3q+IisLjxkfEtCri5MQ5E/7iHNJMgXvfVRjIUCvx3f8\nwyhKSsj/aBaWwMA6ttoZsTJwRKwMBIJLnrBgH6aN68PsJclsScm0Hx87NO4CrypjS8pp5vyeAsCY\nIbH0im1Z+UCLBc2mDbh/8yVufyxBYTA4nc7elXJRvQq8/vc26pR96O57kJJrr6/1PK4ixMARuxg0\nrhkCweXMxW7wfvFnMmkb9pHtHUDvzq0ZfaNrN3lXmfN7CmZL2c8VxODcOTxmf477vK9Qp6UCYJTa\no3tgFOrkfXgsmHfRTWvUmzfhMesjTOERFLz2Vq3nqdE1G+QqlwpiZSAQ2Fm58xjzl6cBMGJQDIkJ\ndRPFcrEbvIfX7eaLueNICYnlneJnoAoxGJ4YXenPlZFToGfhSmtR5Ur/+i0W1Nu34vHNXFjyC97F\nxVi0WopvuxPdAw9h7NUb1cE0vF9+4aKb1igK8vF9bCwoFOTN+hy8vWs9V00QYlAZQgwEArsQ2H6u\nKzG4WM55BZDr6UfHE8l88N1kVMNjMMVVrNzp7+3msmto4cqDdpdS29a+HDmVB8CjiWG4z52Dx7yv\nUKfssw6Ojqbg3pEU33Wv3Y+fk1+M8uHxBBiNZE55A81FNK3xevkFVBlHKXr8SYy9etd6npoiNpAd\nqb+NeoHgP0FOgZ7ZS5IvuGlr2+AN9HVnwm3xNb5G94QwliQMASA47wzNbrwW7d9/Oo35Y9MRRr27\nklHvruSPTUdqNH9zfw++vr4ZP2T+wnV3XIXPc0+hSj2Afsgt5Py0BGQZ3fjHnTZ0d30wl7bJW9gR\n3pW5irY1fk82tMv+xmP+Nxjj4il85oVaz1MbxMrAEeEmEgjsjBgU4+QmcgXHJ2ywumds7pfhidH4\ne7vZN3hry+gb41Bc9Tbmrr+hzM+DYh1+D9xN4fMvUTTxaVAo+HlNmQD8vOYIN10RecE5hydGo9Hr\niNu6goHT/8Y9aTcApjZhFE98Ct3d92EJDrYOLp/9q9MxaOF0DEo1cwaMJrSW4Z+Ks2fxmfQYFq2W\nvE/mgFZbq3lqixADR4QYCAR2EhPaXLRrqLw4uOK2WbQ6jb83HwNgcO823NG/ohBZfP3IGfEgAZ9O\nZ0XHgXQ9spOgd95Av3svpk9nVxj/yHurmXBbPHGRFcMzVSnJhM6by+RFC1Hm52FRKtFffwPFD4yi\npP81oFJVeI0jnp/MwOvsKTZedy+hfROq3Z+oFIsFn8kTUZ7NouCVNzF1iK35HBeJEANHhBgIBDXm\n2+UHOLJsA6f8W9MtIYxesdYn6OGJ0Xy7TLaPKzGYnF5niyoqLjFSWGwErCsQmxAA/L35GEs3H0Oh\nqBjm+U30IMapPqXz0d08fddUnv1zGrF//4bh5gxGPDGV+fvKGiuWGM3MXJzEZ0/3tx7Q6XD7/Vdr\n05htWwAwtWxF4dhxFN97v8uRQI5Na9p9MpUYF3oVVIbbj9/j9ucSSq64Et0j42s1x8UixEAg+A/w\nxZ/JbEyyPqH3iQ+2h2PabshKlYLxt3S0h3nWJPzz5B8rmfnDcxRpPdgQ0wfv0Q8QfffQaoupOUYV\n2XDctLZhwfp8Vj7Ms9AvgH86XsONe/4m/vg+ptz+BuNXzWbgnn+4ZfK9yLdOYZtvhNNcqoNpuH8z\nF/cfv0N5/jwWhYKSAdege+AhSgZdD+rqb4mHT+Ty+hebAfhg46wKTWscI5NsrrELoTx+DO8XnsHs\n5U3+jE+rXYnUF2ID2QGL6GcguEyxCUH5n2035KzzOvvN3/F4dl6x0/HKyAgMI9/NG88SHdcm/8sV\nk+7Hp3MsXm+9Rqussqd8rca1m1yQv2uROMMTo9k3bCQmhZLbty7GpFIxfeB45vQfhfpcNs99NZn+\nKatRmwz0T9vAF/++Q0CfbnjO/hhUaooef5Iti1dzd48nGJXiTfKxXOf3lZnP5E82MvmTjWRk5tuP\nv/nVFrLzignZt5WA5X9UaFpjc41tScm0i0KVmM34PP4oyvw8Ct+aijk8wqX3Xh+IlYETwk0k+G8w\n+ZONFSJ5zGYLs5ck23+2kZ1XzHOzNxLZyg+w3oRPZBXw0aK9mMwWcPPk5x638MD6+WyL7MZ5T3/6\npm3Ec/r7PArcHN6epCsGE3PXGKfraTWVb7Rm5RTTNaY5ew5mY3b4WzRb4OW5mzl+psz9E+jnS0qv\nQcRvXkr3w9vZ1rYnSxKGcjoojJcWvcpTSz/iqaUf2ceX9L3aWi568E2g1TLtvdWUGM0Azm4kLpwP\noTSbeHjVF5hRUPDOtFqXlPb4/BO069da9yjuHlGrOeoKIQaOiBLWgsuUPvHBTisC2xP/hNvi7W6i\n4GYe9s1ePy/nSJYz54s5c956Y3TcELbxZ5cbuWXHEjqc3M/ohz7n88SHuebkLu45s5WWG9bQOv0A\npoXTOda1D4tDrmBrVHeKVFW7T7QaFRq1En25fQZHIQDIztWTPXYCbF7KHVsXsyOyGz0Pb+PRMxsr\nzDnpnvcITLwSgOElFvxrGazz4oO92Dj6OcKzM8i67V7o3NXpvKvJbqoD+/F66zXMzZuT//7Mei1C\n5woKS9O68VmysvKrH1VPaFavxP/OW6whapMmAxAU5ENj2lQZTdEmaJp2CZucmfzJRvvTbqCvu/1p\nNyjIhze/3Fzpjd5Vbt32C6PWfcMPve5kwZX32I/7F+ZwlbyOxJRVRJ85DECh1pMN7fqwMnYAKSEd\nsCjKnqwVCugS3ZyBvcKZ/sMu+5P7hZjy29v0PrSVHA8//HVWd09K6w6sk67kytRNdDyRTEZgG964\n+QVO+7eyv06lAIVSARYLnu4aVEolIS288HLX0De+JV//bd0Ad9w3CaIYc3QMKBSc27QTS/PmNf+w\nSkrwH3wNmqQ95H7zPSWDb6z5HA4EBflctJKIlYFA8B/CthIwmy20CrRW9ByeGE1QkA/DE6M5eDy3\nwqaujfBgb9IzC6qc+68ug7kn6Q9u3vUHSxJuIt+jdEPVy5/fE4bwe8IQ2mQfY0DKavrvX82gff8w\naN8/ZPoGsbpDf1Z16M+JgBAsFtiVdpazucUuCQHAjz1vp8fh7WhMBv7ofAPLOw/iSPMIq12dBzNq\n7dfcvPN33v/uGabeNJm9YZ0AMFlAq1Dg4+Vmf9/nHZLlKs2HeOEFlHm55L8zrXZCAHh+MBVN0h50\nd4+4aCGoK1SvvvpqY9vgyKtFRSWNdnFV+lHcF/2Aoe/VGK6wLie9vNxoTJsqoynaBE3Trv+qTRmZ\n+bw5bwcrth2jfZg/FuCbpTIHT+Ry77Xt2HbgDOmZ+ZzIKiSnoISrE0I5eiKHzcmZFJeYKp0zt/DC\nNptUakrMCnod3orWXcuOkIrZxXmefuwJ78yShCHsC43DolAQnXmIrhl7uWn3X3Q7soMirScZzcPI\nq+J6lT0Cn/MJZFWHfvxwxXC2RPci19Pffs6iULIzIoGC5sH0OLCRxJRV5Hn4kNbSmr9gMlvQ6Y0V\n5gwN8qbKSX4OAAAgAElEQVS71ML+e/KRbL5+cwEDv57KqVZRpL8yjR9WH2aHnEV0iB/uWteerdXb\nt+Lz+DjMoW3Im/c9uF042sgVvLzcXrvYOcTKwBGRZyC4TCi/+Rkd6md3AVX19D9zcZLTU3Ft+LvT\ndQzb/gsDt/zGDx1vIM/Tr9JxFoWSvWGd2BvWic8Sx9L74GYG7F9Dt6M7GffPZ6yX+lZ5jcr+OpUK\nyPRvecExv8cM4NAdwTy/ZCqPrvyc8LMZfD5gNCaV822wY1SA3U00aeZ68otK8PHQoisu4a1/PkeJ\nhZl9HyRl/i4c9tldq4NUWIjP+IfBYiF/5mf2cNSmgAgtdUSIgaCJUlWYY1UU6Q2V/gw4CUGgr3u1\nGbNjh8a6vLdZonFjUc/b8DAUM2rt1wTmn632NXqNG2s69GNhrzsA2BGZ4NrFHDC7+CebEhLLk/dO\n43BQBDfsXcrri1/FR5fnNObJO7swdmgcX/8tk1tYgtkCuUUl9E1aSftTMuvaXcm+NvEuX9MR79de\nRH3kMLpHHsPQp2rBawyEGFSGEANBE8Mx7v/9H3bzyHureeS91SQfya60OJxOX+bq0elN9I2v2KAl\npLknU+7vZk+KcuwQ5sjsJSlOfxLeunyu27sM9xIdAL66PPrK63n2j//x2we3MnbVFwBck7KKOV8+\nQtvMQy69x4H7/gFgRceBLn4qzijNlbu3ypPl24Jnh7/DxujedDq+jw++m0zY2YwK4wqLy0TUQ1/E\ng+vnoVdrmdtvpNM4VwQVQLNyBR5ff4mxQyyFz7/kkq0NiXATOdLIoV0CgSvk68puUjMXJ9G1XVC1\n9X9sUTGOnDhbhL+3W6W1gEZPXVnhyTcgP5uJy2fQNX0PAI/982m1tmrMRl769S2evGca53yqbtvo\nXqLjqtQNZPq2IL15GKHZxzge6HpdpEf/+Ywb9i7FoFSj12jRq92sX5qy7yVqN/u5Pr3assyvBec9\n/WmZm8m075/hfzdNZndUN/ucjnsnd29eiH9hDoXPvchLj93B0zPWYijd3A4L9q42y1hx/hw+T4zH\notFYexRcRInr+kKIgSPCTSRootiigADyCkswmCqPssk8V8iod1c6HStf76c85WsB7ftzI7enbaJP\n2ibaZtWs/DPASf9WrIztz4iN3wMQWHiOl357i+eGv41eU/lNsG/qBjwMxfzc/VZGrpvHwOSVLOpx\nG/OvvAezsvrM5UI3LwD0Gi1nfFvgZijBzajHp7gAN6MercnZVUbKKm5x+NXTUMwD675lR2Q3Hp62\nCouDEoaeO86QXX9Ym9aMe5yoED86RgawK83qAtuVdpZR766ssqgegPezT6LKPE3BlFcwxXeq9v00\nBkIMHBFiIGiiOJZ9Tj6S7VQ3KCSorBNW+TyBuc8lkpGZj9GYjlLh7Ft30yoq7D/csWUR929YUGP7\nDrWIYmmn61gjXYXOzROAhKO7iT25nwOtJNqfknnqrw95Z+izTjkFapOBtpmHeGL5LADWS1fSrOA8\nA1JWc8e2xcSeTOF/Nzx9wVUFwE89hzFo3wo0xhJeGfYKOV7NnM4rzSa0RqtAuBn0uBlLcDPocTfq\n0Rr0uBn1pDcPA8BocviQLBbGrPoCtdnEwadewWJU8PX87WRUEmL79+ZjlYqB8bvvcf/1Z461jad4\n5KP4VxjRNBBi4IgQA8ElQOb5Inv8feZ5a0burtSsKsfPXJxEbmkoq6Mg6EssvPrVNqexR4LK6v4f\nDQwjItvqSzeo1GhMzuGXRRp31ra/mqWdBnEouKLP/Kcet/Lyb/vJ9G2BXq3likNbGLl2HtuiuhN3\nIpmOx5Npf1LG3VgWwfTWopf4qccwfu5+K3dsW0zcif3MmD+JD66fyM4LbCwXuXkxv889jP/3M+7b\nsICZgx5zOm9WqijWelCs9ahyjsrofWgLCem72RHelZnZLYl2KMkd6OteZU6GDeXJE/i/8DQ6jTuv\n9xtH+JojLndfa2iEGAgElxjl21Fq1cpKk7PCg70rdBtTKRWYTVU/7GyP6s6QSb9wtbyOsSvn2I87\nCsGxgFB+TRjKuvZXobvAzXV7VHfSA8O4Mm0jT94zjWf+fI9hO35l2I5fKx1/NDCM4LwzPLz6S847\n5An46fJ47ZfX+bHn7Szoc3eVbqPl8ddy0+6/GLjvX/7ocgNHWkRVaZsraA16Rq+ei1GpYs6A0c4r\nBiA61I+evkFO+y1OmM34PDEObVE+Hw98lNP+rQi/KIvqFyEGjoiqpYJLEGO5/QPbE2t6ZgELVx5k\nwm3xTF+0l8JiA1GtfTl4PIeqEnv9C88zYcXH9Dy8vdLzc696gF963OqSXRaFksU9buXJpdMZmLyS\n1295iYfWzOW0f0vcDMX0O7AOD0MxJ/1b8fHAR9gb1hnfolxu3f4rN+3+q8J8d279idgT+5l241Oc\n8w6ocN6sVPFF/1G8sfhVRq+ey5Q73nAKCukVG8zhk7lk5VT9NP/U8M6sTzpNoc5A3PefEpx3hp+7\n3cKJgBAoKsHTXcXVXUPQ64328tRV7RO4fzUH7ZpVFPYfyL5rboWsIrakZBIS5Flt57XGQISWOmAR\nVUsFlwCOLSjDg70rRP2cz3e+2YUF+9AuzJ8So5kDGVUIgcXCnVsW8e3sB52EQG4ZQ6ZvEADLOl7r\nshDYWCtdxRmfIK5LWk6huxezrh2Pf2EO1yetQGMy8GPP25lw30fsDesMWDOUv7n6AR4a/TmLu99S\nYb6OJ5KZ/u0kuh7dVen1dod3YVtkdzod30fvQ1uczm1JybygECiAuMhAxg6N4+krA7lrxy+c82rG\nD73vtI9ZtfMkk0d0Z+zQuAtGEKkOpuH9+suYAwIonvkJx7PKCuw5tuRsStRYDCRJ+kySpM+rGdNd\nkqT1kiQVSpIkS5J0X+1NbEDsewaNa4ZAcCESE9rQJ97aTcyxVlCgrztatdJJHGJCL5zhqjSbuGH3\nX/z+4a3c57Bx/Fen63l8xAfMGDSBgILzZHsFVIivdwWTSs2v3W7GzVjC9G8n8enX4+knr+NAy3ZM\nvPd9vu07ghJNxZtqnqcfX189khGPfF3hnL8ul1d/fp17NyyoNLdgbr+RGJUqHlzzNWqjocL5qrAA\nn/y6l1HvrmT/XWNQ64vJnfKqfUPcxjd/JF9wnpycQnT3jECh05H5+jTMwRVzPJoiNRIDSZJeBx6u\nZkxzYCmwHegKzAS+lCSpdpkkDYkoYS24RHAsRw1WF8iU+7vh4+lcl3n+8jRGvbuSLSmZFWv6WCx8\n/M3jPLqy7Nnu13738tCzP/PpwEdIbx7OhOWz0JiNfDzwEYpKwzdrSnKotZ9v84JsPPVFfJr4MM/e\n9Q7pQRHVvjbX05+bJy62/16i0gCgxMJdWxbx5k8vE1Bwzuk1xwNC+avzYFrnnuam3X/WyNbtB87S\nKWMvfVM3sr+VRPqAmyrsBfy06sINa05OfomQo/tZ1aEfU4sjGT11pdNnP6xf03MRgYtiIElSpCRJ\nK4GxQHo1w8cAObIsT5RlOVWW5VnAfODpizNVIBBUxdihcXz3j3zB6BbHR5zm+Vm8+ssbhJ4/YT9m\nVKrQ5Jwjv8A6x5Bdf9D+dCprpKvY1rZnjW1SmwwM37yQ6fOftB/7rdtQ/upyg0u5AzbMShUrO/QH\n4M2bX+C3hCHoVVbRiz9udRt1Sd/t9Jrvew8n382b4Vt+xLcot/yUVeLYtGZ24hhm/pJMoG/1CWLJ\nR7J55L3VTHtqDn1/n0uWdyCzB4zh0Mk8zBbrZ69UWEN9m+J+Abi+MugDZADxwNFqxvYF1pY7thq4\nsiaGNQr2/WOxMhA0bWxuIseftx+ovg4QgF9RDjPnTaTb0Z3sDO/CQw/N5r3Bkzjj24Ib9/zN7dt+\nptX5U9y3fgE5Hn58njim+knLEZ+RxC/T77Anntm4dcdvLtUrKs/OCGsDmbZnDvNF/4cY89BnLOl6\nEyUqDf66XF5b/Br3bvjO7jYq8PDh+yuG460v4p5NP7h8nRt3/014dgbL4wfaw2XL92W+fUBZGK2t\nFMiHP+4BnY4Jf36Iymxm0b3P0zGh7SVV1MAlMZBleYEsyyNlWT7jwvBQ4ES5YycBT0mSKoYANCVE\nnoHgEmH0jXHMfS6Ruc8l2pvbV4a3R8WAQb3ajXx3a6KaymyiRK1lTYd+PHPXOwCEZWcwYcUs3Ewl\nfHX1A5jVaoJNBQTkZxOcm0noueNEZB0lIuuIs8/eYiHhyE5+/+AW3v6prPZOlk9z1jpUIf16zmie\n/f1/dEnfjcLiWr+CXeFdMKMg4ehOAM57BzBnwGjGjPqM37vcgEml4q4tP9LVYYXwV+fBHG8WwvV7\nl1WoPeThVvHW51uUyz2bvqPAzYtvr7S2oHzwhvYVxp3J0dnrQNn6HZst8MD6b2lz7jhLut7EMu9o\nxg6N4+EhsSgV1lXBmCGxLr3XxqI+Qks9gfJrVVuwc9MryOGIEANBEyb5SDbTf9prj3cfOzSWXrEt\n7RnJCkXF/7olhoo322KtB5PufZ+Jy2bQ+9BWZnw7ifdueJK9beLRq7T0PrTVPnbSshmwbEaVNp33\n9OdAKwmTUkXfNOdWk/tbSXx6zVhrIptCwY89b2fWtxMB6Ju2kb5pGzntF8zyjgP5J+4azlcSLmoj\nz9OPg8HRdDh5AA99kX1T95xPIJ8nPsziHsPodmQn+0LLhNGkUjO330he/vUtHlozl1eGvWL/G1cq\nlIDzZ3P/+vl464uYPWC0vfT2hqTTjBgUY18dRLT0Zu2u8s+60Dl9D0N3/cGxgFC+6XsfZou1qxzA\nyyN72LukNWXqI7RUB5QPD7D9XlgP16t7hBgImiAzFyc5JT7N+T3FfrzEaHb6b9vz0FYeXPM1HQ7t\nQl2uLo9GraRZmxa8NfR5vug3Cp/ifN746RXu2vwj2d4B5Lt5Ux3GUp9/s6Icrji0xUkI9oXEctuE\nhTxz91Rr4lfpDTg9KIJ/YwcA8GPP21gedw3+hTncv2EBX80ZzfNL3iXhyM4qq4/ujOiK2myi87G9\nFc5l+zRneadBFWofbYvszq7wziSk76bbkR3244XFztnU0acPcu2+fzjSPJy/Og+2H085eo6Edi3s\nq7DgAOdN9OGJ0Vwd4ckzKz/GrFLz0eCJGLVueHuo7VVmbaVDmjr1sTI4BrQqd6w1UCDLcrU7OUFB\njaigAdY/Ak9PLZ4OdjSqTVXQFG2CpmnX5WKTooIDWkFQkE8lx+HqA+voJ69j2I5fKdJ6sDO8K9ui\nurMzqjs5Hr7WuHeFgt+6DeVAK4ln/nyPezeV+fcfeHguLyx5B+l0Gmf8W7I8fiBdDu+k4wmrAKkr\nuWGnBbflvRue5GSzkCrfw+Iew7gmZRWdMpKYfPdUvuj/EP0OrOX6vcvoc3AzfQ5uJtM3iBUdB7Ii\nbqBTTaIdkQncteVHEo7uYnN0b1c/NL7oN4oZ305i9Nqv2RvRFUO5zWuFxczYVXNQYuHzAaOdNrfN\nFvh1w1Emj+gOwLg7uuC2ZB8ADw3tSICvO5O3zoPzZ+C113j/5acAGPXmcgp0VsFRqhRN8v9geepD\nDNYDI8sdSwQ2uPLixmxers4pohlQVKinsNQO0VDddZqiXZeTTTf2CXNKWBozpANZWfk8NqwjMxcn\nYQGiWvkiH8th5qDxZHsHcvPOJXiW6OxuGTMK5Fbt2BbVna1RPUhvHo7cWuKJER8wadl0e8LZoKTl\nJIfEIZ1Oo0XOaUasm1+tfTGZh3ju9/+R1CaepNCO7AuNo8DD+SZ4LLANm9v2pPehrXQ8vo99beJZ\n2vl6lna6jujMQ1yXtJyrD6xlxMbvuXvTQrZFdWdZ/CB2RnQltWUMBW5e1n0Di8XlkvMZzcNZHn8t\ng/cuY9CepfzZ1bnn8ICU1U5Na8qzbtcJu2vIFhaq1xvJzi5A9dMi/ObPx5DQjZwxEyArn5wCPcH+\nHuTm6/Fy1zD+lo71/n+wLsRGYamhS0SSpFVAmizLD5f+rgECgHOyLBskSWoBHAAWAtOBa4FpwHWy\nLK+pZnpLo4rB7p00G9SfokcnUPjaW8DldTOpb5qiXZeTTY+8t9peg0irVvLZ0/0rHff49LX2p9K2\nmQeZsOIT2p45bD9vUihRlW7cnvEJIrVlDM3zzxJy/iQ++qob3gOkBscw69pH7fsAapOBmNMHiT+W\nRKdjSbQ/KeNmshbFM6PgSFAESW3i2dumI8mhcRS5eSGdPMB7PzyHSaHkcIsoUkI6kNK6AykhHcjx\naoZHiY6rDqzj+qTlxGRaY/qzvANZ0XEgnY/tJe7Efh4dOYvjAaEuf3Z+RTnMnjsOk1LFw6M+pbB0\nA91DX8Tsr8bhWVLEIw9+zFmfIKfXla/06kizgnN8/v0k3IwlnF+5HlNba2b47CXJ9mJ2vWKDG6Qw\nXVCQz0XHLdVmZVD+o+kDrAQGAGtlWT4jSdL1wAxgJ9a8hPtcEILGR2wgC+qRPzYdsT/ZD+sXWW/x\n5jYhADgUHM2T90zj5h1LuHfT97gZS0gLjmZLdE8istLpdnQnfdM2YlCqyfQL5kDr9pz0b8nJZq05\n5d+KU/4tiTuRwsRlMwEwqqzjbH8rRpWG/SEd2B/SgR9734naaEA6nVomDqdk2mYd4ZadS+w3/6TQ\njmyL7E7z/Cwiso4Sk3mQm3f+Dlh7Iexv3Z6UkFg+GDwRN4OeQftW0H//Gu7ZvND+vhKO7qyRGOR6\n+rOw1x2MWvcNd29eyBf9HwKsTWuaFeXwbZ97OOsThEqpwFR69x8xKIYfVx6qtAggFguPL5+Fe34u\n+e+8ZxeCS5karwzqmcZdGezZRbNr+1H0yGMUvv42cHk9WdY3TdGupmRT+aYzc59LrNHrbVFDZosF\ni9mCSqXkzsS29qqZE26LJyzYp8J1bLTMOcW4fz6ja8YedBp35l95D391Goy/LodzXgEXTARTmYxM\nWjqdfvI60oLb8sqwV8j3qL6Zu9agp/0pmfhj+4g/nkS7U2lozFaxMimUpDcPo1jjgVGpQoGFyKyj\neOvL4kzy3H1ICenA4aBIAgvOEX3mEG3PHGZx91v4+uqRNfj0wN1sYMZXjxGUf5bH7p+BAgsz5z1B\nlk9zJoycicXNnchSNxtYn+r7xrfko0V77QIxrF8kP685wvV7lzH+n0/ZFd6Z0C1rQFkWi2MLOQXs\nxezqm7pYGQgxcEC9dzfNBl5N0djxFL5hjbluSjcTG03RJmiadjUlmy5WDGw4uoscCfR1Z9q4PhWu\n4xRyarEwYP9qRq+ei29xPqnBMayIH8jZFqEc82tNprt/lb54pdnE+BWfMCj5X44GhvHS7a9VaCJT\nHW6GYjqcPFC6cthHzOk0u8vKqFRxMDiaHE8/jEo1SouZtmcOEZxX1qtBr9KS4+XHxpgr2NOmEwda\nt7e7fCqja0xzjp7K53xpKe+rDm/mmV/fZUtUDzQmAwnpu3nj5he468NJ+Hu7MWXOZk5lW4vKtQr0\n5K0xFTeqt/21mcTRQzCq1Pz56S8k3uTiZnY90lhuossXUcJaUI/YniptP9c1JQYTj7y3GgXO/4MV\nWMNJS4xmUChYFTuAHREJjF49lwEH1tAusyzDVqdx56R/K04EhHDSvzUnmlm/TjZrTaG7N7MGjUev\ncWfI7j95d+EUXrzjtQp+9guh17izO7wLu8O7AOBRoqPDif2VioNBpUZuKbEnrDNGpRqV2UR05kEi\ns45y644l3LpjCWYUpDcPY3/pnkNKSAeyfILsf8u70s46lZNYF9mLwaFx9DpsbeqzI7wruf2utT+9\n24Sg/M92jEZumP0KGPXkzfqkSQhBXSHEwAFRwlpQn9x0hXWf4GJXK479kB3dRFX1RrYAUa19OZCR\nYz+W5+nHBzdM4ucetxB15ggh50/S0ZRNO10WoamplfY+zvHw42Sz1pxs1opsr2aE5Jzk3YUv8OLt\nr3Pav3w0uWvotB7sjEywdzHz1BcSdyLFLg6xJ1LoeMJaJVSv0nKgtcSv3W6mRK1FbTYSczqN9qdk\nIs+mc8PepYB1szklJJY/u9zA/pAOhAV7k19UYhfDeX3vY9oPzwHwy9BHGXtLR5ft9fh4OmzaRPEt\nw9APu6NW77mpIsTAEbGBLKgFdeEjLt/XOC6y6p6/cZGBTpFEiQnWqppjpjq7h5Sl/50jWzkLgY2n\nhnfm/YVwtLTVZatAD05l61BcZyawIJvW508SkX+almdP0PLcCULOn0Q6JRN7cr99juC8LKb89g4T\nHnDOUvbz0qI3GCkuca3chI0iNy+2RfVgW1QPALyKC+h4PNkuDp2PJdH5mPVzKla7sT+kPT/1GIbB\nywff81kM3fU7QQXZ9JPXoTKb2B/SgV1pZ1Epy7wonTP22H9uuX8X/t532X9vFehpXxEolQomf7LR\nvhejStqL1//ehtatKZj6QY3e16WAEAOB4CJZ6NAXF3AplHDR6jSndon/bj9h3weYuTiJB29ob88w\nHjPEWnaiMhyFyNE/pFYp+HyyNdt39pJkDp3Mq/Da8oJzKlsHWDuUnfUJQhMRzt5yrhK1yUBw7hla\nn7eKQ8j5ExwLLNfuEcgtLMHLXU2v2CCnz6amFLp7syW6F1uie6E0m7ji4Ga6pO+h3elUorKO0jV9\nD13T91R43bbIbnx6zVj77yazBa1aSb9AE3dsLSuJPWLDdxjzXsDiay0/MfnurixceZBdqVmUGM32\nDOJpoxLwHT8GhcEAc+diada0y6zVBiEGjoiVgaCBsAmB7Wet2rkyzJzfU+zx7bOXpFBYbLCvABxx\nFCKHh1+UDpvAwxOtVTZzC/QcOmEtAqBUKqqMOrJRmc/cqNJwIiCEEwEhbLvgq60lH8YOjWPf4ewK\n5R9sKCxmYk4fxN2gQ2M0oDUZ0Ni/jGiMJfbfex/cQrRDvsSF6HFkBws+ewCAL68eSVKbeHLbtufh\nDV/hZizhg+ufICj/LPdtWEDRR+9T+PLrAPh7uzF2aByTP9noVA7c6503UB/Yj+7B0Xhcdx00kaCE\nukSIgSOiuY2gFthutuV/rgmO+wATbou3lkR2YP7ytErFwJF2bfw5XLoCmHBbWSatnHGebfutgjFm\nSCy7D2ZX+rQe2sKT42cq2TR1gUBfd8KCvdmVVrE8dceowCpXBzfvWMJDa7+u1TVdpbL5Rz5zJ6ao\ntpj6rsXj80/Q3f8g5oiyTX3Hf4/nwwrweG0Wxqi2FLz8Bh71am3jIcTAEbEyENQC29NkTRjcu42T\nm6j8PsCYIbHMXpJS7Tzlhaiy/QrHVcac31Po0SG4whilQsHro6yRMZWtGAb3bkNuYUmFDmtgLZM9\n4bZ4/t7iXCZaobDuhQxPjCYjM7/SlcZa6Spa5J0hMusoEWeP4q2vXIxO+rUkOO+MPdIIILd1OHM6\nD2N7ZDfiY1ux7XAupnK1N5vnZzHW7yxd0vfgvqisr0FA4pWYAwJQnrN2SfMbcSfn12213wPCgn2Y\nNq4Pivw8mvXvAwoF+bNmg1ftur1dCggxqAwhBoJ65o7+MdzRv+qs1V6xLSksNthLJ48YVPnY8kKU\nkZnvtMKorHTy8MToC7puyjP3uUSnectToDPy+tfbuOdaq43b9lvr+1ss1v2Pz57uX3mYJmUlqAGw\nWAjKP0vE2aPWfgln04nIOkrYueO0zj1d4bVp2kBanz9JJ7WWo7p8Xpp4E+//uJd8nQGNSsnjt5dt\nxBvnfAqLfqCk3wCKb7sT7fq1aDasA6xioE6VCQr2o/jmYRiuvApD36sxtY3Ge8qzqI5lUPjkZIzd\na97t7VJCiIEjl1JbIsFlT2JCm2pdQ+Vv/jMXJ9l93TMXJzFtXB/GDIl12oz293Zj5sSr2ZJy2n78\nqXsT7HN5uKnQ6a1VSW0i5DhvZZgtVldWr9hgVCol5spKODjw6oM97HZL4X7WFYdCQZZvEOE9O7Co\ntGube4mORbPurnSO7kd30r202Q0AXz3Kd8AJ/9a88egsuxAozp7Fc+rbmP38WT3pHT5edxraDGfM\n/Fe5wlOH58wP8Zj/jfV6v/2M+28/V7iWftidNSqOdykixMAR4SYSXGKUv/lXRq/YlhWikXIK9Ow+\nmE2PDsEMT4wmJrI5I19bZp/Lls1cU7akZCK18efIKee9C8eQzVaBnnY3jI3y3dpW7jzG/OVpFGs9\nmHDfh0RkpeNTXIBPcT4+unyiPYx4H00j9HTFfIiQnJMEZ5c1oPF6+zWUebnkv/0/Pl53usxl9sd+\nej2bSMEHM1HodLgv/pGiJ57CFNoGt99+Rru+rHtvQN8emFqHYLjyKrj+WpSdemAOj6jx59OUEWLg\niF0MGtcMgaC2OG58jhwsMXuJNWHLcT8hp0DPW/N2OD3pv/hQc5fmvdDqwIa/jxuf3dvf6VhYsI9d\nDFzp+pWY0Ibh18Xy5peb2ZJSlgsB1oqtL9zXjZ+2ZDhVB70hoSXzv12P0mzi7tHXAtZKxO4L5mHs\nEEvxyNHwfvn27FYKX3wVt79+x+2HBZzbtBPtv8sB0N03EmNsR7Qb1qHZuM6677DoBwIBU5swDFde\nRUmpW8kc4nrhvKaIEANHxMpAcIlRPgrJ8YnbsZQylOU/LFx5sNKbevm5HLHNW1U46rB+kZzIst7s\nK4uoqm3ElW1sicHEviPnMBjNlBjNvPrVNqd9FJvYPf/8zWUvNptxf/ZpFBYL3984nt7FpgouM/vQ\nkFCKxj2O1/tTCbiyO6pTJynpezUF0z4CpZLihx4GsxnVgf0E7NmKfukKNJvW4/7DAtx/WACAKSLS\nKgw2cWhZu6zsxkKIgUBwCVPe3VJTAn3d7TdcV+ZyLPFsw5WCe7WJuCr/uvKx/9+tSOOLZ6u+ttuP\n3+Oxazvr2l3Jj+YQ0lceZOzQuCoT+Ioem4jX+1NRnToJQP6MT52qkaJUYoqNg369ybv7QTCbUScn\nodmwzvq1cQMeC+bhsWAeAMa20Zzp1JPvjSHsjejCgyP6XDCzvLERYuCIWBkILiOqehp3JRy1MjIy\n8/pX4TMAABWpSURBVPF0U5OvM1Q/uB6YcFs8r35VXaqbFUV+Ht5vvIJB68bcfiNdu4C7c/9kc+iF\nN+9RKjHGd8YY3xndI4+ByYQ6aQ+a9evQbFiLZvMmWh/6jqeA037BjPeYU2VDoqaAEANH7IECQgwE\nlz5VPY3X5Cn9iz+T7bkFWjWUlItGbcjgmrBgH8YOrdzNUx7P96aizDrDuSefo21sR9pSvXvK49NZ\nTr+rt2+tWTipSoWxSwLGLgnoHn0Mzw/+h9e0dzCjYEXcNa7P00gIMXBErAwEAicck8zKC4FWrayw\nt1DfVBYZVR5VWioecz7FFBaBZeKTjC33xF/pa1KS8Xr3DcxBLcif9hF+I+/B+6XnyPnr3xornvJY\nBj7jH0a7eSP6oJa8nfgY+8I7NfhnVVOEGDggSlgLBBfG1hugqoQ2V3ElOa5WWCx4T3kGhdFIwRvv\nVHD9VIpej+/4h1GUlJD30SxKrr0e/ZBbcPv9V9x+XoT+tjtdvrzbLz/hPXkSyrxc9DfdTP7703ns\nEilqp6x+iEAg+K/SJz7Y6edp4/owbVyfi75528JUbVVB6wrt33+iXb2SkgHXUHL9DS69xmvaO6iT\nk9DdN5KSa68HoODl17FotXi9+SoUVV+vSZGfh89jY/EdOwqF0Uj+Rx+T9+W8S6q6qVgZOCLcRIL/\nKGt3Hef9BTsA55LZo2+Mq5AQ1mTR6fB++XksajUFb/3PJfeOestmPGZ9hCksgsLX3rIfN4dHoBs7\nHs+ZH+L52SyKnnym6jm2b8X30dGo0o9i6NKV/M++xBRVu4KFjYlYGTgixEDwH+X9BTsxW6xlJWwb\ntPXJhNviCfR1J9DXvc586Z6fzECVkY7u4XGYoquu+2RDUZCP72MPg8VC3qzZWLydVztFE5/C3DwI\nzxkfoDx9quIERiOe772L/5DrUGakU/TEU+T8+c8lKQQgVgbOiBLWAkGDcLH5EeVRHj+G54wPMLUI\npuipqp/iHfF6ZQqq9KMUTZiEsfcVFc5bfHwpfP4lfJ56HK+3X7fmHdiul5GO77gxaLZuxtQ6hPxP\n5mDo07fO3k9jIFYGlSHEQPAf46l7E1AqrA1yLhSy2VTxevVFFDodhS+9hsXHt9rx2hVL8fj2a4yx\nHSl85oUqxxXfcx/nItqh/eE73pr0OX9sOgLffUezAVei2bqZ4qG3cn71xkteCECsDJy5jCsSCgQX\n4uquoXQI9WtsM2qFZt0a3Jf8gqF7T/R33HXBsRmZ+Sz48l9emfMEZo2WvE/mgNsFku5UKt7vPoK3\njr7M48tnkb7zd9i/BoWnF3kzPkU//J7L5r4hxMARsWcgEFxaGI14T3kGi0JBwTvTnMtHOKDIOY/b\nH0vwmPkl7x3ZgxILC699iMTYajbHLRbalrbajMo6SlTWUejZk3MzZmOOalvHb6ZxEWLgiBADgeCS\nwuOrOdbexPeNxNi5q/PJwkLclv2F2y8/oV35DwqDgVhgfyuJlXGJbL9iCBeqqvT/9u48TKryyuP4\nt5peAGloCK2CGME4OYgoCCqLLJFRFHcZUAwoiqJhRMc1iWJcUBQFHBcIjgYTI2oUArIoDCAKUQfX\nqLgdMQEFNYJKL0CvdM8f762mKHqp/Rbd5/M8/TR1+63uH1Vv1am7nRsoKiR34pWMW/tSzbIPRo6n\nx9OzqCpouHvr/saKQSgrBsbsN0IvWrPz5tvcwvJyslevImfhPHL+dxkB7xyBym7dKR0+go0DTmPG\nugJg386soZp9tJ62pwwisNtd5Kc6O5vCeYvo2O9EyMoCrBgYY0xaCF60ZsfdU8n89GNyFswjZ+li\nMgrdm/3uzl0oHT6CsvNGslu6AnAwMK1X/b835y9P0/qaCTW3y848h+IHHqY6r22y/itpwYpBKFsz\nMCbtjZu6miP+9QUznnGtols+OJ2M791lMncf3IFdo0ZTNnwElT17Rbdzt7SU3IlX0nzxwppFjW0n\ncX2sGISyYmDMfmH0G8+QEewuvHs3JRddStnwEVT07Q/NmkX9+zI2beQnJ/SouV3Vvj3bX1xFVZfD\nExU57VkxCFVT/a0YGJPO3j+sJ4Ut2/D6z/szYeZ1kJ0d8+/KfnEJbS4dXXN711X/xc5bbvP2DTQd\nVgxqY2sGxqS1Rb3Prvn3hFgLQWUluRMup/miBTWLChYvd2sXTZAVg1BNYLugMfu7SC6z2ZDA1q20\n776nh1CldKVg6Qqq2+TF/bv3VxEVAxHJAKYAY4FcYDlwlapurWP888AI3PaW4DvsKlUdGnfiJLLr\nGRjT+GWvWEabMRfU3C6+dzql48Y3+Q+Dka4Z3AlcBIwBfgRmA/OBQXWM7w78GvhzyLKyGDOmTs0O\nZH9jGGOSoLqa1mN/Sc7yF2sW/fDWB1R17uJjqPTRYDEQkSzgGmCiqq72lo0CNopIX1VdFzY+GzgC\neLuuNYe0ZUcTGdMoBX78gfZd97zplw8YROHzL0CmbSkPiuSR6Am0AtYEF6jqlyKyCRgIrAsb3xVo\nBnyamIipZy2sjdn/3TP3Hb7YUkTfDeuYtGRqzfLCJ+ZSfubZ9dyzaYqkGHTyvn8dtvwb4NBaxncH\nKoDJIjIMKAHmAXeranpvKmri2wyNaUy+2FLEhFWPcvqHy2uWfa+b9qtLUaZSJNczaAlUqerusOVl\nQG1Xmw62AfwEOB24A7gceDTGjKljm4mMaVSO/+fbAKyVAWzbWmSFoB6RrBmUABkikqGqVSHLc4Cd\n4YNVdZKITFPVAm/RxyJSBTwrIter6vb6/lh+fnwX2o5LRjkAOdnN9srha6Y6pGMmSM9clikyjS3T\nUZ3bcv3o6VQHAnTq1oVBCfz/peNjFa9IisFm73sH9t5U1JF9Nx0BEFIIgtZ73w8F6i0G27YVRxAp\nOQLbd9AeKCutoMjLkZ+f62um2qRjJkjPXJYpMo0x0w2jjgX2tLVO1P8vXR+reEWymegDYAcwOLhA\nRDoDnYG14YNF5DkRWRC2+HjcZqUvYg1qjDEmeRpcM1DVchH5PTBdRH4AtgGzgFdU9S3v0NN2wI+q\nWoE7/+BZEbkOWAT0AqYB01R1V7L+Iwlh+wyMMU1UJGsGALcCTwNPAS8DG4GR3s/6444s6gegqvOA\nS7yv9bhC8N+qenuiQieNFQNjTBMV0RkX3pFEN3lf4T9bgzuvIHTZXGBuIgKmlHUtNcY0UZGuGTQt\ntmZgjGlirBiEspPOjDFNlBWDULbPwBjTRFkxCGEtrI0xTZUVgxDveE1W13/xPW9+8i+f0xhjTOpY\nMQjxx2WfuX9UV/P4kk/8DWOMMSlkxSBEZUYmO3Ja0mPzh5y8fqXfcYwxJmWsGIQYd+4x3H/mrynN\nas7VK2bR6qbroLzc71jGGJN0VgxC9Ol2MFfPvp7yta9T2a07LZ6cAyedRMZ3tv/AGNO4WTGoRVXn\nLmx/cSWlw0fAG2+Qd/IgMt960+9YxhiTNFYM6nLAARTPngMzZpCxbSt5551O8yef8DuVMcYkhRWD\n+gQCcP31FD7/AtW5ueTedC2tbrgGytL76p3GGBMtKwYRqBj0C7avWENF92No8dSfyDt3GBnffuN3\nLGOMSRgrBhGq+ulhFCxdQemIC8h69x3anjyIzHX/53csY4xJCCsG0WjZkuJZj7Fjyn0EfvyBvOFn\n0HzOY9a+whiz37NiEK1AgJLxEyicv5jqvDxyb76RVtdeBaWlficzxpiYWTGIUcWJA91+hB7H0uLZ\nueSdfSoZX2/xO5YxxsTEikEcqjodSsHi5ZRe8Euy3v87bU8ZRNYbr/kdyxhjombFIF4tWlD88GyK\n751OoKCANv9xFi0en237EYwx+xUrBokQCFB62RUULlhKddt2tJr0G3InXgklJX4nM8aYiFgxSKCK\nvv3ZvmotFb1603zeX8g761QyNn/ldyxjjGmQFYMEq+p4CAWLllMy+mKyPnzf7Uf42xq/YxljTL2s\nGCRDTg47HniE4mkPEigups3Ic2gxe6btRzDGpC0rBskSCFA6dhwFC1+iqn0+rW6/hdwJl8GuXX4n\nM8aYfVgxSLLKE/pQsGotFcedQPMF82l7xilkfLnJ71jGGLMXKwYpUHVwBwoWvkjJxePI/Hg9bYcO\nJuvV1X7HMsaYGlYMUiUnhx3TH6T4gUcI7NxJm1HDaTHzIduPYIxJC1YMUqx0zFgKFi2j6sCDaDX5\nd+RecSns3Ol3LGNME2fFwAeVvY9n+8q1VPTpR/NFC2h7+slkbPyn37GMMU2YFQOfVB90EAV/XULJ\nuPFkfvoxbYf+gqzVK/2OZYxpoqwY+Ck7mx1TZ1D00O8JlJbQ5sIRtHhohu1HMMaknBWDNFB24RgK\nFi+nqkNHWk25k9aXXUxgR7HfsYwxTUhmJINEJAOYAowFcoHlwFWqurWO8ccBDwLHAluAu1X1qYQk\nbqQqj+3N9pVraT1+LDlLF9Fsg1L05DPsPvwIv6MZY5qASNcM7gQuAsYAA4FOwPzaBopIe1yxeAdX\nDB4B5ojIyXGnbeSq8/MpnLeIXVdMIFM/I2/oSWSvXO53LGNME9BgMRCRLOAa4GZVXa2q7wOjgAEi\n0reWu4wHClT1WlX9XFVnAnOBGxMZvNHKymLn3fdRNPN/CJSX0XrMBbR84H6oqvI7mTGmEYtkzaAn\n0Aqoab2pql8Cm3BrCeEGAGvDlr0KnBhTwiaq7PwLKVi6gqpDOnHA1LtpfekYAsVFfscyxjRSkRSD\nTt73r8OWfwMcWsf42sa2FJF20cVr2iqP6cn2FWsoHziYnGVLyTttCM02fO53LGNMIxRJMWgJVKnq\n7rDlZUDzOsaX1jKWOsabelS3b0/hcwvZ9auJZG74nOYnDeSucyfx4Pz3/Y5mjGlEIikGJUCGd0RR\nqBygtj4KJd7PwsdSx3jTkMxMdk6+h+nDrqNZ1W5+t+gedr7+lt+pjDGNSCSHlm72vndg780/Hdl3\nc1BwfIewZR2BHapa2MDfCuTn50YQKbXSJdOaIwdXrzlycM3t/PzcgI9xapUuj1UoyxQZyxS5dM0V\nj0iKwQfADmAw8AyAiHQGOrPvjmKA14BLwpYNAV6PMaPxLJlxTtq9+RtjGodAdQStD0TkXtwJZ5cC\n24BZwC5V/Xfv0NN2wI+qWiEiBwKfAc8BDwGnANOAU1XVLgZsjDFpKNKTzm4FngaeAl4GNgIjvZ/1\nxx0t1A/AOyv5NNwJZ+8B/wlcZIXAGGPSV0RrBsYYYxo3a1RnjDHGioExxpgIu5Ymg4jk43ZEnwKU\nA38EblHVBpvwiEge8CHwB1Wd7FcmEckEJuGa+B2M23E+WVUXx5kj7brExpDpAuC3wL/h9inNAaZF\n8vwmK1PYfZcCLVV1SKLyxJJJRA7BHWgxFHeOznzgBlUNP3EzlZmGAPcCRwHfAo+p6rRE5anjbz4K\nZKjqFfWMSWk35AgzJX2eR5spbHzE89zPNYMFwIG4/kbBI5XujPC+s4FD0iDTFOAKXCO/Y3Av5AUi\nMiDOHOnYJTaaTMNwzQkfA47GvVh+A9zsV6awfFcCpyc4S9SZRCQbWAXk4Q7AOB84E7jfx0w/A5YA\ni4HuuOftdhGZkOBMoX9zMu51VN+YlHZDjjBTquZ5xJnCxkc1z31ZMxCRfrijkLqo6lfARyJyE/Cw\niExW1Yp67nsh0IvaT3hLWSYRCQCX47q5vuQtnupNzktw51vEkiPYJXaiqq72lo0CNopIX1VdF3aX\nmi6x3u3PRaQXrkvsqlgyJCDTlcA8VZ3t3d4oIt1wxXWKT5mC9zvCy/BGInLEmWk0cBDQR1WLvPG3\nAQl7440h02m4w8aDz9Mm79PvqbgPYQkjIl1wn6SPAr5sYHjS53kMmZI+z2PIFLxP1PPcrzWDAcCX\n3ptu0KtAa1yX1FqFrFJfzJ5+R35lysAdXrswbHkV0DaOHOnYJTbaTHcB4ZvvqonvcYk3U3BzyZPA\nVODTBGaJNdNQYGWwEHjjn1TV2lrDpyrTNqCdiIwSkYCIdAcGAW8nMFNQf+Ar3KfqTQ2MTVU35Ggy\npWKeR5sp5nnu1z6DujqbguuEWtfEewJ4XFXfFBFfM3mN+1aHLhOR43FnW/8qzhzUkaWuLrHv1TK2\npYi0U9Uf48gSUyZVfTf0toi0xj0myxKQJaZMnltwTReni8jjCcwSa6afAy97q/9jcG8kC4BbVTVR\nH3aizfRX3OsseF5RM+C5kDWFhFHVp72/QwSv51TM86gypWieR/s4QYzzPCnFQEQOw52YVg2Et1Ao\nxW1n22sHmapWikg1dXQ2FZFrcKvUt6VLprDffwTuhbwOt+M5VunYJTbaTDVEpAXwgjcukdtSo8ok\nIr2B64DjEpghrky4tc7LgZeAEbj9YLOAfPZt6ZKqTHm4VjNTgedxn0YfEpE7VPWOBGWKRVp3Q07i\nPI82R8zzPFmbib4GugJHet9Dv47BPal7dTb1jswJUEtnUxHpilsdu7iWSe1LprBxvYG/Ad8DZ8WR\nEdKzS2y0mQAQkZ/gzljviWtHsrmuscnMJCI5wJ9xn7g3JjBDzJk8FcAPuDP031PVJbgX8kUikqhN\nDdFmuh+oUNVJqvqBqgavUvjbBGaKRdp2Q07yPI8mR1zzPClrBqpaCdR5FRYR2QwMC1vc0fte247h\n84EDgNe8HbfgPilMEpGRqnq0D5mC9xuKOzLj78DZEXRmbUgqu8QmK1OwmeEK3PM2UFU/TlCWWDL1\nwRX9+0QkeKRODu5NsgjopqpbUpwJb1mJqoa2AfgE9wGkM7Ddh0x9cGu4od4EsoGfJihTLFIxz6OW\ngnkejbjmuV87kF8DDvd2CAcNAYqA2q7a8jAguMrbw/vajDu6IVGHCEabCREZCCzC7TsYmqBJGdol\nNvh3OlN/l9hBYcsS3SU2qkze+Rqv4DbJ9UvSCySaTG/ijgMPnT8LcfuBerBn31AqM4Fbm+wpIs1C\nlh0NVBLBjsIkZdqCW1MOdTSwG/hHgjLFIhXzPCopmufRiGue+9abSERexz2IV+NO2PoTMFNV7/J+\nfgDQSlW/q+P+G4CnEnzSWcSZvGPE/4FbzT8D9wIOKlPVgjhypF2X2CgzzcMdKTME9+YSVB3JCWHJ\nyFTLfR8HfpaEk86ife4+AlbiNoMeCvwBd4TReJ8yDcOdZ3A7rmX9UcCjwPyQwzoTTkReATYET6by\na55HmSkl8zyaTLWMj3ie+3nS2XnAd7hPJ3NwZzneFfLzG6m/kiWjikWTaTBuNfVo3GFf34R8PR9n\njnTsEhtRJhFpjnscWwFvsecx+Za9XzApy5Tgv5mwTN5zNwj3Yn4XdxDDPNxz6FemZcBw4BzcWsUD\nuGJwQ4IzhQt/PadDN+Q6M6V4nkeUKd5fbF1LjTHGWKM6Y4wxVgyMMcZgxcAYYwxWDIwxxmDFwBhj\nDFYMjDHGYMXAGGMMVgyMMcZgxcAYYwzw/3J2xE9/VAPuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f47aef0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(metropolis_samples[:,0], metropolis_samples[:,2], \".\")\n",
    "plt.plot(unique_samples_met[0:50,0], unique_samples_met[0:50,2], 'r')\n",
    "pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# test speed up ratio for 50-d\n",
    "# 50 dimensional multivariate normal data\n",
    "cov = np.identity(50)\n",
    "X_MVN = np.random.multivariate_normal(np.ones(50), cov, 100)\n",
    "\n",
    "def L_MVN(theta):\n",
    "    cumsum = 0\n",
    "    grad = 0\n",
    "    for x in X_MVN:\n",
    "        cumsum += ((x - theta).T @ np.linalg.inv(cov)) @ (x - theta)\n",
    "        grad += (x - theta).T @ np.linalg.inv(cov) \n",
    "    logp = -0.5 * cumsum\n",
    "    return logp, grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3min 58s, sys: 2min 32s, total: 6min 31s\n",
      "Wall time: 4min 47s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "samples_MVN, burned_in = nuts6_dual_averaging(np.zeros(50), 5000, 1000, L_MVN)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}