{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pymc3 as pm\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "from IPython.core.pylabtools import figsize\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats\n",
    "import pandas as pd\n",
    "import random\n",
    "figsize(12.5, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Landscape Error\n",
    "\n",
    "We model our estimated landscape area using 3 error components:\n",
    "\n",
    "$$\\hat{LA} = LA * \\epsilon_{classification} * \\epsilon_{boundary} * \\epsilon_{landuse}$$\n",
    "\n",
    "Where\n",
    "\n",
    "* $LA$ is the true amount of residential landscape area within the service area of a given utility.\n",
    "* $\\epsilon_{classification}$ is a percent error caused by improper classification in the land cover algorithm.\n",
    "* $\\epsilon_{boundary}$ is a percent error caused by clipping to an improper boundary polygon.\n",
    "* $\\epsilon_{landuse}$ is a percent error caused by clipping to improper residential areas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classification Error\n",
    "\n",
    "Accuracy across the state has been assessed by taking a random sample of pixels and comparing these against reference labels assigned by a researcher. This pixel-based accuracy assessment is then used to generate error bounds on the total amount of estimated landscape area according to the methodology outlined by Olofsson et al. (2013).\n",
    "\n",
    "This results in a standard error of 3.87755%. That is, $se_{classification} = 0.0387755$\n",
    "\n",
    "However, we expect the error to increase when classifying at the level of individual water retailers, and we model this through an increased variance of the classification error at the retailer level. Because we do not have a detailed accuracy assessment at this level, we must make an assumption about our prior beliefs.\n",
    "\n",
    "We do this by assuming that the retailer-level classification error is log-normally distributed with $mu=0$ and a standard error twice the statewide levels $2*se_{classification} = 0.0775510204082$.\n",
    "\n",
    "Our prior error term for the classification is thus\n",
    "\n",
    "$$\\epsilon_{classification} \\sim Lognormal(\\mu=0, \\sigma = 0.0775510204082)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boundary Error\n",
    "\n",
    "The Utility boundary polygons utilized are not guaranteed to be correct, but manual inspection has shown that most of them are, and those with errors tend to be minor. To capture this we model our boundary error according to an assumed prior distribution as before. We err on the side of caution and choose a spread parameter of 0.1 (10%). This is larger than what is typically observed in practice.\n",
    "\n",
    "$$\\epsilon_{boundary} \\sim Lognormal(\\mu=0, \\sigma = 0.1)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Land Use Error\n",
    "\n",
    "The land use dataset used to isolate residential areas is similarly imperfect, but correct in the vast majority of cases. We model this error in the same manner as the boundary error because both result in multiplicative error caused by cropping of land cover classifications. Again, a spread parameter of 10% aligns with prior beliefs based on manual inspection.\n",
    "\n",
    "$$\\epsilon_{landuse} \\sim Lognormal(\\mu=0, \\sigma = 0.1)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ET Error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several possible sources of error when estimating evapotranspiration (ET) using the algorithm described in the methodology document. Three major sources include\n",
    "1. Measurement/estimation error at CIMIS stations\n",
    "2. Algrithm-induced error caused by our estimation methodology (inverse distance-weighted average of 10 nearest stations to utility service area centroid)\n",
    "3. Approximation error caused by attributing a single point estimate of ET (at service area centroid) to all the landscape area in a service area that may span a large and diverse geographic region.\n",
    "\n",
    "The first source is out of our control, and no error bounds are provided in the CIMIS source data. The third is addressable in future iterations through the use of ET data with a higher spatial resolution and the creation of parcel-level outdoor water budgets. For now we do not model these, and instead focus on the second error source: error induced by our estimation algorithm.\n",
    "\n",
    "We can again model the errors as relative\n",
    "\n",
    "$$\\hat{ET} = ET * \\epsilon_{et}$$\n",
    "\n",
    "but this time we can measure the error distribution directly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Error Calculation\n",
    "\n",
    "To briefly recap the current methodology, ET at a given point $p$ is estimated as the inverse distance-weighted average of the ET readings at the 10 nearest CIMIS stations to $p$, with more weight given to nearer stations and less weight given to more distant stations. \n",
    "\n",
    "When estimating utility targets, $p$ is chosen as the centroid of the utility service area. To estimate ET errors, we instead choose $p$ as the location of a CIMIS station, and estimate ET using the 10 nearest stations (excluding the one located at point $p$). This gives an estimated ET from our algorithm, as well as a true ET as measured by the station located at $p$. Relative errors are then calculated as $ET_{est} / ET_{true}$. \n",
    "\n",
    "The distribution of the errors for monthly ET is shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv(\"data/et_estimates_with_truth.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean monthly error:  1.015\n",
      "Std. dev. of monthly error:  0.13\n",
      "5/95 percentiles:  [ 0.84088676  1.23041756]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x1162f4050>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABAcAAAF0CAYAAABfSlCKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3X+cXmV95//Xh1LA1KJ2U0G0selXi+m2uzRRNFXBH92o\ncb2Lq92ENkWTWtsaWJetAe1iE8TtmtgKEtLWbaPUxg4orpHytSH4A2vQyjrj6lqTdLXiVJHUqQqU\nMaDw2T/OGbnnnh/JPTOZa3LO6/l43I8w51z3meua+z0Xc3/uc64TmYkkSZIkSWqvE0p3QJIkSZIk\nlWVxQJIkSZKklrM4IEmSJElSy1kckCRJkiSp5SwOSJIkSZLUchYHJEmSJElqOYsDkiRJkiS1nMUB\nSZIkSZJazuKAJEmSJEktZ3FAkhaIiLgjIh6KiAuO0O7Wut3v9Ww/t97+0WPbU5USEcsiYndEHIqI\n70+Wgyme96667ZEe76zb33qU7bsfDx77n4DaJCK2HG3GJUmzd2LpDkiSfiDrx2zaHe0xphQRrwDe\nBVybmRtmcyzNnYhYBHwIWAJ8BtgDPAj876M8RAJfBvZN02Zs34eAf5hk/yvr49wM3DXJ8aWjEhHn\nAh8Dbs3M503RbNbzmSTp6FkckKTm+DSwDBgt3REdE08Dngjsy8xzZvD8qJ97xIJPZm6b9AARr6z/\n8y2Z+Tcz6IPUj+3AADBSuiOS1AYWBySpITLzMPD3c3ComINjaO49sf73S0V7Ic2NI84zmfkt4Fvz\n0BdJEq45IEmNMd2aAxGxPCKuj4h/jIj7I+LuiPhyRNwQEZ2udncA76Q6lfeVPdeUf7TnmI+IiNdH\nxGBE3BMR90XEFyLiioh49DT9fFZE7ImIb0fEvRFxe0T8Wr1v0mvXu7dHxPqI+GREfKfevqTeviQi\nLo2Ij0TEVyPicP09PhERr46ICW9GIuKJ9TH+ISr/KSI+V4/lzoj444h4TN32pIh4Y0Tsj4jRiPh6\nRFxVn+7ft4h4QUTcVK8fcH99vOsiYkVPu3Mj4iHg2npT9+ty3F3nX6+bcHlE7IuIr9VjH4mIWyLi\nl4/w3DMi4q0R8fk6c/8SEQejWlNh5STtHxER/7nOwLfqTNwRETdGxPk9bcfW/Fgyxfe+NiZZE6R7\ne0T86/r37M6o1oT4vbrNiRGxLiJ21fm5u87QgYh4e0Q87gjj/g91Vr5R/7y+UY/pkog4uacfl05z\nnP9Yt/nb6b5fV/sfzCn1z/JNEfHF+vfjH7raPS0itkXEp7v6eFf9c37+JMf9GPBRqnnmOTF+nuk+\n7rRrDhzt75Ak6eh45oAkNVz9x/mHqOb8zwGfBH4IeDywmqpQfGPd/L3AM4BnUX1C3X19+oGuYz6G\n6o/7fwvcDXwE+B5wLvBfgV+JiOdm5nBPX9YCu6g+Nfw/wBfqfrwzIn7mKMZyNfAa4DbgJmApD1+T\n/GvAm4CvAAfrvj8OWAk8E/h3wHRvPncBvwR8vB77LwC/CTwtIs6hus7+Z4Fbqc7QeDbwn4AnAf/+\nSH3vGccVVD+nh6hej2GqS0J+GXhZRLw6M99VN7+LqjDwJCa+Lsfj9dj/BdhAlafPA9+hWkfhOcDz\nI+Lpmfm63ifVOX4f8Cjgn4APAw8APwmcT/Wz+FRX+ydQvWbLgPuoMvPPVHl7FtVrOdD1LY50fftU\n+8e2PxN4B3AnVYYeAdxbtzkNeHc91v1Uv4c/ApwFXASsjYiVmTlunYeIOBG4Hngp1foSt1P9ri0G\nfgb478B1VPl5O3AB8FsRsS0zJ+vrxrqv26cZ52ROocr9MuBvqNa5+LGu/b9P9fr9HdV6GPcB/x/w\nYuDfR8RrM7P7e/418F3ghVT53tO175td/z3la3IUv0O/kZnX9jlOSWq3zPThw4cPHwvgQfWm9kHg\ngiO0+1jd7vd6tp9L9YfyR3u2f7Ruv3aSY/0ocHbPtlfUx3nnNH24rm5zG/Doru2LqN60PwR8ouc5\njwPuqfuysWffs6jeSD0EPDjJ93uoft63gadN0acVwM9Msv104LP181/Ws++J9bEfonrD/4SufY+h\nKjI8SPVm7lM9Y30i1ZvNB4GVfbzOL6y/333A83r2ra/3HQaW9fu6TPM93zXT507xOpwzi2M8G/jJ\nSbY/meoN3oPAU3v2PYHqjfWDwJuBE3v2LwZ+oevrAP5X3f5DwI/1tD8JeGHPtrHfvyXT/Awn/H52\n/WwfBN48xXMfSVVA6u33D9XjeQj4q0me94f1vi8BPzvJ/ucCP9r19SfqfvzSJG3/dX2su4AfPsrX\namxOeaj+HfrxKdq9ADhtku1Pr1+3w8Djpjj2R6f5/pvrNr1z3Yx+h3z48OHDx/QPLyuQpIXn2pjm\nlnFUf1T347H1v3/duyMz783M2/s5WET8BPByqj/AX52Z3+k63ijwaqo/zH8hIp7R9dRfp3qT9MnM\n3NHTj33AHx/Ft39rZv6vyXZk5mBmfnGS7XcBl1C9YZzqzIEELsrMr3U979t1n4LqjdWGnrF+leps\nA4AJp05P43X199uRmeMu1cjqbIGbgB8GXtvHMY9W76UivY/OkQ8xO5n5icy8Y5Lt/xe4gurn/fKe\n3b8DnArcmJmXZeb3e547kpmf7Nr0Eqpi0TeAl2d17Xp3+wcys/vT6rnw98AbJ9uRmf+SmTdN0u8H\nM/MyqrMNXhgRPzK2LyJ+nIc/6X95Zn5hkuN+LDPv7dr0dqqf38ZJunFRfaw/zczv9Tc0kqqg981J\nd2benJmHJtn+aWAHVZ5/qc/vOZ2Sv0OS1FheViBJC88+pl907kU8/Ib/aNxOdQryX0bE7wN/m5mz\nuVb9HKpLEQYz8+96d2bmnRFxM9Ch+mRz7Prmc6n+oP/LKY77Hqo/+qfz/ul2RsRJwCqqlf0fC5xM\n9WbpR+smZ07x1O8Dt0yy/f/W/w5n5v4p9gdwxvTd/kH/fojqcgWAP5+i2U6qT5mfezTH7FPvpSK9\nhqfZN2fqN8EvAn6e6lP/k+pdY9fe975OL6B+Y3uU3+KFdfu/rAtWx1oCuzNz2ss8IuLfUBWSllJd\nVjD2Ic2J9X8/ieosFahe/5OAz2Tm0d6u8gPAP1JdnvHTmfn39fc9FfhVqrMK/uRoB9Xln3qKLxNE\nxI9RXUbws1Rn3fxwvevJ9b9T/e71ZQH8DklSY1kckKSF588y891T7awX8+qnOPAG4Oeo3jC9CPhu\nRAxRXUP8nsw8MM1zJ/P4+t+vTNPmy1Rvmh/fte0J9b93TPGcqbYfVZv6LIXrgZ9g6mvHT51i+zcy\n86FJtv9L/e9Ub5rHPrU9Zap+9fhXddtk6p/fl+t/Hz/F/tk4qlsZHksR8RKqU/F/jKmv4e99ncbu\n1HC0We23/Vy4Y6odUS1auQs4j8nHHEwcd99jyMwHI+KPqNYiuJBqTQyAV1IVI96fmV8/2uN1uWO6\nnRHxG8Db6u/R7+9ev0r/DklSY3lZgSQ1XGYeysynUX2K9maqT/J/Hvhd4O8iYtN8d6nP7Q83yLx/\nsu0R8QhgN1UBYidwNtWbzxMz84eAp4w1neLQkxUG+tmvoxARZ1CtV/EY4C1UC1o+KjN/qH6dXkj1\nGi3E22ke6W+m706z7y1UhYEvUp1e/3jg5K5xjy2kOBfj/tO6Lxd0XabwGqrfr2tmeMwpxxYRy6nO\nRjgJ2ER1ltIju8b2Wyzc11SS1MXigCS1RGb+TWb+XmY+n+qN829TvWH4/YhY2sehxj55/Klp2vxU\nfezuTynH/vsnp3jOVNuPxjlUZ1MMZeZv1OsP3N11mveTp3nufPpnYKzAMdXPb2z7TD7hXeheQrWK\n///MzN/NzC9k5r907Z/qdRo7c+MpU+yfbXuo7nwAD1+C0uuJU2w/Gr9M9fuwpl574K6e9QcmG/dM\nxkC9vsJ7qMZxQUT8IvDTwN9l5sf77/oRja3jcXVm/mFmHszM7mLCT8/x92v775AkHTMWBySpheoF\n2f4H1a3kTgD+TdfusTdJU1169jdUn6SfFRE/17szIk6n+gQYqjsrdD8vqG47N5lfPbreT2rstmpT\nnf7/a7M49pyp13oYu+b/lVM0+/X6349Osf94dqTXaaoM7KHKzm8c5fcZa39+fVbJ0Rh7I7msd0dE\nnAYsP8rjTGbKcUfEC6jWXej1UarfxRURcVaf3+9qHl6Y8ELqxfv6PMbRmm5sJwP/gcnPCjrSPDMp\nf4ck6dixOCBJDRcRv1PfYaB3+1N4+BPLr3btGlux/2cmO15m/iPV/eZPAN5RL0Q2dsxFVKc1nwLc\nlpl/2/XUncAo8KyIeE1PX55JdSbDTI0tFvj8iBj35i4iXg38R47isoV58odUb9x+OyKe170jIl5J\n9en6A1Rv8ObSQhj/2Ov08rqIBEBEnFDft34lk/fzbVTrO3Qi4oqIGPeGMiJ+vM7QmBupbr13BnBD\nd0br9idHxAsZ78NUr8ulEfGo7mMDf0F1Pf1MjY37op5+nEl1R4wJY67vDDB2t4wbIuJf97aJiOdG\nxIQzHeo7G3yU6ne4Q3UL0b+YRf+ns7/u4ysi4pFdfTuZqv8/OcXzxuaZJ9eLDE5lsjyU+h2SpEZz\nQUJJar7LgLdGxAGqP+S/S/Wm6VlU91n/857V0P+W6tZqPx8Rg8D/Ab4HHMzMP6jbbKRaffzpwJfr\nRRK/T3VHgsVUC4Kt6+5EZn49In4TuBa4pn7T/nddfXkb1TXL/d5mjcz83xHxQao3Qp+NiFuBbwFn\nUZ3W/N/qn0NxmbmnfiN8GXBLRNxG9anrU6g+nf4+8JtT3B1hNgJ4dkS8a5o2w5m5eY6/b7e/Agap\nxvn3EfFxqnvVP53qTgVvAV7f+6TM/MeIeDlVUep3gVdFxKeosvJEqjU03gPcVrfPiHgp1RkELwSG\nI2If1Snpj6da6+DbjD8tfQfVmQnLgYP18X+E6s4XXwU+SLVuwExcXvf9zRGxhir3jwWeTXVGzZ1U\nhZFel1C9ue4An4uIT1MtwreY6vaaZ1Dd+eDeSZ57NfA8qjfX1x7Duza8i2rhw58HvhIRn6C6K8Kz\nqYqEbwf+c++T6tf0M1S3nPxC/d+HgZHMfENX0wlrFRT8HZKkRvPMAUlaWGb76W5OcozXAO+keiN1\nDtVpvj8J3Ayc17t6fX0P9FVUn74+nupU7w3A6q4236K6ndgbgH8A/h3Vbcy+SfVG/KmZOeE048x8\nD9Ublr1Ub+o6wCOBVwHb62Yj04xtOi+nKi4cAJ5Z9+mr9Vh2MvnPhmm2z9X+iU+o3oC/CPgQ1Rua\nX6Z6c3w98AuZOdUt2vr+Xj3P/Snggmke//4ojzOzDlSnhJ8L/D7VJ8fPq78eBJ5BlclJx5iZt1Dd\nJu/tVG/sX0D1xv9RVLe0+5Oe9sPAU4FLgS/Ux38psITqTh2X9LS/myrT766//wupfl5/Um+/e6Zj\nz8wP1OP8MHA61SfbPw78HlUOvjfFmL+XmS8FfoXqVptPpsr5z1EV4F4H3DXFt/0I1Zv0BP5oJv0e\n68Zkfevq491UP+c/onpdXkj1s95D9Ub9c9Mc4z9Q3dr0R6nO7tlQ/3vE7z+L3yFJ0hTiCLfklSRp\nXkTEBVRnFXywfkMkaYYi4lXA/wD2ZObqI7WXJKnvMwci4pERcVVE3BERoxGxLyKe2tPmTRFxZ73/\nloh4Us/+kyNiR0SMRMS9EXFDRPRzz25J0nEoIn6iXtytd/szgT+g+oRwutPeJR1BfQvDN1D9Pv1h\n4e5Iko4TM7msYCfwfKrTTH+W6jS3D0fE4wAi4lKqlXFfTXWf6fuAmyPipK5jXEV1+unLqE5xPQN4\n/wzHIEk6fjwP+HpEDEbEB+ri8GeBTwD/CnhnZt5YtovS8SkiXhcR11Kdyv+TwF9n5keKdkqSdNzo\n67KCiDiFatGbl2Tmnq7tnwE+lJm/FxF3Am/NzCvrfacCh4BXZOZ766+/Caytr8EbW613P/CMzLx9\njsYmSVpg6vn+d6gWKzuNasG371CtLP/OzHxvwe5Jx7V6YdBzqNbt+Cvgd+o1ASRJOqJ+71ZwItXK\n1vf3bP8u1a2pllIttPODKnVm3lOvrrsSeC/VojUn9rQ5GBHDdRuLA5LUUJl5kOrMMklzLDOfW7oP\nkqTjV1/Fgcz8l/rWPm+sb4l1iGoF3ZXA/6UqDGS9vduheh9UnxQ9kJn3TNNmnIj4V1SrEt9BdZsb\nSZIkSZLa6hTqu09l5j/PxQH7PXMAqvtWvxP4OtV9ZIeobkOzYi46NIUXUN2/WJIkSZIkVX6V6v34\nrPVdHMjMrwDPjYhHAKdm5qGIuI7qPtd3AUF1dkD32QOnUV1PSt3mpIg4tefsgdOY+l69dwDs2rWL\nZcuW9dtlac5cfPHFXHnllaW7oRYzg1oIzKFKM4MqzQyqtP3797Nu3Tqo3yvPhZmcOQBAZn4X+G5E\nPIbqk/3XZeZXIuIuqrsZfB5+sCDh04Ed9VMHqc44eD7QvSDhEuBTU3y7wwDLli1j+fLlM+2yNGuP\netSjzKCKMoNaCMyhSjODKs0MagGZs8vu+y4ORMQqqrMDDgJPBrYBXwSurZtcBVwWEV+iqmJcAXwN\n+CD8YIHCncDbIuLbVHc/uBq4zTsVaKG7666pTm6R5ocZ1EJgDlWaGVRpZlBNNJMzBx4F/Hfg8cC3\ngBuAyzLzQYDM3BYRi4B3AI+munf1izLzga5jXAw8WD/3ZGAPsHGmg5Dmy9e//vXSXVDLmUEtBOZQ\npZlBlWYG1UQzWXPgfcD7jtBmC7Blmv33AxfVD+m4sWLFsVx3UzoyM6iFwByqNDOo0sygmuiE0h2Q\njifnn39+6S6o5cygFgJzqNLMoEozg2qiyMzSfTiiiFgODA4ODrrwhyRJkiSp1YaGhsbOYFmRmUNz\ncUzPHJAkSZIkqeUsDkh9WL9+fekuqOXMoBYCc6jSzKBKM4NqIosDUh9WrVpVugtqOTOohcAcqjQz\nqNLMoJrINQckSZIkSTqOuOaAJEmSJEmacxYHJEmSJElqOYsDUh/27dtXugtqOTOohcAcqjQzqNLM\noJrI4oDUh23btpXuglrODGohMIcqzQyqNDOoJnJBQqkPo6OjLFq0qHQ31GJmUAuBOVRpZlClmUGV\n5oKEUmH+T0ClmUEtBOZQpZlBlWYG1UQWByRJkiRJajmLA5IkSZIktZzFAakPmzZtKt0FtZwZ1EJg\nDlWaGVRpZlBNZHFA6sOSJUtKd0EtZwa1EJhDlWYGVZoZVBN5twJJkiRJko4j3q1AkiRJkiTNOYsD\nkiRJkiS1nMUBqQ8HDhwo3QW1nBnUQmAOVZoZVGlmUE1kcUDqwyWXXFK6C2o5M6iFwByqNDOo0syg\nmsjigNSHa665pnQX1HJmUAuBOVRpZlClmUE1kcUBqQ/etkalmUEtBOZQpZlBlWYG1UR9FQci4oSI\nuCIi/iEiRiPiSxFx2STt3hQRd9ZtbomIJ/XsPzkidkTESETcGxE3RMRjZzsYSZIkSZLUvxP7bP96\n4DeBC4AvAk8Fro2I72TmNQARcSlwYd3mDuDNwM0RsSwzH6iPcxXwIuBlwD3ADuD9wLNnNRpJOg4M\nDw8zMjJyxHaLFy/2kwlJkiTNi34vK1gJfDAz92TmcGb+T2AvcHZXm9cCV2TmTZn5BaoiwRnAeQAR\ncSqwAbg4Mz+emZ8F1gPPjIju40gLztatW0t3Qce54eFhzjxzGStWrDji48wzlzE8PDzu+WZQC4E5\nVGlmUKWZQTVRv8WBTwLPj4gnA0TEvwWeCXyo/nopcDrwkbEnZOY9wKepCgtQnW1wYk+bg8BwVxtp\nQRodHS3dBR3nRkZGOHx4FNgFDE7z2MXhw6MTzjAwg1oIzKFKM4MqzQyqifq9rOAtwKnAgYh4kKq4\n8F8z87p6/+lAAod6nneo3gdwGvBAXTSYqo20IF1++eWlu6DGWAYs7/tZZlALgTlUaWZQpZlBNVG/\nxYE1wK8Aa6nWHDgLeHtE3JmZfzHXnZMkSZIkScdev5cVbAPekpnvy8y/y8z3AFcCb6j33wUE1dkB\n3U6r9421Oalee2CqNpNavXo1nU5n3GPlypXs3r17XLu9e/fS6XQmPH/jxo3s3Llz3LahoSE6nc6E\nU3c3b9484Vqi4eFhOp0OBw4cGLd9+/btbNq0ady20dFROp0O+/btG7d9YGCA9evXT+jbmjVrHIfj\ncBwtGMfHPvaxnq17gYnjqE7UGm8hjaMpr4fjcByOw3E4DsfhOBzHQh/HwMDAuPfAS5cuZe3atROO\nMVuRmUffOGIE+N3M/B9d294AvCIzn1J/fSfw1sy8sv76VKpLBi7IzPfVX38TWJuZH6jbnAnsB56R\nmbdP8n2XA4ODg4MsX97/abjSXBkZGWHx4sWlu6Hj2NDQECtWrKBaV2C6+WwIWEHvvGcGtRCYQ5Vm\nBlWaGVRpD/9NyYrMHJqLY/Z75sBfAZdFxOqIeGJEvBS4GPifXW2uqtu8JCJ+Dng38DXgg/CDBQp3\nAm+LiOdExArgncBtkxUGpIVkw4YNpbugljODWgjMoUozgyrNDKqJ+l1z4ELgCmAH8FjgTuCP620A\nZOa2iFgEvAN4NPAJ4EWZ+UDXcS4GHgRuAE4G9gAbZzgGad5s2bKldBfUcmZQC4E5VGlmUKWZQTVR\nX8WBzLwP+C/1Y7p2W4At0+y/H7iofkjHDS9rUWlmUAuBOVRpZlClmUE1Ub+XFUiSJEmSpIaxOCBJ\nkiRJUstZHJD60HvrE2m+mUEtBOZQpZlBlWYG1UQWB6Q+DA3NyV1CpBkzg1oIzKFKM4MqzQyqiSwO\nSH3YsWNH6S6o5cygFgJzqNLMoEozg2oiiwOSJEmSJLWcxQFJkiRJklrO4oAkSZIkSS13YukOSMeT\nTqfDjTfeWLobarHJMjg8PMzIyMhRPX/x4sUsWbLkWHRNLeJcqNLMoEozg2oiiwNSHy688MLSXVDL\n9WZweHiYM89cxuHDo0f1/FNOWcTBg/stEGhWnAtVmhlUaWZQTWRxQOrDqlWrSndBLdebwZGRkbow\nsAtYdoRn7+fw4XWMjIxYHNCsOBeqNDOo0sygmsjigCQ1wjJgeelOSJIk6TjlgoSSJEmSJLWcxQGp\nD7t37y7dBbWcGdRCYA5VmhlUaWZQTWRxQOrDwMBA6S6o5cygFgJzqNLMoEozg2oi1xyQ+nD99deX\n7oJaZv/+/eO+vvTSSxkaGppyvzQfnAtVmhlUaWZQTWRxQJIWpG8AJ7Bu3brSHZEkSVILWByQpAXp\nO8BDHPkWhR8C3jgvPZIkSVJzWRyQpAXtSLco9LICSZIkzZ4LEkp9WL9+fekuqPXMoMpzLlRpZlCl\nmUE1kcUBqQ+rVq0q3QW1nhlUec6FKs0MqjQzqCayOCD14fzzzy/dBbWeGVR5zoUqzQyqNDOoJnLN\nAUmaI8PDw4yMjEzbxlsPSpIkaSGyOCBJc2B4eJgzz1zG4cOjpbsiSZIk9a2vywoi4isR8dAkj+1d\nbd4UEXdGxGhE3BIRT+o5xskRsSMiRiLi3oi4ISIeO1cDko6lffv2le6CFqiRkZG6MLALGJzmccUs\nv5MZVHnOhSrNDKo0M6gm6nfNgacCp3c9/h2QwHsBIuJS4ELg1cDZwH3AzRFxUtcxrgJeDLwMOAc4\nA3j/zIcgzZ9t27aV7oIWvLFbD071WDrL45tBledcqNLMoEozg2qivi4ryMx/7v46Il4CfDkzP1Fv\nei1wRWbeVO+/ADgEnAe8NyJOBTYAazPz43Wb9cD+iDg7M2+f1WikY+y6664r3QW1nhlUec6FKs0M\nqjQzqCaa8d0KIuKHgV8FdtZfL6U6m+AjY20y8x7g08DKetNTqQoS3W0OAsNdbaQFa9GiRaW7oNYz\ngyrPuVClmUGVZgbVRLO5leFLgUcBf15/fTrVJQaHetodqvcBnAY8UBcNpmojSZIkSZLm0WyKAxuA\nv87Mu+aqM0eyevVqOp3OuMfKlSvZvXv3uHZ79+6l0+lMeP7GjRvZuXPnuG1DQ0N0Op0Jtx/bvHkz\nW7duHbdteHiYTqfDgQMHxm3fvn07mzZtGrdtdHSUTqczYbGSgYEB1q9fP6Fva9ascRyOw3E0ZBww\nAEwcB6wBPtOzbS8wcRxw7STbhuq2k90usbf9cN32wISWV1111bivm/56OA7H4Tgch+NwHI7DcRzP\n4xgYGBj3Hnjp0qWsXbt2wjFmKzKz/ydFLAH+ATiva32BpcCXgbMy8/NdbW8FPpuZF0fEc4EPA4/p\nPnsgIu4ArszMt0/x/ZYDg4ODgyxfvrzv/kpzZdOmTbz1rW8t3Q0tQENDQ6xYsYLqjgTTzVPvAdbN\not0m4K1H0W7SXgIrcC7VbDkXqjQzqNLMoEp7+G9PVmTm0Fwcc6ZnDmyguhTgQ2MbMvMrwF3A88e2\n1QsQPh34ZL1pEPh+T5szgSXAp2bYF2neLFmypHQX1HpmUOU5F6o0M6jSzKCaqK+7FQBERACvBK7N\nzId6dl8FXBYRXwLuoLqh99eAD0K1QGFE7ATeFhHfBu4FrgZu804FOh5cdNFFpbug1jODKs+5UKWZ\nQZVmBtVEfRcHgF8EfgJ4V++OzNwWEYuAdwCPBj4BvCgzH+hqdjHwIHADcDKwB9g4g35IkiRJkqQ5\n0HdxIDNvAX5omv1bgC3T7L+f6qMvy22SJEmSJC0As7lbgdQ6vSuTSvPPDKo850KVZgZVmhlUE1kc\nkPpwySWXlO6CWs8MqjznQpVmBlWaGVQTWRyQ+nDNNdeU7oJazwyqPOdClWYGVZoZVBNZHJD64G1r\nVJ4ZVHnOhSrNDKo0M6gmsjggSZIkSVLLWRyQJEmSJKnlLA5Ifdi6dWvpLqj1zKDKcy5UaWZQpZlB\nNZHFAakPo6Ojpbug1jODKs+5UKWZQZVmBtVEFgekPlx++eWlu6DWM4Mqz7lQpZlBlWYG1UQWByRJ\nkiRJajmLA5IkSZIktZzFAakPIyMjpbug1jODKs+5UKWZQZVmBtVEFgekPmzYsKF0F9R6ZlDlOReq\nNDOo0sygmsjigNSHLVu2lO6CWm9L6Q5IzoUqzgyqNDOoJrI4IPVh+fLlpbug1jODKs+5UKWZQZVm\nBtVEFgfQa9IkAAAgAElEQVQkSZIkSWo5iwOSJEmSJLWcxQGpDzt37izdBbWeGVR5zoUqzQyqNDOo\nJrI4IPVhaGiodBfUemZQ5TkXqjQzqNLMoJrI4oDUhx07dpTuglrPDKo850KVZgZVmhlUE1kckCRJ\nkiSp5SwOSJIkSZLUchYHJEmSJElqub6LAxFxRkT8RUSMRMRoRHwuIpb3tHlTRNxZ778lIp7Us//k\niNhRH+PeiLghIh4728FIx1qn0yndBbWeGVR5zoUqzQyqNDOoJuqrOBARjwZuA+4HXgAsA34H+HZX\nm0uBC4FXA2cD9wE3R8RJXYe6Cngx8DLgHOAM4P0zHoU0Ty688MLSXVDrmUGV51yo0sygSjODaqIT\n+2z/emA4M1/Vte2rPW1eC1yRmTcBRMQFwCHgPOC9EXEqsAFYm5kfr9usB/ZHxNmZefsMxiHNi1Wr\nVpXuglrPDKo850KVZgZVmhlUE/VbHHgJsCci3gucC3wd+KPM/DOAiFgKnA58ZOwJmXlPRHwaWAm8\nF3hq/X272xyMiOG6jcUBSTqG9u/ff8Q2ixcvZsmSJfPQG0mSJC0E/RYHfgr4beAPgf9GddnA1RFx\nf2b+BVVhIKnOFOh2qN4HcBrwQGbeM00bSdKc+wZwAuvWrTtiy1NOWcTBg/stEEiSJLVEvwsSngAM\nZuYbM/NzmfmnwJ8CvzX3XZMWnt27d5fuglpvNhn8DvAQsAsYnOaxi8OHRxkZGZllX9VUzoUqzQyq\nNDOoJuq3OPANoPd81P3A2EdLdwFBdXZAt9PqfWNtTqrXHpiqzaRWr15Np9MZ91i5cuWEX869e/dO\nuoLoxo0b2blz57htQ0NDdDqdCX8Eb968ma1bt47bNjw8TKfT4cCBA+O2b9++nU2bNo3bNjo6SqfT\nYd++feO2DwwMsH79+gl9W7NmjeM4DsaxcePGRoyjKa/HQh0HDAATxwFrgM/0bNvL5HcguHaSbUPA\nRmCyN+297Yfr4x6Y2JSbgeVdj6cAW4DR+utlAOzZs6cRr0dTcrWQxvGOd7yjEeNoyuvRxnFcffXV\njRhHU16PNo7jD/7gDxoxjqa8Hk0fx8DAwLj3wEuXLmXt2rUTjjFbkZlH3zjiPcATMvPcrm1XAk/L\nzGfVX98JvDUzr6y/PpXqkoELMvN99dffpFqQ8AN1mzOpigzPmGxBwvpWiYODg4MsX768d7ckFTc0\nNMSKFSuoPnmfbp56D7CuQLt+2g4BK3DOlSRJWpge/tuTFZk5NBfH7HfNgSuB2yLiDVSLCz4deBXw\nG11trgIui4gvAXcAVwBfAz4IP1igcCfwtoj4NnAvcDVwm3cqkCRJkiRp/vVVHMjMz0TES4G3AG8E\nvgK8NjOv62qzLSIWAe8AHg18AnhRZj7QdaiLgQeBG4CTgT1U58pKkiRJkqR51u+ZA2Tmh4APHaHN\nFqoLWKfafz9wUf2QJEmSJEkF9bsgodRqky0eIs0vM6jynAtVmhlUaWZQTWRxQOrDqlWrSndBrWcG\nVZ5zoUozgyrNDKqJLA5IfTj//PNLd0GtZwZVnnOhSjODKs0MqoksDkiSJEmS1HIWByRJkiRJajmL\nA1If9u3bV7oLaj0zqPKcC1WaGVRpZlBNZHFA6sO2bdtKd0GtZwZVnnOhSjODKs0MqoksDkh9uO66\n60p3Qa1nBlWec6FKM4MqzQyqiSwOSH1YtGhR6S6o9cygynMuVGlmUKWZQTWRxQFJkiRJklrO4oAk\nSZIkSS1ncUDqw6ZNm0p3Qa1nBlWec6FKM4MqzQyqiSwOSH1YsmRJ6S6o9cygynMuVGlmUKWZQTWR\nxQGpDxdddFHpLqj1zKDKcy5UaWZQpZlBNZHFAUmSJEmSWs7igCRJkiRJLWdxQOrDgQMHSndBrWcG\nVZ5zoUozgyrNDKqJLA5IfbjkkktKd0GtZwZVnnOhSjODKs0MqoksDkh9uOaaa0p3Qa1nBlWec6FK\nM4MqzQyqiSwOSH3wtjUqzwyqPOdClWYGVZoZVBNZHJAkSZIkqeUsDkiSJEmS1HIWB6Q+bN26tXQX\n1HpmUOU5F6o0M6jSzKCayOKA1IfR0dHSXVDrmUGV51yo0sygSjODaqK+igMRsTkiHup5fLGnzZsi\n4s6IGI2IWyLiST37T46IHRExEhH3RsQNEfHYuRiMdKxdfvnlpbug1jODKs+5UKWZQZVmBtVEMzlz\n4AvAacDp9eNZYzsi4lLgQuDVwNnAfcDNEXFS1/OvAl4MvAw4BzgDeP9MOi9JkiRJkmbvxBk85/uZ\n+c0p9r0WuCIzbwKIiAuAQ8B5wHsj4lRgA7A2Mz9et1kP7I+IszPz9hn0R5IkSZIkzcJMzhx4ckR8\nPSK+HBG7IuInACJiKdWZBB8Za5iZ9wCfBlbWm55KVZDobnMQGO5qIy1YIyMjpbug1jODKs+5UKWZ\nQZVmBtVE/Z458LfAK4GDwOOALcDfRMTPUhUGkupMgW6H6n1QXY7wQF00mKqNtGBt2LCBG2+8sXQ3\nNM+Gh4eP+EfA/v3756k3GwAzqLKcC1WaGVRpZlBN1NeZA5l5c2a+PzO/kJm3AKuBxwD/8Zj0rsfq\n1avpdDrjHitXrmT37t3j2u3du5dOpzPh+Rs3bmTnzp3jtg0NDdHpdCb84b958+YJtygZHh6m0+lw\n4MCBcdu3b9/Opk2bxm0bHR2l0+mwb9++cdsHBgZYv379hL6tWbPGcRwH47j77rsbMY6mvB7zMY7h\n4WHOPHMZK1asmPaxbt267pEAE8cBa4DP9GzbC0wcB1w7ybYh4G4mP3ugt/1wfdwDE5tyVc/Xo3Xb\n8a/Hnj17FtzrMeZ4z9XxPo5XvOIVjRhHU16PNo7jl37plxoxjqa8Hm0cx7nnntuIcTTl9Wj6OAYG\nBsa9B166dClr166dcIzZisyc3QEibgduAf4M+DJwVmZ+vmv/rcBnM/PiiHgu8GHgMd1nD0TEHcCV\nmfn2Kb7HcmBwcHCQ5cuXz6q/ktSPoaEhVqxYAewClk3T8kPAG4FBYLp56j3AugLt+mk7BKzAOVeS\nJGlhevhvVFZk5tBcHHMmCxL+QEQ8EngS8OeZ+ZWIuAt4PvD5ev+pwNOBHfVTBoHv120+ULc5E1gC\nfGo2fZGkY2sZ07+hnq/LCiRJkqS511dxICLeCvwV8FXg8VQ33P4ecF3d5Crgsoj4EnAHcAXwNeCD\nUC1QGBE7gbdFxLeBe4Grgdu8U4EkSZIkSWX0e7eCJwB/SXUR63XAN4FnZOY/A2TmNmA78A6quxQ8\nAnhRZj7QdYyLgZuAG4BbgTuBl818CNL86b1GSZp/ZlDlOReqNDOo0sygmqjfBQnPz8wnZOYjMnNJ\nZv5KZn6lp82WzDwjMxdl5gsy80s9++/PzIsyc3Fm/mhm/nJm/tNcDEY61oaG5uRyHmkWzKDKcy5U\naWZQpZlBNVG/Zw5IrbZjx44jN5KOKTOo8pwLVZoZVGlmUE1kcUCSJEmSpJazOCBJkiRJUstZHJAk\nSZIkqeUsDkh96HQ6pbug1jODKs+5UKWZQZVmBtVEFgekPlx44YWlu6DWM4Mqz7lQpZlBlWYG1UQW\nB6Q+rFq1qnQX1HpmUOU5F6o0M6jSzKCayOKAJEmSJEktZ3FAkiRJkqSWO7F0B6Tjye7duznvvPNK\nd0OtthuYnwzu37//iG0WL17MkiVL5qE3WkicC1WaGVRpZlBNZHFA6sPAwID/I1BhAxz74sA3gBNY\nt27dEVuecsoiDh7cb4GgZZwLVZoZVGlmUE1kcUDqw/XXX1+6C2q9+cjgd4CHgF3Asmna7efw4XWM\njIxYHGgZ50KVZgZVmhlUE1kckCRNYRmwvHQnJEmSNA9ckFCSJEmSpJazOCBJkiRJUstZHJD6sH79\n+tJdUOuZQZXnXKjSzKBKM4NqIosDUh9WrVpVugtqPTOo8pwLVZoZVGlmUE1kcUDqw/nnn1+6C2o9\nM6jynAtVmhlUaWZQTWRxQJIkSZKklrM4IEmSJElSy1kckPqwb9++0l1Q65lBledcqNLMoEozg2oi\niwNSH7Zt21a6C2o9M6jynAtVmhlUaWZQTWRxQOrDddddV7oLaj0zqPKcC1WaGVRpZlBNNKviQES8\nPiIeioi39Wx/U0TcGRGjEXFLRDypZ//JEbEjIkYi4t6IuCEiHjubvkjzYdGiRaW7oNYzgyrPuVCl\nmUGVZgbVRDMuDkTE04BXA5/r2X4pcGG972zgPuDmiDipq9lVwIuBlwHnAGcA759pXyRJkiRJ0szN\nqDgQEY8EdgGvAr7Ts/u1wBWZeVNmfgG4gOrN/3n1c08FNgAXZ+bHM/OzwHrgmRFx9syGIUmSJEmS\nZmqmZw7sAP4qMz/avTEilgKnAx8Z25aZ9wCfBlbWm54KnNjT5iAw3NVGWpA2bdpUugtqPTOo8pwL\nVZoZVGlmUE10Yr9PiIi1wFlUb/J7nQ4kcKhn+6F6H8BpwAN10WCqNtKCtGTJktJdUOuZQZXnXKjS\nzKBKM4Nqor6KAxHxBKr1An4xM793bLokLVwXXXRR6S6o9cygynMuVGlmUKWZQTVRv5cVrAB+HBiK\niO9FxPeAc4HXRsQDVJ/+B9XZAd1OA+6q//su4KR67YGp2kxq9erVdDqdcY+VK1eye/fuce327t1L\np9OZ8PyNGzeyc+fOcduGhobodDqMjIyM275582a2bt06btvw8DCdTocDBw6M2759+/YJpxaNjo7S\n6XTYt2/fuO0DAwOsX79+Qt/WrFnjOByH41jA46huIdh7CuEo0AH29WwfoFpKZcJIgM/0bNtbH6PX\ntZNsG6rbjkyyr7f9cN22dxxQ1Xi7TTWOPUw9DnPlOByH43AcjsNxOA7HMR/jGBgYGPceeOnSpaxd\nu3bCMWYrMvPoG0f8CPDEns3XAvuBt2Tm/oi4E3hrZl5ZP+dUqqLBBZn5vvrrbwJrM/MDdZsz62M8\nIzNvn+T7LgcGBwcHWb58eb9jlKQZGxoaYsWKFcAgMN388x5g3QJudyyOOQSswLlZkiRpfj38Nyor\nMnNoLo7Z15kDmXlfZn6x+0F1q8J/zsz9dbOrgMsi4iUR8XPAu4GvAR+sj3EPsBN4W0Q8JyJWAO8E\nbpusMCAtJBM/TZbmmxlUec6FKs0MqjQzqCaa6d0Kuo079SAztwHbgXdQ3aXgEcCLMvOBrmYXAzcB\nNwC3AncCL5uDvkjH1CWXXFK6C2o9M6jynAtVmhlUaWZQTdT33Qp6ZebzJtm2BdgyzXPup1pVy5U8\ndFy55pprSndBrWcGVZ5zoUozgyrNDKqJ5uLMAak1vG2NyjODKs+5UKWZQZVmBtVEFgckSZIkSWo5\niwOSJEmSJLWcxQGpD733OpXmnxlUec6FKs0MqjQzqCayOCD1YXR0tHQX1HpmUOU5F6o0M6jSzKCa\nyOKA1IfLL7+8dBfUemZQ5TkXqjQzqNLMoJrI4oAkSZIkSS1ncUCSJEmSpJazOCD1YWRkpHQX1Hpm\nUOU5F6o0M6jSzKCa6MTSHZCOJxs2bODGG28s3Q212gZgYWVw//79R2yzePFilixZMg+90XxwLlRp\nZlClmUE1kcUBqQ9btmwp3QW13pbSHejyDeAE1q1bd8SWp5yyiIMH91sgaAjnQpVmBlWaGVQTWRyQ\n+rB8+fLSXdAcGh4ePuJpgUfzqfj8WkgZ/A7wELALWDZNu/0cPryOkZERiwMN4Vyo0sygSjODaiKL\nA5JaaXh4mDPPXMbhw96nePaWsbCKFpIkSeqXxQFJrTQyMlIXBo70qfeHgDfOT6ckSZKkQrxbgdSH\nnTt3lu6C5tzYp95TPZaW69qkzKDKcy5UaWZQpZlBNZHFAakPQ0NDpbug1jODKs+5UKWZQZVmBtVE\nFgekPuzYsaN0F9R6ZlDlOReqNDOo0sygmsjigCRJkiRJLWdxQJIkSZKklrM4IEmSJElSy1kckPrQ\n6XRKd0GtZwZVnnOhSjODKs0MqoksDkh9uPDCC0t3Qa1nBlWec6FKM4MqzQyqiSwOSH1YtWpV6S6o\n9cygynMuVGlmUKWZQTWRxQFJkiRJklqur+JARPxWRHwuIu6uH5+MiBf2tHlTRNwZEaMRcUtEPKln\n/8kRsSMiRiLi3oi4ISIeOxeDkSRJkiRJ/ev3zIF/BC4FlgMrgI8CH4yIZQARcSnVBbGvBs4G7gNu\njoiTuo5xFfBi4GXAOcAZwPtnMQZp3uzevbt0F9R6ZlDlOReqNDOo0sygmqiv4kBm/v+ZuSczv5yZ\nX8rMy4B/AZ5RN3ktcEVm3pSZXwAuoHrzfx5ARJwKbAAuzsyPZ+ZngfXAMyPi7Dkak3TMDAwMlO6C\nWs8MqjznQpVmBlWaGVQTzXjNgYg4ISLWAouAT0bEUuB04CNjbTLzHuDTwMp601OBE3vaHASGu9pI\nC9b1119fugtqPTOo8pwLVZoZVGlmUE10Yr9PiIifBT4FnALcC7w0Mw9GxEoggUM9TzlEVTQAOA14\noC4aTNVGkiRJkiTNo5mcOXAA+LdUawr8MfDuiHjKnPZqCqtXr6bT6Yx7rFy5csI1P3v37qXT6Ux4\n/saNG9m5c+e4bUNDQ3Q6HUZGRsZt37x5M1u3bh23bXh4mE6nw4EDB8Zt3759O5s2bRq3bXR0lE6n\nw759+8ZtHxgYYP369RP6tmbNGsfhOBzHPI+jcm3P18NAh2qq63YdsKln22jddl/P9gGqK6Z6rQE+\n07Ntb32MI/ULYKhuO5txQLX0S7epxrGHqcfRe63lVOPYOKFt03PlOByH43AcjsNxOA7HMZfjGBgY\nGPceeOnSpaxdu3bCMWYrMnN2B4i4BfgSsA34MnBWZn6+a/+twGcz8+KIeC7wYeAx3WcPRMQdwJWZ\n+fYpvsdyYHBwcJDly5fPqr+SBNX/CFasWAEMUq2xOpX3AOsa0K7k9x4CVuAcLkmSNDce/luWFZk5\nNBfHnPGaAz3HODkzvwLcBTx/bEe9AOHTgU/WmwaB7/e0ORNYQnWpgrSgTVblk+aXGVR5zoUqzQyq\nNDOoJuprzYGI+H3gr6nOV/1R4FeBc4FVdZOrgMsi4kvAHcAVwNeAD0K1QGFE7ATeFhHfplqz4Grg\ntsy8fdajkY6xVatWHbmRdEyZQZXnXKjSzKBKM4Nqon4XJHws8OfA44C7gc8DqzLzowCZuS0iFgHv\nAB4NfAJ4UWY+0HWMi4EHgRuAk6kuat04m0FI8+X8888v3QW1nhlUec6FKs0MqjQzqCbqqziQma86\nijZbgC3T7L8fuKh+SJIkSZKkwuZizQFJkiRJknQcszgg9aH3ViTS/DODKs+5UKWZQZVmBtVEFgek\nPmzbtq10F9R6ZlDlOReqNDOo0sygmsjigNSH6667rnQX1HpmUOU5F6o0M6jSzKCayOKA1IdFixaV\n7oJazwyqPOdClWYGVZoZVBNZHJAkSZIkqeUsDkiSJEmS1HIWB6Q+bNq0qXQX1HpmUOU5F6o0M6jS\nzKCayOKA1IclS5aU7oJazwyqPOdClWYGVZoZVBNZHJD6cNFFF5XuglrPDKo850KVZgZVmhlUE1kc\nkCRJkiSp5U4s3QFJUjvs37//iG0WL17sqZqSJEkFWByQ+nDgwAGe8pSnlO6GWu0AcLxl8BvACaxb\nt+6ILU85ZREHD+63QLDAOReqNDOo0sygmsjLCqQ+XHLJJaW7oNY7HjP4HeAhYBcwOM1jF4cPjzIy\nMlKqozpKzoUqzQyqNDOoJvLMAakP11xzTekuqPWO5wwuA5aX7oTmgHOhSjODKs0Mqok8c0Dqg6c6\nqzwzqPKcC1WaGVRpZlBNZHFAkiRJkqSWszggSZIkSVLLWRyQ+rB169bSXVDrmUGV51yo0sygSjOD\naiKLA1IfRkdHS3dBrWcGVZ5zoUozgyrNDKqJLA5Ifbj88stLd0GtZwZVnnOhSjODKs0MqoksDkiS\nJEmS1HIWByRJkiRJarm+igMR8YaIuD0i7omIQxHxgYj46UnavSki7oyI0Yi4JSKe1LP/5IjYEREj\nEXFvRNwQEY+d7WCkY21kZKR0F9R6ZlDlOReqNDOo0sygmqjfMweeDWwHng78IvDDwN6IeMRYg4i4\nFLgQeDVwNnAfcHNEnNR1nKuAFwMvA84BzgDeP8MxSPNmw4YNpbug1jODKs+5UKWZQZVmBtVEJ/bT\nODNXd38dEa8E/glYAeyrN78WuCIzb6rbXAAcAs4D3hsRp1L9dbs2Mz9et1kP7I+IszPz9pkPRzq2\ntmzZUroLar0tpTsgOReqODOo0sygmmi2aw48GkjgWwARsRQ4HfjIWIPMvAf4NLCy3vRUqqJEd5uD\nwHBXG2lBWr58eekuqPXMoMpzLlRpZlClmUE1UV9nDnSLiKC6PGBfZn6x3nw6VbHgUE/zQ/U+gNOA\nB+qiwVRtJGnGhoeHj3gt4P79++epN5IkSdLCN+PiAPBHwM8Az5yjvkjSrA0PD3Pmmcs4fHi0dFck\nSZKk48aMLiuIiGuA1cBzMvMbXbvuAoLq7IBup9X7xtqcVK89MFWbSa1evZpOpzPusXLlSnbv3j2u\n3d69e+l0OhOev3HjRnbu3Dlu29DQEJ1OZ8KnjJs3b2br1q3jtg0PD9PpdDhw4MC47du3b2fTpk3j\nto2OjtLpdNi3b9+47QMDA6xfv35C39asWeM4joNxnHXWWY0YR1Nej8nG8brXva4uDOwCBuvHLqq1\nTz/cte2K+hnX9hxhGOgAB3q2Xwds6tk2Wrfd17N9AJg4DlgDfKZn2976GL16+wUwBJzF5Hcs6G0/\n1TigOumr21Tj2MPU49jds22qcWwEbu3ZNlS37R3Hn0x49kLJVVN+P+ZqHG9+85sbMY6mvB5tHMfr\nX//6RoyjKa9HG8dx4YUXNmIcTXk9mj6OgYGBce+Bly5dytq1ayccY9Yys68HcA3wj8BPTbH/TuDi\nrq9PBb4L/HLX1/cDL+1qcybwEHD2FMdcDuTg4GBKJb3mNa8p3QUdweDgYAIJgwk5zWPXcdruNTM8\n3kIcS++jeu2c6xc+50KVZgZVmhlUaQ//zcvy7PM9/VSPvi4riIg/As6n+sjnvogYO0Pg7sw8XP/3\nVcBlEfEl4A6qj+e+BnywLkbcExE7gbdFxLeBe4GrgdvSOxVogduxY0fpLqj1zKDKcy5UaWZQpZlB\nNVG/aw78FlV14tae7euBdwNk5raIWAS8g+puBp8AXpSZD3S1vxh4ELgBOJnq3NWN/XZekiRJkiTN\nXl/Fgcw8qjUKMnML09yMOzPvBy6qH5IkSZIkqaAZLUgoSZIkSZKaw+KA1IfJVj2V5pcZVHnOhSrN\nDKo0M6gmsjgg9aH3tjXS/DODKs+5UKWZQZVmBtVEFgekPqxatap0F9R6ZlDlOReqNDOo0sygmsji\ngCRJkiRJLWdxQJIkSZKklrM4IPVh9+7dpbug1jODKs+5UKWZQZVmBtVEFgekPgwMDJTuglrPDKo8\n50KVZgZVmhlUE1kckPpw/fXXl+6CWs8MqjznQpVmBlWaGVQTWRyQJEmSJKnlLA5IkiRJktRyFgck\nSZIkSWo5iwNSH9avX1+6C2o9M6jynAtVmhlUaWZQTWRxQOrDqlWrSndBrWcGVZ5zoUozgyrNDKqJ\nLA5IfTj//PNLd0GtZwZVnnOhSjODKs0MqoksDkiSJEmS1HIWByRJkiRJajmLA1If9u3bV7oLaj0z\nqPKcC1WaGVRpZlBNZHFA6sO2bdtKd0GtZwZVnnOhSjODKs0MqoksDkh9uO6660p3Qa1nBlWec6FK\nM4MqzQyqiSwOSH1YtGhR6S6o9cygynMuVGlmUKWZQTXRiaU7IElSt/379x+xzeLFi1myZMk89EaS\nJKkdLA5IkhaIbwAnsG7duiO2POWURRw8uN8CgSRJ0hzxsgKpD5s2bSrdBbVekzP4HeAhYBcwOM1j\nF4cPjzIyMlKqo63nXKjSzKBKM4Nqor6LAxHx7Ii4MSK+HhEPRURnkjZviog7I2I0Im6JiCf17D85\nInZExEhE3BsRN0TEY2czEGk++CmlymtDBpcBy6d5LCvXNQHOhSrPDKo0M6gmmsmZAz8C/G/gNUD2\n7oyIS4ELgVcDZwP3ATdHxEldza4CXgy8DDgHOAN4/wz6Is2riy66qHQX1HpmUOU5F6o0M6jSzKCa\nqO81BzJzD7AHICJikiavBa7IzJvqNhcAh4DzgPdGxKnABmBtZn68brMe2B8RZ2fm7TMaiSRJkiRJ\nmpE5XZAwIpYCpwMfGduWmfdExKeBlcB7gafW37e7zcGIGK7bWByQNKnh4eEjXmd+NCvdqxm8q4Ek\nSdLcmeu7FZxOdanBoZ7th+p9AKcBD2TmPdO0kRakAwcO8JSnPKV0N1ppeHiYM89cxuHDo6W7UtgB\noO0Z9K4GpTkXqjQzqNLMoJrouLpbwerVq+l0OuMeK1euZPfu3ePa7d27l05nwjqJbNy4kZ07d47b\nNjQ0RKfTmfBp5ObNm9m6deu4bcPDw3Q6HQ4cODBu+/bt2yesWDo6Okqn02Hfvn3jtg8MDLB+/foJ\nfVuzZo3jOA7G8ZznPKcR4zgeX4/3ve99dWFgbCX7a6iWLOldyf7sCc+HIaADTHbWwbU9Xw/XbQ/0\nbL+OiXcKGK3b7uvZPgBMHAesAT7Ts21vfYwj9QuqcTyH2Y0DqmVfuk01jj1MPY7dPdumGsdG4Nae\nbVO9HpMtPTPZOMbuarCa8a/9PqpM/BnddzXYtWtX438/5nscv/3bv92IcTTl9WjjOH7913+9EeNo\nyuvRxnH82q/9WiPG0ZTXo+njGBgYGPceeOnSpaxdu3bCMWYrMiesKXj0T454CDgvM2+sv14KfBk4\nKzM/39XuVuCzmXlxRDz3/7V3/zGWlWcBx7/P4Lrjkto2uRQWQ9LU2s0aI3WnpqIptEEX25IVS5tq\nmaTpNgoRKlnTRo1uXGpSsxCoRSE2EUtJVxPSprAxNTT0h90uv8yMoNEpayPk0rosvWKXRjrjln38\n41Yu6l4AAAvcSURBVN6lwzA/zpk5d96Zc7+fZAJz5rlnnt08+9x7n/ue9wD3A6+ev3ogIp4EPp6Z\nn1jk9+wCpqampti1a9eq85XWqtvt+glkIdPT00xMTNB/07dcHzgETLY4rstL71hQ9XzrmeNGiZsG\nJvC5o3n2QpVmDao0a1Cl/fC1MROZOd3EORtdOZCZTwBPA5eeOTbYgPDNwAODQ1PADxbE7KD/avfB\nJvORmuaTgMqzBlWevVClWYMqzRpUG9XecyAizgZeD5y5U8HrIuJC4NnMfIr+etU/johvAk8Cfwp8\nC7gXXtyg8A7gloj4H+B7wK3AUe9UIEmSJEnS+lvNhoRvAr5Cf+PBBG4eHP80sDczb4yIbcAngVcB\nR4C3Z+b/zTvHPuAF4LPAVvoXtl67qj+BJEmSJElak9qXFWTmP2bmWGaeteBr77yYA5l5fmZuy8zL\nMvObC84xl5kfysxOZr4iM9+Tmc808QeShmnhpiTS+rMGVZ69UKVZgyrNGlQbbaq7FUilPf/8qN9G\nT+VZgyrPXqjSrEGVZg2qjRwOSDXccMMNpVPQyLMGVZ69UKVZgyrNGlQbrWbPAUmSNo2ZmZlKcZ1O\nx92nJUnSyHI4IElqqePAGJOTk5Wix8e38fjjMw4IJEnSSHI4INXQ6/XodDql02iVbrdLr9dbMa7q\np7/t1wOswWq+C5wGPgPsXCF2htnZSXq9nsOBCuyFKs0aVGnWoNrI4YBUw969ezl8+HDpNFqj2+2y\nY8dOZmfd1Ke6vYA1WM9OYFfpJFrFXqjSrEGVZg2qjRwOSDUcOHCgdAqt0uv1BoOBKp/sfgHYP/yk\nNrwDpRNotSorVNybwF6o8qxBlWYNqo0cDkg17Nrlp4/DUeWTXS8r6LMGh6P6/gTuTWAvVHnWoEqz\nBtVGDgckSaq8P4F7E0iSpHZyOCBJ0ovcn0CSJI2msdIJSJvJHXfcUToFjTxrUOXZC1WaNajSrEG1\nkcMBqYbp6enSKWjkWYMqz16o0qxBlWYNqo0cDkg13HbbbaVT0MizBlWevVClWYMqzRpUGzkckCRJ\nkiRpxDkckCRJkiRpxHm3AkmSapqZmVkxptPpeLtDSZK0aTgckGrYs2cPhw8fLp3GptDtdun1esvG\nVHmDpYX2ANZgOceBMSYnJ1eMHB/fxuOPz7RyQGAvVGnWoEqzBtVGDgekGq677rrSKWwK3W6XHTt2\nMjv7fOlUWsgaLOu7wGngM8DOZeJmmJ2d5MiRI+zcuVzc5lxhYC9UadagSrMG1UYOB6Qadu/eXTqF\nTaHX6w0GAyu9gfoCsH99kmoNa3Bj2AnsWubn7V5hYC9UadagSrMG1UYOByQN0UpvoLysQG1Vb4VB\nr9fbVMMBSZLUPg4HJNXiXgJSHSsNyCRJkjYGhwNSDffccw9XXHFF6TSKcS+BjeAeYHRrsK2qDNTm\n5ubYunXrinHrsYfBqPdClWcNqjRrUG3kcECq4eDBg619Iqi6IsC9BEo7iMOBNqm+NwGcBbywYtTW\nreN87nOfZfv27cvGrWWI0OZeqM3BGlRp1qDaqOhwICKuBT4MnAc8BnwoM/+pZE7Scs4555zSKQxF\n/RUB7iVQTjtrcHRV3ZvgzMBtpbgjzM39HpdffvmKv3ktGyG2tRdq87AGVZo1qDYqNhyIiPcCNwO/\nDTwC7APui4g3ZObyH19KapR3F5BKqzpwqxLnRoiSJKm+kisH9gGfzMy7ACLiGuCdwF7gxoJ5Sa1S\nbwNBVwRI7VBtI8Sqm4dWvQShSr8peT5JkrS0IsOBiNgCTAAfO3MsMzMi7gcuKpGTtNlUedF8/Phx\nrrzyPczNfX+dspK0OdTZ6+Dl+xicPHmS6enpl56xRr+psi9C0+eD5ocSdc7ZNAcnkqSmlVo50KG/\ns9KJBcdPADsWiR8Hb4+mtTl27Bj79+/n1KlTy8ZFBFdffTUXXHDBy3529OhRDh069OL3Y2NjnD59\nesXfXTWuamyv1+MjH/kDTp2arXRO+CCw3IvmfwXupX/ZwHL/zo4O/rvecSV/90aLOwocqhBXMsfN\nFrcZchxG3GlW7g0A/8Hc3N0v28dgYmJiifiVzrn4+ZbW3Pm2bNnKTTcdpNPpLBlTt79WOSc0+3xR\nJ8cS+dWJW+05Fz4fDzPHzfZ3MwpxGyHH5WqwZH51YjudjnsnbGLz3huPN3XOyMymzlX9l0ZsB74N\nXJSZD887fhC4ODMvWhD/Pl76aliSJEmSpFF3VWb+bRMnKrVyoEf/fkznLjh+LvD0IvH3AVcBTwJV\nPyqVJEmSJKmNxoHX0n+v3IgiKwcAIuIh4OHMvH7wfQBd4NbMvKlIUpIkSZIkjaCSdyu4BbgzIqb4\n4a0MtwF3FsxJkiRJkqSRU2w4kJl3R0QH+Cj9ywkeBS7LzO+UykmSJEmSpFFU7LICSZIkSZK0MYyV\nTkCSJEmSJJXlcECSJEmSpBG3YYYDEXFtRDwREd+PiIci4udXiH9rRExFxGxEHIuI969XrmqnOjUY\nEZdExOkFXy9ExGvWM2e1S0S8JSIOR8S3BzW1p8Jj7IVqTN0atBeqaRHxhxHxSEQ8FxEnIuLzEfGG\nCo+zF6oRq6lBe6GaFhHXRMRjEXFy8PVARPzqCo9Zcx/cEMOBiHgvcDPwJ8DPAY8B9w02LFws/rXA\n3wNfAi4EPgH8dUT8ynrkq/apW4MDCfwUcN7ga3tmPjPsXNVqZ9PfnPV36NfXsuyFGoJaNThgL1ST\n3gL8BfBm4JeBLcAXI+LHlnqAvVANq12DA/ZCNekp4PeBXcAE8GXg3ojYuVhwU31wQ2xIGBEPAQ9n\n5vWD74P+X8itmXnjIvEHgbdn5s/OO/Z3wCsz8x3rlLZaZBU1eAn9f6Svzszn1jVZjYSIOA1ckZmH\nl4mxF2poKtagvVBDNRjSPwNcnJlfXyLGXqihqViD9kINXUT8N/DhzPzUIj9rpA8WXzkQEVvoT0O+\ndOZY9icW9wMXLfGwXxj8fL77lomXlrTKGgQI4NGI+K+I+GJE/OJwM5Vexl6ojcBeqGF6Ff1PZJ9d\nJsZeqGGqUoNgL9SQRMRYRPwGsA14cImwRvpg8eEA0AHOAk4sOH6C/pKcxZy3RPyPR8TWZtPTCFhN\nDR4HrgauBN5Ff5XBVyPijcNKUlqEvVCl2Qs1NINVfH8OfD0z/32ZUHuhhqJGDdoL1biI+JmI+B4w\nB9wO/HpmfmOJ8Eb64I+sKlNpxGXmMeDYvEMPRcRPAvsAN0GSNBLshRqy24GfBn6pdCIaWZVq0F6o\nIfkG/f0DXgm8G7grIi5eZkCwZhth5UAPeAE4d8Hxc4Gnl3jM00vEP5eZc82mpxGwmhpczCPA65tK\nSqrAXqiNyF6oNYuIvwTeAbw1M4+vEG4vVONq1uBi7IVak8z8QWb+Z2b+c2b+Ef0N069fIryRPlh8\nOJCZp4Ap4NIzxwZLeC4FHljiYQ/Ojx/YzdLXYEhLWmUNLuaN9JeVSevFXqiNyF6oNRm8Kfs14G2Z\n2a3wEHuhGrWKGlyMvVBNGwOWukSgkT64US4ruAW4MyKm6E/Z9tHfcOFOgIj4M+D8zDyzLOevgGsH\nuzL+Df2/iHfTn+5Jq1GrBiPieuAJ4N+AceC3gLcB3jZJqxYRZ9P/lCEGh14XERcCz2bmU/ZCDVvd\nGrQXqmkRcTvwm8Ae4H8j4swnYSczc3YQ8zHgJ+yFGobV1KC9UE0b1Ng/AF3gFcBVwCX03/AP7f3x\nhhgOZObdg9uEfJT+8odHgcsy8zuDkPOAC+bFPxkR7wQ+Dvwu8C3gg5m5cIdGqZK6NQj8KHAzcD7w\nPPAvwKWZ+bX1y1ot9CbgK/R3RU76NQbwaWAv9kINX60axF6o5l1Dv/a+uuD4B4C7Bv+/HXuhhqd2\nDWIvVPNeQ/+5dztwkn5N7c7MLw9+PpTXhNG/Y5skSZIkSRpVxfcckCRJkiRJZTkckCRJkiRpxDkc\nkCRJkiRpxDkckCRJkiRpxDkckCRJkiRpxDkckCRJkiRpxDkckCRJkiRpxDkckCRJkiRpxDkckCRJ\nkiRpxDkckCRJkiRpxDkckCRJkiRpxP0/rel3cIhMRzMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x104588c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_conditional = df #[df.et_month == 8]\n",
    "lognormal_et_accuracy = df_conditional.accuracy_ratio[df_conditional.accuracy_ratio < 3]\n",
    "\n",
    "## Landscape Error\n",
    "print \"Mean monthly error: \", np.round(lognormal_et_accuracy.mean(), 3)\n",
    "print \"Std. dev. of monthly error: \", np.round(lognormal_et_accuracy.std(), 3)\n",
    "print \"5/95 percentiles: \", np.percentile(lognormal_et_accuracy, [5, 95])\n",
    "lognormal_et_accuracy.hist(bins = 70)\n",
    "plt.title(\"Histogram of ET accuracy ratio\", fontsize=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Overall Error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The equation for the the estimated outdoor efficiency targets used in this methodology is\n",
    "\n",
    "$$Outdoor Target = Outdoor Standard * \\hat{LA} * \\hat{ET} * Conversion Factor$$\n",
    "\n",
    "Combining the outdoor standard and the conversion factor into a constant $C$, and substituting in our error models gives\n",
    "\n",
    "$$Outdoor Target = C * (LA *  \\epsilon_{classification} * \\epsilon_{boundary} * \\epsilon_{landuse}) * (ET* \\epsilon_{et})$$\n",
    "\n",
    "So, we can calculate the overall error for the outdoor target as the product of the individual error components. This is done using a monte carlo simulation, where random samples are drawn from our prior distributions, or from the empirical error distribution in the case of ET. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with pm.Model() as model:\n",
    "    et_prior = pm.Lognormal(\"et_prior\", mu=0, sd=0.1, observed=lognormal_et_accuracy)\n",
    "#     has_boundary_error = pm.Bernoulli(\"has_boundary_error\", p=0.5)\n",
    "    boundary_prior = pm.Lognormal(\"boundary_prior\", mu=0, sd=0.1)\n",
    "    landuse_prior = pm.Lognormal(\"landuse_prior\", mu=0, sd=0.1)\n",
    "    classification_prior = pm.Lognormal(\"classification_prior\", mu=0, sd=2*.0387755102041)\n",
    "    outdoor_error_dist = pm.Deterministic(\"outdoor_error_dist\", et_prior*boundary_prior*landuse_prior*classification_prior)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "class_samples = [random.choice(list(lognormal_et_accuracy))*boundary_prior.random()*landuse_prior.random()*classification_prior.random() for i in xrange(20000)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean monthly error:  [ 1.028]\n",
      "Std. dev. of monthly error:  0.214\n",
      "5/95 percentiles:  [array([ 0.72775317]) array([ 1.4040262])]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABAMAAAF0CAYAAABWofDwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3X2cXGV9///3JyTcBDSA3IkSQAQMKOCud1G5KxiUmxG0\nGLEqJN5gG1qLNunNrzbR9qvfgIptUq1fDYK2JijFNKYWInJnRAF35EZgjRFCgGBwA5SbTUiye/3+\nuM7CZHZmZ+bsnr12P/N6Ph7zSPbMOWc+Z97nmptrzrmOhRAEAAAAAADax4TUBQAAAAAAgNFFZwAA\nAAAAAG2GzgAAAAAAANoMnQEAAAAAALQZOgMAAAAAAGgzdAYAAAAAANBm6AwAAAAAAKDN0BkAAAAA\nAECboTMAAAAAAIA2Q2cAMIaZ2Uwzu8bM1pvZZjN7wsx+ZWYLzeyg1PUVzcxuMrN+Mzuhavr8bPo/\npKqtkpmdmNVzwyg+5hXZY354tB4Tw2dm52e5XV41fdT3oaGMlzolycwOzmp6IHUtjZjZLDO7w8ye\nzWruM7OpqesaT3jtaw/Z+/zljecEMBx0BgBjkJm93Mxuk7RUUknSY5J+IOkWSQdKmitpjZn9Wboq\nR0XIbvXuc6neF7EqQz03GJ+GnWmT+86o1TMSzGxdtk1DfWkeE7UOxczOkLRE0lGSfiLpiuz2bLqq\nRtcI7Z9jPutmNLlfu2BmC8ZSBz6AF01MXQCAHZnZnpJWSzpEUpekD4UQuivunyDpk5IWSlpkZhNC\nCItT1JqYpS4gsb+R9AXFjiKMf7dJmiapN3UhmWsk/VzS/6YuRI2//D2q+NxtG51ycvtjxe348xAC\nv3jCRadGk5raVjN7s6S/lTRd0sviJDtL8X3ul5K+GUK4tchCgXZDZwAw9vyrpEMl/U7SKSGEpyvv\nDCH0S7rMzLZk837RzH4cQvjN6JeKgjTs6AghbJS0cRRqwSgIIWyRtGYEVjUinWQhhGckPTMS6ypa\nCGG7Rua5K9rAL8Brk1aRVrt34rYrU4PszeyPJF0raSfF9vx7SZMl3SPpdZLOl/SIJDoDgBHEaQLA\nGGJmh0qaqdiD/lfVHQGVQghfk3SXpEmS5lWs47vZ4Xjz6i1rZmdm83TVuO9wM/u6ma3Nxil4ysxu\nNrM/qbOuF87rN7PjzeyHZvZ4di7sh7N59jCzj5nZf5rZmux82WfN7G4z+yczm9Lsc5RX5bnOZrab\nmX3OzO4zs+eqzzU2s13N7NNm9nMzezJ7HrqzsRr2bvFxTzGzRdlYD38wsy1m9rCZLTOzN9SYf52k\nyxX3gQuymgduN1TMN+i82dHOvoXtv8bMNpjZ82a2Mfv7LXXm7zezvuz/s8zs1qyOFw6nbWafq1jf\n+83sJ2a2KXvu15nZEjM7vM7jr8vWM9XM3l2x7KCxK4bY5p3M7C/N7J7seXzczK42s9cOsUzdc/HN\nrMPMrsr2m+fN7H/N7HfZOkuVtau5faeZNtvwcO6sHX3ezH6bbeejZvZNMzuwle2rmOeF7CtrUPwS\nbZLWVW3TCdl8Q44ZYGavyNrgmor9erWZfdzikVbV87+w7WY22cy+kG3jFjN7LGt7g7ZxiO36VrYd\nJ2fbcVPFNlSPyfBGM/te9lwOtJcVZnZqnXW/8DpgZkdn+8kGM9tuLRySneNxa47nUnH/oHFdmt0/\ns3n3MrOvZO1xi5k9lGW4V4Pt2MnMPmEvvm5sznL/56EyG+Y+MlDr2qzWIcfTaHa/zuZ9T9am7rE4\nZtBmM3vA4mvYEXXW3/Q+YWazzeyXFt8He8zsR2Y23Rq0V4unMn7ZXnwPfdrMbjezOWa2U9W8/ZIG\nHnfgdIFa+//nFTsCPh9CmKZ4dNJPQwjvCSEcLukNkn481HMLoHUcGQCMLWcpdtI9JemHTcz/HUlf\nzJYb8C1J71fsRb+kznKzFD+QLamcaGbnSrpS0i6SuiX9t6Qpkt4s6TtmdnII4aNV6xo4/O99kj4h\n6X7FN+y9JT2fzXOspK9LelzSbxQP99tLUqekv5N0rpm9JYTwZBPbPFy7SrpJ8bDiWyTdmdUqKX7I\nkXSdpNdK2iTpdsVfSDsUx2o418xODCE83OTj/ZukV0q6V/H0j+2SXiPpXEnvMbOZIYQfVMz/PUlv\nkfR2xV8QV1fc113x/1qHXY529kMysy9K+pSkPsXMb1H8AFySdJaZfTSEcGWdZf9F0p9J+pmklYpH\nywxsbzP7nMzsSkkfUjx8/BbF/a8jew5mmtl7Qgirqh564DH+StJFku6Q9D+KY3X0qQEzM0lXS3p3\nVstNkp5UfB5vV/wy1DQzO0XSjxTfr+9S/FVsJ0mvkHS64uvFimz2VvedIZ+/BnZWPO/9dYrb2JU9\n7mxJp5vZ8SGE3zW/pTWtVTyn/lzFXwj/Uy+eXx8Ufzkckpm9UfHXxj0lrVcce2WKpJMkvVXS2WZW\nyo4uqBSy+W6VdJCknyr+Qjld0oclnWBmx2ZHUDTy02x975K0n+Lry0DtL2RkZh+T9DXFL4i/knSj\npIMlnSHpTDNbEEL4XI06g6S3Kb7GbpB0s6Td1OSRHcN83KFU39/U/mlm+2X3vVrSE4rvhRMkfUDS\nOxVfS2ttx86Kr1unSNqcbcfTijn/uaTzzGxGCOHOquWGs4/so/jaNkUx519K2lr/KZHU2n59laQt\nku5TbG8TFd+bZkl6n5m9I4Twixp1NdwnzOyriu2/L6v9McX2fLOkf65XfNZZsTzb5nWSVim+b7xJ\n0iLFfebMEMLA6+UVko7LbndmtwGV+8Cx2b8137tCCL+qVxOAYQghcOPGbYzcFL+M9Uu6vsn5j8/m\n75N0cDbNFN+g+yS9qcYyL1P8wL9Z0l4V01+bTXtO0rurljlI8YtIn6QPVt13Y0UNF9ap8xWSTq4x\nfVfFL7B9khbVuP/G7L4TqqbPzx7zH1p4bk/MlulX/MC7b535VmeP+XVJu1dMn6D4IWVQPhXrvqHG\n+kqSptSZvlXxC+ouVfedn63v8iG2Z+B5+3DFtFHNvsHz/bFsG7olHV1139sVz0XfLOmwqvsG9qUn\nJb2xzrqb2ec+kc2zUdLrqu77h+y+TZJeVnXfg9l9WyWd0ez2Viw/J1t+g6QjqvafxRV1X161XM19\nSNIN2fzvr/FYL6nOucl9p5nnr+Z6qtrRbyS9ouK+nSV9P7vvZ822kersa0x/MKt1ap3lDs6WfaBq\n+s4V7WGxpJ0q7jtE0gPZff9YZ9v7Fb9cVr4OTJFUzpb76xb3jZqvZxVtcKtiZ+EHqu47TfELYZ/i\nqWPVrwMDWf5Tjv017+PW3Zbs/pqv0U3unwP70I2S9qiYvqfiOBYD2/vhquX+b8V+eVDF9J0k/b/s\nvrWSJo7wPnJd5T7SwnM/5H6dzXOupN1qTB94fbu7xn0N9wnF959+xdfhN1fd95cVy1e/Hu0vqSfb\nXz5edd9ekq7Plvv7ZvaHqnn+kC372opl6u4n3LhxG5kbpwkAY8u+ij36zZ4L/njVsgohBMVOBVP8\n9aDaBxVPLfivsOMv8X+v+MHo/wsh/FflAiH+Cj47W+df1KnlJyGEr9e6I4TwaAjhxhrTtyj++tun\n+KFnNARJc0IIf6i+w8xOU/wl6FeS/jSE8FxFrf2S/lrSryWdbGZHNfVgIawIIQwahC2EsELxQ+/L\nFA8fHraE2e8g+3V8geJz/f4Qwg6/5IUQVkv6R8Vfky6ss5pLQwh3NHiouvuc4i/7QdJnQwj3VD3+\n5yTdrfjl4mM1lg2Srggh/HeDx6/lL7Pl54cQXjiPPdt/Pq0mfs2usl/27/8MKjKEZ0IIt+eoccBQ\nz18jQdKnQwiPVtSzVbE990p6i9U5FWQUnat4JMoGSReHF3+pVAhhneI+YpL+PPtVudqzkmZVvQ78\nr+KXTpNU8xD6nP5S8Vffa0II3628I4RwneKXWVM8OqmWNZI+k+BxR5SZvVLSOYpfHD8RQnjhSgsh\nhKcUvwTXWm4XxX0vKGb9cMVyfYqD7m5UPMLojysWHe4+slXxS/FzNe4bthDC90MIm2tM/zfFjpGj\nzew1dRYfap/4pOJz9S8hhNuq1v0VxSOiarlY8QiixSGE/1e13JOKR81sVzyqqlX/qfhcLzezDyke\ndQGgYHQGAD5dofhGPzP7kFRp4DDxFw5Xzr68vTP783t11llW/HD8+hofioLiG/mQsnMR55nZ4uxc\ny28pHp66VdK+NgpjB0h6PNQfjfgMxW25JvvytoPsy/Yt2Z9vbfYBs/MrP2pmXzSzb1g8h/hbko7O\nZjmyhfobuUKjm30tr5f0ckm/C1WH5Fa4Ofu33vPYaH+qu8+Z2SskvSr789t1lv+W4gfPeh0xDffn\nGo97oKTDsj//o/r+EMLzis9xK4Oo3Z7N/10ze1v1+bjD0FSbHcJTIYSVg1YaO9muzf48aRjrHwkn\nKW7nshDCoCsNhBCuUTwC5SWKpyxV+2UI4fEa0+/P/n3FCNUpxSMnBjrzahk4ref4rM1WCpKWZ69P\no/m4RThB8bNpOdQYFDeEcJdiR161N0jaQ9ITIYQf1Vhus6RlGtzmT9Lw9pFfhRAearBNw2Jmh2Xn\n4l+WjR8w8P6xfzZLrfePuvtE9hoyPfvzu9X3V0yvlffp2bprvleEEDZI+q3i+/mr625UbZ9WPP3g\nUMX9cY7iqXQ/MrNPmdm+La4PQBMYMwAYW3oU34D3bzRjZr+K/7/wS3cI4UEzu1nxg945ih+CZGbH\nSTpG8VJclQPxvEzSSxXf5B9p8JkvZPNXX9JuXb0FsjfxaxTPYaz1gdWy6S9V8ZcyWzfEfa/Kavkn\nM/unIeYLyo7EaMTM5iuOizDU6+1Lm1lXMxJlX23gi/hh2eBRQ62v3vO4rsFjDDXPwJe0TZW/LFb5\nXdW8eR6/2iuzf3tCCPUuEfhgi+v8W8XzeN+peM75ZjMrK56n/x+h4rKjOawraNkHFdvRK4eYZzQM\nZDvUc/6gYudVrf1gfZ1lBgZ23TVnXbU0qnVgf91VsQ32VN2/LtHjjrSBfaZRZq+rmtZM1rXa/HD3\nkXVDLDcsFgcu/FdJH28wa733j3V1pu+jmGcYYp560wde21c38V6xr1q4ckZ2dMV7zOwYxfeuP1Yc\n2+c0xde/z5rZx0IIy5pdJ4DG6AwAxpYuxUO5O8xsQq1fp6u8Kft3U41fJy5X/NXjAmVfCBUP9w6S\nrqz6xaDyKKErmqiz1iBjgw5lrLBEsSPgZ4rnAd4t6cmBQzLN7FFJB2h0Ljs1VJ0TFJ+f1Xrxg2M9\nNQexqmRm71Hc3mcUB9K7UdKG7PQImdn/kfQ3GvntHu3sqw2s8/eK59MOpeaXi+xX9EaGynK4ilx3\n00K8hOQbs0G7TlVsR29WPKLi78zsb0IIl+ZcfdHb2PR+PUq/Oreq0evvWDIm9tcK7XLkaZHP+18q\nnkb1mOLh+T+XtDE7HUdm9h+KA8bWazvDqa3eUSYDuX5fcYyZoWzK9cAh3C3p7uwl4WBJCxVPUftT\nSVeY2W0hhFY7VQHUQWcAMLb8UNKXFAeperfiqMZD+ZDim3atKw/8p+KASKdkh00/Lum87L4rqubt\nUfzgsKviJQ2fyFN8LWY2WfEXzT5Jp4eq0bez+w8YqccbpoHzTP8rhPDlEVjfuYr5/F0IYUmN+2te\n3m4EpM5+4HncFEKYPQLra9XAeewvM7M96hwd8KqqeUfycfcxs8l1jg44JM+KQwi3KDtFJTtV4wJJ\nX5X0eTO7OsGH40Ma3BcUrwk+YGCE9ZfUWebg4Zc0yEAerxpinkOr5k3lUcU6X6U4cny1gW3YojjC\nfurHbSbPPKctDORwyBDz1LpvYLlDa9w3oFabH8v7yMD7x8frjF+S9/1jk2Kn7s6KOdU6uuiQOss+\nrHiVh4UhhHLOx29JdrrInGw8iTMVT+dbPBqPDbSDdum5BcaFEMIDevGc4kvNrO7h42Y2R/Gw7+2K\nlxesXtdmxcsSTVAc1OcsxcM8V4cQ1lbN268XDx1/3/C3ZAdTFEdzfrq6IyDzIY3OEQHN+B/FWkZq\nMMOBSxYOOtw4O3XiHXWWG/ignavDdgxkf4diJ8NRZjZthNbZtGxQu4EjOy6oM9sFih+0h7wmeI7H\nHbjW/Qeq78++xA98wB/O42zNBu+6WzHjYyruHta+04I9zeyM6olmto9eHIPipoq7XvjSZWa1ajtz\niMfKu003KbbnmbXGujCzcxRHQH9G8aislG5SrPWCOvd/JPv3liaOGBuNxx3Ic1D7NrPdVH8sjkZZ\n3qLYPjrM7Iga6z5WO+7vA36pOK7J3mY2aF8ys10Vf0UPikdoDbhJ6faRRs/FUO8fRyteqq9lIV4i\n8efZn4NepxpMH3iPbPW9YiRel9Zl/47k6TlA26MzABh75ii+6R0q6cbqUevNbCcz+5Skryh+sJkX\nQrh/0Fqiy/XiB72Bw8S/VWfezypej/2LZvbhWoftmtnR2YejVmxUHIBpTzP7YNX63iLp88r35WhY\nX6jq+C/FL7JvygZpGjSasZntZWYXZudzNnK/4vP/cTObVLGOKYoD29Xr7Bn4RbWpKxbUkSz77MPm\nZxXfY5ab2dtqrG+CmZ1sZm9uZp05fFFx+z+TnYNa+difUfwg/aSkb47w434le9wFZvbCwF7Z/vIl\nSQe2sjIz+7SZHVRj+mv04i+DlacIjcS+06wvZUeeDNS0s+LRCrtLui2EMPCFQyGE9YoDi+2peFUO\nVSx3kuL+Us/ANh09xDy1fF/xi9SBki6rHHzRzA5VzGNgRPVG14Yv2j8rduyebWZ/UnmHmc1QPG88\nqEbHb6LHvV5xP5+TDZw5sMxkSd9QvCRpLUPun9lVAH6g2IH8NTN74cgDM9tLcf+qtdzziufXm+J+\nObViuYmS/kXxCLQHJF1dsWjKfaTRfj3w/jGn8jXZzF6u+P4xnMFE/yVb919Uvwab2ScVT0Gs9R57\nqaSnJH0qG9RvUvUMZnZI9b6kJtqwmX3FzGoeJZh1DL0/+/PnteYBkFMYA9c35MaN2443xZHYb1M8\ntL4v+/93FUfa3ah4LutmSRc1sa579eL1kJ+WNHmIed+r+AtIv+IHpGslfUfxWtvrs+nfrVrmRg1x\nvelsnk9m8/QrvpH/h6SfKn4IvUJ1rrdcb91q4prFNWo4UQ2uc17x3Hdlj/uM4vgB31X8AFlW/NLc\nJ2nnRutWPNRyUzb/w4ofPJcrfgl9RPFDc61rcU/K7u/Larkim/evKub5lmpca3s0s2/iOf+/Fbnf\no/gh/7uKv8Y/kd1Xfa3qmteab3Wfy+a7Iptvq+LRD/+h+AG7X/FXxBk1lqm5L7awzZZl3Kd4aPX/\nZNv8O8VzbBerxnXWh9iHnsym36d4+se/Z8/f1uwxqtfTzL7TTJs9v0GdqyXdmj2PKxTHpngku2+D\npFfXWOc5im2+X7EtXaXY+bZdL7bpQdkrXjJuYB++Otueb0g6PLv/4Oz+B2os+wa9eP3yByUtlbRS\n8fKHfdn+PbGZba+4v+7jDWe/VbzM5bZs3b/Msv6pXnwf+PsayzR8HWiirjyPO1HxShd92T76w+y5\n3Kj4ejGc17b9FS+L16d4hNHVivv+E9n05dm6P1y17p0lrcqWey7Leali53p/Vttxo7GPNPm8N9qv\n36T4Pt+XbfcyST/Ktu3u7Dmp9Tw0tU8oXsmnL8v+BsXXx7sVX1u+mK372hrLvT17LvsUx4W5XvG9\nYoVih1+/pFurltlP8f2lL9u3Ls+29YKKeforarlU0k+yepYrvpb2Sfr3vM83N27cat+SF8CNG7f6\nN8VD8a5R/CK5OfvQdafigDpNfVlRvE7ywIe6JU3MPzX7IHBX9iHlOcVfU36SrevQqvlvVPww3+iL\n2VnZh4BNilcMuE3ZF8HsA9j26m2qt27FLw59kj7TwnN5YrbMT5qYd5LiB+TrFc+3f15xEKcuxV/S\nTml23dnz+e1sG3uz53Kx4kjLdbdD8ZezH2QftgY6IG6ouL+ZzoBCs2/yeX9Ltv0PZNv/lOIX8v9U\nPGphStX8/ZK2N1hnU/tcNu/MrP5N2QfKdYpHAxxeZ/6a+2KL2zxBcfCve7Jtfjzb3tcpfokYlEe9\nfUhxrIdvZpn8oWIf+qGks+o8fqN9p+Hz10ydknZT7PBZq/j6tCGr9RVDrPedioeCP5PtYz+T9N6h\nslfsYJmn+MXguYp9+oTs/oOzv39X5zFfofhL6G+zOp9S7Mz4mKQJzW57xf1DPt5w9ltJb1TsJHlU\n8XXnccUjlv6ozvzD7gzI87jZMi9VfD18KGtb6xV/ud9Hw3hty+bZS/Eom4eyzB5SfN3ce6htVmx7\nF2b71VPZsmskXSbp5UNsy4juI00+50Pu19k8R2fP1SPZPN2KR9PtXu95aGWfUDxq7I5s3ZsUOxve\nqjiQcb/qfPnOMl6QLTvwPD+k+B7/GUlH11jmbYoDyvZU5H55xf1vlfRlxR8MHtaLHZ4vZDGcfZwb\nN261bxZCEAAAAACY2eWKHR6fDiF8JVEN8yUdHNIMQgu0jZbGDDCzvzWz283saTPbaGY/qDXAS9Uy\nJ5pZf9Wtz8z2G2o5AAAAACPPzI7KxnionGZm9jHFjoAtiqdMAHCs1VE9j5e0SPGcsomSviBplZlN\nC3H06nqCpCMUDwuME0J4vMXHBgAAADB8cyW9z8x+pXh6yO6Kp3Acong6y5+GEDamKw/AaBjWaQLZ\nSNuPK57ftLrOPCcqDgayVwjh6dwPBgAAAGDYzOw0xTEROhXHAJio+Jl+taR/DiHcnrA8AKNkuNch\n3lPxV/8nGsxnku7MrvP6a0kLQgi3DvOxAQAAALQohHCd4oB+ANpY7iMDsmue/lDSS0IIJw4x3xGK\now//UtIuir2QH5L0phDCnXWWeZmk0xRHfd6Sq0AAAAAAAMa/XRVP47kuhLBppFY6nM6Aryl+YX9b\nCOGxFpe9SdJDIYTz69z/AcXrnQIAAAAAAOlPQgjfHamV5TpNwMwWSzpd0vGtdgRkble83mg96yTp\n3//93zVt2rQcq8d4cPHFF+uyyy5LXQYKRMb+kbF/ZOwfGftHxv6RsW/333+/PvjBD0rZ9+SR0nJn\nQNYR8G5JJ4YQ1ud83OMkDdWJsEWSpk2bpo6OjpwPgbFuypQp5OscGftHxv6RsX9k7B8Z+0fGbWNE\nT6FvqTPAzL4q6TxJJUnPmdn+2V3/G0LYks3zeUmvGDgFwMw+KelBSfcqnuvwMUknS3rHiGwBxq3f\n//73qUtAwcjYPzL2j4z9I2P/yNg/MkYerR4Z8AnFqwfcVDV9lqRvZ/9/uaSDKu7bWdKXJB0oqVfS\n3ZJOCSHc0mqx8OXRRx9NXQIKRsb+kbF/ZOwfGftHxv6RMfJoqTMghDChiXlmVf19qaRLW6wLbaCz\nszN1CSgYGftHxv6RsX9k7B8Z+0fGyKPhl3ugKOedd17qElAwMvaPjP0jY//I2D8y9o+MkUfuSwsW\nycw6JHV1dXUxEAYAAAAAoG2Vy+WBoz86QwjlkVovRwYAAAAAANBm6AxAMrNmzWo8E8Y1MvaPjP0j\nY//I2D8y9o+MkQedAUhmxowZqUtAwcjYPzL2j4z9I2P/yNg/MkYejBkAAAAAAMAYxZgBAAAAAABg\nRNAZAAAAAABAm6EzAMmsXr06dQkoGBn7R8b+kbF/ZOwfGftHxsiDzgAkc8kll6QuAQUjY//I2D8y\n9o+M/SNj/8gYeTCAIJLp7e3V5MmTU5eBApGxf2TsHxn7R8b+kbF/ZOwbAwjCHV6w/CNj/8jYPzL2\nj4z9I2P/yBh50BkAAAAAAECboTMAAAAAAIA2Q2cAkpk7d27qElAwMvaPjP0jY//I2D8y9o+MkQed\nAUhm6tSpqUtAwcjYPzL2j4z9I2P/yNg/MkYeXE0AAAAAAIAxiqsJAAAAAACAEUFnAAAAAAAAbYbO\nACTT3d2dugQUjIz9I2P/yNg/MvaPjP0jY+RBZwCSmTdvXuoSUDAy9o+M/SNj/8jYPzL2j4yRBwMI\nIpn169cz8qlzZOwfGftHxv6RsX9k7B8Z+8YAgnCHFyz/yNg/MvaPjP0jY//I2D8yRh50BgAAAAAA\n0GboDAAAAAAAoM3QGYBkFi5cmLoEFIyM/SNj/8jYPzL2j4z9I2PkQWcAkunt7U1dAgpGxv6RsX9k\n7B8Z+0fG/pEx8uBqAgAAAAAAjFFcTQAAAAAAAIwIOgMAAAAAAGgzdAYgmZ6entQloGBk7B8Z+0fG\n/pGxf2TsHxkjDzoDkMzs2bNTl4CCkbF/ZOwfGftHxv6RsX9kjDzoDEAyCxYsSF0CCkbG/pGxf2Ts\nHxn7R8b+kTHy4GoCAAAAAACMUVxNAAAAAAAAjAg6AwAAAAAAaDN0BiCZJUuWpC4BBSNj/8jYPzL2\nj4z9I2P/yBh50BmAZMrlETvdBWMUGftHxv6RsX9k7B8Z+0fGyIMBBAEAAAAAGKMYQBAAAAAAAIwI\nOgMAAAAAAGgzdAYAAAAAANBm6AxAMqVSKXUJKBgZ+0fG/pGxf2TsHxn7R8bIg84AJHPRRRelLgEF\nI2P/yNg/MvaPjP0jY//IGHlwNQEALdu+fbv6+vqamnfChAmaNGlSwRUBAAAAPnE1AQBjwjPPPKNX\nv3qadt1116Zu++33cj344IOpywYAAABQYWLqAgCk19/frwULFqi7u7vhvM8995weemitpL+WdHSD\nuXv01FOf0tq1a3XooYeORKkAAAAARgCdAUhm+fLlOvvss1OXAUnr1q3TP/7jP8qsU9LeDec3+6BC\n+Iyk3RvM+Y2RKA9jGO3YPzL2j4z9I2P/yBh5cJoAklm6dGnqElAlhEsVwqombt9R444ASVpRdMlI\njHbsHxn7R8b+kbF/ZIw8GEAQgB544AEddthhkm6QdPIIrnm9pIO1atUqveMd7xjB9QIAAADtgQEE\nAQAAAADAiGipM8DM/tbMbjezp81so5n9wMyOaGK5k8ysy8y2mNkaMzs/f8kAAAAAAGA4Wj0y4HhJ\niyS9WdKpkiZJWmVmu9VbwMwOkbRS0k8kHSvpnyV908w4ZhgAAAAAgARa6gwIIZweQvhOCOH+EMI9\nki6QNFW00vcLAAAgAElEQVRS5xCL/amkB0II80IIvwkh/KukqyVdnLdo+DBr1qzUJaBwf5W6ABSM\nduwfGftHxv6RsX9kjDyGO2bAnpKCpCeGmOctkq6vmnadpOnDfGyMczNmzEhdAgp3fOoCUDDasX9k\n7B8Z+0fG/pEx8piYd0EzM0lfkbQ6hHDfELMeIGlj1bSNkl5qZruEEJ7PWwPGt/POOy91CSjcuyX9\nhe666y5NmjSp4dx77723jjnmmOLLwoihHftHxv6RsX9k7B8ZI4/hHBnwVUlHSXr/CNUyyOmnn65S\nqbTDbfr06Vq+fPkO861atUqlUmnQ8nPmzNGSJUt2mFYul1UqldTT07PD9Pnz52vhwoU7TFu/fr1K\npZK6u7t3mL5o0SLNnTt3h2m9vb0qlUpavXr1DtOXLl1a87CdmTNnsh1sx5jbDmmNpJKknqrp8yUt\nrJq2Ppu3u2r6IkkD2/FSTZy4v+bOnauTTz654e31r3+9yuUyebAdbAfbwXawHWwH28F2sB1tuR2d\nnZ2DvgOfeeaZg5YfCRZCaH0hs8WSzpJ0fAhhfYN5b5bUFUL4VMW0CyRdFkLYq84yHZK6urq61NHR\n0XJ9AFrzwAMP6LDDDpN0g6STR3jtT0h6son5HpN0vK699lqddtppI1wDAAAAMD6Vy2V1dnZKUmcI\noTxS6235yICsI+Ddkk5u1BGQ+bmkU6qmzcimo41V95rBo9WS9pZ0WBO3Q9KUiGGhHftHxv6RsX9k\n7B8ZI4+WOgPM7KuS/kTSByQ9Z2b7Z7ddK+b5vJldWbHYv0l6lZktNLMjzezPJP2xpC+PQP0Yxy65\n5JLUJaBwZOwd7dg/MvaPjP0jY//IGHm0emTAJyS9VNJNkjZU3N5XMc/LJR008EcIYZ2kMySdKulO\nxUsKfiSEUH2FAbSZZcuWpS4BhSNj72jH/pGxf2TsHxn7R8bIo6WrCYQQGnYehBAGjY4QQrhFUmcr\njwX/Jk+enLoEFI6MvaMd+0fG/pGxf2TsHxkjj+FcTQAAAAAAAIxDLR0ZAGD86O3t1bve9W6tWfOb\nhvNu3749+x8vCQAAAEA74MgAJFN9bU6MrA0bNuiWW67X739/sn7/+1lD3np6Pibpq5LePsJVkLF3\ntGP/yNg/MvaPjP0jY+TBz4BIZurUqalLaBOzJJ2U6LHJ2DvasX9k7B8Z+0fG/pEx8rAQQuoaBjGz\nDkldXV1d6ujoSF0OMC6tXbtWhx9+uKQbla4zoBWPSDpI1157rU477bTUxQAAAABjQrlcVmdnpyR1\nhhDKI7VeThMAAAAAAKDN0BkAAAAAAECboTMAyXR3d6cuAYUjY+9ox/6RsX9k7B8Z+0fGyIPOACQz\nb9681CWgcGTsHe3YPzL2j4z9I2P/yBh50BmAZBYvXpy6BBSOjL2jHftHxv6RsX9k7B8ZIw86A5AM\nl0BpB2TsHe3YPzL2j4z9I2P/yBh50BkAAAAAAECboTMAAAAAAIA2Q2cAklm4cGHqElA4MvaOduwf\nGftHxv6RsX9kjDzoDEAyvb29qUtA4cjYO9qxf2TsHxn7R8b+kTHysBBC6hoGMbMOSV1dXV3q6OhI\nXQ4wLq1du1aHH364pBslnZS4mmY8IukgXXvttTrttNNSFwMAAACMCeVyWZ2dnZLUGUIoj9R6OTIA\nAAAAAIA2Q2cAAAAAAABths4AJNPT05O6BBSOjL2jHftHxv6RsX9k7B8ZIw86A5DM7NmzU5eAwpGx\nd7Rj/8jYPzL2j4z9I2PkQWcAklmwYEHqElC4BakLQMFox/6RsX9k7B8Z+0fGyIPOACTDlSLaARl7\nRzv2j4z9I2P/yNg/MkYedAYAAAAAANBmJqYuAAAqXXnllVq9enXD+fbdd19ddNFFmjCBPk0AAACg\nVXQGIJklS5boIx/5SOoyUKglkprNeF9NmnS6rr76Vkm3Nph3u7Zte1TTpk3TO97xjuGViGGhHftH\nxv6RsX9k7B8ZIw9+UkMy5XI5dQkoXCsZ76Jt2/5b27ata+J2myRp+/btxZSNptGO/SNj/8jYPzL2\nj4yRh4UQUtcwiJl1SOrq6upiMAwgp7Vr1+rwww+XdKOkkxJXM9IelfRK/ehHP9K73vWu1MUAAAAA\nhSmXy+rs7JSkzhDCiPX8cGQAAAAAAABths4AAAAAAADaDJ0BAAAAAAC0GToDkEypVEpdAgpHxt7R\njv0jY//I2D8y9o+MkQedAUjmoosuSl0CCkfG3tGO/SNj/8jYPzL2j4yRB50BSGbGjBmpS0DhyNg7\n2rF/ZOwfGftHxv6RMfKYmLoAAK2599579fDDDzecb8OGDaNQDQAAAIDxiM4AYBzZsGGDOjreoK1b\ntzQ1/047vUR9fQcVXBUAAACA8YbTBJDM8uXLU5cw7jz77LNZR8D3JT3c8NbXt07SYWmKlSSRsXe0\nY//I2D8y9o+M/SNj5EFnAJJZunRp6hLGsf0kvbKJ296pCsyQsXe0Y//I2D8y9o+M/SNj5EFnAJK5\n6qqrUpeAwpGxd7Rj/8jYPzL2j4z9I2PkQWcAAAAAAABths4AAAAAAADaDJ0BAAAAAAC0GToDkMys\nWbNSl4DCkbF3tGP/yNg/MvaPjP0jY+RBZwCSmTFjRuoSUDgy9o527B8Z+0fG/pGxf2SMPCyEkLqG\nQcysQ1JXV1eXOjo6UpcDjBlr1qzRkUceKelmSSekLiehRyW9Uj/60Y/0rne9K3UxAAAAQGHK5bI6\nOzslqTOEUB6p9XJkAAAAAAAAbYbOAAAAAAAA2gydAUhm9erVqUtA4cjYO9qxf2TsHxn7R8b+kTHy\noDMAyVxyySWpS0DhyNg72rF/ZOwfGftHxv6RMfKgMwDJLFu2LHUJKBwZe0c79o+M/SNj/8jYPzJG\nHnQGIJnJkyenLgGFI2PvaMf+kbF/ZOwfGftHxsiDzgAAAAAAANpMy50BZna8ma0ws0fNrN/MSg3m\nPzGbr/LWZ2b75S8bAAAAAADkNTHHMrtLulPSEknXNLlMkHSEpGdemBDC4zkeG47MnTtXl156aeoy\nUKi5korL+Oqrr9avf/3rhvPts88+uuCCC2RmhdXSrmjH/pGxf2TsHxn7R8bIo+XOgBDCtZKulSRr\n7ZP1H0IIT7f6ePBr6tSpqUtA4YrK+GWaNOntuvLKZvoj+9XX97QOOuggnXrqqQXV075ox/6RsX9k\n7B8Z+0fGyMNCCPkXNuuXdHYIYcUQ85wo6UZJ6yTtKunXkhaEEG4dYpkOSV1dXV3q6OjIXR/gzZo1\na3TkkUdKulnSCanLGScek3SgVq5cqTPOOCN1MQAAAEBLyuWyOjs7JakzhFAeqfWOxgCCj0m6UNJ7\nJb1H0sOSbjKz40bhsQEAAAAAQJXCOwNCCGtCCN8IIfwqhPCLEMJHJN0q6eJGy55++ukqlUo73KZP\nn67ly5fvMN+qVatUKg0ex3DOnDlasmTJDtPK5bJKpZJ6enp2mD5//nwtXLhwh2nr169XqVRSd3f3\nDtMXLVqkuXPn7jCtt7dXpVJJq1ev3mH60qVLNWvWrEG1zZw5k+1gO3Jvh3S9pMHbIc2UtLxq2ipJ\ntcb5nKM49McOW5LN21M1fb6khVXT1mfzdldNX6Q4VkCl3mze6u1YqtHcDvYrtoPtYDvYDraD7WA7\n2A62YyxvR2dn56DvwGeeeeag5UdC4acJ1FnuEklvCyG8rc79nCbQBrq7u/Wa17wmdRnjyvg7TaBb\nUuqMOU2gSLRj/8jYPzL2j4z9I2PfxvNpArUcp/gJHW1s3rx5qUtA4cjYO9qxf2TsHxn7R8b+kTHy\naPlqAma2u6RXSxq4ksCrzOxYSU+EEB42sy9IOjCEcH42/yclPSjpXsUBBD8m6WRJ7xiB+jGOLV68\nOHUJKBwZe0c79o+M/SNj/8jYPzJGHi13Bkh6g+LVAUJ2+1I2/UpJsyUdIOmgivl3zuY5UPHE4bsl\nnRJCuCVnzXCCS6C0AzL2jnbsHxn7R8b+kbF/ZIw8Wu4MCCHcrCFOLwghzKr6+1JJl7ZeGgAAAAAA\nKEKqMQMAAAAAAEAidAYgmepLd8AjMvaOduwfGftHxv6RsX9kjDzoDEAyvb29qUtA4cjYO9qxf2Ts\nHxn7R8b+kTHysBBC6hoGMbMOSV1dXV3q6OhIXQ4wZqxZs0ZHHnmkpJslnZC6nHHiMUkHauXKlTrj\njDNSFwMAAAC0pFwuq7OzU5I6QwjlkVovRwYAAAAAANBm6AwAAAAAAKDN0BmAZHp6elKXgMKRsXe0\nY//I2D8y9o+M/SNj5EFnAJKZPXt26hJQODL2jnbsHxn7R8b+kbF/ZIw86AxAMgsWLEhdAgq3IHUB\nKBjt2D8y9o+M/SNj/8gYedAZgGS4UkQ7IGPvaMf+kbF/ZOwfGftHxsiDzgAAAAAAANoMnQEAAAAA\nALQZOgOQzJIlS1KXgMKRsXe0Y//I2D8y9o+M/SNj5EFnAJIpl8upS0DhyNg72rF/ZOwfGftHxv6R\nMfKwEELqGgYxsw5JXV1dXQyGAVRYs2aNjjzySEk3SzohdTnjxGOSDtTKlSt1xhlnpC4GAAAAaEm5\nXFZnZ6ckdYYQRqznhyMDAAAAAABoM3QGAAAAAADQZugMAAAAAACgzdAZgGRKpVLqElA4MvaOduwf\nGftHxv6RsX9kjDzoDEAyF110UeoSUDgy9o527B8Z+0fG/pGxf2SMPLiaADCOcDWBPDbJ7ACFsF1m\n1nDuPfd8me644xc67LDDRqE2AAAAYGhFXU1g4kitCADGppcphJ9I6lbjvs+n9eSTc3XffffRGQAA\nAADX6AwA0AZOUHNHUmyUNLfgWgAAAID0GDMAySxfvjx1CSgcGXtHO/aPjP0jY//I2D8yRh50BiCZ\npUuXpi4BhSNj72jH/pGxf2TsHxn7R8bIgwEEgcS2bNmi8877sH772982nPf557do7dpuSbdIOr7w\n2trPRkkHaMWKFTrrrLNSFwMAAAAwgCDg1bp167R8+fclnS3pwCaW+KCktxZbFAAAAADX6AwAxoxP\nS3p76iIAAAAAtAHGDAAAAAAAoM3QGYBkZs2alboEFI6MvaMd+0fG/pGxf2TsHxkjDzoDkMyMGTNS\nl4DCkbF3tGP/yNg/MvaPjP0jY+RBZwCSOe+881KXgMKRsXe0Y//I2D8y9o+M/SNj5EFnAAAAAAAA\nbYbOAAAAAAAA2gydAUhm9erVqUtA4cjYO9qxf2TsHxn7R8b+kTHyoDMAyVxyySWpS0DhyNg72rF/\nZOwfGftHxv6RMfKgMwDJLFu2LHUJKBwZe0c79o+M/SNj/8jYPzJGHnQGIJnJkyenLgGFI2PvaMf+\nkbF/ZOwfGftHxsiDzgAAAAAAANoMnQEAAAAAALQZOgOQzNy5c1OXgMKRsXe0Y//I2D8y9o+M/SNj\n5EFnAJKZOnVq6hJQODL2jnbsHxn7R8b+kbF/ZIw8LISQuoZBzKxDUldXV5c6OjpSlwMUqru7W9Om\nTZP0U0lvT11Om9so6QCtWLFCZ511VupiAAAAAJXLZXV2dkpSZwihPFLr5cgAAAAAAADaDJ0BAAAA\nAAC0GToDkEx3d3fqElA4MvaOduwfGftHxv6RsX9kjDzoDEAy8+bNS10CCkfG3tGO/SNj/8jYPzL2\nj4yRB50BSGbx4sWpS0DhyNg72rF/ZOwfGftHxv6RMfKgMwDJcAmUdkDG3tGO/SNj/8jYPzL2j4yR\nB50BAAAAAAC0GToDAAAAAABoMy13BpjZ8Wa2wsweNbN+Mys1scxJZtZlZlvMbI2ZnZ+vXHiycOHC\n1CWgcOMz41tuuUXLli1reFu1alXqUpOjHftHxv6RsX9k7B8ZI4+JOZbZXdKdkpZIuqbRzGZ2iKSV\nkr4q6QOSTpX0TTPbEEL4cY7HhxO9vb2pS0DhxlvGL9GkSUfoi1/8YtNL3HjjjTrppJOKK2mMox37\nR8b+kbF/ZOwfGSMPCyHkX9isX9LZIYQVQ8yzUNK7QgjHVExbKmlKCOH0Ost0SOrq6upSR0dH7vqA\n8aC7u1vTpk2T9FNJb09dDrQ1uzXyB0mv0vLly/Xud7+74JoAAADQrsrlsjo7OyWpM4RQHqn15jky\noFVvkXR91bTrJF02Co8NAC3aObs1srnoQgAAAIDCjMYAggdI2lg1baOkl5rZLqPw+AAAAAAAoMKY\nvprA6aefrlKptMNt+vTpWr58+Q7zrVq1SqXS4HEM58yZoyVLluwwrVwuq1QqqaenZ4fp8+fPHzTw\nxvr161UqldTd3b3D9EWLFmnu3Lk7TOvt7VWpVNLq1at3mL506VLNmjVrUG0zZ85s++3o6elxsR3S\nyOQhXS1pbtW0XkklSaurpi+VNHg7pJmSlldNW5WtY9CWKA79scOWZPP2VE2fr8GDAa7P5q3ejkV6\ncTsG1jPet2NA7e0Yy/tV0e2jurbxuh3V2I4Xt+Od73yni+3wkkcR29HT0+NiOyQfeRSxHQP/jvft\nGMB2DN6OynWP5+2o1K7b0dnZOeg78Jlnnjlo+ZEwGmMG3CypK4TwqYppF0i6LISwV51lGDOgDZRK\nJa1YUXfXaRu+xwwoSfKa8R8k7df2YwbQjv0jY//I2D8y9o+MfStqzIDRODLg55JOqZo2I5uONrZg\nwYLUJaBwC1IXgILRjv0jY//I2D8y9o+MkUfLnQFmtruZHWtmx2WTXpX9fVB2/xfM7MqKRf4tm2eh\nmR1pZn8m6Y8lfXnY1WNc46iPdkDG3tGO/SNj/8jYPzL2j4yRR54jA94g6VeSuiQFSV9SPEH3s9n9\nB0g6aGDmEMI6SWdIOlXSnZIulvSREEL1FQYAAAAAAMAoaPnSgiGEmzVEJ0IIYdDoCCGEWyR1tvpY\nAAAAAABg5I3pqwnAt+qRPeERGXtHO/aPjP0jY//I2D8yRh50BiCZcnnEBsLEmEXG3tGO/SNj/8jY\nPzL2j4yRx7AuLVgULi2IduL70oKecWlBAAAAFG88X1oQAAAAAACMIS0PIAigOZs3b9bzzz/fcL6n\nn356FKoBAAAAgBfRGQAUoKenR0cffawef3xDk0uYpClFlgQAAAAAL+A0ASRTKpVSl1CYxx9/POsI\n+CdJVzdx+6mk16UptlB+M0bkuR0jImP/yNg/MvaPjJEHRwYgmYsuuih1CaPgZElvTV1EQu2QcXtr\nj3bc3sjYPzL2j4z9I2PkwZEBSGbGjBmpS0DhyNg72rF/ZOwfGftHxv6RMfKgMwAAAAAAgDZDZwAA\nAAAAAG2GzgAks3z58tQloHBk7B3t2D8y9o+M/SNj/8gYedAZgGSWLl2augQUjoy9ox37R8b+kbF/\nZOwfGSMPOgOQzFVXXZW6BBSOjL2jHftHxv6RsX9k7B8ZIw86AwAAAAAAaDN0BgAAAAAA0GYmpi4A\nAMazxYsXa+XKlQ3n23vvvfW5z31Ou+yyyyhUBQAAAAyNzgAkM2vWLH3rW99KXQYKNUuS14z31oQJ\nH9dNN90p6e4G827T9u2/0vTp03X22WePRnGjhnbsHxn7R8b+kbF/ZIw86AxAMjNmzEhdAgrnOeOd\n1N//dfX3NzPvJkn7KIRQcE2jj3bsHxn7R8b+kbF/ZIw8GDMAyZx33nmpS0DhyNg72rF/ZOwfGftH\nxv6RMfKgMwAAAAAAgDZDZwAAAAAAAG2GzgAks3r16tQloHBk7B3t2D8y9o+M/SNj/8gYedAZgGQu\nueSS1CWgcGTsHe3YPzL2j4z9I2P/yBh50BmAZJYtW5a6BBSOjL2jHftHxv6RsX9k7B8ZIw86A5DM\n5MmTU5eAwpGxd7Rj/8jYPzL2j4z9I2PkQWcAAAAAAABths4AAAAAAADaDJ0BSGbu3LmpS0DhyNg7\n2rF/ZOwfGftHxv6RMfKgMwDJTJ06NXUJKBwZe0c79o+M/SNj/8jYPzJGHhZCSF3DIGbWIamrq6tL\nHR0dqcsBWnbffffp6KOPlvQzSW9NXQ6S2yRpH11zzTU655xzUhcDAACAcaRcLquzs1OSOkMI5ZFa\nL0cGAAAAAADQZugMAAAAAACgzdAZgGS6u7tTl4DCkbF3tGP/yNg/MvaPjP0jY+RBZwCSmTdvXuoS\nUDgy9o527B8Z+0fG/pGxf2SMPOgMQDKLFy9OXQIKR8be0Y79I2P/yNg/MvaPjJEHnQFIhkugtAMy\n9o527B8Z+0fG/pGxf2SMPOgMAAAAAACgzdAZAAAAAABAm6EzAMksXLgwdQkoHBl7Rzv2j4z9I2P/\nyNg/MkYedAYgmd7e3tQloHBk7B3t2D8y9o+M/SNj/8gYeVgIIXUNg5hZh6Surq4udXR0pC4HaNl9\n992no48+WtLPJL01dTlIbpOkfXTNNdfonHPOSV0MAAAAxpFyuazOzk5J6gwhlEdqvRwZAAAAAABA\nm6EzAAAAAACANjMxdQFoXz09Pdpnn31Sl4FC9Ugi4wE33nijnnnmmYbz7bXXXjrrrLNGoaLhox37\nR8b+kbF/ZOwfGSMPOgOQzOzZs7VixYrUZTQthKCvfe1r+t3vftdw3k2bNo1CRePBbEnjJ+Pi7KFJ\nk16rRYsWNb3EddddpxkzZhRY08gYb+0YrSNj/8jYPzL2j4yRB50BSGbBggWpS2jJ/fffrzlz5mji\nxINltlvD+SdOPFXbtx89CpWNZQtSFzBG7KJt2+6S1NfEvE9K2r+pIwjGgvHWjtE6MvaPjP0jY//I\nGHnQGYBkxtuVIvr7+yVJ27cvk/SWtMWMG+Mr42JNUHPDtEwqupARNd7aMVpHxv6RsX9k7B8ZIw8G\nEAQAAAAAoM3QGQAAAAAAQJuhMwDJLFmyJHUJKBwZe0c79o+M/SNj/8jYPzJGHrk6A8xsjpk9aGab\nzewXZvbGIeY90cz6q259ZrZf/rLhQblcTl0CCkfG3tGO/SNj/8jYPzL2j4yRh4UQWlvAbKakKyV9\nXNLtki6WdK6kI0IIPTXmP1HSDZKOkPTC8NghhMeHeIwOSV1dXV0MhoEx49e//rVe97rXSfq5GEAQ\nxXlS0t66+uqr9d73vjd1MQAAAEisXC6rs7NTkjpDCCPW85PnyICLJX09hPDtEEK3pE9I6lW8oPhQ\n/hBCeHzgluNxAQAAAADACGipM8DMJknqlPSTgWkhHlpwvaTpQy0q6U4z22Bmq8zsrXmKBQAAAAAA\nw9fqkQH7SNpJ0saq6RslHVBnmcckXSjpvZLeI+lhSTeZ2XEtPjYAAAAAABgBhV9NIISwJoTwjRDC\nr0IIvwghfETSrYqnGwzp9NNPV6lU2uE2ffp0LV++fIf5Vq1apVKpNGj5OXPmDBpZs1wuq1Qqqadn\nx+EN5s+fr4ULF+4wbf369SqVSuru7t5h+qJFizR37twdpvX29qpUKmn16tU7TF+6dKlmzZo1qLaZ\nM2e2/XaUSqVxuR3SpRo8Sn5ZUklS9bAZ8yUtrJq2Ppu3u2r6Iklzq6b1ZvOurpq+VNLg7ZBmSlpe\nNW1Vto5qc1T8dgw87njfjgGjux3joX1UtxGvr1ftvB0HHHCAi+3wkkcR21EqlVxsh+QjjyK2Y2C5\n8b4dA9iOwdtRWct43o5K7bodnZ2dg74Dn3nmmYOWHwktDSCYnSbQK+m9IYQVFdOvkDQlhHBOk+u5\nRNLbQghvq3M/Awi2gVWrVmnGjBmpy2gaAwjmsUrS+Ml4bBhfAwiOt3aM1pGxf2TsHxn7R8a+jYkB\nBEMI2yR1STplYJqZWfb3rS2s6jjF0wfQxnjBagdk7B3t2D8y9o+M/SNj/8gYeUzMscyXJV1hZl16\n8dKCkyVdIUlm9gVJB4YQzs/+/qSkByXdK2lXSR+TdLKkdwy3eADw6qmnntLGjdXDswy2yy67aM89\n9xyFigAAAOBJy50BIYTvmdk+kj4naX9Jd0o6LYTwh2yWAyQdVLHIzpK+JOlAxVMM7pZ0SgjhluEU\nDgA+7awJE/bQRz/60ebm3nlX3X77L3TssccWXBcAAAA8yTWAYAjhqyGEQ0IIu4UQpocQfllx36wQ\nwh9V/H1pCOHwEMLuIYR9Qwh0BECSBg2YAY/IuHW7q7//dkkrmrgt1datW/Tb3/42VbG04zZAxv6R\nsX9k7B8ZI4/CryYA1LN06dLUJaBwZJzPNElnNXF7V6oCX0A79o+M/SNj/8jYPzJGHnQGIJmrrroq\ndQkoHBl7Rzv2j4z9I2P/yNg/MkYedAYAAAAAANBm6AwAAAAAAKDN0BkAAAAAAECboTMAycyaNSt1\nCSgcGXtHO/aPjP0jY//I2D8yRh50BiCZGTNmpC4BhSNj72jH/pGxf2TsHxn7R8bIg84AJHPeeeel\nLgGFI2PvaMf+kbF/ZOwfGftHxsiDzgAAAAAAANoMnQEAAAAAALQZOgOQzOrVq1OXgMKRsXe0Y//I\n2D8y9o+M/SNj5EFnAJK55JJLUpeAwpFxsSbKbGede+65mjBhQsPblCl76/777x/RCmjH/pGxf2Ts\nHxn7R8bIY2LqAtC+li1blroESdKzzz6rzZs3N5zviSeeGIVqvBkbGfu1u0K4SdKvFUKjebfo6af/\nQvfcc4+mTZs2YhWMlXaM4pCxf2TsHxn7R8bIg84AJDN58uTUJejxxx/XUUcdo02bNja5xARJU4os\nyZn0Gfs3Pbs18oykvxjxRx8L7RjFImP/yNg/MvaPjJEHnQFoaxs3bsw6Ai6VdEQTS7xc0sj9qgoA\nAAAAKdAZAEiSTpD0ptRFAAAAAMCoYABBJDN37tzUJaBwZOwd7dg/MvaPjP0jY//IGHnQGYBkpk6d\nmroEFI6MvaMd+0fG/pGxf2TsHxkjDwuNh6AedWbWIamrq6tLHR0dqcuBY/fcc4+OOeYYSbeJ0wTg\n2zOSXqqrrrpK73vf+1IXAwAAgCaVy2V1dnZKUmcIoTxS6+XIAAAAAAAA2gydAQAAAAAAtBk6A5BM\nd3d36hJQODL2jnbsHxn7R8b+kbF/ZIw86AxAMvPmzUtdAgpHxt7Rjv0jY//I2D8y9o+MkQedAUhm\n8SuIdWwAAA7QSURBVOLFqUtA4cjYO9qxf2TsHxn7R8b+kTHyoDMAyXAJlHZAxt7Rjv0jY//I2D8y\n9o+MkcfE1AUAAEaDSZI+/el5+tznvtBw7ilTpuj73/+uDjzwwKILAwAAQAJ0BgBAW9hD0uV65JFf\n6pFHGs27TdI3dMstt+j9739/8aUBAABg1HGaAJJZuHBh6hJQODIeW2ZJ+tcmbpc1vUbasX9k7B8Z\n+0fG/pEx8qAzAMn09vamLgGFI2PvaMf+kbF/ZOwfGftHxsjDQgipaxjEzDokdXV1damjoyN1OXDs\nnnvu0THHHCPpNklvSl0OMEY8J2kPLV26lNMEAAAAEiuXy+rs7JSkzhBCeaTWy5EBAAAAAAC0GToD\nAAAAAABoM3QGIJmenp7UJaBwZDyebd26VVu2bBny9sgjj2jr1q2pS0WBeK32j4z9I2P/yBh50BmA\nZGbPnp26BBSOjMennTRhwm46//zztdtuuw15O+igg7TnnnvrrrvuSl00CsJrtX9k7B8Z+0fGyGNi\n6gLQvhYsWFDYuu+//3499thjDed74IEHCqsBkrQgdQHIZVf19/9M0r1NzPtbbd78Od1777069thj\niy4MCRT5Wo2xgYz9I2P/yBh50BmAZIq6UsRDDz2kY499vbZte76p+Xfa6aXq6zugkFrA1UDGr9dn\nt0Y2S/pcwbUgJa7q4x8Z+0fG/pEx8qAzAO489dRTWUfADyQd03D+vr6XStqn6LIA19auXavbb7+9\n4Xx77LGHjjrqqFGoCAAAAEOhMwCOvVLSq1IXATg3URMnvlLz58/X/PnzG85tZrrhhht00kknFV8a\nAAAA6mIAQSSzZMmS1CWgcGTs37e1ffuvJN3dxO0OhRD0yCOPpCoWOfBa7R8Z+0fG/pEx8qAzAMmU\ny+XUJaBwZOxfWfE0m9c1cXttohoxHLxW+0fG/pGxf2SMPCyEkLqGQcysQ1JXV1cXg2GgZXfddZeO\nO+44SXdIekPqcgC8YIuk3fSd73xHH/zgB1MXAwAAMC6Uy2V1dnZKUmcIYcR6fhgzAAAwqsrlsqZM\nmdJwvpe85CWMLQAAAFAQOgMAAKNkoiZNmqbLLrtMl112WVNLrFy5UmeccUbBdQEAALQfOgMAAKNk\norZt65L0dBPzbpU0VZs2bSq4JgAAgPbEAIJIplQqpS4BhSNj/1rNeDdJ+zd5w1jAa7V/ZOwfGftH\nxsiDIwOQzEUXXZS6BBSOjP0rNuMHHnhAd9xxR8P5dt99dx111FGF1tKueK32j4z9I2P/yBh5cDUB\nuMPVBAAP+jRx4qu0ffv6ppf48Y9/rFNPPbXAmgAAAEYfVxNA21u+fLkeeuihhvM9+uijo1ANgGLt\npO3buyRtaGLe7ZI6dfnll+vee+9tOPcee+yhWbNmacIEzpQDAADti84AjAv33XefzjnnHE2YsKvM\ndmo4/6RJR2nbtkNHoTIAxdknuzXSr0mTZuh731uh731vxZBzhhDU39+r3XbbTR/4wAdGpEoAAIDx\niJ9FkMzy5cubnvf555+XJPX3/0x9fc82vG3bdq+klxVUOZrXfMYYr8ZCxhO0bdt1Tb029Pc/I0na\nsmVL4prHj1ZeqzE+kbF/ZOwfGSMPjgxAMgsXLtTZZ5+dugwUaqEkMvZtfGa8fft2bd26teF83//+\n93Xbbbc1tc7dd99dn/nMZzR58uThljem8FrtHxn7R8b+kTHyyNUZYGZzJP2VpAMk3SXpz0MIdYd7\nNrOTJH1J0tGS1kv6PyGEK/M8NvzYd999U5eAwpGxf+Mv4wkT9tCFF16oCy+8sOllJk16XcN5tm27\nR69+9av1kY98ZDjljTm8VvtHxv6RsX9kjDxa7gwws5mKX+w/Lul2SRdLus7Mjggh9NSY/xBJKyV9\nVdIHJJ0q6ZtmtiGE8OP8pWO827Ztm4466jg9/HDjQQH7+rZn/9ul2KIAODdB/f2rJd3d5Py7SzpH\n27ZZk+vuz18aAADAKMpzZMDFkr4eQvi2JJnZJySdIWm2pEtqzP+nkh4IIczL/v6Nmb09Ww+dAW3s\nueee0/333yXpk5Je0cQSByseXAIAw3Fsdht5119/vXp7exvON3nyZM2ePVs77dR4QFQAAIAitNQZ\nYGaTJHVK+vzAtBBCMLPrJU2vs9hbJF1fNe06SZe18tgYP9atW6eenkEHiQzy7LPPZv+7QNJxRZYE\nAIWbNOmduvrqlbr66pUN5+3v79XHP/5x7b//K5uYd7u++P+3d/chcl11GMe/T2J3N4k2CcRutloI\n9aWamhetotU0SdmamhRq1QVfQgxaqrUlDRVBi4jQPywRW7VqCCI1LVWhEMRYkJZ2UyQkaWvajaDb\nDcFd+pKXttuQSEzjmv78496VdZzZvXecnd2z83xgIXv3d2fO8Mszw5w598wPtrF06dIJa+fMmcPl\nlxebNB0YGKCvr2/CuqNHj7Jr1y7Wrl2LNPEKidmzZzN//vxCYxgeHmZwcLBQbXt7O8uWTXy5hpmZ\nmRVTdmXAImA2cKLi+AngshrnLK5Rf6Gk9og4V+WcDoD+/v6Sw7OiIoLt27fz9NM1t3r4j7Nnz3Lk\nyJFJHM3DZFec2MyzF/j5VA/CJpV7PGpk5AaKb6Z4ADjEicpXx6qeYfPmzXWPqxF6enqm9P5HtbW1\n0d3dPWX339/fz6lTpwrXL1y4kHnz5hWqXb58eaHawcFBjh8/XngMixYtoqOjY9yaPXv2sGHDhsK3\nOTAwwMmTJwvXr1q1iq6urgnrhoaGOHbsWOHbLfLYRnV2dtLW1lb4tmea3t7eUvukWHp6e3vZtGkT\nw8PDherb29vZsmULCxYsmOSRWSOMeV9c7EmvIEVE8WKpC3gJuDIinhxzfBuwOiL+Z3WApAHgvojY\nNubYerJ3gHOrTQZI+gLwqzIPxMzMzMzMzGwG2xgRv27UjZVdGfAqcB7orDjeCdSapj5eo/50jVUB\nkF1GsBEYAvxl0GZmZmZmZtaqOoAlZO+TG6bUZEBEjEg6CHQDuwGUXUDYDdxb47T9wPqKY+vy47Xu\nZxho2IyHmZmZmZmZWcL2NfoGZ9Vxzj3ATZK+KOk9wA5gLrATQNJdku4fU78DuFTSNkmXSboF6Mlv\nx8zMzMzMzMyarPRXC0bEQ5IWAXeSLffvA66NiFfyksXAJWPqhyRdR/btAbcBLwI3RkTlNwyYmZmZ\nmZmZWROU2kDQzMzMzMzMzNJXz2UCZmZmZmZmZpYwTwaYmZmZmZmZtZgpmQyQdKukQUlnJR2Q9KEJ\n6tdKOijpdUmHJW1u1litfmX6LGmNpDcqfs5LuqiZY7ZiJF0labekl/JeXV/gHOc4IWV77AynR9Id\nkp6SdFrSCUm/lfTuAuc5y4mop8fOclok3SzpkKRT+c8+SZ+Y4BxnOCFle+wMp0/St/K+jbvhfiOy\n3PTJAEmfBe4Gvgu8HzgEPJJvSlitfgnwMPA4sAL4MfALSR9vxnitPmX7nAvgXWSbUC4GuiLi5cke\nq9VlHtnmobeQ9W1cznGSSvU45wyn5SrgJ8CHgWuAC4BHJc2pdYKznJzSPc45y+l4Afgm8AHgCqAX\n+J2k91YrdoaTVKrHOWc4UfmHp18he+80Xt0SGpDlpm8gKOkA8GREbM1/F9l/8nsj4vtV6rcB6yNi\n+ZhjvwHmR8SGJg3bSqqjz2vIntwWRsTppg7W/i+S3gBuiIjd49Q4xwkr2GNnOHH5ZO3LwOqI2Fuj\nxllOWMEeO8uJkzQMfCMiflnlb87wDDBBj53hREl6M3AQ+BrwHeDZiPh6jdqGZLmpKwMkXUA2o/X4\n6LHIZiMeA66scdpH8r+P9cg49TbF6uwzgIA+SUclPSrpo5M7Umsi57g1OMNpW0D2adJr49Q4y2kr\n0mNwlpMkaZakzwFzgf01ypzhhBXsMTjDqfoZ8PuI6C1Q25AsN/sygUXAbOBExfETZEtYqllco/5C\nSe2NHZ41SD19PgZ8FfgM8GmyVQRPSFo5WYO0pnKOZz5nOGH56q0fAXsj4q/jlDrLiSrRY2c5MZLe\nJ+nvwDlgO/CpiHiuRrkznKCSPXaGE5RP8qwE7ih4SkOy/KaihWaTKSIOA4fHHDog6R3A7YA3tjGb\n5pzh5G0HlgIfm+qB2KQp1GNnOUnPkV0zPB/oAR6QtHqcN4uWnsI9dobTI+ntZJO110TESDPvu9kr\nA14FzgOdFcc7geM1zjleo/50RJxr7PCsQerpczVPAe9s1KBsSjnHrckZToCknwIbgLURcWyCcmc5\nQSV7XI2zPI1FxL8i4m8R8WxEfJts47GtNcqd4QSV7HE1zvD0dgXwVuAZSSOSRoA1wFZJ/8xXdlVq\nSJabOhmQz3QcBLpHj+UPrhvYV+O0/WPrc+sY/zoZm0J19rmalWRLnSx9znFrcoanufxN4ieBqyPi\n+QKnOMuJqaPH1TjLaZkF1Fom7AzPDOP1uBpneHp7DFhG1qcV+c+fgAeBFVF9x/+GZHkqLhO4B9gp\n6SDZLNXtZJtg7ASQdBdwcUSMLmPZAdya75h4H9mD7iGb4bbpq1SfJW0FBoG/AB3ATcDVgL/qZhqS\nNI9shnl0pvJSSSuA1yLiBec4fWV77AynR9J24PPA9cAZSaOfMJyKiNfzmu8Bb3OW01RPj53ltOT9\n+wPwPPAWYCPZJ4rr8r/79ThxZXvsDKcnIs4A/7WXi6QzwHBE9Oe/T8rrcdMnAyLiofyrbe4kW8rQ\nB1wbEa/kJYuBS8bUD0m6DvghcBvwInBjRFTunmjTSNk+A23A3cDFwD+APwPdEfHH5o3aSvggsIds\nV+og6x3A/cCXcY5nglI9xhlO0c1kvX2i4viXgAfyf3fhLKesdI9xllNzEdnzchdwiqxf68bsRu7X\n4/SV6jHO8ExRuRpgUl6PVX3VgZmZmZmZmZnNVM3eQNDMzMzMzMzMppgnA8zMzMzMzMxajCcDzMzM\nzMzMzFqMJwPMzMzMzMzMWownA8zMzMzMzMxajCcDzMzMzMzMzFqMJwPMzMzMzMzMWownA8zMzMzM\nzMxajCcDzMzMzMzMzFqMJwPMzMzMzMzMWownA8zMzMzMzMxazL8BtyaZoLJ7xmQAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1162c04d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "overal_error_dist = pd.Series(class_samples)\n",
    "\n",
    "print \"Mean monthly error: \", np.round(overal_error_dist.mean(), 3)\n",
    "print \"Std. dev. of monthly error: \", np.round(overal_error_dist.std(), 3)\n",
    "print \"5/95 percentiles: \", np.percentile(overal_error_dist, [5, 95])\n",
    "\n",
    "\n",
    "overal_error_dist.hist(bins=70, normed=True, histtype=\"stepfilled\")\n",
    "plt.title(\"Overall relative error distribution for outdoor target$\", fontsize=16)\n",
    "figsize(12.5, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Under the previous assumptions, we can therefore expect 90% of the water retailers to have their true outdoor targets within approximately $\\pm$40% of the estimated target."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "Olofsson, P., Foody, G. M., Stehman, S. V., & Woodcock, C. E. (2013). Making better use of accuracy data in land change studies: Estimating accuracy and area and quantifying uncertainty using stratified estimation. Remote Sensing of Environment, 129, 122-131."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}