fig4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os, glob\n",
    "import datetime\n",
    "from netCDF4 import Dataset\n",
    "from collections import namedtuple\n",
    "from scipy import stats\n",
    "\n",
    "from matplotlib.pyplot import *\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def colorplotLeft(ax):\n",
    "    \n",
    "    almost_black = '#262626' # tweek from Thorsten\n",
    "    \n",
    "    for spx in ['top','right']:\n",
    "            ax.spines[spx].set_visible(False)\n",
    "\n",
    "    # Only show ticks on the left and bottom spines\n",
    "    ax.yaxis.set_ticks_position('left')\n",
    "    ax.xaxis.set_ticks_position('bottom')\n",
    "    \n",
    "    # set axis positions\n",
    "#     ax.yaxis.set_label_position(pspines[0]) \n",
    "#     ax.xaxis.set_label_position(pspines[1]) \n",
    "    \n",
    "    for spx in ['bottom','left']:\n",
    "        ax.spines[spx].set_linewidth(0.5)\n",
    "        ax.spines[spx].set_color(almost_black)\n",
    "    \n",
    "    for ticks in ax.xaxis.get_ticklines() + ax.yaxis.get_ticklines():\n",
    "        ticks.set_color(almost_black)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def anomaly(data):\n",
    "\t\"\"\"calculate anomaly of a 1D array (assumes no missing time slots)\"\"\"\n",
    "\td1=np.ma.copy(data)\n",
    "\tind=np.ma.array(range(len(data))) % 12\n",
    "\tfor num in range(12):\n",
    "\t\td1[num == ind]=data[num == ind]-np.ma.average(data[num == ind])\n",
    "\treturn d1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "## ceres data; from fig. 3 file\n",
    "## based on regression up to 7/17\n",
    "ceresR=-0.17\n",
    "ceresRta=-0.46"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## read in models & calculate correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "## this cell reads CMIP5 model files that have been output by processCMIP.py\n",
    "\n",
    "modelRun = namedtuple('modelRun', ('t', 'ts','t_trop','ts_trop','flux'))\n",
    "modelType = namedtuple('modelType', ('piControl','abrupt4xCO2'))\n",
    "\n",
    "def readFile(modx,runx):\n",
    "\t\"\"\"read global and tropical avg. nc file\"\"\"\n",
    "\tfni=Dataset('cmip5/{0}.{1}.nc'.format(mn,runx))\n",
    "\tovar=[]\n",
    "\tfor ii in ['flux','ta','tas','t_trop','tas_trop']:\n",
    "\t\tovar.append(anomaly(fni[ii][:].squeeze()))\n",
    "\tmodout=modelRun(flux=ovar[0],t=ovar[1],ts=ovar[2],t_trop=ovar[3],ts_trop=ovar[4])\n",
    "\treturn modout\n",
    "\n",
    "# get names of models\n",
    "fn=glob.glob('cmip5/*.{}.nc'.format('piControl'))\n",
    "models=[ii.split('/')[1] for ii in fn]\n",
    "models=[ii.split('.')[0] for ii in models]\n",
    "models.sort()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ACCESS1-0 ACCESS1-3 BCC-CSM1-1 BCC-CSM1-1-M CCSM4 CNRM-CM5 CNRM-CM5-2 CSIRO-Mk3-6-0 "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/anaconda/lib/python3.6/site-packages/scipy/stats/_stats_mstats_common.py:107: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  slope = r_num / ssxm\n",
      "/usr/local/anaconda/lib/python3.6/site-packages/scipy/stats/_stats_mstats_common.py:119: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  sterrest = np.sqrt((1 - r**2) * ssym / ssxm / df)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CanESM2 FGOALS-g2 FGOALS-s2 GFDL-CM3 GFDL-ESM2G GFDL-ESM2M GISS-E2-H GISS-E2-R HadGEM2-ES INMCM4 IPSL-CM5A-LR IPSL-CM5A-MR IPSL-CM5B-LR MIROC-ESM MIROC5 MPI-ESM-LR MPI-ESM-MR MPI-ESM-P MRI-CGCM3 NorESM1-M "
     ]
    }
   ],
   "source": [
    "skip=17*12;length=17\n",
    "\n",
    "feedbacks={};modelOutput={}\n",
    "for mn in models:\n",
    "\tpiC=readFile(mn,'piControl')\t\n",
    "\tabrupt=readFile(mn,'abrupt4xCO2')\n",
    "\tfeedbacks[mn]=modelType(piControl=piC,abrupt4xCO2=abrupt)\n",
    "\tprint('{} '.format(mn),end=\"\")\n",
    "\tsys.stdout.flush()\n",
    "\n",
    "\t## analyze the run\n",
    "\tm1all=[];abrupt=[]\n",
    "\tfor latr,varn in zip(['90','30'],['ts','t']):\n",
    "\t\t\n",
    "\t\tmod1=feedbacks[mn].piControl # control run\n",
    "\t\ts=len(mod1.t)\n",
    "\n",
    "\t\ttvar=varn\n",
    "\t\tif latr == '30': tvar = tvar+'_trop'\n",
    "\t\tcommandStr='tg=mod1.{}'.format(tvar)\n",
    "\t\texec(commandStr)\n",
    "\n",
    "\t\t# flux is the same for all calcs\n",
    "\t\tflux=mod1.flux\n",
    "\t\t\n",
    "\t\tm1=[]\n",
    "\t\tfor ii in range(0,s-length*12,skip):\n",
    "\t\t\tindx=range(ii,ii+length*12)\n",
    "\t\t\tfitx=stats.linregress(anomaly(tg[indx]),anomaly(flux[indx]))\n",
    "\t\t\tm1.append(fitx[2])\n",
    "\n",
    "\t\tm1all.append(m1)\n",
    "\n",
    "\tmodelOutput[mn]=modelType(piControl=m1all,abrupt4xCO2=[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## figure 4 plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcQAAAIfCAYAAADuTqxuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XucVVX9//HXGzJgRMuvIF6QGbxNFiIp+VXLxMs3L9HFvCEjSfJt/GahfVO/XlDzl41YWaJjWmMmFuOFvGUKmmJkmuR1VBCPWTKImobhhQYR4fP7Y62Dh8M5M2dmzpzr5/l4nMc5s/fae68zW1mz1t57vWVmOOecc9WuX7Er4JxzzpUCbxCdc845vEF0zjnnAG8QnXPOOcAbROeccw7wBtE555wDvEF0zjnnAG8QnXPOOcAbROeccw7wBtE555wD4EPFrkA1uPXWW+0rX/lKsatROcaNC+/z5xezFi7NuHhe5vt5caVHuRTyHmIBvPzyy8WuQnlobYW6OujXL7y3tha7Rs65KuI9RFcaWluhsRE6OsLP7e3hZ4CGhuLVyzlXNbxBrDTJ4cRys2ABrF694bKODpgyBa6+esPlbW3hvS++qw/3OVe1Cj5kKukISSbpYynLdpE0R9ILkhZLmi1pWFy3l6QHJCUkPSfpF5JqJE2W9E9JbSmvj0vqJ+lySQslPSPpUUkj475OjMuejuu/FJcfLWmRpHWSxnZS9/+QdK+kv8b3Lfr691U10hvDzpbvtlt4uZIyZ84c5syZU+xqONdjxeghHgc8CEwALpA0ELgL+I6Z/Q5A0gHAUEkAvwEmmNnDCguOBDaL+7rJzL6VunNJxwHbAqPNbJ2k4cC/4/s0YA8ze0vSYGBo3Gwh8BXg513U/SxgnpldLOms+POZPf5N9IVy7eHU1YVh0nS1teX7napMTU1NsavgXK8UtEGMjdCngQOAO4ALgInAw8nGEMDM/hDLfw+4zswejssNuDmuy3aYbYBXzWxd3GZZLD8SeAdYGZevTPm8uIt9Jn0JGBc/XwfMp9QaxHLV1LThNUSAmpqwPN2VV4b3k08uTN1cTq6M5+VkPy8bWLNmDcuWLePdd98tdlUqVv/+/fnoRz/KkCFD6Nev5wOfhe4hfhm428yel/QvSXsAo4DHs5QfRWh4sjlW0mdSft4HmA08KGk/YB4wy8yeBJ4CXgNelDQPuDW1Ec7RMDN7FcDMXpW0VTe3d9kkb5yZNg2WLoURI0JjmOmGmtmzw7v/w1tSZsfz4g3ihpYtW8Zmm21GXV1dLn90u24yM9asWcNrr73GsmXLGDFiRI/3VegG8ThgRvx8Y/y5NzYaMgWWSaoHDoyveZKONrN5kg4FPgUcBFwqaU8zu6CXdciovr6+EWiED57Pcl1oaPA7Sl3Feffdd70x7EOS+PCHP8x2221HIpHo1b4K1iBK2pLQQI2SZEB/wID/B+yfZbNFwJ7Ab7tzLDNbDcwF5kp6jdAznReHXB8BHpF0L3AtYdg2W52vBT4JvGJmhwOvSdom9g63AV7Ptm0ikWgBWgCam5utO/V3zlUWbwz7Xm+GStfvIw/1yNVRwK/MrNbM6sxse+BF4HlgX0mfTxaUdKik3YArgBMk/WfKuuMlbZ3tIJL2kLRt/NwPGA20S9o2DtEmjQEy3MXxATP7mpmNiY0hhOueJ8TPJ9DNhto5B62trdTV1dGvXz/q6upo9QkYXIko5JDpccDFactuIdxUMx6YIWkGsAZ4GjjVzF6TNAG4JF6vWwc8ANwat0+/hngysDlwtaQBcdkjhIZ1WNzPtsC7wD+B/4HwKAjQTLjr9C5JbWZ2SIbvcDEwW9IUYClwdM9+Fc5Vp9bWVhobG+mIN0+1t7fTGCdgaPDhcldkCqOIri81Nzfb1KlTi10NV0Kq9bryggULWJ3h2dIBAwaw9957F6FGfe+73/0u22677UbL6+vri1Cbja1YsYL/+I//YNNNN2XdunUMHTqU888/nylTphTk+Ndddx3f/OY3AXj//fdZs2YNgwYNAuBDH/oQr732GgMGDOhsF+stXryYXXfdNdMqn8vUOVdaMjWGnS13fa+trY0hQ4awcuVKOjo6mD59OieddBLLly8vyPFPOOEEVq5cycqVKznnnHMYP378+p/ffPPNnBvDfPCp21z5ueSS8H766cWtRy9UYiLEJfG8nN7Jeamrq6M9wwQMtbW1Ffk7gdBr6XFvsLU1t0eReqGtrY099vjg9or999+ftWvXsmLFCoYMGbJB2YsvvpjHHnuMm2++ef2yU089FTPj8ssv5wc/+AGXX345b7/9Nttuuy1XXnklBx10ULfqsvvuu/f+S/WQ9xBd+bnzzvByJeXOO+/kzi7OS1NT00Yz2tTU1NCUaQKGapec8L69Hcw+mPA+zzchPfnkk+y5554AvPnmm5x99tnsueee7LTTThuVPe6445gzZw5vv/02AGvXrmX27NlMnDiRRCLBFVdcwaOPPso777zDPffcQ11dXbfq0tbWxpgxY3r9nXrKG0TnXME0NDTQ0tJCbW0tkqitraWlpcVvqEkaN+6D15QpG87cBB9MeJ8skwdtbW1cdtllbL755myxxRa8/vrr3H333RkfFamtrWWPPfbg9ttvB+D++++npqaGvffem/79+7N69WqeffZZ1qxZQ11dHTvuuGPO9Xj77bdZsmRJxgZx0aJF9O/fn2XLlvX8i+bAG0TnXEE1NDSwZMkS1q1bx5IlS7wxzKY7E973+BCrWbx4MU8//TRvv/02N998MwsWLGCTTTbJus3EiRO54YYbALj++uuZOHEiADvttBMzZszgggsuYKuttmLChAm88sorOdflqaeeYrPNNmPkyJEbrfvBD37ApEmTWLx4cTe/Yfd4g+hcNh5Y7Apt/vwPXrW1mcskJ7zPwzXXhQsXMmDAAHbYYQcAjjzySEaMGMEtt9ySdZujjz6a+fPns2zZMm677bb1DSKExvLBBx+kvb0dSZx5Zu5TPbe1tTF69OiNeqZPP/0022yzDYcccog3iM5tZNCg8OpLBbp+U0kGDRq0/nZ5lwdNTWGC+1TZJrzvoSeffJJRo0Zt0Agdfvjh3HHHHVm3GTp0KOPGjeNrX/saI0eOXP+YQyKR4P7772f16tUMHDiQQYMG0b9//5zrku364aWXXsqZZ57Jxz/+8T5vEP0uU1e6urpG0pfP8nUnsLgUFeGOzblz5xb8mBWtOxPe91CyV5bq0EMPpbm5mXfffZeBAwdm3G7ixIl89atf5Yc//OH6ZatXr+ass85i8eLFbLLJJuy77760tLQAcNhhh7HffvtxzjnnZK3LU089xTe+8Y2N6vfQQw8xceJE1q5dy9q1a3v6VXPiD+YXgD+Y30PFfHj9j3/Mvm7/bFPvlpAKfYShHHXysLjrwoQJE7jqqqvYYouQxb7XXnvxyCOPZC3f2wfzC95DjNOk3QrsambPxWW7EFIwdiFM3fYMMDVO3bYXcAlh6jUjhAufAhwD/Ah4OWX3E4Hn4r4OjOXfBY4xsxclnQj8b1zeD5hmZr+VdDRhku9dgb3M7LEsdb+QkIm4jjCx92Qzy/2qseuebP+oX3hheD/vvL47tgcWd9uF8byc15fnxVWNxx9/nEGDBq1vDAEGDhzIG2+8wZZbbtknxyzGNcTjCI3aBABJA4G7gKvMbCcz2xW4ChgqaRjwG+BMM6snNFh3A5vFfd0UJ99Ovp4FjgW2BUab2W7AEcCbkoYD04DPmNloYG/CnKkAC4GvEOZJ7cyPzGy0mY0B7gTO7/Vvw3XfvHnh1ZcKcP2m0sybN495fX1eXNXYc889ufbaazdY9sADD/RZYwgFbhAlDQY+DUwhNoiEXt3DqWG9ZvYHM1sIfBO4zswejsvNzG42s9c6Ocw2wKtmti5us8zMVgBbAe8AK+PylWb2Yvy82My6DNIys7dTftyU0NN0laihAVpaQo9QCu8tLZ7X6FwFK/SQ6ZeBu83seUn/inFMo4DHs5QfBVzXyf7S0y72AWYDD0raD5gHzDKzJ4GngNeAFyXNA25NbYRzJakJ+CrwFnBAd7d3ZcQDi52rKoUeMj0OuDF+vjH+3BvpQ6arzGwZUA+cTbjWN0/SQWa2FjiUkMv4PHCppAu6e0AzmxazHFuBb2UrV19f31hfX/9YfX39YwsXLuzBV3POOVdIBWsQJW1JuNHlF5KWAGcQrvctAvbMslln67Iys9VmNtfMzgAuIvRMk0Ouj5jZdMKQ7ZFd1PlaSW2S5mRYfX1n2ycSiZZEIjE2kUiMHTVqVHe/guvMlluGlyspW265ZU7Xdzwg2JWqQg6ZHgX8ysxOSi6Q9EdCb+1sSZ83s7vi8kMJd49eATwi6S4z+0tcdzxwX7aDxGHYf5jZK5L6AaOBp2Mw8NZm9kQsOgbIcBvhB8zsa2n73tnM/hp//CLhjlZXaJ3MouGKp7PZTZI8INiVskI2iMcREudT3UK4qWY8MEPSDMJjF08Dp8bHLiYQku63IgyBPkB4bAM2voZ4MrA5cLWkZIjWI4SGdVjcz7aERzH+CfwPrH8UpBkYCtwlqc3MDsnwHS6WVB/r0Z7c3rm+UmlBwpkCgjs6OpgyZQpXl8OEBz3w3e9+l379uh6MK5XA4GpWsAbRzMZlWHZ5yo+HZtnuYWC/DKtmxlcmd2dY1k4Yss10jNuA27LsK7Vcp0OsrkDOPju8T59e3Hq4Dfz9738HWD8vZiYeEOxKmU/d5srPww8XuwYFU06hucnebGd19oBgV8p8cm/nXMF4QLArZd4gOucKxgOCS8+KFSuQxODBg6mpqaG2tpZrrrmmYMe/7rrrGDx4MIMHD2bgwIH0799//c8f/ehHCzqc7g2ic66gPCC4tLS1tTFkyBBWrlxJR0cH06dP56STTmL58uUFOf4JJ5zAypUrWblyJeeccw7jx49f//Obb77JgAEDut5JnniD6MrP8OHh5UrK8OHDGe7nJa8K8cxmW1sbe+yxx/qf999/f9auXcuKFSs2KnvxxRdz1FFHbbDs1FNP5ZRTTgFCsv12223HZpttRn19fbfntm1ra2P33XfvwbfIEzPzVx+/Lr/8cnNVYNYss9paMym8z5pV7Bq5EvDss8/2aLtZs2ZZTU2NEeZMNsBqampsVp7/u5o0aZKdffbZZma2YsUKmzRpku255562bt26jcouWbLEBg0aZG+99ZaZmb3//vu29dZb28MPP2zPPfecDR8+3F5++WUzM3vxxRfthRde6FZdRo4cabfcckuPv0snv+uc/q32HqJz+dDaCo2NITLKLLw3NoblzuVo3Lhx619TpkxZP4FBUvKZzWSZfGhra+Oyyy5j8803Z4sttuD111/n7rvvRto4QrC2tpY99tiD22+/HYD777+fmpoa9t57b/r378/q1at59tlnWbNmDXV1dey444451+Ptt99myZIljBkzZoPlc+bMYdy4cey3337su+++fRoS7I9duPKTHJbbaafi1iPVggWQfvG/owOmTIEKfeA83bdfeAGAGYU4LxX6iEaqQjyzuXr1ahYvXsxzzz3HjjvuyC233MKUKVPYZJNNsm4zceJEbrjhBr761a9y/fXXM3HiRAB22mknZsyYwQUXXMCiRYs45JBD+MlPfsK2226bU12eeuopNttsM0aOHLl+2TvvvMP555/PQw89xIABA1ixYgX9+/fv3ZfuRJ/0ECWtjXOAPiXpCUn7pqzbRdIcSS9IWixpdsw9RNJekh6QlJD0nKRfSKrJsP9NJF0s6a+SFkp6RNJhcd2Jkp6R9HRc96W4fKakDkmbpeznMkkmaUj8+ZeSXpfU6WzcuZZzfWTlyvAqJdn+kaqiB87bVq6krdTOS5mZP3/++ldtbW3GMslnNvPx3ObChQsZMGDA+skUjjzySEaMGNHpNHxHH3008+fPZ9myZdx2223rG0QIjeWDDz5Ie3s7kjjzzDNzrktbWxujR4/eoGf6j3/8g5UrV/LSSy8BbBAW3Bf6qoe4ykKILpIOAaYD+6eEAX/HYvSSpAMIYcAQwoAnmNnDCguOJIQBd6Tt/0JC7uEoM1sdG9T9U0KA9zCzt2L+4tCU7V4gJN7PivOcHkCYMzVpJmGat1918f1yLef6QnJIpZR6CXV1YZg0XW1tadWzLyWH8Krl+/axpqamDeZ9hfw/s/nkk08yatSoDRqhww8/nDvuuIMTTzwx4zZDhw5l3LhxfO1rX2PkyJHsuuuuACQSCV5++WU+/elPM3DgQAYNGsS6detyrktbW9tGw6U777wzp5xyCkcddRT9+/dn5syZ7Lbbbj34prkpxDXEzYHk7Uq9DgOOPcavA1PNbHUs+5qZzaaTEODoBkLCBsA44CHg/ZS6PAD8q6svlGs5V0WamiDtgXNqasJy53qgEM9sJntlqQ499FDuvfde3n333azbTZw4kfvuu2+D3uHq1as566yzGDJkCFtvvTWvv/46F110EQCHHXbY+s/ZPPXUUxs1iAAnn3wybW1tnHjiiVx22WXd+Xrd1lc9xEGS2oCBhJ5ccg7R3oQBJ+0ELLUN0+uTugoB/ivwJUlbECYbnwUclsMxnetc8h+padNg6VIYMSI0hv6MneuFhoaGPn1O84orrtho2bhx4/j3v//d6XaTJk1i0qRJGywbPXo0jzzySMbyc+fO7bIujz322EbLnn/+eXbZZRcA3nzzzT6fAq8QQ6b7AL+S1OehgGa2NkZHfQo4iBACvKeZXZBS7FZCFuJ/AidtvJf8qK+vbwQaofISC4ou/g9SchoaqroB3KVUz4srW+effz4vvvgigwYNor6+nubm5j49Xp/fZRqvBw4hXMtbBOyfpWgyDPi36Ssk3UOIb3oMOAUYIWkzM3snw/GMEPn0iKR7gWuBC1KK3Ag8QRieXZfp1uK0Y28PJHuZPzOzn3W6QZRIJFqAFoDm5mbLZRuXo5aWYtfAZdDi58Xl2Y033ljQ4/X5NURJHwP6A28QUub3lfT5lPWHStqNcJPKCZL+M2Xd8ZK2NrNDzGyMmf23mXUA1wCXS/pwLLdNLLttDAhO2igE2MyWEm68uTKX+pvZS/HYY3JtDJ2rVoWYWcW5vtLX1xABBJxgZmuBVZJ6Ggac6lzg+8Czkt4F/g2cD2xClhDgVGb280yVlnQD4WabIZKWAd81s41muc21nOsjMWHde4ql5cADD+RPf/oT778f7lNrb2+nMZ4rn6/UlYM+aRDNLOuTk2b2HN0PA04v9x7wf/GVLlsI8OQsy+tSPh/X1bG7U871keefL3YN+kS5X2t+4IEHCFcsPpCcWeXqKpmcIKlSsx0rnU/d5pzLi/TGMKmQ8T3O9YZP3eZciSj3XsXAgQMzNn7JmVWqmZllnBvU5U93JgHIxnuIzrm8GDlyJP36bfhPSr5nVilHAwcO5I033sjag3a9Y2a89957vPzyy2y66aa92pf3EF35yTCbhSu+Qw45hNraWp577jmWLl3KiBEjaGpqqvobaoYPH86yZcv45z//WeyqVKwPfehDfOQjH2HIkCG92o/8r5a+19zcbFOnTi12NZxzrlrlNF7tQ6auOrW2hgm5+/UL7/68nHNVz4dMXfk5/vjwPmtWz7ZPhvkmUwSSYb5Q1VOv9dbx8bzM6ul5ca7IvEF0vVfo5+fa2np33HIM8y2DuzSXLVtW7Co41ys+ZOqqj4f5Oucy6LKHKGkt8AzhouRa4Ftm9ue4bhdgBrALYRq2Zwg5ha9J2gu4hDAptwEPAqfEuUhT978JIfD3SGA1IQz4u2Y2V9KJwP/G7fsB08zst5JmAscAw5ITfEu6jDDx91AzWy7pl8B44HUzy5q00Y1yBswys0nx5w8BrwJ/MbPxXf0eK1qhey+9DaL1MF/nXAa59BBXxYmtdwfOBqYDSBoI3AVcZWY7mdmuwFXA0Jhg/xvgTDOrB3YF7gY2y7D/CwmZiaNig/QFYDNJwwmTcH/GzEYDexPmPU16AfhSrEs/4ADg5ZT1M8kyRVyaXMv9GxglaVD8+b/SjufKhYf5Oucy6O6Q6ebAivh5IvBwagCvmf3BzBYC3yTEKz0cl5uZ3Wxmr6XuTFIN8HVCr3J1LPuamc0GtgLeAVbG5SvN7MWUzW8Ajo2fxwEPAe+n1CWnVPtcy0VzgWRSx3GxDq7Q9tknvHqqoSFMDF5bC1J4b2nxG2p6aZ999mGf3pwX54osl5tqkskVAwk9ueTk2aOAx7NsMwq4Lod97wQsNbO3M6x7CngNeFHSPODW1MYX+CvwJUlbEBqnWcBhORyzN24Ezpd0JzAa+CU5TEbu8mz69N7vo8rDfPvC9HycF+eKKJcGcZWZjQGQtA/wK0lZr7Xli5mtlXQo8CngIOBSSXua2QUpxW4FJgD/CZxUgDo9LamO0ADP6axsfX19I9AI5Z9i4Jxz1aBbQ6ZxCHQIMJQPEu4zybpO0j2S2iT9gnAdcISkTNcWk0Otj5jZdELDd2RakRsJ1yDvNbMuZ3aVtH08dpukjXIScyx3B+FmoU6HSxOJREsikRibSCTGjhrV538/VJcjjwwvlzf5CPY98sgjOdLPiytj3XoOUdLHgP7AG8D1wNmSPm9md8X1hxJuNLkCeETSXWb2l7jueOA+MzskbZ/XAJdLOsnM3pO0DaFHeD+wtZk9EYuOATa4NdDMlkqaBtyXS/3N7KW4n96U+yXwlpk9I2lcLsd1efbGG8WuQUVpbW2lsbGRjjhRQU+Dfd/w8+LKXHeuIUJ49OIEM1sLrJI0HpghaQbhsYungVPjYxcTCOn1WwHrgAcIQ5zpzgW+Dzwr6V3C3ZznA5vE7bcF3gX+CWzUqzOzn2eqdK6p9rmWSzneMuCybOtd5aiWoe4FCxZsFNvUk2DftjhhQqF+b9UeKeXyr8sG0cz6d7LuObI8shCHV7u84cTM3iMk3/9fhtUHZliGmU3Osrwu5XNOqfbdKDc4w7L5wPxctneuVGUL8PVgX1dtfOo257Kolh5IXV0d7RkmKuhusG+yZ1gtvzdXebxBdOXnoIOKXYOK0tTUtME1ROhZsO9Bfl5cmfMG0ZWf884rdg0qSvLGmWnTpvUq2Pc8Py+uzHmD6JyjoaGh6pPtnfO0C1d+DjssvFxJOeywwzjMz4srY94guu4phaT5VavCy5WUVatWscrPiytjPmTqcudJ8865CuYNYrGU40PfpZI0Hx8Az/o79Nv+nXM9UBFDppK2lnSjpL9JelbSHEm7xNccSS9IWixptqRhkmoktUp6RtJCSQ9KGhz3ZZJ+nbLvD0n6Z0y4SD3mpyStlXRUob9v0XjSvHOugpV9D1GSgNsI+YsT4rIxwDDCvKPfScZGSTqAMDH5icBrZrZbXF5PmHoOUoKAzWwVGYKAJfUHfgDc0+OKl2MvplSS5i+5JLyffnrhjum6NH78+GJXwbleKfsGETgAWGNmP0suMLM2SSeSIcAYQFIjKROFm1kibZ/JIOCb+SAIOHUauqnALYRoqurR1LThNUQoTtK8N4Ql6XQ/L67MVcKQabag4s4CjH8JnCnpYUnfl7Rz2vobgQmSBhKCgP+SXCFpO+AI4GdUG0+ad85VsEroIXZb7EHuAHwOOBh4VNI+ZrY4ru8sCHgGcGYMMM56jIoNCC6FpPnk77Mch50rmM9l6spdJTSIi4BMN7YsAvbPtpGZrSTEUd0qaR1wOLA4pUgyCHgcsGXK8rHAjbExHAIcLul9M7s9df+JRKIFaAFobm627n0l55xzhVYJQ6b3AwMkfT25QNKngBeAfSV9PmX5oZJ2k/RpSVvEZR8GPk5a+DBhWPV7ZvZM6kIzG2lmdTFq6mbg5PTG0Lm+kI9Ue+dcdmXfIJqZEa7p/Vd87GIRcAHwCjAemCrpr5KeBSYDrwM7An+U9AzwJPAY4SaZ1P0uMzMPAnYlIZlq397ejpmtT7X3RtG5/FFoT1xfam5utqlTpxa7GpUjj9cQy+X6bqZUe4ABAwaw9957F6FGG2uLEyaMGTOmT4/j1yhdD2S/4SNFJVxDdNXmmGOKXYOCK4dU+6FDhxa7Cs71ijeIrvycfHLedlUuvY18pdo757Ir+2uIrgp1dGw4OUAVaGpqoqamZoNlPUm170sdHR10VNl5cZXFG0RXfg4/PLyqSENDAy0tLdTW1iKJ2tpaWlpaSirU9/DDD+fwKjsvrrL4kKlzZcJT7Z3rW95DdEEpBP8651wReQ/RefCvc87hDWLpKObzcKUS/JurrgKCM/E7MZ1zXfAG0ZVf8O/WWxe7Bi6DyZMnF7sKzvVKwRtESVsTEiM+BawGlgDfBhLAKWbWHMtdATxmZjMlzSRM1P0WYcaB75jZvFhuPrADUBuncUPS7cDBZjY4Sx2+Cvxf3JeAX5rZJfE4xwDDzOydWPYy4BRgqJktl7QEeAdYC7xvZmPz8ospZg+mVIJ/XVnzBtGVu4LeVJOSbj/fzHY0s48D5xDS7V8HTo2TbWdyhpmNITSe6VmEbwKfjsf4KLBNJ3U4LO7jc2b2CWAPQkOb9ALwpVi2HyGA+OW03RxgZmPy1hgWW1NTCPpNVYzg31wtXx5erqQsX76c5X5eXBkr9F2mGdPtgZeAfwLzgBO62MfDwHZpy24EJsTPXyHEOmVzNnC6mb0Sj/+umaVeKLsBODZ+Hgc8BLzfRZ3KW7kF/x51VHi5knLUUUdxlJ8XV8YK3SB2lmIPcDFwmqT+nZQ5FEiPW5oHfDZuNwG4qRd1+CswNMZDHUdobFMZ8HtJj0tq7GQ/5aWhAZYsgXXrwnupNobOOddHSuqmGjN7UdIjwMQMq38k6YfAVkD69P5rgQcJPbtBZrakszT7HNxKaFj/Ezgpbd2nzewVSVsB90p6zsweSN9BfX19I9AI5ZOo4Jxz1azQPcRFwJ5dlLkIOJON63YGsBNwLnBdhu1uBJqB2akLJTVJapPU1o063AhcCNxrZutSV6QMtb5OuB66V6YdJBKJlkQiMTaRSIwdNWpUF4dzLjMPBXaucArdIGZLt69N/mxmzwHPEsJ9NxAbp8uAfpIOSVv9J2A64Rpg6jbT4g0wyZC26cAP492uSBog6ZS0bZYC04ArU5dL2lTSZsnPwOeAhTl+d+e6xUND4WJhAAAgAElEQVSBnSusgg6ZmplJOgKYIeks4F0+eOwiVRMhyT7bPr5PeGzintTlwCU51GGOpGHAffGuVwN+maHczzNsPgy4LQ7Hfgi43szu7uqYbkO9HUI+4PXXAfhDhQ9FZwoF7ujoYMqUKVxdghMmvB7PSz4vEXi0lSukgl9DjEOOmRJeR6WUeYqU3quZTU7bxy3ALfHzuCzHyfgMYlx3LXBthuWTNy4NZlYXPy4Hds+2X1cYf9hqq2JXoSDKIRQ41VZVcl5c5Sqpm2pcdej1X/0vvRTet9++13UpZeUWCvxSPC/bV/h5cZXL0y5c+Zk0KbwqXDmEAqeaNGkSk6rgvLjK5Q2icyWqHEKBnaskPmTqXAnzUGDnCsd7iM455xzeILruam0N6Rj9+oV3fybOOVchfMjU5a61FRobQ3gwhMioxjidayGH9U47rXDHcjk7zc+LK3OKEYKuDzU3N9vUqVPzt8NiPZC+YEHm0OABA2Dv9Olli6QEH0dwzhVdTpNb523IVNLWkm6U9DdJz0qaI2kXSSZpakq5KyRNjp9nSnoxzjX6lKSDUsrNl7RUKbN0S7pd0spO6vBVSQslLYp1OD3lOB3Jadfissti3YbEn5dIeibW5bEs+99e0h8kLY7HOLUXv7Lyk+2B8EI/KN7R8UEv1ZWMRCJBIpEodjWc67G8DJmmBP9eZ2YT4rIxbBj8+3Mzey/D5meY2c2SDgBagJ1T1iWDfx/sZvDvK5IGAqkPRSWDf2d1EfzbWcLp+8BpZvZEbFwfl3SvmT3byTb5V6xeUF1dGCZNV1tb2Dole8jeGywpJ50UgmFKcdIA53KRrx5iVQT/mtmrZvZE/PwOsDhDnStXUxOkPShOTU1Y7pxzZS5fDWLVBf9KqgM+Cfylq7IVo6EBWlpCj1AK7y0tHibsnKsIBbnLtJKCfwEkDSZMLv5tM3s7U5mKDQhuaPAG0DlXkfLVQ6zI4N94E01bfP1PPO4mhMaw1cyyDuF6QLBzzpWXfDWIFRn8a2YvJY9hZj+LNw9dAyw2s5/k9Jtx+XfuueHl+lxrayt1dXX069ePurq6TsOJzz33XM718+LKWF6GTKso+PfThDtXn0npmZ5jZnO6qp/Lo4MPLnYNqkJrayuNjY10xEdc2tvbaYwTMWSaX/VgPy+uzPmD+QWQ9wfzq11b/FtkzJjOy6WpqGu5BbBgwYKMYcQDBgxg7wwTMaxcGR4RHjw4azZ3UfhjII4cH8z3qdtc+fl2HHjwf+j6VKbGsLPlL7zwAgBjuvmHinOlwhtEVzW8p9A9dXV1tGeYiKG2tjbj7zLZA/ffsytXnnbhnMuoqamJmrSJGGpqamjyiRhchfIG0TmXUUNDAy0tLdTW1iKJ2tpaWlpaPLDYVSwfMnXOZdXQ0OANoKsa3iC68nPRRcWugcvgomo6L62tMG0aLF0KI0aE+Xz9D4ey5w2iKz/77lvsGrgM9q2W81IqQdku7/waois/f/5zeLmS8uc//5k/V8N5mTZt4zzOjo6w3JU17yG68nPOOeHdb+8vKefE81Lxj10sXdq95a5s5NRDlLS1pBsl/S0m0c+RtIuky2NC/TOSHpU0MpZfkpJEvzZOjr1Q0u9i0G9yv5+QdL+k5yX9VdJ5yhJnIWm+pKWp6yXdLmll/DxO0p0ZtvuSpKdjHR6T9JlOvucx8fstknR9ljIDJN0k6QVJf4kxUM65ajFiRPeWu7LRZYMYG6DbgPlmtqOZfRw4hxDJtC0w2sx2A44gJNynWxUnxx4F/Av4ZtzvIOAO4GIz2wXYHdgXOLmT6rxJmE+U2LBuk8N3nAfsHicBPxH4RZbvuTMhZPjTZvYJNp6HNWkKsMLMdgIuBX6QQx2cc5XCg7IrVi49xAOANWb2s+QCM2sD/g28moxRMrNlZraii309zAcJ8xOBh8zs93H7DuBbwFmdbH8jIc8Q4CuEfMONSPqUpCcl7WBmK+2DCVs3JUz6ncnXgZ8mv0OMgcrkS3wQU3UzcFC2Xq1zrgJ5UHbFyqVBzJZEPxv4QhyK/LGkT3a2k5h6fxChVwjwifT9mtnfgMGSNs+ym3nAZ+O+JgA3ZTjOvsDPgC+Z2d/jsiMkPQfcReglZrILsIukhyQtkHRolnLbAS/F+r4PvAVsmaWsc64SNTTAkiWwbl1498awIvT4phozWyapHjgwvuZJOtrM5qUVHRSjkuoIDeC9cXkyoinj7rMsXws8SBiuHWRmS9I6Z7sCLcDnkoG/sa63EeKdPksICM6UU/MhYGdgHDAc+JOkUWaWPgycqTe4UX3r6+sbgUbwlIW8mzGj2DVwGczw8+LKXC49xKxJ9Ga22szmmtkZwEXAlzMUWxWv39UCHyZeQ4z7HZtaUNIOwEoze0fSPbH3mX7N70agmdBDTfcqIYsxY2/VzB4AdpQ0RFJT3H8y13AZ8FszW2NmLwIJYOcs5baP9f0Q8BHCtdENJBKJlkQiMTaRSIwdNWpUpuq4nhozptvRT67vjRkzpuSTLroTeOyqTy4N4v3AAElfTy6I1+j2l7Rt/LkfMBrYeGr8yMzeAk4BTpe0CdAKfEbSwXEfg4DLgR/G8ofEm3H+O21XfwKmAzdkOMybwOeBiySNi/vdKXmNT9IehEb5DTObFvef/D/4dsL1UuIdsrsAf89Q7g7ghPj5KOD+lGuUrhDuuy+8XEm57777uK+Ez0sy8Li9vR0zWx947I2iS+pyyDQm2R8BzJB0FqEHtgS4G/iJpAGx6CPAFV3s60lJTwETzOzXkr4ENEv6KdAf+HUO+zDgkk7WvybpC8BcSScShkC/KmkNsAo4NksDdg/wOUnPEoZmzzCzNzKUuwb4taQXCD3DCRnKuL70/e+Hd09o75a+Hrpvi8HNpdpLzBR43NHRwZQpU7j66quLVKvCqPhnQ/Mkp2uI8XrcMRlWNWcpX5fyeXDaui+kfH6G0GDlUoeM5ZL7N7P5wPz4eSnhph2Av5DDoxGxkfxOfHVW7l3g6Fzq7JwrHd0NPHbVx2eqca5K9HUvodQDgrsbeOyqj89l6pyrCh547LriDaJzrip44LHrivwGyb7X3NxsU6dOLXY1KkciEd7r64tbD7eBRDwv9X5eXOnJaTYxv4boyo//g1uSvCF05c6HTF35+d3vwquctbZCXR306xfeK+BZuN/97nf8rtzPi6tq3kN05efHPw7vX/hC5+VKVYUmrv84npcvlOt5cVXPG0RXfYo9t+yCBZD+7FtHB0yZAuX6gLg/tuAqQEkOmWYLJO7BfmZKejE5F6mkP8flwyTdKemp5P7j8jpJJunClH0MkbRG0hXx5+/EbZ6WNE9Sbb6+t6sS2R4E9wfEnSuqkushpgQSX2dmE+KyMcAw4Pke7PIMM7s5bdn3gHvN7LK4/9Ep6/4OjAfOiz8fTZiIPOlJYKyZdUj6BmHu1WN7UC9XLMXuzdTVhWHSdLW1xa+bc1WsFHuI2QKJn4w9sickPRPnQU326hZLulrSIkm/jxOFd2YbQmpFcv9Pp6xbBSyWlEziOJaUZA0z+0MMMwZYQIiKci53nrjuXEkquR4i2QOJ3wWOMLO3YxrFAknJsOGdgePM7OuSZgNHArPiuh9JOjd+XmRmDcBPgZskfQu4D7g2NT+REDE1QdI/CBN9vwJsm6FOU4C5Pf6mrmd+/eti16B3kjfOTJsGS5fCiBGhMSzjG2oAfl3u58VVvVJsELMRIdbps8A6QnL9sLjuxdiLhNCY1qVst9GQqZndE7MXDwUOI/Q+U0ML7yYECb8G3JSxMtLxhDzH/TOt94DgPrT99sWuQe81NJR9A5hu+0o4L66qleKQabZA4gZgKLBnzCZ8DRgY16XejbCW3GKt/mVm15vZJOBR4LMp694jNKynAbekbxszHKcBXzSzjHdCeEBwH7rppvByRZMpaPemm27iJj8vroyVYg/xfkJP8OtmdjWEQGKgFnjdzNZIOiD+3COSDgQWxBtjNgN2BJamFfsx8EczeyPmCye3/STwc+BQM3u9p3VwvXDVVeH9WL+XqRiSQbsd8TnKZNDuiBEjGDZsGMf6eXFlquQaxE4CiS8ALpf0GNAGPJfjLlOvIQLsReiBXiHpfUIv+Rdm9qikupR6LGLDu0vX7w8YDPwmNpRLzeyLOX9BV3Uqbcg8W9BuIpHg1Vdfrbjvm28eNVW6Sq5BhE4DiffJssn6MUkzuyTl8+Qs5X8UX+nHXZK6r5TlM4GZ8bPHtLuqli1Q14MCXLkryQbRuUpSaT2CbEG7AwYMYMyYMRX3fV31KMWbapxzJSxb0O7IkSOLVCPn8sN7iK783Jw+8ZArpGSg7rRp01i6dCkjRoygqamJQw45pMg1c653vEF05WfIkGLXoOo1NDR40ryrOD5k6srPzJnh5UrKzJkzmennxZUxbxDdhsohuNYbxJLkDaIrdz5k6j5QocG1zjmXC28QC6nUH1gul+DatjhtbXd/n/44gHOuEz5k6j7gwbXOuSpW0B6ipLXAMymLvmxmSyTtRQja3Q54B3gVOMvMnonbNQLfidu8DXzHzB5M2e9QQkTTt8zs5ynLlxDCfJenLBsGXANsD2wCLDGzw7vxHcYAVwGbEyYSbzKz3GY0LvUeSrkE1yZ7hqVUJ+dc2Sv0kOmqmFSxXmygZgMTzezPcdlnCBNuPyNpPHAS8BkzWy5pD+B2SXuZ2T/ibo4mhPUeR5h4uzPfA+41s8visUZ38zt0AF81s79K2hZ4XNI9ZvZmN/dTepqaNryGCKUZXDtnTrFr4DKY4+fFlblSGDL9FnBdsjEEMLMHzez2+OOZhEzD5XHdE8B1wDdT9nEcIappuKTtujjeNsCylGM9namQpB0lLZD0qKTvSVoZyz9vZn+Nn18BXifEUpW/hgZoaQk9Qim8t7SU3g01NTUbJ867oqupqdloBhvnykmhG8RBktri67a47BPAE51s8wlCNmGqx+JyJG0PbG1mjxB6ml1lz/wUuEbSHyRNi728TC4DLjOzTxGGYzcSh3o/DPyti2OWj4YGWLIE1q0L76XWGAJceWV4uZJy5ZVXcqWfF1fGCt0grjKzMfF1RKYCkv4iabGkyzrZj4Dk1PoTCA0hwI2E3mJWZnYPsANwNfAx4Ml4DTLdPsBv4ufrM9RzG+DXwNfMbF36+vr6+sb6+vrH6uvrH1u4cGFnVXLdNXt2eLmSMnv2bGb7eXFlrBSGTBcBeyR/MLP/BM4DPhIXPUvIL0y1R1wOoQGcHG+guQPYXdLOnR3QzP5lZteb2STgUeCzkpqSvdeuKixpc+Au4FwzW5CpTCKRaEkkEmMTicTYUaM2SpRyFSJTcrxzrjyVQoP4U0KDtm/KstQLET8EfiBpS1h/l+dk4EpJ9cCmZradmdWZWR0wndBrzEjSgZJq4ufNCDfvLDWzacneayy6ADgyfp6Qsv2HgduAX5nZb3BVK5kc397ejpmtT473RtG58lT0B/PN7B+SjiU0etsRblJZTrgbFDO7Iy7/syQjPJZxvJm9KukkQuOU6hbC0OmF8eenJSWHNGcTHum4QtL7hD8IfmFmj2ao2reBWZJOI/QG34rLjwE+C2wpaXJcNtnMuuxZVpO+TE2fER/M/3aRJzrIlhw/ZcoUri6liQz6iOceukpT0AbRzAZnWb4A2L+T7a4iPPuXvvyCDMueBj4eP9dl2eWPuq4tLwN7m5lJmkC4kQczmwXMymF7V+GyJcdnW+6cK21F7yGWsD0JPUkBbwInFrk+ZaUQvYe+P0LnsiXH19bWVmXvqRq/s6sspXANsSSZ2Z/MbHczG21mnzWzF4pdJ1dasiXHN5XaRAbOuZx4g+jKzyWXhFeRNTQ00NLSQm1tLZKora2lpaWlaoNzL7nkEi4pgfPiXE/JzLou5XqlubnZpk6dWuxqVA6fy7QkJW+k8qFTV4KUSyHvITrnnHN4g+h6q7U1pGT06xfe/Rk851yZ8rtMXc+1tm6YjtHeHn6G0pwD1TnnOuENYrXoi4fYFyzYODy4owOmTIG+fDD96RhQku/v5Ne+emXQoEHFroJzvVLwIVNJa1MSL9ok1cXle0maL+mvkp6QdJek3VK2a5T0XHw9EjMTU/c7VNKaOHtN6vIlkoakLRsm6U5JT0l6VlLOQW6S/kvS45Keie8H9uT3UBGyPYDe1w+mjx4dXq6kzJ07l7lz5xa7Gs71WDF6iOUeErwc+IKZvSJpFHAP0FUGY/H1Re+nri4Mk6arrfXelnOu7JTKTTWlGhJ8tKSFsSf5QCz7ZAwGhpDUMVDSgFy/aEVpato4qLemJizvSxdeGF6upFx44YVc6OfFlbFiNIjlFBJ8PnCIme0OfDHD+iOBJ82sOievbGiAlpbQI5TCe0tL399QM29eeLmSMm/ePOb5eXFlrCSGTNNJ+guwOfB7Mzs1WzGyhwRfA/wk2/7N7B5JOwCHAocRQoJHmdk/04o+BMyUNBu4Na2OnwB+AHwu0zHq6+sbgUbo2+SHomto8DtKnXMVoVSGTEsyJNjM/gc4F9geaEvJZBxOiJ36qpn9LdP+PSDYFYuHFjvXM6XSIJZkSLCkHc3sL2Z2PuFmmu0lfZSQj3i2mT2Ul2/vXJ54aLFzPVcSzyGWcEjwj2JPU8A84ClgGrATcJ6k82K5z5nZ6738NbgMMg03/79FiwD4biUPRfdQMUOLF8XzUsqXCHyeVdcZn9y7AHxy754r5X9cS9Ef//jHrOv23z9rBnfV8AaxauU0uXdJ9BCdy8b/AeseDy12rudK5Rqic7k7++zwchspZmjx2Wefzdl+XlwZ8x6iKz8PP1zsGpSsZDjxtGnTWLp0KSNGjKCpqakgocUP+3lxZc4bROcqTENDQ0EaQOcqjQ+ZOuecc3iD6Fzp8LBl54rKh0xd+Rk+vNg1yL8KCFseXonnxVUVbxBd6cn12cNKekaxWGHLeTQr+SEf58UfEXFF4EOmzpWCYoUtO+fWK2oPMQYDXwrsDawA3iPMW7oC+C3wYiy63MwOlnQB8HXgn8CmwDPAuWb2bNzffOB0M3usi+OeDvw38D6wFvixmf0qbr8DUGtxCh9JtwMHm9lgSbWE1Iv+wCZAs5n9LA+/Cpeqq97Bt78d3mfM6POqFEwFhC1/O56XGZV0XlxVKVqDKEnA7YRg4IlxWS0hd3AF8CczG59h00vN7JJY/ljgfkm7ZYhuynbc/wH+C9jLzN6W9BHgyylF3gQ+DTwYJ/LeJmXdq8C+ZrZa0mBgoaQ7UgKDXSG0tRW7BvnX1LThNUQoTNhyHrVV4nlxVaWYQ6YHAu+l9rDMrN3MmnPdgZndBPwemNiN454DnGxmb8d9vGVm16Wsv5EPkjK+QkoOopm9lxIGPAAfcnb5UqywZefcesUcMv0E8EQn6/dL5hICvzGzbH8qPwF8LJcDxqinzbJlGEbzgKsl9Sc0jI2EbMbkPrYnxD/tBJyRrXdYNQHBLn88bNm5oiqZu0wl/RT4DOE64hlkHzLdaNPuHAboKt5jLfAgcCwwyMyWhNHdwMxeAkZL2ha4XdLNZvZa+k4SiUQL0AIh7aIbdXTOOVcExRzyW0RIvQfAzL4JHAQM7eZ+PgkszrZS0rWS2iTNicOk/5a0Qxf7vBFoJmQnZhR7houA/bpZX9dbu+wSXq6k7LLLLuxSIeeltbWVuro6+vXrR11dnQcsV4li9hDvBy6S9A0zuyouq+lsg3SSjgQ+B5yWrYyZfS1t0XTgp5KOjTfVbA5MMLOWlDJ/iuVuSDvecOANM1slaQvCzTc/6U6dXR60tHRdxhVcS4Wcl9bWVhobG+mINzi1t7fTGCdJ8DliK1vRGkQzM0lfBi6V9H+ERyn+DZzZxab/K+l4wmMXC4ED0+4wvUvSmvj5YTM7Om37q4DBwKOx3Brgx+l1Ay7JcOxdgR9LMsLw6yVm9kxX39W5UubXuDe0YMECVqc9/9nR0cGUKVO4ukwmSciXasvQLOo1RDN7lQ/u6Ew3P0P5C4ALOtnfuByOaYRnHX+Y6/ZmNji+3wuM7uoYro8lpzSrkB5JpXj++ecByn7YNL0x7Gq5qxwlc1ONczmL//C6/MhXLyDZ0yz3XkVdXR3tGSZJqK2tLfvv5jrnz9E551yKpqYmamo2vJ2hpqaGpjKaJMH1jDeIzjmXoqGhgZaWFmpra5FEbW0tLS0tfkNNFfAhU+ecS9PQ0OANYBXyBtGVnzFjil0Dl8EYPy+uzHmD6EpXaytMmwZLl8KIEWGi64aGykq5qCCecuHKnTeIrjRVQIK8c668eIPouqdQD3F3liD/ne+En3fdtTB1KaQyvq3/+OOPB2DWrFlFrolzPZO3u0wlDZN0vaS/S3pc0sOSjpA0TtJbcT7RNkn3xfIXSHo5LvurpFslfTxlf/Mlje3imDMlvZiy7z+n1OVOSU9JelbSnLi8TpJJujBlH0MkrZF0Rfz5O3GbpyXNixmNybI7x/3+LX7HP0j6bL5+hy5FZwnyyZcrKcuWLWPZsmXFroZzPZaXHmKxwn6jM8zs5rRl3wPuNbPL4r5TZ5f5OzCeDyKdjiZM0p30JDDWzDokfYMwo82xkgYSYp9ON7M74n5HAWOBB7pR3/JWqB5MZwnydXWFrYtzrirkq4dYrLDfbLYB1v+pamZPp6xbBSxO6X0eS0qqhZn9wcySseULgOHxcwNhbtQ7UsouNLOZeaivS9fUFBLjU5VZgrxzrrzkq0HMKew3vqZ1Ui7nsN8UP0rZdzKj5afANXFIc1rMLkx1IzAhplesBTKG/AJTgLnxc1ff0eWTJ8g75wqsT26qKVDYb9JGQ6Zmdk/MPDwUOAx4Mg5vJt0NXAi8BtyUsSIhUWMssH+W9bcBOwPPm9lX0tfX19c3Ao3gaQI9li1Bfp99Cl8X16V9/Ly4MpevHmLBw3672pGZ/cvMrjezScCjwGdT1r0HPE7IUbwlw3EOBqYBXzSz5N0b6d/xCGAy8B+Zjp9IJFoSicTYRCIxdtSoUZmKuJ6aPj28XEmZPn0606dP93BdV7by1SDeDwyMN6Ek9TTs94ZsZczsa2Y2xswO72JfB0qqiZ83A3YElqYV+zFwppm9kbbtJ4GfExrD11NWXQ98WtIXU5Z16zs6V+mS4brt7e2Y2fpwXW8UXTnIy5BpEcN+IVxDPDfl572APYErJL1PaPR/YWaPSqpLqfMiNry7dP3+CAHCvwk3z7LUzL5oZqskjQd+ImkGYbj1HeD7XXxHl29HHhneb9moc99rPrzdc4sWLeJf//oX69at22B5tYbrQvlHYVWbvF1DLFLY7+Qsq34UX+nllwAbjV/GO0Vnxs8Hd3K854BOe6euAN54o+syruDWrFmzUWOY5OG6rhz4TDXOpfC/6Htu3LhxLFiwIGPj5+G6rhx4HqJzLm9Gjhzp4bqubHmD6JzLm2HDhnm4ritbPmTqys9BBxW7Bi6Dg+J58XBdV668QXTl57zzui7jCu48Py+uzPmQqctda2uYWLtfv/Duz5Y55yqI9xBdbkopsPeww8L73Lmdl3MFdVg8L3P9vLgy5Q1iMZXTQ+CdBfYW+oHrtrbw3te/P39MoFtWrVpV7Co41ys+ZOpy01lgr3POVYB8BQQPAy4F9iYEAr9HCNZdAfwWeDEWXW5mB0u6APg6YYq3TYFngHPN7Nm4v/mEIN7HOjnmTEISxVtxUYeZ7Rvrcg2wPbAJsMTMDo/Ttr0IfN/Mzov7GAK8CvzczL4l6TvAfwPvx7qdaGbtuWzbo19cOfVAOgvsLfT3SPYMy+n355wreb3uISpM+Hk78ICZ7WBmexKmcEsG6/4pTsg9Jm1atEvjsp0JEUz3S+puOsYZKfveNy77HnCvme1uZh8Hzkop/3cgNYbqaDacz/RJYKyZjQZuJjTquW5b2Tyw1zlX4fIxZHog8J6Z/Sy5wMzazaw51x2Y2U3A74GJeajPNsCylH0/nbJuFbBY0tj487HA7JSyfzCzeNcIC/igUe9y24pXSoG948eHlysp48ePZ7yfF1fG8jFk2lWS/H6S4l0Q/MbMsnUpngA+1s1jpyZdLDKzBuCnwE2SvgXcB1xrZq+kbHMjMEHSP4C1wCvAthn2PQVIv10u120rMyA4W2BvoZ1+erFr4DI43c+LK3N5v8tU0k+BzxCuI55BGDLN5c9G9eBwZ5jZzakLzOweSTsAhwKHAU9KSk24uBu4kBDfdFPGioRIqrGEa5Sputw2KZFItAAtAM3NzZbrF3LOOVcc+RgyTU+S/yZwENDd64GfBBZnWynpWkltkuZ0tSMz+5eZXW9mk4BHgc+mrHsPeBw4DdgoUE/SwcA0QkDwBrdQdrWtK5Bx48rrkZVeKKf0+XHjxlXOaIirSvnoId4PXCTpG2Z2VVzWrSR5SUcCnyM0NBmZ2ddy3NeBwAIz65C0GbAjsDSt2I+BP5rZGzEEOLntJ4GfA4ea2etZDpFxW+fyLZk+3xEnQ0imzwM+V6hzfaDXDaKZmaQvA5dK+j/C4wr/Bs7sYtP/jUOTmwILgQPN7J8p6++StCZ+ftjMjs6wj9RriAB7AXsCV0h6n9AD/oWZPRofnUjWeRGZ7xD9ETAY+E1s7Jaa2RfTvm+2bV0e5NLDmBEfzP92hfdGMmULlnL6fFs8L6XSS/T8RdddebmGaGavEh61yGR+hvIXABd0sr9xORxzcpZVP4qv9PJLgFEZls8EZsbPB6evz3Vb5/ItW8q8p8871zd86jZXUnL6qz72QCq9B1BXV0d7hskQSjV9flyVnBdXubxBdOXnmGOKXYOCaGpq2uAaIpR2+vwxVXJeXOXyBtGVn5NPLnYNCiJ548y0adNYunQpI0aMoKmpqWRvqDm5Ss6Lq1zeILryk+wxpU8lV4HKKX0+2ZOtqYLz4iqTN4iu/Bx+eHj3a1Ul5fB4XvwaoitXHv/knMtNa2tIPenXL7yX8CQBzvWE9xCdc11rbYXGxg+Gq6dH6DYAACAASURBVNvbw89QGvPbOpcH3iC6rpXIg9brxQfAS65elWzBgo3DoDs6YMoUSE4S4OelelXIMHlJDJlKGibpekl/l/S4pIclHSFpnKQ7U8rcKekpSc8m5zSV1E/S5ZIWSnpG0qOSRmY4xjhJb8X5UJOvgyVtL+kPkhZLWiTp1E7quTLt58mSrsj378O5kpNtMgCfJMBVkKL3EFMChq8zs4lxWS3wRWBFStFk8O9lsczouPxYQgTTaDNbJ2k4Yeq4TDZK3pC0DXCamT0R5z59XNK9ZvZsnr5i+Su1v/5mzgzvkycXsxbVpa4uDJOmq61d/9/HZD8vrswVvUEkS8Aw0CxpXEq5bQghwskyT6csf9XM1sXly+iGOO3cq/HzO5IWA9sB3iCWKv8Ht/Camja8hgjhsZeUSQIm+3lxZa4Uhky7ChhO+ilwTRzenCYpGcw7G/hCHAL9cUysyGa/tCHTHVNXxgnAPwn8Jcv2g1K3J/RaXaEtXx5ernAaGqClJfQIpfDe0rLBDTXLly9nuZ8XV8ZKoYe4gQwBw0D24F8zWyapntDTPBCYJ+loM5uXYfdZw4olDSZkHH7bzN7OUr1VZjYmZZvJhCDhjdTX1zcCjVA6s/9XjKOOCu+lNpRb6RoaOr2j9Kh4Xvw5RFeuSqGHmHPAcLbgXzNbbWZzzewM4CLgy/GmnGRvLmOjlSRpE0Jj2Gpmt8Zl26ds/z/d/VKJRKIlkUiMTSQSY0eN2igow7kulVM4sHOVoBQaxPuBgZK+kbJso7mfJB0oqSZ+Xh/8K2mP5PCppH7AaKDdzG4zszHx9Vi2g8ebeq4BFpvZT5LLzeyllO1/lm175/pCMhy4vb0dM1sfDuyNonN9p+hDpt0IGM4W/HsocLWkAbHcI0C2RyH2i9f+kr4P/AOYBDyTsu4cM5vT6y/nMurtEHI1BASXWzgwlF5AcGd8WNdlUvQGEXILGDazbMG/dwN353CM+cBHsqxWjvUcnPbzTDwg2PUBDwd2rvBKokF01aXXf53fdFPYz7HH9r4yJarcwoEBborn5dgKPi+usnmD6MpPFfyDW27hwOANoSt/pXBTjXPd89JL4VXBGhoaaGlpoba2FknU1tbS0tJS0tmIL730Ei9V+HlxlU1mVuw6VLzm5mabOnVqsatROZI3bZTo0GG1St5MU6pDuq6q5XSfiPcQnXPOObxBdM6VOw8udnniN9U458qXBxe7PPIG0blKU6wH44sREJxLcHEl8+u1eeUNois/p51W7Bq4DE4bPrzwB/XgYpdHJdMgShoGXArsTQgGfg/4Yfx8upmNj2WuAbYHNgGWmNnhcQ7TGYS0CwPeBY4xsxfTjjEO+C2Quvx0IAH8CtgaWAe0JIOIM9RzLfAM4Xf3IjDJzN7s9S/A5e4LXyh2DUpbkXoNRTkrOQQXO5erkripJk6wfTvwgJntYGZ7EqZyS/+T83vAvWa2u5l9HDgrLj8W2BYYbWa7AUcA2RqpP6VM2j3GzO4D3gdOM7NdCQ3yNyV9PMv2q+J2o4B/Ad/s2bd2PZZIhJcrKYlEgkShz0tTUwgqTpUWXOxcrkqiQST07N5LTZUws3az/9/encfJUdV7H/98wxIYQR4koizJNJtzlVwIGCWgCAgoKIjIEnCI5nkiA7J4UUBlj8DgVYggmzKANwojBBEUkEUJBOFC2GRIWBwgkIFAZJdt2PN7/jinSaWneqZnprfq+b1fr35116lTVaeqIGfOqVPnZ2cV5FsLWJTIMy+RvtjMlsT0RWb2cqkHN7PFZvaP+Ps14GFgnRI2vaPEfK6cDjggfFxdOeCAAzig2velhMDFzpWqXrpMNwb+UUK+c4BZkg4BbgT+x8yeAS4DbpO0NTAbuNjM7iuyj8KIF3uY2YL8gqQcsBlwZ38FkbQcIW7jhWnrPUCwc1UyQOBi50pVLy3EZUg6R9L9ku5OppvZDcD6wPnAfwD3SfqomS0CWoCjCM8AZ0vavsjuC7tMk5XhKoRAwYeZ2atFtl85VqgvAh8B/paWyQMEO+dcttRLhfggsHl+wcwOJrS+PlqY0cxeMrPfm9kU4G7gCzH9bTO7zsyOBE4Bvi5p90TU+4n9FUDSCoTKsNPMrohpYxPbHxizvmlmE4BmYEX8GaJrEJ2dneRyOUaNGkUul/NgxG7EqZcK8SZgJUnfTaQ1FWaS9EVJTfH3qsAGwJOSNpe0dkwfBWwC9JjZlYmW4D3FDh4H9VwIPGxmv8inm9lTie1/ndzGzF4BvgccEStT5zKrs7OTtrY2enp6MDN6enpoa2vzStGNKHXxDNHMTNLXgdMl/RB4HngD+FFB1k8DZ0t6j1CZX2Bmd0vaCThf0uiY7y7g7CKHK3yGeDLwL2AKMD+x7mgzu3aAct8n6X7CiNiLSjpZN3zHHlvxQ4y0575z587tE3y4t7eXadOmcX6JL7i//HIYxzbSrl2pfNLz+lcXFSKEkZ6EiiXNnJjnVODUlG2vB64v4RhzgNWKrC5pNnQzW6Vg2V+Kq7Yddqh1CRpOYWU4UHqa1VdfvVzFca4m6qZCdK5k+SnCJkyo2CFG2l/zuVyOnpQX3Jubm0u+Fl3xvkyo4H1xrpLq5Rmic6U77LDwcWXT3t5OU8EL7k1NTbQP4gX3ww47jMP8vrgM8wrROUdraysdHR00NzcjiebmZjo6Omj19/vcCOJdps45IFSKXgG6kcxbiG7k8YCyzrkU3kJ0I4sHlHXOFeEVohu8Wr9n9sor4Xso5RjpAWUr6JRS7ssIG73rssUrRJc9qxV7lbQEHlC2YrYazn1xrg6UrUKU9HrypXVJU4GJZnbIIPaxMG7zQrGAwWZ2ZbFAv2Z2oyQjRLuYEve5PLAYuDMGGW5l6Qw4rwPfNbP7U8oylTAJwNOJ5G8C/6SEYMQNrdZ/5d9+e/jeaqvBb+sBZSvm9nhfthrKfXGuDtRlCzERMPi3ZvbNmNYMfC2R7VYz2yVl8zeA8ZJWNrM3gR1ZtlJ7AtjGzF6WtDPQAWxRpCizCit0SfuyNBjxEknrxmO6ajn66PA9lAqsvX3ZZ4jgAWXL5Oh4X0bapAaucVRllKmkXSXdKek+STfG1h+S1pD015h+HkunTys1YHAx1wFfjb/3BS5J7Of2RPDgucC6gzydYQUjdjXmAWWdc0WUs0JcOREqqQs4MbHuNmCSmW0GXAr8MKafANwW068CxsX0UgIGb508nqQNEusuBfaRtBIh8kWxYL/TCJVnMZMLjrEyIRjxrnF5hqTNBiinqzetrbBwISxZEr69MnTOUd4u03ycQGDpM8S4uC4h0v1ahBiC+edtXwC+AWBmf5GU2tKSdA7weUKr8TMxuViXKWY2T1KO0DpMjVghaTtChfj5fs6pT5cpsEhSC6EV+0VCMOK9zGx2MlNLS0sb0AY++79zzmVBtV7MPws428z+EzgAWCmxzlLylxwwuB9XAaeR6C7Nk7QJcAGwm5m9GNMOTrQE1+5vx2nBiAvzdHd3d3R3d0/s7u6eOH78+EEU2zlXCg9o7MqtWoNqVmPpwJZvJ9L/DrQCJ8cBLvn4MTcBp0j6rpn9Kqb1CRg8gN8Ar5jZ/DgqFQBJ44ArgClm9kg+3czOAc5J5EvdqaTNgX+Z2TOJYMTzBlk2NxxnnFHrErgUZ1TxvuQDGvfGwVH5gMaATz/nhqxaFeJ04A+SniYMZFkvpv8EuETSP4BbgCeh5IDBfQL9mtnl+QUzWwT8MqUsxwNrAOfGSu89M5uYkg/CM8Rkl+pBwIcpPRjxgLw71bnBK0dA43rjo3NrT2ZpPZaunM466yw79NBDU9d5hTh4n46R2e/1gLR15eV4X6oRKPiWW24pum6bbbap+PErwSvEiiopAHxdvoc4kvj/BEOQ/yPCr11dyf9xV43/pssR0Ni5Qh7twjmXOeUIaOxcIa8QnXOZ4wGNXSV4l6lzLpM8oLErN28hOuecc3gL0WXReefVugTD19kJxxwDTz4J48aFycUz3to5rxHuixvRvEJ02dPSUusSDE9n57IRN3p6wjJkulJsyfp9cSOeV4gue66+Onzvuuvw9lOrd0Dnzu0bkLi3F6ZNg3p6qXyQry9cHe/LrsO9L87VSGaeIUp6PX7nJJmkQxPrzo6TiSNppqReSasm1v8ybjMmLn9c0qWSFkh6SNK1kj6R2PdJiW3HSHpX0jKz0UjaM+YtNsuNq5QZM8Inqworw4HSM2LGjBnMyPJ9cSNeVluIzwH/Jek8M3snZf1jwG7AxXG+0e2Ic6nG4MNXEoIP7xPTJgAfA54CHgd2AY6L+9qLMNn4B2Jl+z2Kh5VyWVCrF7hzudBNWqi52ScbcK6GMtNCLPA8MJtlJwpPugSYHH9vC/wv8F5c3g54tyD4cJeZ3RoX3wQeTrT8JhNiICadBPwceGsY5+BGqvZ2KHipnKamkO6cq5msVogA/w0cLmm5lHWPAh+VtDohJuKliXXjgXsH2Hc+wPC6wPvAM/kVMSDwWDO7ZjiFdyNYayt0dIQWoRS+OzoyPaDGuUaQ1S5TzOwJSXcB3yyS5QpgH2ALQgzGwbie0Ap8FpiVT4zdr6cDUwfagQcIdv1qbfUK0Lk6k9kKMToFuJwQV7HQpcA/CM8KlyTiGz4I7NnfTs3sHUn3AocDGwP5YXOrElqYc+L+Pg5cJelrZnZPch/d3d0dQAeEaBeDPzVX1EUX1boELsVFfl9cxmW5yxQz+yfwEGEQTOG6J4FjgHMLVt0EjJa0fz5B0mckFcaMmQH8yMxeTOzzFTMbY2Y5M8sRYjv2qQxdhY0dGz6urowdO5ax/dwXj3Dv6l3WW4gA7cB9aSvMrM/UGTH48O7AGZJ+TBgYsxA4rCDfgxSMLnV1YlbsxZ48uf98rqpmxfsyOeW+eIR7lwUeILgK+gsQ7IbA4yHW5XPprq4uACZMmNBnXVqEe4DRo0czadKkipdtqDy2YsMoKUBwprtMnXPZkFYZ9pfuXC00QpepcyNOPbZc8q3WtLJ5hHuXBd5CdM5VnEe4d1ngFaJzruI8wr3LAh9UUwU+qKbMXnghfI8ZU9tyuGW8EO/LGL8vrv74oBrXoMaMafzKsLMzTAI+alT4zsA7e2PGjPHK0GWaD6px2TNzZvieOrWWpaicjAYQnhnvy9RGvS+u4XmF6LKn0hVird/xy0oA4QIz43uIU/P3pxgfVerqlHeZOldvGjSAsHP1bsgtREmvm9kqknLAw0A3sCJhou2DYrYzgC8CRpgibe8YpWIhMNHMXuhn/ysQIk7sAbwN9AInmNl1cfunzGzrRP4uYHkzG19QJoC5ZnZgIu9mhIm/dzKzG4ocf2os4yEF6QuB1+I5vQx8y8xSor26zKp1CyarAYR9BiGXceVqIS4wswnAJsCngK8TAuuuDWxiZv8J7A78exD7PAlYCxhvZuMJESdWTaxfVdJYAEmfLFam+DmwYN2+wG3xeyi2M7NNgDnAsUPch3PpPICwczVR1i5TM3sPuB3YkFCZLTazJXHdIjN7uZT9SGoC9gcONbO34/bPmlkycv1lhEoXQsV2SYn7FiH801TgS5JWKmW7Iu4A1hnG9s715QGEnauJsg6qiRXZ9sDxwHzgNklbA7OBi80sNSpFig2BJ83s1X7yXA7MBE4jtB5bgSmJ9etJug94FTjWzG6N6Z8DnjCzBZLmAF8hBBMeip2APw1xWzdU115b6xJUXgYDCF87Eu6La2jlqhA3iM/wDPizmV0HIKmF8Azxi8BsSXuZ2ewyHfMl4GVJ+xCeF/Ym1i0GxpnZi5I+DfxJ0saxgt2XEDyY+D2FwVeIN0v6GPAcRbpMW1pa2oA2qM/IBJlW2J3o6kLh1GzOZU1ZnyGa2WZmNj2faGZvm9l1ZnYkIbr914vtQNINkrokXQA8BoyTtGqx/NEs4BwKukvjcV+Mv+8FFgCfkLQcYZDO8XFwzFnAzpJWlXRwPH6XpLUHOO52QDMhXuKJaRm6u7s7uru7J3Z3d08cP378ALtzg3LuueHj6sq5557LucO8Lx5E2NVSxV67kLR5vmKRNIow4KboaEwz+3KsVL9jZr3AhcCZklaM+1hL0n4Fm10J/BxYZqSopI/Gyg9J6wMbAY8DOwD3m9nYGPW+Gfgj8HUzOycxCOeZgc7PzN4kBBX+lqSPlHBJXLlcdln4uLpy2WWXcdkw7ks+iHBPTw9m9kEQYa8UXbVU8sX8NYHzJY2Oy3cBZyfWz5O0JP6+zMx+ULD9scDJwEOS3gLeIDyb/ICZvQb8DCCMlfnAF4ATJb0HvA8caGYvSdqXUIkm/RH4LnBRyjlMlZRs1S4TydTMFku6BDiYMCrWuQE1ahd6PkDwUM8vLYhwb28v06ZN4/w6npCg0jw8VvUMuUI0s1Xi90KgT5+gmV0PXF9k21wJ+38H+GH8DLh9shxm9kdCRVeYZ2pK2lXAVSnpMwmDdgrlCvL5rN3OlYEHEXa15lO3OVdljfoXf38BgkvhQYRdrfnUbc65uuBBhF2teQvRZY+3FurScFtx+WDBxxxzDE8++STjxo2jvb3dgwi7qvEK0TlXN1pbW70CdDXjXaYue047LXxcXTnttNM4ze+LyzCvEF32XHNN+DSazs4Q6WLUqPCdsffvrrnmGq5pxPviRgzvMnWuHnR2QltbCAQMIfxTW1v47V2IzlWFV4gu2xrlJfe5c/sGAO7thWnTICsvpccX8/u9Jz4gytWxQXWZSno9fuckvRnn/XxI0q8ljYqfMyU9IGm+pLslrRe3WShpzAD7X0HSf0t6NO7jLkk7J7a/tSB/l6QHUsrUJenXBXk3k2SSvtzP8afGPNsn0naPaXvG5TmSuiXdH89vwmCuoXOpir187i+lO1c1w2khLjCzCZKWB24iTNw9mqVBgZdIWpcw5VqpkkGB344RJbZJrF9V0lgze6q/oMBF9p0MCnxDkTwQwlbtSwhZBbAPcH9BnlYzu0fS/wVOBXbs96xcea288tLfjdLiyOVCN2mh5ubMnOPKO+8cflx3XW0L4twQDbvL1Mzek5QPCvweBUGBS91PIijwesmgwIRAwHn5oMCnsTQo8BQGkAgKvCNwq6SVzOytItlvBbaWtAKhgt8Q6CqS9w7gyIGO78qsEf/BbW9f9hkihDBXGXop/bpGvC9uRBn2KNNEUOD5hApr19hlOUPSZoPYValBgb8Rf+8KXF2wfj1J90m6JQYmzvsgKDAwhxAUuBgDbgS+DOxGyjynCR4g2JVHayt0dIQWoRS+Ozp8QI1zVTScFmIjBwW+FPgesBpwOHB0wfpOSR8ClgM2T9uBBwiuoJNiYJHjjqttOcqttTXTFeBJ8b4c12j3xY0Yw2khNmxQYDO7ixA5Y4yZPZJShlZgPeD3sSx9eIDgCpo9O3xcXZk9ezaz/b64DCvri/kNFhT4KPq2DJNlf5cQs3FSkQE+zrkK6OzsJJfLMWrUKHK5nAcQdmVT7plq1gSujq9CzCMMsikMCrwofn6Rsv2xwPOEoMAPEJ7PPZ/MYGavmdnPYrzEpC/E/d9PeNZ4oJm9ROguTQsK/M3+TiS2cm8eIM+bwAzgiP7yOefKo7Ozk7a2Nnp6ejAzenp6aGtr80rRlYXMrNZlaHhnnXWWHXqoxxEum/wz2Yy8jjBSlBIPcbjP0+fOnZsaMHj06NFMmjRpWPseCo/TmBkqJZPPVOOyZ401al0Cl2KNKtyXtMqwv3TnBsNbiFXgLUTnyiOXy9GTMoFBc3MzCxcurH6BXFaU1EL0aBfOucxob2+nqalpmbSmpibaMzSBgatfXiG67DnqqPBxdeWoo47iqArfl9bWVjo6OmhubkYSzc3NdHR0eFBhVxb+DNFlzx131LoELsUdVbovra2tXgG6ivAWonMDyXjgXudcabyF6Fx/PHCvcyOGV4iucVRiztisBe719+KcGzKvEF32rLtu9Y7lgXtLtm4174tzFTDkClHS62a2iqQcIfJEN7Ai8HfgoJjtDELUCwPeAvY2syfiJNsTzeyFfva/AiFg8B7A24TIFieY2XVx+6fMbOtE/i5geTMbX1AmgLlmdmDMtxB4DXifEK3iWDP7c8rxp8YyHlKQnt/egJeBb5lZ0flaXQVcfHF6eiVaRw0QuLdaLi52X5zLiHK1EBeY2QRJywM3ESJcjAbWBjYxsyWS1gXeGMQ+TwLWAsab2duSPgZsk1i/qqSxZvZUkcm1F5jZhCL73s7MXoihqv4K9KkQB5Df/ieE+Vf3H+T2LisaIHCvc640ZR1lambvAbcTgv2uBSw2syVx3SIze7mU/cSgw/sDh5rZ23H7Z83sskS2y4DJ8fe+FISDKtGHCa28oboDWGcY27uhOOyw8KkGD9xbssMOO4zDqnVfnKuAsj5DjBXZ9sDxwHzgthi5fjZwsZndV+KuNgSejMF9i7kcmAmcBuxKiFE4JbF+PUn3Aa8SukVvTay7WZKA9YG9SyxTmp0IETlcNXV1Vfd4GQ/cWy1d1b4vzpVZuSrEDeIzPAP+bGbXAcQuyS/Gz2xJe5lZuSKIvgS8LGkfwvPCRJ8Wi4FxZvaipE8Df5K0caKCzXd5bhDLNcfMXh/EsW+OXbjPEbpM+2hpaWkD2mD4M/w755yrvHJ1mS6IQXY3M7Pp+cQYxf46MzsSOIXwbDGVpBti5PoLgMeAcZJWHeC4swgR65fpLo3HfTH+vhdYAHyicGMzWwA8C3xK0sHx+F35IMf92A5oBh4ETkzL0N3d3dHd3T2xu7t74vjx4wfYnXPZ1tnZydy5c7nllls8aK/LrIrNVCNp83zFImkUsAlQdDSmmX05VqrfMbNe4ELgTEkrxn2sJWm/gs2uBH4O3FBw7I9KWi7+Xh/YCHg8pYxrAusBPWZ2Tjz+BDN7ZqDzi8GBDwO+JekjA+V3rlHlg/bmQzB50F6XVZV8D3FN4HxJo+PyXcDZifXzJC2Jvy8zsx8UbH8scDLwkKS3CCNUj09mMLPXgJ8BhEeCH/gCcKKk9wivVxxoZi8l1t8s6X1gBeDHZvZskXOYKinZql0mAqmZLZZ0CXAwYVSsq4ZP9GnsN7R673JPC9rb29vLtGnTOL/OJi/wgL6uPx4PsQo8HqIbjnqvEG+55Zai67bZZpui62rBK8QRq6R4iD5TjXN1rt7/Ee8vaG+9l925JI924bKnrW3pBNuu5jxor2sU3kJ02fPII7UugUvIxyacNm0ab7/9Ns3NzbS3t3vMQpc5XiE654attbX1gwE03k3qssq7TJ1zzjm8QnQuezo7QxSOUaPCt7/v51xZeJepy54JxYKYjACdnctG3+jpWTrAqMbP7CaM5PviGoK/h1gF/h5inajz9/lKMnduenDi0aNh0qS+6Vnizx5d5ZT0HmJVukwlmaSLEsvLS3pe0jVxeaqks+Pv6ZKejnOKPiRp38R2knSspEclPSLpZkkbJ9avIuk8SQskPSjp75K2SCnP1Hj8rsTnU5JGSTpT0gOS5ku6W9J6cZuFkm4t2E+XpAfKf8WcKyKtMuwv3TlXsmp1mb4BjJe0cpwDdEfg6X7yn25mp0naCLhX0uVm9i5hirStgE3NrFfSl4CrYiSLt4ALgCeAjWJQ4vWBtODBALPM7JBkQqx8+wtqPFBQYlcN+8UpbQcbob0RWiC5XOgmLdTcXPPz2y/el4sHe1+cqxPVHFRzHfDV+LukgL5m9ighrNPqMelHhKDBvXH9XwkBiVtjKKctCLEP80GJHzezvwyijAMFNS5HUGI3XIsWhc9I1N4OBS/B09QU0mts0aJFLBqp98U1hGpWiJcC+0haiRD54s6BNpC0OfComT0n6cPAh2LIpqR7gI3jp8vM3i+xPJMLukxXJlR4u8blGZI2K9jmcuAb8feuwNUlHsu58mhthY6O0CKUwndHR80H1DjXCKo2ytTM5knKEVpW1w6Q/fuS9idEtN9pgLwiBCYerD5dpsCiAYIa9xeUeBkeINhVTGurV4DOVUC130O8CjiNgbsaTzezFkL35O8krRSj3b8RnwsmbQ48RAjWu2mMvbiMwQT/LSGocWpQ4kIeINg557Kl2hXib4ATzWx+KZnN7ApCl+i3Y9KphKDBKwNI2gH4PPD72JV6D/ATxeCIkjaStFupwX9LDGqcGpTYVdGWW4aPqytbbrklW/p9KbvOzk5yuRyjRo0il8t54OUKquqL+Wa2CPjlIDc7Efi9pPOBswgDbObHAL//AnaLI1cBvgPMAB6T1Au8CBxZZL+TJX0+sXwQ8GH6D2rcX1BiVy0//WmtS+BS/NTvS9l1dnbS1tZGb5yIoaenh7Y4EYNPnl5+/mJ+FfiL+c7VRtaf38+dO5e3U94xHT16NJMyNBFDHUz4Xj8v5jtXVnvsET6uruyxxx7s4felrNIqw/7S3fD4XKYue158sdYlcClerMP7Ugctk2HJ5XL0pEzE0NzcnPlzq0feQnTOuTrV3t5OU8FEDE1NTbTXwUQMjcgrROecq1Otra10dHTQ3NyMJJqbm+no6PABNRXiXabOOVfHWltbvQKsEm8huoHVW0Da7bcPH1dXtt9+e7b3++IyzF+7qIJMv3ZRGJAWwmTSPn+mcy47SnrtwrtM61k9vEOVFpC2txemTYPzz69NmeqRj/hzLvO8y9T1rx4D0s6bFz6uruy8887svPPOtS6Gc0OWuRaiJAMuNrMpcXl5YDFwp5ntImkqMNHMDpE0HdgfeB5YETjJzC6J2wk4hjBPqhECFh9iZg/G9asQpoHbAXiLOA2cmd0paSHwGvA+8J6ZTazIydZDq6MeA9LmW871cH3cB958882BMzlXx7LYQnwDGJ+f4BvYkVCZFXO6mU0AdgPOk7RCTD8Y2ArY1Mw+AfwUuCrGawS4gBDuaSMz2xiYCoxJ7He7OFl4ZSrDelHHAWmdc66cslghAlwHfDX+LilyvZk9SohfuHpM+hFwqJn1xvV/BW4HWiVtAGwBHGtmUMJwUQAAIABJREFUS+L6x83sL2U9iyzwgLTOuREiqxXipcA+sTW3CXDnQBtI2hx41Myek/Rh4EMxZFTSPcDG8dNlZu8X2Z0Bf5V0r6S2IZ9FVrS2wsKFsGRJ+PbK0DnXgDL3DBHAzOZJyhFah9cOkP37kvYH1gd2GiCvCJXdQD5nZs9IWhP4m6R/mtnfkxlaWlragDbI/oz7dWeXXWpdApdiF78vLuOy2kIEuAo4jYG7S083sxZgMvA7SSuZ2avAG5LWL8i7OfAQ8CCwaQwS3Ec+yLCZPUcIGPzZwjzd3d0d3d3dE7u7uyeOHz9+MOflBnLEEeHTYLIeCPaII47giAa8L27kyHKF+BvgRDObX0pmM7uC0CX67Zh0KnBmfnCOpB2AzwO/j12p9wA/iaNRkbSRpN0kfUjSqjHtQ8CXgAfKeF5uBMoHgu3p6cHMPggEm7VK0bksy2SXKYCZLQJ+OcjNTgR+L+l84CzCAJv5kt4H/gXsZmb5sePfIbx28ZikXuJrF8DHgCtjPbk8oQK9frjnM9IMpxv5jK4uAA6bMKFMpam9tECwvb29TJs2jfMzMgFCV7wvE4rcFw9X5Opd5ipEM1slJW0OMCf+ngnMjL+nF+S7F2hJJP0kftKO8yrhHcY0mw6mzM4NxAPBOld7masQXWMYVmshti4bqcXRCIFgt23A++JGliw/Q3SuYXggWOdqzytE5+qAB4J1rva8y9Rlz95717oEFZH1QLB7N+h9cSOHV4guew46qNYlcCkO8vviMs67TF329PYuG7C4nDo7Q4SPUaPCt78HWLLe3l56K3VfnKsCbyG67PnKV8J3uUczdnZCW9vSyranJyyDz99agq/E++KjTF1WeYXoKqOS87fGF8DLfoy5c/sGPu7thWnToF5fjvfKx7myqUqXqSSTdFFieXlJz0u6Ji5Pjctdkh6Kk3Hn088uss+FkubHbboknRnTJ0m6M6Y9HIME5/dlkrZP7GP3mLZnyv5zkvpMySZppqQn4v7vT+7PZVyxl+D95XjnRoRqtRA/COobp0ZLC+o7K0a5XxN4UNJVJex3OzN7oSDtt8DeZna/pOVYdmaa+YQIGbPj8j7A/YM9GeBIM7tc0nZAB7DREPbR2CrZcsm3DMt9jFwudJMWam72lphzI0A1B9WUFNQ3RpBYADQP8ThrAovjvt43s4cS624FPitpBUmrABsCXUM8DsAdwDrD2N7Vk/Z2KHg5nqamkO6ca3jVfIZ4KXB87CbdhBCtYuvCTDEk0/rAY8CnBtjnzXFiboDfmtnpwOlAt6Q5wPUx/a2Yx4AbgS8DqxFCSK03jHPaCfjTMLZ3QzF1amX2mx84c8wx8OSTMG5cqAx9QE1JplbqvjhXJVWrEEsI6jtZ0ueBt4EDzOylGFGiP326TM3sREmdhLBM34zH2zaR5VLge4QK8XDg6EGfDJwq6eeE1uiktAweILiCKvkPb2urV4BD5BWiy7pqv4fYX1DfWWY2wcy2MLMrkyskLZcYPHPiQAcxswVm9itge0Kg3zUS6+4CxgNjzOyRxDG2SBzjawMc4khCd+uxhGeWfXiA4Ap64YXwcXXlhRde4AW/Ly7Dql0hDiqob158Fjghfo7vL6+kr2pp03Ij4H3g3wXZjqKgZWhmdyaOMeCAHjNbQojHOErSl0s+GTd8e+4ZPq7qOjs7yeVyjBo1ilwut0wA4z333JM9/b64DKvqe4hDDOrbn+QzxHlm9i1gCnB6DOr7HtBqZu8nu1/N7LoS998iaVFi+fvJlWZmkk4GfgjcMNSTcC4LOjs7aWtr+2A2mp6eHtrixAVZnoPVuTyZWa3L0PDOOussO/TQQ2tdjMZRqdcuSj78tjU5bq3NnTs3NWDx6NGjmTRpEl1xwoQJEyZUu2gV57PvZN6AA1LA5zJ1zpUorTLsL925rPGp25wbpJHaWsjlcvSkTFzQ3NzMnDlzPmg5j9Tr47LPK0SXPd/9bq1LMCK1t7cv8wwRoKmpifY4ccF3/b64jPMK0WXP5Mm1LsGIlB84c8wxx/Dkk08ybtw42tvbP0if7PfFZZxXiC57nnoqfI8dW9tyjECtra1FR5Q+Fe/LWL8vLqN8UI3LnilTwqdaPGhwSaZMmcKUat4X58rMW4jO9ceDBjs3YniF6CqrEu/sVSpAcJosBg2uleR98ZGmLoO8y9S5/njQYOdGjIq1ECUZcLGZTYnLyxPiFN5pZrtImgqcSggUvCJwupmdH9MnmtkhKftcCLxGmJ8U4O9m9j1JkwhTwo2On1lmNj3u63+AHcxsdtzH7sAVwF5mdnnB/nPAE8DJZnZcTBsTy31eDGA8HdgfeD6W+yQzS43t6KhMS6GaM9V40ODS1XgGIeeGq5Jdpm8A4yWtbGZvAjsSKr+kWbGSWRN4UNKAk2qTEvKJEHFibzO7X9JyQEti3XxCCKjZcXkf4P5+9v84sAtwXFzeC3iwIM/pZnaapI2AeyVdbmbvllB2Vw6HH169Y7W3L/sMETxocBGHV/O+OFcBlX6GeB3wVeByQqV0CSlBgc3sOUkLgOYhHmdNQisOM3sfeCix7lZga0krEFqPGwJd/ezrTeBhSRPN7B5gMnAZsHZKuR+Nk4ivDjw3xLK7wdp11+ody4MGl2zXat4X5yqg0s8QLwX2kbQSsAlwZ1omSesD6wOPlbDPmxNxC/PRJ04HuiVdKemAeLw8A24EvgzsRojJWGq51yV0zz5TpNybA4+amVeG1dTdHT7V0toKCxfCkiXh2yvDVN3d3XRX8744V2YVbSGa2bz4XG5f4NqULJMlfR54GzjAzF5Khmkqok+XqZmdKKkT+BLwzXi8bRNZLgW+B6wGHE5BLMQU1wMnAc8Cs1LWf1/S/oRKfKe0HbS0tLQBbTByoyNUzAEHhG9/VlVXDoj3xecydVlVjVGmVwGnEbpLC82KAXm3MLMrkyskLZdoCZ440EHMbIGZ/QrYHthU0hqJdXcB44ExZvZI4hhbJI7xtUT+d4B7CZXnH1MOd7qZtRC6U39X0CIFoLu7u6O7u3tid3f3xPHjxw9UfOcaUn8BhZ2rN9V4D/E3wCtmNl/StqVuFJ8FlhRYTdJXgWstBHfciNDN+e+CbEcBbxUc487kMWJrNm8GcIuZvVis1WpmV0j6NvBt4LxSyurcSOEBhV3WVLxCNLNFhFciyuVmSfnXLuaZ2beAKcDpcYDLe0Crmb2frMjM7LrBHMTMHqTv6NI0JwK/l3S+mS0ZzDHcyDESus3zAYLz55oWULi3t5dp06ZxfsYnNfBu4cZUsQrRzFZJSZsDzIm/ZwIzU/Kkpsd1uSLp+xRJL3aMqUXyLyR0rRbdj5lNL1h3L8u+5uGcwwMKu+zxqdtc9hx7bK1LMGgjoUVx4403ArDDDjsAAwcUdq7eeIXosif+g+vqyw4F92WggMLO1Rufy9RlT1fX0omkXd3o6ur64DkihIEzHR0dNDc3I4nm5mY6Ojp8QI2rWwoDM10lnXXWWXbooYfWuhiNw+fMrEv5wTTeHerq0IAvuIO3EJ1zzjnAK0TnBqezM0TAGDUqfPuL5s41DB9U41ypOjuXjXzR0xOWwec3da4BeIXoamM4L6onI7NX09y5fQMD9/bCtGmQ5RfN/Zmfc0CNu0wlmaSLEsvLS3pe0jVxeWpc7pL0UJxQO59+dpF9LpQ0PzFH6ZkxfZKkO2PawzHQb35fJmn7xD52j2l7puw/J+nNRJl+Lcm7nqtpvfXCp9qKvVDuL5oDcMopp3DKKafUuhjODVmtW4hZDSK8wMwmSFoeuAn4OnBFCeVyeVlsleRyoZu0UHNzNs+nzLbaaqtaF8G5YamHlk0+iDAsDSLcR4w5WLYgwmZWGET4s5JWkLQKAwcRzpfpPeD2mN9Vy+23h0+1tbdDU9OyaU1NId1x++23c3st7otzZVIPFWJWgwgjqYkQbmp+KfldmRx9dPhUW2srdHSEFqEUvjs6fEBNdPTRR3N0Le6Lc2VS6y7TrAYR3kBSF6Ei/XNaJA0PENygWlu9AnSuQdVDCxGyF0R4QSzTZoXRL/I8QLAbqTwosMuqmrcQo6wGEXbOJTz77LMeFNhlVl1UiFkNIuzccDVSd3pXVxevvvoqhfMjN0pQ4MHyOV2zp6YVYiMFEXZVdMYZtS6BS7Hhhhty7733pq7zoMAuC+qihejcoEwoqZc8ExqtFeFBgV2W1cugGudKd+ON4ePqyo033si+++5LU8G7mh4U2GWFV4gue04+OXxcXTn55JO54447PCiwyyzvMnXOlVVra6tXgC6TvIXonHPO4RWia1QeyNc5N0jeZeoajwfydc4NgVeIrnyq9ZJ5vqIrdrxGDeSbVIevMJx33nm1LoJzw+IVosuewhBMhTyQb020tLQMnMm5OlaVClGSAReb2ZS4vDwhNuGdZraLpKnAqYTgwCsB55nZ6THvdOB1MzstZb8fB84APkOIhrEQOMzMHpG0ESHk0ycJc5a+CpxgZn+P2+4MnAR8CBBwjZkdEY93ArCRmT0W834f+AXwGTO7R9L1wFqE63crcHCcV3Vkq1ar5eqrw/euu6av90C+NXF1vC+7FrsvztW5ag2qeQMYL2nluLwjofJLmmVmE4DPAcdIGtvfDhUmIb0SmGNmG5jZpwghmz4WYx3+BeiI6z4NHEqIp4ik8cDZwH5m9knCVGyPJ3Y/H0hO9bYnkAwovLeZbRq3+yiwVykXwZXJjBnhU4wH8q2JGTNmMKO/++JcnavmKNPrgK/G3/uSHuoJM3uREAR4rQH2tx3wrpn9OrFtl5ndCrQCd5jZVYl1D8R5SwF+CLSb2T/juvfM7NzEvv9ECBScD0z8CvB8Yl+vxp/LAysS4iK6euGBfJ1zQ1DNCvFSYJ/YetsEuDMtk6RxhG7TeQPsbzyQPpMwbAz8Y4jbQuhefSq2JPcFZqWU8wbgOeA14PIByuqqrbUVFi6EJUvCt1eGzrkBVK1CNLN5QI5QwVybkmWypAcJXZe/NLO3UvIMiaQrJT0g6YpBbHYpodv064Su2WWY2ZcJrdjRwBcL17e0tLS1tLTc09LScs8DDzwwxJI755yrlmq/mH8VcBrp3aWzzGxjYGtgRhww8wFJYxOR6w8EHgQ+XeQ4DwKb5xfMbHdgKvCRxPpi2+ZdTYih+GSii3QZsdK+iti9mtTd3d3R3d09sbu7e+L48R4tyrmBdHZ2ksvlGDVqFLlcjk6fTMFVWbVfu/gN8IqZzZe0bVoGM7tD0kXAfxEi2OfTn2LZyPUCTpG0v5mdH9M+AzQBvweOkvS1xHPE5CiLU4ErJN0WR6SOIoxO/UXieG9K+hHwSLJ8klYBVjWzxXG07FcII01dtVx0Ua1L4FJcNIz70tnZSVtbG73xHdOenh7a4mQKPi+qq5aqVohmtgj4ZQlZfwb8Q9Ip/ezLJO0OnCHpx8BbLH3t4k1JuwC/kHQG8CzhWd/Jcdt5kg4DLpHURBgU85eUY1yacugPAVdJGg0sB9wE/Doln6uUsf0OQM6sbas1sUEdmjt3bp8gwr29vUybNo3zMzKZgsd7zL6qVIhmtkpK2hxgTvw9k0TUejN7Bsh3mU7vZ7/PAHsXWfdPQuut2LbXANekpKcez8y2TSx+pth+XRXMimOcJk+ubTncMp577jkA1lxzzUFvW1gZDpTuXCX4TDUue371q/DdYBVi1lsY+RbuUM4jl8vRkzKZQnNzc+avi8sOj3bhnKu59vZ2mgomU2hqaqLdJ1NwVeQVonOu5lpbW+no6KC5uRlJNDc309HR4QNqXFV5l6lzri60trZ6BehqyluIzjnnHN5CdFnT2QkLFsDTT4eoFu3tPi1bnbj8cp/B0GWbV4guOzo7Q+T7fIDgnp6wDF4p1oExY8bUugjODYtXiK48qvFS+dy5fYP89vbCtGmQkZe3Gzke48yZMwGYOnVqTcvh3FBV7RmiJJM0I7GcD8Y7mH1MlfR8Yk7TLkmfkjRK0plxAu/5ku6WtF7cZqGkWwv20yXpgfh7DUk3S3pd0tn9HHvbeA7TEmmbxbQjBnMeboiKvaTtL2/XhZkzZ35QKTqXRdVsIb4NfEPST83shcFuHOcNhTAJ+CEF6/YF1gY2MbMlktYlBCXOW1XSWDN7StInC3b9FnAcISTUQLNwzwcmAxfG5X2A+wd7Lg2pGi2fXC50kxZqbm7olpdzrjqqOcr0PaAD+H7hCknNkmZLmhe/x8X0mZJ+IelmwvymxawFLDazJRDmTDWzlxPrLyNUZFAQnNjM3jCz2wgV40CeBFaS9LE4ufhOhMDHrhra26Hg5W2amkK6c84NU7VfuzgHaJW0WkH62cDvzGwToBM4M7HuE8AOZnZ4XJ5c0GW6MqHC2zUuz5C0WcH+Lwe+EX/vSgjtNFSXA3sBWxGCEHt/XbW0tkJHB4weHZabm8OyD6hxzpVBtaNdvCrpd8D3gDcTq7ZkaYV1EfDzxLo/mNn7ieU+XabAIkkthEC9XwRmS9rLzGbH9S8BL0vaB3gY6B3GaVwGzAL+g9DS3CotU0tLSxvQBiM7ikHZtbYuHUDj3aTOuTKqxYv5ZwDTCGGUirHE7zeK5kpuYPa2mV1nZkcCpxAi3SfNIrRQ04IT9yFp90QrdGLiOP8C3gV2BGYX294DBFfQtdeGjxtQNYPuXnvttVzr98VlWNVfuzCzlyRdRqgUfxOTbycMULkIaAVuG8w+JW0O/MvMnonBfjcB5hVku5LwrPEGwgCcgcp5Zdwmf4xtE6uPB9Y0s/fDo0RXVYXPEV2qagfdLZyc27msqdV7iDOAZLfn94DfSDoSeB74v/1sO1nS5xPLBwEfBs6PQXsB7iI8l/yAmb1GHJhTWIlJWhj3saKkrwNfMrOHihXAzG7vp3yu0s49N3wfdFBVD5u1ru9qB919+umnAVhnnXXKvu9SeJgoN1xVqxCTQYLN7FmgKbG8kPDsr3CbqQXLM0kEEi5wfZHj5lLSFpJ4xSItT8o2c4gBjQvSpw+0rSuzyy4L31WuELOm2kF3n3/+eaB2FaJzw+Uz1ThXoqy1QKoddHc4AYKdqwce7cK5BuVBd50bHK8QnWtQHnTXucHxLlPnGpgH3XWudDKzgXO5YWlpabkAWFTrckSfBu6tdSHqlF+bdH5divNrk67ersvC7u7umQNl8gpxhGlpabmnu7t74sA5Rx6/Nun8uhTn1yZdVq+LP0N0zjnn8ArROeecA7xCHIk6al2AOubXJp1fl+L82qTL5HXxZ4jOOecc3kJ0zjnnAK8QG56kj0j6m6RH4/fqRfKNk/RXSQ9LekhSrrolrb5Sr03M+2FJT0s6u1ieRlHKdZE0QdIdkh6UNE/S5FqUtVok7SSpW9Jjkn6csn60pFlx/Z0j4f8fKOm6/CD+ezJP0mxJzbUoZ6m8Qmx8PwZmm9lGhPiNff6jjX4HnGpmnwQ+CzxXpfLVUqnXBuAk4JaqlKr2SrkuvcC3zGxjYCfgDEn/p4plrBpJyxFiqe4MfArYV9KnCrJNA142sw2B04mRdRpZidflPmCimW0CXM6ywd/rjleIjW834Lfx92/pGziZ+B/x8mb2NwAze93MeqtXxJoZ8NoASPo08DHgr1UqV60NeF3M7BEzezT+fobwB9RHq1bC6vos8JiZPW5m7wCXEq5RUvKaXQ5sr8YPljrgdTGzmxP/lswF1q1yGQfFK8TG9zEzWwwQv9dMyfMJ4N+SrpB0n6RT419/jW7AaxMDTs8Ajqxy2WqplP9mPiDps8CKwIIqlK0W1gGeSiwvimmpeczsPeAVYI2qlK52SrkuSdOA6ypaomHyuUwbgKQbgY+nrDqmxF0sD2wNbAY8CcwCpgIXlqN8tVSGa3MQcK2ZPdVIf/CX4brk97MWcBHwbTNbUo6y1aG0G184PL+UPI2m5HOWtB8wEdimoiUaJq8QG4CZ7VBsnaRnJa1lZovjP15pzwYXAfeZ2eNxmz8Bk2iACrEM12ZLYGtJBwGrACtKet3M+nveWPfKcF2Q9GHgL8CxZja3QkWtB4uAsYnldYFniuRZJGl5YDXgpeoUr2ZKuS5I2oHwh9Y2ZlaZ6NRl4l2mje8q4Nvx97eBP6fkuRtYXVL+GdAXgYeqULZaG/DamFmrmY0zsxxwBPC7rFeGJRjwukhaEbiScD3+UMWy1cLdwEaS1ovnvQ/hGiUlr9mewE3W+C95D3hdJG0GnAd8zczqfqCeV4iN77+BHSU9CuwYl5E0UdIFAGb2PuEf+9mS5hO6Qs6vUXmracBrM0KVcl32Br4ATJXUFT8TalPcyorPBA8BbgAeBi4zswclnSjpazHbhcAakh4DfkD/I5YbQonX5VRCz8of4n8jhX9I1BWfqcY555zDW4jOOecc4BWic845B3iF6JxzzgFeITrnnHOAV4jOOecc4BWic1UlaaakhUPYboKk6ZI+krLOJE0vR/nqmaSPS7pK0kvxnA+L6btKmi/prZj+fyTNkTRnkPufLqmiw+7jMb5YyWO4ofOZapzLhgnACcDF9J0BZUvCrCGN7njC1F9TgcXAwjgrTCdwO3Aw8A7wGmHKvcG6ALi+LCUt7gSgHbipwsdxQ+AVonMDiBOdK76IXLhudK2no2rwadOSPgncb2ZX5hNifL1VCS+F/z2Rd9AzLZnZIkbGHxauCO8ydQ1D0qaSrpT0oqQ3Y+DSoxLrJen7Mf0dSYslnR3n5EzuxyS1S/qxpCcIrY7/lLRtXPcNSedLeh54tuD4V0l6OR7/fyVtXUK5fyLpH5JekfSCpJskTUqsnwr8T1x8NJbBFIPQpnWZKgRuvSOW4xVJf5LUUpBnjqTbJO0Qj98r6QFJqWGwUsq9nqSLJP1L0tuSHpf0y4I8+0m6P3ZnvhDzr5Wyr/0L8l2Y7x6WlItdmdsS5pXNn/9MYGHcxYUxbU7i3OYUHOOjks6V9FQs71OxPKPj+j5dppKWl3SUpH/GbZ6RNEPSSok8uXjsAxRmaVks6d+Srpa0biJfft/HJM5heinX2lWHtxBdQ1AIQTQHeAz4PuEv/Y2ATRLZ2oGjCEFNryYENT0J2FTSNgXRGqYCjxOmtHuDMGnxanHdWYQwNlOAleLxNwduJQRE3Z8QQPdA4EZJW5nZvf0Ufx1CUNlFwIeA/YC/S5poZvMIE2ifDBwL7MXSVsziItdip7jNTcBkwtRZJwK3SZpgZk8nsm8A/BL4KfACcDhwuaT/MLPHihVY0nrAXfE8TwAeJUz0/KVEnjbCPJazCNd9beAUYAtJm5vZ6zHff8fjnkkIs7VOPN/xkraK57ll3Nf7LO0OXQxcA/wh5v8L8GqR8q5O6Fb9SMw7jxDWajdC6KpirfyLgV0JAX9vJ7RSTwJywB4FeY+Kef5f3PcMQnduPsLDlsAdwMx4LuAt0vpiZv7xT+Y/wN8Jsdmaiqz/CPAWMLMgfT9CyJqvJdKMUAGuXJB327juypT9zybM57hiIm25mPanRNpMYGE/57Ec4Q/VbuCXifSp8dgbpmxjwPTE8j2ECmr5RNp6wLvALxJpc2LaRom0NQmVztEDXO/fAa8Da/dzHs8CNxekfz6W93txORePd3xBvs/FfF9PpN0GzCnIt2HMN7UgfU4yL+EPgveBzfo5p+nhn8QPlreO+/5WQb7WmD4hcQ4G3FKQ74iYvnYizYCTa/3/i3/SP95l6jJPUhPhH9BOWxqdu9AkYDThL/6kS4H36Bun7Xoze7PIvq5MLkhaOW7/B2BJ7GZbnjBJ+o2ESbD7K/8Okm6W9GIsy7uEoM0t/W1XZF8fAjYHZlnimaeZPQH8L33P81GLke9jvucI4Z7GDXCoLwHXmFmfcD9RC6Fy7UwmmtltQE+iHDsSHt105q9bvHZ3Elp7/V67QfgScLeZ3TeIbXYidJf/saBsf43rC8v2l4Ll+fF7oGvp6oR3mbpGsDrhH9X+up/yryss081oZu/FiqjwdYbU7sgi6z5CaBEdFz99SBplKQF0Y1frtYSIAdPivt8njHhcqTB/CVYnVMRp5f8X0FyQlhaz7+0Sjr0GQ7jeiXLk168Zv4t1z5Yr6vwawP2D3GZNQnfq6/3sM6nwWua7YYdyH10NeIXoGsHLwBLCs6di8v9YfRx4MJ8Y/+JfA3ixIH9/76MVrvt3PP45hK7EvhsUjya/B6FV+A0zezdRrtXjfgfr5Vi+j6es+zh9z3OoXqD0651Wjnvi73x5vkQoe6FqlTfNi4Ru9mIDo4q1jl1GeZepy7zYTXobsF/svkwzl/AX+z4F6ZMJfxjeMozjv0EYULMp8A8zu6fw08/mTYQW4QeVrMKL24XdbPnWRrHzS5blXmAvhddF8vtsBrZiGOdZ4K/ALmkjRqNuwjPEZa53HCTTnCjH3wh/TIxLu26xq7dc5f2spE0Hsc31hNbdakXKNpQK8R0GuIeudryF6BrFEYR/ZO+QNIPQnbc+YeDDoWb2kqRfAEdJeoPQTflJwojD2+j7/GewfkAY2HODpAsJXYVjCM/zljOzYgFjrwcOA2ZK+h/Cs8PjgKcL8uXfqztY0m8Jzxnnmdk7Kfs8jnA+10g6lzDK9CfAK4SRj+VwAvBV4HZJpxC6PNcBdjKz/czsfUnHA+dJupjw7HYdwkjfR4mvkZjZAkk/A86Or4XcQmiVjSU8X7zAzG4uQ3lPB75JGPV7MuH53hjCKNMDzey1wg3MbI6kSwijbn9BGFW7hDCI5ivAj8zskUGW4yHgq5KuJ7SInxlixeoqodajevzjn3J9gM0Ir1P8G3gT+CfhH638ehFeyegm/KW+mNDN+eGC/aSOBGTpKNMdihz/k4RBOs8RWnSLgKuAryTyzKRglClwKPBELPPdwA4UjJKM+U4gVJT5FmUuUd7pBXl3Igzxf5NQEf4ZaCnIMwe4LeU8FlIwGrfI+W4AXELojnyb8JrK6QV59iM8u3ub0AV5EbBWyr6mEFrxbxCe2T0MnA1LwHCRAAAAj0lEQVSsm8gz5FGmMW1NoCPe93cIo5J/C4yO66eTGGUa00YB/xXP4a14Le8Hfk5oOcLSUabfKfLfy7aJtM8RWvBvpd03/9T2o3iTnHPOuRHNnyE655xzeIXonHPOAV4hOuecc4BXiM455xzgFaJzzjkHeIXonHPOAV4hOuecc4BXiM455xzgFaJzzjkHwP8HU5Gl/TZI4voAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7ff9a0354048>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure(figsize=(7,9))\n",
    "plt.subplots_adjust(left=0.2)\n",
    "\n",
    "if models[0][0] == 'A':\n",
    "    models=models[::-1]\n",
    "\n",
    "values1=np.array([np.average(modelOutput[ii].piControl[1]) for ii in models])\n",
    "values2=np.array([np.max(modelOutput[ii].piControl[1]) for ii in models])\n",
    "values3=np.array([np.min(modelOutput[ii].piControl[1]) for ii in models])\n",
    "errorbar(values1,np.arange(len(values1))+0.25,xerr=[values1-values3,values2-values1],fmt='o',color='r',label='$R$ vs. $T_A$')\n",
    "\n",
    "values1=np.array([np.average(modelOutput[ii].piControl[0]) for ii in models])\n",
    "values2=np.array([np.max(modelOutput[ii].piControl[0]) for ii in models])\n",
    "values3=np.array([np.min(modelOutput[ii].piControl[0]) for ii in models])\n",
    "errorbar(values1,np.arange(len(values1))-0.25,xerr=[values1-values3,values2-values1],fmt='o',color='k',label='$R$ vs. $T_S$')\n",
    "\n",
    "xlim(-0.7,0.3)\n",
    "xlabel('correlation coefficient',fontsize=16)\n",
    "axvline(ceresR,color='k',ls='dashed') # from Fig. 1 file\n",
    "axvline(ceresRta,color='r',ls='dashed') # from Fig. 1 file\n",
    "\n",
    "# axvline(np.average(values1),color='b',ls='dashed')\n",
    "gca().set_yticks(range(len(models)));\n",
    "gca().set_yticklabels(models);\n",
    "ylim(-0.5,len(models)-0.5)\n",
    "colorplotLeft(gca())\n",
    "legend(fontsize=12)\n",
    "savefig('fig4.eps')\n",
    "# savefig('fig4')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### older code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "if False:\n",
    "    figure(figsize=(9,7))\n",
    "    plt.subplots_adjust(bottom=0.2)\n",
    "    values1=np.array([np.average(modelOutput[ii].piControl[0]) for ii in models])\n",
    "    values2=np.array([np.max(modelOutput[ii].piControl[0]) for ii in models])\n",
    "    values3=np.array([np.min(modelOutput[ii].piControl[0]) for ii in models])\n",
    "    errorbar(np.arange(len(values1)),values1,yerr=[values1-values3,values2-values1],fmt='o',color='k')\n",
    "    ylim(np.average(values1)-0.4,np.average(values1)+0.4)#;grid()\n",
    "    ylabel('correlation coefficient',fontsize=16)\n",
    "    axhline(ceresR,color='r',ls='dashed') # from Fig. 1 file\n",
    "    axhline(np.average(values1),color='b',ls='dashed')\n",
    "    gca().set_xticks(range(len(models)));\n",
    "    gca().set_xticklabels(models,rotation='90');\n",
    "    xlim(-0.5,len(models)-0.5)\n",
    "    colorplotLeft(gca())\n",
    "    # savefig('fig4.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "if False:\n",
    "    figure(figsize=(7,9))\n",
    "    plt.subplots_adjust(left=0.2)\n",
    "    values1=np.array([np.average(modelOutput[ii].piControl[0]) for ii in models])\n",
    "    values2=np.array([np.max(modelOutput[ii].piControl[0]) for ii in models])\n",
    "    values3=np.array([np.min(modelOutput[ii].piControl[0]) for ii in models])\n",
    "    errorbar(values1,np.arange(len(values1)),xerr=[values1-values3,values2-values1],fmt='o',color='k')\n",
    "    xlim(np.average(values1)-0.4,np.average(values1)+0.4)#;grid()\n",
    "    xlabel('correlation coefficient',fontsize=16)\n",
    "    axvline(ceresR,color='r',ls='dashed') # from Fig. 1 file\n",
    "    axvline(np.average(values1),color='b',ls='dashed')\n",
    "    gca().set_yticks(range(len(models)));\n",
    "    gca().set_yticklabels(models);\n",
    "    ylim(-0.5,len(models)-0.5)\n",
    "    colorplotLeft(gca())\n",
    "    # savefig('fig4.eps')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}