Effective N testing

EffN testing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     },
     "base_uri": "https://localhost:8080/",
     "height": 260,
     "output_extras": [
      {
       "item_id": 2
      }
     ]
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 2641,
     "status": "ok",
     "timestamp": 1517625919949,
     "user": {
      "displayName": "Sharan Yalburgi",
      "photoUrl": "//lh4.googleusercontent.com/-ypYfwsB3KsY/AAAAAAAAAAI/AAAAAAAAD0s/JStFE5QSLnE/s50-c-k-no/photo.jpg",
      "userId": "108333366686678045862"
     },
     "user_tz": -330
    },
    "id": "DJ8Nsu6K6IgB",
    "outputId": "f36259dd-91c1-4d36-fb65-a09371fb2d2b"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: pystan in /home/sharan/anaconda3/lib/python3.6/site-packages\r\n",
      "Requirement already satisfied: numpy>=1.7 in /home/sharan/anaconda3/lib/python3.6/site-packages (from pystan)\r\n",
      "Requirement already satisfied: Cython!=0.25.1,>=0.22 in /home/sharan/anaconda3/lib/python3.6/site-packages (from pystan)\r\n"
     ]
    }
   ],
   "source": [
    "# !pip install pystan "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     },
     "base_uri": "https://localhost:8080/",
     "height": 34,
     "output_extras": [
      {
       "item_id": 1
      }
     ]
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1243,
     "status": "ok",
     "timestamp": 1517625921261,
     "user": {
      "displayName": "Sharan Yalburgi",
      "photoUrl": "//lh4.googleusercontent.com/-ypYfwsB3KsY/AAAAAAAAAAI/AAAAAAAAD0s/JStFE5QSLnE/s50-c-k-no/photo.jpg",
      "userId": "108333366686678045862"
     },
     "user_tz": -330
    },
    "id": "0Ol7leKkAKbZ",
    "outputId": "228f8a49-e38d-4c78-900c-7f94af757123"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "import pystan\n",
    "import pystan.chains\n",
    "\n",
    "import numpy as np\n",
    "from collections import OrderedDict\n",
    "\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     },
     "base_uri": "https://localhost:8080/",
     "height": 322724,
     "output_extras": [
      {
       "item_id": 1
      }
     ]
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 4954,
     "status": "ok",
     "timestamp": 1517626121793,
     "user": {
      "displayName": "Sharan Yalburgi",
      "photoUrl": "//lh4.googleusercontent.com/-ypYfwsB3KsY/AAAAAAAAAAI/AAAAAAAAD0s/JStFE5QSLnE/s50-c-k-no/photo.jpg",
      "userId": "108333366686678045862"
     },
     "user_tz": -330
    },
    "id": "62cZLYrKAUG3",
    "outputId": "9146a963-b878-4cf1-90f2-dc29f3c80f35"
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "f1 = 'file1.csv'\n",
    "f2 = 'file2.csv'\n",
    "# import urllib.request\n",
    "\n",
    "# url1 = \"https://raw.githubusercontent.com/stan-dev/pystan/develop/pystan/tests/data/blocker.1.csv\"\n",
    "# url2 = \"https://raw.githubusercontent.com/stan-dev/pystan/develop/pystan/tests/data/blocker.2.csv\"\n",
    "\n",
    "# c1 = urllib.request.urlopen(url1).read()\n",
    "# c2 = urllib.request.urlopen(url2).read()\n",
    "\n",
    "# with open('file1.csv', 'wb') as fx1: # bytes, hence mode 'wb'\n",
    "# # fx1 =  open('file1.csv', 'r+b')\n",
    "#   fx1.write(c1)\n",
    "    \n",
    "# with open('file2.csv', 'wb') as fx2: # bytes, hence mode 'wb'\n",
    "# # fx2 =  open('file2.csv', 'r+b')\n",
    "#   fx2.write(c2)\n",
    "\n",
    "    \n",
    "# read csv using numpy\n",
    "c1 = np.loadtxt(f1, skiprows=41, delimiter=',')[:, 4:]\n",
    "c1_colnames = open(f1, 'r').readlines()[36].strip().split(',')[4:]\n",
    "np.testing.assert_equal(c1_colnames[0], 'd')\n",
    "c2 = np.loadtxt(f2, skiprows=41, delimiter=',')[:, 4:]\n",
    "c2_colnames = open(f2, 'r').readlines()[36].strip().split(',')[4:]\n",
    "np.testing.assert_equal(c1_colnames, c2_colnames)\n",
    "np.testing.assert_equal(len(c1_colnames), c1.shape[1])\n",
    "\n",
    "n_samples = len(c1)\n",
    "np.testing.assert_equal(n_samples, 1000)\n",
    "\n",
    "c1 = OrderedDict((k, v) for k, v in zip(c1_colnames, c1.T))\n",
    "c2 = OrderedDict((k, v) for k, v in zip(c2_colnames, c2.T))\n",
    "\n",
    "lst = dict(fnames_oi=c1_colnames, samples=[{'chains': c1}, {'chains': c2}],\n",
    "               n_save=np.repeat(n_samples, 2), permutation=None,\n",
    "               warmup=0, warmup2=[0, 0], chains=2, n_flatnames=len(c1))\n",
    "\n",
    "# lst['samples'][0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lst['samples'][1]['chains']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     }
    },
    "colab_type": "code",
    "id": "4JymrftZ8gWO"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from pymc3.stats import statfunc, autocorr\n",
    "from pymc3.util import get_default_varnames\n",
    "from pymc3.backends.base import MultiTrace\n",
    "\n",
    "def effective_n(mtrace, varnames=None, include_transformed=False):\n",
    "    R\"\"\"Returns estimate of the effective sample size of a set of traces.\n",
    "    Parameters\n",
    "    ----------\n",
    "    mtrace : MultiTrace or trace object\n",
    "      A MultiTrace object containing parallel traces (minimum 2)\n",
    "      of one or more stochastic parameters.\n",
    "    varnames : list\n",
    "      Names of variables to include in the effective_n report\n",
    "    include_transformed : bool\n",
    "      Flag for reporting automatically transformed variables in addition\n",
    "      to original variables (defaults to False).\n",
    "    Returns\n",
    "    -------\n",
    "    n_eff : dictionary of floats (MultiTrace) or float (trace object)\n",
    "        Return the effective sample size, :math:`\\hat{n}_{eff}`\n",
    "    Notes\n",
    "    -----\n",
    "    The diagnostic is computed by:\n",
    "    .. math:: \\hat{n}_{eff} = \\frac{mn}{1 + 2 \\sum_{t=1}^T \\hat{\\rho}_t}\n",
    "    where :math:`\\hat{\\rho}_t` is the estimated autocorrelation at lag t, and T\n",
    "    is the first odd positive integer for which the sum\n",
    "    :math:`\\hat{\\rho}_{T+1} + \\hat{\\rho}_{T+1}` is negative.\n",
    "    References\n",
    "    ----------\n",
    "    Gelman et al. (2014)\"\"\"\n",
    "\n",
    "    def get_vhat(x):\n",
    "        # Chain samples are second to last dim (-2)\n",
    "        num_samples = x.shape[-2]\n",
    "\n",
    "        # Calculate between-chain variance\n",
    "        B = num_samples * np.var(np.mean(x, axis=-2), axis=-1, ddof=1)\n",
    "\n",
    "        # Calculate within-chain variance\n",
    "        W = np.mean(np.var(x, axis=-2, ddof=1), axis=-1)\n",
    "\n",
    "        # Estimate marginal posterior variance\n",
    "        Vhat = W * (num_samples - 1) / num_samples + B / num_samples\n",
    "        print(\"vhatnew \", Vhat)\n",
    "        return Vhat\n",
    "\n",
    "    def get_neff(x, Vhat):\n",
    "        # Number of chains is last dim (-1)\n",
    "        num_chains = x.shape[-1]\n",
    "\n",
    "        # Chain samples are second to last dim (-2)\n",
    "        num_samples = x.shape[-2]\n",
    "\n",
    "        # Calculate within-chain variance\n",
    "        W = np.mean(np.var(x, axis=-2, ddof=1), axis=-1)\n",
    "\n",
    "        rho = np.ones(2 * num_chains + 1)\n",
    "        negative_autocorr = False\n",
    "        # Iterate over different lags of autocorrelation\n",
    "        for t in range(1, 2 * num_chains + 2):\n",
    "            if negative_autocorr: \n",
    "                print(\"break\")\n",
    "                break\n",
    "            auto_corr = []\n",
    "            for m in range(num_chains):\n",
    "                auto_corr.append(autocorr(x[:, m], t))\n",
    "            print(\"acorrmean \", np.mean(auto_corr))\n",
    "            rho[t - 1] = 1. - (W - np.mean(auto_corr)) / Vhat\n",
    "            negative_autocorr = sum(rho[t - 2:t]) < 0\n",
    "        print(\"rho \",rho)\n",
    "        tHat = 1. + 2 * rho.sum()\n",
    "        neff = num_chains * num_samples / tHat\n",
    "        return abs(neff)\n",
    "\n",
    "    def generate_neff(trace_values):\n",
    "        x = np.array(trace_values)\n",
    "        shape = x.shape\n",
    "\n",
    "        # Make sure to handle scalars correctly, adding extra dimensions if\n",
    "        # needed. We could use np.squeeze here, but we don't want to squeeze\n",
    "        # out dummy dimensions that a user inputs.\n",
    "        if len(shape) == 2:\n",
    "            x = np.atleast_3d(trace_values)\n",
    "\n",
    "        # Transpose all dimensions, which makes the loop below\n",
    "        # easier by moving the axes of the variable to the front instead\n",
    "        # of the chain and sample axes.\n",
    "        x = x.transpose()\n",
    "\n",
    "        Vhat = get_vhat(x)\n",
    "\n",
    "        # Get an array the same shape as the var\n",
    "        _n_eff = np.zeros(x.shape[:-2])\n",
    "\n",
    "        # Iterate over tuples of indices of the shape of var\n",
    "        for tup in np.ndindex(*list(x.shape[:-2])):\n",
    "            _n_eff[tup] = get_neff(x[tup], Vhat[tup])\n",
    "\n",
    "        if len(shape) == 2:\n",
    "            return _n_eff[0]\n",
    "\n",
    "        return np.transpose(_n_eff)\n",
    "\n",
    "    if not isinstance(mtrace, MultiTrace):\n",
    "        # Return neff for non-multitrace array\n",
    "        return generate_neff(mtrace)\n",
    "\n",
    "    if mtrace.nchains < 2:\n",
    "        raise ValueError(\n",
    "            'Calculation of effective sample size requires multiple chains '\n",
    "            'of the same length.')\n",
    "\n",
    "    if varnames is None:\n",
    "        varnames = get_default_varnames(mtrace.varnames,include_transformed=include_transformed)\n",
    "\n",
    "    n_eff = {}\n",
    "\n",
    "    for var in varnames:\n",
    "        n_eff[var] = generate_neff(mtrace.get_values(var, combine=False))\n",
    "\n",
    "    return n_eff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "u1W3GE9RJKTU"
   },
   "source": [
    "## For Scalar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     },
     "base_uri": "https://localhost:8080/",
     "height": 317,
     "output_extras": [
      {
       "item_id": 1
      },
      {
       "item_id": 2
      },
      {
       "item_id": 3
      }
     ]
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1512,
     "status": "ok",
     "timestamp": 1517626474877,
     "user": {
      "displayName": "Sharan Yalburgi",
      "photoUrl": "//lh4.googleusercontent.com/-ypYfwsB3KsY/AAAAAAAAAAI/AAAAAAAAD0s/JStFE5QSLnE/s50-c-k-no/photo.jpg",
      "userId": "108333366686678045862"
     },
     "user_tz": -330
    },
    "id": "yuaPK9GzGi_D",
    "outputId": "e0852c1f-3013-43bd-d654-d85592f353b8"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsXXec3MTZfkbS7t6dewMbFwyYDqaZXkMLoYcvoSYhEEoKIY1USnqAhAAhJAQCBNNDCWA6BoMBU21sinvF3T6Xu7Ov7Uqa7w9ppJnRjKTdW5/PPj2/n723q9FoNBo9884z77xDKKXIkCFDhgzdC8bmLkCGDBkyZOh8ZOSfIUOGDN0QGflnyJAhQzdERv4ZMmTI0A2RkX+GDBkydENk5J8hQ4YM3RAZ+WfIkCFDN0RG/hkyZMjQDZGRf4YMGTJ0Q1ibuwA6DBw4kI4cOXJzFyNDhgwZtihMmTJlDaV0UFK6Lkv+I0eOxOTJkzd3MTJkyJBhiwIh5PM06TLZJ0OGDBm6ITLyz5AhQ4ZuiIz8M2TIkKEbIiP/DBkyZOiGyMg/Q4YMGbohMvLPkCFDhm6IjPwzZMiQoRsiI/8MGTJ0PmY9D2xYtblL0a2RkX+GDBk6F04JePR8YOypm7sk3RoZ+WfIkKFzQV3vc/2izVqM7o6M/DNkyJChGyIj/wwZMmTohsjIP0OGDJ0LSjd3CTKgSuRPCDmJEDKbEDKPEPILxfECIeS//vH3CSEjq3HdDBkyZMhQGTpM/oQQE8A/AHwJwB4AziOE7CEl+xaA9ZTSUQBuAXBjR6+bIUOGDBkqRzUs/4MAzKOULqCUFgE8CuAMKc0ZAMb6fz8B4DhCCKnCtTNkyLDFIZN9ugKqQf5DASzhvi/1f1OmoZTaABoBDJAzIoRcRgiZTAiZXF9fX4WiZciQIUMGFapB/ioLXu7a06QBpfQuSukYSumYQYMSdyHLkCHDlohswrdLoBrkvxTAcO77MADLdWkIIRaAPgDWVeHaGTJk2OLAyD9TfjcnqkH+HwLYmRCyAyEkD+BcAOOkNOMAXOj//RUAEyjNuv8MGbo3MgrYnOjwBu6UUpsQcgWAlwGYAO6llE4nhPwOwGRK6TgA9wB4gBAyD57Ff25Hr5shQ4YtFJnd1yXQYfIHAErpCwBekH67jvu7DcBXq3GtDBkybOnIZJ+ugGyFb4YMGTJ0Q2TknyFDhs5FJvt0CWTknyFDhk5GRv5dARn5Z8iQIUM3REb+GTJk6Fxksk+XQEb+GTJk6GRk5N8VkJF/hgwZMnRDZOSfIUOGzkUm+3QJZOSfIUOGDN0QGflnyJAhQzdERv4ZMmToXGSyT5dARv4ZMmToZPjkn23mt1mRkX+GDBk2D7IRwGZFRv4ZMmToXGSk3yWQkX+GDBk6GZns0xWQkX+GDBkydENk5J8hQ4bORSb7dAlk5J8hQzdHyXHx3oK1nXjFjPx1+PFj0/DcJ8s75VoZ+WfI0M1xxxvzce5d7+Gd+Ws2d1G6Pf730TJc8fDUTrlWRv4ZMnRzrNnYDgCYuWJD51wwk326BDLyz5Chm2Pb3jUAgFVNbZ10xYz8uwIy8s+QoZujT20OALChzd7MJcnQmcjIP0OGbg5mh3ea230m+3QJZOSfIUN3h0/GnbfkKiP/roCM/DNk6ObodMs/gxK0k0dEGflnyNDNwTiHdJbtn8k+SridXC0Z+WfI0M3BLE6d5W87LpY1tFbziv5nNtTg4XQy+2fknyFDN0cSFf/h+Zk4/IYJwXqA6l85AwC4meyTIUOGzQGiMf3fnFMPAGhoKVXnQpnso0RG/hkydEGsay7i8BsmYPbKTloF24nofC7OZB8VtijNnxDSnxAynhAy1//sp0n3EiGkgRDyXEeulyHD5sJrM1dhWUMr7nxz/uYuStWR1tsn8wbatNjSNP9fAHiNUrozgNf87yr8BcDXO3itDBk6jLaSgxWNlU9edppHTLUw5xXg/jNizfu0LoZVGyFkso8SW5qr5xkAxvp/jwVwpioRpfQ1AFvfeDnDFofvPzIVh14/oewXbYv1hf/vBcCCNwA7ebJW27FV/Z4z8ldhi5J9AGxLKV0BAP7nNh0v0ubBmo3tVXZny9AVMX7GKgAVDLEDX/gtDazEcZa/n3LLu7mtCl1O9iGEvEoI+Uzx74xqF4YQchkhZDIhZHJ9fX21s4/FmD+8isNvmNCp19xccF2KiXPqO32Y2RVgGh7DlZxyLf8ttK4Yo8fJPuiE8A4bVgKvXAO4Tib7aNDlZB9K6fGU0r0U/54BsIoQMgQA/M/VHSkMpfQuSukYSumYQYMGdSSrDDF4+IPFuPDeDzDu487ZMagrwfTJsOS6ZZ235VrHXcTyf+7HwDt/Bxa8vgkvsmXD4cj/8gcmb/LrdVT2GQfgQv/vCwE808H8MnQClqxvAQAsb+is+O1dB8zyt8u0/Bm2uAnfwPLXd3bhfEYn3Fspk1Z14FWf1RuqvaAuio6S/w0ATiCEzAVwgv8dhJAxhJC7WSJCyFsAHgdwHCFkKSHkix28bga7HXj+KqBlXfnnVnF0+dwny7HU70w2BVyX4v0q7i8byj5lWv5VK0FnIwX5p57P6EAtWHnv027femWfWS8A0x6u+HSXY//OMDE6RP6U0rWU0uMopTv7n+v83ydTSi/h0h1JKR1EKa2llA6jlL7c0YJ3e3z6OPDhv4HXfrtZi3HFw1Nxxu2TNln+d721AOfc9V6wyrSj8Lm/fPLfUmWfVJZ/J8xmW95uYZ7XkYL8m1YAv+kDfPbkJizEJsaj5wFPfyf6+/rPgbtPAFrXx57Or/DtjFHY1rvC9/GLPI2xDPzIehy4eU/gg3+nO2H9os1nxbDrOpXvvlSt9rW2uVidjBRYUL8RALC8oRVoXgusmt6h/CzTa/Llyj7BpOgWR/7+Kx4zx9EpUT2tgvdpa6TG1f5z/eiBTVeGzYW3/gos/QCY/lRssi7n7bPFYvr/PO+CMvAFYxrQtBRYODE58dLJwN/2ASbfW2EBOwj2UsdYdFsDmEzjUgB3HQ3ccViH8jNIZbJPiC2O/b2PWNmnE0jH9MnfKWoMpi2tXstBunAWPPd3xjPZesm/AtTBn2RJMym11l/mv/jdTVegOKQYzm9qdEYDZcNfh1KgcUmH87MqdfXcUmXqwNlH304Y6RgabqoKLQuWf5yGtqVWdAyCoVU83fKyT2fUQrcj/+Z2Gy1FtVRSSxj5p/CCMUzv03WqVLIy0QUs/04xGH2CcKs0JK5ownfVDBw+7Wcw4Wx5sk8Ky5+RTtK9HX/zm5UXQ9D8fWyxPWqZSDlhxJN/ZyhA3Y789/z1y9j3d+PDH5pWAOsWAKCc5Z/CeyUg381F/n7nsxnJvzNC0IayT3XJ3y7Hz/+py7DDypewG1my5YkTaSZ8O4ODzQRvnzS96vt3Ao+cV91ydQpSyj7cI8pkn02Eou3Xcut64ObdgNv2w6HGDNTBt/hLrcDi94DWBn0mm93yr1z2qVp8rirlEwd2m9WaDBNW+DavAX7TF1iU3ltpa7T8O0XzNxjVUMTKPnFlefFnwOwXql2yTQ/e8l8+DRh/HVBsjiRzndBxojMMq62f/Kfcpz/GkfteZCEKxJeDWtcD934RePQCfWM0LO9zCyT/IIsOFqEzZZ9qXUuQfZZ8AICW5RW2NS/ySlPHFXcUgdeR7n3Z/PVaclwsqN+4CTpDTvN/9dfApL8By6d6HQCbO1y3AHvdsxPONN4GEOucVTVs/eT/7A/0x7iGeHWOW5zRssb7/PxtvXdJQP6Vu1p2CF1A8+8U68TgJnyrANaZ2A7tEpPmmx7pNX9dDfM+5xU/Br69xmay+eYB/vzSLBz714mYtiRmxF8JgronQONS70/XAR45F/j7/t731TMBAKeY73mHM8t/E4MR9+hzcKt9Fv5SOhvOHmeJhL56hvrcQHPvvhO+nQHmmsnLPh2Jx28IE77JcW9kdGnZZ8OqMHgaQxU0f1oVLxS+4mI0/804CfzpskYAwJqNVV63wss+NsubAgu5CXSfT4xO7Py6Ofn7e5Ludiputb+Cfzhnwt7+yHTnGknD2E2MoPPZfC+L7tLVXKxiBJwQ5nnKbW9XnB+jIJdSrgMtg/wrvnIn4NkrPQlLWKeSxvL3k5Qh+xRtF22lctq+n/kmkH0cl+KuN+ejtdixdzFvee9UefeVBpzs4/hOJVSY3Q3aognv98zy39RgFj6TcAAUdzw+3bmb29umA5Z/tTRNVZjjlz5biZ1+9UIVwzF4pLChPRyNrevAimLWZ3uGf7myTxd3TWSrZ/nnm+Ie735rgZdEc3+C7ON/nvmPSTjixjIidLLr62SfDgypxs9YhT+9MAs3vDiz4jwAIO+v/q46+fOyD3N15euA0qBhGgH5V7cIKnRz8vcfspkLfppcn7JKWGPdbJq/4qW224H2TbRh2pyXgfniy65qoFMXe/FLZq5oqsplmUxz58QFVckv0PxdlyOcpDctJKZOiXxZKZT+5Mnkb5fBNOwSM1Y0Yc3GhMiTzVxAPnYidRFb3xUYJj0KniE2e1XH2n7B2lTkzz0Xh5N9wgSc7JNZ/p0Dx5d9mNsmgIvu/zj2lPn1G3HSrW9iY5v/EDtZ9lm0phmTF62DUq/+97HA9cM2zYUfPht4QNylUzWCYI3W1C0XjcPsl4D7ThUIQJ1NB14MYQ4hnhjrN7RXNaLopofCn7yc0U0a2Sdt3V8/HPjLjsDKz8Tra8tReada8OWaVU3twFs3ewHi7PJHhyH5V3s0z9WZ0vJ3OdmHRg5vKljJSbZiKGSfJPxjwjzMWrkBUz9vxpFAp0/4HnPTGwCARReJcw4/fmwabl71WVl5dbR9qc5nC2crspAfPd+rT9cORmOmIp+OTIqxzsRxeW8fdX6n3PYWVm9ox6KR4W9d2fCvxPIvdzJXrirXpcHoTEC7P/JbOxcYvJdo+XfU24dS4R7ZHFPRdoGJf/Z/bA/DSKdEfpNb/kbIF3wdLHk/MEAN4vqHM8t/0yIg/1D2STRYg3fJf5k214SvZEn976NlZWfR0YlZFZ8wy78Swz/MmLP8FRmRKvicCJa/Jr9zWx7B6cY70vldmf0Z0lv+iTGOJv8HuzpztIcTd0QLDKsE2aecXlV651ibc1wazntUIMeyttZmbyLNn79H/nncdwow7zUAwMHGLOxJFlbNtTkO3Zb8563eEDYizvLvUYgfBRjBA2Tku7nIn3+ZKju1w+SveInZi1hZ1tF7Mqpu+XPrBhK8fX6cewK35W+v+Fqdjgosfz7MhdLafO6H+EfzVZFLBOcndR7MsEqUfTQXUKZRk793L6JX0cX3fYhrn043Imb331rcVLJPjLvrxlXBn88Xrq7y9dXotuR/8/g5gavnksb0W6YFj49WZvk7LsWGtlJZ5yiR9mVKKEuHiqA4neVZSSA2V9GhsQ5mx4E9gt/YpFhy+Sj+/NIsfL42XEovhIso09WTgHZt2SdJ87/vVOCJi4UzeMs/qIZpDwPrFmquINZVMvlb4fXZZ0dDOktWPWtzwijGKWHG8iZMmLUaD7z3eWKWJcdFux/2pay4T2mgiuopv7ec0wl/yqbE1kn+mpp77pNww3LLMIJGNO5T0S3x8uIPtVmHhn9llvdvxk3H3r95JYwvVCk6Qvr8MLkjRVD8Flr+FeTNTuHujXUilZDuwjXN+Ocb83HCLW9idVObnw834ZtyMpQVy4TbtUWfOMvfKQGL3orslGX7kzS90QzCLOqnv+M5DyhIMGL5JxGlKZG/60Ddcqj0GYM42Sf40cbJt72VnJePA34/PpBOq0W87y9Yi6a2kvq5yBcxRPLvDHQb8ndciisenhp8twwSkL9LTCHty+5BaOm5vTLrUPNllkz6ljJjeVNghZS7mYjTvA7uTbthb+K7PHbA8mfeDB3RFf/4/Az8+62o+yXjgkryNojUoa6dD9fPcH59aL3rLP9Xpq/E4rVhRFZGCkXbxcHXv+Zfwy+fS7ngWoqy8qTmv7RGFSz/tpKzSfc8joCVt32j8rDtUtSiDZ/UXIrjlvw9bM+t6zi3xBByTSW6iQaWP5NjNFq8ynPMpXh99uqoHCXJPuxVEjoit7zRdVNbWK7UHk0x2NBWwjl3vYdvPzAFqWSfMpxOqoWtlPyj5NAsxfA3DRJsgejCjKSnmo0XwuCE5ZPv1U9/GvxdLvk/9vjDMDauwPesZwAgIMVKJpybWkt+HpU38n+/tRB3vDE/8ntg+XdkVEFdoH428Pf9MWbxPZHDOs3/sgem4PhbwtWt8joaQAoX8fDZ0YQMfFhvZriBdtjP/9L7J5e3OIrh1d96LoxcWdtKDh5873OOHFWyD/MKU5NuyXHRB14nuMf618T2rCJ/qa4S27Es+7i2ur7Z8cXvervkAXj4g8W46D8f4qmpkjODNNpgFr9o+Vc+F1eNBVasLJ8ua9RY/lK9GSIHZbJPpVCRf7vY+C0zlH1s3vJnBgrR9cQExxhTUWhbq72WDhbnuVIsk/w/WLiWL16HdMlGn/wral9/3R1NK6Okz+AEsk8lmTPQYNeu7RqnRY7GefvwcprKKmXkLxxb8l40I478Wcq0cw1xeGvumspOfOc275PbZe6GF2fhmqc/w/gZ/mShcseo+MWItkOR96PZOiSfTP6K82Mhr0R3bSRawncfBwBYut671xWN0uZK2glfLo/Vla/23bnpPeCG7YGpD6VbL6DYo0B0qlBp/rLlHzVANzW6EfmLDSZO9gGgFZkt2Lgv/xccPOP32mvpwHuulL+NoFd+pjo7bOP2CizRhlavQVdknW9Yjv/c/ofgqwUbo0nYGSitsHJBKfDBv9mXyOEcbKAx2bW1XTGvknp/AEW8dQNUqflTSqu205gWVq1frlC+WeuHuWhlfukql8JgXkNtCdsuRR6eMWAbedFi5sifdXwyZy1vbMW9b0uTw/ICJu8P6Tv3W8Pist6jaZ+LHSh7lkLZHvt66vxkfGn5P4G2BuCZ7wJv3xL8bjsuxs9YJY5+GpYAf9gmEjre5t8DPrwDQ8Tyz2Sf6iCV5U8CXdBRrHXThvGPkFH6l94yOfIvd8LXt/Sp/8ienFz+fra242LhmmY0tdp+XpWBcBbwz6z/Ylzh2sDSYpZgh5anUxeY85L/dzSf63IPALfsAbSsi82mPWaxzuJ1Cbo7R/7sZTcIVTql3PjSbOz4qxeCydM0ULlVvjF7tX6BUc7fBpEj/3BNhbxeQTHhq7P8XRcFRv4kp7X8c7DFS/i4/P4p+N1zM7CgnptTkIOW8b+5tpjH9KeBW/cO/Ny9Ipv+pdST/d994EPhe7VDIRiUq6vWsI39dfwcXHr/ZLzLr/pe5xs+0/8n5CEYQcqJeKnMKgN0E6Pbkr9BSGDlOJzmTxV/iefJ10rf8EwjrO5yNX9mWbKrTf18XfA9rdV5w4uz8IWb3giIL6no331oCm54cVbkd74D3IMs8v7YsBJAeF8ds/zDuhnZ9KEwsgCAE4wp3h/t8fGDRGmNAuN/jeFFb5L6iSlL48sgbOXpkxCocpHX2HcWecUpo0OX637WyiZ88z8f4rpnOJ/05jXA/WcCD301tMiVnZK0Ulll+UsaeEvRhuNS2A4NyD8i+3CSRx7MYBBHOSzgnvC8BfKXHCN42afUAjx+off3ck7eq+svlFWuc7YKlqGaUWQBwOTJP98z+NMLqyJB4y7MZFnbpWgr2dE0lEI5N9OJ6Dbkv7FdYfn4L4SrqHidlhmx/Mshf+5Zl6v5M2vb9R9Z+AKQ1J41k+aLcWqSvBpe+HQl/jVxPmatFEmWr4Og4yy1ApP+BtePXdKh91F6fuMK12IYCd1xTUSH0SpLup2L0dILrcCkW/Gr1VdF0inhk2w7DV3wTLhKlY3N5STq3xxka7WxxSPgRWu4Tmf1TGDB68DcV0ILlCN/Nu0Tlimd5U8pxR7XvYxrnv4MJccNNH/bSLb8KVV7colGrULaYb85Gi8c3juHbfaugYVNS/4GuI6yEJL/St9dmMUS8qAOEcKX6YMFnkz1POdqLoeokBtWFt6hUij0zRYp1nePvBk0SFdRDbr3OEr+6UlctPzLe7i6QAQU6Ru/3KDStq+TbhX9pfmJT5uR/wd3AuOvwwmNTwCoguwj4aX8z7nrRy1cVR3wljjrPHvSjahFWySteH0KTPDmNdqQC5R+3YQvk/PK6dDl0roq3lbp9JzLJuu8g9ForOUfkj9re49PXiJq/iTnPUcGh7f8QycBVV0LjzvO8qdU3fD4TiEYyUSTAYBF1BO+MhbVnI8rzf8pj8XB0Fj+JVshaWoiw6qcDZ6ayo8243lk01P/Vkv+0apTWv7+y6WcxtM0KBNyPly6/37dc8fTQPD2UUkE79zunc+9CE1tJfxm3PTAwyWw/Lnrpg3JWy1t1ISL0WQ+LNhwWBPySanW8cLqzl21Aa/NXKXLIh4K8u9JQsI2FJa/oDow/34nJAmTI+7tSEKkzqblwDLP3bAVBe666glfFsG0HPKXn0WgbwuJFORvh94+bsD1cZo/Sxy2W/dTr4M2DYKS44qyz4RwMp8PD54jzPKnye0oTvMHhZLaeMvfH4kHdyPdTo6I58dV+yWWv+H7098D7j4httgMOtmHlUMY4THVgA+quG4BaGt0K0jBS01u45thX5CtlPyTNX+bm4XnZR9GHLoGftGc7+qvNXNcbLH4MMdKzX/ijd4nN7S/Zfwc3PfOooDsWal4AkwbLkK22CodWo4iyzCucC1+ZT0cWv7+m0H8l/j12fX41tjJFeWf9CIE90548g/v5R1f3uI7WF4qcMpo9m00D2vVJ97l5EVey6YAnz4RbgrfAc1fnUiRH2cYsOcXqj6sN+BPiGr+Nc9cCsAzRlwXqIFn4dvyKtP7Tw/+PN74yMuGqi1/O1HzT1gU6XDvpyTDyl2ZKY3A4jqj4Mi0B4GlH3h/L/kA2Lhae47Jj7gUOp94Pf94W2P40237YeDbv+FS+CM0jvxXNbaKpkRG/lVCCvJ3XBpq/jRaDU3WAGXW27ZKPu6qhte6XnkuT/7N61WNLzqEtB2KHGz8Nf8v/4ho+VMQHHr9BOX1ZMhFrVQqHUi8hr63sSAkf7/sNU4znspfF65Env0ioLCC4gsa/yKYwUrg8Ab4F5J1hjwhWZyOq5L5Ajz7Q8+TyAffURj8qm7X9UIgPPktL1QIypR9VPod4rRzVnheu/fLJROUwE1s67KogWAaBC6l6OGPqvZseENb3lq0+1lTVdSHFBO+nOWvemcUlr8OljRyiCd/gsONT8Uf7zkBuGln7TmC7MPdC6tlsaNTe2f1mftUmJ8idMVfXp4ZuygtW+RVIVraowszXpkhShCC5a8g/4eGXo1Wmsccd2i8pq56QR84S5mUyT7HGNNw4vOHAgveEBNoXPZ7gQ9Z4FvYMS6nUz4XO5/VG9pw8t/eChbNhGdQXPfMZ8E2fmlhcBKFbPnvUZqO/Yx5uDr3EAahAXjkXODxb5aVv5OWRLm655+RKmqpRTjZIy5Cz5T/CF9tzhPMYN4+b90E/K5f8Hsg+5Rh+UdlH/8aPJGrZJ91C4BxVwIL3ghGWeFUUpRk4nacs0zvSfZEa+SYjMDa5iZ8c6ZmJKua8A0mfjWyj6OQfeQ6suoAAPdY1wPPXBGeGvN+9iXNeCh/vfa4CqaG/BkE7zqNC61rensJELg43Jzu/x0i8v5mln91cObt0YBOs1aKW7zxiy+oYmi3Fn3xpjsaDsz4zR1UD235R7jw3g8ijZfFCx9jzPZ+WPKhfCYA4AePfBTsU0uIOMx1QdAbG/GHnEdSVDGk/Oq/xBj0T0xZihkrmiKWKaXA23PX4EPJhe2TpQ1YH7NPbkj+FE7QcXr31oN62v9GWgOTWdv1s7V5qVDflDL2DR8AjjfG/E9em+VlH0rTL4zjh+ZBbJ/JYgdRieYf6bqVlr+i3b1xPfDRWOD+M/Dl9V45AldIflI1gJ78meXfiyTXt8k2GUFItjWWiTFkFv6Xvw5OiYuMq4qrkRQOpeFzrsiS7MMW5vXYBgAwgDQBUx8IjvPk7xrlbeDC0IsL5W4hnvx5y5/qvJf8ej/VeI/7hQp/x8k+1YgvlIQOkT8hpD8hZDwhZK7/2U+RZl9CyLuEkOmEkE8IIed05JppsKE1+kDYLj0MtusGFa7SgNtLjm/n0QT/bfVDmjinPnJesj8+Cc7lV03yXiYUBKPI8sg5fMOSL6NzQaTwNuPg069qasPpt0/Cibe+GVPOkPxLkuXfyyf/VhTCek0IskUpxc+f+CT4PnVx/OKt8ESO/JMsf072OWjHgehdk25FJV+vOm8fRv4b29RWoAq6CV8xUXx8mmElr41EXT1Vlr+cF0XOIKCUem6wCeBX+LKyF3Im/pi7F/sb82A1LOKyTpjwTdI0Astf/JmRv4zeG+cH5SvWDky8FwEbV3vOCjp7gJur2NudiZ5oEdrVzx77SH2aX+8FErZ9kfw112FftwDZ5xcAXqOU7gzgNf+7jBYA36CU7gngJAC3EkL6dvC6sVC9pPLQkNf8+R6YpWq3XVB4Dyze8k//lAb38fyXB/ZkHiTSi0BCIuevyRMX9ZcayeWN2+BEt6CMUs99jRHnZ8sacfCfvJWW9Rv0exzwDdeRNH82SmmmNeGIRREjhke77eK/3Irlja0p91/lZR+uHhmR8hZajrPmTKLZelABvu513j5MzvvGvR+kyhOINptwrjb9JCBzXzViXSPVlr8BCtP0Jnx5WVEH9iwpaDCKrsmFNeLwnUuSq2eSVSutdmV14uZ6RtOumYez3/8KfmI95qUtd7OYm3YG7jpGLwSy/FZ+gjtLV+P3uf8IXNLYHN9xttLQW2xrk33OADDW/3ssgDPlBJTSOZTSuf7fywGsBjCog9eNBZFcwWzHFa1Ag3jWMHUjDY21jXbbgetb/kviQgHEPDSVrpszCQYw8n/9j8Bfd+VLHvzPbyXHr2iUyV9l+cvQrymgKDluUM759eqwvzIYkQqav3T9ZtRy5B9v+cudU3sppQUtyD4qyz88zktnJqhiklRdR/x5Bly/gxbPldeQpIEsCbL2mTjhy2F0m7fSOSR/hUcNy2++6BRgwIVlGHApFaxTHUyuPW/VAAAgAElEQVR4oZ9RbMYlYz25sqGlFDz1RWs04R0ErZ/7Hgf/fiIppTob+YvnsXSJN181xpijTBOBqk7XztWfxtI3ewu19iHzcc/b4RyZ7HkUwruHVoQylCG0JWmstwWS/7aU0hUA4H+qx2U+CCEHAcgD0IeFrAJkK7hNkl9ypuFr/g5A1DazZ/kTGKA45y5F1EeGmMYmyy8u9Ugnz68Q3Bj1hTdAg5j7BFHNn8+W/R1H/rqYM67rES/rG9KuAwjJn4bkL1n3LSiEnZbknvjMtGXC5KjcObXb5ZE/pRSvzQy9p1huN70S7j3LW/4GUQy7V34q/wJAXFAkd7sMdTnvNSprG1pZmlNJginDEsfKPqzES8VRCYFnBFGki1ZqwMWL+V9im/sOh0EIatGGdtsJvM/ufpN7pZNcPVPKPsHXYGAjntcXG2C+8Ufx3KRtVfkREOdeql0lP/tFL3ibH+6DguCjxaH3Gj8yFIrh13ubhvyTLP9OUH2SyZ8Q8ioh5DPFvzPKuRAhZAiABwBcRKm6myOEXEYImUwImVxfX69KkgpyY5Zlm7xlhN4+mlCqUxc3AL7lHw/9cVlqcl2P/LWKQyD7uEKZTUnzFyko1N910C0Co6AoOaHsk1bBYuEAdiQrcClbRLNsinhNaobl5jT/8TNW4QePTsPfJ8wNfpO9ZIrFlOQ/8UYsferX+OkTn+CX/1OTNwP/khKiiMt/55HA4vdjz9Nt5jJvlUcGR++SfkDLW/6uS3Hp/ZP9svGyT0ryDzMVP70MlecYcGH6mn+aPZEtOBhprILZvBJn5j/AzJqLcTqdGJxpQif7qCZ8K9P85dHSb3NjMUQK+c3LPvVUseBSIP9Q2qyhrXgqfx12I4vF9PNeBe48CiiKo39WliTLn39XTYn8N7eff+KsF6X0eN0xQsgqQsgQSukKn9yVKycIIb0BPA/gGkqp1oymlN4F4C4AGDNmTMWdX8TyV5B/oPkTQ0t6qTbt0zy0AoqRxupSb3JQznXSvDU4dMcBMIIwAqHlL9+Pp5jymn8a2Uev+dtuKPukJX9mRfeO8RKxYUZfjIbFqP/csxCXN4Qrdhskjb+U1vKf/hSGAXii7SDhZ0pppO55C14Xo0fwOAnS8uSvrsfeaME69BZiCSWB749bS/zogkNKQhi44CnAHopYy18CAYVlErg0fr6IoQ8JFx7e6N4MADi3x2TQdi//vM5DRg7pDArbpfHEo/Pzl4rJ1h6IacJrt9F89PZ58rfD8/ejs7CfOQ+/sh6K5tm6DlPmLcMB3E8NLSX065GPhJoIiuFfl38H9jHC0RHxZxTDcon5DHBTOj10AB2VfcYBuND/+0IAz8gJCCF5AE8BuJ9S+ngHr5cKMhHKXjd50/C9fShATK1blSf7lDmB5GMnslwp+3iSsVjtF9z9Pu6dtDCw0ky4aLfVlr9XrhAb2rx0cS+wztvHpZ7kEsRD1+YgIh8JcaHIG0a0TLfujQs+8KaFeHI+/dYJ6MF5nKSWfWIgPxbegh/VPgNXOA9ARkt7VPvmXUQNoo7qeVfeI0P+mSWWT1PbQqeUUvbZ5Z2rgEfPC276b6/OCTet11r+FKav+aeRfQYhulBvRJ9ccBfCvAFP/o99HZj7avBb0XZw0X+iIywBGvLXCAbiqdyD5yWXAHyd2qEB0uoH8NPNf7z40Tzh+8wVXrBDvezDFmOGZb7QGh/8nRQj7Cbnz8p8q4mOkv8NAE4ghMwFcIL/HYSQMYSQu/00ZwM4CsA3CSHT/H/7dvC6sUhj+XsTvk6koVHJyg7W3OrMYk2DbEdOK/uo3seFa5qDqxmEBlYkIUQg/1bUCJb/jBWe50Ws5a/Z9YuFnWW3ptP8ZW8QFuQrDi5IzJBYvNb/8r/G9JpvBd9XdHCPW94dkYF/Sb+++i+40H06ct4nS6MEZ0mWv+rZscnGSlf48nozn/2vn/kEyYg+s3fmrcFPHvvYz1D9ihNQGMQzAJKlTaAf8dpZcyGc1rNoMVgt3Q/cOhr5nXjjTxz5u8n7HgSxfeRyid/V3Rpn+UOxKTpP/sGGQUCb68m/OsOmThplzF/jda468gc3ilcfjWhawtd2VcdVZXSI/Cmlaymlx1FKd/Y/1/m/T6aUXuL//SClNEcp3Zf7F92br4qQSUeeUPQsf6b5G1qTt3/PGvSr8wao2uBpfqN+Z764u9CJxhT0ekGMA6STfbwyutjg7zZG4ArePlHpIZpDJZY/09oD8tFkcXPun8L3Akm2zH+ZeySW/PusmQrcvAfQ1oS9jEXCsQWrN6hPSgmqoI2cpN2rT4z+Lmv+cShP9uGMDO40XvN3nGTLXzX9TwjFupZikEKFIDxISs2/xu/wF7eGrov9103DPobn+fLP/G1cMaT8lk0B5rziH0oxmpE0f1YnaWJRES5/teXPtd23bw7+LPpCVEFj2NQRj/zZyMBxXHy2rDG2jQP6jjXJ26edKDquKmOrXOErD2NlKzCi+Qtp+ZThhK82Vrvf2O5+a6Hw889zj6Jm5hORchgEIAprbPaqjWgphRJOa1Et+xhwtSKV8I1SzPb9sXX7/QbkH8g+0Zz3JItwgqleyJIEK0Ye+tLqu4CmZcCKjyPHyt0rV04vW/6HjxqAfjVpXHGSyT8ul0o3c3E1ln+aehDIJwj0RtHg7w+gk33YGgFKk6/jwEDBD/7WhB5hHjoiV42Gixv8QynERUWZN7bbmL0ifvMe7wJh/vxeDAE04RjYFXWjWhYGnNW3Q4FT//52orePrnNI8vYpqkYtVcZWSf5yxcoGg+ft40Lt5y96SrC8io6rjt3iDyP5OCc6eJo/AZEahAkHJd+1FPBIhicS/uU04SopSLbenvtkBb5465t4ZfpKrZ8/kynYPasMq/PN16I/poTOihKOWYXIsTQyBA/5hZW9CR+65BBcf7o+kFfcdXOSq6dLESEn2w9xUY7mzxO+QP6EeIVvawqI42NzL7iUYHXffYDR5wr5GArLH+AjjMZb/m4Ky9+GhRrf4m2iIfljr6+oT4jR5pc3tCY/X0VDfHLKUhQj9ct3mn7Hx0/4Jln+HFhd68ifyT7sOmwNie5e2ISvvmOVzvvkUeFre0b+lUFuzPJwMS/7+XPHRZ3eCF4d23GD7fpEeOlrc8l7cLouYBISsRYKKAl6sQEXrSUHd7+1AOtbitGFRgrIjXCRr0lO+Xy9VmNtl2QfVTNO5fGkQdzioTjyTyND8FC9sDJ/kFJyCIMkdycCiltenRP5nYWxKCewm2j5h3/3rrWAF64CbhiOk01vYvS6uquxc/v9eOHAscAxPxfyEWPEh5Z/IFNqLX9G/smdrQ2Ls/zrwgN9h4sJV/hzFDHkn6pjl84n8DbMSdcSk8hf3UGz9yqvkTRriUj+bGGfXvYxhPTRo1RpxDF0ec2/K4IqvBfk6rdMIvj5a2Ufzkq3XYqluolI10VtPgX5U4ptsBb7LxU9TQooouSEFj0BRf2Gdvzh+Zl4ZtryIKgWECVGVaxwINy8prG1pJ2vYC6grgv84/V54v6xrMwdIf9Yy9/XpOUY8iif/H9gPSV8p1BIWGnIP+G6uheZhbgoZztBnezjujQIhDeSrPQvbMKB6T1Hqb5E8uHJvxzLP77TaqNm8CwF3/m+24sJ7zwKePnqaLRaDt4qFX09bey1E7DkPWD8r8P9Coi3MFNeua+KlcN7+ygJVCNVpbf8PTDy17XV2vZ6fMN8Of2Er4RM9qkAthsdxkbc/gzia/6u3qcY3hC8R9sqLKo5Hy1FB2PfjfqBAwDcEvJmclU6lOJkd2LkZSug5Ms+HuTy8+m/br0qHA+Gu8IZFHe+6U3E1eRMrZ8/s1RdSvGXl2cr5aHY2PcJiCV/NipQWInyBt1J+Jb1ovDdpTS6V0EVyV/uWJhXSTn7I+hknzUbi2jxjc9av4Okfht1KQVMkRQMheYPJG8TyuYv0kz42jAD2edZ51B8t3glrup5PTBwFyklBd69HXhJFeJLvK4OKxv95/TObcLvOTN+wSUb4BDZz1+GTvbx21xvTZwjtucBMwZfn+UtaQrq/7ArI+f8LjdWOzJIsvyLmeVfPmwnalnIE76mYXCxffSLvPiHM2XeMv1FnWKqTdQp9SwIFyaecI4Kfq8hRRQdUfMXyhtpQOHxMLyyuHowOE6IXvOXJnyVZe6A5a+zoq62HsQw4ntHKci/XM1fhYhnSCmF+2ii7MPSeR8raH8AwHPOoQDK2yZTN9p8e96aYD+GwL2QsJEFIpa/SP5Mh+YLHS/7zF65EQaJL3eJW5LVjBq84B6Cmfm9gW12jz0v7rrlHreMeKryPIKo8B6Uo/mzerQ0hgdbUFbjRweeu3qjcB56DdHkq76fy63nAg8iFdppuqizHcHWR/6uG7GsVzSKVl9g+cfE9gEgjAq2eff3+os6pVRWn0spCijCNkWdu4CSoBfL5Y9+Vwx3Nccd19XKPmyeIY60qqX5z+PcN4OQEIByGJ7kPpcICtCGJfiy8VZIBiksf7nDWEbF3dyCsMZ+nbT4ERtZHZVD/hGphwNbWBaQsk98jusChV5A/52CtJ/UXMbfAQCJQBuWQAUCio8WN+DeSQsTCblEQ0nT9jsCyzSAuv7A2Q/gsQHfiz1fhhXTpPiy0OA3X/aJSJ4hDERHBmo/fzX5P5z/U0yJw1GYfI1A7TXUZC07dzD0I/FBFDPNvwI4CtmnqVV84KZJQm8fTWwfAILVNGJ9zKrExe/CtJMtS8f1Iig6Rl5oFAWU0NRmQxeqwZIakKB1Eqb5qy1/26XaCd9Q9tGXuVqyz7/fXKhOpLD8Te3CGTVkkqagqH3pR7glfwd+d6CfVwrL35VcYh1ph7dgwV/wXZxvKUf20Wn+AGDKIQN4y9/KA1d+hBeG/0SRqVgu2rgUaFqqvD7/jiR1tjsYYfDBom+R5liAqj1Ox2e9Dos9nwcBRdz0WFAu6oYT2IQkyj4GUSzuVMo+5UdgBcRItgCwDdbjqfx1GGL6e/dqeKTc+SuGNrLpyX/Tjy06GY5LI8NY2fJllv8nS9Zhm9Yi3AGaSRnO8o8bouHR83G1Uac/7rrA1PtB3D2QRwmOURAahayNJ8k+coAo75MrN/e349LERV4L1zRHjt1/8UH44/Mz4a7tCPmHMXu+tkJjWSnWIMidXRLedz35gcDFzbk7cOJrC5BrWwEAGJT3y8At5ddBXlTFb+Ho5e+TKhVJP9zopJwJX17zF4/xC9JcSgIjhEmL81ZvxPLG5P0WnIZl2hdcDC+cvr5ZnZT4zXMS9tzlMYKsxonGFO30Cl+WwS1zAHh7YFimERmDivKmgvxV1vPsF6K/pUCOMPL3yneh9TL2M+ZhNF0AGwYsTR2YZc5fNdAe6EuaswnfSjCgZwHfPmoH4Td5WG0anrfP0nXNaGpztDts8b/2SNjtqMaNsSw/fgR49gc4ft0jnuxj5IWGy+QRtr2gARcEbkCe8svJhycmEgHxvwFex1dyXYwiS/F+4bsYjLXBsbhwBJbpxR2sircPMbD32hfViapg+W+L9SBwMYKsxpfNSejhEz8AWK5P+gkbygBi/H+A36gmKCz3f1jnjHTK8fYRtoGVLX9hs3kSjBL+NdELDHb8zRMxf62qPYqyj2Pr75n3nCnHOmXkf/xuYZgHopE8VLCIi6/QV7THecNNnMeKny/wyF98fkryl+S/V5wDomkUyHNhzIHQvddbd2NoLf9y56+YG3gm+1QIU6pwneVvwoUDQ/vS8pZ/D9IuWLKpQSnQ6kXoq3MaUUDRl31CyPkaoLg1909ML1yMWrQFlv5815tUOtSYEZaRvfBchnyDc1xvw5avmhOxLWnAWWa4v3GcX7pleBPhHZl6DTX/mA7kiYsiP+UU5P+RO0qbxeHmdPzQehLbkbWRY7mA/FPEI5LI35Zejz6kGb2xMSBjRlQEFAeO7Fee7MNr25E4RGE5HC4QQNF2cbsfClspx0kbpgj76kqQ95NNixJM/OtrB+D7x4WL5kicdFomeMOGlWq7te/CLDZG0h5rhlFiCInKQiXVuIczAjbSGnzkJi/+A6IxfPj6d0Eii0UZypV9WNtXrk6uMrZK8pcnWZTePv7cAIWhn6iTPCW2J9GNVxJRagVKbX45DORpCbYhTvjek/8rjjI+FrTkM8x3YBEXQ8i6wKJ5z90DAPCj3JNhEYNPNZE7LoXjUCym2/r3EEbd1rmAAkDB8uqlI5p/ECQrThZoinpRqZbMP2hrI4sDAE4yPsQ2WB/Ny/UJUDPRx0MOPSBb/ldaT2NC4argObFO2fR3xQLSSz98PyM/Bj4shgtDSMs2qFGPyEQJ0CnqpS6elFSyz8wxageHEqxg28rgqmVY/kkoSqTXD004Ycq3ceCzx+MAY67mLE8Sk+VRZR1Jo6G16J2qXGyClq0j4OeDCihp2/hhxvRU+TOwzi9b4VshZCJ0bNHqswyCQfYqHGTMhAlHP1yXHuh+xjx1uhi4N44EXv8DAO/VLKAIh+QjVso11oOcq2dY/rPNiRjtB88qKiwZlvYbh4wIiy1b/i5FM/W0015cDP44S7UmZ8Jx05H/n0vnKH8PRjTSBu6tqok4DqoY6UkaqAtDiDnPUJblLz2TIYqRxEDSFJV9iBcbH0g/6ctb/sz4YBvB8yS2nA7ADEVMG5fqLf80so9I/tFCyxIYgw0zuFeGaln+4/f5m9DGDdcOnmm+Pdqx81CtV1C6KUvy3xrVhi8xCGUf6Z41dXCG+U5Z+bPRZEb+FUKWfYgt6nymQXAfvRp9SAt2NZbCoYhYMzyo71t9njlBPLD9EYllMbjdggxawj72JzBo1AqtQVHw82+k3gTyt61ncYn1IhxKsJ72ipy3vzEPABXumTXQ3jUWbNcLocsIVQ5Ny4Mf6NTkDDiUehOOMbCpgeFnXKM8plvk1ZqgZ6os/1JEfxfhguB3ubHRvByf/N1k8pcnnwcQTXRRJvtwayxY+0nr7qny9mE+5BYcrNr+VIxpuwMnF69XF1VFbBuW++XxvY9iZB95P1kZKoPo2tI3AZCoz32VLH9iWGIcK1pCj5j2KhRB8vEHot5a3o9ifmtpOss/KKOk+TM0tnV8Dwoe2YRvpZBjg0ieHgXLECQC16WB1SWe6FfPqOMAAIOJtLtOoWdZxapzPStmQ2FwdCEa9ygMQoG+IwEAD9jH49T2P+Do4q1Yhyj5A8CXjA8welUYn57dSZ+6nK/5h0Pi2hivpUcuPST4u/+0OzDYXpE44WsRV7s2ShfbJ2kyS03+8QSjK2fv5kXAb/oAC99MHMWk2SwEUE/4mga3CleBFY2tQqRW1YbzNX58KIs4sEkea9BHW1cy+fAIaqKY7N46BGtF2af3UADAxsLg4KfvFH+AHxS/iwecE73ybSLLnxiGMOoxqR1E00yCBTdq+aukGKcEfg6qWuT/4Acxi0ArQGdM+G51rp6AIjyCI1r+vf1OdbE7CHc6p8FxKSyDRGwMNuFLcp4VPphIQ88yfYaZMbW074HAqtfFY9wiFQKK2poCprUdgGvbLubuS01OXzUnYsxybvIrsPxzsF1vpy5GqDUxk9bb9an1zkMzer31e9yGgXgQxyrTUjMP4g+hdbtSHWFEYwUBvq94TJ9SGfmryXDwmnfDNEYehqsnk7R6fdTVkwZRXVVqCaUUh14/AUdxe/wKm7n4DSMgfzhoT+qoYiowsPyLURmMwQDFnmQRni/8SjxwyHeAYQdi6artcGL7jdiLLMSL7sFCEnmUXI63TxwMwxC8fUxaxAVWuqiytaaLP+buEX5zVKNWux3I1QGlZhDQVJr/MjoAQwMJ0K9b6fkUE8JpMPyo+B3cmLsLec32jwxKT6UqY6u0/CMrYm2R1vvmPIt0rHMiHnKOh0Opt2JRRhDkW+PDn3KD7SB58Bn1WeY3ZvfK74JKHgS6xThyGAUCilHb9ETB37eAIiR/tvmGCmw0z+pvABr0vgpm2Dh1nDmIRD00gOQhbSXkryPDfCnUy11FEDkhD+lGfl66VHMtD8KEr6m3/Nkiwzfn1IdlUcT2YcEBTbiRNQYy0szF0JiFbQZc7EBWKA5YwIhDUHIo5tDh+J97VCRJZJRcNcvfFGUft4Qvm5NSnbtdLwunmB8Iv2k1/1xN8LUN0aiyMp71w3cA3GS6VP+yZ5gO/HseB3nie1Ng6yR/aZEXlTw9+uY8cmEP3vUt/wjYsDGnIf8UHiRiQXyrQRFvpTdpETR/L0CVmE7nNiYTAQEwtG8tLD+GEaUioaq8YgAvDhAQupvlYQclOK94tXgrfoiKZloApRT39FATpQqqiWseKvIvJsY6UddN3g51eyeR/MXOVTsxHcTQCTX/nEEwkqyA2xadnF3RFPXJ50vLRoQ1uVDzt2k8ocbLPv56hBjZZzeyWN3B+AZHnCeYTP5E4+Yow1Zp8Hw+hoFldGDwfWOzfuQiwyhE31EHBkrUxGral/uxCFi1QrpvFn8amzdveLDAcXJrS7o3HmnIP5vwrRBEtr4ckaR7Wd53tvzboWrNn7DqsQrqIaTGI0ILv1hUkdcg0oidDM8SM0AB14lY/jqCk1cfe5uOePfkWf7iwiHe44cHqwPe15pZt2wVbXjAa5yPO0cDAMbVfRkfuzsq85WxyHc71UFJ/nKH0W+k8LVXwiI8wJN94iBb/tFFXh6I3+mbnOafh403Cj9BzbjLIunXboxKbVTQ/NmEbyj7lDok+/ifMZb/3/L/VM+T+AZAHPkbkvFiptjICIjvsACP/P9Q+lr43SnhJefAVHkj3yPykwsDu7aPxZfauUlzu12w/L1yxXdeouHhe2ZFFpQljyDY2TUx+1wwZIu8KkRkJy/JQu9tepXPvE70lr//m1mAYSosz7JlH98i8yMQ6kDggtAo+ess/x6RSbHQddB2XVBKhYVDtRoPCkKAa4Z/gvO2CUNX/zjnbUUZIQpf9qEgGNavLvg7DR6pOQ+32Wdqj+cUemhE9tn7qwCA9dSbdG+E+PIX++8WyaOmGLputpCopSj7+eukF9sRN/IYSBrx29mner8tnxJJr/KcEVf4ep9M9umo5R+sRE4IZqd6Xk9OXYGdfvVCYkhoHkaMpxyPJKnKIAY+pqNwXPtfAACfr16fegEasaNt2vP/McTn6JQiln9Sp9SCsLNgdyqTf2tq8k9XV51h+W+VE76RBU++5V+T81atDqr1GhR7YA6lSusl0B+tgjepJbsLlh0kyn8pExqbCS/onOytoCP/GsiWv6cjW4av+VMxWJjO3dMgBJfU36AplVQ//g5cFARf2G0b3Pra3HQN+7Tb8P1+p+LafxdxpfW0MokqvEOE/P26edA5Hicb72MjFV/oZedPQONtR2JfY36YL/f8Wkkd6qhoGcsT1zodt1QqASRsZwcac4INpIzWdcDNewrpD7IdTCqI1v/AxwtwDIKi7WIfAkwqFFGz3ERbwUENKaVybdUhsPwTvH1yij2W1yyZBcc9AC1FvaQpW/7ydx2S5jGIP+nU6HfohsKDh2FtbjAGlLzNbt5w9sFRTrRNuyC44gujcP/r3D7RTjtQ209I5yR0tJPcvXCrfRa2QQPOzU0EgYvLreeENOknaDPy36SIuHz5lv+4K47ALtv2Qvu8NwGED8xx1fHCA/nBzCkntZrb2hAdbOrBRvoUCREKQb3NsSXLX3eOKuAV9dcueL76FJbBkT9pUw48TOkl/m7pR/hn7hZ1IX3LnyehVHGADrgQPZc3xlpbqWQfrmNsRg2OMj8VDruUopXqrbFWow6yjSCv0LI1rwcrn+yxsY72RH+yEZQAZIdjwt8bWzFp7hoh7TGDB2HakgY0tJRw+KgBmDRvLYb3qcOSdS1wQdDe7wTEjQ7jyd+fi7DjLX/VOgw2Sb8yJnCc/KrI7UaHJAubucvyK911Hm5za/fFgNJLADzDgNjREMkODIwYUCfWlWsDUkh1J6Hd2jBwq/0V/Nh6DIRSHG18HHFm4OeH6mkf3GN/Cb/IPSpnlTrYQ2eEd9gqyZ/Ifv6uGI61QL2GzTR/WfbZY0hvzFjRBJMtxjILSvJ/f5WBY8twdGAPPokkb8vdDtDtUlv+8gvCyN80CNZuLIJCJNQ4y59He9wkVkD+UowTCY6Rh+mKVm+vQi52lKAi/4h1xtxwQbGADsFoLBQOU0pjF5O1kdrIb1ReJKQhK53X1Xy6HfqTOZhR2A87nRLuRPXZnHr8bKYoBz1y2CH4yf2TscG28Y/998fPZn2Es0cMw2OrvRDMlxV2AKD3zIqTUFKTv0J77gdvgvzJj9ShoL38pfAOKSktWfM3AThBO4ojf5uE5FiEqYza6oKgZ8ES25rrCruhXXvqHnjueX3YCCCUajyjzVXOEbRwhkY7clivWZNDQfD94hX4e/722Gt+59iobFltbJWaf2QbQF+bJ4zcfC1UkH048mc+21awyXg+sopxKnbDT0rfLq9gwUbp8ZtR15KiJylJJlZq/ROe5Vt0XKxobPM6Au4lOnHnXthpUHTMIu+RGudeyXy7gzOoels6h4tj9Kh9DACgV41VtuUfSU9CgrjNPiuSnlJgA/Rhth0SvTcqTeDraltVvmkDTkYT9ep04nKC3a59Kfh3zdPR9Q78hG/Rn0Ngfv4AhAVhacvAwNwjSUIkU9Waj5vss2PPAbwImjxO22e7xHOA5IlVI7D8Q683ncHDP78SLG8iVwKFgR4FSzRKqCO4KQ/rV5vYKYXt2vs8/8g9I2l4Q8Oh8eHcmlLoBZcfow9kWC1sfeRfbMHw2eIyfyb7BIbtes9KDPZedamwapH5bAuWvyTBPGQfi/Upg0IxNPubqtMUw2SjrSFyzXqo45AoZR8AR+0cLiriyeLMgcvQvy46rDSlsBNJvvWAOH+hIn/b9CbLltKB+IXtecJEXkgJqqie0bzD76qX16XAOikcxsy9f86drXo9xd9U4TQAtbeGTfKYS73VsTPdEetNnOQAACAASURBVDh9n+3ws5N2xbB+tVi9QUVM4S20lfzV1xz5T5i1OnIOj54xK19PML1RBklwRVbJPtPpDrHnAJwR5YNN+AMAjpEWjXFIIllGmTnf6ylus3d+/qBIc1rLv0feFEZJJbuEl2eFE/9500jslNzA8vew87bR97DZ5cgfRhjUUAIFSayHzkLXKEU1UWpBr4aZ4m/+hK9BCNC4DBh/HQAEmrBn+YdVwSSggAytQsTy1waDG3mkvmj+OZTGa/4AQIobIoHlHneOxudkaDStZLETeJbliP7hS2nCRQPtgSZaB0y5D/uWPorkIy+GO3r30KI776ARsAucv7Q/mhI0f4VM5Jg1kd/ylgErZmGQasJX98IQ6djvSxfgrJ4PgoJGyHvt4MO5a0TlBOq6mOzugiZai7Pbr8UMOlJbRhk2yeEG+zzs2nYfnnUPw8l7D8F3jxklPAMe/CKv9lLU8l/WEC/Z9CTxx0eT+ahpWR6b5pe5R2KP6yDbLj3yJn5b+rq3FuSYn6tPQjL5l2yvHnoUPMOEQL+qnXKEXYKlHOU4vuXPGw6NzW2CUWOZyWRMQUBI6KKt2qClzTWF9LoJW4rkeugsdI1SVBOqeB6+Vw4BgMZwT9PWYMJX1PyZ9l2y/Ng9Tcsjmr/Kyn2115nACb/TFy1w9Yyu8FVBdvWkMDDBiG6ZF93blMKlQCEX1kUONmyY+LN7AQBgsB1d3UmkMBiXHRPqjn1qc7C+9x6Q84asxCd//spyd+ZSguaawcHfPExLP6pQuXrqZB9AnA9YT3thPXrBdRGJhcTvz2CoZBPqgoDiY3cnfEDL25y8rq4HABL4Z7O6z6lWjkOcW27z91XgyZ/hoB28TeKH9q3Fx9edGPyetOHNuMK1qG3W6/YdgTw3RAjBs7VnYF6P/YXfS9I8jTLQGp/eN9J61jDy13v7OFwZirCUezRTX/PnDZQc7MB5gMBz9Ehancs21GH55I0o+dsu0HjYL7yywcBTzhF4QBGGnILE1sNMpFsrUw1sfeSvsCgJ5WSf1nB1K3tRS44raP5sEPDR8AuBPsOB3U6J5Kuy3KlhRs0iBQo5Sz+c5RuGyjpWdG7yC+JN+NKATIaRegwn9XBg4El6DACCPm5DJB9TcpczrNB6sQwC9B4C7Ogt6kKw0lHv7bNr+1hsqB3qH5MmrxXeVcG1Usk+Hgio0DHYML21Df5EsHBNgfyjL/Ca9et9ySyd9wqPfI04gVzwI3TmLfV9UsHyZ+QfTXvBwWGo7j6cVPeYc0yqcq0ubJ8qXTlQufVPuOoYvPaTo4XfZOs3ydWzWPRkqLq8d961uYdwmDlDmZaXGz1LPvo+LacDUJMzBccJCw5K3KKtnEmSA/5BnIsoGNFrNbaWcOdEb4KeEgM2LExxd1HkFT/SuNi6HviVIuzGJsDW5+2jIEemfRqESOTvNbKSEwbmCtIBaO4xDPiRP1nHrPBRJ2D6nDl4zd1PdXH1yEPCsbsPxoRP1Mfm0+1QoCWMNFZFNH/A72AiV9VY/j7x3JG7BXsbizDbHebpm3UDcHLzUzgwP144z7pfGqpyE2PBC3/UT4GlH3rhrFeK7pUyaZZgCXrpD48Pd00yLQuaiM8pZR8S/M8fc2DC8UNaTHVHoaXHMNT5FjDhnrGpWKD3LfIsQICJzmh1wWLgSu6DAflrLf+wtm551dugpVZh+es2iWlHHh9tfzH2//ze2HKtLOyIPm3LtRFWK4HKr793TVTmaEdOmJtIIlkm+xTyybTEjyR1EstSOgimQYQRleXv3sdgGiSxU6LUa/8B+SuSb2y3ATP0CgL0nV1cPVAjp48lVmVsfZa/gjBHt7yHV/I/xeAHjgTG/5pPDAAo2a6g+TPk+N+Y5r/tnjileD2aEA3n7HUgyVZjzjJxwPb9lMcoCJZQf6JW0ZGorFIV+VMABX/irB/ZiAnOvji/eLXHwsdejSk1h2IK3UX4h5GHA/me3gT3tnsBvYdibP8rMc45NFzFOXR/4KfzgB4DItdWlY39Jjd409S/cKoJX5n81zQXlcccGCj5C9taUYN3TuH2YODahs5dky9zOaCGTP7etXKa0Acvf7YKTVIMeLbCl4dlhp1n5Bh3D3S7/RUpgBKxUgcdqzbkAH5uQgygvYd6Ml0a7yG+Pia7uwrHLiv+CN8o/hwOTBgEqM2FnYkJJyDlFhRgGsmW/1Un7e6NIPzvBTN+vo7JOiryp5rfGVIulq4Ktj7LX2EZH976OlxC0DbwNOTyJpDviTs2HAH4Rn274wqaP/NmEJatMyK2pAnMq+ai/rnfYtCsB3HQyH6pLH+AYGAPtQ86BcFq+B2DokMiMbLPdHd77Gl87lEXpYGMUEARK2l/rEUf5ECBMRfjX9P2wKT14k5VZ511SiTvN3qdhteXH4IfR6y9aCtVvURskowi3IgcACxVuAx2LIXs859Ji/BTM3rd/UYOwAerKE67/W2vlNxpRpLmH1yrfFCd5a+Rff47eUnktzqFxcs6D1XkVIubeLQdV2n/2jCrPsGYNpxDO80JzSTJq2ZonwIW3XAsZi1Js11qWCFrpYn91bQfplHPVZIQgnYuVIVFXDSjBn8snY9X3QPwd0ISg7J9ef/h+MnLa0D99zFH1C0ktPhF8m+nOWE/67jOJm3dVgMdahWEkP6EkPGEkLn+Z8ScJYRsTwiZQgiZRgiZTggp0zm+3EJpNlImFA2n3AV89T7gjNux32EnBMdkzZ8E53AZMJnAklaN9twGg0buDQCosxBL/oE2z6VhO3YFlwEJF4xw93LNKd4EpEr2YZjhevpu/zoLV31x18D6LKAUrGZmJBK/0oAvj192TXKhihTHQ9lHHJ8oyf+YX3rHUlj+obcVFSypQr4AmwtKFkgU2+0Pws/rxGzcUonl314nBqtjpK+b8JVx3al7YLfBUddS1crz4BhXT81r1ZO7JX+bk2oiLT/VjdhH+K7cXEVI4D3TXrXJoRJY8Mbxzv6Qw6XwoRYMAmyQRlg1KOLfzqlYSIfAIAR/Okssp+JqcCm38lgT0ytskV4F9aj1DMV2zsZO8vZJGyqjGuhoq/gFgNcopTsDeM3/LmMFgMMopfsCOBjALwgh6VaFVIIYcuQr9pAdB+A3p3kborNQCAwsmfAgmM+0bPkDACMypxQ74RsEVCMGWFP5WUmMAknhhSuQ72W3wb25c6X78of/7CV//sojcOTOgwLrswbFQBc9fV+v6uOiNgrlYZ2FfF/cIisGpeUfyD7i+UUrlM2aWFyebbznoSJ/FwYm/l/onsrylSd8TcsS3XAJgGvqgW+NF0ZNcbJPqjAVElp7iROrrOOdXx8NO6DCxUfsIHhnMYSyT7RrtUh83QNAicYvqKsEaQ2HOYf+RfjuJsXC8ueYhvZNo3lTfKnuYXyn9EPh/iY6ozGLhpPk7B2e74aT/8KeAQbBsXtG3acFEANjtu8XGgWa9RPsuAMDPztpV9xynheRlJe/kiZ8O1P26WirOAMAW1E1FkAkVCOltEgpZW4khSpcMx4x5CsfylvcRJDCQhPSO5zPvwwWJ94txVr+J/qLb0AMYIS3ZeJiKbwxBQnWH/CdT+AXrsifkeWAXj6J+mnzlgEDLvLECUJZ/Pn/vMnM9xeui+SjQmj5x8s+FGrSbLe9HFwQQbswLTMYqdxqfwV/7XsNsMsXvftRuHq6IHBzoWUcxn+RtsA0TNgc+RuEeCu0TSsIHAYALZZ+4+6yLf/DrkRrH3FFJrP8py2JelXpUFDMg7CRg0r24TuwgmZnqBLMzWb5mzW9cFnxR8H3Gt0MPwCcdCMw0t8TO410SilQ6I3DdhmC352xd/DzS64YApq12/9y3lEGaBiw10hxPULwn4sOxGVH7eR9lzqxST08F1w+fEvBMgNvOTku1dZi+W9LKV0BAP7nNqpEhJDhhJBPACwBcCOlVLn6hBByGSFkMiFkcn19vSpJhyBXa4HTYy2F7EPSWv5sMth1FFfRFOSwKzH5tFcxk4oWowsSWP45Gk5qBu9+DPkfu8cQIbVlEhT8JfxsKKzcsSwGzMskXZuMJvp0eVNwjOcvfhhdhIX3ag4PZC516AIiZM8TtC1Y/jmJ/LkcuLqbOPjCsu5Dh6V0IHDi7wXX1VNHD0E/3y2Trd5NA9X8AGuX6gnfsJ4MzQb1JaS3/J91DklOBMUoUIO8RfCKeyCafWMmbv9oHPJtrpGly980gPsvPggH7RhuABPZ2IgAuw/pLbSX7foUgk6VEALU9AF2O1V/IWKgV00OA3r6hp8k+7zd05OQeeeGnEkCo5B3LU20/LuS5k8IeZUQ8pni3xlpL0IpXUIpHQ1gFIALCSHK3TwopXdRSsdQSscMGjRIlaRDkBstP8xWb+bCYbdTvAVOgxQBl1igqLQ7exEDIATtfXaIHKIgIVFz5B9v+fubijCZyE+bM4xgCX9HQ8Sm0fxVlr9O9nFdGvzmwAs9ze4tLm6NCrzmW997L0H24SUKwrmuEisa2C0oWxnkH2zowrWtm766T9DWeG+fp757GB6+RNwPl4duchgILf/nrzwivDY34WtSNfkXYSW6MgLA70tfw/dLVyamA9IaAuGohSXvlbYJprhAjraHz5Z7JyjEd9kgBI9/+1B8/7jQzfi4XQcELrgmIZ68eu5DMeUxxE/J8pcnjB1qePfuy8GlsjR/fTGqjURvH0ppdJmaD0LIKkLIEErpCkLIEACxAUkopcsJIdMBHAngibJLWwFcSoJtHeU2VcPJPnxY2sDbhz/hlJu8fwCA58WMGOk6JaTyFYkZZlKEVnqe22w88PNWTGgHfvGsHCs/BYobYQDYy1gEoOMbQkeGo8F3r1y9a3KJ5M9LF7ZLBY3UoWGeOvJXvReRxXLSlpv8y7TNgDA8RUObvqMuR/Zh+jFfPfwocuxFB+H8u98H4M3bFMvYHhEAetcyxgzrOUjPWaCGhvxtakQmRFUo655Tsj8jf9Y+4xb2CUgh+wwpLQZhCqxA/gQ1loFmPzCeQYCeBQuo5eRa6gadsqrO9eVimr/YPncZ3AdYLLZ1yyCBIiCSP1GGQWHYkmSfcQDY+PlCAM/ICQghwwjx4uf63kCHA5jdweumwlf7P4HLS6HmKFcsb/nz7TKUfdT59pfdNA3O8pfEWbXVRbj/o8fYPqZ9WxYFvwbZqvYdYBYgewke+zrwryOAfx2BB/PeFnaNtJydBxSlStD8f37SbhgxIHoNfijMT1pSGlr+Lgyvc/OvsafxeWIZ+AlfHpFNebivQwcNCP7WxmYCMLiPQtrTgJE/TyL836O2DSe2CYF6xzgNrj5592CtCXv+w/qFI5aGnb8c/D2l78nKPBaua4et2oJUgmrF+Wn7bIcv7TU48nvaWwgINpibkE785gve53BpNBRDgOtzg9FK83i511lhbtw74VJDkDYNlZREKSf7pLiRoHPxE0uyz1fGDMdDlxzMaf6GN4rzlYCiRP5xI7G0klo10FHyvwHACYSQuQBO8L+DEDKGEHK3n2Z3AO8TQj4GMBHATZTST5W5VRm2VSfsiBTV/PmHwFv+6vQML/5ACt7GZB+F5W/kFBPEhC0Xj8IFwUfuzkC+Jz4b6fWrx+++bRgLJm6hTLDhfA/gnAeBcx7E5cUf4ZvFn2G8e4D+vBgw0kl64fce1gd7btc38jvlPvvV8S54BKyGXUpiyZiBL4IutUyuQofPxXGPmwTdd7h6AZ4KpoL8+ReYd9U0DSJEjwWAvnU5PHKpWmvv3yMf8Frw+AnBqG28DqV1xDEY2fYwRrY9hIl91dtiOjAxlw4DAHyn+APtfSjJf/QQ3PG1aLtJ6+3DCDY0TqQEQ/cHrl4VdgIpUDTrsFdpLD6uPZh7UUXZp90OyTnkfp783diJ9AgC2Udt+RNiYv8R/QDOmClYps8HUdkn1s+/K8k+caCUrgVwnOL3yQAu8f8eD6D89fJVgGUQYTcm2fLnY6mIh4jitxCRoVlv31Vs8F6R1mRYBUDeSDtW9vGDg/1qGRZNWghgBob2rQk0f2GR14hDgcXvchfzO4b+OwK7nwYAeNntWP/OrHXdcFQgDSnNF/fcFu5sdn0D3zxsZHDMMokg+1Q63JXPkldq63KNnYfl6/jyN4GGJcB/L1AmZeSvWqDllYcbERAS2fXqG4dsj0N3Ckck2/YuYFWTNzH65f2GYmVTNFQx6yjDjoWgTTNN4sDA1aWL8axzKF529Zuhq8hf+8xTNqloncijMtPzxFKAEiOyKZNXJvhbk3JCFVcgF4YwyR50xA4ni7lO4OwRJ8NxmYjXkef2DM/SD/378zh4l0HARu+9n+TuhTHGHK6M+rZelgzVQWxat8vNjEjcDqleRcs/Ct0QLDJ0HzIauHSCOpa5qbL89bIPr73yqvoxuw7CWfsPxZcPCH2Y8X/3APyyfjYqqAn3GahWW4rko7oHiRXO2HdoQCn9etYIw3HT8HZFAjw/6IP96JXxZeCvpr6xvCacgoyiG5OOEIy74nA89/0jgCH7ALvrPUHYwr2eqoAvECd8DYNEvDlqJYJ8/1fHBxvtGAYJ7pmP7cPIvy5v4muHeO2hzVG/yg4M/PT/jsaz7mFqi3PnE4X7EO5Nww5pO+pBvQp4+FJO0onMG8XRj87Y8NqM41LOqufJn6BXQdER80ELqYsjdvak1Z6qtJGLSrKPvF6BGDCNsPsc1r+nF6qj/46Y9dWJuFXYbCg+mPuWJPt0aZgGEcLKyvXKu3r2qeUWYjArW5Ov0h1r6AH+7L70aFXrAhImfFUoWCZuPntf9OvB6dHEADauFr8Dnuuaj45aEtpFXn7t7Mpp2vJ99a4Jt2uU78sySGA1tyOHn560K5Kgei/kV0n2mGm31ZbdNr31+/sCBKOH9cVeQ/VrARjKsfxV6FkTPe/5K4/Ep7/xSVmSfYCQ/E2D4Jdf8lZ+r9eY/g4MYWQhgBiB95qKknREVE6TOmyngdw36cSYBZl68vdgu7zlL8q35x40PHoib/lTB9eeugde+8nR2LZ3dH5n8ckPqK/K6kNe4eu3+92H9Pa/huVx++0IecOjuMn1LWmRV5eGaYi79ERln/DYAG4S16HxUkfspJ0sIpqKYW0M+ScuyOHPJUTsXJivd4G3/NO1poM0lneo+avzOXoXvUturxqLI3/xvkyDgIXK2nFw/8RRWKRcmglfmfxbS2pS/PqhO+ozL8P6YuTfQ0P+ceEZAOCgkdF6r8mZ6OV79YSWf3icSYCGQYKiTlsmSYs+DGIoVw4HmfrtSWn5a6W+KjFUTD3rQkGwctqOG3ZOkuW/olGxy5ktWv4508BOg6LBGQFgxPCRUjmlsigsfwBwAkMpbMtyYL9kzT+z/KuCvEkSJnxFr4Aj/aEgI3fdc4i15kzJmVll+ad8eYjixResJWKI39u9zbcrsfz3GBK/JWXU8FeZ4WJz8iIhskkwMf0Xdt0m8JS55vR9U5WRz2EDPK+XJmmfXrkTaZPJ/5yHgG8+j0IuZrifVtRG6Gtfp5F9kqp/sMLyFM/Xyz4mCWUhXeROyyRidFoZksuueO2EU3Q49Vbga09WcGJyWvar4+o0fxKJ5eOdwO30FRPXSX1RacI3Yvmbfpn879IkPw8K4u3he9KNyktlln+VkDMNQfOPc/U0DYJ/XrA/nvj2oejh64C6Xji2dx64i/dg9/c9YGMsf1UI30R/6zjdu92PI8Np/vIEY6pshfKwUZDmRL5nMsSOr29dLhLsiuGyo3bEqIE+8SlWTSt907ksnnSOwu9KX8c/bNHLRY6f3y7P7O5+qhdGoAK9mYd7jrcoiBl2ugBuSRpuD02nwcCeH0/No4d5XlV8+9FFzLRMM+JhJKCP5wkURJLloLPwE63TMRcBo9S7WKWFrt6Y5l9y1Jo/BcGGNsWaB8Hy14ir+30tyEUqjHidJMvf4C1/sV0EHm+HqONbZpZ/lXDJkTsI5B/V/LmOwSDoVZPDmJH9A7dK3XOIlX0I8R4sWwncWxHDzs943+F9cf1ZewuHaJJPNpEsf74hM8ufk32S2pJuj9nIZbV+/jT6k49te9cEm13LLz4hJAxMpvD4cP2O5H13NxzZfouffZiHCwP3Ol8SY8YPPzgi++w9TKfbJzzDBBh10f0M/nH+/rjx//bWnaJEYrgNRVH+du6+ePp7h6NPbS4oqt7yN+Ijix5wMXDOQ3jCOSpyqGLLX4syyF+xITufg2f5S1o8PMtfuVKat/wPvlx90dNvB65TxLxKMeEL8LKP3vL/w5l7YdwVh0OHjPyrhAO27y/42MqQ3fAYLj9qR/SusXD4qIGq09LF3zjoMuDMO4ADLooeY54yhOC8g0YIh3h5RNkOeItVHs4rZJ8kP+ZvHLp97HGt5q8qHH+xOq/uRgzwOhelJc9eIoVHFCP/Hxa/hyVc8Lsvt/8Wlxd/GEl/45i3gYteFF78n5ywSzAJF0Gc5Z9G9mHxnLh7PmX0EJxz4AjNCSHqFCM+bVEUfWyPgoV9h3vWfyj7aOYcTCPGWKFeG9r9VKiIWT/hWylBpXGqjwfrbEuuGxaZs7S/MmYEbjlHISOyCd8z7wB2iHZ0XuZEPQkd8fOXXT192YeK3wFERl29a/PByE2Fdxes1R6rNra+zVwksAiZQHyj5Y2jvYb2wSe/+WLHLmxawL7nAwveiB6LIZdh/etQuy7FQi5VPsdeDbx9K7BnuPKzzVZPeL7646PR2FrE1MXxUSfDqJ6xyfzEfuqTbgQOYAu/1d4+ADjyjwZ9cY3ob4QAU+nOysxc0wIMUyD/2E46lsBS3GzQ8ZZPaL1qLLQU08Uv0ivy4nFdvBiXUvW8zwEXAaPPib22rvoq1aUNTRz8csDIf0F9c2iwce/BCXsMBvoo4jYxV0+VDJt40UBf8j7k++jlrYKu9dueye1VEZnw70TLPglbPfmvh7i8XodNNtwyFFUcQ/7bD+iBmT86SfhNiOUuT/jyx/oMBy54TDjXNAhKTpQ62CrRJPKHzvIPXgQubzaRVtMHyNVyZdQsbGEvkaKOXBL9Lc0T4jV/ed9bMbOOyT7K55oSPQsWViEmwqUCunsJLX+1wbB9/x5qC/60WxOvqes8K/ZFT7WcNh4beo0C/G24567257jSjNRsX/ZROmCkhE7zz3vrMs7cfxjwKtCvR3iN6HxL1yH/rVr2ARC7wpfHJltZpwzHEEc8eu8k70fJ8udfKMVLkLSTFJMEdBuNs45H7+1DhdTKMurgSgHp+JwU58WRDtN/+TAR/9/emUfNTV0H/Hf9rd7w8hnv/rxggxewsf1hY1bHBsdsJlCzuBD2YzBQaKBhCQmEpOcUOC0ESktDaQJpehIIkJhCgQABkh4CxTQEHJYChYKL2QmU9cP26x/zNJ9GI81II2lGn+b+zpkzI+nN01ukq6v77ruvcsSImJp/lfVoK+FMLLpmzfxQaQ+ZN54fnOA/O7c4gTXgVh5VxZuoEknfEfctuBaWnFnz/9f1ns1T8y8tP+C+VoI8eYqafwjhH/SQ8pp99vs2LOobPxjS0WaL01eeMpObav6NoVKzpyb8fTX/Cufy83xxX4veAd8q+Y4e2uHv+mY5elE3r773KWcum+573BGqod6M+maElZXJ93b6k3+Eh6+AweXLQPilD6Osu/uxZkWzwvrCRSpOUKqMM8g7sK16HiLC31Z4SIjXJOHBWbf4L1bsyDsf9cJ/UuaVVT3veDw+6jB2e+fnvD94Ouy1P/z22pryuXvbYlYNGQt4lgNx3xNBwt95E20P5+Dgj/O2a88x/zgY7DOBzlWeavM83LgXvqkHTSH81/RexIdmEHf6XMvtLQPo3botRbNPxJertiqaWolgbamqZf/whEXc+fTrXHGPfyDVzrYWLrbLWVY6XWDzmGqaf4V2nbpP8OCb5Wen7c6Jt2/mhbc+qvjwPnhuYSGbnScMY/ywTl7/4DPfpQ9DlatCrP++/4cX/vvsuD27TiyfexFxXZ2qvMgkplO6MLxTzTOX2Xj2S58LbfpISh+6q/vrHLHpCC5OIC/fN9kSzT+gzw++uhAKZXKwp01VvIHdvNeQ81CQCpp/hav4lxViL6VBPs0+B10Jh/59cfO32+bwBzPVV5MZaqfX11Xzr6SS+q0U5sYr7KsI/+6uQZy+1F+rD4PTLsHLOPrY/N2pihFMw7Xv2x2TeGnbOB7Z+bswfgETJ00t9lEQIwa1MWd8n3Bd3VOY3l9Z869QnmoPYAhnZ7b86KRFnLOiL3yFM+sz6TguZwz4Vtm+shj6242DQdXjKEF93Q7D0ubnxhnG7DO4C/Y+J6bZxTPgGzTzd0CpC3lpFtlp03wK/91Ohvn+URi9OItl1HXAt2wSiasbXFrZ/rPH2uBdLndM74BvWPt6jTjtEhTYrQSf1cZ8ZylX4NqZP2Zl7+Vs7loMax+Elr6JYkFd5BWi1TxkKmYGITX/2q8Xb+TRpIgSoiFUfhkSVA7ecAlAaV9Enb0bBW9UT+/95pgL3XMKyjNJvFi1kk/hHwHH57quA75eSbj24b7fLs1/7LBOnvnOylJfda9JpeQCTL4OfWafgLzddfnSRTBxUcnsTse/P7T5fUArX9DKNtdobd8p/Mvg9YTpG4uu0eifsOZflr291rb6eGHF4W2Gw6xVJftCr57lQ9Kyv1jb0x+FY3zCP1ThtnV7+DsmlAj/ZNvU9zyO8PeO+wy0b1Sf+EwU8+YBcNIv+e/D7gpOmzJNI/yXTPOPbLibDaw1c+zQdE4cRvPvcAWYqmaPLRnwlfKHQcL0af4hbJejpsMp95WElxhuFzJ3Lz9YCech4+epExyCwluyCoPMrlQl7HVO3+9Qmn/tt46jaGwJsYBNFLYyAI4qjUgpMcqZmuI/ehbMCFwdNpCFk0dU9V5LYiJZ1VkV28pt+wAMtCEyPn0/3Gm6F7NtdPB4W9o0xYAvUIzX4+W8lTtx3JLJTAoZ5iAyHT6RA8teTV133aiXJgAAEDpJREFUWRSbv+evlQTSwskjfJfkq4bzQhQsByrfbGO6Cg/d7qHhbsoBReHv0vz98t2ug92ndbH+yddL3hKgT2hVVAIrCcWUNX/H9zvM6mVR8Mst5PIGvmTT7FOl3e29dfPa3dn0/qfJntyr+XuvgeF2dreNlxSQSclWe0vtXmNxaRrhHxRzZVB7K9MCQrsmwtCxsPYhzr72Fq5ut4PQXqkkEYS/91UzpOZ/27o9qpfV73TBhvbCd7XXbBvbSP7vjVDn883WZ32Fx76xH6+99wnrn3zdR/O3f4vi7VPSB2lr/oX/btmWrH16+cwxZfskhtkni8K/vbVKmazwXzyti8WVUwYzMGBA3GmPp35qtz1t27UDHLceJvSU7L7qqHl9q5t72rTd4+77yAXLailxTTSN2adrSIyZfXEZP98TNTFIXBHC7OMVWq4u3Brsz18r4qOJR2KYXeLy/ZdDJR9nF08f6VpfoW/A1zOwGzCyG/a5FEjamr9j9knY5u+n4HzYFS5cth9JDYMl+QwJq/nHYsRkWPeIzwFXRbqX+IYlYdrSsrf9w+a73wRKG8Ndn5ljhzJ+eAjFIyGaRvg3ml7jeskqk/2uC6Ja2IAys4/dnr5faBe+KDgCoGbhP3o2IKFndp6451SuWTOfwxdMKO67+uj5HNUziZ3HlwZpK8a69+QhAfs9qUo3u11vRmE0/ziTvJwB34TNPsUIoXbuxBm9Z/HZ8B1qzq+eSwqGpS7CH2DMHFh+CWw/q2+fuz3mHllbvp42bWsLDgKXNir860RJ6OHtd/QcLZ8RG4jXe8gR/vuen+qAb7CcqiLAWjvgW+/Air8Mdb6WAcKqeeNLBM/UUYO5fPXcsvDHfgudgMs0Gza2z1m/Kx2ATFnzd8JoDxsUbhA8MkMLprZOeiuHH3fxo5MWle2r58IigRxydclmh5+fv5skvX32PgfOeNS1w9UgYcJEhMBt80/YCliVprH5N5pe29Rbt5tIS4fHsyiKr36Q5p+Sf7PjmVImSOetgVf+HZZeGCKTdC4zRzh5b3e/pQ/Lcd3IHZ6wzynb/Nct3YEZY4ayYna5jT4RbCiDgfJ5aBfmfXyW48yE5j92Lpz7fLEz3eZAX1L183cL/xqigxYyKdlqb22c5p974f/7i1dUHvirkTvO3DPSurOO5i/bfFYZKgnZUEWolJkbHAmYzkXfp0V7DnQMgSNvSuWcYQmaQNbn5p+sn3+vaaFdAqb2R6C1ZQAra/C8Ck1b4c2ioPnHGfD1ZNviHyE2LDX1h0gxZHJhU/iHYxdy2o+fCDhJHSZ5gb+9P1QepY3qngG8dp9pteVZI7kX/mm9WldakMEPx+YvvoOyEvDbL2lAMLeULvog00qjOHHPKTy16QPArfnX4upZaYavj/A//xUWXHo/GztPsf/PsMXUBi8bxOeRQ0u58Xr7/Oa8Zbz5of8KW6nh084VH5wTFtanLAlp/m4Onuuz6l+K5F74Z4Wizd/4CP8omn+g2Scd4Rx7wDdhLjlkTvF3kHmnpklebtp8zD4DR/CRe7H4LAt/u8TkVgbEuiy8z8exwzoZO6z2ENG1FSJCO0/eE8bX7t1UFfd1Uavwz4AlzUGFf51wbP7+7pgJCP8UTFvgGvCt82BUGIK9fQrfiWr+Zf/PsPDf7RRuffR5bnhjXxYGrOQWhqRs/kELwYf7c8h2/ubbsTywQtHhWg+6VrNPUFuMr762Q9Jk+ArOF0Xh713/Ezyaf0Szz5IzCt8pTROP7eefItWee6HHejo9i7z7af6BJ88gLW3c13Usn9POZ1/U/tRuSUj4F2eJ15Jf2HZuba+D8Hc5atRs9vHh1N/AV3+RXH4hUc2/TlQc8I2i+Xsv8J0OgG9/EK9wlU4Xd8JUivRp/l6bfwhvn2ETCy6RB19ZrsWFceOLsZJXPRg+sCCc4jy0kwo++mfLZ/Bx7xb+dFH1xe3LyNJD1hWzqnazj88DcNzc2vKKSYZaNt9spYVbtuwLx/9r+cE4Nv+UcQTs1gxK/6AHUyj9srUDzn228PAsyzhEG2dJKPnwjYNmsW7pDnx5TnivorOXzyjZTiq8w7CBbfzV4XMZ2F7DAzNL7Vyi+Sds9mkAGWrZ/HPellNhyl7lB6KEZa6zxjnbzqqdOKJ+087DUt3mn2Z432zfOsMGtnH+ypkhomD28bX9d+SVyw4qbidl9olFltq5zTXgn6Tm3yDU7JM1otr8U+a4JZNZ0D2CXSYOq564zvhFAAV3YLcUydBNnBaOD/rg9hY+7q194DgWWWpnEegcDp/9MReafyzhLyIjgZuBKcArwJHGGN9g1iKyHfAs8HNjTLhAL81Chs0+IpJJwQ/BZh9n8ZsF3SNIjSwJpZRw2veJb+3fwFJkrJ0HjSwI/zQnk9WJuJLkAuABY8wM4AG7HcR3gYdjni+nRPD2SdujoR8RZJNePK2LRy5YxlfmT/A9XjnTlOLt9EOc0BCdbS10tjXousuS2QcKq9UBDBld2/8zpDTENfscCiy1v28CHgLO9yYSkYXAGOAeoMd7vOmJpPln5+JpNE5T+MV7qSk07nkv68PVRSbi+WdN+O+yuvCpmQy0qSWu8B9jjNkMYIzZLCJlj0MprCP3N8BXgeWVMhORtcBagO7uGtzC+i1RhL8KJwcR4YrVc1k8NaFQ1lVCYn/n0DlwbzKn6g9kXfgPqsV7qNFkoU0tVYW/iNwP+PmLXRTyHKcD/2aMea3aJA9jzPXA9QA9PT3Z8y2MwQE7j+XujQGrWWXY5p91juyZVLdzHbdkSlMJ/7ARQVMl4Hr/9de/xJDO/uivkoE2tVRtPWNM4ErLIvKmiIyzWv844C2fZEuAvUXkdGAI0C4iHxljKo0P5I7rjq0UcCpGYDdFSYksyP6g6727K6U1t9OmP2n+VbgDOB64zH6v9yYwxhzj/BaRE4CeZhP8VYkSz7/aSl+KkhCZiOevyk5qxJUklwG3iMjJwKvAEQAi0gOcZow5JWb+zUEUs09rgjFFlOic+mv45N1Gl6J5yJ3wz8AD1RJL+Btj3sVnENcYswEoE/zGmBuBG+OcM59EcPVUGsu4eY0uQXORl/tBBhTmBmSoPmpDyAJRonoCHHMrbKnzohpKOdP3g50ObHQp8k1eNP8BrbC1l9xo/kpSRDD7AMxo5IxLpcixtzW6BPknQ5pyLBzhn6H65OSx2s+RCN4+itJM5Ebzz97M8Zy0bH8nouavKM1CXu6H4szx7Ch3OWnZfk4UV09FaSbycj84wl/NPkoJUQd8FaVZyI3wz97wavZK1IxE8fNXFA/XrJnP0I6c3cojd4D3XiJLZpJYOMI/Q8pdzq6YHKDCX4nIqnnjG12E5DnhLnjt0fxManSEf4aWQ1VJkxmsRqDCX1Fgu3Ew57BGlyI5HOG/rUErovmgkiYrTKgU+E1RlH5NUfhvaWw5XKjwzwoZsgUqipIwKvwVRVGaEMfVU4W/EkiGBoQURUkI1fyVYNTsoyi5ZfnF0DkcRs9qdEmKqKtnZlCNX1Fyy7R94YL/aXQpSlDNX1GUEnYaM1T9D5oA1fwzg95tSja4++y9G10EpQ6o8FcUpYQBmVi5XUkbNftkDrX9K4qSPir8FUVRmhAV/oqiKE2ICv/MofZWRVHSR4V/5lCbv6Io6aPCPyuoY7WiKHVEhX9WaO0sfGs8f0VR6oD6+WeFw6+Hx2+ACT2NLomiKE2ACv+sMHQsLPtmo0uhKEqToDYGRVGUJkSFv6IoShOiwl9RFKUJiSX8RWSkiNwnIi/Y7xEB6baKyJP2c0eccyqKoijxiav5XwA8YIyZATxgt/341Bizq/2sinlORVEUJSZxhf+hwE32903AV2LmpyiKotSBuMJ/jDFmM4D9Hh2QrlNENojIoyKiDwhFUZQGU9XPX0TuB8b6HLoownm6jTGvi8g04Fci8rQx5iWfc60F1gJ0d3dHyF5RFEWJghhTeyAxEXkeWGqM2Swi44CHjDE7VfnPjcCdxphbq6R7G4iz4vEo4J0Y/++PaJ3zT7PVF7TOUZlsjNm+WqK4M3zvAI4HLrPf670JrAfQJ8aYz0VkFLAncEW1jMMUvhIissEY01SxErTO+afZ6gta57SIa/O/DNhfRF4A9rfbiEiPiNxg08wCNojI74EHgcuMMc/EPK+iKIoSg1iavzHmXWC5z/4NwCn29yPALnHOoyiKoiRLnmf4Xt/oAjQArXP+abb6gtY5FWIN+CqKoij9kzxr/oqiKEoAuRP+IrJSRJ4XkRdFJCjcRL9DRCaJyIMi8qyI/EFEzrb7feMrSYFrbDs8JSILGluD2hGRFhH5nYjcabenishjts43i0i73d9ht1+0x6c0sty1IiLDReRWEXnO9veSvPeziHzNXtcbReQnItKZt34WkR+IyFsistG1L3K/isjxNv0LInJ8reXJlfAXkRbg74ADgNnAGhGZ3dhSJcYW4FxjzCxgd+AMW7eg+EoHADPsZy1wXf2LnBhnA8+6ti8HrrJ1fh842e4/GXjfGDMduMqm649cDdxjjJkJzKNQ99z2s4hMAM4CeowxOwMtwNHkr59vBFZ69kXqVxEZCVwCLAYWAZcEBdSsijEmNx9gCXCva/tC4MJGlyuluq6n4F77PDDO7hsHPG9/fx9Y40pfTNefPsBEe1MsA+4EhMLkl1ZvnwP3Akvs71abThpdh4j13Q542VvuPPczMAF4DRhp++1O4Mt57GdgCrCx1n4F1gDfd+0vSRflkyvNn76LyGGT3Zcr7GvufOAxguMr5aUtvgecB2yz213AH40xW+y2u17FOtvjH9j0/YlpwNvAD62p6wYRGUyO+9kY87/AXwOvApsp9NsT5LufHaL2a2L9nTfhLz77cuXOJCJDgNuAPzfGfFgpqc++ftUWInIw8JYx5gn3bp+kJsSx/kIrsAC4zhgzH/iY4FDpkIM6W7PFocBUYDwwmILZw0ue+rkaQXVMrO55E/6bgEmu7YnA6w0qS+KISBsFwf8vxpjb7e43bVwl7Pdbdn8e2mJPYJWIvAL8lILp53vAcBFxJii661Wssz0+DHivngVOgE3AJmPMY3b7VgoPgzz3837Ay8aYt40xXwC3A3uQ7352iNqvifV33oT/48AM6yXQTmHQKBcrh4mIAP8EPGuMudJ1yImvBKXxle4AjrNeA7sDHzivl/0FY8yFxpiJxpgpFPryV8aYYyiECVltk3nr7LTFapu+X2mExpg3gNdExAmQuBx4hhz3MwVzz+4iMshe506dc9vPLqL2673AChEZYd+YVth90Wn0AEgKAyoHAv8FvARc1OjyJFivvSi83j0FPGk/B1KwdT4AvGC/R9r0QsHz6SXgaQqeFA2vR4z6L6UQDRYKdvH/AF4EfgZ02P2ddvtFe3xao8tdY113BTbYvv4FMCLv/QxcCjwHbAT+GejIWz8DP6EwpvEFBQ3+5Fr6FTjJ1v1F4MRay6MzfBVFUZqQvJl9FEVRlBCo8FcURWlCVPgriqI0ISr8FUVRmhAV/oqiKE2ICn9FUZQmRIW/oihKE6LCX1EUpQn5fwVbBS9I9y3UAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fad8d129048>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "n_eff = [\n",
    "    466.099,136.953,1170.390,541.256,\n",
    "    518.051,589.244,764.813,688.294,\n",
    "    323.777,502.892,353.823,588.142,\n",
    "    654.336,480.914,176.978,182.649,\n",
    "    642.389,470.949,561.947,581.187,\n",
    "    446.389,397.641,338.511,678.772,\n",
    "    1442.250,837.956,869.865,951.124,\n",
    "    619.336,875.805,233.260,786.568,\n",
    "    910.144,231.582,907.666,747.347,\n",
    "    720.660,195.195,944.547,767.271,\n",
    "    723.665,1077.030,470.903,954.924,\n",
    "    497.338,583.539,697.204,98.421\n",
    "]\n",
    "\n",
    "ess = []\n",
    "for i in range(len(n_eff)):\n",
    "    ess.append(pystan.chains.ess(lst, i))\n",
    "\n",
    "pd.DataFrame(data=dict(Target=n_eff, PyStan=ess), columns=['Target', 'PyStan']).head()\n",
    "lst.keys()\n",
    "slst = lst['samples'][1]['chains']\n",
    "slst.keys()\n",
    "slst['d'].shape\n",
    "sim = lst\n",
    "m = sim['chains']\n",
    "\n",
    "ns_save = sim['n_save']\n",
    "ns_warmup2 = sim['warmup2']\n",
    "ns_kept = [s - w for s, w in zip(sim['n_save'], sim['warmup2'])]\n",
    "\n",
    "n_samples = min(ns_kept)\n",
    "plt.plot(lst['samples'][0]['chains']['d'])\n",
    "plt.plot(lst['samples'][1]['chains']['d'])\n",
    "trace_values = [lst['samples'][0]['chains']['d'], lst['samples'][1]['chains']['d']]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scalar input generated by model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     },
     "base_uri": "https://localhost:8080/",
     "height": 69,
     "output_extras": [
      {
       "item_id": 1
      }
     ]
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 1222,
     "status": "ok",
     "timestamp": 1517626800094,
     "user": {
      "displayName": "Sharan Yalburgi",
      "photoUrl": "//lh4.googleusercontent.com/-ypYfwsB3KsY/AAAAAAAAAAI/AAAAAAAAD0s/JStFE5QSLnE/s50-c-k-no/photo.jpg",
      "userId": "108333366686678045862"
     },
     "user_tz": -330
    },
    "id": "qrbFuE8i55bu",
    "outputId": "c979dc9c-2b23-4d01-d9c2-10b00ed2e808"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:pymc3:Auto-assigning NUTS sampler...\n",
      "INFO:pymc3:Initializing NUTS using jitter+adapt_diag...\n",
      "INFO:pymc3:Multiprocess sampling (2 chains in 4 jobs)\n",
      "INFO:pymc3:NUTS: [mu]\n",
      "100%|██████████| 2000/2000 [00:00<00:00, 2737.87it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Vhatold  [ 0.95052935]\n",
      "rho  [ 0.41958492  0.20705641  0.07990572  0.05096835  0.00974834 -0.00199782]\n",
      "vhatnew  [ 0.95052935]\n",
      "acorrmean  0.420976441285\n",
      "acorrmean  0.20913573551\n",
      "acorrmean  0.0805219196838\n",
      "acorrmean  0.0503138208127\n",
      "acorrmean  0.00946224350183\n",
      "rho  [ 0.44361314  0.2207471   0.08543953  0.05365924  0.01068153]\n",
      "Vhatold  [ 0.95052935]\n",
      "rho  [ 0.41958492  0.20705641  0.07990572  0.05096835  0.00974834 -0.00199782]\n"
     ]
    }
   ],
   "source": [
    "import pymc3 as pm\n",
    "# D = [1]#,2,4]#,8,16,32,64,128,256,512,1024]# number of dimensions, ie shape of the gaussian\n",
    "# D = [1024]\n",
    "eff_n_new = []\n",
    "eff_n_old = []\n",
    "# for id in D:\n",
    "with pm.Model():\n",
    "    mu = pm.Normal('mu', 0, 1, shape=1)\n",
    "    tr = pm.sample(1000, tune=1000, chains=2)\n",
    "\n",
    "## convert to scalar\n",
    "tr = [tr['mu'][:1000].reshape(1000,), tr['mu'][1000:].reshape(1000,)]\n",
    "\n",
    "eff_n_new.append(effective_n(tr))\n",
    "\n",
    "eff_n_old.append(pm.diagnostics.effective_n(tr))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     },
     "base_uri": "https://localhost:8080/",
     "height": 34,
     "output_extras": [
      {
       "item_id": 1
      }
     ]
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 892,
     "status": "ok",
     "timestamp": 1517626736093,
     "user": {
      "displayName": "Sharan Yalburgi",
      "photoUrl": "//lh4.googleusercontent.com/-ypYfwsB3KsY/AAAAAAAAAAI/AAAAAAAAD0s/JStFE5QSLnE/s50-c-k-no/photo.jpg",
      "userId": "108333366686678045862"
     },
     "user_tz": -330
    },
    "id": "JcEITzzo5-n5",
    "outputId": "3ca58dab-69bc-4acf-dbb4-12bc3c4d2ab0"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "eff_n_new  [835.20763909944083]\n",
      "eff_n_old  [840.0]\n"
     ]
    }
   ],
   "source": [
    "# eff_n_old[0]['mu'].mean()\n",
    "print(\"eff_n_new \", eff_n_new)\n",
    "print(\"eff_n_old \", eff_n_old)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "vhatnew  [ 0.00331131]\n",
      "acorrmean  0.482815610574\n",
      "acorrmean  0.332760953691\n",
      "acorrmean  0.205459325913\n",
      "acorrmean  0.157693900831\n",
      "acorrmean  0.0919225955553\n",
      "rho  [ 145.80893142  100.49311379   62.0486065    47.62366737   27.76103506]\n",
      "NEW EFFN  2.60257154686\n",
      "Vhatold  [ 0.00331131]\n",
      "rho  [ 0.48345803  0.33414833  0.20791086  0.16368281  0.10153236  0.10553312\n",
      "  0.07091255  0.07571874  0.07063633  0.04890764  0.05859839  0.00773918\n",
      "  0.01506161  0.00163325  0.03169346  0.01217826]\n",
      "OLD EFFN  436.0\n"
     ]
    }
   ],
   "source": [
    "print(\"NEW EFFN \",effective_n(trace_values)) ## facing problems here\n",
    "\n",
    "print(\"OLD EFFN \",pm.diagnostics.effective_n(trace_values))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "colab": {
     "autoexec": {
      "startup": false,
      "wait_interval": 0
     },
     "base_uri": "https://localhost:8080/",
     "height": 34,
     "output_extras": [
      {
       "item_id": 1
      }
     ]
    },
    "colab_type": "code",
    "executionInfo": {
     "elapsed": 992,
     "status": "ok",
     "timestamp": 1517626586502,
     "user": {
      "displayName": "Sharan Yalburgi",
      "photoUrl": "//lh4.googleusercontent.com/-ypYfwsB3KsY/AAAAAAAAAAI/AAAAAAAAD0s/JStFE5QSLnE/s50-c-k-no/photo.jpg",
      "userId": "108333366686678045862"
     },
     "user_tz": -330
    },
    "id": "h7N7FUqw6G9U",
    "outputId": "15e02d32-7895-474d-f912-6f600a115be7"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1000, 1)"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tr['mu'][:1000].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [],
   "default_view": {},
   "name": "pymc3.ipynb",
   "provenance": [],
   "version": "0.3.2",
   "views": {}
  },
  "kernelspec": {
   "display_name": "Python (pymc4)",
   "language": "python",
   "name": "pymc4"
  },
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}