model_evaluation_notebook.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model Evaluation Tutorial "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first part of this tutorial describes a python algorithm/script that is used to compare observed and modelled datasets. The user provides two arrays (i.e. obs and sim) and an estimate of the uncertainty around the observed data (constant) to the comparison function. When called by the user the function performs the comparison and returns a summary statistic. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following block of code creates the sim_vs_obs comparison functions. The first function is named \"compare\" and the second \"compare_daily\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare(sim, obs, unctype, unc) :\n",
    "\t\t\n",
    "\timport pandas as pd\n",
    "\timport numpy as np\n",
    "\n",
    "\tsim = np.array(sim)\n",
    "\tobs = np.array(obs)\n",
    "\t\n",
    "\tfit = pd.DataFrame(columns=['oneorzero','obsday'])\n",
    "\n",
    "\tif (unctype == 'pc') : \n",
    "\n",
    "\t\tfor ct in range(len(sim)):\n",
    "\n",
    "\t\t\tif ( np.isnan(obs[ct]) ) : continue \n",
    "\n",
    "\t\t\tif ( (obs[ct] >= 0) and (obs[ct]*(1-unc) <= sim[ct] <= obs[ct]*(1+unc)) ) :\n",
    "\t\t\t\t fit = fit.append({'oneorzero':1,'obsday':ct},ignore_index=True)\n",
    "\t\t\t\t continue\n",
    "\n",
    "\t\t\tif ( (obs[ct] < 0) and (obs[ct]*(1-unc) >= sim[ct] >= obs[ct]*(1+unc)) ) :\n",
    "\t\t\t\t fit = fit.append({'oneorzero':1,'obsday':ct},ignore_index=True)\n",
    "\t\t\t\t continue\n",
    "\n",
    "\t\t\tif (ct not in [0,1,2,3,len(sim)-3,len(sim)-2,len(sim)-1,len(sim)]) and (obs[ct] >= 0) and \\\n",
    "\t\t\t\t (      ( obs[ct]*(1-unc) <= sim[ct+1] <= obs[ct]*(1+unc) )  \n",
    "\t\t\t\t\t or ( obs[ct]*(1-unc) <= sim[ct+2] <= obs[ct]*(1+unc) )  \n",
    "\t\t\t\t\t or ( obs[ct]*(1-unc) <= sim[ct+3] <= obs[ct]*(1+unc) )  \n",
    "\t\t\t\t\t or ( obs[ct]*(1-unc) <= sim[ct-1] <= obs[ct]*(1+unc) )\n",
    "\t\t\t\t\t or ( obs[ct]*(1-unc) <= sim[ct-2] <= obs[ct]*(1+unc) ) \n",
    "\t\t\t\t\t or ( obs[ct]*(1-unc) <= sim[ct-3] <= obs[ct]*(1+unc) ) ):\n",
    "\t\t\t\t fit = fit.append({'oneorzero':1,'obsday':ct},ignore_index=True)\n",
    "\t\t\t\t continue  \n",
    "\n",
    "\t\t\tif (ct not in [0,1,2,3,len(sim)-3,len(sim)-2,len(sim)-1,len(sim)]) and (obs[ct] < 0) and \\\n",
    "\t\t\t\t (      ( obs[ct]*(1-unc) >= sim[ct+1] >= obs[ct]*(1+unc) )  \n",
    "\t\t\t\t\t or ( obs[ct]*(1-unc) >= sim[ct+2] >= obs[ct]*(1+unc) )  \n",
    "\t\t\t\t\t or ( obs[ct]*(1-unc) >= sim[ct+3] >= obs[ct]*(1+unc) )  \n",
    "\t\t\t\t\t or ( obs[ct]*(1-unc) >= sim[ct-1] >= obs[ct]*(1+unc) )\n",
    "\t\t\t\t\t or ( obs[ct]*(1-unc) >= sim[ct-2] >= obs[ct]*(1+unc) ) \n",
    "\t\t\t\t\t or ( obs[ct]*(1-unc) >= sim[ct-3] >= obs[ct]*(1+unc) ) ):\n",
    "\t\t\t\t fit = fit.append({'oneorzero':1,'obsday':ct},ignore_index=True)\n",
    "\t\t\t\t continue \n",
    "\n",
    "\t\t\tfit = fit.append({'oneorzero':0,'obsday':ct},ignore_index=True)\n",
    "\n",
    "\t\n",
    "\tif (unctype == 'qt') : \n",
    "\n",
    "\t\tfor ct in range(len(sim)):\n",
    "\n",
    "\t\t\tif ( np.isnan(obs[ct]) ) : continue \n",
    "\n",
    "\t\t\tif ( (obs[ct] >= 0) and (obs[ct]-unc <= sim[ct] <= obs[ct]+unc) ) :\n",
    "\t\t\t\t fit = fit.append({'oneorzero':1,'obsday':ct},ignore_index=True)\n",
    "\t\t\t\t continue\n",
    "\n",
    "\t\t\tif ( (obs[ct] < 0) and (obs[ct]+unc >= sim[ct] >= obs[ct]-unc) ) :\n",
    "\t\t\t\t fit = fit.append({'oneorzero':1,'obsday':ct},ignore_index=True)\n",
    "\t\t\t\t continue\n",
    "\n",
    "\t\t\tif (ct not in [0,1,2,3,len(sim)-3,len(sim)-2,len(sim)-1,len(sim)]) and (obs[ct] >= 0) and \\\n",
    "\t\t\t\t (      ( obs[ct]-unc <= sim[ct+1] <= obs[ct]+unc )  \n",
    "\t\t\t\t\t or ( obs[ct]-unc <= sim[ct+2] <= obs[ct]+unc )  \n",
    "\t\t\t\t\t or ( obs[ct]-unc <= sim[ct+3] <= obs[ct]+unc )  \n",
    "\t\t\t\t\t or ( obs[ct]-unc <= sim[ct-1] <= obs[ct]+unc )\n",
    "\t\t\t\t\t or ( obs[ct]-unc <= sim[ct-2] <= obs[ct]+unc ) \n",
    "\t\t\t\t\t or ( obs[ct]-unc <= sim[ct-3] <= obs[ct]+unc ) ):\n",
    "\t\t\t\t fit = fit.append({'oneorzero':1,'obsday':ct},ignore_index=True)\n",
    "\t\t\t\t continue  \n",
    "\n",
    "\t\t\tif (ct not in [0,1,2,3,len(sim)-3,len(sim)-2,len(sim)-1,len(sim)]) and (obs[ct] < 0) and \\\n",
    "\t\t\t\t (      ( obs[ct]+unc >=  sim[ct+1] >= obs[ct]-unc )  \n",
    "\t\t\t\t\t or ( obs[ct]+unc >=  sim[ct+2] >= obs[ct]-unc )  \n",
    "\t\t\t\t\t or ( obs[ct]+unc >=  sim[ct+3] >= obs[ct]-unc )  \n",
    "\t\t\t\t\t or ( obs[ct]+unc >=  sim[ct-1] >= obs[ct]-unc )\n",
    "\t\t\t\t\t or ( obs[ct]+unc >=  sim[ct-2] >= obs[ct]-unc ) \n",
    "\t\t\t\t\t or ( obs[ct]+unc >=  sim[ct-3] >= obs[ct]-unc ) ):\n",
    "\t\t\t\t fit = fit.append({'oneorzero':1,'obsday':ct},ignore_index=True)\n",
    "\t\t\t\t continue \n",
    "\n",
    "\t\t\tfit = fit.append({'oneorzero':0,'obsday':ct},ignore_index=True)\n",
    "\n",
    "\n",
    "\treturn round(float(sum(fit.oneorzero))/float(len(fit)), 4)\n",
    "\n",
    "\n",
    "# same-day-only comparison\n",
    "# sim/obs : simulated and observed data\n",
    "# unctype/unc : type of uncertainty info values (amount or %) and uncertainty info\n",
    "def compare_daily(sim, obs, unctype, unc) :\n",
    "\t\n",
    "\timport pandas as pd\n",
    "\timport numpy as np\n",
    "\n",
    "\tfit = pd.DataFrame(columns=['oneorzero','obsday'])\n",
    "\n",
    "\tif (unctype == 'pc' ) : \t\n",
    "\n",
    "\t\tfor ct in range(len(sim)):\n",
    "\n",
    "\t\t\tif ( np.isnan(obs[ct]) ) : continue \n",
    "\n",
    "\t\t\tif ( (obs[ct] >= 0) and (obs[ct]*(1-unc) <= sim[ct] <= obs[ct]*(1+unc)) ) :\n",
    "\t\t\t\t fit = fit.append({'oneorzero':1,'obsday':ct},ignore_index=True)\n",
    "\t\t\t\t continue\n",
    "\n",
    "\t\t\tif ( (obs[ct] < 0) and (obs[ct]*(1-unc) >= sim[ct] >= obs[ct]*(1+unc)) ) :\n",
    "\t\t\t\t fit = fit.append({'oneorzero':1,'obsday':ct},ignore_index=True)\n",
    "\t\t\t\t continue\n",
    "\n",
    "\t\t\tfit = fit.append({'oneorzero':0,'obsday':ct},ignore_index=True)\n",
    "\n",
    "\tif (unctype == 'qt' ) : \t\n",
    "\n",
    "\t\tfor ct in range(len(sim)):\n",
    "\n",
    "\t\t\tif ( np.isnan(obs[ct]) ) : continue \n",
    "\n",
    "\t\t\tif ( (obs[ct] >= 0) and (obs[ct]-unc <= sim[ct] <= obs[ct]+unc) ) :\n",
    "\t\t\t\t fit = fit.append({'oneorzero':1,'obsday':ct},ignore_index=True)\n",
    "\t\t\t\t continue\n",
    "\n",
    "\t\t\tif ( (obs[ct] < 0) and (obs[ct]+unc >= sim[ct] >= obs[ct]-unc) ) :\n",
    "\t\t\t\t fit = fit.append({'oneorzero':1,'obsday':ct},ignore_index=True)\n",
    "\t\t\t\t continue\n",
    "\n",
    "\t\t\tfit = fit.append({'oneorzero':0,'obsday':ct},ignore_index=True)\n",
    "\n",
    "\treturn round(float(sum(fit.oneorzero))/float(len(fit)), 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create 2 arrays of example modelled and observed Net Ecosystem Exchange (NEE) data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "\n",
    "sim = np.array([3.20077443e+00,  3.89910202e+00,  3.89860972e+00,  3.59352206e+00,\n",
    "        3.50346302e+00,  2.85976351e+00,  2.27231692e+00,  1.87589841e+00,\n",
    "        1.47556985e+00,  1.70785929e+00,  8.45365624e-01,  9.08928307e-01,\n",
    "        2.20131278e-01,  4.71314923e-02, -6.64671327e-01, -7.67143037e-01,\n",
    "       -1.15083273e+00, -1.54285147e+00, -2.08870210e+00, -2.04016013e+00,\n",
    "       -2.14789965e+00, -2.71875385e+00, -9.00387997e-01, -3.59870503e-01,\n",
    "       -9.65612313e-01, -1.47855614e+00, -1.21594273e+00, -1.92575484e+00,\n",
    "       -1.39136095e+00, -1.41483321e+00, -1.15581224e+00, -6.62026835e-01,\n",
    "        7.08218222e-02,  4.21816310e-02, -2.24205982e-01,  4.01284725e-02,\n",
    "       -3.14518533e-02,  1.78466466e-01,  3.58637157e-03,  4.23179849e-01,\n",
    "        5.94252198e-01,  4.27781032e-01,  5.76851040e-01,  8.33181359e-01,\n",
    "        1.07823158e+00,  1.28288443e+00,  1.32597672e+00,  1.36958171e+00,\n",
    "        1.41083922e+00,  1.25352988e+00])\n",
    "\n",
    "obs = np.array([1.31088054e+00,  1.98014745e+00,  2.57458262e+00,  8.87110118e-01,\n",
    "       -5.34521438e+00, -2.83797040e+00, -5.19935729e-01,  7.09602213e-01,\n",
    "        3.06549885e-01, -8.91126422e-01, -1.88046196e+00, -1.91680980e+00,\n",
    "       -1.57424605e+00, -1.81943746e+00, -3.71028183e+00, -1.74877815e+00,\n",
    "       -6.00830765e+00, -5.53215022e+00, -6.47411172e+00, -6.23745443e+00,\n",
    "       -6.00962183e+00, -6.33951073e+00,  8.74391606e-01, -1.90858384e+00,\n",
    "       -2.82342206e+00, -3.42489523e+00, -1.51987052e+00, -4.27835205e+00,\n",
    "       -3.36605751e+00, -3.24301280e+00,  2.68654946e+00, -6.20634775e-01,\n",
    "       -2.79786651e-01, -2.85994399e-01, -8.37594278e-01, -1.26462998e+00,\n",
    "       -8.58365794e-02,  2.14254504e-01,  6.47751816e-01,  1.01457125e+00,\n",
    "        9.53587635e-01,  5.77756800e-03, -4.49344337e-01,  5.81501931e-01,\n",
    "        1.28945119e-01,  7.06964359e-01,  1.00089889e+00,  3.03851273e-01,\n",
    "        4.86533088e-01,  1.45005787e+00])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the two datasets and show a 1gCm-2 uncertainty around the observed data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/vm/miniconda2/lib/python2.7/site-packages/pandas/plotting/_converter.py:129: FutureWarning: Using an implicitly registered datetime converter for a matplotlib plotting method. The converter was registered by pandas on import. Future versions of pandas will require you to explicitly register matplotlib converters.\n",
      "\n",
      "To register the converters:\n",
      "\t>>> from pandas.plotting import register_matplotlib_converters\n",
      "\t>>> register_matplotlib_converters()\n",
      "  warnings.warn(msg, FutureWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1122706d0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEsCAYAAADdO/TjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXd8G/X5x9+nPWx5O05iO3ZCErKcRRIgAQIUKKVNy4aWskpDGWUUKBRaCqWlQH8FCpRVZiFsyi6jjBSShkD23sN2vG1ZkrVP9/39cZESx7ItyfJK7p2XX0ms093XsnTP91mfRxJCoKGhoaGh0Rm6/l6AhoaGhsbARjMUGhoaGhpdohkKDQ0NDY0u0QyFhoaGhkaXaIZCQ0NDQ6NLNEOhoaGhodElmqHQ0NDQ0OgSzVBoaGhoaHSJZig0NDQ0NLrE0N8LSAf5+fmirKysv5ehoaGhMahYvnx5kxCioLvjDgpDUVZWxrJly/p7GRoaGhqDCkmSdidynBZ60tDQ0NDoEs1QaGhoaGh0iWYoNDQ0NDS6ZMDmKCRJ0gPLgD1CiO/393o0NDT6hnA4THV1NYFAoL+XctBgsVgoLi7GaDSm9PwBayiAa4GNgKO/F6KhodF3VFdXk5mZSVlZGZIk9fdyBj1CCJqbm6murqa8vDylcwzI0JMkScXAacBT/b0WDQ2NviUQCJCXl6cZiTQhSRJ5eXk98tAGpKEAHgR+DSj9vRANDY2+RzMS6aWnr+eACz1JkvR9oEEIsVySpLldHDcfmA9QWlqa3kXU1MDHH0N0TGz0RZYkyMuD004D3UC1sRoaGhrpZSDe7WYD8yRJ2gW8ApwgSdKLBx4khHhSCHGEEOKIgoJuGwuT48IL4dJL4Wc/U78uvVT9uuQSmDcPfvxj8PvTe00NDY0BzWWXXcaGDRvScq6ysjKampq6PObuu+9O+rzPPfccV199darL6pQBZyiEEL8RQhQLIcqA84DPhRAX9NkC/vtf+OwzuOsu2LWr49e998Jrr8EJJ0BDQ58tS0NDo3956qmnGD9+fJ9dLxVD0VsMOEPRrwgBv/sdDBsGN9wAI0Z0/Pr1r+GNN2D1apg1C9K0w9DQ0Bg4eL1eTjvtNCZPnszEiRN59dVXmTt3bkwqKCMjg5tvvpnp06fzne98h2+++Ya5c+cycuRI3n33XaDj7v773/8+Cxcu7HCtH/3oR0yfPp0JEybw5JNPAnDLLbfg9/uZMmUKP/nJTwB48cUXmTlzJlOmTOHyyy8nEokA8OyzzzJmzBiOO+44Fi9e3Cuvx4DLUeyPEGIhsLDPLvjpp/DVV/DII2C1dn7cGWeonse8eXDUUarhOOmkPlumhsYhw3XXwapV6T3nlCnw4INdHvLRRx8xbNgwPvjgAwBcLhePPfZY7HGv18vcuXO59957Of300/ntb3/Lf/7zHzZs2MBFF13EvHnzEl7OM888Q25uLn6/nxkzZnDmmWdyzz338Mgjj7Bq78++ceNGXn31VRYvXozRaOTKK69kwYIFnHTSSfz+979n+fLlZGVlcfzxxzN16tQUXpSu0TyKKFFvoqQELrus++NnzIClS1Uv49RT4R//6P01amikSENbA/6wlldLlEmTJvHpp59y880389VXX5GVldXucZPJxHe/+93YsccddxxGo5FJkyaxa9eupK710EMPMXnyZI488kiqqqrYunVrh2M+++wzli9fzowZM5gyZQqfffYZO3bsYOnSpcydO5eCggJMJhPnnntuyj9zVwxoj6JP+fe/1Rv/k0+C2ZzYc0pLYdEiOO88mD8f/vxnmDsXjj9e/bukpDdXrKGRMJ6QB5PBhNXYhac8EOlm599bjBkzhuXLl/Pvf/+b3/zmN5x88sntHjcajbGSU51Oh3nvPUOn0yHLMgAGgwFF2VfhH6+PYeHChXz66acsWbIEm83G3Llz4x4nhOCiiy7iz3/+c7vvv/32231SSqx5FKB6E7ffDiNHwsUXd3lotbsaT9Cz7xsOB7z7Ljz+uOrSvvOOWjVVWgqHHaZ6J3F2CBoafUkwEiQcCff3MgYNNTU12Gw2LrjgAm688UZWrFiR9DnKyspYtWoViqJQVVXFN9980+EYl8tFTk4ONpuNTZs28fXXX8ceMxqNhMPq7+zEE0/kjTfeoGFvAU1LSwu7d+9m1qxZLFy4kObmZsLhMK+//nqKP3HXaB4FwNtvw4oV8Nxz0IUWijfkZadzJyNzRpJpztz3gMEAl1+ufikKrF0LX3wBCxfCK6+o5162TOu90Og3ZEUmIGvaSYmydu1abrrpJnQ6HUajkccee4wbb7wxqXPMnj2b8vJyJk2axMSJE5k2bVqHY7773e/y+OOPU1FRwdixYznyyCNjj82fP5+KigqmTZvGggUL+OMf/8jJJ5+MoigYjUb+/ve/c+SRR3LHHXdw1FFHMXToUKZNmxZLcqcTSUSbygYxRxxxhEh5cJGiwOTJEArB+vXqTT8OQgjWN67HF/aRYcpgfEGCZXIvvQQ/+Qk8/7zqaWho9AMra1eSYcpgdN7o/l5Kt2zcuJFx48b19zIOOuK9rpIkLRdCHNHdc7Ut7uuvw7p1cMcdnRoJgGZfM23BNrLMWXhCHhSRoLrIeeepie/bbgOfLz1r1tBIklAkRFAO9vcyNAYph7ahiERUAzFhAnRRLRCOhNnZuhOHxaEmjgSJf+h0Ovi//4Pq6n5LzGkc2ihCQQihhZ40UubQNhQvvQSbNsGdd3aZP6h2VwNg0O3zOJL60B17LPzoR2pVVH19ysvV0EgFRSjoJB0REUncE9bQ2I9D11AIAXffrVYqnX56p4d5Q17q2+pxmPeNxTDqjbiD7uSud++9EAioHoyGRh8SUdTkpoSkVT5ppMShaygkCT74AJ56qlNvQgjBztadWI3WdrXKJr2J1mBrctcbMwauuEJtzNNkPzT6EEUoiL1/ZEXu7+VoDEIOXUMBat/E9OmdPtzsa6Yt1NahScmkNxGUg8l/6G6/Hex2uPnmVFaroZESilBiGx3NUGikwqFtKLoglsA2dz6JNenkYH6+Wv30/vvw+ec9XKGGRmJE8xKaR9Ezdu3axcSJE/t7GR3YX6ywt9AMRSfES2C3Q5Cads4116hd2zfeqPZwaGj0MhERQQiBUWfUKp8GGFG5j4GOZiji4A/7qffWd+lNmAwmXEFX8ie3WNTqp5Ur4cUO85g0NNJO1KPQ6/SaoUiC+++/n4kTJzJx4kQe3FvaLssyF110ERUVFZx11ln49vZG3XLLLYwfP56KiopYB3djYyNnnnkmM2bMYMaMGTEJ8DvuuIP58+dz8sknc+GFFzJr1izWr18fu+7cuXNZvnw5Xq+XSy+9lBkzZjB16lTeeecdAPx+P+eddx4VFRWce+65+PtgiJom4RGHKncVJr2pS7Ets96cfOVTlPPOU3sqbrlFLZt1dG6QNDR6SkSJIEkSBp2BYGRwNd098s0jbGvZltZzHpZ7GFfP7HoK3PLly3n22WdZunQpQghmzZrFcccdx+bNm3n66aeZPXs2l156KY8++iiXXnopb731Fps2bUKSJFpb1UKXa6+9luuvv545c+ZQWVnJKaecwsaNG2PnX7RoEVarlQceeIDXXnuNO++8k9raWmpqapg+fTq33norJ5xwAs888wytra3MnDmT73znOzzxxBPYbDbWrFnDmjVr4kqDpBvNoziAtlAbLf4WMkwZXR6n1+mRFZlQJJT8RXQ6ePhhqKvTymU1ep1wJIxO0qGX9Fp3doIsWrSI008/HbvdTkZGBmeccQZfffUVJSUlzJ49G4ALLriARYsW4XA4sFgsXHbZZfzrX//CZrMB8Omnn3L11VczZcoU5s2bh9vtxuNRBUXnzZuHde/Mm3POOScm5vfaa69x9tlnA/DJJ59wzz33MGXKlJiqbGVlJV9++SUXXKAO/ayoqKCioqLXXw/NoziAytZKLAZLwscH5AAmvSn5C82apUqTP/SQqljbB79sjUMTX9hHfVs9o/NGE1JCCCH6RJo6HXS38+8tOtPAO/B1kyQJg8HAN998w2effcYrr7zCI488wueff46iKCxZsiRmEPbHbrfH/j18+HDy8vJYs2YNr776Kk888URsDW+++SZjx47tdh29jeZR7Ic76MYVdGEz2hI6XkLCG/KmfsG774acHLW/Qktsa/QSX+z6gru+vAtvyKs23Sla0113HHvssbz99tv4fD68Xi9vvfUWxxxzDJWVlSxZsgSAl19+mTlz5tDW1obL5eJ73/seDz74YGwq3cknn8wjjzwSO+eqLib1nXfeedx33324XC4mTZoEwCmnnMLDDz8cM1orV66MrW3BggUArFu3jjVr1qT/BTgAzVDsRQjBrtZdXYacZEVmU9Om2P/Nhh7kKQByc+G+++B//1PVZTU0eoFGXyMREaHF34IQWolsIkybNo2LL76YmTNnMmvWLC677DJycnIYN24czz//PBUVFbS0tHDFFVfg8Xj4/ve/T0VFBccddxwPPPAAoE6uW7ZsGRUVFYwfP57HH3+80+udddZZvPLKK5xzzjmx7/3ud78jHA5TUVHBxIkT+d3vfgfAFVdcQVtbGxUVFdx3333MnDmzd18MNJnxGC2+FrY0byHXltvpMQvWLuCjbR9x94l3U+IoQREK7qCbGcNmpO4KKoqqBbV5s/qV2/n1NTRS4dbPbuXTHZ9y09E3UZpVyriCcV1W9PU3msx476DJjPcQRSjsdu0mw9y5N1HlruKT7Z8AsL1lOwA6SYcQomeVJDodPPooOJ1w662pn0dDoxOi4dEWfwugdWdrJI9mKIAmXxNBOdhpUloIwfOrnsdqsGIz2Njh3NHu8R7XpldUqI14Tz4JccYlamj0BG/YiyRJtPhb0Ou0yqeDiaiEfG9zyBsKWZGpclXhsHTuii+pXsLm5s2cO+FcRuaObFfXbdAZ2s/QTpU77oChQ9XEdi+MMtQ4NBFC7PMoAi3opcHRdHcwhMT7gpAcSkg6vqev5yFvKBq9jchC7lSqwxf28fLalxmZM5K5ZXM5LOcwqt3VsV2Z2WDGFUihQ/tAHA544AF1vnYXSS8NjWRQhIJfVjt3nX6n2nQ3wD0Ki8VCc3OzZiy6QVEUZNF9GFEIQXNzMxZL4mX/B3JI91EIIah2V+Mwde5NvLXxLVxBF9cfdT2SJDEyZyQCVX788PzDMelNOP1OIkoEvU7fswWdfbYqe/6b38Axx2i9FRo95kBDodfpB3x3dnFxMdXV1TQ2Nvb3UgY0siIjKzImvQmd1PWe32KxUFxcnPK1DmlDAaCgdHqDr3JV8cmOTzih/ARG5owEiP29w7mDw/MPjx0bjASx6RLrv+gUSYKnn4ajjoLvflctmy0r69k5NQ5pFKHgC6t6RC1+NfTUFmkb0E13RqOR8vLy/l7GgEYIwYraFQQjQQ4rPIxMc2avXm/AhZ4kSSqRJOkLSZI2SpK0XpKka/tjHUIInl/9PDajjbPGnxX7fpYli3xbfoeEdkpKsvEoKYGPPwa/H045BZqa0nNejUMSWZEJhAMYdAZ8sk/1JoSqKKsxePGEPMhK5yHzdDPgDAUgAzcIIcYBRwJXSZI0vq8X8b+q/7G5eTPnTTivQxPeqJxRbHduj/3fqDemJ08RZcIEeO89qKyE004Dbw+6vzUOaXxhHwoKQzOGAmr4SSC0kaiDnEZvI0a9sc+uN+AMhRCiVgixYu+/PcBGYHhfrsEX9vHyupcZlTOKY0cc2+HxkTkjafI1xYyDxWDpWYd2PObMgZdfhmXL1NxFWPtgayRPW6gNgGKHGp/WeikGP7Ii0+xvxm60d39wmhhwhmJ/JEkqA6YCS+M8Nl+SpGWSJC1Ld9LrrU1v4Q66uXjKxXHjuPvnKYCYfHPaP3w/+hE89hh8+CH8/OegVYFoJEl0AxM1FM6AEzi4DEVEibCjZQdbm7dS66mlNdCKP+wnogzM8Jo/7KfV35ry810BV5/nmAZsMluSpAzgTeA6IUSH7boQ4kngSVAlPNJ13SZfE5/u+JRjRhxDWXZZ3GPKs8uRkNjh3MHUoVNj3/eH/elPKs2fD7W1+/os/vzn9J5f46Am6lEMyxwGqB6FTtKlJo8/AFGEwnbndpx+JxaDhdZAa7u+AovBwqjcUd2ODegrWv2tbGnZgoTE1KFTU8ox1LXVYTV2VKTtTQakoZAkyYhqJBYIIf7Vl9f+18Z/ISFxxuFndHqM2WCm2FHMNue+xjtJkvCFfb1TfXD77aqxuOceGDYMfvnL9F9DIy6yIhNRIpgN5v5eSkq4g24kJLIt2WSYMmLd2YOh6a47FKGww7kDp99JjjUHoMMNtC3URq2nltF5o/tjiTGEENS11bGrdReZ5ky8IS8NbQ0McwxL6jwBOYAn5CHX2reacAMu9CSp/tTTwEYhxP19ee0qdxWLKhdx0qiTyLPldXnsYbmHsdO5M9YUZNAZem+XJknw97+roahrr4U33uid62h0wBvy0uBt6O9lpEx0XK/NaCPHkrOv6W6A91J0hxCCHc4dNPuaY0YiHnajnRZ/S796UBElws7Wnexq3UWONQeT3oTD7KDaU510UYHT70TXD7ftAWcogNnAT4ETJElatffre31x4TfWv4HFYOEHY37Q7bEjc0biDXup99YD6myKXq0k0evhpZfUHosLLoAvv+z+Odu2QWvqsVANddfqDQ/eqrOovIzVYCXXmhvrpRjMHoUQasNrk7epSyMBewf8SNAa6J/PQSgSYlPTJhp9jeTZ8mKNcXqdHgSx+0ciRL0Su6nvkthRBpyhEEIsEkJIQogKIcSUvV//7u3rbmnewoq6FXx/zPcTimeOyhkF7Eto6yRdQu30PcFvAP+br0J5OcybB+vWxT/Q54Nf/QrGjIELL+zVNR3s7N+wNhjxBD1IkoTVuNdQBFoGhYxHZ0TnxjR4G7ocCbA/dqOd2rbaXl5ZR3xhH2vr1xKMBMmxdDRoDouDGk9Nwt6ON+wlGAn2aVlslAFnKPoDIQSvrX+NbHM2p4w6JaHnDMschklvikmOR2do9yZ7PHtYr9QRfP8dsNng1FOhqqr9QV98AZMmqbpR48fD++/Dzp29uq6DGVmRCcmhAVtB0x1t4TYMkgGX34XdaMcddBMREQRi0P1MQggqXZXUtdXFvfF2hklvIhAO9GwaZZIoQmFbyzYMOkOnG0+dpENCor4tMa+iydeU2tjlNKAZCmB1/Wo2N2/mh4f/sNukZVTWV6/TU55dHmu8k5B61VAoQolVrGzOCBB+/11wuVRj4XSq/54/H044QZ1xsXAhfPSR+u+9M3g1kkdWZMIiPGjLST1BD1ajFYPeEBtW5PQ7EUIMupGo9W311HhqyLXmJl0aqtfpafY19+j6yYTrGr2N+EK+bquTMs2Z1HhquvXwIkqEBm9DwmOa080hbygUofDa+tcYYh/C3LK5XR4rhKDF3xKrTR+VO4rdrt3IioxO0vXqDs0X9iGEIMOUQUgOsa3ETuRfb8KWLXDyyar38PTTcNNNsGYNHHccFBerIaqnn4bA4I1J9yeyIiNH5EF3U43iC/uwG+3odXqyLFnA4Gy6a/W3srN1J9mW7JT6B+xGO/Xe+pQ/o7Iis6FxQ0LGJigH2e3aHXu9u0In6dBLeura6ro8zh10I4RoJ/7Xl+q6h7yh+HbPt1S5qzhr/Fnd1jS7g26yLdkI1F/QyOyRyIpMpatSzVH04gfPFXDFxAsdFgfukJud08oRzz+vdm/n5cHSpeoMbut+u5grr1T1orRKqZSI/k4Hq+RFW6gNm9GGEIJsSzYw+Jru/GE/W1q2kGnOTFmhWa/TowgFTyi12TGeoIeAHGBby7ZuQ1i7XbvRS/qE15ppzqSurS6+VyEEwYZa6mq3YjG0lwlfUbeCu7+6u0+q8gZkH0VfEY6EeW/Le5RllTFr+Kwuj5UVGYGgNKuUDY0bANWjADWhPTJnJIpQUITSreRvKjT5mrAa9hmAHEsOTb4mzN+bQ8m2baqYoClO/PKEE9Sk9qOPqtVSGkkREREMOsOgrBISQrTr7YnG9Vv8LUhSL1fppYlwJMyW5i2Y9KYex+fNBjP1bfUxg5kMDd4G7EY7kiSxpXkLEwsnxk0qt/pbafY1d1teHyMcxrSnjuwdW3E3fElBs1/NO1ZWolTuhqoqzF4fYy1m3Of8iNZfXERojHrfWVS5CFfARZ41wWv1gEPaULy35T1a/C38YvovunVnPUEP5TnlZJgyMOgMRJQIedY8HGYH21u2852R3wHoFUMRlIME5ECHUsAcSw5V7iqMhSMpimckQM1RXHEFXH89rFwJU6fGP04jLhElgklvGpSVTxERwS/7KcwoBNQuZbPeHOulGOjGL9pQF46Eu5xAmSg2ow2n30koEkrK6IQjYVoDrbGwlyfoYYdzB2PyxrS7b8iKzA7njnZNt5LXh7G6BkN1Dcaq6N971O9V1WCoa0BS2k+oE0VFhIYX4R05jNAx09CVjsC0ZQeO194m+5+v4j1hDpU/P5dV3lXMLp3d8zk4CXDIGgohBO9vfZ8xeWOYWDixy2MDcgCr0Uq+LR8Ah9kRc+kPyzksViIrISU0ljBZ2kJtEMeOSZJErjWXHc4dmA3mzmvKL7oIbr1V1Y168sm0r+9gJmoo0iYj34coQsEf9mM1WNFLeiJEyLXm4gw4B0UvRZWrCqffmXAZbKI4/U6GZAxJ+Phod3vUKGSaM3H6nezx7IlpaAHUemqRhUxm2Ej2358k5/F/Ymhor0MnDAbkYUMIFw/Hd+xRyMXDCJcMI1wyHFdhFlJJKV69jIREpjmz3aaz+bbryXr+FbKffolNd/4Sw/EWTmuogBl+6OV5FIesoZAkib+f+nf+W/nfLr2JqPs+oWBC7JeWZc6ixdeCzWijPKecFXUrYjvO3jAUTb4mLPr4Ywx1ko5McyZ73Hs6NxQ5OfDjH8OCBWoOIzt51/tQJSIiGPXGwWsoZNVQGA1G9IqebEt2TMZjIPdSNHobYxVO6cRuslPXVkehvTDhpHiDtwGLsf3nL9uSTZWrCrvRTo41B1/Yxx53NcWffUvhHX/BtLsa74nH4DzqpzFjIJcMRx5SoDbPxsEoBN6wF4fRETcqEcnPpeWGK3FefRlvvXAZozZu5pjb/4Hv3OsguzD5FyMJDulkdrQJqSs8QQ+F9sJ27qTVaI3t8KONdzudaq9Cug1FRInQGmjtkMjaH7PBjDfs7bpx58or1Ua8f/4zres72IkOh1GEMmiSv1FCcohgJIjFYMGoM2I2mGOGYiCHnlp8LWxr2ZZyhVNXRL3DREOJQTmIO+ju8PmTJAmH2cHWlq34wj5qF33EhPOuofiSaxFWK9VvPMOeV5/Ced3leM76AYFZ05GHFXVqJKLnzDBldBu6rg42sjlPMOUXd7Du05cQw5PTi0qFQ9pQdIesyCgo7dxLoN2bpjxHHdm43bm9V5qYovIRiXxguqzGmDYNZs1Sk9qaXHnCREufe7tPpjeIKsdajVYMOgNmvRmH2RGTs4iISK94wKkihGCPew+bmjeRZcnqtdi7UW+kyZfY5Eh30B037Bs9j7XVi++SCxh54lnYNm2j/r7fs3vh2/jmzk7jitvzVeVX6CQdR5fOJnBYWa9dZ380Q9EFnqCHsuyyDokvg86A2WAmFAmRYcqgKKMo1niX7g9eq781ISlis8FMs7+bGu+rroLNm+Hzz9O0uoObdr9LafCVyEZLQU16E3pJj9VgJcucRUREcAVcva9PlgQRJcJ253YqXZXkWnOTkt+WFZkqV1XCHpLdlHhPRV1bXacDguwffc64E84h79X3aJn/U3Z+8wmuS38Mht6L6CtC4X9V/2PykMmxBsq+4JDNUXTHgQnsA8k2Z9Psb8akNzEqZxTrG9YjEGk3FM3+5oS0560GKy3+FiJKpPOd2Nlnq9VPjz4KJ56Y1nUejOz/uxyMnczRxlCr0YpJb8JitMRuLi2BFnItuciKjJn+lVAPRUJsbd6KN+xNvKx0L/6wn/v+dx/bWlTJ/wJbAcWOYoodxQzPHM6o3FEUZRS1e45O0sV6KroqlQ3IAXxhX4fcn87TRsFv7yZrwZsEJh5O9RvPEBo/NqH1rmtYx1MrnuKs8Wcxp3ROUj9r9PnOgJMLSvu21F0zFJ3gC/uYWDix03ihw+yIdVOOzBnJ4qrFtPpb02oo/GE/oUgoIbVISZI61M13wGKBn/0M/vpXqK5WO7dBDUXV1MDataAo8L0+Eesd8Oz/uzToDGpCu2/nxfSIqEdh0Vsw6AwYdIZYL4XT74wZiv7EF/axqXETSCTd3xCQA/x1yV/Z6dzJ+RPPJ6yEqXZXU+2uZm3DWmRFRi/puWPuHR2GkNmNdnY6d3baDwHEvK79sS5aStEvf4NhTy3N1/+C5puuit+/dABCCP6z4z8sWLsgpgYxc/jMpHtDFlUuwm60M7Wob8vcNUMRh7ZQG3m2vC5VZPfPU0RHo1a5q9pNvEvHOpLBoDPg9Du7Hp50+eXwl7/AddepQ5DWrFENREvLvmMWL4ajj05x1QcPQTnIP1b8g59W/JRsS/ag66WIlnWa9Ka4hkLkiH41FN6Ql/UN67EYLV0Wa8QjFAnxwJIH2NK8hatmXMWs4vYNs7IiU+Op4c6Fd7Jw10IunnJxu8fNBjPBQJAdrTsYkzsmbg6wzluHzaRqK0n+APl/eoCcx58jNLKMqg9eIjAjsc+6rMg8v+p5Fu5eyLSh05g7Yi73f30/n+38jFMPOzXhn9kX9rG8ZjnHjDimzxVktRzFAQghCEfCFGcWd3mc2WAGST0+WjnlDXnTGp5o8jUlNfLQarTS7G/uWgNm5Ej4wQ/gzTfh2WchGISzzoKHH4ZPP4XCQvjtb9Ow+sFPjaeGlXUrWVm3Er1Oj18eXCWybUF1o2Ez2mKGItow2uJvwagz9mvlU5OvCYPekJqR+PoBNjZt5IojruhgJEDdNJVmlTJ92HS+rv46bkWgw+LA6XdS46np8Jgv7MMf9qsd4es3UXrSWeQ8/hzOyy5g9xdvJWwk3EE39yy6h4W7FzJv7Dyum3UdU4dOZULBBN7f8n5SJcrf7PmGkBLimNJjEn7Dj3jBAAAgAElEQVROutAMxQG4g24K7YXd3qB1ko5MUyahSCgmreGX/WlLDsqKjCfkwaxPPH4cnbLX7Yf/5ZdV6XGXC5YsUdVlr75azVvceqsqVf7ZZz38CQY/0YqzJl9Tv99UU8EVVEMnVqMVnaRTE8SS2tHvDDj7dSSqEKKDLE0ihCNhHlr6EOsa1nHZtMs4quSoLo8/pvQYvGEvK2tXxn0825LN7tbdtPrbDzZq9beiQyL7yX9SevLZ6FtaqX79aRrv+R3CnpiC6+7W3dz+xe3sdO7kyhlXcvb4s2Oey5njzsQddPOfHf9J6Fyghp2GZgyNRTD6Es1Q7EdUqyk6iL47ssxZsTp1CYlQJJQ2V94b8iKEiL2xqtxV3Pb5bWxt3trl8ySkWBKzU2w2KCtT5T0O5PLL1dzFbbcd8mW00XLjRm9jLAc0UKqEEsEb9oKkFjroJF3MWORYcmK9FJ2NRFWEgjvoxh104wl68AQ9tIXa8Ia8aWnU84V9ag4hiRJYWZF55NtHWF2/mkunXMqxI47t9jkTCieQZ83jy8r4EyF1kg6HxcGWli0xoymEoHnXBsZdfBOFt/4J39zZ7P7yXXzHJ5583tW6iz98+QcUoXDbsbdxVHF7gzY6bzSTh0zm/S3vJ9TMWd9Wz+bmzcwpnZP23pJE0AzFfrgDboY7hnc7kyJKhimDiBJRJ4gZrATkQNr6KJwBZ7s45OLKxVS6Krn/6/u7lCS2Gq0J14jHxWKB229XlWjffz/18xwERHMSjT5VhkGSpEFV+RT1SKNGAlB7KSyOmIxHvJt+RImwtXkrGxo3sKlxExubNrKxcSMbGjawvnE9q+tX97hT3R10J3XDE0Lw5PInWVG7ggsrLuT48uMTep5O0jG7dDZr69fG5NUPxKQ3YdQZ2dK8hYgSIfje2xx+4jnYF39D/X23U/PiY0TyE+8QV4TCMyufwWqwcufcOzv1AM4cdybesJePt3/c7TkXVS5CQkqpUiodaIZiL9Eb/hB74how+8dWrUYrftmfFo9CCEGzr7mdW76qbhUljhJ0ko6/LP5LpzOAzQYzbaG2nu18L74YRo2C3/1OrYI6RImGnhp9jWreRwweaW5Q+4DsRjsSUmznbjaYVQmavUOwQkqoXU5LVmS2NG+hNdBKrjWXbGs22ZZs9e+9/zboDOx07uzRPIQGX3JDeN7Z/A5Lqpdw9vizOWnUSd0eb6jaQ+Gv76Rsxsn88K1NEAqxuHJxp8fbTXbEli14L70Ayw/PIFyQR+Wnb+K69CcgSXy751s+3PZhQmv9fOfn7GzdyU8m/aTLmd7lOeVMHzqdD7d+2GXhihCCxVWLmVAwIe2SJomiGYq9uINuShwlSVUTmPSmmM69zWhTPQrRc48ianCiH+5GbyN7PHs4ZsQx/OqoX9EabOX+Jfd3GV9OtmKqHUYj3HEHrF59SM+xiIaeQpGQOjgGQUhObL7xQMAX9sWqdqIehdWoNt2FIiG8Ya/adLfXSwpHwmxu2kxbqK3LG1yGKQNX0BXztJIlKAcJhAMJl4Yuq1nGmxvfZE7JHH4w5gddHmvcup0hV99C+YyTyXrhdeQhBUx4cAEzPlzNt289gmHH7vZPiESwf/gZw8/+GZOPPZuMBa9T+7Nzqfz4dUKHjwbUG/ULa17gpbUvsXDXwi6v7/Q7eW39a0wsmMiRxUd2+7OdOe5MfLKvSyO0uXkzjb7GfvMmQDMUgPoBMegMFNgLknpeVO8lKAdjoad0xLA9QU+7+u3V9asBmDxkMqNyRnH1jKvZ1bqLh795OO4O16Q3dd+l3R3nn69Ozbv9dpAHzy46nexfDtvoa1TFAQdR5VNbqA270Y5AtAs9Rcunnf59A4xCkRAbmzYSkAMJTWbLsmSxq3VXSvkKT9DTqSzGgVS5qnh82eOMyhnFJVMv6TRcZV6zgaGXXkPZ0aeR+c6HtP7sx+xc/inV7y9g15IPOXb4bOqbd+P7wSkM/dl1WBd/Q86DT1A+/TsM/+mVmDZtpemWa9i28nNq7/o1Ots+b2dz82acASc5lhyeX/18TIUhHgvWLkBWZC6aclFCobWSrBJmDZ/Fx9s+7pBbFELwdfXXPLXiKSwGC0cMOyKxF60X0AwF6gdqRNaIlLRlogltq9GKP+xXB9f3MAnc4m9pV3W1um41hbZChmYMBWDq0KlcMuUS1tSv4dlVz3a4ntVoxel39qz5T6+Hu+5SJT9efDH18wxiojtuUL06g84waAyFIhQCckAN70igl9T3tklvioUvWvwtCCHwhrxsbNqIrMhd9+Dsh0FnQC/p2d26O+n3e7O/uVM15P1xB93c//X92Iw2rp11baceSN6fH2TECadj+2IxLdddzs6VX9D4p9tUET4gfFg5Y+54mMjMGbx10Uxsn39FyQ9/SsEf7ydcXkrNcw+zc+XntNx4FQwt6tA/9XX115j0Jn5/3O/JseTw0NKH4haMrKlfw9I9S5k3dl6HbvCuOGPcGYQiIT7Y+gGgGohlNcu47fPb+Pu3f8egM3DNzGvi5k77ahyq1nAn1F1Wqpr30TGTNqONurY6JEmdSRH9YKaCL+yLdWOHIiE2NG3guBHHtduhHF9+PC3+Ft7e/Da5llzOHH9m7LGoRIE35E34gx+X00+H6dPhzjtVmfIEOlAPJrwhLwadAYGg0de4rzt7EBBRIvjCPrXiCV3svWPQGWId0C3+FkqzSql0VcaUS5Mh05xJs6+ZfH9+wp8fWZFxBV1kmbv2WmRFVm/IATe3HXtbp6Ewx4uvk/fXx3CddzqNf7oVJSu+/pHVaGVG+Rw+syznjBWfkPufrwhOmRSbFtfVOpbuWcrUoqnk2fK4ZtY13PXfu3jkm0e4efbNsc1lKBLi+dXPU5RRxGmjT0vgldjHsMxhzC6Zzac7PqXUUcpH2z5il2sXRRlFXHnElcwqnhVXISJaFZlMr1WqHPIehVFvZET2iJSn0kUT2laDNbbb7MlOXlZkIkoktp6NjRsJRUJMKZrS4dgzxp3BsSOO5e3Nb7OqblW7x/Q6Pa6gK+V1ACBJ8Mc/wq5d8PTTPTvXIMQX9pFhysBhdsQ8ikAk0KdD7VMl6lFElWOjGHQGHGYHEhLOgBO7yR5rxEsFh9nBduf2riXu9+PAsu94CCF4ftXzbG7ezM+n/7zTqiHrl0sYcuMdeI+fQ/0Dd3VqJKIcM+IY/LKfZW1b8Jzzo26NBMCGxg20hdpi/Rpl2WVcMvUSNjZt5NX1r8aOe3fzuzR4G7hkyiUpdU3/6PAfISsyjy9/HF/Yx/xp87nnxHs4quSoTu9N3rA3aQHFVDmkPQpJkhiWMSylGbpRjHojZr0Zs8GcluFFB37gVtevxqQ3cXj+4R2OlSSJS6Zcwqq6Vfyv6n/tjInVYKXJ29RBIj1pTjkFZs9Ww1Bnnql2bh8i+MI+LAZVSC9WciwgrIR7PL+5t4koEfxhP2aDuV1INdqhHZ1LYdKboAdq3ka9ESksUemq5LDcw7o9vtnX3O2N9JPtn8Q6mTtLCBu3bmfYJdcQGlVG7dMPqgUY3TAufxwFtgK+3P1lt416UZZULcFutFNRWBH73pzSOexw7uDDbR8yMmckpVmlfLD1A2aXzGZ8wfiEznsgQzKGMH/afGQhM7tkdkI3/1AkRKG9bz6PA9KjkCTpu5IkbZYkaZskSbf05rWKMot63MDiMDsw6ozIiqzOsOihoRCoO1YhBKvrVjM+f3ynNyaDzsDkIZNZU7+mXQ+HUW8kGAn2vPNWkuDee6GhQW3Su/56VUDwECBqKArsBTT4GmLfHwxNd37ZjyxkrAYrJt2+945ep0dCihmKdOCwqB7Xgd3NB6IIhZZAS6dlsUII3vvfs7z2zp+Ys0fP2SXxdZB0zS0MP/9yhMnInpefQHEkFl6VJLUPYX3jepp93Rd7hCIhltcuZ/qw6R2M248n/ZgxuWP4x4p/8Ni3j2HWmzl/4vkJraMzZpfO5rgRxyVkJCJKpEeeYLIMOEMhSZIe+DtwKjAeOF+SpNTMdB/hMDtiv1xfyNejEll/2B9zNWvbamnwNTC5aHKXz5laNBVv2MvWlo5d256gJ+W1xJg9G9avh3POUTWhysvhiivUkNRBTHTedKGtkBZ/C4pQkKTBMcAophxrsHQo0jAbzGRZsnAGnGm7XqY5k+3O7V0aUV/Y1y6suj/hNg/PPXgx7791H6f8r5E/37WEUdNPJPeBx9F59pV6S8EQwy+8CkNtPTUvPIpcmpzHPKd0DgLBV5VfdXvs6rrV+GV/h65qUDdoV8+8GrvRzi7XLs6ZcE5C1WLpoi3URlFGUcoh82QZcIYCmAlsE0LsEEKEgFeAH/bzmrrEarS203vqiUfhDXtju5doWWy8/MT+TCyciEFnYGVdez0bm9FGoze1WvcOjB0Lzz0HW7fCJZfAM8/A6NHqvxvTdI0Bhl/2x2aSyIocKydNNB7fn+xvKA7coVoMFrLN2bGfpwOhEMPPvYySk89m2HnzGXL1LeTffi85f3sSx4uvY9xZ2eEpUY93W8u29u//QAD++1+oqcEVcHW8sQmB/O47/P3Xx/F19df8REzk/Ac/pe6DV/FPm0z+nx6gfMrx5N37MDpnK0OuvRXr0hXUPXIvgSO6/lzEo9BeyLj8cXxV+VW3uaavq7/GYXYwLn9c3MdzrDnccNQNnDX+LI4vS6xTPF1ERKRPm+8GYo5iOFC13/+rgY7ykAMIi8ESqzwIysGeGwqdaihW1a6iOLO40+FJUaxGK4fnH87K2pXt3F+zwYzT7yQoBxOWJemW8nJ4/HFVYfb//g8ee0xVoH3ppfScfwDhl/1YDJbY69/oa6TYUTwo5MZjY1AN1g5hy6iMhzfsjfvecLz6DvbPvsJ35BEYGpowb9qKvsWJzre3WMNqoeG+3+M+/4x2z8s0Z9Lqb2V3627KssuQNm9WvdC1awEoLMzHUTGe0JRJBCZPIJKfi/vB+7grcwXOPBtXzfkNFadeTASIlAyn5pUnMa9aR95fHyXvL4+Q87cn0IXCNN16HW2nt5+ZEr3pJxJGPnbEsTyx/Ak2NG5gQuGEuMf4w35W1a3iuLLjuiybH5E9ghHZI7q9ZjoJRULYjLakOtt7ykA0FPF+0x1MvyRJ84H5AKWlpb29pi7R6/RkW7Jjg4NSNRTRSpUcSw7+sJ8tLVs4ZdQpCT13atFUXljzAnVtdR1quF1BF4WGNCe9iovhwQdVbaj77lPlPsbF33kNVqLlpdGEYYO3gfLs8kFRIusKqBVvZr25g0dhNVhxmNQKIWfA2f79Isvk/u0JAlMmUv3ei2qOai+Sz49hTw1DbrqTol/+BuuipTTcezsiY99grSxLFnVtdeS8+h7ZN/xGFaB89lnCLU20Lv4Mx/ptZHy+CElRWFIMfzjRgGHU4dxwzl8pz++YDA9OmUjNC49iWr+J3If+gTykkJbrf9HumFAkhCfoQadTZ5tHq6qMOqOq43RAfmHGsBm8ufFN/rHiH9x1/F1xS8iX1y4npITihp16ijvgJiIiXXa/d4U37GVkdt8qyA7E0FM1ULLf/4uBDtlTIcSTQogjhBBHFBQk11HdGxTaClGEgj/iT1kYMBwJxxq81jeuR1ZkJg/pOj8RJRqeOrBM1m6yU9tWm9J6EuLGG9WbwR/+0HvX6AeEEKoysNFCni0PCYlGb+OgmUsR7e43GzoaCqPeGBtgdGBCO/PN9zHtqqL5V1e0MxIAwmYlPHoU1W8+S/NNV+N47R1KTzoL04bNsWN0Xh/jbriH7CuuJXzENFUG5uKLab38Inb87Q52L/6AbTuXs/Sth7npFyPJ/s5p/PaS5+Iaif0JTTicuif+StMfbm63rlAkhDfkZWLhRGYOm8nkIZMZVzCO8uxyHGYHYSWMO9C+Oc5sMHPNzGtwB908+u2jcTd2X1d/Tb4tP6FKrkQRQtDiayHTnInFaEmpqz2qOdaX+RAYmIbiW2C0JEnlkiSZgPOAd/t5Td1SmFGIQBCUgylXxewf+15VtwqrwcrovNGJXd9eyPDM4R109016E/6wv/fCJfn58MtfwquvwoYNvXONfiCiRAhFQmrvhBwg15qrDtrZO/Mj3bPR0020iCEqMb4/Rr0xdqNpZygiEfLuf4zghLF4T+1iprpeT/PNv6T6zWfRuzyUnnw2jhdex7x2IyNOPIOsN9+n/qYrWfnCffgK1NLzJl9TLFQi7DYWFQYJD8nn2rm3pLyzjhqJcQXjyDRnIkmqYXSYHRTYCxiVO4rD8w9HFh2LD8pzyrloykWsa1zHGxva65l5gh7WNazjyOFHpk3SW1Zkmn3NDHcMZ0zeGMqyymKik8ngl/3kWHL6vDx7wBkKIYQMXA18DGwEXhNCrO/fVXVPpkl1X0ORUMpS1EE5iBACIQRr6tcwqXBSUs00U4umsrl5cwejoJf0aSuFjMsNN4DdflB5FdEPsV5Sh/vk2/LbieAN9BJZV9CFJEnYjLYOMfb9R6Lu/77IePdjTNt3xfUm4uE/9ih2L3wb/5FHUHT9byk98Qx0Xh/Vbz2P6+ZrMZtsbG5S34/uoLtdLmRl3UrKs8tTTsiGI2Hagm2MzRuLw9x5o53FYMFhcsQtEz9uxHEcX3Y87215j2U1y2Lf/7bmWyIikpCoXyIE5SDugJux+WMpySqJacTZjfaky9eDcpAhGYkrXKeLgZijQAjxb+Df/b2OZIjGOXsiNR6teNrt2o0z4Oy2LPZAphRN4f2t77O2fm278ZB2k536tnqGZw7vnaEnUa/innvUXMWE+AnCwURUOdasN6NDR6G9kHUN62KPh5UwZtJUINALtIXa0Em62DyK/THoDJj0JuxG+z5DoSiqNzFmFPXfncOLy55Qpch1OgySod3fRxcfHZsNHynMZ89rT5HzyFOYN22j8Q+3xGY3WI1WPEEP21q2tXvfuYNutrVs44xx7ZPhiRKOhPGEPIzNH0u2tftm2aKMIra2bI07cvWnFT9lt2s3Ty5/kmGZwxiWOYwl1UsYljGM0qye5z7bQm0IIZg4ZGK7ngdJkhiRPYINDRsSHgUbLS3ukSxPigw4j2KwYtKbMOlNBOVgyn0UvrAPo84YK4utGFLR7vG2YFuXfRGj80ZjN9o75CkMOoO6A+uJ9Hh3HGReRdSjMOqNSJJEga2A1kBrLI8UWbMaWrtuMOtPvCEvNoMqCHigodBLepAg15obK5HN+PAzzBu3sOvai7hn8X18vedrZEUmIAdUSfG9UvcbGjfwyLePUN9Wv++EOh3Oa+ZT9+h9HQb8ZJozY+rKUVbWrUQgui37joesyLiDbsbkjkk4ZOWwqB5HvHChUW/kmpnXYNQb+dvXf6PGU8Pmps0cVXJUjzdVnqAHk97ExMKJcRvjMk2ZZJozEy6O8Ia9DLEP6bPeif0ZkB7FYEQn6bAY1ARVyh5FyEuGOYPVdaspzy7vIC0SVsKxwUTx3ng6SceUoimsql+FIpR2byiTwUSjt7H3diN5eXDNNfDnP6texcSJvXOdPmL/8lK9Tk+ONQeBoNnfTO62GhynXggzZ8JXXyUUpulrvGEvNpMNSZI6CFRKkoRZbybbkq023QlB7v89Su3hxfwmZzmNnmaum3VdXI+2xd/CLZ/ewjMrn+GWObckdDM9MPG6snYludZcRmQlV1YaUSK4Ai7G5o1NSsTToDNQaC/EGXDG/dzk2fK4asZV3LPoHv701Z8QiB6HncKRMEIIxuaN7VSyRJIkSrJK2NC4ISFhPzkik2fL69G6UkXzKNKETtJh0VvwhX0p5Sii0h++sI9tLds6VDvJioxJb1KbfwSd7kKmFE2hLdTGtpZt7b5vM9po8jelbVRrXH71K8jIOCi8imiex2KwYJAMsVh6o7uO8hvuUg9avBgWLOivJXZJdBYFdPQoQA2pRSfd2f/zX1w7NnDVhXk0BVq48egbOw175lpzOX/i+Wxo2sB/d/836XWFIiHWNqxlatHUpEehtgZaGZkzMiWl5wJ7QZd5pfEF4zl3wrm4g27Ks8vjyoS7g+6EixjcQTflOeXd6lo5zA4cJke3xSahSAir0dqnvRP7oxmKNKHX6bEa1eFFciR5jyJa8bS2fm1ct9wf9pNnzcNsMHN4weEEI8G4HcKTCiehl/Qdwk86SYcQItax2ytEvYrXX481WQ1WojkKi8GC0WCMeXf+N17GvnoDex65G2bMgF//Gjy9+JqmgBACn+yLSeDHNRQGtenOFXThefRBrj7TirMoi1vm3MK4gq77YeaWzWVc/jheWvtS0kUSUTXkaUOnJfU8V8DF0IyhKSdy7UY7Jr2pS2PxvdHf4+zxZ/PjST/u8FgoEooZq+7whrxkW7ITTtSXZJUQCHed1PaGvBRl9FyXLlU0Q5EmdJJOnXIXCaQUeoqKAUaTbgdKK8uKHEvc2Yw2xuWPoy3Y1uFadpOdsXljO5TJgnpzaBdb7g1+9SvIzBz0XoVPVnd4ZoMZu8FOhikDYyBE28fv4T7tJOpPm6vqXtXWqsq6PcAVcKU1f6QIRdWp2hvOiGcoLAYLWeYs9K1ubizdSNvoUn5z7G0J9Q1IksTPpv6MiBLh+VXPJyW7vqJ2BRaDpVNZjHh4gh4yzZmUZqeeXJYkiaKMotgGoLNj5o2dF1ep2RvyUpJVQq41t0Nfxv4oQiEoB9XO9ARv6pnmTLKt2XHXJoTAFXChl/Q9UrnuKZqhSCM2oy1lqfGoGGCtp5ahGUPbvcmiH8RoKAHUN9eYvDG0+ls7XGvK0ClUe6pp8Da0+77VYKU10JqSVlGdp45qd3X3oavcXLj2WnXW9po1SV9noLB/6MlitKBTBCWrd1ObY6TxL3cgiwiRGUeoWlcPPqhOAkyRZl9zesQb9xIR6tAim8GGQWeIe8OK9hsYq/ZgMpi5+eyHKMsuS/gaQzKGcOb4M1lRt4Jv9nyT0HOEEKysW8nEwokJz2yIfi4Oyz2sx0ncXGtuTJk5FXIsOYzMGamOxO0k9OsKuijOKk56mFCJo4RgZF8DnhACd9CN0+8k35bPpCGT0ifDkwKaoUgjdpMdX9iHhJS0oWgLtWHUG6ltq2VY5rB2jwXkANmW7A718Lm2XMpzymMjLaNMLVJLF1fXrW53vCRJINHljuhAhBBUu6vZ0bqDPe49rGtc1+WuDFClyB0OOP54+MlP1Dh+U1PC1xwI+EL7DIVJb6LoH69QvNvJruOnEClUtZ/CSlhN3lutqnFMcaCRJ+xJa0Nk1KOIJwgYxagzMn15LWcs9XDH2MsZll+W9HVOGXUK5dnl/HPNPxPyiHa17sIZcDKtKLGwUzgSJiAHGJvfeUI4GaLGMVkJlqAcJMOUEetyH503Om7RSigSwqQzxUYWJ4PdZCfXkktbSK1sbA20kmfNY8rQKZTnlPerkQDNUKQVu9EeexMmayh8sg9FUWj2N3dIpAUiAfKs8asdhmYOpdhRHNP2AbVufGjGUFbUrYi7xkQlPaJGospVRZ41Ty1HFOps4Bp3Tec/Y24ufPghnHYa/Oc/cMEF6sCjI49Ux6pu2pTQ9fuTtnAbEhImvQnL9kqK//IEOcWHsadon1cXjoRhyBD1Z/r4Y3g3eQGBiBLBH/KnVT8qokRi8iOdjeQ11NYz/td/4XLfeKwX/zyl6+h1en4+7ed4Q14WrOk+qb+ybiUSUkL9QYpQcAVcjMkbk9YEblFGUdJNbj7ZxxD7vtyIzWjjsNzDcAVc7TZonqCHkTkjuxQR7IrhjuGEI2GyLdlUDKmgPKc84R6L3kYzFGnEbrTHZMaTMRRRMcAmv7rrPtCjQNBlWevQjKEdXOopRVPY1LSpww3IpDfhDXm7vTEJIdjduptqdzW51txY+MJqtJJjzaHSVcmGxg2dn+foo+Gf/4S6OvjmG7jjDrWM9M47EdOmEvky+YqZvsQb8mLSm5AUgfWKq1EsJixnnY8n5IndaGLVbVddBePHq56UP8ndaiSIQKRVP8oT9CAQqnKsIY7UQySC+ZKfowsEqX3yrwhz6nIQJVkl/GDMD1hUtYg19V2HGlfWrmR07uguO6mjOP1OyrLLUpb36AyH2RGba58wYl8vRpRcWy7FjuJYctsddFNgL+iRBpPdZGfa0GmMyh3VJ3Owk0EzFGlk/+7sZMpQo01c0Z3+/q5rOBLGYrB06Xoa9cZYY1OUqUVTkRWZ9Y0d1U/0On3nswhQDdfO1p3UeevaGYkoOklHri2XcCTMmvo1XU8L0+nU6qDbb4clSxCVlQSHFaGb90NYtarz5/Uz0el2Q557HWnRYmr+cBPZRWWAKjeu1+n3VaoYjWpie+dOVXo9CYJyEJ2kQxFK2gYiRSvbzAZzfI/ivvvQLVzI7rtuIDy65yqk88bOY3jmcJ5a8VSnVUHNvmZ2uXbFOro7QxEKLb4WhmQMYWhm8iGc7tDr9BTaC7sPn+7FH/bjMDviaisNdwzHYXbg9DtRhEKJoyTOGZIjHSG23kAzFGkk05QZUx1NZscSTS7XemqRkNqVAPplf0JNNgW2gnZx7miX9vKa5R2OtZvs1LXV0eJvwR104wv7CEVCRJQIilDY4dxBQ1tDXCNx4HnsJju1nsTVaYNF+Wxc8CDCkanO496yJeHn9iW+sA+roqfknsfg1FMJ/Pic2O620duIQWdo7wWccAKcdZaas9i9O+HrRHNTkD79qP2bBTvkKJYuVRsizzkH1/lnpKWvxqg3csURV+ANe3nw6wfjFktEh2pF82fxiA6HKs0qpTy7vNdKQfNt+Qn3OgXkQKdzqXWSjlG5ozAZTJRll/V7HqE36dZQSJJ0kiRJ/5Akacre/8/v/WUNTuwmNX7tDyc35S4qBljbVkuBraDd7iWiRBIqizuw49SgM3DEsCNYVruswwfXoDOABFubt7KpcRPr6texqnYVy2qWsaxmGS2+loSbmow6I17Zm/DPG5SDuIdk433/LTX5e9JJUF2d0K8MVFEAACAASURBVHP7El/YR0arF70/ADffjNVoiwnpNfmaMOqMHRO4f/2r+vf11yec2PaEVJkHSZJSFpM8EHdQLVYw682xIVgAuFxw/vnqLJEnnsBstKTNixmRPYIrjriC7c7tPLXiqQ4lsyvrVjLEPqRjWHUvQTmIJ+hhbN5YhjmG9Wq/gN1ox6K3dFv9F51r0VU4yaQ3MbFgIgW2/h910Jsk4lFcCdwEXCBJ0glA8gIthwjR2KtPTm54UVQMsNZTS1HmvkR2dEZzIsk8q9GKWW9utyudXTKbgByI61VkmDLIseaQbc2OfeVYc8gyZyUktBZFkiQkISVcchv1XsKjR8JHH6l6SSedNOCqonxhH5kte8MTU6diNVgxG8yY9KoUSrREsl0Io7RUzcW89VZCHdtCCNWj0BkRQqTNo3AFXUhIarNgNJQhBFx5JVRWqtMIs7Mx681pnf99xLAjOHv82SypXsI7m9+JfT8gB9jQuIGpQ+N3Y7eF2ghFQkwonJBS13WySJLE0Myh3YafopLe3Sk4R/XADmYSMRSNQohWIcSNwMnAjF5e06DFYXYgEATCgaQNhUEyqKWxGft2XNE5CInWj+fb89uFQw7PP5w8ax6LqxYnvJaU3vASCVeStIXaMBv2GrRp0+C992DXLvjud8GdeNlub+OX/WQ0tBIcOQIcDkwGE3pJT6GtkAaf2p9i1BvbSY8DqjjinDlqgruy42zp/YnWzUuS1DGU1QOiPRkWg2VfBc4LL6gG4ve/VwsN9j6eqoBlZ/xgzA+YUzKHNze+GeuvWNewDlmR45bFtvpbMeqMnQrn9RZ5tjzMenOXRR1BOUiB/eD2FBIlkTvQB9F/CCFuAf7Ze8sZ3EQb4gJyIKkwgi/kwxPyEIqE2nkUQTnYaVlsPLLMWe1izpIkcXTJ0axtWNuufDbd6CRdwslBT8iD1WDd54Ece6zanLd6NcybB4HkShd7C3/Yj6O2heBkVTLdqDOCpGoGNXpV42A32WnwNrSP8+v1arWXEHDRRaB0vmGIhhxBDQemq0TWE1Kn28WGFq1bpxquY4+FW2+NHWcxWNKu/SVJEpdOvZTRuaN5YvkT7HTuZEXtCuxGe4chXO6gmyxrFuMKxvV5fN+gMzAmf4wquRPHq4r+XqJzZg51ujUUQoh3ACRJyt/7/4d7e1GDlWgJXTCS+JS7cCSMIhTq2uoA2nkUsC/vkQh2k71D6d/RJUejCIWle5YmfJ5kMelNCXUWy4qMHJEx6o3tQ1WnnabeXL/8Uo2hy+kLh6RK0O8ms7mNYIWqgmvUG0GoRQNNvqaYhpKiKB2NcHk5PPQQLFwIDzzQ6TV8YR/GVg+Zr76NIZK+Ell30A0S2IxWLM++qKrc2mzw4ouqIduLUW+Mjd5NJ0a9keuOvI5MUyYPfP0AK+tWUjGkokMIJ6JEGJ45PKnhXOkk2g/R6m/tkFPxhr3k2fJS7ok42Eim6umZXlvFQULUdQ7JoYRjv1HPI1oaG222C0VC2Iy2pEYe6iQdudbcdjvTYkcxZVllLKpclPB5ksWkN+EJebrV/AnKQZDUeQgdchrnnw9/+xu8/TZccUXKXc7pQAhBuLUFqwzhyaqhiEph5NvU8F50XoXNZIvfwHjRRXD66eoOvhOBxPCXXzDxlJ8y9KqbyX/hTQJyICndpM7whryYFB2HX3k7lit/qYaaVq2Ckvblm715g3aYHdxw1A34ZT9tobYOZbERJYJep28nS9Mf5NnyGO4Y3qGsNxwJd1rtdCiSjKE4uLM1acBqsKJDl9SUu1hpbFstVoM1VuHkC/vIt+UnvYY8a16Hm/DRpUezs3UnNZ6apM+XCFFl2u7CbcFIMLYTj5v8/uUv4bbb4Kmn1BLOfiKiRJA9LmxhiEzd10VsNphjv5+ojpbFYMET9HQMG0kSPPEE5OSoMibBfT0uRCLwxz9SMu+nYDISmDSOvL89Cf7UBCUPJLR9M4VL1pD98X8J/+ku+OQTGNqxJ6ErQxGKhHD6nbQGWmn1t+L0O3EGnDj9zlhVVXeUZJXw/+29eZhcZ3Xn/zl3q7V39aKtW7Zl2ZZk4xUL8CKHJSYYnEAgBA+YMUsMiQMmkxDHz2+IAwyDhzGYJexkhglLIEyAZ0hCMOCFRQZJtvFuY1tLW2otLfVW+/L+/rh9r6q6q6qrW1VdddXv53nqUau2/nb17Xvu+55zvuemF97EBUMXcP7gHDfkfIrecO3y6+Vifdd6OkOdfhWbN8tlOXMm7c5iLilad4kXEEzDJGyHSefTde/9pnIpDGPWDLDjhBmgUqquDta5VDq4t63dxjce/gY/3/9zXr/59Yt+z7qYTWjXWgFNZ6axTRvLsKp7G33wg3D4MHz4w67tx5//eXP01iCVTyGJJHZnD0bviRxRxIr4vRRHk0d9h1/TMBlPjbPOXlf+Rv398OUvwzXXuIHv9tvhwAF485vhJz9h/NpXMPGJjxB68BHW/8H1DHz9u+RuuXjhpqtf/crN6XR3u4HIu3V3w6c/Tein3yK2LsLj//eLnPl717lNjxWwDKuiSV4mnyGVT3Fm35mYYqJQ/kpHoTg8c7jq8Ky5nDd43rxJjeBesbdqCM9cvH6Ihw89TLaQJVvI0h/rb8kkuXZlMYGi9aG/zfGsxjOF+qfczWRnsAy34mnzqs2Ae0VjGuaSPG5s0ybuxMnkM36CsCfSw5aBLfxy/y/5w3P+sGlXcV4XazW80ZBeHmXuFD7AvRL/7Gfdctn3vMc92f7xHzdFbzVmsjOYM0lk7aayq+5oSS+Fl9AGNziPzYyxpmPN/J/nVa+CG290O7a7u12n2ZkZ0p//DM++8iK6o3FSl11K8kWXsOYzXyV1019BT43f+6OPwuWXQ7Z6OfLRd29AXXYxiQu21DzZGWJgGVbZ7yGdT5PJZ9jSv6Vqfixmx3jo0EPki/klbV95/QmLyb81G8d02NS3iUePPIpSalFFJCuBxYTMW5qm4hTBFJOwFSaVq3/rKZlLUigWOJY65lc8JXPJk1qWr4qumrcVctn6yziSPMKT40u3w66FYzo1tySUUiTyCb8BTESqr7pM0y3lvPJKeMtbXMO9ZWTm6AGMdAYZHi47EUbsCI7pELNjZWWxpmFSKBaq//wf+xiccYa7rTY0BLt2MX3d6/0RqpOZKcbffxPO4XHML36purB8Ht76VteZ9+GHXRv3e+5x8zr/8A9wxx0Uv/svHLjiAsIxt1R7oavi0l6KZC5JrpBjy0D1IAHuFtyG7g1LrqRL59N0hbpalsSuRkeog9O6XSO+dgpi7UDdgUIp9UgzhZwKGGL4U+7qKY/1zADHU65XklfxlCvkTqrxqDPUSZHyssyL1lyEYzr8Yv8vFv1+Simen3q+ZqLVMZ2aVtOZQgZUSZ+GonYNfzgM3/sebNkCr3vdss62SD88O/TpjDPKTrS24TZW9Uf75/VPhK2wX7k2j1jM/Vk++lHXQuOcc/xtuCeOPsFN/3YTT5+3lukXX0T4jk9WNxa8/XbYuRM+8xl3Jvm557olr9de6waQm2+m+OprSOfTRC13VbJgoLDcQJHIJigWi2zu31zXSrY/2k9XqKvusuhS0oX0kvJvy8FAbIAtA7VXYiuRRX0aInKxiPyLiOwWkd+IyMMiEtzpNA1GRIhaUddBluKCFSzZQhYR8b2SSnsoItbS3SPDVnhe123YCnPJmkv41fO/WlQHsFKKr/7mq/z1j/+a239xe9XRl5ZhkS1Ur/YqNSwEd697wVVXV5fbvR2PuyfCXGM6lxci/djsIX1m+bAcv0S2pJfCI2JHmExPzvs5fTZvdsemRtzf63TO3YZ7avwpFIr9k/s5/Jd/inXoCHzuc/Nf/8gjbtf3618Pb3hDVe2Foju0yFv9LISXjBeELQNb6nYtFRFO6znN9whbFKpyLq0dEJFFVRquFBYbNr8G/APwOuDVwDWz/2pm8abc1TO8yKv8OTB9wDUDjA26wUU4qYNVROiL9M3bfnrx+heTyCV46NBDVV5ZTlEV+dIDX+KuZ+/iotUX8fT409xy1y38fN/PKwZBEal6okzmkvOu0uo6wQwNuVfQDzywaGfWpZJ54lGKjo3VNzBvReE13Xm9FKUIUtcM6UKxQDqXxjZt9k263dtHkkfIvORSpi97obvySJYk+3M5N1B2d7ufRQ28VWpZV3YNIlaEznDnkpreInaE4a7huqugwD3mw3ZtN2RN+7HYQHFEKfV9pdRzSqm93q0pygJKzIn5NfYLBop8FhRlZoC5Yo6IFTnphHNPpGfe9teW/i10h7rrsvTIF/N8dudnuXfvvbz27Nfynkvfw4d/58Os61zH53Z9jk/96lPzThBKqardxZ75XbaQ9beo6raPeN3rXGfW226Dxx+v7zUnQeq5pyjGo4TMUFmgMA0TA4NVEdd9dG7tfTwU5+D0wYX7SUpGXu6ddP98PEfafe+7AQ4dchP6HrffDrt2wd//vZvcr/Xe+Yx/DNWTA+iJ9LClf8uST9yD8UEidqTurvJkLslAVPcnBI3FBooPiMiXROSPReS13q0pygJK3Im77rHFhYcXJXIJLNOtePK893OFXEOakKJ21O9v8DANk23rt/Hg2IM18wm5Qo5P/epT7BjdwRu3vJE/OOcPEHHtz2+94lbeuOWNPDD2ALf8+BZ2HzwxRc82bWZyld/Xq3j6p0f+iQ/d+yEMMRZngvfpT7t7/W97m9uH0CxmZsiMjVKIRQlZoXkBO2yFT9iNz8lTWIZFtpj150FUw1t1pXIpDs0cKnuvxCUXUHzp77irikTCTVrfdhv80R+5wXIBprJTvs4y59gqmIZ5UkllQwzO6DnDH9i1EMVi8aSG+2haw2IDxX/GdY+9GnfLydt+aggi8j9E5InZ/Me/iEj9NqZtQtyOU1AFcsXcglfMiZxbBTQ2M+YHinwx35D9W9Mw6Qn3zLOFeMn6l5Av5n3DtrlkC1k+vuPj7D64m+tfcD2v2vSqsscNMXjVplfxd9v/ju5wNx/f8XHuH3XtQapVPuUKOb8T96nxpziUOFS96a4ag4Nu5/Yvf+kOCWoWDz5IyoJChRUFuNstXSH3RDc3TwFuFdGhxKGa38KbQTE6NYpCuVVU3nsJZP6/v4EjR1wbkOuvd3skPv3puuSXGgIuV1VRzIm543gztaugCsUCtmmfVP5N0xoWeyS9QCl1blOUuPwIuEUplReRj+KW5L6/id+v4cSdOCi3BHDBFUU2QaaQIVvI+lPtFKphc3JXRVdxbPxYWRXLSNcIazvW8o1HvsHde+6mJ9zj2o2HXZvxn+37GU8efZK3X/h2rhy5sup7r+9az23bb+M9//4efnP4N1y67lIc0/GnfZWeYDMF17qjUCzw/PTz5It5soVs2RZMXVx3HXzzm64txqtf7ZacNprdu0nZUIzFTpjqlRC1o3539jzn2NnHxxPjrOtYVzUxPJWZwjEdf9vp/KHz2TG6w83ZKMhcehGRV7zCLadVCr7zHVhVX5WQt5pZzkAB7lTG8dQ4qVyq6s/tuQ20Qze2ZnEsdkWxQ0Q2N0UJoJT6D6WUVwqzA1hX6/ntSDwUr8tqPF/Mo5Tytx5KB7o0ahxipVpwEeGt57+VF659IR2hDo4kj3D/8/fznce/w5d2f4mnxp/iXZe8q2aQ8LAMi5GuEfZOlKep5ia0SwczeXmTRDaxuBWFK96tCLJteMc7ajqzLpndu0n0xrEiUTcnMSdQhK0wlmkxEB3g6WNPV5AoOJbjJ6nnopTyV5L7JvcRs2Oc1XcWBeX20vgFAbfd5gaJN74RXlv/7u5Mxt36C5mhZQ0UpmGyqXeTf+FTiXwx3/AZ2JrlYbFH0mXA9SLyHJDB7dZWSqn5Pfonzw3APzXhfZuKZ0u80J5ttpBFoeaZASqlGlae55gOXeGueVd5Z686m7NXnT1Pz2R6EsuwFvXHPNI1wr8/8+9+l64gpPPpsu/nbbXsP7rfvy+ZTy5tUM+6dW710zvfCV/8IvzJn5Q/rhSMjbkltR1LsIjetYupiwYJmaGKDWu2aSNKeNH6F/H9J7/PeHJ8nhVF3IkznhzneOr4vM/S87sSEfZO7mWka8Q3nzuSPMJI14hrb7Jtm5vA3ry467KJzIQ/tGi5yzwjdoRNvZt44ugTdIW7yqqu/G7sFpsAapbGYlcUVwNn4g4wWlJ5rIjcJSKPVLhdW/KcW4E8bjlutfd5p4jsFJGdR47M3wJoFR2hDhRqwUCRK+QQEQ5MHyBqRf1ZEpZhNfRKcCA6UNdQIcd06I/1L/qKb6R7hHwx7xsOmobpV315eIns0qvsJa0oPN7+dnjpS+Ev/9LNW7zvfW7T2datbsJ7zRro7YWXvcx9/Nln63vfZBIee4ypNav8QDcvUMyWyF4xcgUKxX377qv4Vh2hDvZM7JlXAuyttgrFAqOTowx3D/vDcY4kZyfneXmlCy90Gw8XgbeiCNvhljSNdUe6Oa3nNCbS5dbd3rQ4bdsdTOo6kkRko4i8pLQkdrYsdphFrkqUUi9TSm2tcPPmXlyPG4CuUzXqDJVSX1BKXayUurh/gZLB5aTTcb2O0vk0+UL1hjLvhOGNP/VmJjf6isubkdEI++pKjHSNALBnYg8wP6FdVEXfumP/1H4/UT+dnfb9nhaNiLuaUAre+153O+qZZ9ycxY03usnu973PNeB773vd+7duhVtucRvXqvHww1AsMjPQXXXrxjbdsaUDsQG29m/l7j13V/wZHNMhW8z6W4seXj/J2MwY2WKWka4ReiO9CMKRxBFMMU9qgNF0dtpv/GxVd/FgfJDV8dVl5cOZfKZtu7E1C1PvkfQJoFLNX2r2sYYgIlfjJq9fo5SqYi/a3ngnwkwhU9PGI5FL+GaApdYdjfaYsQyL/mj/vKv8RjEYd7dpvNWCZ+XhBaZsIetbd+yb3MeWfndi3FRmChFZuq32aafB00+7wSCRcAPA974Hd9wBf/ZnbnnpY4+5z7njDteJ9mMfg0sugd/+tvJ77nJni0/1xd2GNZl/9WuIgW3aFIoFrtxwJeOpcR49/GjFt+sKdTE6NVq2ovP6SbzPa6RrBMuw6Iv2cSR5BNMwyRfzS548VzrdrpVX78Pdw3SHu8sGWmn/pOBSb6DYoJSaZ9WhlNoJbGignk8DHcCPRORBEangZdDeeIGi2ohFj1Q+RaFY4Hj6eFlpbDP+mFbFVi19m2cBDDEY7hr2VxQiglLK/37eymk6M83x9HFO7zmdiBVhOuMOOjqpUZxDQ+6chVpVNBs3ws03w09+4q46bBve/e7Kg5F274ZVq0iEjZqdzRErQr6Y56LVFxGzY9yz956KzzPEwDKtsi03bxtu7+ReLMPyf/f90XJbkMWM0i1lOjNNxHYbNlvpV2SIwek9p2MZFhPpCeJOXFtjBJh6j6RaG6UNK4pWSm1USq1XSp0/e7uxUe+9XNQbKBLZBEdTR4ETFU9Cc3xm4k4cy7AaPh/ZY7hrmH2T+/xVhCB+6au31bJ/yk1kr+9cT0eow6+5r7s7uyFCh+EjH4Ef/Qi+8Y35j+/eDRdeuKAFRtgKkyvmsE2by4cvZ9fBXVVtLLzE9mR6klwhR76YxzRM9k7uZV3HOn97a65/1GIT/Zl8hsOJw0xmJola0bqcY5uNbdps6tuEQultp4BT75H0axF5x9w7ReRtwK7GSgo2tmnjmA6ZfKbqSdArjT08405J8yqe4OQ8nqphiMFQfKhp208jXSOk8im/r0BE/MFE3laL1zMw0j1CZ6jT30tvVvCqyo03ujOkb74ZjpX4MmUy7vbVhReSzCUJWdXLS6N21L8IuGLDFeSL+Zq2KHEnznPHn/M9wJRS7Jvcx3D3sP+c/mg/E5kJfyVWz4oiW8hyNHGURw8/yoNjD/Lc8efKVqWtDhTgVkKdO3Bu2wwp0iyNeo+k9wL/WUTuFpH/OXu7B3g78J7myQse3vCidD5d9aowW8iCzDcDVKLqsl1YCr2R3uatKGZPeKV5Cu8K29tq2T+5n65QF52hTjqdTqYyU+7W03KuKMCddfGFL8D4OLy/pJfzkUcgl0OVrigq5CjAteb2Vk/rO9ezsXcj9+y5p2rBQMgKkSlk/FLoifQEU5kpvxAAOFEim3DzFOlc9Uq1QrHAU+NP8cDBB3hm4hmKqkhPxG2czBQyRO0oIlJV/3ITtaNtN3tCszjqChRKqUNKqRcDtwF7Zm+3KaVepJSqYsK/MjENk4gVqTk3O1fI+Q1o/bF+bNNumBlgNSJWhLAVbkquYn3negwx/MY7x3SYzkyXbbXsn9zP+q71gFs6Op2ZXrzfU6N4wQvcqqgvfQnumy1v3e16VhXOfwGZfAbHdKqe3Ly5FB7bR7bz/PTz/PZYlSQ57oyQI4kjZY6xw13lKwpwR6xahkUyX72WYzI9ybHUMTc4hHvKDP1SuRRRuz22njSnDos9ki4A7Nnb+SLyPhF5m4icv8DrVgyGuInQTL76ONRMPuP3UJRVPDWxGUlE3O2nJQyaWQjHdFgTX+OfAL3KHc94MF/MMzo96p8Yva0nQ4y6ejyawgc+ACMjbsNeNutWPHV3k16/miLFmp3N3lwKj0vXXUrYCnP33rurfjvTMOmL9hFzYv42XGmg8PbwjyRdF9laJbJjM2NVhwslc0nXegRDW2VoGsZiA8XFwI3A2tnbO4HtwBdF5K8aKy2YeFPuaq0oErkEppgcmjnkDyvKFXNNH+bieRQtBaVUzV6M4e5h9kzuKbvPG5V5cOYg+WKe4U73xNgR6vD9nppVjbUgsZhr2/34466N92wiO1lwT9ALrSiUnPgswlaYbWu3sWN0R80TvPd++yb3MRAd8E/2x1PH3UIGw+Fw4jC2YVcNoOl8munsdFU/sGQu2ZKubM2pzWIDRR9woVLqL5RSf4EbOPqBK4C3NlhbYPGGF0HlmRSpfIrp7DTZYtZfUSjVODPAaoSsEB1Ox6Kv4gvFAuPJcY6ljlUdTDTSNcKx1DG/bl6QEy6pk6PAiStorykxkUssuQy0Ifze77kT4z70IXjoIbjwQn/FVWtFISKEjPIJgts3bCdbyLJjdMeC33bv5F5Guk/kJwQBcVcVR5JH/BLjShca48nxqisFpZQ/3a5d8hOaU4PFBophoPQSMAeMKKVSuN5PGiBmx2pOuUtkE34ppFdHv1wjGAfjg4vq/PV6Pc7oPYOtA1vd/6eOz1tdeIlZb1vFNm1msjMVewY6Q7OBIptoTY6ilDvvhFDI3X668EK/MqySxXgpETtSdiI/ved01nWuq7n9BCdmUJRuO7mOaXNKZGV+iWxRFTk4c7DqyjOdT/vuwzp5rGkkiw0UX8d1kP2AiHwA+DnwDRGJAY81XF1AiTtxf0Uxt9LIK40dm3FrADx7cWhOaexcusJd/hXrQuSLeSbSE2zs3chgfJCOUAdbB7cyGBvkWOpY2baRd4VcWvnkWXfsm9zHmo41/snLO9HNZGdat/XksXq127FtWfDiF/srirAVrrnHH7bCZYFCRNg+sp1njz9b1TkW8GdQeIE1X8xjieWPWPUDhZpfIjudmSZfyFcNAp5HlOdwq9E0ikUFCqXUB4F3ABPAJHCjUurvlFIJpdR1zRAYROJOnEzBTWbPXVF4pbEHZw4Ss2N0hjp9M8DlsFywDIveSK8fyKqRL+aZTE+yqW+TX7rpvX5DzwY2928mk8/4eYi4E6cv0le2olgdX42IMDo16ucnwL2qVkoxk52hqIrL30sxl3e8wy2XHRnxP5eQVXtFEbWj8674XzL8EizD4u49d1d9XWk/CbiFDZ4fV3+kn2Q+6X4uFN1RuSWMzYwRtqtvT3qroXqn22k09bLo+jml1C6l1J1KqU/MWnho5uA1PGUKmXmBwi+NnT7IUPyEGWDUqlzF0gz6o/1Vcw2exqnMFGf1nVW1Uaor3MV5g+fRE+nxzd/mzqawDIupzBTH08f90thCsUBXuAuFak13djU6T+RNYOGtp0qrv7gTZ9vabdyz556K0+/AXXHFnTg9YdelN1vI0hXqImSF6I30Am6JrGM6ZdMJM/kMx1PHq1Y7Af6Wot560jQaXWjdBOJOHKVUxSl32UIWEeFQ4pDfkd0MM8BadIQ6CFkhJtITTKQnOJ46znRmmkQ2QTKXZDozzdl9Z9Mb7a35PrZpM9I14m9jDXcNc2D6QNl20tyegWwhS9xxTfe8xHfLVxQleCuKStPtSrFN201Cz+EPN/8hCHz9ka9XfJ03g6J0WytshYlYEbojs5PzPBfZkkBxLHVswb4IfzVkhvSKQtNQdKBoArWGF3nbLeOpcT9QNMsMsBqGGGwd2MqW/i2c1XcWG3s3MhQfojPUSdgMs3lgs3/SWgjbtH0fqZHuERTK93UC/K9LA0XEitDhdPhOp22xopjF6/1YaOvJG2w09/fbF+3j2rOuZeeBnfzmULmPpj+DoiSRLSKErBARK+KvMg4nDpf1Ung5rYWOkbIVhc5RaBqIDhRNoDPUicItVZx7IknlU/6+/mBs0L9/uevebdMm5sToCnfRF+1jbedazug9g7P7z/arkuqlM9RJtpD1T4D7Jk4kc/dN7KM71O2/Z1EVCVthYk7MtfmQNltRZJP+hLhagcI0TAZiAxVzPa/c+EqG4kP8n9/8n7I8RukMCnA/C1NMHNMhakdxTIeYHfMHGKULaZRSbil1IbvgiFxv2yxqt24WhebURB9NTaAj5K4o0rn5DrLJbJLx1DgAA3E3Sdws19jlwgsU/dF+olbUT9gC7Jva5+cnwP1ZY06MuB33/Z6WPJOiCSRyCUzDxDKsBU+2/bH+iuW9tmnz5vPezNjMGD985of+/f42XHf5Npz3Gpi1G0+WVz4dThyua466F7QWCnIazWLRR1MT8LaeMsVyG498MU9BFTiccF1jPTNAhEDvKUfsCEVVREQY7h72A4U3ItVbaXjGh1E76tp4zPo9tbxEwNVUsgAAIABJREFUtoRELkHYDNfllRSzY4Ttyv5Z5w2ex4WrL+S7T3yXYynXpXbv5F5sw/ZLojP5jL/S8n7/c+3Gk7kk48nxuuxdUrkUlmFhG7ZuuNM0FB0omoBf9TTH78lPZM8cIu7EiTvxppsBLgch84Qp3UjXCPsm97nNYdOz1h2zgcKr7rJNm7gT9/2e2ilQJHNJvwS1UrK6FBFhdXx1Vfv26869jqIq8o1H3NkXeyf3sq7zxAyKoir6VUz+iiLWz9HkUb9A4GjyKAh1HR/JXNLddjIMvaLQNBR9NDUBb0Uxd3iRVxo7NjPm5ydyhVzNkscg4JiO38Q30jVCtuDOivaHFc1uPWULWTpCHZhi0hHqoKAKZPKZtgoUqXyKsOlu3dRzcu4J94CqPJN8IDbANZuuYcfoDh478pg7g6K0Ixt851dvq6s/0k+umGMiPYGIMJGeqLmamMnO8NPnfspH7vsIP93zU3+FogOFppHo0ogm4M0ymBsoSktjz1p1FrA8ZoDNRkTocDrIFrJ+I9meyT3sm9znWnfMbrXkCjnidhzTMP2fOZlPttbvaQ7JbHJRfQi2adMX7WMqM1Xx93jNpmv42b6f8cXdXyybQeFtOZauxkpLZA8nDrOucx2JXGJecUEql2LXwV3sGN3BI4cfoaAKDMWH+P2zf5/Lhy9HKW0xrmksOlA0AcuwCFvheYFiJjuDUopjqWMMxdzS2OUwA1wOOkIdHJo55Ft17JvYx97JvaztWFt20g3bYX9FAa4tRTbSZiuKGmNQKzEYG+Ro4ihUqEdwTIc3nfsm7rz/ToCybbiYFStbtUTsEyWyR5JHOGvVWfNWm0opPnjvB9k/tZ++SB9Xb7yabeu2lfVmHE8dX5Yuf83KQQeKJuBZjafz6bKr5VQ+xWRmEoViMN660thmEHfiPF98HsuwWNexjj2TexidHOXcwXPLnhcyQ4gI3aFuUO6+elutKHJJOkOdi+psjjtxQlaIXCFXsTrpotUXce7AuTxy+JGyfpJVkfI50hEr4lvBH02689Tnbn89fvRx9k/t583nvZmXn/7yqttjekWhaSQ6UDSB0nGopT0CyWzSPwF4gWK5XGObTekWynD3MPeP3k+mkCmz7rAMyz+R9kX6UCimslO+31M7XAWncqma87Ir4Q2F2j+5v2Kjoojw7kvezd6JvUTsCOCuKLxVlUfYDmMYBj3hHr8ybi5377mbmB1j+4btFYOE9zlrNI1EX3Y0AW/KXenWU6XS2EKx4JYytsEJ8mQpHcc50jVCpuB6SXlmgF4i26M32otCMZ1pn+5spRTpQpqItfh5Dr2R3oqW8h5xJ86WgS3+/0VJ2WcGsyNWlZT3UpQwk53h1wd+zYvXv7jqxUVRFXVXtqbh6EDRBETEnZudS1GkiFKKXCHnl8bG7JhfGhv0iicPb7uttEMb5lQ8OScCRcSKELEibdWdrVALzsuuRsgK0RPpWdCVt/R7la7CwE2MK9S8XgqPX+z/Bflinu0bttd83yD35GjaEx0omkTciZPIJfzhRdlCFqUUhxKH/G2n5TYDbDZeh7ZX2dMT7imz7igNio7pEHfc7mxUezjIZvIZCqpAyAwtaZU3FB8inVt4emCukKuYMPeb7qL9HEsdKyuEUErx0z0/5fSe0+eV2JZSVEW99aRpODpQNImYHfPnLhRUwS+NHZsZKzMDPFVWFOD2j+QKOSJ2hNXx1Wzo3lD2eOlWi2M6xJwY05lpFO1h41FqMb6UzuaOUAeWaS24Oiq17ijFsw7pjbjbcuPJcf+xZ48/y+jUKNtHttd870Kx4A5C0mgaiD6imkTMiZHMJ1HKdRgtLY1tpRlgMykt8715283+z+b1DJT+rLZp0+l0MpGZQKHaYuvJc471HHEXiyEGQ/EhDk4fpCvcVfV52UK2qvFixI7QF3FngBxOHPZXn3fvuRvHdNi2bltNDQpVly+URrMY2nJFISL/RUSUiKxa+NntSdyJU1RFcoUcRVUsK431VhRBNwOcS+mKYXXHan/okdczUFqyaYrpbz1ZhtUW3dneGNSQtbStJ3CruWoltQHfnbYSpU13XoVcKpfil6O/ZNu6bX7VVDX01pOmGbRdoBCR9cDLgeqDhwNA3ImDwh9elMql/K2EU8UMcC6WYeGYzrxtpGwhSzxUvtViGm7TnVf11A6Bot6hRbWI2BF6Ij3+6qQiwryKJ4+oHSXuxLEMy698uv95t9R4oW0n0OWxmubQdoEC+DjwV8B885wAEXfiKBSpfIpsIUu+mOdQ4hDg9lDkim5CM8hmgJXoCHXMO+l71h2leCsKhSKVS7VFoKh3DOpCDHcN+yvJuZTOoKhEyAohCL2RXr+U+u49d7O2Yy0bezfW9f11oNA0mrYKFCLyGuB5pdRDrdZysnSF3D3qVD5FOp9GRDicOOyXxp5qiWwPz/OpFBHxHVk9TMP0y2WT+WRbBIpGzXMIW2HWd673Z4KXMrefZC62YSMiDEQHOJI8wv7J/Txz/Bmu3HBlXRcVIqK7sjUNZ9kvPUTkLmCowkO3An8DvKLO93kn8E6A4eHq5YKtonTKXSafOeEaW1oaW8eMgaBRKfgpNb9nwFtRgJtEbodA4eUoGjH4ZzA+6E60K2TLVg+ZfIb+aH/V19mmDcq1G995YCf37L0Hy7C4bP1ldX9vHSg0jWbZjyil1MuUUlvn3oBngdOAh0RkD7AO2C0ilYIKSqkvKKUuVkpd3N9f/Q+vVXj9EdlC1l9RHEoc8iuevJGgpxohK1RmuT3XusNDxN1eQbkn6Hbwe/JmTp/s1hO4K6bTek5jJlOeq5jbTzIX27BRouiP9jOdnea+vfdx0eqLaq5CStHOsZpm0DZHlFLqYaXUgFJqg1JqAzAKXKiUGmuxtCVROpPCa7YbT46fqHg6RTye5uKYDoYY/v58ra2W3vCsjUd2Gmh9d/ZMbgZB/DzBydId7qY70j0vsV0tkQ3ucREyQvRGewF3W65WJ/a81yOnhCWMpr1om0BxquFtq6QLaTKFzAnX2NkVhVKnbr2716EN8607Sin1e2qH7uxENoEhRl3zsutBRNxBTvksRVWsOIOiEhE7Qm/YDRT90X629G+p+fy56BWFptG07RE1u7I42modS6V0HGq2kPVLY4fiQyilEJFTqjS2lNKEdq2tlrAVJmpF3aSv0PLu7GQu6a4mGpgQjtgR1nWtYyoz5Q6psuMLJqUjVoTeSC+WYXHVhqsWfH6+mGcqM8Xx1HGidlRXPWkajj6imoRt2oTNWQfZQp7DyVnX2FO4NNYj5sTKSkOrbbWEzBAdoQ6/a73VW0+JXIKI5Ta0NfKqfCg+xKGZQ8xkZ1gTX7Pg8yN2BMd0+OjLPsqqaOWe00KxQDKXJF/MYxkWQ/EhesI9RO3oKXtcaVqHDhRNotRNNebEylxjk7lk4Mef1sJLaFey7ijFMR1idsw1BqT1W0/JXNKveGrkydYyLE7rPo3fHP5NXUlp7/MaiA1UfDxXyDGdmfa732N2TAcHTVPRgaJJeDMpUvkUXeEuDs2ccI3NF/OnZGmsR8gMgVS27ijFMi06Qh0cSx0DWp/M9raemrF10xPpYbhruK5Kt4VyVzPZGU7rOa1sSqJG00zaNkcRdPwpd7O202OJMT+RXSgWTsnSWA8RIWbFmMnO1LyCNsVtupvOTmMaZs1eitKS22aRyqWWbDG+ECLC6T2nL+jVBLVtXZRSCEJPpKeR8jSamuhA0SS8FUUylyRXyDGeHPcDxalaGltKZ6iTRDZRc4vNNNymu+mMWx5bLVCMTo1yYPpAU3SWks670+1anQy2DAsRqRgck7kkPdGeU/740bQXOlA0kagdJZFLcDhxuMw19lQujfWIh+JuQr/Gyqkev6eJ1AT7JvZxcObggq6sJ0sqn3IHCi1hFkUjEXHdZSs1IWYKmTKbeo1mOdCBoonEQ3FSuVSZGeCpXhrrEbbCxJ14zStfy7D8uQye1UkpmXyG3x77LZ3hTvKFvL/yaBapXOqkLMYbScSKzMvZ5It5HMOp2pei0TQLHSiaSNyOk8qfCBRD8SHyxfwpXRrrETJDdIe7a66cTMP0k/qJXLmNR1EVefb4sxiGgWM6hO2w76baDIqqSDqfXvJ0u0YTtaPz+kpmMjOs7lh9yh87mvZDB4omEnNiZAoZDkwdIGbHiNkxcsXcKekaOxfTMDm95/QFn1fq91S69XRg+gCT6Uk/xxGxIhxPHW+aeWC2kCWv8k2relosETsyL1AolPt5aTTLjA4UTcSbwbBnYg+D8UFE5JQvjV0s3nzo6ew0IkKhWPDzEqWVPSKCEsVEeqIpOjznWMd02iJQeHbjHqmcW2ZdyydKo2kWOlA0kXgojlKK0enRFVMau1g887vJ9CQoN1fx9LGn6Qx3zttiidkxDs4cbEqprDe0qG0CxazduEc6n/aLITSa5UYHiibSYbtJx3wxv6JKYxdDxIoQtaNMZ6dRKJ6beA7TqDwBzjEdUtmUP2CokfhDi8xwWySzbcPGM7AtFAuYhukn/jWa5UYHiibSFe5CzV4Wel20K6E0djE4pkPcifs2HslsbXsT27R9g8XFki1k2Te5r2Kg8baeQtbJz6JoBKZhYuDatc9kZxiKD7WFLs3KpPVr7FOYeOjECW8luMYuhdJA0RXuWvBkGHNijCXGWNu5tu4rf28WyHMTz1FQBcZmxtjcv7ksIHlbT42YbtcowlaYfDFPoVigL9LXajmaFUx7/EWcopSeiAZjgyumNHYxlHZn13OC9oYiecOOFiKdT/PE0Sd4+tjTxJwYvZFeonaURw8/ykTqRGK8NJndiKFFjSBqR30blHqsPzSaZqEDRRPxAoXnGrtSSmMXg+f3NJWdqvs1ESvC2EztwYdFVeTg9EEeGnuIVD5FX7TPT1I7pkPMifHE+BP+Npafo2ijFUXEjpDIJVgdX91qKZoVjt56aiKdjpt8HIgN6NLYKpiG6c+k8JK2CxGxIxxLHiOTz1QsF01kEzxz/BlSuRSdoc6K7+mYDp2hTp48+iRn9J7hB4p2yVGAG7S6nC66wl2tlqJZ4bTHX8Qpite85ZU16tLY+ViG5VtSzJ0tXQtDjHk9FYVigdGpUR4+/DBKKXoiPTUDj2VY9ER6eObYMxxKHEIQIlakrQLFcNdwW1RhaVY2ekXRRAwxuGrDVVy85mJAl8ZWwjMGBPyEdj3EnBgHpw/6q7VENsFvj/+WTC5Dd7i77pO9aZj0Rns5njru/37aJVBE7IjOTWjaAh0omoghBq8957V0h7vdO9TCQ2lWGqZh+iuKehPU4H6OM9kZpjJTTGen2T+5n6gdpTvSvWgNhhi+71Qj52VrNKcK+i+iiZSecJRSKFG6NLYCPZEelFJMZibnPVYoFth5YOc83yNwg8UTR5/g+ann6Yn0nNTVt2cxrlA6UGg0c9B/EU3EEMO3m8gX80SsiC6NrUBv2LXxqGQj/u3Hvs2d99/Jj5790bzH4k6cmBOjJ9Jz0if3TD5DyAxh0Nh52RrNqYAOFE3E28YoqiK5Yo6IpfebK9EXdZvJvO5sjwfHHuQHT/8AQwzu3XNvRY+nRuV80vl028yi0GjaDR0omoxjOhRVUZfG1iBkhYg5sbIVxXhynM/t/BzDXcO8aeubGJ0e5dmJZ5umIZ1Pt8V0O42mHdGBosmYYlJURQrFgq5gqYJv4zHbdJcv5vnMrz9DQRW46YU3cfnI5TiGw31772uahnYaWqTRtBs6UDQZy7BOeDzpiqeKOKZD3I77K4p/fuyfefrY09xwwQ0MxYeI2lEuWnMRO0Z3NG1wkZfMbgeLcY2m3Wi7QCEiN4nIkyLyqIjc3mo9J4tpuCsKVOP20081PL+nqcwUDxx8gB88/QNeetpLedG6F/nPuXLkShK5BLsO7mqKhnQ+jWM6Okeh0VSgrQKFiFwFXAucp5TaAnysxZJOGsuwKKiCLo2tgWVYdIY6GU+N8/ldn3fzEue+qew5m/s3syq6inv33tsUDV6g0CsKjWY+bRUogHcB/10plQFQSh1usZ6TxjEdsoWsLo2tgSmu31O2kPXzEnNXXyLCZesv49HDjy55HkWhWPBdYkspqiLZQrZt5mVrNO1GuwWKTcDlInK/iNwjIpdUe6KIvFNEdorIziNHjiyjxMVhGZYfKDSVMQ2T7pDbUe3lJSpx+cjlKBQ/2/+zRX+Pp8af4taf3Mp7f/heHj38aNlj6XwagJCpA4VGU4ll/6sQkbuASmeCW3H19ADbgEuAb4nI6apCAb1S6gvAFwAuvvjixg9RbhBeoNClsdUxxeSStZeweWAzG7o3VH3eQGyAc1adw7177+U1m15T1wotmUvyrUe/xY+f+zF9kT76In187Jcf408v+VPfg8sPFHpFodFUZNn/KpRSL6v2mIi8C/i/s4HhVyJSBFYB7btkWABDDNeVVJfGVsU0TGzTZiQ6suBzrxi5gs/v+jxPjj/J2avOrvncnQd28tWHvspEeoKrz7ia121+Hflinjt+eQefvP+T3HDBDWzfsN0PFLZp62S2RlOBdtt6+i7wOwAisglwgKMtVXSSGGLgWI4ujV0Ax3AoqMKCz7tkzSVErEjNpPbx1HE+ef8nufP+O4k7cT6w/QNcd951hK0wcSfO+1/yfs4dPJcvP/BlfvDUD0jlUwCEzXDbTLfTaNqJdltnfwX4iog8AmSB6yttOwUJQwxCZkiXxi6AY7kd7AsRskK8cO0L2TG6g7e84C1l8z2UUvxs38/42sNfI1fI8YbNb+CVZ75y3nZSyApx87ab+cKuL/DNR7/J5lWbgfaabqfRtBNtFSiUUlngP7VaRyMxDRPHcHRp7ALYhu1Omatj5+eKkSu4Z+893P/8/Vw5ciUAx1LH+MoDX+GhQw9xVt9ZvP3Ct1dNioObO3rXxe8i5sS469m7AB0oNJpqtFWgOBUxxCDmxHRp7AI4plP3PIoze89kKD7EfXvv44rhK7h377187eGvUVAF3nzem3n56S+v6/MWEd5y3lvocDr416f/tSEutBrNqYgOFE0mZIZY17mu1TLaHsd0KBQXzlGAe4K/YvgKvvXYt/jwfR/myfEnOWfVObztgrcxGB9c1PcVEV57zmu59qxrmc5O60Ch0VRA/1U0GRFd8VQPIStUV47C47LhyzDEYM/EHq5/wfXcctktiw4SpXjVTjpQaDTz0SsKTVtgill1u6hQLKBQZUnpnkgPt15+q9sbMTvPohHoQKHRzEcHCk1bYBomVKhvU0oxkZ5AEHqjvWWPberb1FANSukxqBpNJfRfhaYtMMVEVYgUE5kJVsdX45hOxbnZjUYHCo1mPvqvQtMWVOqITmQTxOwY67vWMxQfqmjo1yi8dh1dnabRzEcHCk1bMHeyXK6QI1fMsbF3I6Zh0hvtratze6nMzYFoNJoT6EChaQtMw8QQA6UUSimmMlNs7Nnod16HrTCdTqfvy9RoiqqoA4VGUwUdKDRtg2O6fk8TmQnWdKyZl7weig+RyqWa8r2LqqjnZWs0VdCBQtM22KbNdGaamB2r2KTYGe4EWFS/Rb0opbeeNJpq6EChaRs8PywvLzEXy7AYjA+6nlANpqiK2mJco6mCvoTStA1d4S5WRVeVOcLOZVV0FWMzY0t6/5nsDJl8Bsuw6Ap3lT2mcxQaTXX0ikLTNgzEBuiJ9NR8TsyOETbD5Aq5ut83V8gxnhwnZsc4b/A8OkOdHE8fL3uOQukchUZTBR0oNIFCRFjdsZqZ7MyCz1VKMZWeIpVLsalvE5v6NhFzYmzs3UhvuJdjyWN+/4TeetJoqqPX2prA0R3uZqF5VtlClqnMFEPxIdZ3ri+bMGgaJmf0noFlWByaOURPpEcnszWaGui/DE3gCFkheiI9JHNJona07DGlFJOZSWzDZuvAVjpDnRXfwxCDDd0bsAyL/VP7McTQgUKjqYLeetIEksH44Lzmu0w+w7HUMQZjg5w7eG7VIOEhIqzvWs9p3aeRLWS1z5NGUwV9CaUJJB1OB6aYFIoFDDGYyEwQMkJsHdhKR6hjUe+1umM1ITOk54ZoNFXQgUITSEzDZDA2yOj0KIKwtmMtazvXLjkhPbcLXKPRnEAHCk1gWRVbxVRmipHukUWvIjQaTf3oQKEJLFE7ytbBra2WodGc8ujsnUaj0WhqogOFRqPRaGqiA4VGo9FoatJWgUJEzheRHSLyoIjsFJEXtlqTRqPRrHTaKlAAtwO3KaXOB/7r7P81Go1G00LaLVAowGun7QIOtFCLRqPRaGi/8tj3Aj8UkY/hBrEXt1iPRqPRrHiWPVCIyF3AUIWHbgVeCtyslPqOiLwB+DLwsirv807gnQDDw8NNUqvRaDQaWciueTkRkUmgWymlRESASaVUbWc393VHgL1L/LargKNLfO1yEyStECy9QdIKwdIbJK0QLL0nq3VEKdW/0JPabevpAHAlcDfwO8DT9byonh+0GiKyUyl18VJfv5wESSsES2+QtEKw9AZJKwRL73JpbbdA8Q7gThGxgDSzW0sajUajaR1tFSiUUj8DLmq1Do1Go9GcoN3KY1vBF1otYBEESSsES2+QtEKw9AZJKwRL77Jobatktkaj0WjaD72i0Gg0Gk1NdKDQaDQaTU10oNBoNBpNTXSg0GgChIj0iIie+7rCWe7jQAeKWUTk31qtoRQRuVhEfioi/ygi60XkRyIyKSK/FpELWq1vLiJyQ8nX60TkxyIyISK/EJFNrdQ2l9nP85sicp+I/I2I2CWPfbeV2iohImtE5KuzzgVHgUdFZJ+I/G2p9nZARK4u+bpLRL4sIr8Rka+LyGArtVUiSMdCK4+DFRUoROTCKreLgPNbrW8Of49rs/4D4BfA55VSXcBfzz7WbvxZydd3AN8CeoH/AXy2JYqq8xXc7v+bgNXAPSLSN/vYSKtE1eAfga/M/v5fD3wHOAe3D+ozrRRWgf9W8vX/BA4CrwZ+DXy+JYpqE6RjoWXHwYoqjxWRAnAPIBUe3qaUiiyzpKqIyANKqQtmv96nlBqu9Fi7ICK7lVIXzn794OxMEe+xttJbQd9/Am4BXgN82/s52gUReUgp9YKS/+9SSl00+/UTSqmzW6eunAWOg7L/twNBOhZaeRy0VWf2MvA48CdKqXkeUiKyvwV6apEWkVfgzuVQIvL7SqnvisiVQKHF2iqxTkQ+iRuE+0XEVkrlZh9rq+0RwBaRsFIqDaCU+kcRGQN+CMRaK60iR2ZPYD8BXgfsAZg1zmy3XYEBEXkf7nHQKSKiTlyNtptWCNax0LLjoB1/cc3kb6n+M9+0jDrq4UbgL4AbgN8FrhKRCdxtpz9vpbAq/CWwC9gJ/A0QBxCRIeD7LdRViS8Bl5beoZS6C3c5/0hLFNXmBtwr3P/A1e1t8/XiXv22E18EOnB///8b193UOw4ebKGuagTpWGjZcbCitp40Go1Gs3hW2ooCEfldEfmsiHxfRL43+/XVC79y+ami9XdbrWuxiMh/bbWGucx+tm8TkQ1z7r+h8itaS4nekTn3t53eIP2N1aIdj9tqNFvrilpRiMgngE3AV4HR2bvXAW8BnlZKvadV2uYSJK0LMTcZ32pE5CPAS4DduBU5n1BKfWr2sd3tlMAEEJH/BlxGAPTq47Y1NFvrSgsUTyml5tX0zyaDnlJKndkCWRUJklYAEZmq9hAQUUq1TeGEiDwMXKCUyotIN/B14Eml1M3tVqEFwdKrj9vm0UqtK23rKS0iL6xw/yW4g5LaiSBpBZgAzlRKdc65deDW0rcTllIqD6CUmsC9Su8UkW8DTkuVVSZIevVx2zxaprVtouUy8Vbgs+K2vnvL4vXA1Oxj7cRbCY5WcLcaRoBDFR77+jJrWYhnRORKpdQ9AEqpAvA2EfkQbtlhuxEkvW9FH7fNomVaV9TWk8dsqd5a3CXbqFJqrMWSqhIkrUFBRCIASqlUhcfWKqWeX35V1QmaXtDH7anGStt6AkApNaaU2qWU2onbr9C2BEnrXETkb1utoRJKqdTck66ntR1PukHTC/q4XS6WS+uKDBRzeE2rBSyCIGmFYOkNklYIlt4gaYVg6V0WrTpQVPZ9aleCpBWCpTdIWiFYeoOkFYKld1m0rsgcRSkiYiiliq3WUQ9B0grB0hskrRAsvUHSCsHSu1xaV/yKwvuQg9CFGSStECy9QdIK7a13btd7ida26yKHYOltldYVHyhKeHurBSyCIGmFYOkNklZoM72zXeS3AucCPxaRUrPNP6v8qtYRJL2t1Lqitp50F2bzCJLeIGmFYOkNUhc5BEtvK7WutBWF7sJsHkHSGyStECy9Qeoih2DpbZnWlRYovM7GSrRrF2Yl2k0rBEtvkLRCsPQ+I+5wLcDtIldKvQ14EndsZ7sRJL0t07qitp40Gk1zCVoXeZD0tlLrigsUItIFXI1rL6CAA8APZ5dybUWQtEKw9AZJKwRLb5C0QrD0tkrritp6EpG34Hr6bweiuDNxrwJ2zT7WNgRJKwRLb5C0QrD0BkkrBEtvK7WuqBWFiDwJXDo3+opID3C/quCj3yqCpBWCpTdIWiFYeoOkFYKlt5VaV9SKArecsFJkLNJ+bftB0grB0hskrRAsvUHSCsHS2zKtbVN/vUx8GNgtIv8B7J+9bxh4OfDBlqmqTJC0QrD0BkkrBEtvkLRCsPS2TOuK2noCf5n2u5R45eMmg463VFgFgqQVgqU3SFohWHqDpBWCpbdVWldcoJiLiFyjlPp/rdZRD0HSCsHSGyStECy9QdIKwdK7XFp1oBDZrZS6sNU66iFIWiFYeoOkFYKlN0haIVh6l0vrSktmV6LdEla1CJJWCJbeIGmFYOkNklYIlt5l0aoDBfxJqwUsgiBphWDpDZJWCJbeIGmFYOldFq0rLlCISKeInOH9Xyn1q9n7z2udqsoESSsES2+QtEKw9AZJKwRLb6u0rqhAISJvAJ4AviMbTmNKAAADqklEQVQij4rIJSUP/6/WqKpMkLRCsPQGSSsES2+QtEKw9LZUq1JqxdyAB4HVs1+/cPZDf+3s/x9otb6gag2a3iBpDZreIGkNmt5Wal1pDXemUuoguEs2EbkK+H8iso7KHY+tJEhaIVh6g6QVgqU3SFohWHpbpnVFbT0B03P29w7iGmxdC2xplagqBEkrBEtvkLRCsPQGSSsES2/LtK60FcW7mBMclVLTInI18IbWSKpKkLRCsPQGSSsES2+QtEKw9LZM64psuBORQUr83JVSh1osqSpB0grB0hskrRAsvUHSCsHS2wqtKypQiMj5wOeALsCbBrUOdybxu5RSD7RK21yCpBWCpTdIWiFYeoOkFYKlt6VaW53JX84bbtXApRXu3wY81Gp9QdUaNL1B0ho0vUHSGjS9rdS60pLZMaXU/XPvVErtwJ0W1U4ESSsES2+QtEKw9AZJKwRLb8u0rrRk9r+JyA+Ar3LCz3098Bbg31umqjJB0grB0hskrRAsvUHSCsHS2zKtKypHASAir8QtJyv1c/++UupfWyqsAkHSCsHSGyStECy9QdIKwdLbKq0rLlBoNBqNZnGsqByFiHSJyH8XkcdFZHz29vjsfd2t1ldKkLRCsPQGSSsES2+QtEKw9LZS64oKFMC3gOPAVUqpPqVUH3AVbnnZt1uqbD5B0grB0hskrRAsvUHSCsHS2zKtK2rrSUSeVEqdtdjHWkGQtEKw9AZJKwRLb5C0QrD0tlLrSltR7BWRv5rtbATcLkcReT8nqgjahSBphWDpDZJWCJbeIGmFYOltmdaVFij+COgD7hGR4yJyDLgb6KX9fF2CpBWCpTdIWiFYeoOkFYKlt2VaV9TWE4CInI3b9r5DKTVTcv/VSqm2qpsOklYIlt4gaYVg6Q2SVgiW3pZpXY7W83a5AX8OPAl8F9gDXFvy2O5W6wuq1qDpDZLWoOkNktag6W2l1pb/8Mv8QT8MxGe/3gDsBN4z+/92m2YVGK1B0xskrUHTGyStQdPbSq0rzcLDVLPLNaXUHhHZDvyziIzgdjm2E0HSCsHSGyStECy9QdIKwdLbMq0rLZk9Jq5VLwCzH/o1wCrg3JapqkyQtEKw9AZJKwRLb5C0QrD0tkzrikpmiztbNq+UGqvw2EuUUj9vgayKBEkrBEtvkLRCsPQGSSsES28rta6oQKHRaDSaxbPStp40Go1Gs0h0oNBoNBpNTXSg0Gg0Gk1NdKDQaDQaTU10oNBoNBpNTf5/0cot0ETTdEgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt \n",
    "import pandas as pd \n",
    "%matplotlib inline \n",
    "\n",
    "DF = pd.DataFrame()\n",
    "DF['obs'] = obs\n",
    "DF['sim'] = sim\n",
    "DF.index = pd.date_range('2002-01-01', periods = len(DF), freq='D')\n",
    "\n",
    "plt.plot(DF.sim,c='red',label='simulated')\n",
    "plt.plot(DF.obs,c='green',label='observed',alpha=0.75)\n",
    "plt.fill_between(DF.index, DF.obs + 1 , DF.obs - 1 , \n",
    "                color='green', alpha=0.15)\n",
    "plt.ylabel('gCm$^{-2}$')\n",
    "plt.xticks(rotation=90)\n",
    "plt.legend(loc='best')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can right click on the figure above to save it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Call the compare_daily() function specifying that the uncertainty is provided as a quantity (i.e. \"qt\") and is equal to 1 gCm-2 for all the observed data points. The value returned expresses the fraction of modelled data points that were within a +/- 1 gCm-2 range around the observed data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.36\n"
     ]
    }
   ],
   "source": [
    "print (compare_daily(sim,obs,'qt',1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fraction of modelled data points that were within a +/- 1 gCm-2 range around the observed data points when considering a +/- 3 days/timesteps window can be returned by calling compare() instead of compare_daily(). The time lag window can be reduced or increased by editing lines 28-46 of the compare() and lines 65-83 of the compare_daily() functions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compare(sim,obs,'qt',1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The lines below repeat the two function calls performed just above this time using a percentage (\"pc\") instead of a quantity (\"qt\") in order to express observational uncertainty. Both with and without time-lag consideration the observational uncertainty is set equal to 25% of the observed value (i.e. 0.25 expressed as a fraction)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.08\n"
     ]
    }
   ],
   "source": [
    "print (compare_daily(sim, obs, 'pc', 0.25)) # no consideration of timelag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.32\n"
     ]
    }
   ],
   "source": [
    "print (compare(sim, obs, 'pc', 0.25)) # +/- 3 days timelag considered"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can repeat all the above comparisons using your own data by saving them as a .csv using \"sim\" as a header for the simulated data and \"obs\" as a header for the measured data and then running the lines below to load the file. Make sure you change the \"/Users/vm/Desktop/laiDF.csv\" into the location of your .csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "ename": "IOError",
     "evalue": "[Errno 2] File /Users/vm/Desktop/gcel tutorial/neeDF.csv does not exist: '/Users/vm/Desktop/gcel tutorial/neeDF.csv'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIOError\u001b[0m                                   Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-8-fc06c9a147b5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mDF\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mread_csv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"/Users/vm/Desktop/gcel tutorial/neeDF.csv\"\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# will load your data and store them as a dataframe called DF\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/Users/vm/miniconda2/lib/python2.7/site-packages/pandas/io/parsers.pyc\u001b[0m in \u001b[0;36mparser_f\u001b[0;34m(filepath_or_buffer, sep, delimiter, header, names, index_col, usecols, squeeze, prefix, mangle_dupe_cols, dtype, engine, converters, true_values, false_values, skipinitialspace, skiprows, skipfooter, nrows, na_values, keep_default_na, na_filter, verbose, skip_blank_lines, parse_dates, infer_datetime_format, keep_date_col, date_parser, dayfirst, iterator, chunksize, compression, thousands, decimal, lineterminator, quotechar, quoting, doublequote, escapechar, comment, encoding, dialect, tupleize_cols, error_bad_lines, warn_bad_lines, delim_whitespace, low_memory, memory_map, float_precision)\u001b[0m\n\u001b[1;32m    700\u001b[0m                     skip_blank_lines=skip_blank_lines)\n\u001b[1;32m    701\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 702\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0m_read\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilepath_or_buffer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    703\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    704\u001b[0m     \u001b[0mparser_f\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__name__\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/vm/miniconda2/lib/python2.7/site-packages/pandas/io/parsers.pyc\u001b[0m in \u001b[0;36m_read\u001b[0;34m(filepath_or_buffer, kwds)\u001b[0m\n\u001b[1;32m    427\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    428\u001b[0m     \u001b[0;31m# Create the parser.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 429\u001b[0;31m     \u001b[0mparser\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mTextFileReader\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilepath_or_buffer\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    430\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    431\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mchunksize\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0miterator\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/vm/miniconda2/lib/python2.7/site-packages/pandas/io/parsers.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, f, engine, **kwds)\u001b[0m\n\u001b[1;32m    893\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'has_index_names'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwds\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'has_index_names'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    894\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 895\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_make_engine\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mengine\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    896\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    897\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/vm/miniconda2/lib/python2.7/site-packages/pandas/io/parsers.pyc\u001b[0m in \u001b[0;36m_make_engine\u001b[0;34m(self, engine)\u001b[0m\n\u001b[1;32m   1120\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_make_engine\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mengine\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'c'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1121\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mengine\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'c'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1122\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mCParserWrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptions\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1123\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1124\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mengine\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'python'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/vm/miniconda2/lib/python2.7/site-packages/pandas/io/parsers.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, src, **kwds)\u001b[0m\n\u001b[1;32m   1851\u001b[0m         \u001b[0mkwds\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'usecols'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0musecols\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1852\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1853\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_reader\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mparsers\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mTextReader\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1854\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munnamed_cols\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_reader\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munnamed_cols\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1855\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/parsers.pyx\u001b[0m in \u001b[0;36mpandas._libs.parsers.TextReader.__cinit__\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/parsers.pyx\u001b[0m in \u001b[0;36mpandas._libs.parsers.TextReader._setup_parser_source\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mIOError\u001b[0m: [Errno 2] File /Users/vm/Desktop/gcel tutorial/neeDF.csv does not exist: '/Users/vm/Desktop/gcel tutorial/neeDF.csv'"
     ]
    }
   ],
   "source": [
    "DF = pd.read_csv(\"/Users/vm/Desktop/gcel tutorial/neeDF.csv\") # will load your data and store them as a dataframe called DF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DF.sim # will print your modelled data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DF.obs # will print your observed data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can add a date index to your dataframe (i.e. DF) so that you have temporal information related to your data points. In the lines below a daily date index, with starting date being 2002-01-01, is added. The obs and sim data are then ploted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "DF.index = pd.date_range('2002-01-01', periods = len(DF), freq='D')\n",
    "\n",
    "# plot sim and obs data\n",
    "plt.plot(DF.sim,c='red',label='simulated')\n",
    "plt.plot(DF.obs,c='green',label='observed',alpha=0.5)\n",
    "plt.ylabel('gCm$^{-2}$')\n",
    "plt.title('NEE')\n",
    "plt.xticks(rotation=90)\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(compare_daily(DF.sim, DF.obs, 'pc', 0.25)) # compare obs and sim without considering +/- 3 days timelag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(compare(DF.sim, DF.obs, 'pc', 0.25)) # compare obs and sim while considering +/- 3 days timelag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (\"2002 : \" + str(compare(DF.sim['2002'], DF.obs['2002'], 'pc', 0.25)))\n",
    "print (\"2003 : \" + str(compare(DF.sim['2003'], DF.obs['2003'], 'pc', 0.25)))\n",
    "print (\"2004 : \" + str(compare(DF.sim['2004'], DF.obs['2004'], 'pc', 0.25)))\n",
    "print (\"2005 : \" + str(compare(DF.sim['2005'], DF.obs['2005'], 'pc', 0.25)))\n",
    "print (\"2006 : \" + str(compare(DF.sim['2006'], DF.obs['2006'], 'pc', 0.25)))\n",
    "print (\"2007 : \" + str(compare(DF.sim['2007'], DF.obs['2007'], 'pc', 0.25)))\n",
    "print (\"2008 : \" + str(compare(DF.sim['2008'], DF.obs['2008'], 'pc', 0.25)))\n",
    "print (\"2009 : \" + str(compare(DF.sim['2009'], DF.obs['2009'], 'pc', 0.25)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DF = DF.resample('W').mean()\n",
    "print (\"2002 : \" + str(compare(DF.sim['2002'], DF.obs['2002'], 'pc', 0.25)))\n",
    "print (\"2003 : \" + str(compare(DF.sim['2003'], DF.obs['2003'], 'pc', 0.25)))\n",
    "print (\"2004 : \" + str(compare(DF.sim['2004'], DF.obs['2004'], 'pc', 0.25)))\n",
    "print (\"2005 : \" + str(compare(DF.sim['2005'], DF.obs['2005'], 'pc', 0.25)))\n",
    "print (\"2006 : \" + str(compare(DF.sim['2006'], DF.obs['2006'], 'pc', 0.25)))\n",
    "print (\"2007 : \" + str(compare(DF.sim['2007'], DF.obs['2007'], 'pc', 0.25)))\n",
    "print (\"2008 : \" + str(compare(DF.sim['2008'], DF.obs['2008'], 'pc', 0.25)))\n",
    "print (\"2009 : \" + str(compare(DF.sim['2009'], DF.obs['2009'], 'pc', 0.25)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DF = DF.resample('W').mean()\n",
    "print (\"2002 : \" + str(compare(DF.sim['2002'], DF.obs['2002'], 'qt', 1)))\n",
    "print (\"2003 : \" + str(compare(DF.sim['2003'], DF.obs['2003'], 'qt', 1)))\n",
    "print (\"2004 : \" + str(compare(DF.sim['2004'], DF.obs['2004'], 'qt', 1)))\n",
    "print (\"2005 : \" + str(compare(DF.sim['2005'], DF.obs['2005'], 'qt', 1)))\n",
    "print (\"2006 : \" + str(compare(DF.sim['2006'], DF.obs['2006'], 'qt', 1)))\n",
    "print (\"2007 : \" + str(compare(DF.sim['2007'], DF.obs['2007'], 'qt', 1)))\n",
    "print (\"2008 : \" + str(compare(DF.sim['2008'], DF.obs['2008'], 'qt', 1)))\n",
    "print (\"2009 : \" + str(compare(DF.sim['2009'], DF.obs['2009'], 'qt', 1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}