{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic statistics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find out where we are and save data files\n",
    "- use os.getcwd() to know your current path\n",
    "- create a data folder at this location, for example *data\\*\n",
    "- unzip the data file downloaded at:\n",
    "https://perso.ensta-paris.fr/~monchaux/doc_M1IPP_102/data_in_zip.zip\n",
    "- by now, you should have in the data folder several files:\n",
    "    - **meteo_orly.csv**\n",
    "    - **signal_cylindre.txt**\n",
    "    - **covid-19-pandemic-worldwide-data.csv**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Directory already exist\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import csv\n",
    "\n",
    "os.getcwd()\n",
    "\n",
    "if os.path.exists(\"data\\\\\"):\n",
    "    print('Directory already exist')\n",
    "    print()\n",
    "else:\n",
    "    os.mkdir(\"data\\\\\")\n",
    "    print('Directory created')\n",
    "    print()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- read the data file **meteo_orly.csv** and print the first lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-2-806c48438e08>, line 4)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-2-806c48438e08>\"\u001b[1;36m, line \u001b[1;32m4\u001b[0m\n\u001b[1;33m    nblines =\u001b[0m\n\u001b[1;37m              ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "with open(\"data\\\\meteo_orly.csv\", \"r\") as csv_file:\n",
    "    csv_reader = csv.reader(csv_file, delimiter = ';')\n",
    "    linecount = 0\n",
    "    nblines = \n",
    "    print(f\"Overview of {nblines} first lines:\")\n",
    "    for lines in csv_reader:\n",
    "      if linecount <= nblines:\n",
    "        print(lines)\n",
    "        linecount += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading all data\n",
    "- choose a data column to work with (temperature, wind direction or wind speed)\n",
    "- isolate the desired column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data\\\\meteo_orly.csv\", \"r\") as csv_file:\n",
    "    csv_reader = csv.reader(csv_file, delimiter = ';')\n",
    "    # create empty list for data\n",
    "    timedata = []\n",
    "    ?? = []\n",
    "    for cnt,lines in enumerate(csv_reader):\n",
    "        # do not take header into account\n",
    "        if cnt > 0:\n",
    "            timedata.append(lines[0])\n",
    "            tempdata.append(float(lines[??]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Average\n",
    "We generate our own algorithm to estimate the average of a list and test it against the built-in numpy function:\n",
    "- write a for loop to calculate the average of the whole dataset\n",
    "- compare the output value to the built-in output\n",
    "\n",
    "Change your program to do the same only for the N first values\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-4-16250ee1c854>, line 5)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-4-16250ee1c854>\"\u001b[1;36m, line \u001b[1;32m5\u001b[0m\n\u001b[1;33m    for :\u001b[0m\n\u001b[1;37m        ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# calculate average with a for loop\n",
    "temp_average_ih = 0\n",
    "for :\n",
    "    \n",
    "\n",
    "\n",
    "# calculate average with built-in python function\n",
    "temp_average_bi = np.mean(tempdata)\n",
    "\n",
    "print(f\"In house calculated average temperature is {temp_average_ih}\")\n",
    "print(f\"Built-in calculated average temperature is {temp_average_bi}\")\n",
    "\n",
    "\n",
    "# calculate average with a restricted for loop\n",
    "lendata = 10\n",
    "\n",
    "temp_average_ih = 0\n",
    "for :\n",
    "    \n",
    "\n",
    "# calculate average with built-in python function\n",
    "temp_average_bi = np.mean(tempdata[0:lendata])\n",
    "\n",
    "print(f\"In house calculated average temperature over {lendata} data points is {temp_average_ih}\")\n",
    "print(f\"Built-in calculated average temperature over {lendata} data points is {temp_average_bi}\")\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Conclusion:\n",
    "Our calculation is equal to the built-in function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Standard deviation\n",
    "Do the same work for standard deviation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-5-324c38fdd394>, line 5)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  File \u001b[1;32m\"<ipython-input-5-324c38fdd394>\"\u001b[1;36m, line \u001b[1;32m5\u001b[0m\n\u001b[1;33m    for :\u001b[0m\n\u001b[1;37m        ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "# calculate standard with a for loop\n",
    "temp_average2_ih = 0\n",
    "temp_average_ih = 0\n",
    "\n",
    "for :\n",
    "    \n",
    "\n",
    "# calculate standard deviation with built-in python function\n",
    "temp_std_bi = np.std(tempdata, ddof=1)\n",
    "\n",
    "print(f\"In house calculated standard deviation of temperature is {temp_std_ih}\")\n",
    "print(f\"Built-in calculated standard deviation of temperature is {temp_std_bi}\")\n",
    "\n",
    "# calculate standard with a for loop rstricted to lendata points\n",
    "lendata = 10\n",
    "temp_average2_ih = 0\n",
    "temp_average_ih = 0\n",
    "\n",
    "for :\n",
    "    \n",
    "\n",
    "# calculate standard deviation with built-in python function\n",
    "temp_std_bi = np.std(tempdata[0:lendata], ddof=1)\n",
    "\n",
    "print(f\"In house calculated standard deviation of temperature over {lendata} data points is {temp_std_ih}\")\n",
    "print(f\"Built-in calculated standard deviation of temperature over {lendata} data points is {temp_std_bi}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Work with second file\n",
    "We now work with a hot-wire signal acquired in the wake of a cylinder in a wind-tunnel\n",
    "\n",
    "The file to load is **signal_cylindre.txt**\n",
    "\n",
    "From now on we work only with built-in functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data\\\\signal_cylindre.txt\", \"r\") as csv_file:\n",
    "    csv_reader = csv.reader(csv_file, delimiter = ';')\n",
    "    # create empty list for data\n",
    "    speeddata = []\n",
    "    for cnt,lines in enumerate(csv_reader):\n",
    "        # do not take header into account\n",
    "        if cnt > 0:\n",
    "            speeddata.append(float(lines[0]))\n",
    "\n",
    "# calculate average\n",
    "speeddata_average = \n",
    "# calculate standard deviation\n",
    "speeddata_std = \n",
    "\n",
    "print(f\"Built-in calculated average velocity is {speeddata_average}\")\n",
    "print(f\"Built-in calculated standard deviation of velocity is {speeddata_std}\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting data and statistics\n",
    "- import matplotlib and plot the data as time series\n",
    "- superimpose the statistical values on the plot as horizontal lines\n",
    "- zoom on the x axis by setting the axis limits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 4.5), tight_layout=True)\n",
    "\n",
    "# data points\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "ax1.plot(, '+')\n",
    "ax1.set_ylim([0, 16])\n",
    "\n",
    "# statistics as lines\n",
    "ax1.plot([0, len(speeddata)],[speeddata_average, speeddata_average],'r')\n",
    "ax1.plot([0, len(speeddata)],[XX, XX], 'k')\n",
    "ax1.plot([0, len(speeddata)],[XX, XX], 'k')\n",
    "\n",
    "# labels\n",
    "ax1.set_xlabel('Time (AU)')\n",
    "ax1.set_ylabel('U (m/s)')\n",
    "\n",
    "# Same zoomed in\n",
    "# data points\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "ax2.plot(, '+')\n",
    "ax2.set_ylim([0, 16])\n",
    "\n",
    "# statistics as lines\n",
    "ax2.plot([0, len(speeddata)],[speeddata_average, speeddata_average],'r')\n",
    "ax2.plot([0, len(speeddata)],[XX, XX], 'k')\n",
    "ax2.plot([0, len(speeddata)],[XX, XX], 'k')\n",
    "\n",
    "# labels\n",
    "ax2.set_xlabel('Time (AU)')\n",
    "ax2.set_ylabel('U (m/s)')\n",
    "\n",
    "# zoom\n",
    "ax2.set_xlim([XX, XX])\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convergence of the estimators\n",
    "- build a loop to estimate the average and the standard deviation on an increasing number of successive points\n",
    "- build a loop to estimate the average and the standard deviation on an increasing number of non successive (distant) points\n",
    "- in both cases, plot the estimated values as a function the number of data points used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# work with consecutive data points\n",
    "ave_list1 = []\n",
    "std_list1 = []\n",
    "nbpt_list = [10, 30, 30000]\n",
    "for nb in nbpt_list:\n",
    "    ave_list1.append(XX)\n",
    "    std_list1.append(XX)\n",
    "\n",
    "# plot results\n",
    "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.5), tight_layout=True)\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "ax1.plot(nbpt_list, ave_list1, 's-')\n",
    "ax1.set_xlabel('Number of points')\n",
    "ax1.set_ylabel('Estimated average')\n",
    "ax1.set_xscale('log')\n",
    "\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "ax2.plot(nbpt_list, std_list1, 'o-')\n",
    "ax2.set_xlabel('Number of points')\n",
    "ax2.set_ylabel('Estimated standard deviation')\n",
    "ax2.set_xscale('log')\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# work with distant data points\n",
    "ave_list2 = []\n",
    "std_list2 = []\n",
    "nbpt_list = [10, 30, 100, 300, 1000, 3000, 10000, 30000]\n",
    "for nb in nbpt_list:\n",
    "    data = np.random.choice(speeddata,XX)\n",
    "    ave_list2.append(XX)\n",
    "    std_list2.append(XX)\n",
    "    \n",
    "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.5), tight_layout=True)\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "ax1.plot(nbpt_list, ave_list1, 's-b', label = \"succesive\")\n",
    "ax1.plot(nbpt_list, ave_list2, 's-r', label = \"distant\")\n",
    "ax1.set_xlabel('Number of points')\n",
    "ax1.set_ylabel('Estimated average')\n",
    "ax1.set_xscale('log')\n",
    "ax1.legend()\n",
    "\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "ax2.plot(nbpt_list, std_list1, 'o-b', label = \"succesive\")\n",
    "ax2.plot(nbpt_list, std_list2, 'o-r', label = \"distant\")\n",
    "ax2.set_xlabel('Number of points')\n",
    "ax2.set_ylabel('Estimated standard deviation')\n",
    "ax2.set_xscale('log')\n",
    "ax2.legend()\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusions:\n",
    "- \n",
    "- \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Density probability function and histograms\n",
    "Calculate and plot histograms and/or probability density functions for:\n",
    "- the wake signal\n",
    "- the temperature or wind signal\n",
    "\n",
    "Explore the parameters:\n",
    "- number of bins\n",
    "- density or not"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'np' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-6-e296e33fcc51>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mhist10\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbin_edges10\u001b[0m  \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mhistogram\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbins\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m10\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdensity\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mhistXX\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbin_edgesXX\u001b[0m  \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mhistogram\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbins\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mXX\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdensity\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mfig1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0max1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0max2\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfigsize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m12\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m4.5\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtight_layout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'np' is not defined"
     ]
    }
   ],
   "source": [
    "hist10, bin_edges10  = np.histogram(data, bins = 10, density = True)\n",
    "histXX, bin_edgesXX  = np.histogram(data, bins = XX, density = True)\n",
    "\n",
    "\n",
    "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.5), tight_layout=True)\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "ax1.plot(bin_edges10[:-1], hist10,label = \"nbin = 10\")\n",
    "ax1.set_xlabel('Velocity (m/s)')\n",
    "ax1.set_ylabel('PDF')\n",
    "ax1.legend()\n",
    "\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "ax2.plot(bin_edges10[:-1], hist10,label = \"nbin = 10\")\n",
    "ax2.set_xlabel('Velocity (m/s)')\n",
    "ax2.set_ylabel('PDF')\n",
    "ax2.set_yscale('log')\n",
    "ax2.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discovering data fitting\n",
    "Polynomial fit on data with more and more parameters\n",
    "\n",
    "- We define polynomial functions of increasing order\n",
    "- we perform fitting of data points with higher and higher order poynomials\n",
    "- we estimate the R-square and plot the data points along with the fitting polynomial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\romain\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\minpack.py:808: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  category=OptimizeWarning)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAI4CAYAAABndZP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeZyN5f/H8ddlkiWUUNHwHWU3ZoYGlSK7soUKlbKU9KVQlqm+5Ns6otBCkaUFYaxFJJJ9GSVfS0YioTJkGQxmxvX749b8osHMOOfc55x5Px8PD+Oc+5z7faZzPn3OdV/3fRlrLSIiIiLBJJfbAUREREQ8TQ2OiIiIBB01OCIiIhJ01OCIiIhI0FGDIyIiIkHnCrcD/F3RokVtWFiY2zFEJJvWr19/wFpbzO0cmaWaIxL4LlR3/KrBCQsLIz4+3u0YIpJNxphf3M6QFao5IoHvQnVHh6hEREQk6KjBERERkaCjBkdERESCjl/NwclISkoKe/bs4eTJk25HkRwib968hIaGkjt3brejiAtUc8TXVHO8w+8bnD179lCwYEHCwsIwxrgdR4KctZaDBw+yZ88eSpcu7XYccYFqjviSao73+P0hqpMnT1KkSBEVGvEJYwxFihTRt/ccTDVHfEk1x3v8vsEBVGjEp/R+O9ewhQluR/A5vQfEl/R++ydP1J2AaHBExD0jFm13O4KI5DCeqDtqcLJo0KBBDB069KLbzJo1iy1btng1x759+7jvvvsuud1rr73m1Rwi4l2qOSLZE7QNjpvD6r4oNiVKlCAuLu6S26nYSHYMW5hAWMxcwmLmAqT/nBMPV2WWao5DNUeyy9N1J2gbHE8Oq7/66quUL1+eBg0asG3btvTbx4wZQ/Xq1YmMjKRNmzacOHGClStXMmfOHPr27UtUVBQ7duzIcLvzDRo0iA4dOlCvXj3Kli3LmDFjAGeGfd++fQkPD6dKlSpMmTIFgF27dhEeHg7AhAkTaN26NU2aNKFs2bL069cPgJiYGJKTk4mKiuKhhx7y2O9Dgl/vhuXYFduUXbFNAdJ/7t2wnMvJ/JdqjmqOXB6P1x1rrd/8ueWWW+z5tmzZ8o/bMuNf/b/I1uPOFx8fb8PDw+3x48ftkSNH7M0332yHDBlirbX2wIED6du98MIL9u2337bWWvvoo4/aadOmpd93oe3+7sUXX7QRERH2xIkTNjEx0YaGhtq9e/fauLg426BBA5uammp///13W7JkSbtv3z67c+dOW7lyZWuttePHj7elS5e2hw8ftsnJybZUqVJ29+7d1lprr7rqKo/8HnKa7L7vglFWPktAvPWDWpLZP6o5qjn+QjXnXJ6oO35/HZysGLYw4ZxvUX8Nc/WsXzbbHeCyZcto1aoV+fPnB6BFixbp923atIn//Oc/HD58mGPHjtG4ceMMnyOz27Vs2ZJ8+fKRL18+6taty9q1a1m+fDnt27cnJCSE66+/njp16rBu3ToiIiLOeWz9+vW5+uqrAahUqRK//PILJUuWzNZrFvm7nvXLuh3Bb6nmqOaId3ii7gRVg9O7Ybn0ohIWMzd9mOtyXegUvo4dOzJr1iwiIyOZMGECS5Ysuaztzt+PMQanOb20PHnypP8cEhJCampqph4ncik6LHVhqjkO1RzxNE/UnaCdg+MptWvXZubMmSQnJ5OUlMTnn3+efl9SUhLFixcnJSWFiRMnpt9esGBBkpKSLrnd+WbPns3Jkyc5ePAgS5YsoXr16tSuXZspU6aQlpZGYmIiS5cupUaNGpnOnzt3blJSUrL4qkXELao5Ip4RtA2Op4bVq1WrRtu2bYmKiqJNmzbceeed6fe9/PLL1KxZk4YNG1KhQoX029u1a8eQIUOoWrUqO3bsuOB256tRowZNmzbl1ltvZcCAAZQoUYJWrVoRERFBZGQk9erV44033uCGG27IdP6uXbsSERGhCX8iXqaa41DNEX9hMjsc6QvR0dE2Pj7+nNu2bt1KxYoVXUrkO4MGDaJAgQL06dPH7ShCznnfeZoxZr21NtrtHJmlmqOa4y9yyvvOGy5Ud4J2BEdERERyrqCaZBzIBg0a5HYEEclBVHMk2GkER0RERIKOGhwREREJOmpwREREJOh4tcExxvQ2xmw2xmwyxkw2xuT15v5ERFR3RAS82OAYY24EngairbXhQAjQzlv78yd33XUX5596er7OnTtz3XXXpS9el5HExERq1qxJ1apVWbZsGffccw+HDx/m8OHDjBw5Msu5lixZwtVXX03VqlWpUKFChqeHbtmyheuuu4677777nCuT/vrrr9StW5eKFStSuXJlRowYkeX9Z2T+/PmUL1+eMmXKEBsbe8Htpk6dSqVKlahcuTIPPvhg+u0hISFERUURFRV1ziXtH3roIcqXL094eDidO3dOv/DYoUOH0q/zUaNGDTZt2uSR1yH+IafWnUvVnMx+flVz/t+Fas7u3btp1KgRFStWpFKlSuzatQuAxYsXU61aNcLDw3n00UfTX8uRI0do3rw5kZGRVK5cmfHjx3vkdUgmZLRAlSf+ADcCvwLX4pyt9QXQ6GKP8eTCd26qU6eOXbdu3UW3+fbbb+369evTF6/LyOTJk+0jjzzyj9v/vuhdVnzzzTe2adOm1lprT5w4YcuXL2+XL1+efv/evXtt5cqV7bJly+wzzzxjO3XqlH7fvn377Pr166211h49etSWLVvWbt68OcsZ/i41NdXedNNNdseOHfbUqVM2IiIiw+dMSEiwUVFR9s8//7TWWvvHH3+k33ehhf3mzp1rz5w5Y8+cOWPbtWtnR44caa21tk+fPnbQoEHWWmu3bt1q69Wrl+HjA/F95w9webHNrNadnFJzMvv5Vc1xXKzm1KlTx3711VfWWmuTkpLs8ePHbVpamg0NDbXbtm2z1lo7YMAA++GHH1prrX311Vdtv379rLXW7t+/3xYuXNieOnXqH/sMxPedv7hQ3fHaCI61di8wFNgN/AYcsdZ+df52xpiuxph4Y0x8YmKit+Jk265du6hQoQKPPvooERER3HfffZw4cQKARYsWUbVqVapUqULnzp05derUOY8dO3YsvXv3Tv/3mDFjeOaZZwDncuzXXnvtBfe7YcMG+vXrx7x584iKiiI5OZmwsDAOHDhATEwMO3bsICoqir59+2brdeXLl4+oqCj27t0LwNGjR2nbti2jR4/mjjvu4M0336RYsWIMHDgQgOLFi1OtWjXAuSx8xYoV0x+bXWvXrqVMmTLcdNNNXHnllbRr147Zs2f/Y7sxY8bQvXt3ChcuDMB11113yee+5557MMZgjKFGjRrs2bMHcL4t1q9fH4AKFSqwa9cu/vjjj8t6HeI/MlN3cmLNycznVzXn/12o5mzZsoXU1FQaNmwIQIECBcifPz8HDx4kT548lCvnrJ/UsGFDpk+fDjjreyUlJWGt5dixY1x77bVccYWu0OILXvstG2MKAy2B0sBhYJox5mFr7ad/385aOxoYDc5VRS/6pL16wYYNng0aFQXDh190k23btjF27Fhq1apF586dGTlyJD169KBjx44sWrSIcuXK8cgjjzBq1Ch69eqV/rh27doRERHBG2+8Qe7cuRk/fjwffPBBJmNF8dJLLxEfH8+77757zn2xsbFs2rSJDZfxuzh06BDbt2+ndu3aABQqVIhly5ads83gwYMzfOyuXbv4/vvvqVmz5j/umzhxIkOGDPnH7WXKlCEuLu6c2/bu3XvO6sOhoaGsWbPmH49NSEgAoFatWqSlpTFo0CCaNGkCwMmTJ4mOjuaKK64gJiaGe++995zHpqSk8Mknn6QPb0dGRjJjxgzuuOMO1q5dyy+//MKePXu4/vrrM3ytElgyU3eyUnN69ep1WZ+zjERFRTHcxZpzoc+vas7/u1DNSUhI4JprrqF169bs3LmTBg0aEBsbS9GiRUlJSSE+Pp7o6Gji4uL49ddfAejRowctWrSgRIkSJCUlMWXKFHLl0vk9vuDN33IDYKe1NtFamwLMAG734v68pmTJktSqVQuAhx9+mOXLl7Nt2zZKly6d3rE/+uijLF269JzHXXXVVdSrV48vvviCH3/8kZSUFKpUqeLz/H+3bNkyIiIiuOGGG2jWrFmW1pgBOHbsGG3atGH48OEUKlToH/c/9NBDbNiw4R9/zi80QIarFme0inJqairbt29nyZIlTJ48mccee4zDhw8DzvHw+Ph4Jk2aRK9evdixY8c5j/33v/9N7dq109fziYmJ4dChQ0RFRfHOO+9QtWpVfZsKLkFRd7xVcy71+fWGYKo5qampLFu2jKFDh7Ju3Tp+/vlnJkyYgDGGzz77jN69e1OjRg0KFiyYXlcWLFhAVFQU+/btY8OGDfTo0YOjR49m6Xcg2ePNyr4buNUYkx9IBuoDF595eymX+NbjLed/AIwxGX5QMvLYY4/x2muvUaFCBTp16uSNeOd47733GDNmDADz5s2jRIkS59x/55138sUXX5CQkMAdd9xBq1atiIqKytRzp6Sk0KZNGx566CFat26d4TZZ+TYVGhqa/i0HYM+ePf/I+9d2t956K7lz56Z06dKUL1+e7du3U7169fTtb7rpJu666y6+//57br75ZgD++9//kpiYeM432EKFCqVP8rPWUrp0aUqXLp2p1y8BwaN151IjLd7ijZqTmc9vduSkmhMaGkrVqlW56aabALj33ntZvXo1Xbp04bbbbksfkfrqq6/SR4HGjx9PTEwMxhjKlClD6dKl+fHHH7O0Qrtkjzfn4KwB4oDvgP+d3ddob+3Pm3bv3s2qVasAmDx5MnfccUf6/I2ffvoJgE8++YQ6der847E1a9bk119/ZdKkSbRv394jeQoWLEhSUlKG93Xv3j39G0xGH9y/lCtXjueee+6CQ8Lns9bSpUsXKlasmD6PKCNZ+TZVvXp1tm/fzs6dOzl9+jSfffbZOWdC/eXee+/lm2++AeDAgQMkJCRw0003cejQofQ5CAcOHGDFihVUqlQJgA8//JAFCxYwefLkc4aDDx8+zOnTp9O3qV27ts++yYr3BUvd8XTNyezn90JUc5yaU716dQ4dOsRfc7cWL16cXnP2798PwKlTpxg8eDDdunUDoFSpUixatAiAP/74g23btqU3SOJdXj0QaK190VpbwVobbq3tYK09delH+Z+KFSvy0UcfERERwZ9//smTTz5J3rx5GT9+PPfffz9VqlQhV65c6W/o8z3wwAPUqlUrfcIaQPv27bntttvYtm0boaGhjB07NtN5ihQpQq1atQgPD8/2hD+Abt26sXTpUnbu3HnJbVesWMEnn3zC4sWL00/LnjdvXrb3DXDFFVfw7rvv0rhxYypWrMgDDzxA5cqVARg4cCBz5swBoHHjxhQpUoRKlSpRt25dhgwZQpEiRdi6dSvR0dFERkZSt25dYmJi0otNt27d+OOPP7jtttvS5xaAs2Jv5cqVqVChAl9++aXHTj0V/xEMdcfTNedyP7+qOU7NCQkJYejQodSvX58qVapgreXxxx8HYMiQIVSsWJGIiAiaN29OvXr1ABgwYAArV66kSpUq1K9fn8GDB1O0aNHLeh2SOSazw56+EB0dbc+/loPbS8jv2rWLZs2aXdb1Upo1a0bv3r3Tz94R/+f2+y5QGWPWW2uj3c6RWao54i/cft8FsgvVHU3l9qLDhw9Trlw58uXLp0IjIl6nmiPy/3T6yCWEhYVl+5vUNddckz7RTEQkM1RzRDwjIEZw/OkwmgQ/vd9E7wHxJb3fvMPvG5y8efNy8OBBvQHEJ6y1HDx4kLx5tT5jTqWaI76kmuM9fn+IKjQ0lD179uCPl1SX4JQ3b15CQ0PdjiEuUc0RX1PN8Q6/b3D+utCSiIgvqOaIBAe/P0QlIi5auBBWrnQ7hYhIlqnBEZF/OnwYHnsMGjWC115zO42ISJapwRGRc82dC+HhMH489O8PGVzyXkTE36nBERHHwYPQoQM0awaFC8Pq1RAbCzq7Q0QCkBocEYEZM6ByZfjsMxg4ENavh+rV3U4lIpJtfn8WlYh40f790KMHTJsGVavCggUQGel2KhGRy6YRHJGcyFqYPBkqVYLZs+HVV2HNGjU3IhI0NIIjktP89hs8+aTT2NSoAePGOYenRESCiEZwRHIKa2HCBGfUZsECGDrUucaNmhsRCUIawRHJCXbvhieegPnz4c47YexYKFvW7VQiIl6jERyRYGYtfPCBc12bpUvhnXdgyRI1NyIS9DSCIxKsfv4ZHn8cFi+GevXgww9BayyJSA6hERyRYHPmjDNSU6UKrFsHo0fD11+ruRERv7d9+3aGDx/OTz/9dNnPpREckWCSkACdO8OKFXD33c7hqZIl3U4lIpKh06dP8+233zJ37lzmzp2b3tjkyZOHMmXKXNZzq8ERCQZpafDWW85ViPPmhY8+cpZdMMbtZCIi5zh06BDz5s1j9uzZzJ8/n6SkJPLmzUvdunXp1asX99xzD6U9MOLstQbHGFMemPK3m24CBlprh3trnyI50ubNzqjN2rVw770wciQUL+52Kleo7oj4p7179zJr1ixmzJjBt99+S1paGjfccAPt27enefPm1KtXj/z583t0n15rcKy124AoAGNMCLAXmOmt/YnkOCkpMHgwvPQSFCrkXJm4bdscPWqjuiPiP3bt2kVcXBwzZsxg1apVAFSsWJF+/frRsmVLqlevTq5c3psK7KtDVPWBHdbaX3y0P5HgtmEDdOrk/N22rTOpuFgxt1P5G9UdER/79ddfmTp1KlOnTmXt2rUAVK1alVdeeYXWrVtTsWJFn2XxVYPTDpic0R3GmK5AV4BSpUr5KI5IgDp1ylk36vXXoUgRZxXwVq3cTuWvMqw7qjkinrV//36mTZvGpEmTWLlyJQDVqlVj8ODB3Hfffdx0002u5DLWWu/uwJgrgX1AZWvtHxfbNjo62sbHx3s1j0jAWrvWmWuzebMzgXj4cLj2WrdTncMYs95aG+0HOTJVd1RzRLLn2LFjzJw5k0mTJrFw4ULS0tIIDw+nffv2PPDAA5d9BlRWXKju+GIE527gu0s1NyJyAcnJ8OKL8OabzuThuXPhnnvcTuXvVHdEPCwtLY1vvvmGjz/+mOnTp3PixAlKlSpF3759efDBB6lSpYrbEc/hiwanPRc4PCUil7BiBXTpAtu2OVclHjIErr7a7VSBQHVHxEMSEhIYP348n376KXv27KFQoUI89NBDdOjQgVq1anl1ovDl8GqDY4zJDzQEnvDmfkSCzvHj8MIL8PbbUKoULFwIDRq4nSogqO6IXL5jx44RFxfHuHHjWLZsGbly5aJx48YMHTqUFi1akC9fPrcjXpJXGxxr7QmgiDf3IRJ0vvkGHnvMWUuqRw9nQnGBAm6nChiqOyLZY61l/fr1jB49msmTJ3Ps2DHKlStHbGwsHTp0oESJEm5HzBJdyVjEXyQlQb9+8P77UKYMfPst1K7tdioRCXJHjx5l0qRJjB49mu+//578+fPTtm1bunTpwu23344J0GtrqcER8QcLFkDXrrBnDzzzDLz8Mnj4qp4iIn+3ceNGRo4cyaeffsrx48eJjIxk5MiRPPjgg1wdBHP91OCIuOnQIXj2WRg/HipWdCYV33qr26lEJEidPn2a6dOn895777FixQry5s1L+/bteeKJJ6hRo0bAjtZkRA2OiFvmzIFu3WD/fnj+eRgwwFkoU0TEw37//Xfef/99PvjgA37//XfKlCnDm2++SceOHbnWz66n5SlqcER87cAB6NkTJk2CiAj4/HO45Ra3U4lIEFq3bh0jRoxg6tSppKSkcM899/DUU0/RqFEjvz2921PU4Ij4UlwcdO8Of/4JgwbBc8/BlVe6nUpEgkhaWhqzZ8/mrbfeYsWKFRQsWJAnn3ySHj16ULZsWbfj+YwaHBFf+OMPp7GZPt0ZrVm40Bm9ERHxkGPHjjF27FhGjBjBzp07KV26NMOHD6dTp04UKlTI7Xg+F9zjUyJusxYmToRKlZxDUa+9BqtXX1ZzM2xhggcDikig+/3333nhhRcoWbIkvXr1okSJEkyfPp3t27fTs2dPjzQ3gVh31OCIeMvevdCyJTz8MJQrBxs2OIekrri8gdMRi7Z7KKCIBLKEhASeeOIJwsLCeP3116lfvz6rV69m+fLltG7dmpCQEI/tKxDrjg5RiXiatTBhAvTuDadPw1tvwdNPgweLjYjkXN9//z2vv/46cXFxXHnllXTs2JFnn302R82vyQw1OCKe9MsvzgX7vvrKuQrxhx+CB4rOsIUJ53yDCouZC0DP+mXp3bDcZT+/iPi/pUuX8vrrrzN//nwKFSpETEwMPXv25Prrr/fK/gK97hhrrdsZ0kVHR9v4+Hi3Y4hk3Zkz8MEHzlIL1sIbbzjXuPHCaZhhMXPZFdvU48/rCcaY9dbaaLdzZJZqjvg7ay2LFi3ipZdeYtmyZRQrVozevXvz73//26dXGw7EuqMRHJHLtWOHszjmkiXQsCGMHg1hYW6nEpEAZq1l/vz5vPTSS6xevZobb7yRt99+my5dupBfy7hkiiYZi2RXWhqMGOGcEfXddzBmjLOmlJebm571dZxdJFhZa5k3bx41a9bknnvu4bfffmPUqFHs2LGDp556yrXmJhDrjkZwRLLjxx+hSxdYuRKaNnVWAA8N9cmuA+HYt4hkjbWWBQsWMGjQINasWUNYWBgffvghjzzyCLlz53Y7XkDWHY3giGRFaioMHgxRUbB1K3zyiXN9Gx81NyISfBYvXkytWrW4++67+e233xgzZgwJCQl06dLFL5qbQKURHJHM+t//oHNniI+HVq1g5Ei44Qa3U4lIgFq9ejUvvPACixcvJjQ0lFGjRtG5c2eu1PItHqERHJFLSUmBl15yllj45ReYOtVZckHNjYhkw8aNG2nRogW33XYbmzZtYvjw4Wzfvp1u3bqpufEgjeCIXMx33zmjNj/8AO3bO5OKixVzO5WIBKCdO3cycOBAJk6cyNVXX82rr77K008/TYECBdyOFpTU4Ihk5NQpePlliI11GpqZM+Hee91OJSIBKDExkVdeeYVRo0YREhJCv3796N+/P4ULF3Y7WlBTgyNyvrVroVMn2LIFHn0Uhg0DFSIRyaITJ07w1ltvMXjwYJKTk+ncuTMvvvgiN954o9vRcgTNwRH5S3Iy9O0Lt90GR4/CvHnOmlJqbkQkC9LS0hg3bhxly5ZlwIABNGzYkM2bNzN69Gg1Nz7k1QbHGHONMSbOGPOjMWarMeY2b+5PJNuWL4fISBg61Lkq8aZNcPfdbqeSbFDdETctWLCAqlWr0qVLF0qVKsWyZcuYMWMG5cuXdztajuPtEZwRwHxrbQUgEtjq5f2JZM3x49Czp7MwZkoKfP21s6aUD9d4EY9T3RGf27JlC3fffTdNmjTh+PHjTJ06lZUrV3LHHXe4HS3H8lqDY4wpBNQGxgJYa09baw97a38iWbZ4MVSpAu+8A0895Vznpn59t1PJZVDdEV87cOAA3bt3JyIiglWrVvHmm2+ydetW7r//fowxbsfL0bw5gnMTkAiMN8Z8b4z50BhzlRf3J5I5R47AE084zcwVV8C33zqnf+tUzWCguiM+kZKSwvDhwylTpgwffPAB3bp146effuKZZ57RtWz8hDcbnCuAasAoa21V4DgQc/5Gxpiuxph4Y0x8YmKiF+NIRoYtTHA7gm99+SWEh8OHH0KfPs71be680+1U4jmXrDuqOe4L9Lrz1VdfERkZSe/evalZsyYbN27k3XffpWjRom5Hk7/xZoOzB9hjrV1z9t9xOIXnHNba0dbaaGttdDFdQM3nRiza7nYE3zh0CDp2hHvugUKFnEUyhwyBfPncTiaedcm6o5rjvkCtOzt27KBly5Y0btyY06dPM2fOHObPn0+lSpXcjiYZ8FqDY639HfjVGPPX1PH6wBZv7U/kgmbPhkqV4NNP4YUXnKsT16zpdirxAtUd8YYTJ04wcOBAKleuzOLFi4mNjWXz5s00b95c82z8mLcv9PcUMNEYcyXwM9DJy/uTTBi2MOGcb1BhMXMB6Fm/LL0blnMrlucdOOBMHv7sM+cU8HnzoGpVt1OJ96nu+KFArDvWWmbNmkXv3r355ZdfePDBBxkyZAglSpRwO5pkgrHWup0hXXR0tI2Pj3c7Ro4SFjOXXbFN3Y7hWdbCtGnQowccPgwDBkD//qCJf15njFlvrY12O0dmqea4IxDqzvbt23nqqadYsGAB4eHhvPvuu9SpU8ftWJKBC9UdLdUgweX33+Hf/3bWjoqOdk4FDw93O5WIBIiTJ08SGxtLbGwsV155JcOGDaN79+7kzp3b7WiSRVqqIYfrWb+sz/bl1TMnrIVPPnHm2sybB4MHw6pVam5E/JC/1p2/Rmv++9//0qpVK7Zt20avXr3U3AQoNTg5nC+PfXvtzIk9e6BZM3jkEahYETZsgH79nGvciIjf8be689tvv9G2bVuaNGlCSEgICxcuZPLkyRQvXtwHCcVb1OBI4LLWuZ5N5crwzTcwfDgsXQoVKridTEQCwJkzZ3j//fepWLEis2fP5qWXXmLjxo00aNDA7WjiAfqKK17ltTMndu2Cxx931o666y6n0bn55ssLKyJBITN1Z9OmTXTt2pVVq1ZRr1493n//fcqW9d2hM/E+NTjiVb0blksvKB45c+LMGRg1yjkryhjn565dIZcGI0XEcbG6c/LkSV599VViY2O55ppr+Pjjj3n44Yd1PZsgpAZHAsdPP0GXLs5hqMaNYfRoKFXK7VQiEiCWL1/O448/zo8//sgjjzzCm2++qeUVgpi+9vqZQF+j5WKyfeZEWhq89RZERDhrR40f76wpdYnmJph/lyKeFMyflZ71y3L06FG6d+/OnXfeycmTJ1mwYAEfffSRx5ubYP49BiI1OH4mUNdoyYxszbnZsgVq1YJnn4UGDZx/d+zoHJ66hGD+XYp4UjB/Viqm/Ux4eDijRo2iV69e/O9//6NRo0Ze2Vcw/x4DkRoc8U+pqfD6687SCtu3w8SJzppSukS6iGTCoUOH6NSpE3fffTcFChRg5cqVDBs2jAIFCrgdTXxESzX4gfNn/P/Fn9do8aqNG+jqRDsAACAASURBVKFzZ1i/Hu67D959F66/PlMP1e/SXVqqIXAE82dlzpw5dOvWjf379xMTE8OAAQPIkyePV/YVzL/HQHHBumOt9Zs/t9xyi83p/tX/C7cjuOfUKWtffNHaK66w9rrrrJ027bKeLkf/Ll0CxFs/qCWZ/aOa4wiWz8qff/5pO3ToYAEbGRlpv/vuO5/uP1h+j4HmQnVHZ1GJf1i/Hjp1gv/9Dx56yLlon85uEJFMmjt3Lo8//jiJiYm8+OKLPP/881ypBXZzNM3B8TO+XKPFL5w8Cc8/DzVrwoEDzjybTz/1SHOT436XItkUyJ+VI0eO0LlzZ5o1a0bRokVZs2YNgwYNcqW5CeTfYzDSHBxxz6pVzlybH390Rm/eeguuucbtVHIZNAdHfGnRokV06tSJffv2eX2ujfivC9UdjeCI7504Ac8845z+feIELFgA48apuRGRTElOTqZXr140aNCAfPnysXLlSl555RU1N3IOzcER3/r2W+dqxDt2wJNPQmwsFCrkdioRCRDr1q3jkUce4ccff+Tpp5/m9ddfJ3/+/G7HEj+kERzxjaQk6N7dWRjTWmf175Ej1dyISKakpqby8ssvc9ttt3Hs2DEWLlzIiBEj1NzIBWkER7xv4UJn5e/du6FXL3jlFbjqKrdTiUiA+Pnnn3n44YdZtWoVDz74IO+99x7X6JC2XIJGcMR7jhxxGptGjSBPHli2DIYNU3MjIplirWXChAlERkayZcsWJk6cyMSJE9XcSKaowRHvmDsXKld2Jg/36wcbNjiTikVEMuHQoUM88MADdOrUiVtuuYWNGzfy4IMPuh1LAogaHPGsP/+ERx6BZs2gcGFYvRoGD4Z8+dxOJiIBYunSpURGRjJr1ixiY2NZtGgRpUqVcjuWBBg1OOI5M2dCpUoweTIMHAjx8VC9utupRCRApKSkMGDAAOrWrUuePHlYtWoV/fv3JyQkxO1oEoC8OsnYGLMLSALSgNRAugCYZMH+/fDUUzB1KkRFwfz5zt8iLlDdCUy7du2iffv2rF69mk6dOvH2229r5W+5LL44i6qutfaAD/YjvmYtTJniNDdHjzpnR/XrB7lzu51MRHUngEybNo3HH38cay2fffYZbdu2dTuSBIFLHqIyxvQwxhT2RRgJIL/9Bq1bQ/v2cNNN8N138MILam7EI1R3cobk5GS6devGAw88QIUKFdiwYYOaG/GYzMzBuQFYZ4yZaoxpYowxWXh+C3xljFlvjOma0QbGmK7GmHhjTHxiYmIWnlpcYS189JEz12b+fHjjDVixwjljSsRzvFZ3VHP8w5YtW6hevToffPAB/fv3Z9myZZQuXdrtWBJELtngWGv/A5QFxgIdge3GmNeMMTdn4vlrWWurAXcD3Y0xtTN4/tHW2mhrbXSxYsWyll5869df4Z57oGNHCA+HH36Avn3hCl0vUjzLm3VHNcd9EyZMoHr16iQmJrJgwQJiY2PJrdFf8bBMnUVlnSXHfz/7JxUoDMQZY964xOP2nf17PzATqHFZacUd1sLo0c4ozbJl8M47zppS5cq5nUyCmOpO8Dl+/DgdO3akU6dO1KhRgw0bNtCoUSO3Y0mQyswcnKeNMeuBN4AVQBVr7ZPALUCbizzuKmNMwb9+BhoBmzySWnzn55+hQQN44gnnlO///Q969IBcusKAeI/qTvDZvHkz1atX5+OPP2bgwIF8/fXXFC9e3O1YEsQyc2yhKNDaWvvL32+01p4xxjS7yOOuB2aePXR+BTDJWjs/20nFt86cgffeg5gYCAmBDz5wll3I0lQIkWxT3Qkin3zyCU888QQFCxbkq6++okGDBm5Hkhzgkg2OtXbgRe7bepH7fgYis5lL3JSQAF26wPLl0KSJc3iqZEm3U0kOoroTHE6ePMnTTz/NmDFjqF27Np999plGbcRndJxB/l9aGgwdCpGRsGkTjB8P8+apuRGRLPv555+5/fbbGTNmDP3792fRokVqbsSndPqLOLZsgU6dYO1aaNECRo2CEiXcTiUiAejzzz+nQ4cOGGOYM2cOzZs3dzuS5EAawcnpUlLgtdegalXYscNZR2rWLDU3IpJlaWlpvPDCC7Ro0YKbb76Z7777Ts2NuEYjODnZDz84ozbffw8PPOCc/n3ddW6nEpEAdODAAdq3b8/XX3/NY489xjvvvEPevHndjiU5mEZwcqLTp+HFFyE6Gvbtg+nTnTWl1NyISDasW7eOatWqsWzZMj788EPGjBmj5kZcpwYnp1m3Dm65BV56yVlHavNmZ00pEZFsGDt2LHfccQchISGsWLGCLl26uB1JBFCDk3MkJ0P//nDrrfDnn/D55/Dxx1CkiNvJgsKwhQluRxDxqVOnTtGtWzcee+wx6tSpQ3x8PLfccovbsXIU1Z2LU4OTE6xc6UwifuMNZ87N5s3Q7GLXSpOsGrFou9sRRHxm37593HXXXekLZX755ZcU0Zcln1PduThNMg5mx4/DCy/A228717L56ito2NDtVCISwFauXEmbNm1ISkpi6tSp3H///W5HEsmQGpxgtWSJczXin3+Gf/8bYmOhYEG3UwWVYQsTzvkGFRYzF4Ce9cvSu6EWIpXgM2bMGLp3706pUqVYuHAh4eHhbkfKcVR3Mk8NTrBJSnLm2owaBTff7DQ6deq4nSoo9W5YLr2ghMXMZVdsU5cTiXhHSkoKvXr1YuTIkTRq1IjJkydz7bXXuh0rR1LdyTzNwQkmCxZAeDi8/z707g0bN6q5EZHLkpiYSIMGDRg5ciR9+vRh7ty5am4kIGgEJxgcPgzPPgvjxkGFCrBiBdx2m9upcpSe9cu6HUHE43744QdatmzJH3/8waeffspDDz3kdiT5G9Wdi9MITqD7/HOoXBk++giee865KrGaG5/TsW8JNjNnzqRWrVqkpqaybNkyNTd+SHXn4tTgBKqDB+Ghh5yFMYsUgTVrnDWldPVQEbkM1lpeeeUVWrduTeXKlVm3bh3R0dFuxxLJMjU4gSguDipVgqlTYdAgiI93rk4sInIZTpw4Qfv27RkwYAAPP/ww3377LcWLF3c7lki2aA5OINm/H7p3dxqcatWc69pERrqdSkSCwG+//UaLFi1Yv349gwcPpm/fvhhj3I4lkm1qcAKBtTB5Mjz9tHMa+OuvQ58+cIX+84nI5fv+++9p3rw5hw8fZtasWbRo0cLtSCKXTYeo/N3evdCypTPfpmxZ2LABYmLU3IiIR8ycOZM77riDXLlysWLFCjU3EjTU4Pgra2H8eOcMqa+/hrfeguXLoWJFt5OJSBCw1jJ48GBat25NeHg4a9euJVKHvCWIaBjAH+3eDY8/7syxqV0bxo6FMmXcTiUiQSIlJYUnn3ySsWPH0rZtW8aPH0++fPncjiXiURrB8SdnzjhXIa5c2blY33vvwTffqLkREY85fPgwTZo0YezYsfznP/9h0qRJam4kKHl9BMcYEwLEA3uttc28vb+AtWMHPPaYs3ZUgwYwZgyEhbmdSiTgqOZc2M6dO2natCk//fQTEyZM4NFHH3U7kojX+OIQVU9gK1DIB/sKPGlp8O678PzzzsThMWOcVcB1eqZIdqnmZGDNmjU0b96c1NRUvvrqK+666y63I4l4lVcPURljQoGmwIfe3E/A2rbNmWPTqxfcdRds3uyM4qi5EckW1ZyMzZw5k7p161KwYEFWrVql5kZyBG/PwRkO9APOXGgDY0xXY0y8MSY+MTHRy3H8RGoqvPGGc5G+rVvh44/hiy8gNNTtZCKBTjXnPMOHD6dNmzZERkayevVqypcv73YkEZ/wWoNjjGkG7LfWrr/Ydtba0dbaaGttdLFixbwVx39s2uQshtm/P9xzD2zZAh06aNRG5DKp5pwrLS2Nnj170rt3b1q1asXixYsJ5tcrcj5vjuDUAloYY3YBnwH1jDGfenF//i0lBV5+2VliYdcumDIFpk+HG25wO5lIsFDNOSs5OZn777+ft99+m969ezN16lSdKSU5jtcaHGvtc9baUGttGNAOWGytfdhb+/Nr338P1avDwIHQpo0zavPAAxq1EfEg1RzHwYMHqV+/PrNmzWLEiBG89dZbhISEuB1LxOd0oT9vOnXKGbWJjYVixWDmTLj3XrdTiUiQ2rlzJ02aNOGXX35h2rRptGnTxu1IIq7xSYNjrV0CLPHFvvzGmjXQubMzWvPoozBsGBQu7HYqkRwhJ9ac9evX07RpU06fPs2iRYuoVauW25FEXKUrGXtacjL07Qu33w5Hj8K8eTBhgpobEfGaBQsWUKdOHfLmzcvKlSvV3IigBsezli1zTv0eOtS5ns2mTXD33W6nEpEg9umnn9KsWTPKlCnDqlWrqFChgtuRRPyCGhxPOHYMnn4a6tRxrnGzaBF88AFcfbXbyUQkiA0dOpQOHTpw55138u2331K8eHG3I4n4DTU4l2vxYoiIgHfegR49YONGqFfP7VQiEsTOnDnDs88+S9++fXnggQf48ssvuVpfqETOoQYnu44ehSeegPr1nTWkli6Ft9+GAgXcTiYiQSwlJYVHHnmEt956i6eeeorJkyeTJ08et2OJ+B01ONkxfz5Urgwffgh9+sCGDXDnnW6nEpEgd/z4cVq2bMnEiRN57bXXGDFiBLlyqYyLZETXwcmKQ4fgmWecs6IqVoSVK6FmTbdTiUgO8Oeff9KsWTPWrFnDmDFjeOyxx9yOJOLX1PpnwrCFCTBnjjNq88kn8PzzztWJ1dyIiJcMW5iQ/vPevXupXbs23333HXFxcWpuRDJBDc6lHDjATT27QsuWztWI166FV18FHfMWES8asWg7AAkJCdSqVYvdu3czf/58WrVq5XIykcCgBudCrIVp06BSJe7etgJeegnWrXMWyxQR8YENGzZw5513cuLECZYsWcJdd93ldiSRgKEGJyO//8722k3ggQf4IeQamnccTtjxaoQNXHjOsLGIiCcNW5hAWMxcwmLmcnLPZqrdegcHky1dBn9CNX25EskSTTL+O2th4kTo2ZOyx4/D4MFEPvMM2/6zgF2xTd1OJyJBrnfDcvRuWI758+dzT/OBlLs5jIULF1KyZEm3o4kEHI3g/GXPHmjeHDp0gPLlnVO/+/VzrnEjIuIj06ZNo0WLFuQuEsrSpUvV3Ihkkxoca2HsWOcMqcWLnVW/ly2Dv63n0rN+WRcDikhOMW7cONq1a8ett97Ki+9P4brrrnM7kkjAytnDE7t2weOPw9dfO+tIjR0LN9/8j816Nyzn+2wikqOMGDGCXr160bhxY2bMmEH+/PndjiQS0HLmCM6ZMzByJFSpAqtXOz8vXpxhcyMi4k3WWl555RV69epF69atmT17tpobEQ/IeSM4P/0Ejz0G334LjRrB6NHwr3+5nUpEciBrLf3792fIkCF06NCBcePGcYXm/Yl4RM4ZwUlLc+bXREQ4E4jHjXPWlFJzIyIuOHPmDD169GDIkCE8+eSTTJgwQc2NiAfljE/Tjz9C586wahU0awbvvw833uh2KhHJodLS0ujatSvjxo2jT58+vPHGGxhj3I4lElSCewQnNRViYyEqCrZtc9aRmjNHzY2IuCYlJSX9cNTAgQPV3Ih4SfCO4Gzc6IzarF8PrVvDe+/BDTe4nUpEcrDTp0/Trl07Zs6cyWuvvcZzzz3ndiSRoBV8IzinT8N//wvR0bB7N0ydCtOnq7kREVedPHmS1q1bM3PmTIYPH67mRsTLvDaCY4zJCywF8pzdT5y19kVv7Q+A776DTp2c0ZsHH4QRI6BoUa/uUkT8hyt1JxOSk5O59957+eqrr3j//fd54okn3I4kEvS8OYJzCqhnrY0EooAmxphbvbKnkyfh+eehRg1ITITZs501pdTciOQ0vqs7mXT8+HGaNWvGwoULGTt2rJobER/x2giOtdYCx87+M/fZP9bjO1q92plrs3WrM3rz5ptQuLDHdyMi/s9ndSeTkpKSaNq0KStWrODjjz/m4YcfdiuKSI7j1Tk4xpgQY8wGYD+w0Fq7JoNtuhpj4o0x8YmJiZl/8hMnoE8fqFULjh1zrmkzbpyaG5Ec7lJ1J9s1J4uOHj1KkyZNWLlyJRMnTlRzI+JjXm1wrLVp1tooIBSoYYwJz2Cb0dbaaGttdLFixTL3xEuXQmSkM1rTtSts2gSNG3s2vIgEpEvVnWzVnCw6cuQIjRs3Zu3atUyZMoV27dp5ZT8icmE+OYvKWnsYWAI0uewnW7TIWRjzzBln/ahRo6BQoct+WhEJLh6tO1lw+PBhGjVqRHx8PFOnTqVNmza+3L2InOW1BscYU8wYc83Zn/MBDYAfL/uJ77rLGbnZuBHq1r3spxOR4OG1upNJfzU333//PXFxcbRq1cpXuxaR83jzQn/FgY+MMSE4jdRUa+0Xl/2sISHwzDOX/TQiEpS8U3cy4dChQzRq1IgffviBuLg4WrRo4YvdisgFePMsqo1AVW89v4jI+dyqO3+N3Pzwww9Mnz6d5s2b+zqCiJwneJdqEBHxgSNHjqQ3NzNmzKBZs2ZuRxIRgnGpBhERH/nrbKkNGzYwffp0NTcifkQjOCIi2fDXdW7Wr19PXFycDkuJ+BmN4IiIZFFSUhJ33313+qngLVu2dDuSiJxHIzgiIllw/PhxmjZtypo1a5gyZYpOBRfxU2pwREQy6cSJEzRv3pwVK1YwadIkXcRPxI+pwRERyYSTJ09y7733smTJEj7++GPatm3rdiQRuQg1OCIil3Dq1CnatGnDwoULGTdunBbOFAkAmmQsInIJAwYMYN68eXzwwQd06tTJ7TgikgkawRERuYSYmBiqVaumVcFFAohGcERELuHaa69VcyMSYNTgiIiISNBRgyMiIiJBRw2OiIiIBB01OCIiIhJ01OCIiIhI0FGDIyIiIkFHDY6IiIgEHTU4IiIiEnTU4IiIiEjQUYMjIiIiQUcNjoiIiAQdrzU4xpiSxphvjDFbjTGbjTE9vbUvERFQ3RGR/+fN1cRTgWettd8ZYwoC640xC621W7y4TxHJ2VR3RATw4giOtfY3a+13Z39OArYCN3prfyIiqjsi8hefzMExxoQBVYE1vtifiIjqjkjO5vUGxxhTAJgO9LLWHs3g/q7GmHhjTHxiYqK344hIDnCxuqOaI5IzeLXBMcbkxikyE621MzLaxlo72lobba2NLlasmDfjiEgOcKm6o5ojkjN48ywqA4wFtlpr3/LWfkRE/qK6IyJ/8eYITi2gA1DPGLPh7J97vLg/ERHVHREBvHiauLV2OWC89fwiIudT3RGRv+hKxiIiIhJ01OCIiIhI0FGDIyIiIkFHDY6IiIgEHTU4IiIiEnTU4IiIiEjQUYMjIiIiQSdgG5xhCxPcjiAiOYhqjkhgCdgGZ8Si7W5HEJEcRDVHJLAEbIMjIiIiciFeW6rBG4YtTDjnW1RYzFwAetYvS++G5dyKJSJBSjVHJHAZa63bGdJFR0fb+Pj4TG0bFjOXXbFNvZxIRLLCGLPeWhvtdo7MUs0RCXwXqjs6RCUiIiJBJ2AbnJ71y7odQURyENUckcASsA2Ojn+LiC+p5ogEloBtcEREREQuRA2OiIiIBB01OCIiIhJ01OCIiIhI0FGDIyIiIkFHDY6IiIgEHb+6krExJhH4JZObFwUOeDGOW4L1dYFeW6DKymv7l7W2mDfDeJJqTjq9tsCk1+bIsO74VYOTFcaY+EC6JHxmBevrAr22QBXMry0rgvn3oNcWmPTaLk6HqERERCToqMERERGRoBPIDc5otwN4SbC+LtBrC1TB/NqyIph/D3ptgUmv7SICdg6OiIiIyIUE8giOiIiISIbU4IiIiEjQCbgGxxjTxBizzRjzkzEmxu08nmKMKWmM+cYYs9UYs9kY09PtTJ5mjAkxxnxvjPnC7SyeZIy5xhgTZ4z58ex/v9vczuQJxpjeZ9+Lm4wxk40xed3O5BbVncCkmhN4PFl3AqrBMcaEAO8BdwOVgPbGmErupvKYVOBZa21F4FagexC9tr/0BLa6HcILRgDzrbUVgEiC4DUaY24EngairbXhQAjQzt1U7lDdCWiqOQHE03UnoBocoAbwk7X2Z2vtaeAzoKXLmTzCWvubtfa7sz8n4bxhb3Q3lecYY0KBpsCHbmfxJGNMIaA2MBbAWnvaWnvY3VQecwWQzxhzBZAf2OdyHreo7gQg1ZyA5bG6E2gNzo3Ar3/79x6C5MP4d8aYMKAqsMbdJB41HOgHnHE7iIfdBCQC488OhX9ojLnK7VCXy1q7FxgK7AZ+A45Ya79yN5VrVHcCk2pOgPF03Qm0BsdkcFtQnedujCkATAd6WWuPup3HE4wxzYD91tr1bmfxgiuAasAoa21V4DgQ8HM0jDGFcUYpSgMlgKuMMQ+7m8o1qjsBRjUnMHm67gRag7MHKPm3f4cSRMPmxpjcOEVmorV2htt5PKgW0MIYswtneL+eMeZTdyN5zB5gj7X2r2+9cTjFJ9A1AHZaaxOttSnADOB2lzO5RXUn8KjmBCaP1p1Aa3DWAWWNMaWNMVfiTD6a43ImjzDGGJxjqluttW+5nceTrLXPWWtDrbVhOP/NFltrg2I0wFr7O/CrMab82ZvqA1tcjOQpu4FbjTH5z7436xMkExmzQXUnwKjmBCyP1p0rPBbLB6y1qcaYHsACnNnV46y1m12O5Sm1gA7A/4wxG87e9ry1dp6LmSRzngImnv2f389AJ5fzXDZr7RpjTBzwHc6ZNt8T3JeFvyDVHfFDQVdzwPN1R0s1iIiISNAJtENUIiIiIpekBkdERESCjhocERERCTpqcERERCToqMERERGRoKMGR0RERIKOGhwREREJOmpwxOuMMdWNMRuNMXmNMVcZYzYbY8LdziUiwUk1R0AX+hMfMca8AuQF8uGso/K6y5FEJIip5ogaHPGJs5cUXwecBG631qa5HElEgphqjugQlfjKtUABoCDOtyoREW9SzcnhNIIjPmGMmQN8BpQGiltre7gcSUSCmGqOBNRq4hKYjDGPAKnW2knGmBBgpTGmnrV2sdvZRCT4qOYIaARHREREgpDm4IiIiEjQUYMjIiIiQUcNjoiIiAQdNTgiIiISdNTgiIiISNBRgyMiIiJBRw2OiIiIBB01OCIiIhJ01OCIiIhI0FGDIyIiIkFHDY6IiIgEHb9abLNo0aI2LCzM7Rgikk3r168/YK0t5naOzFLNEQl8F6o7ftXghIWFER8f73YMEckmY8wvbmfICtUckcB3obqjQ1QiIiISdNTgiIiISNBRgyMiIiJBx6/m4GQkJSWFPXv2cPLkSbejSA6QN29eQkNDyZ07t9tRxCWqOeJrqjve4fcNzp49eyhYsCBhYWEYY9yOI0HMWsvBgwfZs2cPpUuXdjuOuEQ1R3xJdcd7/P4Q1cmTJylSpIgKjXidMYYiRYrom3sOp5ojvqS64z1+3+AAKjTiM3qv/dOwhQluR/gHY8w4Y8x+Y8ymv912rTFmoTFm+9m/C1/G83smqEgm6P3mHQHR4IiIe0Ys2u52hIxMAJqcd1sMsMhaWxZYdPbfIpJDqcHJokGDBjF06NCLbjNr1iy2bNni1Rz79u3jvvvuu+R2r732mldzSHBb+stSUs1+t2P8g7V2KfDneTe3BD46+/NHwL0+DeUlqjkSLM6cOcXWrR1ITJzuk/0FbYPj5rC6L4pNiRIliIuLu+R2KjaSHcMWJvCvmC+oN+4+/sw9mrCYuYTFzPXLw1V/c7219jeAs39fl9FGxpiuxph4Y0x8YmKix3aumuNQzZELOXRoEX/88SmbN9/Pr78O9/r+grbB8eSw+quvvkr58uVp0KAB27ZtS799zJgxVK9encjISNq0acOJEydYuXIlc+bMoW/fvkRFRbFjx44MtzvfoEGD6NChA/Xq1aNs2bKMGTMGcGbY9+3bl/DwcKpUqcKUKVMA2LVrF+Hh4QBMmDCB1q1b06RJE8qWLUu/fv0AiImJITk5maioKB566CGP/T4k+PVuWI6ZPUuQliuR/Gk12RXblF2xTendsJzb0S6btXa0tTbaWhtdrJjnls1SzVHNkYs7cGAmISEFKVr0Xnbs6M1PP/XB2jPe26G11m/+3HLLLfZ8W7Zs+cdtmfGv/l9k63Hni4+Pt+Hh4fb48eP2yJEj9uabb7ZDhgyx1lp74MCB9O1eeOEF+/bbb1trrX300UfttGnT0u+70HZ/9+KLL9qIiAh74sQJm5iYaENDQ+3evXttXFycbdCggU1NTbW///67LVmypN23b5/duXOnrVy5srXW2vHjx9vSpUvbw4cP2+TkZFuqVCm7e/dua621V111lUd+DzlJdt9zwWbg4oE2139z2dD+EzP9GCDe+qheAGHApr/9extQ/OzPxYFtl3oO1RzVHH8R7HXnzJlUu3x5Ubt5c3t75kyqTUjoYb/5Brt5czublnbysp77QnXH76+DkxXDFiac8y0qLGYuAD3rl832N89ly5bRqlUr8ufPD0CLFi3S79u0aRP/+c9/OHz4MMeOHaNx48YZPkdmt2vZsiX58uUjX7581K1bl7Vr17J8+XLat29PSEgI119/PXXq1GHdunVERESc89j69etz9dVXA1CpUiV++eUXSpYsma3XLAIwe9tsbi95O61Do92OkllzgEeB2LN/z/b2DlVzVHMkc44cWU5KygGKFm2NMSGUKfM2efKU5Oef+3PmzGnCwz0/LyeoGpzeDculF5WwmLnsim3qkee90Cl8HTt2ZNasWURGRjJhwgSWLFlyWdudvx9jzF/fTC8pT5486T+HhISQmpqaqceJZGTX4V388McPDGk4hN63+99hKWPMZOAuoKgxZg/wIk5jM9UY0wXYDdzv7RyqOQ7VHLmUxMQZ5MqVl2uvdU5+NMZQqlQ/Tp7czW+/jcbaNIwJ8eg+g3YOjqfUrl2bmTNnkpycTFJSEp9//nn6fUlJSRQvXpyUCUcq2gAAIABJREFUlBQmTpyYfnvBggVJSkq65Hbnmz17NidPnuTgwYMsWbKE6tWrU7t2baZMmUJaWhqJiYksXbqUGjVqZDp/7ty5SUlJyeKrlpxuzrY5ALQo3+ISW7rDWtveWlvcWpvbWhtqrR1rrT1ora1vrS179u/zz7IKCKo5EmystRw4MIPChRtxxRUFzrmvQIEIrE3h1Kl9Ht9v0DY4PeuX9cjzVKtWjbZt2xIVFUWbNm2488470+97+eWXqVmzJg0bNqRChQrpt7dr144hQ4ZQtWpVduzYccHtzlejRg2aNm3KrbfeyoABAyhRogStWrUiIiKCyMhI6tWrxxtvvMENN9yQ6fxdu3YlIiJCE/4kS2Zvm02FohUoV8T/Rm/8lWqOQzVHzpeUFM+pU3soVqz1P+7Lm9dZnuLkyZ0e36/J7HCkL0RHR9v4+Phzbtu6dSsVK1Z0KZHvDBo0iAIFCtCnTx+3o+R4OeU9dyGHkg9RbEgx+tzeh9gGsVl6rDFmvbU2YCbtqOao5viLYH7f/fzz8+ze/Qa1au0nd+5rz7nvxImfWLu2LOXLj6d48Y7Zev4L1Z2gHcERkeyZt30eaTaNluVbuh1FRAKctZbExOkULlz3H80NQN68pQDjlRGcoJpkHMgGDRrkdgQRAOYkzOH6q66nZmhNt6OIF6nmiC+cOLGV5OQEQkN7ZXh/rlxXkidPSa80OBrBEZF0p1JP8eX2L2lerjm5jMqDiFyexMQZgKFo0QuvnJI3b2mSk3/2+L5VwUQk3ZJdS0g6nUTLCjo8JSKX78CBGRQqdBt58hS/4Db58pUOvBEcY0xvY8xmY8wmY8xkY0xeb+5PRC7P7G2zyZ87//+xd+dhVZb5H8ff9zkshx0UVBQRFBBFcdfMFnMvU6dya1ptsZpKs+ZXVlPN0jpTmY1WU5k1LZqZmaPlrpWTueS+gSKogAsCsu/n/v2BMpobyzk8Z/m+rotrEh6e53MQ7/mee2Vg9ECjowghnFxJSSqFhVsJDb3pktdZLNGUl2dSVVVq0+fbrcBRSrUCJgE9tdadADMw3l7PE0I0jNaaRUmLGNpuKD6ePkbHEUI4uZMnFwIQFna5AqctAGVlh2z6fHsPUXkAPkopD8AXsP1OPg6of//+/Hbp6dlKS0vp3bs3Xbp0ISEhgRdeeOGC1+3bt4+uXbvW7G1x5ZVXAtWH3n3xxRd1zvXxxx8TFhZG165diY+PZ9q0aedds3btWvz8/LjvvvvO+fy2bdvo27cvCQkJJCYm1hzA11CffPIJsbGxxMbG8sknn1zwmu3bt9O3b186d+7MiBEjyM/PB6p/Dj4+PnTt2pWuXbvy4IMP1nzPl19+SWJiIgkJCTUHAQK89957dO7cma5du3LVVVfZ/QRmZ/Lr0V/JKMiQ1VNO6HJtDsC0adNISEigU6dO3HrrrZSWnv9u2R3anGHDhhEcHMyNN9540WvKysoYN24cMTEx9OnTh7S0tJqvvfLKK8TExNC+fXuWLVtW8/mlS5fSvn17YmJiePXV/22vkJqaSp8+fYiNjWXcuHGUl5fb5HU4g7y8/2KxtMPHp90lr/Pxqd4Lx+bzcC50QJWtPoDJQCGQBXx+kWsmApuBzZGRkecdouWMB5Bde+21etOmTRf9utVq1QUFBVprrcvLy3Xv3r31+vXrz7vulVde0c8///x5n1+zZo0ePnx4nXPNnj1bP/zww1rr6sP4mjZtWnNAntZa79y5U3fo0EHv3LlTjx07Vv/5z3+u+VpSUpJOTk7WWmudkZGhW7RooXNzc+uc4WzZ2dk6OjpaZ2dn65ycHB0dHa1zcnLOu65nz5567dq1WmutZ82apf/0pz9prfU5h/+d7eTJk7p169b6xIkTWmut77zzTr1y5UqttdZ5eXk113377bd66NCh532/M/7O2cKTy5/UHn/10CeLTl7+4ougEQ/btMWHLQ/bNNLl2pz09HQdFRWli4uLtdZajxkzRs+ePfu861y9zdFa65UrV+pFixZd8vXMnDlTP/DAA1prrefMmaPHjh2rtdZ69+7dOjExUZeWluqDBw/qtm3b6srKSl1ZWanbtm2rU1JSdFlZmU5MTNS7d+/WWlf/rOfMmaO11vqBBx7Q77zzzgWf6Yy/d5ezfn203rVrzGWvKy3N0GvWoNPTZ9brORdrd+w5RBUCjAKigZaAn1Lq9gsUWO9rrXtqrXuGhYXZK069paWlER8fz1133UViYiKjR4+muLgYgFWrVtGtWzc6d+7MPffcQ1lZ2TnfO2vWLKZMmVLz5w8++IDHH38cpRT+/tXbVVdUVFBRUXHemTDfffcdb731Fh9++CHXXXcdQM33TJ06lZ9++omuXbte8B1RbTRt2pSYmBiOHj0KQEZGBvfeey8LFy6kU6dOfPHFFyQlJfHRRx8BEBcXR2xs9U6tLVu2pFmzZmRlZdXr2WcsW7aMwYMH06RJE0JCQhg8eDBLly4977qkpCSuueYaAAYPHszXX1/6ULaDBw8SFxfHmd+nQYMG1XxPYGBgzXVFRUUXPfPH3Wit+XL3lwxqO4imvk2NjuPW7NHmAFRWVlJSUkJlZSXFxcW0bNnynO91hzYHqg8JDQgIuOQ13377LXfddRcAo0ePZtWqVWit+fbbbxk/fjze3t5ER0cTExPDxo0b2bhxIzExMbRt2xYvLy/Gjx/Pt99+i9aa1atXM3r0aADuuusuFi5c2ODX4AwqKnIpLU3F37/7Za/18mqBUt42n2hsz31wBgGpWussAKXUAuBK4LP63vCxpY+x7dg2G8Wr1rVFV94a9tYlr0lKSmLWrFn069ePe+65h3feeYdHHnmEu+++m1WrVhEXF8edd97Ju+++y2OP/W+t//jx40lMTOTvf/87np6ezJ49m3/9618AVFVV0aNHDw4cOMDDDz9Mnz7n7jlyww038OCDD15wp9FXX32V119/ncWLF9f7dR8+fJjS0tKaE4JbtWrFhg0bar5uNpsv2iW9ceNGysvLadfu/G7Hf/zjHxc8++aaa67h7bffPudzGRkZ55w+HBERQUZGxnnf26lTJxYtWsSoUaP46quvOHLkSM3XUlNT6datG4GBgbz44otcffXVxMTEsG/fPtLS0oiIiGDhwoXndAvPnDmTN998k/LyclavXn2xH5Fb2ZCxgUN5h/hL/78YHcWh7N//GIWFtm1z/P27EhvbuG1Oq1at+OMf/0hkZCQ+Pj4MGTKEIUOGnPNMd2hzauvstsnDw4OgoCCys7PJyMjgiiuuqLnu7Dbrt23Zhg0byM7OJjg4GA8Pj/Oud3WFhVsBCAi4fIGjlAkfn2hKSmxb4NhzDs5h4AqllK+qfps8ENhrx+fZTevWrenXrx8At99+O+vWrSMpKYno6Gji4qrP6rnrrrv48ccfz/k+Pz8/BgwYwOLFi9m3bx8VFRV07twZqP7HvG3bNtLT09m4cSO7du1qlNfy5ZdfkpCQQNu2bZk8eTIWS90Wth09epQ77riD2bNnYzKd/+vzf//3f2zbtu28jws1NPoCx4RcqEflo48+YubMmfTo0YOCggK8vLwACA8P5/Dhw2zdupU333yT3//+9+Tn5xMSEsK7777LuHHjuPrqq4mKiqppYAAefvhhUlJSeO2113jxxRfr9Ppd1Ze7vsTL7MXv4i++V4VoPLZuc3Jzc/n2229JTU0lMzOToqIiPvus3u8168SR2pzauljbZKvPu4OCgi0A+Pt3q9X1Fks0paW2nYNjtx4crfUGpdR8YAtQCWwF3m/IPS/X02Ivv/2FvNgv7oXcd999vPzyy8THxzNhwoTzvh4cHEz//v1ZunQpnTp1anDWZ599liVLlgDVk/R+a9y4ccyYMYP169czfPhwrr/++lofpJefn8/w4cN58cUXz3kXc7a6vJuKiIhg7dq1NX9OT0+nf//+531vfHw8y5cvByA5Obnm9Xl7e+Pt7Q1Ajx49aNeuHcnJyfTs2ZMRI0YwYsQIAN5//33MZvN59x0/fjwPPfTQ5V+4i7NqK/P2zGNYzDCCLEFGx3Eol+tpsRdbtzkrV64kOjq6Ztj25ptv5ueff+b228+bNVBnztTm1FZERARHjhwhIiKCyspK8vLyaNKkSc3nz0hPT68Z6rvQ50NDQzl16hSVlZV4eHicc72rKyzcird3a7y8ajf1xGKJJj9/vU0z2HUVldb6Ba11vNa6k9b6Dq112eW/y/EcPnyY9eurf/Bz5szhqquuIj4+nrS0NA4cOADAp59+yrXXXnve9/bp04cjR47wxRdfcOuttwKQlZXFqVOnACgpKWHlypWXPPH3twICAigoKLjg11566aWadzCX0rdvX+644w6mT59eq2eWl5dz0003ceeddzJmzJiLXleXd1NDhw5l+fLl5Obmkpuby/Llyxk6dOh51504cQIAq9XKiy++WLNaKisri6qqKqB63s3+/ftp27btOd+Tm5vLO++8U7NCY//+/TX3XbJkSc0Yvztbd3gdmQWZjE+QXRwcha3bnMjISH755ReKi4vRWrNq1ao6HezoKm1ObY0cObJmVef8+fMZMGAASilGjhzJ3LlzKSsrIzU1lf3799O7d2969erF/v37SU1Npby8nLlz5zJy5EiUUlx33XXMnz8fqF41OmqUe6xSLCzcUqv5N2dYLNFUVp6iouKUzTLITsa10KFDBz755BMSExPJycnhoYcewmKxMHv2bMaMGUPnzp0xmUznLFM+29ixY+nXrx8hISFAdZfrddddR2JiIr169WLw4MGXXLL4W4mJiXh4eNClS5d6T/gDeOqpp5g9e/ZFG66zzZs3jx9//JGPP/64Zln25Rq0y2nSpAnPPfccvXr1olevXjz//PM0aVJ9GNt9991Xs+x1zpw5xMXFER8fT8uWLWvelf74448kJibSpUsXRo8ezXvvvVfz/ZMnT6Zjx47069ePqVOn1nTrz5gxg4SEBLp27cqbb7550aXp7uTLXV/i4+HDiPYjjI4iTrN1m9OnTx9Gjx5N9+7d6dy5M1arlYkTJ9Y6j6u0OQBXX301Y8aMYdWqVURERNQs9X7++edZtGgRAPfeey/Z2dnExMTw5ptv1iz7TkhIYOzYsXTs2JFhw4Yxc+ZMzGYzHh4ezJgxg6FDh9KhQwfGjh1LQkICAK+99hpvvvkmMTExZGdnc++99zb4NTi6yspCiouTajX/5gwfn+o3pzadaHyhpVVGfTjiks2LLUWui+HDh9csUxaOz+jfucZUUVWhm/2jmR49b7RN7ocsE28waXPck9G/d7Z06tQ6vWYNOivrP7X+nvz8LXrNGvSJE/Pr/LyLtTvSg2NHp06dIi4uDh8fHwYOlK3vheNZm7aWE0UnGJcwzugowgakzRGO4MwE44CA2k0whuohKsCmK6nsuUzcJURFRdV7hVNwcDDJyck2TiSE7Xy560v8vfy5IfYGo6OI06TNEc6usHALnp7N8PKq/YRqT89gPDyCbTpE5RQ9OLqWqweEaCh3+l2rqKpgwb4FjGw/El9PX6PjOBR3+j0QxnO137eCgq0EBHSv85L46qXiblTgWCwWsrOzXe4XQDgerTXZ2dl13qfDWa08uJKckhwZnvoNaXNEY3K1dqeqqpTi4t11WkF1hsXS1qbnUTn8EFVERATp6ek22aJbiMuxWCxEREQYHaNRzN09lyDvIIa2O39pvjuTNkc0Nldqd4qKdqF1ZZ1WUJ3h4xNNdvZitLaiVMP7Xxy+wPH09CQ6OtroGEK4lKLyIhbsXcDYjmPx9vA2Oo5DkTZHiPorLDyzg3F9enCi0bqM8vJjeHs3fENEhx+iEkLY3vw98yksL+TurncbHUUI4UIKCrbg4RGMxRJV5+89s5LKVvNwpMARwg3N3jabmCYxXBV5ldFRhBAu5MwOxvU5c+vMZn+2mocjBY4QbiYlJ4UfDv3AhK4T3ObgPyGE/VmtFRQW7qj1AZu/5e3dBpAeHCFEPX287WNMysSdXe40OooQwoUUF+9D67J6TTAGMJsteHm1lAJHCFF3VdYqPtn+CYPbDiYi0DVWbQghHENDJhifYbFEyxCVEKLuVqeu5kj+ESZ0nWB0FCGEiyko2ILJ5Ievb2y97+Hj01Z6cIQQdTd722yCLcGMih9ldBQhhIupnmDcFaXM9b6HxRJNWVk6b62o33ElZ5MCRwg3car0FN/s+4bfd/o9Fg/X2DVVCOEYtLZSWLit3vNvzqheKq75Yv36BmeSAkcINzF311xKK0u5p9s9RkcRQriY0tJUqqoK8ffv0qD7+PhU74UT5nOswZmkwBHCTczeNpvOzTrTPbxh77CEEOK3ioqqh5T8/DrX+x7TViQz4K3q+TdhvseJmrqEqKlLmLYiuV73c/ijGoQQDbf7xG42ZmzkzSFvyt43QgibKyzcCYCvb8d632PK4DgeGxRDcXEv7l+eRNqrwxuUSXpwhHAD72x6B2+zN7cn3m50FCGECyoq2onF0hYPD/8G3UcpE35+CVRYG35GnhQ4Qri4vNI8Ptn+CeM7jSfML8zoOEIIF1RUtAs/v042u9/kgfVfan6G3QocpVR7pdS2sz7ylVKP2et5QogL+3jbxxRVFPFo70eNjiKEcEFWaxnFxUkNmn/zW1MGxzX4Hnabg6O1TgK6AqjqRfEZwDf2ep4Q4nxWbWXGphn0jehLj5Y9jI4jhHBBxcVJQJVNe3BsobGGqAYCKVrrQ430PCEEsDxlOQdyDkjvjRDCboqKqicY+/vbrgfHFhqrwBkPzLnQF5RSE5VSm5VSm7OyshopjhDu4Z8b/0kL/xbc0vEWo6MIIVxUUdEulPLEx6fhw0q2ZPcCRynlBYwEvrrQ17XW72ute2qte4aFyQRIIWzlQM4Bvt//PQ/0eAAvs5fRcYQQLqqwcCe+vvGYTJ5GRzlHY/TgXA9s0Vofb4RnCSFOm7lxJmaTmQd6PGB0FCGEC7P1CipbaYwC51YuMjwlhLCPwvJCPtr2EWM6jiE8INzoOI1KKTVFKbVbKbVLKTVHKSUHbwlhJ5WV+ZSVHbLpCipbsWuBo5TyBQYDC+z5HCHEuT7b8Rn5ZfluN7lYKdUKmAT01Fp3AsxUzwEUQtjB/45ocLweHLse1aC1Lgaa2vMZQohzVVmreOuXt+ge3p0rIq4wOo4RPAAfpVQF4AtkGpxHCJdlizOo7EV2MhbCxSzYu4Ck7CSe6veU2507pbXOAF4HDgNHgTyt9fKzr5GVm0LYTlHRTsxmfyyWNkZHOY8UOEK4EK01L/30Eu2btueWDu63NFwpFQKMAqKBloCfUuqcA7hk5aYQtnNmgrEjvpmSAkcIF/L9ge/Zfnw7U6+aitlkNjqOEQYBqVrrLK11BdXz/640OJMQLklrTWHhToccngIpcIRwGWd6byKDIrmt821GxzHKYeAKpZSvqn5LORDYa3AmIVxSeflxKiuzHXKCMUiBI4TL+OHQD/x85GeevPJJPM2OteFWY9FabwDmA1uAnVS3ce8bGkoIF3XmiAZH7cGx6yoqIUTjeemnl2ju15x7ut1jdBRDaa1fAF4wOocQrs6Rl4iD9OAI4XSmrUg+73MbMzay8uBKnuj7BD6ePgakEkK4m6KinXh6NsfLyzEn60uBI4STmb5q/3mfe/mnlwmxhPBgzwcNSCSEcEeOekTDGVLgCOHkth/bzrdJ3zKpzyQCvAOMjiOEcANaWykq2o2/v2POvwGZgyOEU5i2IvmcnpuoqUsAmDwwluUnphJiCWFyn8lGxRNCuJmSkoNYrcUO3YMjBY4QTmDK4DimDI4DqoubtFeHA7Dy4Eoe/3kpbwx5gxCfECMjCiHciCMf0XCGDFEJ4aSs2sqTK54kKjiKh3s9bHQcIYQbObNE3Ne3o8FJLk56cIRwMpMHxgIwZ+ccth7byuc3f463h7fBqYQQ7qSwcDsWSzs8PPyNjnJR0oMjhJOZMjiO0spSnl39LN3DuzO+03ijIwkh3ExR0Q78/bsYHeOSpMARwgnN3DiTQ3mH+Mfgf2BS8s9YCNF4qqqKKCk5IAWOEMK2ckpyePGnF7k+5noGRA8wOo4Qws0UFu4EtBQ4QgjbeunHl8grzeO1Qa8ZHUUI4YaKirYD4OcnBY4Qwka2H9vO9A3TubfbvXRu7rjLM4UQrquwcDtmcyAWSxujo1ySFDhCOIkqaxX3/+d+mvo25bXB0nsjhDBGYeF2/P0TUUoZHeWSpMARwkm8s+kdNmVuYtrQaTTxaWJ0HCGEG6o+omGHww9PgRQ4QjiF9Px0nln9DEPaDeHWTrcaHUcI4aZKS1Opqip0+AnGYOcCRykVrJSar5Tap5Taq5Tqa8/nCeGqHv3+UaqsVbw7/F2H7xYWQriuwsIdAE5R4Nh7J+PpwFKt9WillBfga+fnCeFyFu5byMJ9C3l14Ku0DWlrdBwhhBsrLNwOmBz6kM0z7FbgKKUCgWuAuwG01uVAub2eJ4QryivN45HvHiGxeSKP933c6DhCCDdXVLQdH59YzGbH76+w5xBVWyALmK2U2qqU+lAp5ffbi5RSE5VSm5VSm7OysuwYR1zItBXJRkcQF6G1ZuLiiRwrPMYHIz7A0+xpdCQhhJurXkHl+MNTYN8CxwPoDryrte4GFAFTf3uR1vp9rXVPrXXPsLAwO8YRFzJ91X6jI4iLmLV1FvN2z+Nv1/2N3q16Gx1HCOHmKivzKS1Nxd8/0egotWLPAicdSNdabzj95/lUFzxCiMvYfWI3k76fxKC2g3jqqqeMjiOEEDUTjJ1hiTjYcQ6O1vqYUuqIUqq91joJGAjssdfzRO1NW5F8Ts9N1NQlAEweGMuUwXFGxRKnlVSUMG7+OAK8A/j0pk/lME0hhEM4c0SDswxR2XsV1aPA56dXUB0EJtj5eaIWpgyOqylkoqYuIe3V4QYnEmebsmwKu7N2s/S2pbTwb2F0HCGEAKrn33h4hODtHWF0lFqxa4Gjtd4G9LTnM4RwJfN2z+Nfv/6Lp/o9xdCYoUbHEUKIGoWFO/D37+I0e3FJ37ebmzwwttGeJSu2Lm1jxkbuXng3fSP68rfr/mZ0HCGEqKF1FUVFO51m/g3Yf4hKGKSwvJA9WXvYdWIXe7L2kFOSQ35Zfs1HpbUSfy9/ArwD2PR1AIHegbQJakNMkxhimsTQrkk7Ar0DbZpp+qr9MsfnIlJzUxkxZwQt/FuwcPxCWRIuhHAoJSUpWK3FTjP/BqTAcRmF5YUsT1nOoqRF/HT4Jw7mHqz5msXDQqhvKIHegQR6BxJsCcZsMlNQVkB6fjoFZQXkleVxsvjkOfeMCIygT6s+XBFxBVdEXEH38O74ejr+5k7OJqckh+s/v55KayXf3/Y9zfyaGR1JCCHOUb2DMU6zRBykwHFqReVFzNk1hwV7F7AqdRXlVeWEWEIYED2Au7vcTefmnenUrBPRwdGYTebL3q+grICDuQc5kHOA/Tn72XF8B7+k/8LXe78GwMPkwZWtr2RYu2FcH3s9XZpffixWVmxdWlllGTd9eROpp1JZecdK2oe2NzqSEEKcp7rAMePrm2B0lFpTWmujM9To2bOn3rx5s9ExHN6xwmPM2DiDdze/S05JDu1C2jGq/ShGth9Jv8h+eJhsW7eeKDrBhvQNrDu8juUHl7Pt2DYAWvi34IaYGxjdcTQD2w7Ey+x1yfvIiq1zWbWV2xfczpxdc5h7y1zGdRpndKQGU0r9qrV2moUF0uYIUTs7d46gpOQgvXvvNjrKeS7W7kgPjhNJyUnh5Z9e5rOdn1FRVcHv4n/HE32f4MrWV9p1Vnszv2aMaD+CEe1H8BqvcazwGMtTlvP9ge+Zv3c+H237iCDvIEbFj2JMxzEMaTfkssWOu6u0VjLh2wnM2TWHVwe+6hLFjRDCdRUWbico6CqjY9SJFDgOZtqK5POGbkorS3lt3Wu8su4VTMrEfd3u47ErHiO2aeOtgDpbC/8W3NnlTu7scidllWWsPLiS+Xvns3DfQv69/d809WnK+E7jubPLnfRq2aum+GrMFVtw4Z+lIyirLOPWr2/lm33f8NKAl2SnYiGEQysvz6Ks7Aj+/t2MjlInskzcwfz2bKjlKcvp/G5n/vzDn/ld/O84MOkAM4fPNKy4+S1vD2+Gxw1n9qjZHP/jcZb8fgmD2w1m1tZZ9PmwD/Ez43n5p5fJyM9o9GLDEc/ZKq4oZuTckXyz7xumD5vOM1c/Y3QkIYS4pPz86hOXAgP7GJykbqTAcVB5pXnc+vWtDP1sKArF8tuXM3f0XFoGtDQ62kV5mb24IfYG5twyh2NPHGPWyFm08G/Bs6ufJfKtSG784ka+2fsNFVUVRkc1RF5pHkM/G8rKgyuZNXIWk/pMMjqSEEJcVkHBBqqsJgICehgdpU5kkrED+O1Ko3KVRpbXy1jNJ3j+mj/x1FVPYfGwGJiwYQ7kHGD21tl8vP1jMgsyaebXjLu63MV93e8jrqlte3V++7M8w+hVW3uy9jB63mj25+zn85s/Z2zCWMOy2JNMMhbCdZxpT//Y8zkCvU7x/M//BIxvT3/rYu2OFDgOJuzZpyj2m0GgdyDzRs/j6jZXGx3JZiqtlSw7sIxZW2exKGkRVbqKa9tcy/3d7+eWjrfYvIhzlFVbn+/4nImLJ+Lv5c+cW+YwIHqA0ZHsRgocIVyL1lbWrWvCytS+/OXO742Oc0EXa3dkiMpBVFRV8Piyxznp9Xe6tejGlolbXKq4gep9dIbHDWfBuAWkP57OKwNfIT0/ndu/uZ2Wb7Tk0e8eZfux7UbHtJnSylIeXPwgt39zOz3Ce7D1ga0uXdwIIVxPcXESVVV5pOQ53x5dUuA4gNLKUm6edzPTfpnG1eF3sPqu1YQHhBsdy65a+Ldg6lVTSX40mVV3rmJYzDA+2PIBXf/VlZ7v9+TdTe+SW5LboGc09qqts23O3MyVs66sOThz9V2rHXr+lBBCXMjgTuDgAAAgAElEQVSZCcb94ocYnKTupMAxWHFFMaPmjmJx8mLeueEdfpz4b7faQ8akTAyIHsAXt3xB5hOZvD3sbSqtlfzhuz/Q4o0WjPlqDP9J+k+9JiYbMUacVZTF/Yvup/cHvcksyGTR+EW8OuhVm2++KIQQjaGgYANmcyAPDhpmdJQ6k1bXQEXlRYyYM4K1aWv5aORHTOg2wehIhmri04RH+zzKI70fYeuxrfx7+7/5YucXzN8znzDfMMYljGNswlj6RfbDpByrNq+0VvLe5vd4bs1zFJYXMuWKKTx/7fMEWYKMjiaEEPWWn7+BgIBeKAdrc2tDChyD5JflM/yL4fx85Gc+velTbku8zehIDkMpRffw7nQP784/Bv+DpQeW8u8d/+bDrR8yY9MMwv3DGd1xNKM7jubK1lca2juSW5LLrK2zmLFxBofyDjGo7SDeHvY2HcI6GJZJCCFsoaqqmMLCHURGOudmpFLgGKC4ophhnw1jU+Ym5t4ylzEJY4yO5LA8zZ41x0QUlheyOHkx83bP44MtH/DPjf8k2BLMkHZDuCHmBobFDKO5f3O7Z9JasztrNzM3zuTfO/5NcUUx17S5hrevf5sRcSPsemyGEEI0loKCLUAVgYFXGB2lXqTAaWRWbeWOb+7gl/RfmD92Pjd3uNnoSE7D38uf8Z3GM77TeArKClh6YCnfH/ie7w98z7zd8wBICEvgytZX1nzENom1ScFRUlHCD4d+4Lv93/Hd/u9IyU3B2+zN7zv/nkl9JtG1RdcGP0MIIRxJfv4vgPPtYHyGFDiN7OmVT7Ng7wKmDZ0mxU0DBHgHMCZhDGMSxmDVVrYf2873B75n3eF1fLXnKz7Y8gEAwZZg4kPjad+0Pe2btieuaRzhAeEEW4IJ8g4i2BKMj6cPJRUlFFUUUVxRTGF5IQdzD7I3ay/7svexN2svO47voKSyBB8PHwZED2DKFVMYkzCGZn7NDP5JCCGEfRQUbMBiicLLyznbOSlwGtGHWz7k7z//nYd6PsTkPpONjuMyTMpEt/BudAuvPgjOqq3sO7mPn4/8zObMzSRnJ7Pi4Ao+2f5Jne8d7h9OfGg8E3tMZFjMMK5tcy0+nj62fgnChpRSwcCHQCdAA/dordcbm0oI55Ofv4HAwCuNjlFvUuA0klUHV/HQkocY2m4ob1//tszTsCOTMtExrCMdwzpyX/f7aj5fUFbA/pz9ZBVlcar0FHlleZwqPUVJRQm+nr74evri5+WHr6cvbYLa0D60PcGWYANfiain6cBSrfVopZQX4Gt0ICGcTVnZUcrKjjjt8BTYucBRSqUBBUAVUOlMW7jb0r6T+7hl3i3Eh8Yzb8w82RPFIAHeAXQP7250DGFHSqlA4BrgbgCtdTlQbmQmcb5pK5Id6iwjcb7/nSDunBOMoXF6cK7TWp9shOc4pNLKUsbNH4en2ZPFty4m0DvQ6EhCuLK2QBYwWynVBfgVmKy1LjpzgVJqIjARIDIy0pCQrqaiIofy8mNYraWnP8oAsFja4O3dGpPJ85zrp6/aLwWOgyso2IBSnvj7dzM6Sr1JV4KdTV05lR3Hd7D41sW0CW5jdBwhXJ0H0B14VGu9QSk1HZgKPHfmAq31+8D7UH3YpiEpnVxFxSny8n4gN3cNp06toahoxyWuNmOxRGKxtMXfvwtBQf0I9CputKyifvLzf8Hfvwtms20PQW5M9i5wNLBcKaWBf51uWM7hyu+mvtv/HdM3TGdS70kMjzP+VGsh3EA6kK613nD6z/OpLnCEDZSWHiYt7S8cO/YJUIXJZCEwsB/R0S/i4xODyeSDyeSNyWRBayulpWmUlh6kpOQgqcd24ZX9TzzNb/L2AJi75EmScxOIbDGM319zh9Ou1HFFWldRULCZ5s3vMjpKg9i7wOmntc5USjUDViil9mmtfzz7Ald9N3Ws8Bh3L7ybxOaJvDb4NaPjCOEWtNbHlFJHlFLttdZJwEBgj9G5nF15+UkOH36FjIyZALRq9TBhYbcQGNgHk8m7Vvfo2BGs1jIKCrYwdc4H/KFfDhFBP1JZuZKff/4j/v7dCAkZQkjIIIKC+mE2y2pFoxQV7aGqqtCpJxiDnQscrXXm6f89oZT6BugN/Hjp73J+Vm3lroV3UVBewNpb1mLxcN4uPiGc0KPA56dXUB0E3PuQtwbQWpOePo20tL9QVVVIixZ3ERX1ZyyW+vW2m0zeBAX15fu0HN59cPjpnoIt5OYuJydnOenpb3DkyGso5U1Q0JUEBw8gJGQA/v7dnXqoxNk4+wZ/Z1y2wFFKPQJ8rrXOrcuNlVJ+gElrXXD6v4cAf61fTOcy/ZfpLE9ZzrvD36VjWEej4wjhsOrbvlyK1nob4JYrNm3Jai0jKel+jh//lCZNhtOu3d/x87NNezZ5YCwASpkJDOxFYGAv2rR5lsrKQvLyfiQ3dzWnTq0iLe050tKeOz3ZtQsBAX0IDOyDn19nfH1jMZv9LvssWbFVd7m5q/DyaoGPT6zRURqkNj04LYBNSqktwEfAMq11bYaSmgPfnN7vxQP4Qmu9tN5JncT+7P08veppRrUfxQM9HjA6jhCOrr7ti7Cj8vKT7N59E3l564iK+itt2vzJpnt3Xazg8PDwp2nTG2ja9IaaHHl56ygo2EB+/gaOH/+EzMyZNdd7ebXC1zcOi6UtXl7N8PQMxdMzDE/PMMxmP0wmC19v2MAD/UApL5Qy13yA+fR8IR9MsnVHDa2ryM1dTmjoKKffr03Vpi1R1a9yCNVdvT2BecAsrXWKLcP07NlTb9682Za3bFRaa4Z8NoRNGZvY98g+Wvi3MDqSEI1KKfVrXfe7aqz25UKcvc2xh6KifezcOZyysgw6dPiEZs3GGR2phtZVFBfvo6hoLyUlyRQXJ1NSkkxpaRoVFVloXVmv+yrlicnkg4dHEN7eEac/WuPt3ZqAgB4EBPRymyGyvLxf2Lq1Lx07znWov/tLuVi7U6uyVWutlVLHgGNAJRACzFdKrdBaP2nbqM5r7q65rDy4kpk3zJTiRohakvbFceTnb2b79kGYTN507bqWoCDH2uRNKTN+fgn4+SWc9zWtNZWVeVRUZPHpz1tYuOUAXuYKPE3leJgq8DBVMKRDGIM6hqJ1FVpXnd6zpwSrtYSqqhIqK3MpK8ugsHA72dlLsFqLTz/Xi4CAXgQFXUWTJoMJDr4OpUyN/fIbRU7OUsBESMggo6M02GV7cJRSk4C7gJNUn++yUGtdoar/dvdrrdvZKowzv5s6VXqK+BnxRAZFsv7e9ZhNZqMjCdHo6tqD05jty4U4c5tja6WlR9iypTdKedOt2w9YLK6xb1fU1CWkvVr3bTq01lRUZJGfv4G8vHXk5f1EQcFmtK7AYmlLePj9hIdPwMuruR1SG2fLlr5orenR4xejo9RaQ3pwQoGbtdaHzv6k1tqqlLrRVgGd3TOrniGrOIvvbvtOihshak/aFwcwfcUWrgyeQFVVMd27r3SZ4qYhlFJ4eTUjNHQEoaEjAKiqKuHkyYVkZv6L1NSnSUt7jtDQm4iKegE/vwSnn9BcUZFNfv5G2rR57vIXO4HL9rFprZ//beNz1tf22j6S89mQvoH3Nr/HpN6T5KwjIepA2hfjaV2Fzn2QoqLdJCTMu+DwjzM7s2LLFsxmH5o3v5Vu3dbSu/c+WrWaRG7uCjZv7sqBA0/wr7XbbfYsI+TmrgSsNGkyzOgoNuGag4iNqNJayQOLH6BlQEv+ep1brIIXQriQAweeoGuzTcTG/pMmTYYaHcfm7NWj4uvbnpiYN+jdez8tWkwgPX0ar179AMeOfYazLgTMyVmKh0cIgYG9jI5iE1LgNNB7m99j+/HtTB82nQDvAKPjCCFErUxbkcyEmQ+TkTGd5Wkj6ffPSKKmLmHaimSjozmVmT/kMHT2KP7y8xtkl4axb98dvPN1N/65Yr3R0epEa01OzlJCQoacXkbv/GTxfwMUlBXw1x/+Sv+o/tzc4Waj4wgDOfvYu3A/E6+0stlrFiEhw5iz9N56TcQV1T1E1f/2hxM9NYZ1j2ZgMk3B03Mc+fkLCAx0jj0ni4p2UF5+zGWGp0B6cBrkjfVvkFWcxWuDXnP6DZFEw0xftd/oCELUmtZV7Nt3N2azH+3bz0bjGu/YjaYx0arVg3Trtg5QbN16FUePfmx0rFqpXh4OTZoMMTiJ7UiBU0/HC4/zxvo3GN1xNL1b9TY6jhBC1NqRI29QULCB2NgZeHu3sOlEXHd25ucYENCDHj02ExTUj6SkCSQnP4LVWmFwukvLyVmKn18i3t4tjY5iMzJEVU8v/vgiJRUlvDTgJaOjCINMW5F8Ts9N1NQlQHUjJ8NVwlEVFe0lNfV5QkNvplmz8YD9JuK6m7N/jl5eYSQmLuPgwamkp79BeXkmHTvOxWTyMjDhhVVWFpCX918iIqYYHcWmpMCph5ScFN779T3u634fcU2lYXBX/xt7r/9mYkI0Jqu18vTQlD9xce/I0LqdmUwexMS8jsUSyYEDk9m162YSEuY73LEPp06tQesKl5p/AzJEVS/PrXkOL7MXL1z7gtFRhBCi1o4ceZ2Cgo3Exc10uR14HVlExCTi4t4jJ2cJu3aNoqqq2OhI58jJWYrJ5EdQUD+jo9iUFDh1tOXoFubsmsNjfR4jPCDc6DjCQcgcBuHoiouTSEt7gbCw0YSFjTU6jttp2fIB2rf/iNzcFezceSOVlYVGRwJAayvZ2d8REjLQIYfPGkIKnDp6ZtUzNPFpwpP95AxA8T8yh0E4upSUP2IyWYiNnSFDUwYJD59Ahw6fcurUD+zcOZyqqhKjI5Gbu5qyskOEhY02OorNSYFTB5syNrEsZRlP9XuKIEuQ0XGEEKJWcnJWkp29mDZtnpWhKYM1b34bHTp8Rl7eT+zZMx6rtdLQPJmZ7+Lh0ZSwsDGG5rAHKXDq4JV1rxBsCeahng8ZHUUIIWpF6ypSUp7AYommVatJRscRQPPmtxIbO4Ps7EUkJd2H1lZDcpSVZXDy5LeEh9/jcBOfbUFWUdXS3qy9fLPvG/509Z/kSAYhhNM4enQ2RUU76Nhxnkv+n5izatXqD1RUZJGW9mc8PUNp1+4fjT50mJn5AVBFy5YPNOpzG4sUOLX02n9fw8fDh0l95B2QEMI5VFYWkJr6JwID+7nkHAtn16bN81RUnCQ9/Q28vMKIjHyq0Z5ttVZw9OgHhIQMxcenXaM9tzFJgVMLh/MO8/nOz/lDzz8Q5hdmdBwhhKiVw4dfpaLiOJ07L5KJxQ5IKUVMzHQqKrI5eHAqnp6hhIff2yjPzs7+D+XlmcTFvdMozzOCFDi18PrPrwPwxJVPGJxECCFqp7T0EEeOvEGzZrcRGCjHyTgqpUzEx39MRUU2SUkT8fRsRmjoCLs/NzPzXby9W9OkietuUCqTjC8jqyiLD7d8yO2JtxMZFGl0HCGEqJWDB59FKUXbtq8YHUVchsnkRULC1wQE9GDPnrHk5f1s1+cVFyeTm7uS8PCJmEyu289h9wJHKWVWSm1VSi2297PsYfqG6ZRWlvJUv8YbGxVCiIYoKtrDiRNf0KrVJCyW1kbHEbXg4eFP585L8PZuzc6dN1JUtNtuz8rMfA+lPAgPv89uz3AEjdGDMxnY2wjPsbn8snxmbJzBTR1uIj403ug4QghRK2lpf8Vs9qN16/8zOoqogzMHdJpM3uzYMYzS0iM2f0ZVVQnHjn1MaOhNeHu3sPn9HYldCxylVAQwHPjQns+xl4+2fkReWR5T+001OooQQtRKYeEusrLm0arVJLy8Qo2OI+rIxyeazp2/p7Iyn+3bB1JWlmHT+x89+j6Vlbm0bPkHm97XEdm7B+ct4EngorsYKaUmKqU2K6U2Z2Vl2TlO7Vm1lZmbZtI3oi+9WvUyOo4QQtTKoUN/wWz2p3VrWRThrAICupKY+B3l5UfZtq0/paXpNrlvYeFOUlKeokmTYQQHX2uTezoyuxU4SqkbgRNa618vdZ3W+n2tdU+tdc+wMMdZgr0iZQUHcg7wSO9HjI4ihBC1Uli4g6ys+URETMbTs4nRcUQDBAX1IzFxOeXlx08XOQ0brqqqKmbPnlvx8AgmPv5jt9g2wJ49OP2AkUqpNGAuMEAp9Zkdn2dTMzfNpJlfM27pcIvRUYQQolbS0v6M2RxIRMTjRkcRNhAU1JcuXVZQUZHFtm3XUlp6qN73Skl5guLi3XTo8G+3OY/MbgWO1vpprXWE1joKGA+s1lrfbq/n2VJqbiqLkxdzf/f78fbwNjqOEEJcVkHBNk6e/IaIiCl4eoYYHUfYSGBgH7p0WUllZS7btvWnoGBLne+RlfUNmZnv0br1H2nSZIgdUjom2QfnAt7b/B4mZeKBHq55PocQwvVU994EERHxmNFRhI0FBvaiS5eVWK3lbNlyBYcOvYLWVUxbkXzZ7y0tPUJS0r34+/cgOvqlRkjrOBqlwNFar9Va39gYz2qo0spSZm2dxaj4UbQOkv0jhBCOr6BgK9nZ39K69eN4egYbHUfYQUBAD3r12klo6E2kpj7Dtm39+fznny56vdaa3NxV7Nw5Aqu1nI4d52AyeTViYuO57haG9fTlri/JLsnm4V4PGx1FCCFq5fDhVzCbA2nVSg4DdmWenk3o2HEux4+PYP/+h/lbv0c5eDCFwMC+BAb2xsurOVprcnK+49ChF8nP/wUvr3A6dvwCX99Yo+M3OilwfmPGphl0CO3AdVHXGR1FCCEuq6hoH1lZ84mMfFp6b9zAWyv3M31VCE0tb3F3wgy80l7DbKreicXbuw1msy/FxXvx9m5DbOw7tGgxAbPZYnBqY0iBc5aNGRvZnLmZGdfPcIsldEII53fkyGuYTBaZe+MmpgyOY8rgOACipjYj5aXrKCjYQkHBRvLzN1Benkn79rNp3vw2TCZPg9MaSwqcs8zcNBN/L3/u6HKH0VGEEOKySksPcfz4Z7Rs+TBeXo6zj5hoPGazL8HBVxEcfJXRURyOrKI6La80j692f8XtnW8n0DvQ6DhCCHFZR468DijZtdhNTR7ofvNq6kIKnNO+2vMVJZUlTOg2wegoQghxWeXlxzl69EOaN79TTgx3U2eGqsSFSYFz2sfbPqZDaAd6tTz/3Kna7DUghBCN6ciRaVit5URGPmV0FCEckhQ4wP7s/fz3yH+5u+vdF5xcPH3VfgNSCSHEhVVU5JKZ+Q7Nmo11y+W/QtSGFDjAJ9s/waRM3J7oFCdJCCHcXEbGTKqqCoiMfNroKEI4LLdfRVVlreKT7Z8wtN1QWga0rPn8tBXJ5/TcRE1dAlRP6pJxTyGEUaqqikhPf4umTW/E3z/R6DhCOCy3L3DWpK0hPT+dN4a8cc7nz91rYAlprw43Ip4Qoo6UUmZgM5DhLEfE1EVm5gdUVmYTGfms0VGEcGhuP0T18baPCbYEM7L9SKOjCCFsYzKw1+gQ9mC1lnHkyOsEB19HUNAVRscRwqG5dYGTV5rHgr0LuLXTrVg8Lr6Vtew1IIRzUEpFAMOBD43OYg/Hjn1KeXkGkZHPGB1FCIfn1gXOmb1v7u569yWvkzk3QjiNt4AnAevFLlBKTVRKbVZKbc7Kymq8ZA1ktVZy+PCrBAT0JCRkoNFxhHB4bl3gXGrvGyGEc1FK3Qic0Fr/eqnrtNbva617aq17hoU5z/EGWVnzKS1NITLyGTkrT4hacNsC50DOAf575L/c1eUuaSyEcA39gJFKqTRgLjBAKfWZsZFsQ2vN4cMv4+vbgdDQUUbHEcIpuG2BM2fnHBSK2xJvMzqKEMIGtNZPa60jtNZRwHhgtdbaJTa3ys5eQlHRTiIjn0Ypt222hagTt/2XMm/PPK6KvIqIwAijowghxEVV9968hMUSRbNm442OI4TTcMsCZ0/WHnad2MXYhLFGRxFC2IHWeq2r7IFz6tRq8vN/oXXrJzGZPI2OI4TTcMsCZ97ueSgUt3S4xegoQghxUVpr0tL+jJdXK8LD7zE6jhBOxe12MtZaM2/3PK5pcw3hAeFGxxFCiIvKzV1FXt46YmNnYjJ5Gx1HCKditx4cpZRFKbVRKbVdKbVbKfUXez2rLnZn7Wbvyb2MSxhndBQhhLioM7033t4RhIffa3QcIZyOPXtwyoABWutCpZQnsE4p9b3W+hc7PvOy5u2eh0mZuLnDzUbGEEKIS8rNXUl+/n+JjX1Hem+EqAe79eDoaoWn/+h5+kPb63m1obXmy91f0j+qP839mxsZRQghLkprzc/bnsLbu7XMvRGinuw6yVgpZVZKbQNOACu01hsucE2jbZu+4/gOkrOTGdtRVk8JIRxXbu4K/E1biYx8RnpvhKgnuxY4WusqrXVXIALorZTqdIFrGm3b9Hm752FWZhmeEkI4rDNzb7JLwggPn2B0HCGcVqMsE9danwLWAsMa43kXycCXu79kQPQAwvyc5/wZIYT7mLYimRFv/I38/PX85+BY2j6zkqipS5i2ItnoaEI4HbtNMlZKhQEVWutTSikfYBDwmr2edzlbj20lJTeFqVdNNSqCEEJc0mODYrg65FsqKiL5KX0Qaa8ONzqSEE7LnquowoFPlFJmqnuK5mmtF9vxeZd0ZnjqpvibjIoghBCXdPz4FxQWbiE+/lOqtOxaLERD2K3A0VrvALrZ6/51obXm671fM7DtQJr6NjU6jhBCnKeqqoTU1Gfw9+9B8+a/Z/LAA0ZHEsKpucVRDftO7uNAzgF+1/53RkcRQogLSk+fTlnZEdq1ex2lTEwZHGd0JCGcmlsUOIuSFgEwov0Ig5MIIcT5ysuzOHz4ZZo2HUlISH+j4wjhEtyjwEleRPfw7kQERhgdRQghzpOW9heqqopp29awdRhCuByXL3BOFJ1g/ZH1jIwbaXQUIYQ4T3FxEpmZ79Gy5QP4+cUbHUcIl+HyBc6S5CVoNCPbS4EjhHA8KSlPYTb7EhX1gtFRhHApLl/gLEpeRERgBF1bdDU6ihBCnCMnZxnZ2d8SGfk0Xl7NjI4jhEtx6QKnpKKE5SnLGRk3EqWU0XGEEKJGVVURyckP4usbT+vWjxsdRwiXY8+N/gy3OnU1xRXFMjwlhHA4qakvUFqaRteuP8qBmkLYgUv34CxKWoS/lz/9o/obHUUIIWoUFGwhPX0a4eETCQ6+2ug4Qrgkly1wrNrKf5L/w7CYYXh7yLsjIYRjsForSUq6Hy+vZrIsXAg7ctkhql8zf+Vo4VFZHi6EcCgZGdMpLNxCx45f4ekZbHQcIVyWy/bgLEpahEmZuCH2BqOjCCEEACUlqaSmPk/TpiMIC7vF6DhCuDTXLXCSF3FV5FVyuKYQwiFobSUp6V6UMhEbO1NWdgphZy5Z4KSdSmPH8R0yPCWEcBhHjrzBqVNriIl5C4ultdFxhHB5LlngfL//ewBujLvR4CRCCAEFBVtJTX2W0NCbadHiHqPjCOEWXLLAWZayjDZBbYhrGmd0FCGEm6uqKmbv3tvw9Ayjffv3ZWhKiEbicgVORVUFq1NXM7TdUGlIhBCGS0l5kuLivcTHf4ynp8wJFKKxuFyBsz59PQXlBQyNGWp0FCGEm8vO/o7MzJlEREyhSZPBRscRwq24XIGz7MAyzMrMwOiBRkcRQrix8vIs9u27Bz+/zkRHv2x0HCHcjstt9LcsZRlXRFxBkCXI6ChCCDeltSY5+UEqK3Pp0mUFZrPF6EhCuB2X6sHJKspiy9EtDG0nw1NCCOMcP/45J08uIDr6b/j7dzY6jhBuyW4FjlKqtVJqjVJqr1Jqt1Jqsr2edcaKgyvQaJl/I4QwTGlpOvv3P0JgYD9at37C6DhCuC17DlFVAk9orbcopQKAX5VSK7TWe+z1wGUpy2ji04Qe4T3s9QghhLgorTVJSfegdQXx8R+jlNnoSEK4Lbv14Gitj2qtt5z+7wJgL9DKjs9jecpyBrcdjNkkjYoQovFlZr5Lbu4K2rV7A1/fGKPjCOHWGmUOjlIqCugGbLjA1yYqpTYrpTZnZWXV+xk7ju/gWOExmX8jhDBEcfF+UlL+SEjIUFq2fMDoOEK4PbsXOEopf+Br4DGtdf5vv661fl9r3VNr3TMsLKzez1mWsgyAIe2G1PseQghRXykpj6OUF/Hxs2STUSEcgF0LHKWUJ9XFzeda6wX2fNaylGV0ataJVoF2GwUTQogLystbT3b2YiIjn8TbW9ogIRyBPVdRKWAWsFdr/aa9ngNQVF7EusPrZHhKCGGI1NRn8fRsRqtWk4yOIoQ4zZ49OP2AO4ABSqltpz9usMeD1qatpbyqXAocIUSjy81dxalTa4iMfBoPD3+j4wghTrPbMnGt9TqgUQail6csx8fDh6vbXN0YjxNCCKB69ebBg8/i7R1By5YPGh1HCHEWlziqYWXqSq5pcw0WD9kOXQjReLKzF1NQsIG4uH/JcQxCOBinP6rheOFx9mTtYUD0AKOjCCHciNZWUlP/hMXSjhYtJhgdRwjxG05f4KxNWwtA/6j+huYQQhivMY+Iycr6iqKiHURF/RmTydNejxFC1JPTD1GtSVtDgFcA3cO7Gx1FCGG8Rjkiprr35gV8fTvSvPmttry1EMJGnL4HZ03aGq5pcw0eJqev1YQQDdRYR8QUFGyipCSJ1q3/T86bEsJBOXWBk1mQSXJ2MtdFXWd0FCGEg7nYETG2OB4mK+trlPIgNHRUg3MKIezDaQucaSuSWZO6BoDroqXAEUL8z6WOiGno8TBaa7KyFhAcPBBPzxAbJRZC2JrTFjjTV+1nTdoagi3BdGnexeg4QggHYe8jYoqKdlBamkJY2M22vrUQwoactsCB6vk317a5FrNJxsCFEI1zRExW1gLARGjo7+xxeyGEjTjVzNxpK5KZvmo/AJXqBIdyD5J7fCDTmiYzZXCcwemEEA7gzBExO5VS205/7hmt9Xe2ekBW1tcEBV2Nl1czW91SCGEHTlXgTBkcV1PIhD77OABrJ/6eylIAAAZZSURBVD1CYnMpboQQ9j8iprg4ieLi3cTETLfXI4QQNuK0Q1Slpp009WlKp2adjI4ihHAT1cNTEBoq82+EcHROW+B4+e3l2qhrMSmnfQlCCCeTlfU1AQF9sFgijI4ihLgMp6wOUnNTyS3LkP1vhBCNprT0EIWFv8rqKSGchFMWOGvSTu9/IwWOEKKRnBmeCgu7xeAkQojacNoCp5lfMzqGdTQ6ihDCTWzf/zl+fl3w8WlndBQhRC04XYGjtWZN6hr6R/WnessLIYSwr7KyY/iqLTI8JYQTcboC50DOATIKZP6NEKLxnDy5EJPSMjwlhBNxqn1wzpjQdQKD2g4yOoYQwsWd2Vy0c+hxejYfwt1LU4E0Jg+Mlc1FhXBwSmttdIYaPXv21Js3bzY6hhCinpRSv2qtexqdo7bq0uZETV1C2qvD7ZxICFFXF2t3nG6ISgghhBDicuxW4CilPlJKnVBK7bLXM4QQorFMHhhrdAQhRB3YswfnY2CYHe8vhBCNRubcCOFc7FbgaK1/BHLsdX8hhBBCiIuROThCCCGEcDmGFzhKqYlKqc1Kqc1ZWVlGxxFCCCGECzC8wNFav6+17qm17hkWFmZ0HCGEEEK4AMMLHCGEEEIIW7PnMvE5wHqgvVIqXSl1r72eJYQQQghxNofayVgplQUcquXlocBJO8Yxiqu+LpDX5qzq8traaK2dZqxZ2pwa8tqck7y2ahdsdxyqwKkLpdRmZ9oSvrZc9XWBvDZn5cqvrS5c+ecgr805yWu7NJmDI4QQQgiXIwWOEEIIIVyOMxc47xsdwE5c9XWBvDZn5cqvrS5c+ecgr805yWu7BKedgyOEEEIIcTHO3IMjhBBCCHFBUuAIIYQQwuU4XYGjlBqmlEpSSh1QSk01Oo+tKKVaK6XWKKX2KqV2K6UmG53J1pRSZqXUVqXUYqOz2JJSKlgpNV8pte/0319fozPZglJqyunfxV1KqTlKKYvRmYwi7Y5zkjbH+diy3XGqAkcpZQZmAtcDHYFblVIdjU1lM5XAE1rrDsAVwMMu9NrOmAzsNTqEHUwHlmqt44EuuMBrVEq1AiYBPbXWnQAzMN7YVMaQdsepSZvjRGzd7jhVgQP0Bg5orQ9qrcuBucAogzPZhNb6qNZ6y+n/LqD6F7aVsalsRykVAQwHPjQ6iy0ppQKBa4BZAFrrcq31KWNT2YwH4KOU8gB8gUyD8xhF2h0nJG2O07JZu+NsBU4r4MhZf07HRf4xnk0pFQV0AzYYm8Sm3gKeBKxGB7GxtkAWMPt0V/iHSik/o0M1lNY6A3gdOAwcBfK01suNTWUYaXeck7Q5TsbW7Y6zFTjqAp9zqXXuSil/4GvgMa11vtF5bEEpdSNwQmv9q9FZ7MAD6A68q7XuBhQBTj9HQykVQnUvRTTQEvBTSt1ubCrDSLvjZKTNcU62bnecrcBJB1qf9ecIXKjbXCnlSXUj87nWeoHReWyoHzBSKZVGdff+AKXUZ8ZGspl0IF1rfeZd73yqGx9nNwhI1Vpnaa0rgAXAlQZnMoq0O85H2hznZNN2x9kKnE1ArFIqWinlRfXko0UGZ7IJpZSiekx1r9b6TaPz2JLW+mmtdYTWOorqv7PVWmuX6A3QWh8Djiil2p/+1EBgj4GRbOUwcIVSyvf07+ZAXGQiYz1Iu+NkpM1xWjZtdzxsFqsRaK0rlVKPAMuonl39kdZ69/+3d8e6NkRhFIDXiuYKGp1Sr9XQeY8r0fMCaoknIQqVQqmUiOQSovUcCrIV57RC4t4zzu/7yqn+as2aPTN7bzzWebmb5DTJ57Yf99cer7VebzgTf+ZRkmf7m9/XJA82nuevrbXetX2Z5Cy7P20+ZPa28L8kd/gHjcuc5Pxzx1ENAMA4x/aKCgDgtxQcAGAcBQcAGEfBAQDGUXAAgHEUHABgHAUHABhHweHCtb3d9lPbk7ZX2n5pe2vruYCZZA6Jjf44kLZPkpwkuZzdOSpPNx4JGEzmoOBwEPstxd8n+Zbkzlrrx8YjAYPJHLyi4lCuJ7ma5Fp2T1UAF0nm/Oes4HAQbV8leZHkZpIba62HG48EDCZzOKrTxDlObe8n+b7Wet72UpK3be+ttd5sPRswj8whsYIDAAzkGxwAYBwFBwAYR8EBAMZRcACAcRQcAGAcBQcAGEfBAQDG+QmekHWMhc1GwAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Data points\n",
    "xval = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8 ])\n",
    "yval = np.array([1.5, 3.5, 5.5, 5.9, 5.7, 5.5, 4.5, 6.5, 8.5])\n",
    "# plotting coordinates\n",
    "xabs = np.linspace(0, 8, 50)\n",
    "\n",
    "# polynomial definitions\n",
    "def poly8(x, a0=0, a1=0, a2=0, a3=0, a4=0, a5=0, a6=0, a7=0, a8=0):\n",
    "    return a0 + a1*x + a2*x**2 + a3*x**3 + a4*x**4 + a5*x**5 + a6*x**6 + a7*x**7 + a8*x**8\n",
    "def poly7(x, a0=0, a1=0, a2=0, a3=0, a4=0, a5=0, a6=0, a7=0):\n",
    "    return a0 + a1*x + a2*x**2 + a3*x**3 + a4*x**4 + a5*x**5 + a6*x**6 + a7*x**7\n",
    "def poly6(x, a0=0, a1=0, a2=0, a3=0, a4=0, a5=0, a6=0):\n",
    "    return a0 + a1*x + a2*x**2 + a3*x**3 + a4*x**4 + a5*x**5 + a6*x**6\n",
    "def poly5(x, a0=0, a1=0, a2=0, a3=0, a4=0, a5=0):\n",
    "    return a0 + a1*x + a2*x**2 + a3*x**3 + a4*x**4 + a5*x**5\n",
    "def poly4(x, a0=0, a1=0, a2=0, a3=0, a4=0):\n",
    "    return a0 + a1*x + a2*x**2 + a3*x**3 + a4*x**4\n",
    "def poly3(x, a0=0, a1=0, a2=0, a3=0):\n",
    "    return a0 + a1*x + a2*x**2 + a3*x**3\n",
    "def poly2(x, a0=0, a1=0, a2=0):\n",
    "    return a0 + a1*x + a2*x**2\n",
    "def poly1(x, a0=0, a1=0):\n",
    "    return a0 + a1*x\n",
    "\n",
    "# First fitting with linear function\n",
    "p0=[0, 0] # initial values (2 expected)\n",
    "bounds=(-np.inf, np.inf) # bounds for initial values\n",
    "pars, cov = curve_fit(f=poly1, xdata = xval, ydata = yval, p0=p0, bounds=bounds) # fit itself\n",
    "[a0, a1] = pars # get the result\n",
    "# calculate residual and R-squared\n",
    "f = poly1\n",
    "residuals = yval- f(xval, *pars)\n",
    "ss_res = np.sum(residuals**2)\n",
    "ss_tot = np.sum((yval-np.mean(yval))**2)\n",
    "r_squared1 = 1 - (ss_res / ss_tot)\n",
    "\n",
    "# Do the same for poly2\n",
    "p0=[0, 0, 0]\n",
    "bounds=(-np.inf, np.inf)\n",
    "pars, cov = curve_fit(f=XX, xdata = xval, ydata = yval, p0=p0, bounds=bounds)\n",
    "[b0, b1, b2] = pars\n",
    "f = XX\n",
    "residuals = yval- f(xval, *pars)\n",
    "ss_res = np.sum(residuals**2)\n",
    "ss_tot = np.sum((yval-np.mean(yval))**2)\n",
    "r_squared2 = 1 - (ss_res / ss_tot)\n",
    "\n",
    "\n",
    "# Do the same for poly3\n",
    "p0=[0, 0, 0, 0]\n",
    "bounds=(-np.inf, np.inf)\n",
    "pars, cov = curve_fit(f=XX, xdata = xval, ydata = yval, p0=p0, bounds=bounds)\n",
    "[c0, c1, c2, c3] = pars\n",
    "f = XX\n",
    "residuals = yval- f(xval, *pars)\n",
    "ss_res = np.sum(residuals**2)\n",
    "ss_tot = np.sum((yval-np.mean(yval))**2)\n",
    "r_squared3 = 1 - (ss_res / ss_tot)\n",
    "\n",
    "\n",
    "# Do the same for poly8\n",
    "\n",
    "\n",
    "# Plot all fits along with data points\n",
    "fig1, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(8, 8), tight_layout=True)\n",
    "\n",
    "ax1 = plt.subplot(2,2,1)\n",
    "ax1.plot(xval, yval, '+', label = \"data point\")\n",
    "ax1.plot(xabs, poly1(xabs, a0, a1), '-r',label = f\"poly1 fit - R^2 = {r_squared1:.4f}\")\n",
    "\n",
    "ax2 = plt.subplot(2,2,2)\n",
    "ax2.plot(xval, yval, '+', label = \"data point\")\n",
    "ax2.plot(xabs, poly2(xabs, b0, b1, b2), '-k',label = f\"poly2 fit - R^2 = {r_squared2:.4f}\")\n",
    "\n",
    "ax3 = plt.subplot(2,2,3)\n",
    "ax3.plot(xval, yval, '+', label = \"data point\")\n",
    "ax3.plot(xabs, poly3(xabs, c0, c1, c2, c3), '-g',label = f\"poly3 fit - R^2 = {r_squared3:.4f}\")\n",
    "\n",
    "ax4 = plt.subplot(2,2,4)\n",
    "ax4.plot(xval, yval, '+', label = \"data point\")\n",
    "ax4.plot(xabs, poly8(xabs, d0, d1, d2, d3, d4, d5, d6, d7, d8), '-y',label = f\"poly8 fit - R^2 = {r_squared8:.4f}\")\n",
    "\n",
    "\n",
    "ax1.set_xlabel('x')\n",
    "ax1.set_ylabel('y')\n",
    "ax1.legend()\n",
    "\n",
    "ax2.set_xlabel('x')\n",
    "ax2.set_ylabel('y')\n",
    "ax2.legend()\n",
    "\n",
    "ax3.set_xlabel('x')\n",
    "ax3.set_ylabel('y')\n",
    "ax3.legend()\n",
    "\n",
    "ax4.set_xlabel('x')\n",
    "ax4.set_ylabel('y')\n",
    "ax4.legend()\n",
    "# ax1.set_yscale('log')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Curve fitting on actual data\n",
    "We work on **covid-19-pandemic-worldwide-data.csv** file recensing the cumulative number of dead people from covid-19 in France over the year 2020.\n",
    "\n",
    "We have to:\n",
    "- load the data\n",
    "- plot them\n",
    "- calculate the daily number of dead people\n",
    "- fit a reasonnable model to the data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "with open(\"data\\\\covid-19-pandemic-worldwide-data.csv\", \"r\") as csv_file:\n",
    "    csv_reader = csv.reader(csv_file, delimiter = ';')\n",
    "    linecount = 0\n",
    "    nblines = 10\n",
    "    print(f\"Overview of {nblines} first lines:\")\n",
    "    for lines in csv_reader:\n",
    "      if linecount <= nblines:\n",
    "        print(lines)\n",
    "        linecount += 1\n",
    "    \n",
    "with open(\"data\\\\covid-19-pandemic-worldwide-data.csv\", \"r\") as csv_file:\n",
    "    csv_reader = csv.reader(csv_file, delimiter = ';')\n",
    "    # create empty list for data\n",
    "    deaddata = []\n",
    "    for cnt,lines in enumerate(csv_reader):\n",
    "        # do not take header into account\n",
    "        if cnt > 0:\n",
    "            deaddata.append(float(lines[1]))\n",
    "\n",
    "deaddata2 = XX\n",
    "\n",
    "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.5), tight_layout=True)\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "ax1.plot(deaddata,label = \"cumulative dead people\")\n",
    "ax1.set_xlabel('day')\n",
    "ax1.set_ylabel('cumulative dead people')\n",
    "ax1.legend()\n",
    "\n",
    "ax2 = plt.subplot(1,2,2)\n",
    "ax2.plot(deaddata2,label = \"daily dead people\")\n",
    "ax2.set_xlabel('day')\n",
    "ax2.set_ylabel('daily dead people')\n",
    "ax2.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now fit a model to the data:\n",
    "- What kind of a model do you expect?\n",
    "- Try and fit these models\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "def myfunc(x, a, b):\n",
    "    return a*x + b\n",
    "\n",
    "\n",
    "# first wave between day XX and day XX\n",
    "d01 = XX\n",
    "df1 = XX\n",
    "dead_wave1 = deaddata2[d01:df1]\n",
    "ydata = dead_wave1\n",
    "xdata = range(d01,d01+len(ydata))\n",
    "\n",
    "pars, cov = curve_fit(f=myfunc, xdata = xdata, ydata = ydata, p0=[0, 0], bounds=(-np.inf, np.inf))\n",
    "[a, b] = pars\n",
    "print(f\"a = {a}, b = {b}\")\n",
    "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4.5), tight_layout=True)\n",
    "ax1 = plt.subplot(1,2,1)\n",
    "ax1.plot(deaddata2,label = \"daily dead people\")\n",
    "ax1.plot(xdata, XX, '-r',label = \"my  fit\")\n",
    "ax1.set_xlabel('day')\n",
    "ax1.set_ylabel('daily dead people')\n",
    "ax1.set_yscale('log')\n",
    "ax1.legend()\n",
    "\n",
    "\n",
    "# # second wave between day XX and XX\n",
    "\n"
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}