NDVI, NDWI, PSRI for komotini parcels.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from datetime import datetime\n",
    "from scipy.signal import savgol_filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NDWI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndwi = pd.read_csv(\"C:/Users/antonkout/Documents/Work/antonkoutpersolnal/komotini/index/featurespace_object_NDWI.csv\", header=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndwi.columns.values[1:] = [1,6,16,26,31,36,66,71,86,91,96,111,116,121,136,141,146,151,156,161,166,171,176,181,186,191,\n",
    "                          196,201,206,211,216,221,226,231,236,241,246,251,256,261,266,271,281,291,306,314,316,336,341,361]\n",
    "obj_ids = ndwi.iloc[:,0]\n",
    "del ndwi['id']\n",
    "ndwi = ndwi.replace(-9999999, -999999)\n",
    "ndwi = ndwi.replace(np.nan, -999999)\n",
    "new_days = [\"2018-01-10\", \"2018-01-20\", \"2018-01-30\", \"2018-02-10\", \"2018-02-20\", \"2018-02-28\",\"2018-03-10\", \"2018-03-20\", \n",
    "            \"2018-03-30\", \"2018-04-10\", \"2018-04-20\", \"2018-04-30\", \"2018-05-10\", \"2018-05-20\", \"2018-05-30\", \"2018-06-10\", \n",
    "            \"2018-06-20\", \"2018-06-30\", \"2018-07-10\", \"2018-07-20\", \"2018-07-30\", \"2018-08-10\", \"2018-08-20\", \"2018-08-30\", \n",
    "            \"2018-09-10\", \"2018-09-20\", \"2018-09-30\", \"2018-10-10\", \"2018-10-20\", \"2018-10-30\", \"2018-11-10\", \"2018-11-20\",\n",
    "            \"2018-11-30\",\"2018-12-10\",\"2018-12-20\",\"2018-12-30\"]\n",
    "ndwi = ndwi.T\n",
    "ndwi.index.name = 'doys' #Create decadal dates to have a fixed timestep\n",
    "ndwi.reset_index(inplace=True)\n",
    "nd = list()\n",
    "for i in range(len(new_days)):\n",
    "    nd.append(datetime.strptime(new_days[i], '%Y-%m-%d').timetuple().tm_yday)\n",
    "ndwi=ndwi[1:ndwi.shape[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Number of total outliers found: 15 ---\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\antonkout\\anaconda3\\envs\\e-shape\\lib\\site-packages\\numpy\\lib\\nanfunctions.py:1667: RuntimeWarning: Degrees of freedom <= 0 for slice.\n",
      "  keepdims=keepdims)\n",
      "C:\\Users\\antonkout\\anaconda3\\envs\\e-shape\\lib\\site-packages\\ipykernel_launcher.py:8: RuntimeWarning: Mean of empty slice\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "ndwinn = ndwi.drop(['doys'], axis=1)\n",
    "ndwinn.index = ndwi.doys\n",
    "ndwinn = ndwinn.replace(-999999, np.nan)\n",
    "ndwinn = ndwinn.replace(-9999999, np.nan)\n",
    "l=np.zeros(shape=(ndwinn.shape[0],2))\n",
    "for i in range(ndwinn.shape[0]):\n",
    "    l[i,0] = np.nanstd(ndwinn.iloc[i,:])\n",
    "    l[i,1] = np.nanmean(ndwinn.iloc[i,:])\n",
    "l = pd.DataFrame(l)\n",
    "l.columns = ['std_l','mean_l']\n",
    "\n",
    "#Find 99.7 of the data that are within +- 3std\n",
    "l['lower_limit'] = (l.mean_l-3*l.std_l)\n",
    "l['upper_limit'] = (l.mean_l+3*l.std_l)\n",
    "l = pd.DataFrame(l)\n",
    "l.index = ndwinn.index\n",
    "outliers = []\n",
    "for i in range(ndwinn.shape[0]):\n",
    "    for j in range(ndwinn.shape[1]):\n",
    "        if((ndwinn.iloc[i,j]<l.iloc[i,2]) or (ndwinn.iloc[i,j]>l.iloc[i,3])):\n",
    "            outliers.append(ndwinn.iloc[i,j])\n",
    "            ndwinn.iloc[i,j] = np.nan\n",
    "print('--- Number of total outliers found: %s ---'%len(outliers))\n",
    "ndwinn.insert (0, 'doys', ndwinn.index)\n",
    "ndwinn.index = range(0,ndwinn.shape[0])\n",
    "ndwinn = ndwinn.replace(np.nan,-999999)\n",
    "ndwi = ndwinn.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndwi = ndwi.replace(np.nan,-999999)\n",
    "#Weighted Average Interpolation, search +-5 days to find the reference days \n",
    "#from which the new dates will be created\n",
    "l = np.zeros(ndwi.shape[0])\n",
    "for i in range(ndwi.shape[0]):\n",
    "    for j in range(len(nd)):\n",
    "        x = ndwi.iloc[i,0]\n",
    "        y = nd[j]\n",
    "        if (abs(x-y) <= 5):\n",
    "            l[i] = y\n",
    "ndwi = pd.concat([pd.DataFrame(l),ndwi.reset_index(drop=True)], axis=1) #\n",
    "ndwi.columns.values[0] = 'ref'\n",
    "del new_days,nd\n",
    "colnames = ndwi.columns.values\n",
    "def interp_weight (df):\n",
    "    dfn = np.zeros(shape=(df.shape[0],df.shape[1]))    \n",
    "    #colnames = df.columns.values\n",
    "    i = 0\n",
    "    while (i <= (df.shape[0]-1)):\n",
    "        if (i == (df.shape[0]-1)):\n",
    "            dfn[i,:] = df[i,:] #When i becomes the last row then import the last row\n",
    "            i=i+1\n",
    "        elif (df[i,0] != df[i+1,0]):\n",
    "            dfn[i,:] = df[i,:] #If value of ref at i is different at value of ref at i+1 then import the i line \n",
    "            i=i+1        \n",
    "        else:\n",
    "            if (df[i,0] == df[i,1]): #If value of ref at i is equal to value of doys at the same line i then import i line\n",
    "                dfn[i,:] = df[i,:]\n",
    "                i=i+2\n",
    "            elif (df[i+1,0] == df[i+1,1]): #If value of ref at i+1 is equal to value of doys at the i+1 line import i+1 line \n",
    "                dfn[i,:] = df[i+1,:]\n",
    "                i=i+2\n",
    "            else:\n",
    "                diff = df[i+1,1] - df[i,1]\n",
    "                dist = abs(df[i,1] - df[i,0])\n",
    "                w = (diff-dist)/diff\n",
    "                df[i,:][df[i,:]==(-999999)] = df[i+1,:][df[i,:]==(-999999)]\n",
    "                df[i+1,:][df[i+1,]==(-999999)] = df[i,:][df[i+1,:]==(-999999)]\n",
    "                new_row = w*df[i,:] + (1-w)*df[i+1,] \n",
    "                dfn[i,:] = new_row\n",
    "                i = i + 2    \n",
    "    #dfn.columns = colnames\n",
    "    return pd.DataFrame(dfn)\n",
    "ndwi = np.array(ndwi)\n",
    "ndwi = interp_weight(ndwi)\n",
    "ndwi.columns = colnames\n",
    "ndwi = ndwi.replace(-999999.0000000001, np.nan)\n",
    "ndwi = ndwi.replace(-999999.0, np.nan)\n",
    "ndwi = ndwi.drop(['doys'], axis=1)\n",
    "del colnames\n",
    "d = list()\n",
    "for i in range(ndwi.shape[0]):\n",
    "    if (ndwi.iloc[i,0]==0):\n",
    "        d.append(i)      \n",
    "ndwi = ndwi.drop(d, axis=0)  \n",
    "del d\n",
    "ndwi = ndwi.drop(ndwi.index[16])\n",
    "ndwi.index = ndwi.ref\n",
    "ndwi = ndwi.drop(['ref'], axis=1)\n",
    "ndwin = ndwi.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "polynomial = np.zeros(shape=ndwin.shape)\n",
    "obj_id = np.array(ndwin.columns)\n",
    "#Apply Polynomial Interpolation second order, replace timeseries with less than 3 values with nan to remove them afterwards\n",
    "time_with_less = []\n",
    "for i in range(ndwin.shape[1]):\n",
    "    if (len(ndwin.iloc[:,i].unique())>3):\n",
    "        pd.DataFrame(polynomial).iloc[:,i] = pd.DataFrame(np.array(ndwin.iloc[:,i]).flatten()).interpolate(method ='polynomial', order=2).ffill().bfill()\n",
    "    else:\n",
    "        time_with_less.append(ndwin.columns[i])\n",
    "        polynomial[:,i] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 0 timeseries with less than 3 values\n",
      "--- Number of total outliers found: 5 ---\n"
     ]
    }
   ],
   "source": [
    "print('Found %s timeseries with less than 3 values' %len(time_with_less))\n",
    "polynomial = pd.DataFrame(polynomial)\n",
    "polynomial.index = ndwin.index\n",
    "###Apply again Outlier Detection with Normal Distribution\n",
    "l = np.zeros(shape=(polynomial.shape[0],2))\n",
    "for i in range(polynomial.shape[0]):\n",
    "    l[i,0] = np.nanstd(polynomial.iloc[i,:])\n",
    "    l[i,1] = np.nanmean(polynomial.iloc[i,:])\n",
    "l = pd.DataFrame(l)\n",
    "l.columns = ['std_l','mean_l']\n",
    "\n",
    "#Find 99.7 of the data that are within +- 3std\n",
    "l['lower_limit'] = (l.mean_l-3*l.std_l)\n",
    "l['upper_limit'] = (l.mean_l+3*l.std_l)\n",
    "l = pd.DataFrame(l)\n",
    "l.index = polynomial.index\n",
    "outliers = []\n",
    "for i in range(polynomial.shape[0]):\n",
    "    for j in range(polynomial.shape[1]):\n",
    "        if((polynomial.iloc[i,j]<l.iloc[i,2]) or (polynomial.iloc[i,j]>l.iloc[i,3])):\n",
    "            outliers.append(polynomial.iloc[i,j])\n",
    "            polynomial.iloc[i,j] = np.nan \n",
    "print('--- Number of total outliers found: %s ---'%len(outliers))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---Drop 0 from 45 timeseries with less than 3 values---\n"
     ]
    }
   ],
   "source": [
    "new_doys = list(range(0,365))\n",
    "old_doys = np.array(polynomial.index)\n",
    "new_doys = np.array(new_doys)\n",
    "ndwin = np.zeros(shape=(len(new_doys),polynomial.shape[1]))\n",
    "ndwin[ndwin == 0] = 'nan'\n",
    "\n",
    "for i in range(ndwin.shape[0]):\n",
    "    if (i in old_doys):\n",
    "        ndwin[i,:] = polynomial.loc[polynomial.index==i,:]\n",
    "#ndwin[:,0] = range(0,365)\n",
    "ndwin = pd.DataFrame(ndwin)\n",
    "#ndwin.columns = obj_id\n",
    "print('---Drop %s from %s timeseries with less than 3 values---' %(len(time_with_less) ,ndwin.shape[1]))\n",
    "\n",
    "#ids = np.setdiff1d(np.union1d(obj_id, time_with_less), np.intersect1d(obj_id, time_with_less))\n",
    "ndwin = ndwin.dropna(axis='columns', how='all')\n",
    "\n",
    "linear = np.zeros(shape=ndwin.shape)\n",
    "\n",
    "for i in range(ndwin.shape[1]):\n",
    "    linear[:,i] = ndwin.iloc[:,i].interpolate(method ='linear', axis=0).ffill().bfill()\n",
    "\n",
    "linear = pd.DataFrame(linear)\n",
    "linear.columns = ndwin.columns\n",
    "savitzky = np.zeros(shape=linear.shape)\n",
    "for i in range(savitzky.shape[1]):\n",
    "    savitzky[:,i] = savgol_filter(np.array(linear.iloc[0:,i]).reshape(linear.shape[0]),mode='nearest', \n",
    "                                  window_length = 51, polyorder = 4)\n",
    "savitzky = pd.DataFrame(savitzky)\n",
    "savitzky.index = range(0,savitzky.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndwi = savitzky.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PSRI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "psri = pd.read_csv(\"C:/Users/antonkout/Documents/Work/antonkoutpersolnal/komotini/index/featurespace_object_PSRI.csv\", header=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "psri.columns.values[1:] = [1,6,16,26,31,36,66,71,86,91,96,111,116,121,136,141,146,151,156,161,166,171,176,181,186,191,\n",
    "                          196,201,206,211,216,221,226,231,236,241,246,251,256,261,266,271,281,291,306,314,316,336,341,361]\n",
    "obj_id = psri.iloc[:,0]\n",
    "del psri['id']\n",
    "psri = psri.replace(-9999999, -999999)\n",
    "psri = psri.replace(np.nan, -999999)\n",
    "new_days = [\"2018-01-10\", \"2018-01-20\", \"2018-01-30\", \"2018-02-10\", \"2018-02-20\", \"2018-02-28\",\"2018-03-10\", \"2018-03-20\", \n",
    "            \"2018-03-30\", \"2018-04-10\", \"2018-04-20\", \"2018-04-30\", \"2018-05-10\", \"2018-05-20\", \"2018-05-30\", \"2018-06-10\", \n",
    "            \"2018-06-20\", \"2018-06-30\", \"2018-07-10\", \"2018-07-20\", \"2018-07-30\", \"2018-08-10\", \"2018-08-20\", \"2018-08-30\", \n",
    "            \"2018-09-10\", \"2018-09-20\", \"2018-09-30\", \"2018-10-10\", \"2018-10-20\", \"2018-10-30\", \"2018-11-10\", \"2018-11-20\",\n",
    "            \"2018-11-30\",\"2018-12-10\",\"2018-12-20\",\"2018-12-30\"]\n",
    "psri = psri.T\n",
    "psri.index.name = 'doys' #Create decadal dates to have a fixed timestep\n",
    "psri.reset_index(inplace=True)\n",
    "nd = list()\n",
    "for i in range(len(new_days)):\n",
    "    nd.append(datetime.strptime(new_days[i], '%Y-%m-%d').timetuple().tm_yday)\n",
    "psri=psri[1:psri.shape[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Number of total outliers found: 4 ---\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\antonkout\\anaconda3\\envs\\e-shape\\lib\\site-packages\\ipykernel_launcher.py:8: RuntimeWarning: Mean of empty slice\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "psrinn = psri.drop(['doys'], axis=1)\n",
    "psrinn.index = psri.doys\n",
    "psrinn = psrinn.replace(-999999, np.nan)\n",
    "psrinn = psrinn.replace(-9999999, np.nan)\n",
    "l=np.zeros(shape=(psrinn.shape[0],2))\n",
    "for i in range(psrinn.shape[0]):\n",
    "    l[i,0] = np.nanstd(psrinn.iloc[i,:])\n",
    "    l[i,1] = np.nanmean(psrinn.iloc[i,:])\n",
    "l = pd.DataFrame(l)\n",
    "l.columns = ['std_l','mean_l']\n",
    "\n",
    "#Find 99.7 of the data that are within +- 3std\n",
    "l['lower_limit'] = (l.mean_l-3*l.std_l)\n",
    "l['upper_limit'] = (l.mean_l+3*l.std_l)\n",
    "l = pd.DataFrame(l)\n",
    "l.index = psrinn.index\n",
    "outliers = []\n",
    "for i in range(psrinn.shape[0]):\n",
    "    for j in range(psrinn.shape[1]):\n",
    "        if((psrinn.iloc[i,j]<l.iloc[i,2]) or (psrinn.iloc[i,j]>l.iloc[i,3])):\n",
    "            outliers.append(psrinn.iloc[i,j])\n",
    "            psrinn.iloc[i,j] = np.nan\n",
    "print('--- Number of total outliers found: %s ---'%len(outliers))\n",
    "psrinn.insert (0, 'doys', psrinn.index)\n",
    "psrinn.index = range(0,psrinn.shape[0])\n",
    "psrinn = psrinn.replace(np.nan,-999999)\n",
    "psri = psrinn.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "psri = psri.replace(np.nan,-999999)\n",
    "#Weighted Average Interpolation, search +-5 days to find the reference days \n",
    "#from which the new dates will be created\n",
    "l = np.zeros(psri.shape[0])\n",
    "for i in range(psri.shape[0]):\n",
    "    for j in range(len(nd)):\n",
    "        x = psri.iloc[i,0]\n",
    "        y = nd[j]\n",
    "        if (abs(x-y) <= 5):\n",
    "            l[i] = y\n",
    "psri = pd.concat([pd.DataFrame(l),psri.reset_index(drop=True)], axis=1) #\n",
    "psri.columns.values[0] = 'ref'\n",
    "del new_days,nd\n",
    "colnames = psri.columns.values\n",
    "def interp_weight (df):\n",
    "    dfn = np.zeros(shape=(df.shape[0],df.shape[1]))    \n",
    "    #colnames = df.columns.values\n",
    "    i = 0\n",
    "    while (i <= (df.shape[0]-1)):\n",
    "        if (i == (df.shape[0]-1)):\n",
    "            dfn[i,:] = df[i,:] #When i becomes the last row then import the last row\n",
    "            i=i+1\n",
    "        elif (df[i,0] != df[i+1,0]):\n",
    "            dfn[i,:] = df[i,:] #If value of ref at i is different at value of ref at i+1 then import the i line \n",
    "            i=i+1        \n",
    "        else:\n",
    "            if (df[i,0] == df[i,1]): #If value of ref at i is equal to value of doys at the same line i then import i line\n",
    "                dfn[i,:] = df[i,:]\n",
    "                i=i+2\n",
    "            elif (df[i+1,0] == df[i+1,1]): #If value of ref at i+1 is equal to value of doys at the i+1 line import i+1 line \n",
    "                dfn[i,:] = df[i+1,:]\n",
    "                i=i+2\n",
    "            else:\n",
    "                diff = df[i+1,1] - df[i,1]\n",
    "                dist = abs(df[i,1] - df[i,0])\n",
    "                w = (diff-dist)/diff\n",
    "                df[i,:][df[i,:]==(-999999)] = df[i+1,:][df[i,:]==(-999999)]\n",
    "                df[i+1,:][df[i+1,]==(-999999)] = df[i,:][df[i+1,:]==(-999999)]\n",
    "                new_row = w*df[i,:] + (1-w)*df[i+1,] \n",
    "                dfn[i,:] = new_row\n",
    "                i = i + 2    \n",
    "    #dfn.columns = colnames\n",
    "    return pd.DataFrame(dfn)\n",
    "psri = np.array(psri)\n",
    "psri = interp_weight(psri)\n",
    "psri.columns = colnames\n",
    "psri = psri.replace(-999999.0000000001, np.nan)\n",
    "psri = psri.replace(-999999.0, np.nan)\n",
    "psri = psri.drop(['doys'], axis=1)\n",
    "del colnames\n",
    "d = list()\n",
    "for i in range(psri.shape[0]):\n",
    "    if (psri.iloc[i,0]==0):\n",
    "        d.append(i)      \n",
    "psri = psri.drop(d, axis=0)  \n",
    "del d\n",
    "psri = psri.drop(psri.index[16])\n",
    "psri.index = psri.ref\n",
    "psri = psri.drop(['ref'], axis=1)\n",
    "psrin = psri.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "polynomial = np.zeros(shape=psrin.shape)\n",
    "obj_id = np.array(psrin.columns)\n",
    "#Apply Polynomial Interpolation second order, replace timeseries with less than 3 values with nan to remove them afterwards\n",
    "time_with_less = []\n",
    "for i in range(psrin.shape[1]):\n",
    "    if (len(psrin.iloc[:,i].unique())>3):\n",
    "        pd.DataFrame(polynomial).iloc[:,i] = pd.DataFrame(np.array(psrin.iloc[:,i]).flatten()).interpolate(method ='polynomial', order=2).ffill().bfill()\n",
    "    else:\n",
    "        time_with_less.append(psrin.columns[i])\n",
    "        polynomial[:,i] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 0 timeseries with less than 3 values\n",
      "--- Number of total outliers found: 3 ---\n"
     ]
    }
   ],
   "source": [
    "print('Found %s timeseries with less than 3 values' %len(time_with_less))\n",
    "polynomial = pd.DataFrame(polynomial)\n",
    "polynomial.index = psrin.index\n",
    "###Apply again Outlier Detection with Normal Distribution\n",
    "l = np.zeros(shape=(polynomial.shape[0],2))\n",
    "for i in range(polynomial.shape[0]):\n",
    "    l[i,0] = np.nanstd(polynomial.iloc[i,:])\n",
    "    l[i,1] = np.nanmean(polynomial.iloc[i,:])\n",
    "l = pd.DataFrame(l)\n",
    "l.columns = ['std_l','mean_l']\n",
    "\n",
    "#Find 99.7 of the data that are within +- 3std\n",
    "l['lower_limit'] = (l.mean_l-3*l.std_l)\n",
    "l['upper_limit'] = (l.mean_l+3*l.std_l)\n",
    "l = pd.DataFrame(l)\n",
    "l.index = polynomial.index\n",
    "outliers = []\n",
    "for i in range(polynomial.shape[0]):\n",
    "    for j in range(polynomial.shape[1]):\n",
    "        if((polynomial.iloc[i,j]<l.iloc[i,2]) or (polynomial.iloc[i,j]>l.iloc[i,3])):\n",
    "            outliers.append(polynomial.iloc[i,j])\n",
    "            polynomial.iloc[i,j] = np.nan \n",
    "print('--- Number of total outliers found: %s ---'%len(outliers))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---Drop 0 from 45 timeseries with less than 3 values---\n"
     ]
    }
   ],
   "source": [
    "new_doys = list(range(0,365))\n",
    "old_doys = np.array(polynomial.index)\n",
    "new_doys = np.array(new_doys)\n",
    "psrin = np.zeros(shape=(len(new_doys),polynomial.shape[1]))\n",
    "psrin[psrin == 0] = 'nan'\n",
    "\n",
    "for i in range(psrin.shape[0]):\n",
    "    if (i in old_doys):\n",
    "        psrin[i,:] = polynomial.loc[polynomial.index==i,:]\n",
    "#psrin[:,0] = range(0,365)\n",
    "psrin = pd.DataFrame(psrin)\n",
    "#psrin.columns = obj_id\n",
    "print('---Drop %s from %s timeseries with less than 3 values---' %(len(time_with_less) ,psrin.shape[1]))\n",
    "\n",
    "#ids = np.setdiff1d(np.union1d(obj_id, time_with_less), np.intersect1d(obj_id, time_with_less))\n",
    "psrin = psrin.dropna(axis='columns', how='all')\n",
    "\n",
    "linear = np.zeros(shape=psrin.shape)\n",
    "\n",
    "for i in range(psrin.shape[1]):\n",
    "    linear[:,i] = psrin.iloc[:,i].interpolate(method ='linear', axis=0).ffill().bfill()\n",
    "\n",
    "linear = pd.DataFrame(linear)\n",
    "linear.columns = psrin.columns\n",
    "savitzky = np.zeros(shape=linear.shape)\n",
    "for i in range(savitzky.shape[1]):\n",
    "    savitzky[:,i] = savgol_filter(np.array(linear.iloc[0:,i]).reshape(linear.shape[0]),mode='nearest', \n",
    "                                  window_length = 51, polyorder = 4)\n",
    "savitzky = pd.DataFrame(savitzky)\n",
    "savitzky.index = range(0,savitzky.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "psri = savitzky.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NDVI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi = pd.read_csv(\"C:/Users/antonkout/Documents/Work/antonkoutpersolnal/komotini/index/featurespace_object_NDVI.csv\", header=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi.columns.values[1:] = [1,6,16,26,31,36,66,71,86,91,96,111,116,121,136,141,146,151,156,161,166,171,176,181,186,191,\n",
    "                          196,201,206,211,216,221,226,231,236,241,246,251,256,261,266,271,281,291,306,314,316,336,341,361]\n",
    "obj_id = ndvi.iloc[:,0]\n",
    "del ndvi['id']\n",
    "ndvi = ndvi.replace(-9999999, -999999)\n",
    "ndvi = ndvi.replace(np.nan, -999999)\n",
    "new_days = [\"2018-01-10\", \"2018-01-20\", \"2018-01-30\", \"2018-02-10\", \"2018-02-20\", \"2018-02-28\",\"2018-03-10\", \"2018-03-20\", \n",
    "            \"2018-03-30\", \"2018-04-10\", \"2018-04-20\", \"2018-04-30\", \"2018-05-10\", \"2018-05-20\", \"2018-05-30\", \"2018-06-10\", \n",
    "            \"2018-06-20\", \"2018-06-30\", \"2018-07-10\", \"2018-07-20\", \"2018-07-30\", \"2018-08-10\", \"2018-08-20\", \"2018-08-30\", \n",
    "            \"2018-09-10\", \"2018-09-20\", \"2018-09-30\", \"2018-10-10\", \"2018-10-20\", \"2018-10-30\", \"2018-11-10\", \"2018-11-20\",\n",
    "            \"2018-11-30\",\"2018-12-10\",\"2018-12-20\",\"2018-12-30\"]\n",
    "ndvi = ndvi.T\n",
    "ndvi.index.name = 'doys' #Create decadal dates to have a fixed timestep\n",
    "ndvi.reset_index(inplace=True)\n",
    "nd = list()\n",
    "for i in range(len(new_days)):\n",
    "    nd.append(datetime.strptime(new_days[i], '%Y-%m-%d').timetuple().tm_yday)\n",
    "ndvi=ndvi[1:ndvi.shape[0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Number of total outliers found: 22 ---\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\antonkout\\anaconda3\\envs\\e-shape\\lib\\site-packages\\ipykernel_launcher.py:8: RuntimeWarning: Mean of empty slice\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "ndvinn = ndvi.drop(['doys'], axis=1)\n",
    "ndvinn.index = ndvi.doys\n",
    "ndvinn = ndvinn.replace(-999999, np.nan)\n",
    "ndvinn = ndvinn.replace(-9999999, np.nan)\n",
    "l=np.zeros(shape=(ndvinn.shape[0],2))\n",
    "for i in range(ndvinn.shape[0]):\n",
    "    l[i,0] = np.nanstd(ndvinn.iloc[i,:])\n",
    "    l[i,1] = np.nanmean(ndvinn.iloc[i,:])\n",
    "l = pd.DataFrame(l)\n",
    "l.columns = ['std_l','mean_l']\n",
    "\n",
    "#Find 99.7 of the data that are within +- 3std\n",
    "l['lower_limit'] = (l.mean_l-3*l.std_l)\n",
    "l['upper_limit'] = (l.mean_l+3*l.std_l)\n",
    "l = pd.DataFrame(l)\n",
    "l.index = ndvinn.index\n",
    "outliers = []\n",
    "for i in range(ndvinn.shape[0]):\n",
    "    for j in range(ndvinn.shape[1]):\n",
    "        if((ndvinn.iloc[i,j]<l.iloc[i,2]) or (ndvinn.iloc[i,j]>l.iloc[i,3])):\n",
    "            outliers.append(ndvinn.iloc[i,j])\n",
    "            ndvinn.iloc[i,j] = np.nan\n",
    "print('--- Number of total outliers found: %s ---'%len(outliers))\n",
    "ndvinn.insert (0, 'doys', ndvinn.index)\n",
    "ndvinn.index = range(0,ndvinn.shape[0])\n",
    "ndvinn = ndvinn.replace(np.nan,-999999)\n",
    "ndvi = ndvinn.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi = ndvi.replace(np.nan,-999999)\n",
    "#Weighted Average Interpolation, search +-5 days to find the reference days \n",
    "#from which the new dates will be created\n",
    "l = np.zeros(ndvi.shape[0])\n",
    "for i in range(ndvi.shape[0]):\n",
    "    for j in range(len(nd)):\n",
    "        x = ndvi.iloc[i,0]\n",
    "        y = nd[j]\n",
    "        if (abs(x-y) <= 5):\n",
    "            l[i] = y\n",
    "ndvi = pd.concat([pd.DataFrame(l),ndvi.reset_index(drop=True)], axis=1) #\n",
    "ndvi.columns.values[0] = 'ref'\n",
    "del new_days,nd\n",
    "colnames = ndvi.columns.values\n",
    "def interp_weight (df):\n",
    "    dfn = np.zeros(shape=(df.shape[0],df.shape[1]))    \n",
    "    #colnames = df.columns.values\n",
    "    i = 0\n",
    "    while (i <= (df.shape[0]-1)):\n",
    "        if (i == (df.shape[0]-1)):\n",
    "            dfn[i,:] = df[i,:] #When i becomes the last row then import the last row\n",
    "            i=i+1\n",
    "        elif (df[i,0] != df[i+1,0]):\n",
    "            dfn[i,:] = df[i,:] #If value of ref at i is different at value of ref at i+1 then import the i line \n",
    "            i=i+1        \n",
    "        else:\n",
    "            if (df[i,0] == df[i,1]): #If value of ref at i is equal to value of doys at the same line i then import i line\n",
    "                dfn[i,:] = df[i,:]\n",
    "                i=i+2\n",
    "            elif (df[i+1,0] == df[i+1,1]): #If value of ref at i+1 is equal to value of doys at the i+1 line import i+1 line \n",
    "                dfn[i,:] = df[i+1,:]\n",
    "                i=i+2\n",
    "            else:\n",
    "                diff = df[i+1,1] - df[i,1]\n",
    "                dist = abs(df[i,1] - df[i,0])\n",
    "                w = (diff-dist)/diff\n",
    "                df[i,:][df[i,:]==(-999999)] = df[i+1,:][df[i,:]==(-999999)]\n",
    "                df[i+1,:][df[i+1,]==(-999999)] = df[i,:][df[i+1,:]==(-999999)]\n",
    "                new_row = w*df[i,:] + (1-w)*df[i+1,] \n",
    "                dfn[i,:] = new_row\n",
    "                i = i + 2    \n",
    "    #dfn.columns = colnames\n",
    "    return pd.DataFrame(dfn)\n",
    "ndvi = np.array(ndvi)\n",
    "ndvi = interp_weight(ndvi)\n",
    "ndvi.columns = colnames\n",
    "ndvi = ndvi.replace(-999999.0000000001, np.nan)\n",
    "ndvi = ndvi.replace(-999999.0, np.nan)\n",
    "ndvi = ndvi.drop(['doys'], axis=1)\n",
    "del colnames\n",
    "d = list()\n",
    "for i in range(ndvi.shape[0]):\n",
    "    if (ndvi.iloc[i,0]==0):\n",
    "        d.append(i)      \n",
    "ndvi = ndvi.drop(d, axis=0)  \n",
    "del d\n",
    "ndvi = ndvi.drop(ndvi.index[16])\n",
    "ndvi.index = ndvi.ref\n",
    "ndvi = ndvi.drop(['ref'], axis=1)\n",
    "ndvin = ndvi.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "polynomial = np.zeros(shape=ndvin.shape)\n",
    "obj_id = np.array(ndvin.columns)\n",
    "#Apply Polynomial Interpolation second order, replace timeseries with less than 3 values with nan to remove them afterwards\n",
    "time_with_less = []\n",
    "for i in range(ndvin.shape[1]):\n",
    "    if (len(ndvin.iloc[:,i].unique())>3):\n",
    "        pd.DataFrame(polynomial).iloc[:,i] = pd.DataFrame(np.array(ndvin.iloc[:,i]).flatten()).interpolate(method ='polynomial', order=2).ffill().bfill()\n",
    "    else:\n",
    "        time_with_less.append(ndvin.columns[i])\n",
    "        polynomial[:,i] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 0 timeseries with less than 3 values\n",
      "--- Number of total outliers found: 14 ---\n"
     ]
    }
   ],
   "source": [
    "print('Found %s timeseries with less than 3 values' %len(time_with_less))\n",
    "polynomial = pd.DataFrame(polynomial)\n",
    "polynomial.index = ndvin.index\n",
    "###Apply again Outlier Detection with Normal Distribution\n",
    "l = np.zeros(shape=(polynomial.shape[0],2))\n",
    "for i in range(polynomial.shape[0]):\n",
    "    l[i,0] = np.nanstd(polynomial.iloc[i,:])\n",
    "    l[i,1] = np.nanmean(polynomial.iloc[i,:])\n",
    "l = pd.DataFrame(l)\n",
    "l.columns = ['std_l','mean_l']\n",
    "\n",
    "#Find 99.7 of the data that are within +- 3std\n",
    "l['lower_limit'] = (l.mean_l-3*l.std_l)\n",
    "l['upper_limit'] = (l.mean_l+3*l.std_l)\n",
    "l = pd.DataFrame(l)\n",
    "l.index = polynomial.index\n",
    "outliers = []\n",
    "for i in range(polynomial.shape[0]):\n",
    "    for j in range(polynomial.shape[1]):\n",
    "        if((polynomial.iloc[i,j]<l.iloc[i,2]) or (polynomial.iloc[i,j]>l.iloc[i,3])):\n",
    "            outliers.append(polynomial.iloc[i,j])\n",
    "            polynomial.iloc[i,j] = np.nan \n",
    "print('--- Number of total outliers found: %s ---'%len(outliers))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---Drop 0 from 45 timeseries with less than 3 values---\n"
     ]
    }
   ],
   "source": [
    "new_doys = list(range(0,365))\n",
    "old_doys = np.array(polynomial.index)\n",
    "new_doys = np.array(new_doys)\n",
    "ndvin = np.zeros(shape=(len(new_doys),polynomial.shape[1]))\n",
    "ndvin[ndvin == 0] = 'nan'\n",
    "\n",
    "for i in range(ndvin.shape[0]):\n",
    "    if (i in old_doys):\n",
    "        ndvin[i,:] = polynomial.loc[polynomial.index==i,:]\n",
    "#ndvin[:,0] = range(0,365)\n",
    "ndvin = pd.DataFrame(ndvin)\n",
    "#ndvin.columns = obj_id\n",
    "print('---Drop %s from %s timeseries with less than 3 values---' %(len(time_with_less) ,ndvin.shape[1]))\n",
    "\n",
    "#ids = np.setdiff1d(np.union1d(obj_id, time_with_less), np.intersect1d(obj_id, time_with_less))\n",
    "ndvin = ndvin.dropna(axis='columns', how='all')\n",
    "\n",
    "linear = np.zeros(shape=ndvin.shape)\n",
    "\n",
    "for i in range(ndvin.shape[1]):\n",
    "    linear[:,i] = ndvin.iloc[:,i].interpolate(method ='linear', axis=0).ffill().bfill()\n",
    "\n",
    "linear = pd.DataFrame(linear)\n",
    "linear.columns = ndvin.columns\n",
    "savitzky = np.zeros(shape=linear.shape)\n",
    "for i in range(savitzky.shape[1]):\n",
    "    savitzky[:,i] = savgol_filter(np.array(linear.iloc[0:,i]).reshape(linear.shape[0]),mode='nearest', \n",
    "                                  window_length = 51, polyorder = 4)\n",
    "savitzky = pd.DataFrame(savitzky)\n",
    "savitzky.index = range(0,savitzky.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi = savitzky.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndvi.columns = np.array(obj_ids)\n",
    "ndwi.columns = np.array(obj_ids)\n",
    "psri.columns = np.array(obj_ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "validation = pd.read_csv(\"C:/Users/antonkout/Documents/Work/antonkoutpersolnal/komotinivalidation.csv\", header=0)\n",
    "validation['Length'] = validation['Harvest DOY'] - validation['Seeding DOY']\n",
    "validation.columns = ['parcelid','area','seeding','emergence','harvest','length']\n",
    "validation.index = np.array(validation['parcelid'])\n",
    "del validation['parcelid']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "ndwi.index = np.array(range(0,365))\n",
    "psri.index = np.array(range(0,365))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArQAAAIXCAYAAABtts2aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd1gVZ9rH8e8cOhwB6SBVwYK9xtiNBUvEErusstk1xqibWGOyMclmWzQR0zTGWBKzMVhib6hJEEtUMPZGE2wUY8FOPe8fxPOG2A7IYc6B+3NdcyFTnvnNgHrz8Mwzik6n0yGEEEIIIYSZ0qgdQAghhBBCiKchBa0QQgghhDBrUtAKIYQQQgizJgWtEEIIIYQwa1LQCiGEEEIIsyYFrRBCCCGEMGtS0AohhBBCCLNmqXYAU2BjY4O7u3uFnS83NxcbG5sKO19lIPesbMrzvt27cgUAW1fXcmnPVMn3WtnIfSs9uWdPkJVV/NHTs8RquW+lV1nu2eXLl8nNzX3oNiloAXd3dy5cuFBh54uJiSEsLKzCzlcZyD0rm/K8b5vCwwHovX59ubRnquR7rWzkvpWe3LMnqF+/+OOJEyVWy30rvcpyz3x9fR+5TYYcCCGEEEIIsyYFrRBCCCGEMGtS0AohhBBCCLMmY2iFEEIIUWUVFRWh0+nUjmF0hYWFakd4LEVR0GjK3s8qBa0QwiCV/WEwIYSJ+cPDYOWtqKiI9PR07t27Z9TzmAJ3d3cSExPVjvFEtra2BAQElKmwlYJWCCGEEFVOdnY2Go2GkJAQFEVRO45R3bhxA0dHR7VjPJZOp+PixYtkZ2fj5eVV6uOloBVCGOTS7t0A+LRrp3ISIUSVEBNT/NEI003pdDquX79OYGAglpaVvxTSaDRYWFioHeOJPD09SUtLw9PTs9Q/ZFT+r6IQolwcmjULkIJWCFFBJk0q/miEoQc6nQ6dToeVlVW5ty3KzsrKSv+1KW1BK7McCCGEEKJKqQoPgZmzsnx9pKAVQgghhBBmTQpaIYQQQghh1qSgFUIIIYQwcZ06dUJRFLZt21ZifUREBJGRkQAEBgaiKAqKouDg4EC9evUYN24caWlp+v3nzJmDv79/iTZ27NiBoihERUWVWN+nTx9925GRkURERJT7dZUXKWiFEAbR+vmh9fNTO4YQoqqoVat4EXq2tra89dZbj91n9uzZZGRkcOLECaKiokhJSaFdu3YkJCQA0KFDB86fP1+iyN21axe+vr7ExcXp1xUVFbFnzx46dOhglGspb1LQCiEM0nHuXDrOnat2DCFEVbF+ffEi9P70pz9x8uRJ1q5d+8h9HB0d8fLyIjAwkJ49e7JlyxZatGjBuHHjAGjSpAmOjo4lite4uDimTJnC7t279Q9kHTt2jGvXrtGxY0fjXlQ5kYJWCCGEEMIMeHh48OqrrzJjxgyKiooMOkZRFEaPHs2BAwfIzMzEwsKCtm3b6gvavLw89u/fz7Bhw9BqtZz4bZq0uLg4fHx8qGUmveQyD60QwiCJy5YBUHv4cJWTCCGqhPu/EfqtZ9HY/vp1POlX7lTIuQACXO1ZOKplqY+bOnUq8+bNIzo6muEG/ntcu3ZtANLT0/Hy8qJDhw4sXrwYgPj4ePz9/fHw8KBdu3bExcXRoEED4uLizKZ3FqSHVghhoKToaJKio9WOIYSoKubNK15ECc7OzkyZMoV33nmHgoICg465P4zg/ssKOnbsSFJSEpmZmcTFxdG+fXsAfUELxeNqzWX8LEgPrRBCCCFEmXpL1fLqq6/y8ccf89VXXxm0f2JiIlA8CwJAixYtsLe3Jy4ujri4OH1Pb/v27Xnvvfc4c+YMWVlZZlXQSg+tEEIIs5RfmM/BrIMcvXyUa/euqR1HiAqj1Wp54403eO+998jNzX3svjqdjoULF9K6dWs8PDyA4lfMPvvss/z000/s3btX30MbGhpKbm4uixcvxt3dndDQUKNfS3mRHlohhBBmJb8on7mH5hJ9Jprb+bcBsFAs6BHUg7GNxxLgGKByQiGMb+zYsURFRbFx40aGDBmiX3/jxg0yMzO5d+8ep06d4uOPP+bgwYPExsaWOL5Dhw588MEHODs763tuFUWhTZs2zJs3j+7du1fg1Tw96aEVQghhNnJyc/hLzF9YdHwR3g7ejG8ynonNJ9LCswWbUjcxaMMgNqRsUDumEEZ3f07ae/fulVg/efJkvL29CQ0NZeLEidSqVYtdu3bRrFmzEvt17NiRW7du6Xtn72vfvj23bt0yq+EGID20QggDebVpo3YEUcUVFBUwdedUDmUfIqJeBJNaTMJKYwXAiw1e5EDGAd7Y9QZv7n6T7DvZ/KXhX1ROLJ6KmfUQGtsfe1gBxowZw5gxY/Sf//5lCb+Xk5PzwLqOHTvqHxb7vWnTpjFt2rQH1hs6XlctUtAKIQzSfPp0tSOIKu7TQ5/yc8bPvBDyAtNaTtM/sX1fK+9WLO+znJe2v8RHv3xEoa6Qlxq9pFJa8dTmzFE7gTAjMuRACCGEyUu8lshXJ76igWsD3nzmzQeK2fvc7NxY1H0RdV3q8umhT/n88OcVnFQIoQYpaIUQBjn4/vscfP99tWOIKkin0zHzwEx0Oh1vPfsW1hbWj92/um11FnZfSD2Xesw7Mo+5h+c+9FerwsRNnFi8CGEAGXIghDBI5t69akcQVVTs+VgOZB7ghZAXqO9a36BjnGyc+LL7l4zZPob5R+aTW5jLa81eQ6NIP47Z2LZN7QTCjMjfbCGEECZLp9Ox8NhCrDXWjGtSulegOtk4saD7Ahq5NWLJ8SWM/2E8N/JuGCmpEEJNUtAKIYQwWQlZCRz99Sj9gvvhbu9e6uMdrR1ZFLaI8Frh7Lq4i6Ebh5J4LdEISYUQapKCVgghhMladGwRGkVDZP3IMrdha2nLv9r+izefeZOMWxkM3zSc705/J+NqhahEpKAVQghhks7fOM+eS3vo4t8FP0e/p2pLURSG1R3Gkh5LcLNz4z/7/8OfY/5M0rWkckorhFCTPBQmhDBIyNChakcQVczalLUADAwZWG5tNvFowqo+q4g6GMWqxFW8sP4FegT2ICI0goZuDR85HZhQwSuvqJ1AmBEpaIUQBqk9fLjaEUQVUlhUyLrkdXg5ePGM9zPl2rbWWsvbz77NgJABfHboM7akbWFL2hZ8tb609GpJM89mNPdojm81Xylw1TSudA8BiqpNhhwIIYQwOXsv7SXrThb9gvthobEwyjkauDVgfrf5rOyzkoh6EejQsSZ5DTP2zKDXml50XN6RKTunEJMWQ0FRgVEyCGGoTp06oSgKiqKg1Wpp1aoVMTExAGRkZPDnP/+ZGjVqYGtri7+/P4MGDSI7OxuAXbt26Y/VaDT4+PgwatQoLl++rG8/LS0NRVFITk5W5fqelvTQCiEMsvO33pKOc+eqnERUBVvTtgIQXivc6Oeq61KXuq3q8nqr18m8ncnh7MP8kv0LRy8fZVvaNmLSYgh2DmZay2k86/Os0fOI34T/9rVfv17dHCbktdde4/XXX+fmzZvMnj2bvn37cvz4cUaOHImNjQ0rV67Ex8eH9PR01q9fz+3bt0scf+HCBTQaDWlpaYwdO5axY8eyatUqla6mfElBK4QwyK3z59WOIKqI/MJ8fjr/E/Vc6uFX7ekeBistLwcvegT1oEdQDwAu37lM9Jlovj7xNS9tf4lOfp2Y0mIKAY4BFZqrSkpJUTuByXFwcMDLywsvLy8+++wzvvnmG1atWsXPP/9MQkICzZs3ByAwMJCOHTs+cLynpyeWlpZ4e3vzl7/8hbfffruiL8FoZMiBEEIIk3Ig8wA3827SLaCb2lFwt3dnQtMJrO+3nh6BPYg9H0u/df2ISojidv7tJzcghJFYWlpiZWWFg4MDDg4OrFu3joICw4bGXL58mTVr1mBt/fjXSJsT6aEVQghhUranbwega0BXlZP8Px+tDx90/IChdYcy88BMlpxYwqbUTUxtOZWwwDB5eKwyWDYUrp2tuPNVD4Lh0WU6ND8/n9mzZ3Pz5k3at2/Pl19+yZgxY5gzZw6tWrWiS5cujBo1iho1apQ4ztnZGZ1Ox507dwCYOXPmU1+GqZAeWiGEECajsKiQH8/9SLBzMEFOQWrHeUBzz+Z81/s7ZrSewb3Ce0yNm8ro7aNJzUlVO5qoAmbNmoVWq8Xe3p4PPviAzz//nCZNmjBs2DAuXbrEt99+S8uWLfnqq68IDQ3lyJEjJY4/ePAgCQkJvP/++3Ts2JGJEyeqdCXlT3pohRBCmIyTV05yLfcaL9R+Qe0oj2ShsWBwncF0DejKx798zOqk1byw/gVGhY7ipUYvYW9lr3ZEURZl7C2tSKNHj2bixIlotVq8vLxKbNNqtYSHhxMeHs4///lPmjZtyuzZs1m6dKl+n1q1amFpaUm9evU4fvw4b7/9Nv/9738r+jKMQgpaIYRBmk6bpnYEUQXsvrQbgLY+bVVO8mQuti78o80/6B/cn//s/w+Lji9iY+pGxjUZR59afbDUyH+xTyUqSu0EJqd69eoEBwc/cT8rKytq1qz5wCwHvzd9+nSaN2/O+PHjHxiaYI5kyIEQwiA+7drh066d2jFEJbf34l4crBxo7NFY7SgGa+LRhO96f8cbrd4gtzCXt/e+zYD1A9iWto0iXZHa8cxXWFjxIh4rKyuL7t27s3z5ck6ePElSUhJz5sxh8+bN9OnT55HH1a9fn2effZZZs2ZVYFrjkYJWCCGEScjJzeHor0dp7d0aK42V2nFKxUJjwfB6w9kyYAsvN36ZrNtZTN45mb5r+7I6aTV5hXlqRxSVlKOjI02aNOG///0vrVu3pnnz5ixdupR58+YRGRn52GMnTJjAl19+SVZWVsWENSL5fYgQwiCbfpvkvLdMci6MZF/GPop0RbStYfrDDR5Fa61lXJNxDK0zlKUnl7LyzEre2fsOnx36jBH1RjC4zmCqWVdTO6Z5qF+/+OOJE+rmMBGxsbEPXW9nZ8esWbMe29Pavn17dDrdA+sHDBign/EAeOg+5kJ6aIUQQpiEfRn7AGjj00blJE/P1c6Vic0nsm3gNqa0mIKiKHz0y0d0XdmVmQdmcvHWRbUjClGpSA+tEEIIk5CQmUANbQ1qaM3/AZX7tNZaRtUfxfC6w9l0dhNLTy7lf6f+x7LTy2hg0wCvy140djef8cJCmCrpoRVCCKG6y3cuk3YjjeaezdWOYhRWFlb0C+7H932+54tuX/Csz7McvXeUiM0R/HXbXzn+63G1Iwph1qSHVgghhOoOZh0EoIVnC5WTGJeiKLTxaUMbnzZ8vfFrEp0T2Zi6kWGbhtHVvyuvNnuVQKdAtWMKYXakh1YIIYTqErISAGjp1VLlJBXHy8qLf7f7N6vDV9PFvws7zu3ghfUv8MWRL8gvylc7nhBmRXpohRAG6fDpp2pHEJVYQmYCXg5elWr8rKFqOdfio84f8UvWL7z383t8dvgzdl3cxcwOM6vk/dBbt07tBMKMSA+tEMIg1QICqBYQoHYMUQldu3eNlJwUWni2QFEUteOopplnM1b0WUFk/UiOXD7C0I1D9UMxqqTg4OJFCANIQSuEMMjN9HRupqerHUNUQkcuHwGgqUdTlZOoz9rCmsktJvPpc5+SV5jHX7f9lR3pO9SOpY7k5OJFCANIQSuEMEjchAnETZigdgxRCd0vaBu5N1I5ieno5NeJpT2X4mTtxJSdU9iatlXtSBWvb9/iRQgDSEErhBBCVUcuH8HO0o5gZ/n18u/VcanDkh5LcLV15Y1db7A/Y7/akYSKOnXqhKIobNu2rcT6iIgI/StuAwMDURQFRVFwcHCgXr16jBs3jrS0NAByc3NxcXFh3rx5Dz1Hy5YtGT9+vL6thQsXGu16yptJFrQxMTE0b96cpk2b0qBBA77++msAsrOz6dGjByEhITRo0IDdu3frj7lz5w7Dhg0jODiY2rVrs3r1arXiCyGEMFBBUQHHfz1OQ7eGWGrkOeU/CnIKYn63+dhZ2PHaT6+RfE1+BV+V2dra8tZbbz12n9mzZ5ORkcGJEyeIiooiJSWFdu3akZCQgI2NDUOHDmXp0qUPHHf69GkSEhIYNWqUseIblckVtDqdjuHDh7NkyRIOHTrExo0bGTNmDDdv3mT69Om0bt2apKQklixZwogRIygoKADgww8/xMbGhuTkZGJiYnjllVe4du2aylcjhBDicZKvJ3O34K68LesxQqqH8FHnj7hXcI/XYl/jRt4NtSMJlfzpT3/i5MmTrF279pH7ODo64uXlRWBgID179mTLli20aNGCcePGATBq1Cj2799PUlJSieOWLl1KvXr1aNnSPKfOM7mC9r7r168DcOPGDVxdXbGxsWHFihX6L0jLli3x9PTU99IuX75cvy0oKIgOHTqwTqb8EEIIk3YkW8bPGqKVdyumtJxC+o10/r7r7+h0OrUjCRV4eHjw6quvMmPGDIqKigw6RlEURo8ezYEDB8jMzOSZZ56hTp06JXppdTod3377rdn2zoIJzkOrKAorVqxgwIABODg4cO3aNVavXs3NmzcpKirC3d1dv29gYCDnzp0D4Ny5cwT8bkqh328TQghhmo7+ehSQgtYQw+sO5+jlo2w+u5noM9EMqztM7UiVyoQfJnD+5vkKO59fNT8+7VL6+b2nTp3KvHnziI6OZvjw4QYdU7t2bQDS09Px8vJi5MiRfPnll7z33nsoikJsbCwXLlwgIiKi1HlMhckVtAUFBfz3v/9l3bp1tG3blvj4ePr168fRo0cfmJ/wjz+h/n774356jYqKIioqSv/59evXiYmJKacreLJ79+5V6PkqA7lnZVOe983yt9+AVPavg3yvlU1Z79uBXw9QXVOd+J3xRkhl2spyz9oUtWGfxT5m7Z9FXkoeXpZeRkpnAu7/P/2He1Ref0fd3d25ceMGGk3xL6vzC/IpLCp86nYNlV+QT05OjsH7FxQUkJubi6IojB8/nhkzZhAWFkZ+fj5FRUXk5ORQVFTE3bt3H2j3fm/u7du3ycnJoW/fvsyYMYMtW7bQtm1bFi5cSMeOHdFqtfpjH9WWMd0/544dpZ+qzuQK2sOHD3Pp0iXatm0LFA8t8PHx4ejR4p/iL1++rO+lTU9Px9/fHwB/f3/S0tJKbOvVq9dDzzFp0iQmTZqk/9zX15ewsDCjXdMfxcTEVOj5KgO5Z2Uj96305J6VTVnu292Cu0xbNo3Ofp0J61z17nlZv9cCsgMYtXUU29jG0m5LsdBYGCGd6SqPv6OFhYUkJibi6OiIhUXx/ZsfNr884hmNpaUlNjY2ODk58frrr/PFF1+wZs0arKyssLS0xMnJCY1Gg52dHU5OTiWOvT9etkGDBjg5OeHk5ETnzp1Zs2YNnTt3ZsOGDXz++ecljntUW8ZUWFiInZ0dXbt21X9dDGVyY2j9/Py4cOECZ86cASA5OZmUlBRq167NoEGDmDt3LgDx8fFkZmbSrl07gBLbzp49y86dOwkPD1fnIoSohC7t3s2l380sIsTTSrqWRJGuiHou9dSOYlaaeDQhol4ER389ynenv1M7jvHExDzQOyuKabVa3njjDd577z1yc3Mfu69Op2PhwoW0bt0aDw8P/fpRo0axcuVKvvuu+Huof//+Rs1sbCZX0Hp6evLFF18wcOBAGjduzIABA5g3bx41atRg5syZ7N27l5CQECIjI/nmm2+wtCzuZJ46dSp3794lODiYsLAw5s6di4uLi8pXI0TlcWjWLA7NmqV2DFGJnLpyCoB6rlLQlta4JuOooa3BJ4c+IfN2ptpxjGPSpOJFPNTYsWPR6XRs3LixxPobN26QmZlJWloaW7ZsoWfPnhw8eFDf6XffgAEDKCwsZPLkyQwaNAh7e/uKjF/uTG7IAcCwYcMYNuzBwe6enp4PTCh8n4ODA8uXLzd2NCGEEOXk5NWTAIS6hqqcxPzYW9nz92f+zis/vELUwShmdZAfNqua+3PSvvzyyyXWT548mcmTJ2NnZ4e/vz9dunRh5syZNG5ccmo8BwcHBg4cyFdffcXIkSMrMrpRmGRBK4QQovI7deUUHnYeuNm5qR3FLLX3bU8n305sObuFwbUH08KrhdqRhBHFxsY+sG7MmDGMGTNG//n9N4L90aMe7FqyZAlLlix56LZHtWWqTG7IgRBCiMovrzCPpOtJMtzgKU1rOQ1LjSVRB6NkblpRpUlBK4QQosIlX0+moKhACtqn5Ofox6Dagzj26zFiz8eqHUcI1UhBK4QwiNbPD62fn9oxRCWhfyBMZjh4ai81eglbC1s+PfwpRTrD3h5lFmrVKl6EMICMoRVCGKTjH56QFeJpnLpaXNDKA2FPz83OjWH1hrHk+BK2nt1Kr5oPn4Pd7Kxfr3YCYUakh1YIIUSFO3XlFNVtquNp76l2lErhxfovorXSMvfwXPKL8tWOI0SFk4JWCGGQxGXLSFy2TO0YohIoKCrgzLUz1HOt98ArzUXZONs6M7L+SM7dPMf65ErSszl3bvEihAGkoBVCGCQpOpqk6Gi1Y4hK4GzOWXILc2X8bDkbGToSJxsnFh5bSGFRodpxnt68ecWLEAaQglYIIUSFuj9+VmY4KF8OVg4MqzuMC7cusOPcDrXjCFGhpKAVQghRoe7PcBDqIg+ElbdhdYdhY2HD4uOLZV5aUaVIQSuEEKJCnb56Gq2VFt9qvmpHqXRcbF3oF9yPk1dOEp8Zr3YcUUE6derEW2+9ZdRzpKWloSgKycnJRj1PWUlBK4QQokKlXE+hlnMteSDMSEaGjkSjaFhy4uGvNBXmqVOnTiiKgqIoaLVaWrVqRUxMjFHOtXDhQgIDA0us8/PzIyMjg6CgIKOc82lJQSuEMIhXmzZ4tWmjdgxh5q7cvcK13GsEOwerHaXS8nf0p6t/V3Zf3M2Zq2fUjlN23bsXL0LvtddeIyMjg0OHDtGsWTP69u1bYT2mFhYWeHl5YWFhUSHnKy0paIUQBmk+fTrNp09XO4YwcynXUwCo6VRT5SSV258b/BmAr098rXKSpzBnTvEi9BwcHPDy8iIkJITPPvsMCwsLdux48AHA1157jZo1a2Jvb0/9+vVZvXp1ie2BgYHMnj2bQYMG4eDgQL169fjxxx8BiI2NZfTo0aSnp+t7hGNjYx8YcvDVV1/h6+vLqlWrCAoKwtnZmRdffJHc3Fz9edLS0ujUqRO2trY0adKEFStWoCgKaWlp5X5vpKAVQghRYVJyigta6aE1rgZuDWju2ZytaVv59e6vascRRmBpaYmVlRX5+Q++SMPV1ZXo6GiOHz/OhAkTGDNmDMeOHSuxz6xZs+jTpw+HDx+mffv2REREkJeXR5s2bZg9eza+vr5kZGSQkZFBm0f8du7KlSt8/fXXrF+/njVr1rBu3ToWLFig3x4REQHA/v37+fDDD5kxY0Y53oGS5NW3QgiDHHz/fQDppRVP5X4PbS3nWionqfyG1x3OwayDfJ/4PWMaj1E7TulNnFj8sQJ7aTeFhz90fYdPP6VaQAA309OJmzDhofv0/u1VvZd27+bQrFkPbNf6+elfIZ64bBm1hw8vc878/Hxmz57NzZs3ad++Pd9//32J7b8vHF9++WVWrVrFmjVraNiwoX79Cy+8wMiRIwH4xz/+wZdffkliYiINGjTA0dFRP8TgcXJzc1m4cCGensVv/Bs4cCA7d+5kwoQJnDhxgj179pCUlERwcPEPsJMnT2bMGON8L0oPrRDCIJl795K5d6/aMYSZS76ejNZKi4e9h9pRKr3O/p3xsPdgReIK83wd7rZtxYvQmzVrFlqtFnt7ez744AM+//xzmjRp8sB+X3/9NS1atMDNzQ2tVsvOnTs5f/58iX1+X9zeL1yzs7NLlcfd3V1fzN5v534bSUlJODo66otZgBYtWpSq/dKQHlohhBAVJvV6qsxwUEGsNFYMqTOETw99yk/nfqJ7oDxg9ST3e1kfpVpAwBP38WnXDp927R67T1l7Z0ePHs3EiRPRarWP7D3dtWsXo0eP5oMPPqBDhw5Uq1aNsWPHPjA0wcrKSv/n+38fi4qKSpXn923cb+d+GzqdrkL/nksPrRBCiAohMxxUvBdCXsBKY8Wy08vUjiLKQfXq1QkODn7sUID9+/cTGhrKq6++StOmTalZsyZnz54t1XmsrKwoLHy61yfXrl2bnJwcUlJS9OsOHjz4VG0+jhS0QgghKoTMcFDxXO1cCQsM42DWQfOewksYrFatWpw5c4aNGzdy5swZJkyYUOqhBAEBAWRlZZGQkMCvv/760AfPnqR+/fq0bduW0aNHc/ToUX744QeioqIAjNJzKwWtEEKICiEzHKhjWN1hAHx3+juVk4iK0K9fP0aPHs2f/vQn2rRpQ7Vq1ejRo0ep2ujQoQNDhw6la9euuLu7s2fPnjJl+eabbygsLKRVq1ZMmjSJ6b89VGxjY1Om9h5HxtAKIQwSMnSo2hGEmZMZDtTRyL0R9V3rs/nsZqa0mILWWqt2JMO88oraCUxKbGysQdsUReGTTz7hk08+0a/LycnByclJ//nD5oHV6XT6P2s0GpYuXfrYfSIjI4mMjCyx/d133+Xdd9/Vfx4UFMTOnTv1n//vf/9Dq9Xi4VH+D4VKQSuEMMjTTDEjBMgMB2oaVHsQ7/78LpvPbmZwncFqxzHMuHFqJxBPafv27RQUFFC3bl1OnjzJm2++yfDhw9Foyn+AgAw5EEIIUSFkhgP19Azqib2lPSsTV5boZRPCmO7du8fEiRMJDQ1l7Nix9O/fn48++sgo55IeWiGEQXb+1ltyf2JwIUrj/gwHzzk/p3aUKsneyp7eNXuzMnElJ6+cpL5bfbUjPdn9lxw8YZosYbr69OlDnz59KuRc0kMrhDDIrfPnufWHibmFMJTMcKC+gbUHArAycaXKSQyUklK8CGEAKWiFEEIYncxwoL5Q11BCXUPZfHYzt/Nvqx1HiHIlBa0QQgijkxkOTMPA2gO5W3CXzWc3qx1FiHIlBa0QQgijS7meIjMcmIBeQb2ws7Tj+8Tv1Y4iRLmSglYIIYTRpVxPkRkOTICDlQO9gnpx4soJTl45qXYcIcqNzHIghDBI02nT1I4gzJTMcGBaBld505cAACAASURBVNUexPdJ3/N94veEPhuqdpxH++01qUIYQgpaIYRBfNq1UzuCMFOpOamAzHBgKkJdQ6nnUo9NZzcxucVk7K3s1Y70cGFhaicQZkSGHAghhDCq5OvJgMxwYCoURWFg7YHczr/NlrNb1I4jDNSpUyfeeuutB9ZHRkYSERGhQqInKygoQFGUx762t7xIQSuEMMim8HA23Z/oXIhSkBkOTM/9h8NWJ69WO8qj1a9fvAijyc3NVTtCuZGCVgghhFHJDAemR2utpVtAN45ePqr/gUOYv0WLFtGkSRMcHBwICAhgxowZFBQU6LdHRkYyYsQIpk+fjpubGwMHDqRv37688sorJdpJSEjA0tKSrKwsAFJTU+nTpw9arRYfHx/Gjx/PnTt39Pt/9NFHBAUFYWNjg6+vL++++y4AwcHFv5Xp3LkziqIQGRlptGuXglYIIYRRyQwHpmlAyAAAVieZcC+tKJWioiI+/PBDjh8/zvz581m4cCELFiwosc+6deu4e/cue/bsISoqiqFDh7Jq1SoKCwv1+6xYsYJOnTrh6elJXl4eYWFhhISEcPDgQdatW0d8fDyTJ08GID4+nnfeeYf58+eTlJTEihUr9IXsvn37APj+++/JyMjg448/Ntq1y0NhQgghjEZmODBdzTyaEegYyIaUDbzW7DWsLKzUjqS+Rw1xWLcOgoMhORn69n34PidOFH+MiYFJkx7cXqsWrF9f/Oe5c2HcuFLHmzVrFh999FGJdbm5uQwZMgSA0aNH69cHBQXx6quvsmrVKkaMGKFf7+npyZw5c9Boivs0fXx8+Otf/8pPP/1E165dAVi5ciV///vfAYiOjsbJyYmo3806MWfOHDp37sxnn33GuXPn8PLyokuXLlhaWuLv70+bNm0AcHNzA8DFxQUvL69SX29pSA+tEEIIo5EZDkyXoij0C+7HtdxrxF6IVTuOMMDo0aM5fPhwiSX8d8827N27l+7du1OjRg20Wi3vvvsu58+fL9FG48aN9cUsgIODA71792b58uUAHDhwgIsXLzJgQHEP/rFjxzhy5AharVa/dOvWjby8PC5evEjXrl1RFIVatWrx8ssvs2nTJnQ6XQXcjZKkh1YIIYTRyAwHpq1vcF8+PfQpq5NW0y2gm9px1He/l/VRgoOfvE9Y2JP3KUPvLED16tX1v86/r1q1ahQUFHDz5k169+7N4MGDee+993BxcWHZsmV89dVXJfa3t39wmrYhQ4bw0ksvMW/ePJYvX063bt1wcXEB4NatW3To0IEvvvjigeO8vb2xsrLi6NGj7Nixg61bt/Liiy/yzDPPsP5+b3QFkYJWCGGQDp9+qnYEYYZkhgPT5mbnRnvf9sRdiCPzdiZeDsb9tXCprFundgKzcubMGa5fv87MmTNxdnYGeKB39lF69+5Nfn4+27dvZ9WqVfzzn//Ub2vcuDHr16/H19cXW1vbhx5vbW1Nr1696NWrFxERETzzzDNkZ2fj7u6ORqMpMT7XWGTIgRDCINUCAqgWEKB2DGFmZIYD0zcgeABFuiLWJZtYARkcXLwIg/j7+2NlZcW8efNITU1l/vz5rF271qBjbW1tCQ8P5/XXXycrK4t+/frpt40YMQJra2uGDBlCfHw8ycnJbNiwgSlTpgCwceNG5s6dy7Fjx0hNTWX58uW4ubnh6uqKoij4+fnx448/kp2dza1bt4xy7SAFrRDCQDfT07mZnq52DGFmZIYD09fetz1udm6sSV5Dka5I7Tj/Lzm5eBEG8fDwYMGCBcybN4+GDRuybds2pk+fbvDxQ4cO5fjx4/Ts2RNHR0f9+mrVqhEbG4u1tTXdunWjcePGvPXWW3h7ewPg7OzM8uXLad++PY0aNeLAgQNs3LgRCwsLoPhBtm+//RZvb2/Gjx9fvhf9OzLkQAhhkLgJEwDoXcHjooT5unrvqsxwYAYsNZaE1wpn8fHFHMg8QGvv1mpHKnZ/NoEnjUetIh71tq3fj5GNjIx8YK7XqVOnkpOT88C+f/T8888/8mGugIAAVq5c+dBt7dq1Iy4u7pHtDh48mMGDBz9ye3mRHlohhBBGcX/8rMxwYPr6B/cHZE5aYb6koBVCCGEUMsOB+Qh0CqSZRzN+SP+BnNwcteMIUWpS0AohhDAKmeHAvAwIGUBeUR6bUjepHUWIUpOCVgghhFHIDAfmpVtAN7RWWtYkr1E7ihClJgWtEEIIo5AZDsyLvZU9PYN6cvrqaU5eOal2HKOS70nTVpavjxS0QgiD9F6/XmY4EAa7P8OBDDcwLwNCil93ahIPh504YbQZDhRFQVEU8vPzjdK+KJv8/Hz916a0ZNouIYQQ5U4/ftZJClpzUt+1PiHVQ9icupkpLaZga/nwN0OZO0VRcHZ2Jisrixo1alT6HtuioqIKeVvX09DpdGRlZeHs7CwFrRDCeC7t3g2AT7t2KicR5kBmODBPiqIwIHgAM+NnsuPcDp6v+bx6YWJiij+GhRmleQ8PD9LT00lKSjJK+6bk7t272NnZqR3jiWxtbfHwKNuYeylohRAGOTRrFiAFrTCMfg5aZ5mD1tw8X/N5og5GsSZpjboF7aRJxR+NNOxAo9EQFBREUVHRI18oUFns2LGDrl27qh3jsRRFQaMp+0hYKWiFEEKUu/szHHjae6odRZSSs60zz/k/R0xaDOdvnMfP0U/tSEb1NEWUObn/KtrKqmp8FYUQQlQomeHAvA0ILn44TKbwEuZCClohhBDlSmY4MH+tfVrj7eDNuuR1FBQVqB1HiCeSglYIIUS5khkOzJ9G0dAvuB/Zd7PZe2mv2nGEeCIpaIUQBtH6+aH1q9xj6UT5kBkOKod+wf1QUNSbk7ZWreJFCAPIQ2FCCIN0nDtX7QjCTMgMB5WDj9aH1t6t2Xl+J7/e/RU3O7eKDSAvchGlID20QgghypXMcFB5DAgZQIGugI0pG9WOIsRjSUErhDBI4rJlJC5bpnYMYQZSc1JlhoNK4jn/53CycWJ18uqKn6t17tziRQgDSEErhDBIUnQ0SdHRascQJu7qvatcvXdVZjioJKwtrOlTsw9nc85y5PKRij35vHnFixAGkIJWCCFEuZEZDiqf/iH9AdR7OEwIA0hBK4QQotzcL2hlhoPKo3b12jRwbcDWtK3czr+tdhwhHkoKWiGEEOXm/pRdMsNB5dI/pD93C+6y9exWtaMI8VBS0AohhCg3MsNB5dQzqCe2FrasTpZhB8I0mWRBm5uby/jx4wkJCaF+/fpEREQAkJ2dTY8ePQgJCaFBgwbs3r1bf8ydO3cYNmwYwcHB1K5dm9Wr5S+dEOXJq00bvNq0UTuGMHEyw0HlVM26Gt0Du3P08lH9sBKj6969eBHCACb5YoXp06ej0WhITExEURQyMjL061u3bs3WrVuJj49n4MCBpKSkYGlpyYcffoiNjQ3JycmcPXuWZ599ls6dO1O9enWVr0aIyqH59OlqRxAm7srdK1y9d5XOfp3VjiKMoH9wf9anrGd10mqmtpxq/BPOmWP8c4hKw+R6aG/fvs2SJUv4z3/+o/8J39vbG4AVK1Ywbtw4AFq2bImnp6e+l3b58uX6bUFBQXTo0IF169apcAVCCFE1JV5LBCCkeojKSYQxNPdsToBjABtSNpBfmK92HCFKMLmCNiUlBVdXV/71r3/RokUL2rdvzw8//MCVK1coKirC3d1dv29gYCDnzp0D4Ny5cwQEBDx0mxDi6R18/30Ovv++2jGECUu6lgRAiLMUtJWRoij0C+7HtdxrxF6INf4JJ04sXoQwgMkNOcjPzyc1NZXQ0FDef/99jhw5QteuXTl+/PgDY7L++NaS329/3BtNoqKiiIqK0n9+/fp1YmJiyukKnuzevXsVer7KQO5Z2ZTnfbu+bRsAvzZtWi7tmSr5Xiube/fuEXsiFoALhy9w7eg1dQOZAXP8XnMqdEKDhgV7F1B0usio52r727Mwe3r0KLHeHO+b2qrCPTO5gjYgIACNRsOIESMAaNy4MUFBQZw6dQqAy5cv63tp09PT8ff3B8Df35+0tLQS23r16vXQc0yaNIlJkybpP/f19SUsLMxo1/RHMTExFXq+ykDuWdmU533b9NsrKCv710G+18omJiaG2xa3ccedAT0HqB3HLJjr91rcD3HEXYyjcbvGeDl4Ge9EWi3w4L855nrf1FQV7pnJDTlwc3OjS5cu+p8k0tPTOXv2LHXq1GHQoEHM/e0/1fj4eDIzM2nXrh1AiW1nz55l586dhIeHq3MRQghRxRTpiki5niLjZ6uA/iH9KdIVsS5ZnlMRpsPkClqA+fPnM2vWLBo2bEjfvn1ZsGAB3t7ezJw5k7179xISEkJkZCTffPMNlpbFncxTp07l7t27BAcHExYWxty5c3FxcVH5SoQQomq4UniF3MJcGT9bBbT3bY+bnRtrktdQpDPusAMhDGVyQw4AatasSWxs7APrPT092fbbOL4/cnBwYPny5UZOJoQQ4mEuFVwCoLZLbZWTCGOz0ljRP7g/Xx77kt0Xd9PBt4PakYQwzYJWCGF6QoYOVTuCMGGZ+ZmAzHBQVQyuM5jFxxez7PQy4xW0r7xinHZFpSQFrRDCILWHD1c7gjBhGQUZWCgW1HSuqXYUUQG8HLx4zv85tqdvJy0njUCnwPI/yW9zywthCJMcQyuEEMK8ZBRkEOAYgI2FjdpRRAUZVncYANFnolVOIoQUtEIIA+0cN46d0mMiHuJO/h2uFl6VGQ6qmBaeLQh2DmZd8jpu598u/xOEhxcvQhhAClohhEFunT/PrfPn1Y4hTFDK9RR06GT8bBWjKArD6w3nVv4tNqRsKP8TpKQUL0IYQApaIYQQTyXxWiKA9NBWQb2DelPNuhrfnf7usW/oFMLYpKAVQgjxVJKuJwFQu7pM2VXV2FvZMyB4AKk5qey5tEftOKIKk4JWCCHEU0m6loSNYoOP1kftKEIFI+qNwEKx4OsTX6sdRVRhUtAKIYQoM51OR+K1RLwsvdAo8l9KVeSt9aZ7QHf2ZezjzNUzascRVZT86yOEMEjTadNoOm2a2jGEibl89zLXc6/jZemldhSholH1RwGw9OTS8ms0Kqp4EcIA8mIFIYRBfNq1UzuCMEGnr54GoIZVDZWTCDXVd6tPc8/mbD67mVebvYqHvcfTNxoW9vRtiCpDemiFEEKU2ckrJwHwtfRVOYlQ26jQURQUFbDs1DK1o4gqSApaIYRBNoWHs0kmORd/cPLKSSwUC7ytvNWOIlTW0a8jgY6BrEhcwZ38O0/fYP36xYsQBpCCVgghRJmdunqKms41sVKs1I4iVKZRNIysP5KbeTdZmbhS7TiiipGCVgghRJlcvXeVzNuZhLqEqh1FmIjwWuG427nz9YmvySvMUzuOqEKkoBVCCFEmp68UPxBWz7WeykmEqbCxsGFU/VFcvnuZdSnr1I4jqhApaIUQQpTJyavFD4SFukoPrfh/g2oPwtHakcXHFlNQVKB2HFFFSEErhBCiTE5eOYmCQp3qddSOIkyIvZU9EfUiuHDrAjFpMWrHEVWEzEMrhDBIh08/VTuCMDFHLx8luHow9lb2akcRJmZ4veEsObGEhccW0jOoZ9neIrdOhiwIw0kPrRDCINUCAqgWEKB2DGEism5nkXUni0ZujdSOIkyQk40TQ+oMIfl6MjvP7yxbI8HBxYsQBpCCVghhkJvp6dxMT1c7hjARx349BkAjdyloxcONDB2JlcaKL499iU6nK30DycnFixAGkIJWCGGQuAkTiJswQe0YwkQcvXwUQHpoxSO527vTP7g/x349xoHMA6VvoG/f4kUIA0hBK4QQotSOXD6Cg5UDQU5BakcRJiyyQSQWigULji5QO4qo5KSgFUIIUSoFRQWcvHKSBm4NsNBYqB1HmDC/an70CurFgcwD/JL1i9pxRCUmBa0QQohSSbyWyL3CezLcQBhkdKPRaBQN84/MVzuKqMSkoBVCCFEq93vamno0VTmJMAdBTkH0COzBzxk/czj7sNpxRCUlBa0QQohSSchKQKNopKAVBhvTaAwKCvOPSi+tMA55sYIQwiC9169XO4IwAUW6Ig5mHaSeSz201lq14wgzUdO5JmGBYWxN28qxy8do6N7wyQedOGH8YKLSkB5aIYQQBku9nsr13Os092yudhRhZl5q9BJAhfTSdurUibfeeuuB9ZGRkURERBj9/GVRUFCAoijExsaqHcUsSUErhDDIpd27ubR7t9oxhMoSshIAaOHZQuUkwtyEVA+he0B34i7EGTaWNiameDERubm5akcQjyEFrRDCIIdmzeLQrFlqxxAqS8hKQEGhmWcztaMIMzSuyTg0ioZPDn3y5LeHTZpUvBjJokWLaNKkCQ4ODgQEBDBjxgwKCgr02yMjIxkxYgTTp0/Hzc2NgQMH0rdvX1555ZUS7SQkJGBpaUlWVhYAqamp9OnTB61Wi4+PD+PHj+fOnTv6/T/66COCgoKwsbHB19eXd999F4Dg317z27lzZxRFITIy0mjXXhlJQSuEEMIghUWF7M/YT12XujjZOKkdR5ihms416VOzD/GZ8ezL2KdqlqKiIj788EOOHz/O/PnzWbhwIQsWlHwBxLp167h79y579uwhKiqKoUOHsmrVKgoLC/X7rFixgk6dOuHp6UleXh5hYWGEhIRw8OBB1q1bR3x8PJMnTwYgPj6ed955h/nz55OUlMSKFSv0hey+fcX34/vvvycjI4OPP/64gu5E5SAFrRBCCIOcuHKC67nXaVejndpRhBl7pckrWGos+eQXA3ppn8KsWbPQarUllm+//Va/ffTo0XTt2pWgoCB69uzJq6++yqpVq0q04enpyZw5c6hTpw4hISGEh4dz+/ZtfvrpJ/0+K1euZOjQoQBER0fj5OREVFQUderUoWXLlsyZM4fFixdTWFjIuXPn8PLyokuXLvj7+9OmTRv9mF43NzcAXFxc8PLywslJfmgsDSlohRBCGGT3xeIx1FLQiqfho/VhcO3BHL9ynB/P/2i084wePZrDhw+XWMLDw/Xb9+7dS/fu3alRowZarZZ3332X8+fPl2ijcePGaDT/Xyo5ODjQu3dvli9fDsCBAwe4ePEiAwYMAODYsWMcOXKkRBHdrVs38vLyuHjxIl27dkVRFGrVqsXLL7/Mpk2bjFrUVyVS0AohhDDI7ou7qWZdjUbu8oYw8XRGNxqNnaUdnx36jMKiwicfUAbVq1cnODi4xFKtWjUAbt68Se/evQkKCuL777/nl19+Yfr06eTn55dow97e/oF2hwwZwurVq8nPz2f58uV069YNFxcXAG7dukWHDh1KFNFHjhwhKSkJb29vnJycOHr0KJ9//jnW1ta8+OKL9O3b1yjXX9XIPLRCCINo/fzUjiBUdPXeVY7/epxuAd2w1Mh/HeLpuNm5EVEvgi+Pfcna5LW8UPuFB3eqVcto5z9z5gzXr19n5syZODs7AzzQO/sovXv3Jj8/n+3bt7Nq1Sr++c9/6rc1btyY9evX4+vri62t7UOPt7a2plevXvTq1YuIiAieeeYZsrOzcXd3R6PRlBifKwwnPbRCCIN0nDuXjnPnqh1DqGTXhV3o0MlwA1Fu/tLwL7jYuvDpoU+5nX/7wR3Wry9ejMDf3x8rKyvmzZtHamoq8+fPZ+3atQYda2trS3h4OK+//jpZWVn069dPv23EiBFYW1szZMgQ4uPjSU5OZsOGDUyZMgWAjRs3MnfuXI4dO0ZqairLly/Hzc0NV1dXFEXBz8+PH3/8kezsbG7dumWUa6+spKAVQgjxRFvTtmKpWPKc/3NqRxGVhIOVA+ObjufKvSssOraoQs/t4eHBggULmDdvHg0bNmTbtm1Mnz7d4OOHDh3K8ePH6dmzJ46Ojvr11apVIzY2Fmtra7p160bjxo1566238Pb2BsDZ2Znly5fTvn17GjVqxIEDB9i4cSMWFhZA8YNs3377Ld7e3owfP758L7qSk98bCSEMkrhsGQC1hw9XOYmoaDm5Oey7tI/WPq1lui5RrvoH92fZqWUsPbmUwXUG4+Xg9f8b7/9GaNy4Urf7qLdtffXVV/o/R0ZGPjDX69SpUx+67x89//zzj3yYKyAggJUrVz50W7t27YiLi3tku4MHD2bw4MGP3C4eTXpohRAGSYqOJik6Wu0YQgU/nPuBAl0BYYFhakcRlYylxpKpLaaSW5jLx7/8Yd7VefOKFyEMIAWtEEKIx9p6diuWGhluIIyjTY02tK3Rlo2pGzn+63G14wgzJQWtEEKIRzp/4zz7MvbRoUYHHK0dn3yAEGUwpfkUNIqGD+I/kHlZRZnIGFohhKgE7uTf4XD2YVJyUrDUWBLoGEhTj6bYWj586iBDRZ+JRoeOYfWGlVNSIR4UXD2YgSEDWZG4gh3ndtAtoJvakYSZkYJWCCHM2O3823x14iu+PvE1dwvulthma2FLeK1wIutH4udY+nmE7+TfYU3yGmo61eQZr2fKK7IQD/VKk1fYdHYTH8Z/SLsa7bBTO5AwK1LQCiEM4tWmjdoRxB+k5aQx4ccJpN1II8gpiBdCXiDUNZRCXSEnr5wkJi2GFYkrWJ28mpGhI3mp0Us4WDkY3P6a5DXczLvJ35r+DUVRjHglQoCrnSvjmoxjVvwsFh1dyPj69SEnB3bvhrZtQb4HxWNIQSuEMEjzUszRKIzvxJUTjN42mtv5t5ncfDIRoREl3uDV2rs1f67/Z+Iz4/kg4QMWH1/MhpQNTGs1jR6BPZ7Yfk5uDp8f+RwPew/6BsurOUXFGFZ3GLv3fkvv56eiu1KEYm0DXbpAUBDExEBAgNoRhYmSh8KEEMLMpOakMnb7WPIK85jbZS6RDSIf+jpaRVFo5d2K6N7RzGg9g/yifKbunMq/9/2b/ML8h7T8/+YfmU9Obg6vNXsNO0v55a+oGJaKBR/NPI1fVi5KXj7cugV5eZCSAj16gDwwJh5BClohhEEOvv8+B99/X+0YVd7NvJuM/2E8N/NuMrvjbINeRWuhsWBwncGs7buWll4tiT4TzYsxL5J9J/uh+/947ke+PfUtDd0a0rtm7/K+BCEebc8e7M5nYln0h/UFBZCaCnv2qBJLmD4paIUQBsncu5fMvXvVjlGl6XQ63tn7Dudvnmdqy6l09OtYquNd7VxZ0G0Bo0JHcfjyYQauH8jG1I0lpknafXE303dNp7ptdWZ2mIlGkf8mRAVKTgYrq4dvs7Yu3i7EQ8gYWiGEMBPfnf6O7enb6R7QnWF1yzaNlqXGkiktp9DIvRH/2vcv3tj1Bp8f/pymHk3JuJ3BgcwD2Fna8clzn+BXrfQzIwjxVIKDi4cYPExeXvH227crNpMwC0/1o/evv/4qEyALIUQFOP7rcT5I+AD/av78o80/nnrWge6B3Vnfbz0jQ0dyr+Ae61LWcTDrIN0CurG271oauzcup+RClELbtsUPgFmW7G8rsrCAmjWLtwvxEAYXtO+99x7Lli3Tf96vXz88PDzw8vIiISHBKOGEEEJAbmEur8e9jgYNszvNRmutLZd2nW2dmdpyKtsHbWf30N388qdfiOoUhY/Wp1zar0hFRTru5hVKJ4u5U5Ti2Qxq1QJFQaco5FkqXPSyI3/zRpm6SzySwUMOFi5cyIoVKwDYvHkz+/fvZ//+/URHRzN58mR27txptJBCCFGVfXn0S87dPMek5pOo61K33NvXKBqcbJzKvV1ju5NXwJI9aWw+lsGJSzcAcLa3opGvM70aeNGrkTeOto8YjylMV0AAnDoFNWui5OWx+d/DmMFWXrmxjbGMVTudMFEGF7TZ2dn4+voCsGHDBoYMGULLli1xdXWlSZMmRgsohDANIUOHqh2hSjqbc5ZFxxcRUj2EiNAIteOYjH2pV5i84ggXr9/Fyc6K7qGe2FtbcO7qHfalXCEu8TLvbTzJC818GdUmkGCP8unVFhVEUWDKFAB6/+klvt6YwoKjC3jO7zmVgwlTZXBB6+Pjw+nTp/Hx8WHLli3Mnz8fgDt37mBpKc+WCVHZ1R4+XO0IVdLHv3xMQVEBM1rPwEojvY0AP53OZsz/DmKlUXi3TygjWgdgZfH/I+hu3Mtny7EMlv6czjf7ipdWgS50DfXgubqe1HJ3kDefmYNx4wCwAv7V7l+M2DSCN3e/yYvWL6qbS5gkgyvRv/3tbwwcOBBvb2/s7Ozo2rUrAHFxcTRs2NBoAYUQoqo6nH2YH879QBf/LjT1aKp2HJNwMP0qY745iIONBf/76zPU93lwqISjrRVDWvozuIUfCenX+GpvGj+cyuJA2lX+s/k0ga72PFfXk1ZBLjTxc8bLyVaFKxGlUd+1Pi81eonPj3zOdoft9EbmRzZY/l3s716CjKPgFgJWlfNFKQYXtK+99hpt27bl/PnzdO3aVd8rGxAQwH/+8x+jBRRCmIadv/WWdJw7V+UkVcfHv3yMRtHwt6Z/UzuKSbh+J48Jyw6h0fDIYvb3FEWhZaALLQNduJdfyM8pV/jhdBY/nspm8Z6zLN5zFgCPajb4ONvhbG+Fs50VtlYWWGgUrCw0ONlZEeBqT4fa7rhpbSriMsV94eHFH9evB2B0o9HEno/lx6s/cvTyURq5N1IxnBm4fg5iZ8LJtbTPuwWHAY0lBLSFdq9Bzc6V6iG7Uo0VaNmyJS1atODixYs4ODhgYWFB797yU5IQVcGt8+fVjlClHMo+REJWAn1r9aWmc02145iEN9cc41LOPWa+0PCJxewf2VpZ0LmuB53reqDrqyMp+xaHz13n8IXrHLuQQ2bOPU5n3uBe/h9fUVVMo8AzQa70buRNn0Y+ONnL8A+jS0kp8amVxop/t/s3g9YP4u+7/87y55djb2WvUjgTd/x72PAa5N4An2ak4k/NkDpw6TCk/ABnd0Kd3hD+CTi4qZ22XBhc0N69e5cpU6awaNEiCgoKSExMpGbNmkydOhUfHx8mTpxoiMVGlgAAIABJREFUzJyiqrhzFTIOg6KB6oHFixBV0OJji1FQeLGhjBcE2JV0mc3HMgmr78ngFk/3wgdFUajtWY3antUY3LJkW/fyC8ktKKKwSEdBYRFXbudxOvMGW45lEpt4mZ9Tr/CvTSfp16QGEa0DaFDD/GaHMGch1UPoWa0nG29sZGb8TP7R5h9qRzI9h5fB2rGg9YSBSyCkK0kxMdTsHFa8/Vo67HgXTqyGzw/C8GjwMf8hTQYXtNOmTeOXX37hxx9/pHv37vr1HTp04O2335aCVjydCwfhp38X/+T4e74tof1k5C3NoipJupZE7IVYuvp3paaT9M7mFxbx3oaT2FhqeKt3qFEf6LK1ssDWykL/uYejLfW8Henf1JdbuQVsOZbBt/vPER1/nuj487Su6cIrnYJpH+ImD5pVkA72Hbha7Sqrk1bTxqcNYYFhakcyHYkxsG4cOPnDi1vAyffBfaoHwKAlULsHrJ8AS3rDkKUQ3LXi85Yjg6uENWvW8Nlnn9GmTZsSf2lDQ0NJlncri6exbz4s7AIpP0K9PtBzFvT6EBoPh8zj8N1QGiTPhTx53aGoGv536n8A/LnBn1VOYhpWJlwgKfsWYzrUxM9FvV8xa20sGdTCj7Xj2rJhfDsGNKtBQto1Ri4+QPhne9h6PIOiInmxg7FpFA3/bvdvXGxd+Mfef3Dp1iW1I5mGGxmw5mWwdYKRax9ezP5e4yEwagNYWsN3wyHlp4rJaSQGF7Q5OTk4OT34q5Vr165hZSVjiUQZxb4PW18H11owdi8M+R88MwZajYb+n8PfDkFwN2pc3gnfDIB7OWonFsKocnJz2Jy6mUZujeShF4p7Z+f+lIyzvRUvdayldhy9hr5ORA1uQuzUTv/H3n3HVVX/Dxx/3cvegiggU5YIKC7c2xR3maPc5azMhuX6VeaovlZqZVGZmppmbnMLbnPjFhQB2SKyh2y49/fHUUpFWXdxPc8e53HljnPe3C73vu/nvD/vD+M7OBNxP5e3Nlymz/cn2Xs9SUxslczayJpFnRaRW5LL7JOzKZGVqDsk9ZLLhTKDggwY/KPwmVoVTu1g7N+gZwh/jYSYf5QbpxJVOaHt0aMHK1asKP9ZIpFQXFzMF198Qe/evZUSnEjLXd0Ix/8Hdi1g4iGw8X76PuZ2MHorMXYDIeEc/PEKFD1QfawiWs6aRctZs9Qdhtb7O+pvCssKed1LXMgCYMflRO5mFTC5iyumBprX89zB0pgFL/tyanZP3u7uxr2sAt7deIUBP57iyK374lK8tbFsmbA9Q1eHroz1HsvV1Kv8cOkHFQamgcJ2QvQxaDlGONNZHY1awNidoKMHG1+DuLPKiVHJqpzQ/vDDD+zatYu2bdtSVFTElClTcHNz48aNGyxZskSZMYq0UeIl2P0eWDjCqC1gbPXs+0okRDiPhc4zIOkybJ8IZaWqi1UEQKPOnWnUubO6w9BqMrmMLbe3YGlgSR+XPpU/QMvJZHJ+Pn4HCyM9xnVwVnc4z9XAzIDZfb04OasHU7q6Ep36gInrLvLqL2c4eydd3eHVTQEBwvYcH7b6EL8Gfqy7uY7g2GAVBaZhSgrg0OdgYAEvLazZPuxbw5jtQhuvP4dDQohiY1SBKie0jRs3JiwsjOnTp/PBBx/g4eHBp59+yrVr13B0rN2MU9ELpjgPdkwW/nBe3whmNpU/RiKBXvOg+esQcRAOf678OEUiFbt8/zLxufG84v4KBjpiz9MTEanEpeczpr0TZoZ1o7StvqkB/9e/KSdn9WBse2dC72YzcuU53vnzEgkZ+eoOT+vo6eixtNtSrAytmHdmHjHZMeoOSfXO/QLZ8dB9NpjUr/l+HNvC6K0gL4MNQyHpiuJiVIFqTR3X09Nj7NixfPPNN/z8889MnToVU1PlrY+9YMECJBIJoaGhAKSkpNC3b188PDzw9fXl1KlT5ffNz89n5MiRuLu74+npyY4dO5QWl6iWgj+FjDvw0nywq0aNoEQi1AY5toOzP8GNbcqKUFSBfYMHs+9Ro3ORUuyJ3gPAILdqnjLUUuvOxqIjlTC6nWaPzlbExtyQRa/4cmRGdwJ8bNh/I5mXlp3g+8MRFJdW3OtW9AQfH2GrhI2JDd90/YaC0gJmHJ9BfskL9MWhOF/4PLRwAv/Jtd+fc0cYtRnKioQSv3vXa79PFalyQdLvv//+3NsnTFBsr8TLly9z7tw5nJycyq+bM2cO7du35+DBg4SEhDBs2DDu3LmDrq4uS5YswcDAgKioKGJiYujQoQM9evTA0tJSoXGJainuLFz8HVy6QLu3q/94XX0Yvg5+6ya0G2ngBba+io9TJFKxwtJCgmKDaGrVFA9LD3WHo3axaXkcv51KXx9bGtWru0t1OtU3ZsXYNpyKTGPBnjC+PxxJUNh9lg73U3doWqWdXTumt5zOD5d/YPY/s/mhh7DKnta7/Afkp0P/ucLnoyI07iqcPf1rJKwdACP+ALceitm3ElU5oV20aNFjP5eUlJCcnIyhoSENGzZUaEJbVFTEtGnT2LhxIz16/PskbtmyhZgY4XSCv78/NjY2nDp1iu7du7N582bWrl0LCOURXbt2ZdeuXbzxxhsKi0tUS6XFsPcD0NGHgd+DtIZvNuZ2wh/Y2gGweQxMOQZG4hcXUd12LOEYeSV5DHYTR8EBNoUIK9Npeu1sVXX2sGbfe1346VgUgceiGPzTKQY4S+gtkyOViv1rFWGi70SisqLYF72P7y9/z4zWM9QdknKVFsOZ5WDSUJgMpkjuvWDMNtg0Bv4cBi8tgA7TNHqp3ContI8Syf9KTU1l8uTJjBo1SqFBzZs3jzFjxtC4cePy69LT05HJZDRo0KD8OhcXF+Lj4wGIj4/H2dm5wtuetGzZMpb9Z+ZkVlYWQUFBCv0dnqewsFClx9MUzkl78EoNJ8phOHcu3QHuVPqYRyp6zhydx+Mds5rUFUO47DVbWF1M9BhFvtYePBC6S2j7a1ddf5+/Z/yOFCmGsYYEJdS951iRz5tMLmfTuTLqG0J21EWC7mjuh2h1+QCzW0lZe6uMXTFyYpYeYIK3FCNd7fkdFaXTw/ec00+8rp73Wuss70yoXihrQteQH5+Pv7G/0uNUF9u0M/jl3CXCaSQxR08+9741/fs08fqcFreXYBr8CSkXtnPLdRKFBhq6VK68lq5fvy53dnau7W7KnTlzRt6jRw+5TCaTy+VyubOzs/zGjRvytLQ0ubGx8WP3HTZsmHzdunVyuVwuNzU1laekpJTf9vHHH8sXLFhQpWPa29srKPqqOXjwoEqPpxEepMnlXznK5ct85fLigmo/vMLnTCaTy3e+LZd/bi6XH/lCAUFqH0W+1vYOGiTfO2iQwvanqdTx95manypvvq65fNrhaSo/tqIo8nk7FZkqd569V/7twXCF7VPTFBSXykcs2yd3nr1X3mvpcfmdlFx1h6R5vL2F7QmVvdZS81Plvbf2lrf4o4X8XNI5ZUWnfr/3l8sXWsvlD1IrvWut/j4Lc+Ty7VOEz9ov7OTyQ/Pl8tyUyh+nBM/L12o9pJWWlkZOTo4icmsATpw4QXh4OI0bN8bFxYXExEQCAgK4cOECIIwKPxIXF1deY+vk5ERsbGyFt4k0wInFUJQNL30uNHBWBIkEBiwFOz84+Q2E71PMfkUiFdsXvQ+ZXCZOBntox+W7AAxpZa/mSJTHUE+HN5tK+XRAU6JTH/BK4GkuxWWoOyytYG1kzU+9fsJIx4j3jr5HWHqYukNSvJRbEHcKvF8BEyWPmBqYwasrhLZe5o3g1DJY5iUsdnT7oHKPXQ1VLjmYN2/eYz/L5XKSk5PZsWMHQ4cOVVhAc+bMYc6cOeU/u7i4sHfvXnx9fRk+fDiBgYHMnz+fkJAQkpOT6fywL+aj29auXUtMTAwnTpzg119/VVhcolpIjYCQ1eDgD76Ke60AoGckrC72W3fYMRUmH4UGnoo9hgiArj/+qO4QtNaeO3sw0zOju2N3dYeidoUlZRwMvYefgwVuDZTXRUcTSCQSJnVxxdPGjLc2XGLMqgusGt+GTu4aekpX1XbtqvFDPS09+bHXj0w9NJW3D73Nun7raGzRuPIH1hUX1wiX/hNVd0z3l2Daebi1B67+KSyV69YD6Ku6GJ6jyiO0//zzz2PbmTNnyM3NZdGiRQQGBiozxnJff/01Z86cwcPDgzfeeIP169ejqyvk5DNnzqSgoAB3d3cCAgIIDAzEyuo5zfpFqnNontDXrs+Xyikor+cEw9ZASZ7QOy9HXNdbGcycnTFz1o4JOpokKjOK25m36ePSR+w9C5yKTCOvuIwBze3UHYrKdPVswIZJ7dDTkTBp3UVCYsWRWgDc3YWthlrbtGZJtyXkFOcw9dBUkvOSFRicGpUWw40tYN1EaGOpSlId8HlF6Fc7MxJajVPt8Z+jyiO0x44dU2Ycz/TfMgIbGxuCgyteCcTExITNmzerKCpRlUWfgIgD4DNEWDNaWVy7CT1qd00Teue9eaB2DaZFT8mNiwMQk1oFOxR3CIB+jfupORLNEBQmJB0BPrZqjkS1WjlZsm5CW8asOs+ENSFseasDTe3M1R2WekVFCZe1SGq7O3ZnYaeFfHLqEyYETeD3gN+xNanjr62oQ1CQCR3fU2/XAQ3rLiROCxcpj6wMgj8R2nS9NF/5x2s5BgK+grTb8OdQKFRcbbdWkMshJVxYkCJkFVzbJNRhVXGt+ZPTp3Ny+nQlB/niCY4LxtLAktY2rdUditqVlsk4fOs+XrZmONc3UXc4KtfSyZJV4/0pLC1j4toQUnIL1R2Ser38srDV0mC3wczrMI+E3AQmBU8iJT9FAcGp0bVNgASaj1B3JBrluSO0jo6OSKqY/T+rRZboBXZtEyTfEL5FWrqo5pgdpkFBljBJ7M9hwoonGvYtUuXKSoR6q5BVQrL/JAtHaDcVWr8JBtpds6hporOjicqKYqjHUHSlVT5hprUuxGaQmV/C+I4u6g5FbTq41Wfxq835aOs1pq6/xOYpHdDXFceeamu453DkcjmLzi1iYtBEfg/4nQbGDSp/oKYpyBSWf2/cBSwc1B2NRnnuO+gXX3yhqjhE2qY4D44sBOP60OUj1R67x/8Jy/ad/gHWDICxO8Csjp9iqiGL3Aj4+RNIjxQSe//J4NRemBVbkAUJ54UR2+BP4fwKGPS9UPhfF5UWQWIIpN8RvkA5dwIdzU4SD8cdBqCPcx81R6IZDt28D7x45QZPGtragdv3c/ntZDRfHwzns4He6g5JK4xoMoIyeRlfnf+KCUETWNN3DdZGdWwC3q29UFYMzV9TdyQa57nv9uPHj1dVHCJtc+ZHeJAM/ZeAUT3VHlsigd4LwcgKDn8Oq/vA2J1Q3021caiTXA5nA2kb9jno6EGPT6DjdKErxH/5vAK9PodLa+DoF8Kkuq6zoPvcmq/kpmpyOVzfLMSfnfDv9aY2MHSVsIyjhjoUdwhzfXP87bS3+Xt1nLidip2FIV62ZuoORe1mBjThSnwmq0/F0NGtPr2a2qg7JK0w0mskMrmMxRcWMyFoAr/1/q1u1dTe3AVSXfAaoO5INE61P7Hkcjn37t0jPj7+sU2kQqXFkHABzv0KxxfDySUQthMepFb+WFXIuSeMjlp7Cqex1aXzBzBouZDkrOoFUYfVF4sqyWRwYBYEf0K+gQ1MPgbdZj2dzD6iZwjt34Z3zkKjlkK5xrY3hdeZpisrEZZT3jkVih9A9/8T1iDvNgdKCmD9ELi+Rd1RVig+J57wjHB6OvVET6qn7nDULj49n+i0PLo3aVDlUjdtpqcjZfnIlpgZ6jJnxw0y8+rA32MdMbrpaOa0nUNMdgxj9o8hOita3SFVTUEmRB8H1+5iKV0Fqnw+Lj09nXfeeYe///6b0tLSp24vKytTaGCiCuTcg3OBcOVPKKigrYtECq49oOtMcO6g+vgeOTwfSvKh9yL1n/JtPV4oN9g+GTYMg/bvQM9PQF9LJ5zI5bD3fbj8B7h04XyDifSyqeLpynpOQneInW/Bzb+F/4cj1ituIQxFK8yBrePhzlHhDX742n/f5L0GCKPPG4bBrnfBxkfYNMij7ga9nXurORLNcCJCmKjTzbMO1jUqiZ2FEQsG+zBjyzUW7Anj+9dbqjskrTG66WgsDCz47NRnjDs4jp96/kSLhi3UHdbz3T4AshJhMQXRU6o8Qjt9+nTu3bvHiRMnMDIyYt++faxdu5YmTZqwbds2ZcYoKi2CY1/B8pbCqXwjS+GU8Nid8NZpmHQUBv8EHgHCt7c1fWH7JKGOVdXizsL1TUIdpmeA6o9fEc8AmHIM7FsLXwh+7iA0hNZGx74Uklm3XjB6G6W6xtV7vJ4RDPsdWoyGyGDYMVnoVgEM2L2bAbt3KyHoGshOhN/7Cslsq3EwetvTIxYNm8KIP4QeyNsmCiO2GuRQ3CHM9Mxob9de3aFohBMRqehKJXQUFxV4zJCW9vTyasjfV5M4F52u7nBUKyxM2JRkoOtAfuz1I8VlxUwOnszBGM1Z9apCYrnBc1U5oT18+DDff/897du3RyqV4urqytixY1m2bBmLFy9WZowvtruXYEU3OPE1WDrD63/B9EvQfQ649QRbX3BoDa3GwqhN8G4IePaFG1th1UuQfVd1sZaVwv6ZINWDvl+rtz/ek+q7wcRg6LsY8lJh/SuwabTQhUFbnP8NTn4LjVoJiVxNR1alOsIXJJ9X4dZu2Dejyq29VCLpKqzsBSlhQv3voOVCnXBFHFoLfyupt4RJbxri7oO7hKWH0d2xO/o6+uoOR+2KSss4cyedVs6WmBuK5Rf/JZFImDfIG31dKfN3h1FaJlN3SFqls31nVvdZjZm+GTNPzuSHyz9QJtPAM85FD4SBGOdOYCwuGlWRKie0JSUlWFoKIyANGjTg7l0hUXJ3d+fGDS1KCjRFSSEc+lxIStOjoMen8NYp8Or//ESxvhuM3AS95kHKTVjTDzLjVBPzpTVw/4Yw+ci65o2wlUaq87BW9Bw0HQThe+HXzrB5bN1PbEN3CHWz9d2FFVxq235LKoUhK4QSlktr4diXJJ06RdKpUwoJt8aSb8Afg4VasmFroMuMyr84dXwPLJzg9Pca05v4UXcDsdxAcDU+i/ziMrqIo7MVcq5vwpQuroQn5/LXhRdozkpQkLApWbMGzdg0cBPNrJux6sYq3jv2HtlF2Uo/brVEHxO694ijs89U5YS2ZcuWXLp0CYCOHTuyYMECdu/ezdy5c2nSpInSAnwhxfwDK7oIH8B2fjD1JHSb+exRqCdJJEKrrEHLISse1g0S6m+VKS8Nji4Cc3vo+rFyj1Vbls7w2gbheW3SXxiF/LUzbHxNWHigrok+DjumCLP6x+wQWnIpgq6+8Dw1agUnv+XKgrlc+eYbxey7JtLvCJO8SgqEsxG+r1btcboGwqS4gkw494tyY6yiQ3GHMNY1pqN9R3WHohHO3BFOpYvlBs/2Tg83bM0NWRIcQcaLMkFsxgxhU4GGxg1Z03cNg90GczLxJK/ufpWQ5BCVHLtKbh8QLj37qjcODVZpQvtoAtiXX36JhYUFAIsXL0ZXV5fRo0cTGRnJqlWrlBvliyA/A65tFpZtXTcQMmOF06kTD0NVJ/U8qfV4GPQDZMUJiUC+EtcHD/4UCrOhzxd1Z8KVnR+M/AumnADvl4Vm1b90hH0fQ14dqVVLugqbxoCeMYzZLiTrimRgKtSnWrnBg/vCaS91yE6EP16G/HShxtetZ/Ue7zcSLBvD+V+Esx9qlJyXzLXUa3Rz6IaBjoFaY9EUZ++kY6KvQ3MHC3WHorGM9XX5vwFNyS4oYWlwBQukiGrNQMeALzp9wcKOC8ktzmVi0ESWX15OSVmJegOTlQmfTza+in+P1yKVJrS2tra899576Onp0bu3cHrM3t6ew4cPk5uby9WrV2ndWlyyscpCd+B076AwUnT0S2Fk7eeO8I0r7JwCMSfAd6hQC9tlRu27BLQeL3QbSL0Ffw5XTkISvh+u/SVMBPMZovj9K1ujFkLN6aSj4NAGQlbCjy3hbKBmt67KiBZWQysrFhJzW1/lHMekPozZBhIdyLkLiZeUc5xnyUsTvuhlJ8DLgUK5SHXp6ELbKcIo7c1dio+xGo7EHwGgt4tYbgBQUFzGlYRM2ja2Qk+njvQ+VpNBze3wd7FkU0gCMWlqmPT7ApBIJAzxGMKWgVtoWr8pK2+sZPie4VxJuaK+oBJDhC/zTfqpL4Y6oNJ3jwULFnDhwgXatm2Lj48P3377LUlJSaqITTudWkbT2N/h4Byh3+f1zcLIps8QYTT1owhhBEqRS8V2ek8oQbh7ETaNUuwI1YNU2PM+GFrA4B81ayJYdTm0hglBwvNvYAFB/wc/txdO9WjSpCiArITHRyxdOin3eFau/y6zuHGEcPpfFQqzhbML6ZHCRMMWo2q+L7/XQdcQLv6uuPhqIDg2GCNdIzo1UvL/szriYlwGJWVyOrjVV3coGk8ikTCrrxdlMjk/HI5QdzhazcXChQ39NjCtxTTic+MZd2AcC84uUE9tbaTQ4g8PDekcpKEqTWinTZvGuXPnCA8PZ8iQIfzyyy84OzsTEBDAX3/9RWGhek/f1TmDlnPB+3N4Y78wOWl2HMwIg+FroPUbYKqkHow9P4M2E4UR4O0ThY4EtSUrE/aVlwL9l4J5o9rvU90kkocj5BeE5yw3Gf56XUgek0PVHZ0gJ0koS8mKF75ENB2omuPqGQk10gUZsG6wUBajTIU5wsplydeFBRPav1W7/RlbCZ0bEs7BfeW1Anqe+3n3uZJyhS72XTDWq2ZLNS1VXj/rJtbPVoW/ixVdPKzZdS2JyPu56g5Hq+np6PGW31vsGLwDf1t/tkVsY9DOQWyP2I5MrsJuE1GHhZUv7Vup7ph1UJXP73h6evLFF18QHR3N4cOHcXR05J133sHGxoZJkyYpM0btYt+KTAsfYUStYVPVLQsrkQjL0PoOE2b373lPWFGqNo59KSTIrd+E5sMVE6em0DMSJre9dxlajIGYk8JEvT3vq3dFttxkYZJfZqwwot9yjMoOberoiKl7U6H7Qc5dIY6shMofWBNFuUI5RWIIdPpAmNSlCK0fLud97S/F7K+aDscfRo6cABdxpOWRkJgMzAx1aWpnru5Q6owZvT2Ry+H7w5HqDkW53NyETc1cLFxY3Wc1X3X+Ch2pDvPPzmf0vtHcSFVBd5wHKXDvqjBvQKqj/OPVYTUqWOrWrRuBgYF899136OjosGbNGkXHJVIGqRSG/AoefeDqnxA0t+an0i+ugX+WCjPg+2pxH2IzW3glEKYcB6cOQgurH1sJS/uqur42IwbW9BfauPVfIozoq1C3wEC6BQZC8xFCLWt5Bw0FlyAVPRDqvRPOCy3gXpqvuFIWx3ZCC6+wv2v/ha4GHpUbdHHoovJja6LCkjKuJ2bT2tkSHWkdLldSsZZOlvT0asi+G/e4maQZreiUYvduYdMAEomEQW6D2PPKHsZ7jyc8I5xR+0fx+ZnPyShU4oTrO0eFS/eXlHcMLVHthPbkyZNMmjQJGxsbPvjgA4YMGcKxY1q66pI20tGD4euE5sznf4W/34HqzuC8tklotm/pAqO2aO7SqIrUqAW8sU+YPGZkCYfmwS8dIPKwao5/9zKs7i1MBBuwFNpOVs1xn6XlaBj4PWTGCElt7n3F7DfnHqwdAPFnof00YUKjIuuyJRJhSdzsBGH0V4Xu593ncsplujp0xUjXSKXH1lTXE7MpLpPh7yI2iq+uGb09AfhOrKVVKVN9Uz72/5htg7fRzq4dOyJ3MHDnQLZGbEWujLkWUQ8/Y6rb2eUFVKWENiIigs8++4zGjRvTs2dP4uPjCQwM5N69e6xevZquXbsqO06RIukbC833PQLg2kahHrIqfWplMjj1HeycCmaNhJ6nyqr51UQSidDea9rD+tqcJPhzKGx8XUg0lSUiGNYOFE7Dv7Ye/NVT4hOxcSMRGzf+e0WbN6Hft8KI8bpBQmut2og+ASt7CKfXus2BgC+VM8nQd6hwGbpd8ft+jkNxwsQOsdzgXyGxwsiWmNBWn6+9BX28bTh087721tIGBgqbBnKr58bK3itZ2m0pxrrGLDy7kLcOv0VyXrLiDiIrg6gjYNsczGwUt18tVWlC265dO5o2bcr27dt56623iIuLIzg4mNGjR2NkJI4y1Fn6JvD6n9DubYg/I4w2hqx69mSxtCgheTs8H6ybwMQgYVWyF5GeoVBf+26IkBxFHIDAdnBkIRQrsJWOrAyOfSV0FdDVh3G7a9aySkEiN20ictOmx69sN0UoOUm7Laxql3S1+jvOS4d9HwkT74pyYehq6DFXeR0z7PyEvro3VVt2EBwnlBt0tu+ssmNquouxGejrSMX+szU0tZsrAL+fjlVvIMry88/CpqEkEgl9XPqw8+WdDHEfwpmkMwzZNYQ9d/Yo5gD3rgqTcMVygyqpNKH19/fn3Llz3Lx5k9mzZ2Nvb6+KuESqoKMH/RYLq0HpGgpJxfe+EPQJ3NgG4fvgwkr4ayQE+gu1PC3HwuSj/7ZwepFZOAgts97YB/U9hJriH9sIz11tTz3dDxNKDE58DbbNYNIRcGqnmLgVrf3bQhKany7EfPZnIRmvTNJV2Psh/OAnfJlyag9v/QPNhik3XolE6Azx4D7cU01vyeS8ZK6kXKGbQzex3OAhmUzOxbhMmjlYYKgnTnapiVZOlvg51mPH5UQyX5TVwzSQmb4ZCzstJLBXIMa6xvzfqf9j/pn5FJUV1W7HUULPajGhrZpKE9qffvoJf39/VcQiUpemg+Ddi9DzU5DqwdmfhHZcm0bB/o+FFUpce8CbB+Hln4TVo0T/cuksLKPbfwmU5AvP3aqXhFKB6ia2mbGwe7qwFG/SFejwLkw8pPmj4c2GCT18LRyEyYa/dILrW58esc7PgJDVsKIcu1yWAAAgAElEQVQr/NZN6Alr6QLD1sCbB4R+t6rwaPnI2wdVcrjDcUIdXB+XPio5Xl0QmfKA3MJS2jhbqjuUOksikTChkwtFpTI2XohXdzgvvK4OXdk6eCsd7DqwPXI7Y/aPISGnFp1gog6Dvhk4tlVckFqslstQibSGgSl0nQmdP4LUcEi5CbJSYQKUY1vhUvRsOrrCRC2fV+H4/+DyOtg4HOq7C621vAYJSWlFp9Fz7gmj3zd3QWQwIBcm7fX5om71HbRvBVP/EeqszwbCjkmgYyA8BwamwghuepRwX31TaDVe6NTQqKXqF+RwePiajjgIPT9R+uGCYoPEcoMnXE3IBIQZ+6Ka69/Mjv/tD+ePs7FM7uKKvq642po6WRla8ctLv/Db9d/45dovjNw/kh96/EBrm2quqFqQKUxcbdJfOJsqqpSY0IoeJ5WCjbewiarPpD4MWAKdPxS6SFxZL9QdH54PJg2ggZdwiVxI8NKiIPdh2yuJFDwDoN1b4Nq9bq66ZmAKvT4TyhBCdwgJeloE5KeBgZnQB9m1u7AynjpH+nV0wb033NgC2XfBQnmlVMl5yVxNvUpfl75iucF/XE3IAqClk4p6cWspPR0p4zo6883B2+y/cY9XWoplgeqmI9Xh7RZv42Ptw8wTM5kcPJlFnRYxwHVA1XcSfRzkMrHcoBrEhFYkUgYLe+izSCjjiD4BkUFw95Kw6lXhw6UT9c3AqjF4vCT0uPUIEBJiDWXbsWPV72xiLUwYazdFeQHVVpO+QkIbGQRtJijtMGJ3g4pdTcjG1twQG/MXoO2fko30d2L5kUh+Px3Dyy0aIamLX4Yr0qdul+h0dejKun7rmHZkGnP+mUPSgyQmN69iy8U7D9uhiu26qqzKCW1WVhb16lX8TToqKgp3d3eFBSUSaQ1dA/DsI2yPlJUIo7F1bNWX1nPmqDsExXLrKfx/uHNMqQmtWG7wtPziUm4n59DbW2xFpAiWJvq82sqBjefjuZaYTQtHLRn1/u47dUdQa15WXmzsv5FpR6ax/MpyisqKmNZiWuVfOqKPg2VjsHRWSZzaoMrFNj4+Puzfv/+p67///ntatGih0KBEIq2mo1fnklmtZGQp1O/GnKxaV4YaSM5L5lrqNbo7dMdQVxyJfCT0bg4yObRwFOtnFWVUWycAtlxU0nLUohqzMbFhdcBqfOv7suL6Cn688uPzF2HIiIGsOHDtprogtUCVE9qZM2cyfPhwJk2aRG5uLjExMXTr1o0lS5awZcsWZcYoEok0wKXFi7m0WMuWOXbtDoVZQr9HJQiKDQLE7gZPejQhzM9R7D+rKD6NzGlqZ86eq0kUFCvnC5rKffihsGkBCwMLVvRZQXPr5qy8sZIfLv/w7DvHnBAuXburIjStUeWE9oMPPuDy5cuEhYXh7e2Nn58fTk5OhIaG0r9/f2XGKBKJNEDymTMknzmj7jAUy7W7cBl9XCm73xe9DzM9M7o4dFHK/uuqawnZSCTQ3EFLTo1rAIlEwmttHMgtKuVAaBVWfqwLgoOFTUuY65uzovcK/Br4sTp0NWtD11Z8x+jjgARcxFVYq6Na/T3MzMwwNzcnLy+PwsJCmjZtirm5ubJiE4lEIuVyaAu6RkpJaKOzo7mVcYuXnF/CQMdA4fuvy64mZOHR0BRTA3FesiK90tIefR0pm0PEsgNNZapvSmCvQNzrubP00tKnVxWTyYQyKNtmGj1JWBNVOaHdsGEDPj4+WFhYEBERQXBwMKtWraJDhw7cvn1bmTGKRCKRcugZgnMHiD8HJQUK3fW+6H0A1WvV8wJIyS3kblaB9kxc0iD1jPXp42PD+ZgMYtMUuAy3SKEsDCz45aVfsDWxZd7peZy+e/rfG++HCi0dXburK7w6q8oJ7fvvv09gYCBbtmzB2tqa7t27c+PGDVq0aEGrVnWo+btIJBL9l0tnKCsW2qopiFwuZ1/0PhoaNaSNTRuF7VcbXE8Q2tb5iQmtUrzm7wjA1kviKK0mszWx5deXfsVYz5gPj39IRGaEcMOjs0XihLBqq3JCGxoayqhRox67zsTEhBUrVvD3338rPDCRSCRSCedOwmWc4uqDr6Ve4+6Du/Rr3A8dsaPFYx4tqOAn1s8qRSc3a+zrGbHtUiKlZTJ1hyN6Drd6bizvuZwSWQnvHX2PzMJMYUKYjr7Qm1xULVVOaO3s7Hjw4AEbN27kiy++ICtLeFO6desWfn5+SgtQJBJpBo/XX8fj9dfVHYbiNWoJuoYKTWj3Ru8FxHKDilxLzMJQT0oTWzN1h6KVpFIJw1o7cD+niNN30tUdTu28846wabHWNq35tN2n3H1wlxnHPqAk7gw4tgN9E3WHVudUuSL/xo0b9OnTB3Nzc6Kjoxk1ahT16tVjw4YNJCYmsm7dOmXGKRKJ1MzziTM0WkPXABz8IeGCsOhFLZXISgiODcbVwhUvKy8FBKg9ZDI5VxOy8G1kgZ5OteYki6rhlZb2/HAkkt1Xk+jm2UDd4dTctGnqjkAlhnoOJSIzgo3hG/na3IBPG4vlBjVRrRraSZMmcfv2bQwN/20QPnDgQI4fP66M2EQikUg1nDtCSR7cu17rXZ1NOktmUSYDXAdozxKkChKTnkduYak4IUzJGlub0NzBgqCwZApLtKQnrZab6T+TdgY2bDY3Y6ehWKZUE1VOaC9evMibb7751PV2dnbcv39foUGJRCLNc2LaNE5o64jJo3q1uNPPv18V7IraBUD/xmJ/7ifdSBQmhDVzEBdUULbBfo14UFTKsfAUdYdSc4MHC9sLQFeqy7f5UmxLy/gyciO3M8TuUdVV5YTWwsKC5OTkp66/fPky9vb2Cg1KJBJpngcJCTxI0NKZ045tQaIDCedrtZvMwkyOJhylnW07HMwcFBSc9gi9KyS0vvZiQqtsA5s3QiKB3deS1B1Kzd25I2wvguJ8LO9eYYl+Y8rkZcw4PoPc4lx1R1WnVDmhfeONN3j//fe5efMmEomE7Oxs9u3bxwcffMDkyZOVGaNIJBIpl74JNPSudeuuPXf2UCor5VWPVxUUmHYJS8rBRF+HxvXFCS/KZmthSLvGVhwJTyGnsPa14SIlSwyBsmL8Gvfm4zYfE58bz+dnPkcul6s7sjqjygntggUL6NevH/7+/jx48IA2bdowbNgwRowYwZw5c5QZo0gkEimfQ2vIvYdBUc1mhsvlcnZG7cRc35xezr0UHFzdJ5fLCUvKpqmdOVKpWFusCoP97CkulREcJpYFarzYU8KlSxdGeY2ij3MfDsUdYsOtDeqNqw6pckIrlUpZuHAhGRkZhIaGcvbsWVJSUliyZIky4xOJRCLVsBcWQKj3ILJGD7+edp2orCgGuA4Ql7qtQGJmATmFpfg0EpdLV5V+vrbo6UjqdtnBiyL2FOibgp0fEomEBR0X4GLuwrKLy7iaclXd0dUJ1e6bYmBggLe3N23btsXMTOwjKBKJtISDkNBaPIiq0cN3Ru4EYKjHUIWFpE3CkoT6WR+xflZlLE306erRgNNRaaQ9KFJ3OKJnKc6HuxeFyak6QjdVU31TlnZfiq5Ul49PfExGYYaag9R8z+1D26NHjyq3nTl69KhCAhKJRJqp5axZ6g5Buaw9Qd8Mi9zqj9Dml+RzIOYA3vW9aWLVRAnB1X1hSTkA4gitig1u0Ygj4SkcCE1mbHtndYdTPcuWqTsC1XhYP4tL58eu9rT05LMOn/HJqU+Y+89cfu71s7jy4HM8d4S2c+fOdOrUiU6dOtG6dWsuXrxIQUEBzZo1o1mzZhQWFnLx4kVat26tqnhFIpGaNOrcmUadO1d+x7pKqgP2LTHPi4ay0mo9NCg2iPzSfF51FyeDPUtYUg56OhI8Gopn9lSpp1dD9HWlHAy9p+5Qqi8gQNi03X/qZ5802G0wQz2GcibpDL9e/1XFgdUtzx2hXbRoUfm/33jjDWbNmsWnn3762H2+/PJLwsPDlROdSCQSqZJ9G3RjTkLKTbBrXqWHyOVy/gr/CyNdI/q59lNygHVX6N1sPG3M0NcVVwhTJTNDPbp6NOBo+H3SHxRR31Ss79Y4/6mfrcjcdnO5mX6TFddW0Ny6OV0cnk58RdWood2+fTuvvfbaU9ePGDGCv//+W6FBiUQizbNv8GD2aXuT84d1tNy9WOWHXEm5wq2MWwxyHYS5vng6vSKpuUWk5Bbh20isn1WHAc1tkckh+GYd63bg4yNs2qyC+tknGegYsKz7Mkz1TZnzzxzuPrir4iDrhiontPXr12fr1q1PXb9161asrKwUGpRIJBKpxcNOByRWvR/to7Y6o5uOVkZEWuHfCWFiwq8OvZraoKcjYf+NOlh2oO2eUT/7JAczBxZ3WUxOcQ4zjs+gqEyc5Pek55Yc/Ne3337L6NGj2bt3L23btkUikXDhwgVCQkLYsEHskyYSibSAmQ0FBg0wquIIbdKDJI7EH6FTo0641nNVcnB1lzghTL3MDfXo4tGAExGpZOYVY2mir+6QRI/E/iNcVlA/+6SuDl2Z2nwqK66v4H/n/8f8jvOVG1sdU+UR2uHDh3Pr1i26du1KXFwcMTExdO7cmVu3bjFixAhlxigSiUQqk2XqDqm3oTCn0vtuCt+ETC5jjPcYFURWd4UlZSORgJetmNCqSz9fW8pkcg7dqmNlB9qukvrZJ73t9zYdG3Vke+T28laBIkGVR2gB3NzcWLx4sbJiEYlEIrXLNvXALv0sJF0G1+7PvF9+ST7bIrfhYu5Cx0YdVRZfXRSWlIOrtQkmBtX6yBEpUG9vG3SlEg7cuMeINo7qDkcEQv1s4kXhfeYZ9bNP0pHqsLjLYkbsHcGX57+kaf2meFl5KTXMuqJa7y5paWlcuHCBlJQUZDLZY7dNmDBBoYGJRCKROmSbugv/ePRB8wx77uwhtziX91q+h1Qiztx/lpzCEuLS8xns10jdobzQ6hnr09HdmlNRaWQXlGBhpKfukESJF0BWUmn97JMsDS1Z1m0Z4w+O58NjH7Jp4CYsDMQJl1VOaDdv3sybb76JVCrF2tr6sQUXJBKJmNCKRFqu648/qjsElcgxcQWpLty9/Mz7lMpKWXdzHWb6Zgx20/LOD7V0S6yf1Rj9fW05GZHK4Zv3GdraQd3hVG7XLnVHoFyP+s82rn4brmYNmjGn7RwWnVvEJ6c+YXnP5S/8F+sq//Zz5sxh9uzZZGdnExsbS0xMTPkWHR2tzBhFIpEGMHN2xsy5bqw0VFpaikQi4fjx4wAcP34ciURCaWnlCybIdPShgRfcu/bM+xyIOUBCbgJjm47FWM9YUWFrpdCHCa2vuOSt2vXxsUVHKuFAXVlkwd1d2LRV7CnQNwPbqtXPPmm453AGuQ7iROIJVt9YreDg6p4qJ7Tp6emMHTsWHR1x2TWR6EWUGxdHblxcjR5779493nzzTezt7TE0NMTJyYnhw4eTkpKi4Cgr1rFjR+7du4eubhVPStn5QU4i5KU9dVOZrIzfrv+GiZ4Jo5qOUnCk2qe8ZZc4Qqt2Vib6tHe14mREGrmFJeoOp3JRUcKmjR7Vzzo/u/9sZSQSCZ91+AwPSw9+uvoTZ5POKjjIuqXKCe2oUaPYu3evMmMRiUQa7OT06ZycPr1Gjx06dCixsbFs3bqV8PBw1q9fj5OTE3l5eQqOsmL6+vrY2tpW/QGPZhxXMEp7KO4QsTmxjPIaJdatVcHNpBzs6xlRz1hsFaUJ+jezo7hMxtFw1XyZrJWXXxY2bfSofta5U612Y6RrxHfdv8NY15jZJ2eTnJesoADrniontBYWFnz++ecMHDiQuXPnMm/evMc2kUgkqkhWVhZnz55lyZIldOzYERcXF7p168bSpUtp3LgxAJcvX6Z79+4YGRnh4uLC559//lh5QGpqKqNHj6ZevXpYW1szevRo0tPTy2/Pzs5m6NChGBkZ4enpSXBw8GMxPFlyMH/+fDp37sxPP/2EnZ0d1tbWzJo1C7lcLjzArgVX7pXRqv94DA0N6dy5MytXrkQikbDi+gqMdI0Y6z1Wyc9c3VdYUkZkygO8xdFZjdHH2xapBHGRBXWrRf3sk5zNnfmi8xdkFmXy0fGPKCmrA6PvSlDlce4LFy7QokUL8vLyOHfu3GO3/XeCmEgkEv2XiYkJJiYm7Nq1Cz8/v6dO+6enp9O7d29mz57NqlWrSExMZPLkyRgbGzN79mwAhg0bhoODA//88w8SiYTZs2czZswYDhw4AMAHH3xAWFgYR48eBeC9996rNK7r16/j5ubG0aNHuX37NiNGjKBLly7o6+tTau3Fq1vyadvUmg1/X+bWrVvl+4zKiuJNnzexNLRU5NOklW4n51Imk4tL3mqQBmYGtG1sxfHbqeQVlYqt1NSllvWzT+rl1IsJvhP4PfR3vgn5hk/af6KQ/dYlVX4lHzt2TJlxiEQiLaWnp8fKlSuZOnUq3333HW3btqVXr16MHz8ee3t7AgMD6dGjB7NmzQLA3d2dBQsWMG/ePGbPns3Jkye5ffs2R44cKU+GV65cib29PYmJiZibm7Nhwwb27NlDhw4dAFi0aBH9+vV7bly6urqsWLECQ0NDmjZtSo8ePThx4gS9e/cm6NgpUvIlrHzZFHNvb7y9vQm5GMLXi7/GUMeQcT7jlPukaQlxhTDN1L+ZHeeiMzgansIgsZ2a6j2qn3XrUeP62YpMbzmd0LRQNt3ehF9DPwa6DlTYvuuCF7vHg0gkUomRI0eSlJTEn3/+ib+/P2vXrsXb25tr165x48YNdu/ejampafk2ceJEYmNjkclk3Lhxg9TUVOrVq1d+u6enJwDR0dFER0dTWlpK27Zty4/3338/i4eHB4aGhuU/29ralk9Si4yMxN2uHuYFCVCQCYDcQShHGNV0FNZG1gp7bpQi7gxsHgPfuMLqADj/GzzRO1wVyieE2YsJrSYJ8LFFIqHudDvQNjXsP1sZXakuX3f9moZGDVl4diGRmZEK3b+mq/SrwbhxVRuJ+OOPP2odjEgk0l6mpqYMHjyYwYMHs2jRIlq2bMnSpUt58OABr7/+eoW1+FKplAcPHuDu7s6+ffueut3e3p7w8HCg+qVPenqPN5aXSCSUlZUBIJfLkegZAbmQfIMCR3/2R+8HYFKzSdU6jspd3wo7Jgv/tvWFlJtw4BxEBsHQ1WBUT2WhhCblYGWij625YeV3FqmMjbkhbZwtORaeSn5xKcb6YtmBSsX8I1wqOKEFsDayZmn3pbx58E3eP/Y+G/tvpJ6h6v7m1anSV7HYpkskEgEM2L1bYfvS09PD1dWVvLw8/Pz8OHz4MO7P6Dfp5+dHfHw85ubmNGzY8Knb3dzc0NXV5cKFCwQEBAAQEhJSq/g8PT2JvJtOTpEe5veusSrrGvcj7gNgpm9Wq30rVUQw7JwKZnYw7m9o0ASKHsCB2XB1A2wZC6O3g67yOw6Ulsm4nZyDv4uVOM9CA/XztSMkNpPjt1Pp38xO3eFULCxM3REoR+wpMDBXWP3sk1o0bMHcdnNZdG4RM07MYMVLK5RyHE1TaUK7Zs0aVcQhEom01P379xk7diwTJ06kWbNm6OnpsXfvXvbv38+qVavo1asXK1asYPLkybz77rsYGhpy7do1IiIi+PTTT+nTpw/NmjXj1Vdf5X//+x/29vbcuXOHrVu38ttvv2Fubs6oUaP48MMPWbNmDXK5nM8++6xWMQcEBNCwYUOm7k3izYZBLIsJI/dsroKeESUpyoU974OBKbyxF+q7CdcbmMLLP4GeIYSsgn0zhJ+VLDotj8ISGT7ihDCN1K+ZLQv33uRAaLLmJrTaqDgP7l4Ct54KrZ990ogmI7iTdYeN4Rv58vyXtJe3V9qxNIVYQysSiaok6dQpkk6dqvbjzM3NadGiBf/73/9o3749rVu35o8//uDnn3/mjTfewNHRkZMnT5KQkECnTp3w9/dnyZIlODk5AULZwcGDB2nSpAmvvvoqPj4+TJ8+nXr1/j2N9v3339OkSRO6devG6NGja53Q6urqsn3HTsIzdek3az8pB1N4/6P3MTAwqNV+ler4YshNgl6f/5vMPiKRQN+vwf0luLIerm9Rejihd8UFFTSZnYURLZ3qcfTWfQpLytQdTsWCgoRNmyQop362IjP9Z9KxUUe2R27nn/x/lH48dRMLZ0QiUZVc+eYbABp1rt4bsZGREd988w3fPHx8RZo1a8bBgwefebuVlRWrVz97aUdLS0t27tz52HXlPWWB7t27P/bz/PnzmT9//mP3X7t2LQBBDz9AW7Vqxftfd2Np6V1ec38Vjunh5vZEoqgpMuPg/K/QqBW0fqPi++jowpAV8Esn2DsDHPzBqrHSQhI7HGi+/r52XInP4p/INHp726g7nKfNmCFcalPpwaP+sypIaHWlunzb7VvG7B/Dnuw9dI3tSoBLgNKPqy7iCK1IJBJV4Ps13/PFjUgaJuXhdEnO0qVLGTNmjLrDqtj5FSArhZ6fgvQ58x5MrGHIL1CcC4eUuyBOWFI2Jvo6uNQ3UepxRDXX11dYPe+AuMiC6pTXzzZXyeHM9c355aVfMJOaMfefuVy4d0Elx1UHjUxoCwsLeeWVV/D09KRFixb07duX2NhYAFJSUujbty8eHh74+vpy6j+nQPPz8xk5ciTu7u54enqyY8cONf0GIpGoLssozGDl+ZUkrE7k1LwYFn77KzNmzCjvlatRCrPh8h/Q0Fuoy6uMW0/wGgi3dkNC7SbPPYtcLicsKQfvRuZIpeKEME3laGVMM3sLDt26T1GphpYdaJNH9bNOHZRaP/ske1N7JllOwkDHgPeOvUd4RrjKjq1K1U5o09LSSEtLU0Ysj5kyZQq3b9/m6tWrDBw4kClTpgAwZ84c2rdvT2RkJGvWrGH06NHly1kuWbIEAwMDoqKiCAoK4p133iEzM1PpsYpEIu1RIi9hxvEZSDtI2XRsJUWfmBP14zA+++wzzez6cvkPYcS1wzShVrYqes0DiRQOz1dKSAkZBeQWlooTwuqAfs1syS0s5UxUeuV3FtVOwnmhflYBy91WVyO9RizvuZzismKmHprKnaw7Ko9B2aqU0KalpTFlyhSsrKywsbHBxsYGKysrpk6dqpTk1tDQkP79+5e3emnfvj3R0dEAbNmyhWnTpgHg7++PjY1N+Sjt5s2by29r3LgxXbt2ZdeuXQqPTyQSaaeSshLWZa3j0v1LjGk6hiG+48DcAe5dU3doFZPL4fJ6MLKCZsOr/rgGTcBvFMSdElYsUrDyBRXE+lmN189X6HAgLrKgAiqsn62Iv60/S7otIacohwlBE7Quqa00oU1OTqZNmzYEBQUxa9Ystm/fzvbt25k5cyYHDx6kbdu23L9/X6lBLl++nEGDBpGeno5MJqNBgwblt7m4uBAfHw9AfHw8zs7OFd4mEolqx9TREVNHR3WHoTS5xbm8e/RdwovCecX9FWb6zxRusPODlFtQUqDeACty7yqk3QbfoaBbzQ4MHd4RLs/9ovCw/p0QJo7QarrG1iY0tTMn+OZ9SspUv5rcc7m5CZu2UHH9bEV6OvVkSfd/k9qozCi1xaJoEvl/p/5WYPLkyURHR7N///6nWtYUFhYycOBAXF1d+e2335QS4FdffcWePXs4cuQIBQUFODk5kZeXV3778OHDGTRoEOPGjcPMzIzo6OjyhHfmzJmYmZk9tQLRsmXLWLZsWfnPWVlZbN++XSnxV6SwsPCxJTdfFPmyfDLKMsiT5ZEnyyNflk+eLI8ieRGllFIqL6VMXkYpQgmJ5D//ycpk6OroggSkD//TlegiRYqORAcddKp1+bzHG0oNMZQYoiPRwNPL1aQJr7VH/1+lEmn5/ztNa7SfXJrM+qz13C+9Txv9NoywHIFUInzfd03chkfCFs75fkW2WcWLP6hLk9i1uNzbzznfL8k286j249vcXIRlzk1OtgykyMCqVrH897X247UywjPlLO+qg45YQ/tMmvD3CbAvVsbuGBnv+0nxttLIqTWP0ZTnrTp0ygrpGfImafX8uOI1R+XHf/I5Cy0MZX3Wegwlhky0nIiTvpPKY6qJiRMnkpiYWOFtlVYl79+/n23btlXYf9HQ0JCFCxcybNiw2kdZgSVLlrBjxw4OHz6MsbExxsbGAKSmppYnrXFxceX9Kp2cnIiNjX3stv79+z+13xkzZjDjUTsQwMHBoXyFIVUICgpS6fFULb0gnbD0MMLSwojJjiE+N56E3ARyinNqt+MSxcRXVUa6Rpjpm2Gub465vjlm+mblm4WBBZYGllgZWmFlaIWloSWWhpbUM6iHrlRzuuEp7LUml8Pp0xAVBe7u0KnTY/Wacrmc6OxobqbfJDwjnNsZt7mff5/0wnRyix9fkEBPqlf+vNU3qi9cGtanoXFDPCw98LD0wMqwdslVVZWUlbDh1gYCrwZSKitlRusZ2CXa0bdv33/vFAFs3EJ7Z0Pw16C/27JSWPYuWLnRfmg16mf/y0UGf71Od9No6DmyVuH897X26cXDeDUyoH8/1dcK1iWa8lnglpLL7mUnSTV0ICCgmbrDqZSmPG/VcucoXCijof+rBHRUfexPPmcBBNAmoQ0fn/iY37J/45uu39DDqYfK41KkSj9509PTcXzOaUYnJycyMjIUGhQIo6h//fUXhw8ffqyB+vDhwwkMDGT+/PmEhISQnJxM54d9MR/dtnbtWmJiYjhx4gS//vqrwmMT/UsulxOVFcX5e+e5nHKZ0LRQ7uU9Xov1KFFxNHPEwdQBKyMr6hnUo55BPSwMLDDTM0NPRw99HX30pHroSfWQSqTI5DLkyJHL5QQfCqZnr57IkVMmL6NMVkaprJQSWcljl4/+/eT1T97+rOuKy4rJL80ntzi3fMspziEuJ46c4hxKZM/PqiVIMDcwfyzZtTGxwcHUAQczBxzNHLE3tcdQtw6NLsTFQUAAERkZoKODZ0YGNG7Mgz3bOaETy5mkM5xNOktqQVabGO4AACAASURBVGr5Q4x0jbA3tcfbyhtLQ0sMdAyE0XdZKfml+WQUZJBRmEF0djRFZUVPHbKhcUNa27SmrW1b2tm2w9FcsaUORWVF7L2zl99Dfyc+Nx5HM0e+6vwVLRq2IOjuE43c7R4uT6lpdbRxpyAvBfwn1SyZBfDoA2aN4Pom6D4XpLUfnUvJKSQ1t4heXk8vUyzSTO4NzfBoaEpQaDKLXvbVnFH1wEDh8uHcmDpNzfWzFenm2I3VAauZfnQ6Hxz/gDlt5zDSq3ZfbNWp0oTWycmJixcv4uDgUOHtISEh5SOkipKYmMhHH32Eq6srPXoI3xgMDAw4f/48X3/9NWPHjsXDwwN9fX3Wr1+Prq7wa8ycOZMJEybg7u6OVColMDAQKyvVjPS8SNIL0jmZeJJz985x/t550guF2bESJLjVc+Nlt5dpZt0MX2tfXOu5YqRrVKPjPDrlC6An0cNYz1gh8ddGUVkRucW5ZBZmklmYSUZRxr//LhSStEc/x+bEcjX1KjL503VpjmaOeFp6PrY5mDk89jtrBLkcAgLgzh0iH/6deyYnUxYVQWq31sz5yh0kEtzruRPgEkCLhi3wsvLC0cyxSr+LXC4vT3Dv5t0lMjOSiMwIQtNCORBzgAMxBwDh+erUqBOd7Tvjb+tf49dCcl4yf0f9zV/hf5FRmIGZnhnvt3qfsd5jMdB5Rg2qmS2Y2mheQhu+T7j0HlzzfUh1oPlwOP0DxJ8Fl061DktcUKFu6tfMjuVHIrkQk0EHt/rqDkfw88/CpTYktDH/gIGFWutnK9K8QXM29NvA20fe5qvzX3Er/Rb/1+7/6tagy0OVJrSjR49mxowZNGnShKZNmz52282bN/noo48YP368QoNycHDgWaW9NjY2BAcHV3ibiYkJmzdvVmgsIsHdB3c5EneEI/FHHkvSPCw96O/an/Z27WnVsBWm+qZqjlS5DHQMMDAywNrIukr3L5OVkVqQSkJuAom5iSTkJhCXE0dUVhTHEo5xJP5I+X2NdI3wsvLCp74Pvta++Fr74mTmpN5609OnkcfEIHnYGu8RnTI5Dqkl/KD/Gj4vT8HGpGarDEkkEkz0TDDRM8HR3JH2dv+uN55RmMHF5Iucu3eO03dPs+n2Jjbd3oSeVI/WNq3p1KgTXvW9cDF3oaFxw6cS6JKyElILUrmVfotLKZe4fP8yYenCikN2JnbM8p/Fqx6vYqJXhcb/dn4QfRxKi0FXv0a/q0LJ5UJCa+UKDbxqt6/mrwsJ7bW/FJTQCh0OvMUJYXVKP19blh+J5EDoPc1JaLVFYbbQf9az7/MXPlETR3NHNvTbwJxTc9gZtZNbGbdY1m2Zws+MKVulCe3cuXM5e/YszZs3p3fv3nh5CW+et27d4vDhw/Ts2ZO5c+cqPVCR6qUXpHMw9iC77+zmZvpNQEi6ejr2pKdTTzo26kh9I/GN73l0pDrYmthia2KLv63/Y7cVlhYSnR1NRGaEsGVEcDPjJldSrpTfx0zf7N8Et74vPtY+2BjbKD3JzSjM4HjCcfJ2BzJUUkpF46F6hsb0LHaEGiazlbEytKKPSx/6uPRBLpcTkxPD6bunOX33NCHJIZy7d678voY6hpjqm6In1UNXqktucS5ZRVmP7c/CwIK+Ln3p17gfXRy6oCfVq3owdn4QGQyp4WCnASMsSVcg5y50eLfm5QaP2HgLv1/Y39D/W9Cr2RmVR8KScpBIoKmdWe3iEqmUl60Zja1NOBiazPxBPuKCGIoUexrkZeDaXd2RPFM9w3r83OtnVlxfwS9Xf+G1va/xaftP6e/69DwkTVVpQquvr8+BAwdYv34927Zt4+DBg8jlctzd3Vm1ahVjx45FqoC6K5FmKCor4njCcfbc2cOpu6cok5dhomfCINdB9HLuRcdGHWtcQiB6nKGuId71vfGu711+nUwuIyE3gRtpNwhLCyM0LZSrKVcfS97qG9bHx9qnPMH1qe9T6y8Wcrmc2JxYTt89zZH4I1xOuYxMLqONfiGvP2sBoeJiYYKYCkgkElwtXHG1cGWs91gKSgu4mnKV6Oxo4nLiiM+JJ780n5IyoX66gVEDGho3xNrIGtd6rrRu2BoXC5eal3SU19Fe1YyE9lG5gddAxezPdxgc+gzuHAOv2n2AhSZl49bAFGN9zZkcKaqcRCKhn68tPx+/w+X4TNq4iOV6ChN9TLh00+xJV1KJlLf93sbP2o+5p+Yy+5/ZHE88ziftPsHCQPPPuFTpHUcikTBu3DjGjRun7HhEaiCXy7maepXdd3YTFBNEbkkuUomUjo06MthtMN0du4tJrIpIJVKczZ1xNndmoKuQrJTKSrmTdYfQtFChe0R6GGeSznAy8WT542xNbMtHchubN8bezB57U/sKj1FUVkRyXjLRWdHE5MQQlvb/7N13fFRl9vjxz0x6741UkpDQQid0AkhTqg0Bxb4WWFdl7btrXV3WdfHn10Vlde0iKlhoGhCkdwgJJJSQ3gshIb3MzO+PC0hnkkxyJ5Pz5nVfFzIz956ZhMmZ557nPMkcKDpwvhba3sqeMUFjGB86ntF3jMJmxVBIu6QBt7U1hIcr3Q5U4GDtwLAuwxjWZVj7nNDcJoadiAdHLwiONc3xek5XEtqjq1uV0FbUNpJTVsuMfl1ME5doVzfFBPDe5jTWHS6UhNaU0jeDayB4mVfbv6sZHjicH2b8wMs7X+bnjJ85WHSQvw39G3HBcWqHdk3XTWi//fZbZs6cia2tUjeWmZlJcHDw+SUga2pq+M9//mOea5yLa8o6k8XqtNWsSV9DXlUeANEe0UyLmMaU8ClG14mKtmWttSbaM5poz2hu5VYAGnQNnDh9guRSJcE9curIZTW5AFZY8cqyV3CwdqDJ0ERlQ+VlnRo0aIj2jGZy18kM9h/MsIBhF0+6io+HSZPwr6xUZsHb2irJbHx86y93dxRuwcpqXPmH1I4EKoug6LCyMpip6vE8wpTJKsfXga4RrJpRjnGBFJkQ1qH16uJKkIcDvxwp4G9Te6jfL3riRHXPbwoVeVB6Avrd1aHeLz3tPXln7Dv8cPIH3tz3Jn/c9EfGBY/judjnCHAOUDu8K7puQjtnzhwKCgrw9VVasPTp04dDhw4RHh4OQGVlJc8//7wktB1EWV0Zv2T8wtr0tSSVJgHg4+DDPT3vYVrENKI9o1WOUBjD1sr2/MSxc2qbajledpzsymzyKvPIrcolLTcNF08XaptqsdJYne+j62XvRbi7cgk/3D0cV9trJCChoXD0KAOv0YfW4mk00KUfZO1U+r9aqXg5/fzlyxtMe9we0+G3v0PmNogY16JDnJsQ1lsmhHVIGo2Gm2IC+O/WdBJzK+gX7H79B7Wlt99W9/ymcO7/a/gYNaNoEY1Gwy3dbmFk4Eje3Pcm8Znx7CrYxT297uHunnfjYmtedfLXfVe+tNvAdRYWE2aotLaUzTmb2Zi9kd35u2kyNOFg7cD0iOlMDZ9KrH8sVmY481I0j4O1A/18+9HPt9/5r8XHxzNpogmaeGs0MHKksnVWAf2U5uglx8C/9/Xv31ZOnh2FN3U9Xs+zCe3R1a1IaJUR2p4yQtth3djbn/9uTefnwwXqJ7SWIH2zsg8378v11+Lr6MtbcW9xc+TNLNq7iA8SP2DZ0WXc1/s+5nafaxYtNcHIGlrRsRgMBjIqMtiWt41N2ZtIKE7AgAFrjTVDugxhWvg0xgaPNZsfQtExHFi0CICBz7X/so1m4cI6WrUSWr1eSar9eiv9cU3JO0opPTj5q9IWrAUj8Mn5FQS6O+DuaAatzUSL9At2p4ubPeuOFPDcjd3VLTt48kll31FHavV6JaH16w3OHX+hkRGBI/hhxg+sSV/DB4kf8M7BdyisLuSvQ/+qdmiAJLQWo6i6iL2Fe9ldsJvd+bspri0GlFG78aHjGRcyjlGBozrETEVhngp37lQ7BHV1OTvyXXAI+t+pTgxFh6GmFPrNNf2xNRqIHA/7PoJTJ8G7W7Me3qAzcLK4ivE92qaNm2gfGo2GSb39+WRHJsn5Z+gdqOLvjKv0nO8wilOgugT63KF2JCZjrbVmZuRMpnSdwg8nf2BYQDtNzDWCUQntP//5T5yclObjDQ0NvP3223h4eABQXV3ddtGJy+j0OnKrckkvTyetIo3k0mQOlx6mqKbo/H0i3CKYEDaBYQHDGBIwpEOu+CGE2XEPBXt3dSeGnS83aFlJwHWdS2hP/trshDavCvQG1E2AhEncFBPAJzsy+flIgXw/W+N8/ax5t+tqCRsrG2ZFz1I7jItcN6EdPXo0Bw8ePP/v4cOHc+TIkcvuI4xzpPQIqfWpuOa7YjAYOPdHb9Cj0+uoaaqhpqmG2sZaappqOF13mpLaEopriimpKaG4tpgm/e+rNmk1WiLcIxjeZTiD/AcxNGAovo4d/9KGEGbn3MSw7D2g16mz4k/aJrB2gJA2GhUJGwVaGyWhHfposx6aXaXMr5AOBx3fwBAPfF3s+PlwIU9NjFa/20FHlb4ZrGwh1HxGMS3ZdRPazZs3t0MYnccru17h2OljsMH4x2jQ4OXghY+DD908uhHqGnp+dnq0R7TUwgrRXs4tgVt6Anx7XPfuJlVfBdm7ldnSNm101cXOWfnlm7kdGmubtWpYduW5hFZG9Do6rVbD5N7+fL4ri+NFlXT3N+JDSl0FnM4EG0elzV1b/Yx2FI21ygphwUPA1ojltUWrSQ1tO3sw5kF2JuykV89eaDQatGjRaDRo0GCltcLR2lHZbBxxsHbAw94DT3tPrLXyrRJCdQFn62jzD7V/Qpu5HfSNbVducE7keMjYClk7lL8bKbPSgJ+rHf5unTyRsRA3xQTw+a4sVifmXz2hbaqHhC/g4BdKbfk5VnbQ+1YY+QT4dNJWkBnboKkWokzQZUYYRbKkdjYpbBIch0nR8kMuOpZus2erHYL6LpwY1m9O+5477Wz9bKSJ+89eKnI8bHhRqdc1MqGtaWgirwom9pI2T5YiNswTf1d7VicWXLns4MR6WPcUlGeBrYuyfLJPNDTWKP2aE5dB8vcw6Q0YdH/L+lbPn2+aJ6OGE78o+6jJ6sbRiUhCK4QwStTcNphZ39F4dAU7N3WWwD23fKZ3VNuex7cnuAQodbT8w6iHHMk7gwHoK31LLYZWq2FqnwA+2p7BoZxy+ocoE8FpqocNL8Ge98HOFca/ArF/uPyyetYu+OEhWLsQqoph7PPND2LBgtY/ETUYDJC6Xnm/6CDL3VoCrdoBCCFEh6HRQEAfKEhSJoa1l8oipW636+i2X6FNo1FGgUtPwOksox6SmFMOQL8gSWgtybS+XQBYnVigfOFUGnw0XklmQ4bB/F1KWcGVakRDh8HDW5UynS2LYPcH7Ri5yopToCJHGZ2VCXXtRhJaIYRRtixYwJaOOmJiSl36QWO10qu1vWTtUPahI9rnfOdKDc6VOVzHoZxyNEBMkEwIsyR9gtwI9XJkTVI++kPLYeloKDwMcc/CPWvALejaB3DwgLu+V64qxL8AuQeaF8D06crW0Rxfp+yjJqobRycjCa0QwihVOTlU5eSoHYb6LpwY1l7OJbRh7ZTQho8Bjfb3vrfXcSinHH9HcLG3adOwRPvSaDTc0cOel+r+ifbHh8HOBe5ZBWNfACsjKxadvOD2T0FrDSsfgPpK4wNIS1O2jib5J6VnddgotSPpVCShFUKI5jiX0LZnHW3mDnDpotTktQcHD+gyADK3Xbe0oqSynrzyWsJc5dKqRTEY4PAKHj4yhylWezniGgePbFfKXprLrxdMeBVOZ8Cm100fqzk5laas6Nd9KljJB7z2JAmtEEI0h2e4Mqu7oJ1GaKtPQclRZXS2PevxwuOU3qLXeZ5JuUr9bFdJaC1HyXH48lZY+QBWWited3yWu6r+SIOdZ8uPGfsQBMXC3qXqTKpsLyk/KvteM9WNoxOShFYIIZpDq1UWWChIAr2+7c/X3vWz53SNU/bpW655t0M5ktBajJoy+OV5eH+4Uj/ddy4s2IP74FmU1zSy/WRJy4+t1cLUxcrf1/5ZGQG2RMk/gr3b7/9/RLuRhFYIIZqrSz9oqISy9LY/1/n62ZFtf64LBQ9RGuRnXD+htbPWEiiLIXVchYdh3TPwdm/Y/Z7yge3BjXDz++DkzbQ+SreDVYfyW3ce/xgY/CDk7oOjq00QuJkpS4fCJKXcwNpW7Wg6HelDK4QwSv9nnlE7BPMR0FfZFxwC7zbuM5m5A5x827+fpY09hAxRltttqgdru8vuYjAYSMwpp3egG1baZkz2EcapKVNaQBUlQ8kxpY2aQacsL+saCK5dwD0E3EPBIxScfK5elmIwQF05VJVAdbGSfBUkQtpvUHZ24pVfjNKGq9ctyojqWSFejvQLdmd9ShG1DTocbK1a/pxGPw0JX8HGVyH6pmtPLlu8uOXnUUPy2XKDnlJuoAZJaIUQRukysp1HCM3Z+U4HCRBzW9udp/Y0FB1R6vHU6GfZNU5ZBjdnL3S9fMZ25qkaztQ10TfIHZCE1iQq8uDAJ3Dk+98TzXOsHZSRv4Zq0Ddd/lgrW6UTga0z2DhAY62ycldDjbLnCpf5XYOUlbz6zIbg2Kv+nE3v24VX16Sw6VgxU/oEtPz5OfvCsAWw9U1IWg7977r6fSd1sBU1U35UFl4JH6N2JJ2SJLRCCNFcXpFK0tDWk1uydgGG9q+fPSd8DGx6TSk7uEJCeyjnNAB9g92gSFq6tUpjHWz+B+z5AJrqwNkf+s5ROgT49QKfHuDiryScep2y+taZPCjPVpafPZ0FlQVQXwX1Z5RjOLgrI7m2jsqorr2bklA6+So9ZP17K6O7RnxYmtIngNfWprAqMa91CS3A8Mdgz1LY/rbyHLWtGPE1F2UZyvtB3zlSbqASSWiFEEZZe7bB+ZRVq1SOxAxoteDfR/kFptdfdHnWpNSqnz0noJ+yvGn6Fhj318tuTsypAKB/sAdHi9o7OAtSlg7f3q3UsvrFKK91twlXT/S0VuAaoGxBg9olRD9Xe4aFe/HbsRJOVzfg4dSKpM3eFWIfhG3/Vmppr9YRoFcvZZ+c3PJztZcUKTdQm0wKE0KIlggcoIyEXXpZ2JQyt4OjF/h0b7tzXIuVtZJM5x2AujOX3bw/qwwvJ1uCPR1UCM5CFB+DjydD4RGIew4e3gLRk81y1PLWAUE06PSsTmrl5DCAIY+CtT1sX9zxOx4YDHDoa6V/c8RYtaPptCShFUKIlggcoOzzmrmcp7HqKpQZ06HD1V0PvmucMhEpa+dFX66qbyIl/wyDwjzQyHr1LeJYWwif3qRM/pr1OYx93iwT2XNujPHHydaKFQdyW38wZx/oP0+5ypG9q/XHU1PeASg9DjGzrjh5UrQPSWiFEKIlAgcq+7ZKaLN2gkEPYS1YmcmUws/208zYetGXD2adRm+AwWGtaLbfmdWW0//YImXi36zPoOd0tSO6Lkdba6b0CSApt4LjhSaYBDjkYWW/54PWH0tNCV8q+/53qhtHJycJrRBCtIR7qFIO0FYJbcY2ZX+FyVjtyqe7Monokn60+zLLAIjtKgltsxkM8P1DONflw6Q3oPsUtSMy2m0DgwFYedAEo7Te3SBiHBxdAxUmOJ4aGmuVjhR+Mb+38xOqkIRWCCFaQqNRRmkLDyt9Wk0tcys4eqtXP3uORgNdRyvtw6p+Xylqb0YZjrZW9AxwVTG4DmrfR5AaT573aBjyiNrRNMvgMA9CvRz5/mAeTToTrJQ35BGlpGX/J60/lhqOroH6ChmdNQOS0AohjDL63XcZ/e67aodhXgIHgq5BmdBjSjVlyjHDRqpbP3vOubKDTKXsoKFJz6GccgaEeGBtJb9GmqU0Fdb/FTzCONr1fvP4/jaDRqPhtgFBlFbVszW1FUvhnhM5AdyC4dBXSjuyC/30k7KZs0NfgtZGqZ8VqpJ3IiGEUVxCQ3EJDVU7DPNyvo52v2mPm7UTMKhfbnDOuXXp05Wyg8N55dQ36aV+trkMBlj3lNIjduYH6Kwd1Y6oRW4ZGIRGg2kmh2m10O9OpYdu2qaLb4uMVDZzVZ6j/J+IngxOXmpH0+lJQiuEMEplVhaVWVlqh2FeziW0OXtNe9zMs/Wzak8IO8cjFDzCztfR7jh5CoBhEfJLvFlSfoT0zcrqWKHD1I6mxQLdHRge4cWGlCJKKk1QbtP/TkADBz+/+OsnTyqbuUr8GjBAv2usdibajSS0QgijbH3sMbY+9pjaYZgXR0/wjoKcPaY9bsY2cPZTJs2Yi65xcDoTTmex42QpDjZW9At2VzuqjqO+En55AezdYfwrakfTancOCaVRZ+CbfdmtP5h7iFLWcvxnqC79/eszZiibOdI1wYHPlBXdIserHY1AElohhGidkKFQkQMVeaY5XnUpFCdD2Cjzqq88W0fbcHIzCdnlxHb1xNZafoUYbcubUJkPN7wITt5qR9NqE3r64edqx1d7sk0zOaz/PNA3QtI3rT9WezjxM5zJhUH3KQuQCNXJu5EQQrRG8FBln7PbNMfL3K7szaV+9pyz5Q/lyb/SoNMzIlLKDYxWchx2vwdd+sPAe9WOxiRsrLTMjQ2loKKOX02x7nH3qcro9cEvOsbKYXs/BK21xXw/LYEktEII0RohZxPabBOVHZyvnzWzhNbZB3x74Zi3AzAwPKLjjzK2m02vgb4JbvyXWa8E1lxzYoOxsdLwv+0ZrT+YjT30mQUlRyHvYOuP15ZKjiv15D2mg4u/2tGIsyShFUKI1vAMV/rFmmqENmMbuAYqxzU34XE4N55ioEOR9J81Vn4CHF2tjEAGD1Y7GpPydbXn5v6B7Ms8zYGsstYfsP/ZyVUJX7T+WG1p30fKPvYP6sYhLiIJrRBCtIZGo4zSFh5WJv60xpl8ZU14c6ufPeu0vzIzf65PBlqt+cVnlja9Dmhg7AtqR9ImHhqtfPBauiW99QcL6KusuJX8PTTWtf54baG+Eg59DX69IaTjdqqwRJLQCiGMMmXVKqasWqV2GOYpbCQY9Gf7x7bCyV+VvZnOmt5YHUmTQctwrYkXkrBU+Yfg5AbodTP49VI7mjYR6evC+B5+bDhaxLHCM60/YN/ZUFehTLpKTlY2c5K4HBoqYfCDZvmhszOThFYIIVrrkoUHWix1AwY0rKruwVvxx/l0RwaHcsoxmMkkmfi0GhINEfiX7VfaFolr2/GOsh/5hLpxtLHHb+iGwQCL159o/cFibgeNlZI4mhuDQSk3sHNT6n2FWZFeE0IIo+RvV2bfdxk5UuVIzJBvD3DyVZrmt1BtbR2aExs5qo/kTz9dvIBFpK8zT46P4qYYfzQqjQrVNerYcbKUyS6DGFj9NRQkQtBAVWLpEMrSlYUUIsYpl9ItWEyQG5N7+fNLciGJOeX0bU1/Yhc/iLwBUjfAj9+AgztMmmS6YFsjYwuUHIOh88HWSe1oxCVkhFYIYZSEN98k4c031Q7DPGk0Sp/W4mSoKm72w6vrm1j04afY66pJdorlndn9iH9iNF8+MIR7h4dRVFHHgmUHmbV0F4k55W3wBK5vd/opahp0WEeOUb6QsVmVODqMnf9RylBGWPbo7DkLJ0ah0cDr6462/opC39lg0MHCJ2HhQtMEaAp7/gtolHIDYXYkoRVCCFMIH6PsM7Y262ENTXoe/Gw/AcXKCPjsufczo18g0f4ujOzmzcvTe7H56THcNTSEA1mnmbFkBy/8cJgzdY2mjf861h0uAKDP0PFgbd/68gpLVlUCh75S+s52NZPli9tYlJ8Lc2JD2JtRxqrE/NYdLPomsHNVamnNxekspa632wTwilA7GnEFktAKIYQpnKujTdvUrIe9s/EEu9JLmeV4AINLF6wDB1x2Hy9nO/4+M4ZfnhjN0HBPlu3JZuLirfyaYoKG9kZo1OmJT1ZadXUN8Fa6OuTsMd+Z6Grb8wE01Smjs51o4tDTE6Nxd7ThjXVHqWzNBy4bB+g1U3kNm+pNF2Br7PtIGXGPfVjtSMRVSEIrhBCm4B4Mvr2U9eiNnDC1P7OM9zenMcPvFJ4N+Wh6Tgft1d+Wo/xc+PoPQ/nHLTFU1zfx4Of7eXZFEnWNOlM9iyvacbKUitpGpvQJUL7QNU5JNnJMtJiEJWmoVpIfzwjoMU3taNqVh5Mtz07uTtGZel5eldK6g/Wdq+zNYZS2oQYOfq58TyPGqR2NuApJaIUQwlR6TIXaMsjedd276vQG/vLDEWystLwUcXZ2eM+Z132cRqNhTmwIv/45jlHdvPlmfw63vLeTrFPVrY3+qtYmKeUGU2LOJrThY5R9M0ejO4Wkb6CuHIY+alGrghlr9uBgxkb7sPJg7vkylRYJGQpaG6g/A/q2/cB2XUdWKN/T2Ieu+YFTqEu+M0IIozgHB+McHKx2GOat+1Rlf2ztde+68kAux4sqeXBkGJ6Z68DZH4KHGH0qP1d7Pr0vlifHR3G08AxT393eJiUINQ1N/HKkkN6BroR5n53ZHdAPnHwgdb3Jz9ehGQzKxCE7V+g7R+1oVKHRaPjnbX3wdLLlmRVJHC9s4WIjGg2EdwU3PaT/Ztogm+Pc99TWGfrNVS8OcV2S0AohjBK3ZAlxS5aoHYZ5848B9xA4tkb5RXgVtQ06/r3hOJ5OtsyPOqO0eLpOucGVWGk1PD6+G5/fH4u1VsMfvtjPf7emmbRv7ZrEAirrm7hjcMjvX9RqodtEKE6B8hyTnavDy9wGJUeh351g56x2NKrxdbHn/TsHUN+k4/5P91F0poW11mvWwhxHdXvSZu+CosPKBxR7We7ZnElCK4QQpqLRQI/pUJEDufuuerdv9mVTdKaex8ZF4nT4c+WL/ee1+LSjuvnw04KRRPo488a6YzyzIomGJn2Lj3ehr/Zm42BjxYx+XS6+odtEZZ8ab5LzWIQ9SwENxP5B7UhUNyTcizdujiGvvJY7lu4iv7y2+QfxDFeWlz26bXY7GQAAIABJREFUBupMsApZS+xZquxjH1Ln/MJoktAKIYxyYtkyTixbpnYY5u/cZcn9H1/x5iadng+3ZeDlZMucPm5weCUEDoKAPq06bYiXI9/PH86YaB++O5DLXR/toay6oVXHTM6vIDGnnOl9u+Bqb3PxjRFjQWsNJ6TsAIDybDi+Tto6XeD2QcG8fnNvMk/VcPN7OziYfbp5B1iyBI77QFMtHFWW3dbpDZTVGUjOr+BwbgW5p2to1Jnmw9tlKvLg6GoIHws+UW1zDmEyslKYEMIoqcuVy35Rc6WO7Jr8ekHwUDjyPUx6Axw9L7p57eEC8sprWTghCvuU75Rf1oPuN8mpXext+N89g3l97VE+3pHBzCU7+N89g+jm59Ki4y3dkg7AnUNDLr/R3k0ZPcvYoswC7+ykrdMV3TkkFGc7a55dmcQdS3fxaFwEj46JxMHWiAlz772HTqeDubZkb/yIP20P43hhJQ06Pezafv5udtZaBoR4MLm3PzP7B+LmYHONgzbD/o+VBR6GyPe0I5ARWiGEMLVB94OuHg5dPKJtMBj4aFsGDjZW3D3YD3a+Cw4e0Otmk53aSqvhxWk9eePmGPLLa7nlvZ1sPt781ctSiypZnZRPXJQPfYKuspRp1GSlfVfmtlZG3cGda+vkFSltna5gRr9Afpg/gnBvZ/5v00lG/+s3/r3+OIdzK2i6YHRVrzeQX17LusMFvLI6mayyGtJO1fJz4wC6ViVAeTbDI70YG6jhzxOieHpSNPeP6MrgME8O5ZTz0qpkhv1jI/+3MZWaBuNa511VYx0c+BTcQ38vrxFmTUZohRDC1HrOgPgXYOf/wcB7wE4ZIU3MreBwXgVzh4TgfuQzpdZ24utg62jyEOYOCSHMy5FHvzrI/Z/u48WpPblneBgaIxv9v7MxFYMBnhjf7ep3ipoE6/8CJ+LBZpKJIu+ADn8HtadhzPPS1ukqegS4svZPI/lydxYfbsvg3U0neXfTSay1GjydbAGoqG2k/oLa77t0BpztrHGInQcHdrNqdC6auLnEx8cz6YaLfy5rGppYk1TAB5vTWLzhBCsP5vL2Hf0YEOLRsoCTf4CaUhj5RKdsv9YRyf88IYQwNRt7GPdXqCqCrf86/+Uvd2cBcE8fZ9j6FriFtOkEouGR3vy4YARhXk68vDqFx75OoKLm+is4bTlRwpqkAsZ196X/tRICr0jw6Kq07zJhZ4UOxWCAvWfbOnXSVl3GsrbScu+Irmx9ZixfPjCEBWMjGBPtS6CHA8GejgyL8OLuYaEsuiWGjX+OI9zHiS5u9txw02xw8kVz6CvQX7le1tHWmlmDgol/cjR/uakHBRV13P7BLj7fldn8QA0G2LsUbByh/12tes6i/cgIrRBCtIUBdys1eLveg+5TKffqx+rEfIaEOBG9+RGlUfvUxWBt16ZhdPV24of5I3h2ZRJrkgpIyC7ng7sGEhPkdsX7l9c08MyKRJztrHl1Rq9rH1yjUUZp93yAc00nbd+VtROKjii1s9LWyShWWg0ju3kzspu3kQ+wVhLL7YshbeM172pjpeUPo8MZHeXDI18e4MWfkkkvqebFqT3Rao1chjhnD+QnwMD7lJIg0SHICK0Qwij+w4fjP3y42mF0HFormPr/wMoGPp/JkbXv019/hHf1b0DObhj2R+h9a7uE4uZow/t3DeCft8ZQUlXPrR/s5JMdGej1F4+qnq5u4N5P9lF0pp6XpvUkyMOIUogopdTAp/xAW4Ru/vZ8oOylVZfpTZyobKDUpWu0sPdDox4a7e/CD/OHMzTck093ZvLUd4kX1ete067/KPuh81sQtFCLjNAKIYwy8Lnn1A6h4wkaCPN+hK9uZ2Tyi4y0BUpRRpsmvNquoWg0Gu4YHEKPAFcWLDvIK6tT+CEhj3lDQ4nwdSa1qJL3N6eReaqG+WMiuG1gkHEHDh0Bdm74ndrbtk/AHJXnKKvCRY4H72vUGouWefvt3//uHgzRN8GxtTg4TzHq4e6Otnx6XywLvjrI9wl52NloeePmmGvXkZelK31vu02SVl0djIzQCiFEWwoZQu7cTTzV+DDx3nfDIztgxhLVJpr0CXLnl8dHc98IpQXS0yuSuOW9nTy78jAllfW8OLUnz0zubvTkMaztIPpG3KrT4HRmm8ZudvZ9dLat0yNqR9I5xD4EGAguNL73sb2NFR/MG8iobt58vTeH9zanXfsBe5YCBhgmo7MdjYzQCiGMcmDRIkBGaltiZaqBFbo4pk4YDP6+aoeDk501L03rxRM3RLE+pZDymkZcHay5KSYAl0sXUDBGr5mQtBxSfoIRj5s+YHPUUAMHPzvbqusGtaOxTE8+qezPjdR2HQ3e0QQV/6a8/kZ2B7Gx0vLenQO4/YNdvLX+ODGBboyO8rn8jrXlkPAl+PWGrnEmehKivcgIrRDCKIU7d1K4c6faYXRIqxLz8Ha2ZWSkkZNg2ombow23DwrmD6PDuWNwSMuSWYDwsTRZOSgJbWdxrlVX7MPSqqutrF+vbOdolGWFbXTVyuvfDC72NiydNxBnW2ue+OYQhRV1l9/p4OfQUAXDFijnEh2K/C8UQog2lFZSRVpJNRN7+WNtZaFvuTb2FHsMhLwDUJahdjRtz2BQLk3bukA/adXVrvrOVj487f1vs1vFhXo58a/b+1BW3cDTKxIxXPj4pgble+rs126TNYVpWei7qxBCmIcNKUUATOjpp3IkbavAe5Tyl6Rv1Q2kPWRuh+JkZXKfXcuWFRYtZOdCru9YpVVa6oZmP3xy7wBuHxjEttRSvt57Qau5xK/hTK6yzG0bt9ITbUMSWiGEaEMbUopwsrVieISX2qG0qVPufcDJV0kMLH2RhT0fABpp1aWSzIDpYGULW99s0c/a36b1pIubPa+vTaHoTB3oGmHbW2DvDoPle9pRSUIrhBBtpKSynoPZp4mL9sHO2rKXzzRorKDPLDidATkW3MLrdBYcXwfdJoJXhNrRdEr1dp7K6HjuPkj/rdmPd7W34bWZvalu0PHGuqPKh7DybBj+R1kcowOzuIQ2NTWV4cOHExUVRWxsLCkpKWqHJIRF6DZ7Nt1mz1Y7jA5l07EiDAbLLzc4r+/Zn49DX6kbR1va9xEY9MqladG25s9XtisZ8QRobeDXl6+6HO613NDDjxu6+xJ/KIP6X18HB09lgp/osCwuoX344Yd56KGHOHHiBM888wwPPPCA2iEJYRGi5s4lau5ctcPoUNYnF2Gl1TA2Wv1WXe3CPwYC+sHhFVBXoXY0ptdQrbTq8o6CiHFqR2P5FixQtivxCFU+VBQkQtI3LTr8S9N68YjNOuxqCtGNflZGZzs4i0poi4uLOXjwIHfddRcAt956KxkZGWRmZqobmBCi06lpaGL7yVJiwzxxd7RVO5z2M/hBaKyGxJYlGWbt0DIlUY99SNo6mYPRT4GDB2x8FerONPvhIdanmW+zhjR9AJ82jG2DAEV70hgMllO9f+DAAebNm3dRmUFsbCxvvfUWo0ePPv+1xYsXs3jx4vP/Li8vZ+XKle0WZ11dHfb29u12Pksgr1nLmPJ1O/O//wHgauFXPUz1miWU6PngiJ5ZkVpuCLaosYMrOve6aXX1jDnwCPW27uzou9hiEj+NQcfIhMex1tWyZcB76K1aPxNe3teurf9LLwGQ8MorF339wtctqGgDvdI/JMd3PCkRDxl/cIOBAccW4VOewJ94mvWN/XltqBWutpbx83opS/lZe+CBB8jNzb3ibRa3UtilyzVeKV9fuHAhCxcuPP/voKAgJk2a1OaxnRMfH9+u57ME8pq1jClft7VLlgBY/PfBVK/Z+m8TgVweuyWOYE/jVjTqyC563bS7sNn9HpO62VrOpfnDK2B3MYx5ngljppvkkPK+dh1nf09f+hpd9LrpJ8AXxwnO+JXgiY9C5Hjjjp24HHYnQJ87mNztAVZ9dZBkQzAvT+plymdgNjrDz5pFDRsEBweTm5tLU1MToCSzOTk5hISEqByZEKIz0ekN/Ha8mO7+Lp0imb3MkEdAYwXb31Y7EtMwGGDH/wMbR6XcQJgPrRam/wfsXGHFA1Caev3HFB6BNU8qbeYmL+LG3v70DXbnqz1ZZJ+qafuYRZuwqITW19eX/v378+WXXwKwcuVKwsLCCAsLUzcwIUSnkpxfQVl1A+O6d5LJYJfyCIWY2yBjK+QeUDua1kvbBIWHYcDd4OipdjTiUh6hcPsnUF8JX92utFa7mjP5sHyO0nv2ji/A0RONRsNzk7vTqDOweMPx9otbmJRFJbQAS5cuZenSpURFRbFo0SL+d7buTwgh2su21FIARnXzUTkSFY14Qtlv+7e6cZjCjneUEeehV2khJdQXOR6mLobTmfDhODi58fL7FB6Gj8YrPWenvg0hQ8/fNCzCi7goH348lE9yvgV26OgELK6GNjo6ml27dqkdhhCiE9t6ogRHWysGhLqrHYp6/HpCj2lwdDXk7IPgwWpH1DJ5ByFjC8TMUkYChfkaeC84esP3D8GXt0DIMAgfC3bOkJ8AR1YqH0xu+VBZBOQSz07uztbUEt785Tif3R/b/vGLVrG4hFYI0Tb6P/OM2iF0CFX1TRzIOs3oKMtfHey6xv0Njq1Vmt/fu6Zjdjz47Q1lP/IJdePojC7oRmS0HlNh/k7Y8i9IWg7ZFwxwdR0NN7wMQQOv+NCeXVyZ0bcLPx7KZ1faKYZZ+HLVlkYSWiGEUbqMHKl2CB3C7rRTNOkNjOrmrXYo6vOJhn53QsIXcOIXiL5R7YiaJ2snnNwAvW8DP8uc/W7WWjor3yMMZi6BKf+G4mRoagDXAOXr17FwQjRrkgr4V/wxVj46/LLOScJ8WVwNrRBCqGlbagnQyetnLzT2L2DrDD8/C421akdjPINBadivsYKxL6gdjWgJG3sIHAihw4xKZgFCvByZExvCwexyNh4tbtv4hElJQiuEMMra6dNZO900/Tct2bbUUrq42RPh46R2KObBNUBJCMuzYFsLLiGrJfkH5XL1wHvAK0LtaDqnXr2UrZ09Ni4Sexstb60/jl5vMWtPWTxJaIUQwkRyympIL61mVDcfuVR5odiHwa+30su19KTa0VxfQw2s/xvYuykjzKJT8XW1597hXTlWWMnqpHy1wxFGkoRWCCFMZPvJs+26oqR+9iJW1ko9o64B1v1ZuZxvzrb+C87kwpjnwUm+l53RI3HhuNhbs3jDCRp1erXDEUaQhFYIIUxkW2oJGg2MjJQk6DIhQ6H/XZC+GRK+VDuaqytIUvrO+sfA4AfVjkaoxN3RlodHh5N1qoZv9uWoHY4wgiS0QghhAjq9ge2ppfQJcsfd0VbtcMzT+FfBLRjWLoTsPWpHc7mmBvhpgfL3GUvAykbdeISq7hvRFW9nW/5vYyp1jTq1wxHXIQmtEEKYQFJuOWfqmhgt7bquzskLZi9TOgd8PVtZtMCcbHwFCpNg1EII6Kt2NEJlTnbW/HFsJMWV9Xy2M1PtcMR1SEIrhDDK6HffZfS776odhtmS5W6NFNAHZn0OjTXw2TRIWaV2RIoT8bDrPxAUC3HPqh2NAPjpJ2VT0ZwhIQS6O/D+ljTO1DWqGou4NklohRBGcQkNxSVUlv68mm2pJTjZWtE/pBMvd2usqIkw70flkv638+Dbu6H4mHrxFB+FFQ+AvTvc+pGUGpiLyEhlU5GdtRVPjO9GeU0jH21NVzUWcW2S0AohjFKZlUVlVpbaYZilyrpGDmaXMyzCGxsreVs1SugwmL8Huk+FlJ/gvSHwyRTY/T7kHoDGuvaJ40wBLLsDmmrhji/AQz60mY2TJ5VNZbcMCCLS15mPtmdQWlWvdjjiKmTpWyGEUbY+9hgAU1aZySViM7Iz7RQ6vYHR0q6reVz8YPZXkLMPdr2rXPbP2q7cprUB3x5K/1q/XuDfG4IGg60JF6yoLILPpiqLPkz/D3Qdbbpji9abMUPZJyerGoaVVsOfJ0Tx6FcH+b+Nqbw6o7eq8Ygrk4RWCCFaSZa7baXgwRD8OdRXKatz5ScoE8YKDkHist/vp7WBsBHQdw70mNa65LbkOHx1u5LM3vQWDJjX+uchLNbk3v4MCHHnqz3Z3DM8jAgfZ7VDEpeQhFYIIVppW2opQR4OhHk5qh1Kx2bnDN0mKNs5NWVQlKwkt+lblD626ZthzULoNVPpbRsyDJqzMtuRlbD6SWishmn/pyxvK8Q1aDQa/jKlJ7e+v5N/rDvGR/cMUjskcQlJaIUQohWyTlWTdaqGObEhstxtW3D0hK6jlG34Y0qCe2QlHFoGh75SNq9IJbHtOwdc/K9+rFNp8OtLcHQ1OHrDrBUQMbb9novo0AaGejClTwBrkwrYlXaKYRFeaockLiAJrRBCtMK5dl1xUj/bPhw9IfYPylZ8VFl1LPFr+PVl2PgahA6H8Dil9tbZV1lmtzQVjq2BE7+Avgl6TIcpi8FZSkRE8zw7qTsbkot4fV0KqxaMRKuVD7HmQhJaIYRohW2pJWg1MCxCEtp259sDJr0ON7wEJ35WRm3Tt0DmtivcWaOMxo5+Wkl6hWiBEC9H7hkeyofbMvgpMY+b+wepHZI4SxJaIYRRpLvB5Zp0enaePEW/YHfcHKR3qWqsbaHnDGVrrIOCRCg5BrWnAQN4dIWgQeAmyUeHonJ3g6v549hufHcgl3/9cpwbewdgb2OldkgCSWiFEKLFDuWUU1nfxEjpbmA+bOwhZIiyCdEG3Bxt+NO4bry6JoX/bc9gwVh1F38QCukALoQwSv727eRv3652GGZl6wmlXVdclCS0QphcfLyymaG7hoYS6uXI+5vTOF3doHY4AklohRBGSnjzTRLefFPtMMzKlhMluDnY0DfITe1QhLA8Cxcqmxmytdby1MRoquqbWCpL4poFSWiFEKIFyqobSMqrYGSkN9ay3K0Qnc6UmAC6+7vw2c5MSiplSVy1ybuwEEK0wLbUEgwGKTcQorPSajU8OSGK2kYd729OUzucTk8SWiGEaIGtJ5T+s6Ok/6wQndbEnn7EBLrx5Z4sCipq1Q6nU5OEVgghmslgMLA1tYRoPxcC3BzUDkcIoRKNRsPCiVE0NOlZ8ttJtcPp1CShFUIYxTk4GOfgYLXDMAtHCyopqaxntIzOCtF2IiKUzcyNifJhQIg73+zLIa9cRmnVIgmtEMIocUuWELdkidphmIUt59t1+aociRAWbNUqZTNzGo2GP93QjUadgY+3Z6gdTqclCa0QQjTT1hMl2NtoGRTmoXYoQggzEBflQ3d/F77em01FTaPa4XRKktAKIYxyYtkyTixbpnYYqquub2J/VhnDwr1kyUsh2tKSJcrWAWg0Gh6Ji6CmQceXe7LUDqdTkoRWCGGU1OXLSV2+XO0wVLfjZCmNOgOjpV2XEG3rvfeUrYOY0ieAQHcHPtmRQV2jTu1wOh1JaIUQohk2pBQBML6Hn8qRCCHMiY2VlgdGdqW0qoHvD+apHU6nIwmtEEIYSac3sOlYMd39XQj2dFQ7HCGEmbljcDBuDjZ8uC0dnd6gdjidiiS0QghhpITs05yqbmBCTxmdFUJczsnOmnlDQ8korWbTsWK1w+lUJKEVQggjSbmBEOJ67hoaipVWwzKZHNaurNUOQAjRMfgPH652CKrbkFKEn6sdMYFuaocihOWbOFHtCFrE382eG7r7suFoETllNVKe1E4koRVCGGXgc8+pHYKq0kqqSC+t5s4hIWi1GrXDEcLyvf222hG02J1DQ1mfUsTXe7N5ZnJ3tcPpFKTkQAghjHC+3EDqZ4UQ1zEq0ptgTwe+3Z9DQ5Ne7XA6BUlohRBGObBoEQcWLVI7DNVsSCnCydaK4RFeaociROfw5JPK1gFptRrmxoZSWtXA+pRCtcPpFCShFUIYpXDnTgp37lQ7DFWUVtVzMPs0cdE+2FnL6mBCtIv165Wtg5o1KAgbKw1f7pbJYe1BElohhLiOX1OKMBiku4EQwnheznbc2DuA3ellnCyuUjsciycJrRBCXMfqpHxsrbRSPyuEaJY5sSEAfHcgR+VILJ8ktEIIcQ3FlXXsSjvFmGgfXO1t1A5HCNGBDOnqSZCHAz8m5MnKYW1MElohhLiGnw8XojfAtL5d1A5FCNHBaLUabhkQRNGZerafLFU7HIsmfWiFEEbpNnu22iGoYlViPg42VtzQw1ftUIToXObPVzsCk7h1QCD/tzGVlQdyiYvyUTsciyUJrRDCKFFz56odQrtLL6niQNZpZvTrgqOtvF0K0a4WLFA7ApMI9XIiNsyT+ORCztQ1SulSG5GSAyGEuIoVB3IBmDUoWOVIhBAd2a0DA6lv0rM2qUDtUCyWJLRCCKNsWbCALRYyYmIMnd7AyoO5BLo7MCxcFlMQot1Nn65sFuCmmADsbbSsPPshWZieJLRCCKNU5eRQldN5Ws9sTS2h6Ew9tw0MQqvVqB2OEJ1PWpqyWQAXexsm9/Jnf9ZpMkqr1Q7HIklCK4QQV/DFriy0GrhtYJDaoQghLMCtZ99Lvj8oo7RtQRJaIYS4RGZpNb8dL2Z8Dz+CPR3VDkcIYQGGR3jj72rP9wfz0EtPWpOThFYIIS7x+a4sDAa4d0SY2qEIISyElVbDLQMCySuvZXfGKbXDsTiS0AohxAVqGg18tz+HaD8XmQwmhDCpc2UHK2RymMlJY0UhhFH6P/OM2iG0i425eirrDTwcF45GI5PBhFDN4sVqR2ByET7O9At25+fDhbw2owknO0nDTEVGaIUQRukyciRdRo5UO4w2VVHbyMZcA2FejkyXpW6FUNekScpmYW4dGERto45fjhSqHYpFkYRWCCHO+t+2dGqb4I/jumFtJW+PQgjTm9YnABsrDd8nSNmBKclYt2iWRp2e/PJaCivqaNQZsLXW4utiR5CHgyQAFm7t2QbnU1atUjmStpFTVsPSrekEOMKMfjI6K4TqevVS9snJ6sZhYu6OttzQ3Y/4lELyy2vp4u6gdkgWQRJacV3FZ+pYk1TA+pRCDuWUU9eov+w+9jZaendxY3SUD1P6BBDh46xCpEK03BvrjlLfpOeOXlps5MOZEKIN3TowiF+SC/nxUB7zx0SqHY5FkIRWXFVaSRX/2XSS1Yn5NOkNONhYMTDUg54BrgS4OWBno6W2QUfRmTqOFVaSkF3O/qzTLN5wgiFdPZk7JITJvf2xs7ZS+6kIcU3rDhfw85FCJvfyp4dnqdrhCCEsXFyUD55Otqw8kMujcREyAdUEJKEVlzlT18hb8cf5ak82Or2BoeGezBsaxrjuvjjYXj05bdTp2ZtRxjf7cvjlSCF7MsoIcLPnj+MiuX1gMLbWMuolzE9+eS3Pf38YTydbXp3Ri4RdW9QOSQhh4WyttUzv24VPd2aSlFtB32B3tUPq8CShFRfZfLyYZ1YkUVxZz8BQD164qTsDQz2NeqyNlZYRkd6MiPSmrLqBb/bl8OG2dP7ywxHe35zG4zd045YBQVhp1fskajAYqKxvoqyqgVPVDVTWNdLQpKe+SX/BXnfxv3V66ht1Z/d66s/uG3RXuG+TnvomHQ1Nepr0BpztrHG2s8bVwQZ/N3uC3B0I8nSkh78LPQJcpWWLymoamnj0q4NU1Dby0d2D8HW1VzskIUQnccuAQD7dmcn3B3MloTUB+W0qANDpDfy/X0/w7qaTuNhb889bY5g1KLjFl0E8nWx5dEwEdw8L5bNdmSzdks7TK5L4cFs6T0/qzvgevm16iaWksp4jeRWkFJwhs7SarFM1ZJfVUFbdQIPu8hrg5rLWarCz1mJrrcXO2ursXouLvTVe1rbYWWvRajRUN+iorm8i61QNB7JOX3QMjQYifZwZ3NWT4RFeDA33wtvZrtWxCeM0NOn547IEEnPKWTA2gvE9/dQOSQjRicQEutHN15lVifn8ZUpPuYrZSpLQCkoq63l8eQI7004RE+jGe3cOMNn69U521swfE8ldQ0P575Z0Ptqezh8+38+gUA+eu7E7g8KMG/29Fr3eQHL+GbamlnAop5wjeRUUVNRddB93RxtCPB3p1cUVTydbPJ1t8XKyxdXeBjubs0mplRY7Gy22VpcnqnYX/NvWWtuiUea6Rh0FFXVkl9WQnF/BkbwKEnMqWLYnm2V7sgHoH+LOtD5dmNInAD8zGy0c/e67aocAKKUtpVX1VNU14eVsh4ejTbM/HJ2pa+TRLw+w4+Qpbh0QxFMTo9soWiFEi/30k9oRtCmNRsMtA4L45y/H2Hy8mIm9/NUOqUOThLaTSy+p4u6P95J7upZ5Q0P569QebTKJy9XehqcmRXP3sFDe2ZjK8n053PbBLkZGevPgqK6M7uaDthlJ4qmqerallrLlRAlbT5RwqroBUNbK7ubrzMhIb2KC3OjVxY0IHyfcHW1N/pyay97Giq7eTnT1diIuyuf813PKatiVdoqtqSVsOlZMQnYKr61NIS7KhwdHhjMi0sssJgy4hIaqdu6M0mp+TMhj07FijhdWXjTK7uVky4hIb0Z182Z0lM91PwgczD7Nn79NJKO0mjuHhPDK9F5m8foKIS4Rafmz/2f278Kb8cdYeTBXEtpWMruE9oUXXuCHH37A1tYWOzs7Fi1axLhx4wDQ6/U8/vjjrFu3Do1Gw8KFC5k/f/75x/7973/nk08+AWDu3Lm89tprqjyHjuJQTjn3f7qPitpG3ry1D7MGB7f5OX1d7Xn95hgeGNmVdzamsiapgO0nSwl0d2BCTz+GRXjRM8CVLu4O50dBm3R6SmoNrErMJyH7NPszT3MkvwKDQTlmjwBXbh8UzOgobwaEeGBv07G6KgR7OhLs6ciswcHUNuj47XgxPybk8evRIjYfL6G7vwuPjolgWp8uzUr6Ta0yKwto38Q2p6yGN+OPsyYpH4MBPBxtGB7pRZiXE8521pRW1ZOUW8GqxHxWJeYDEO3nwqhu3gyP9CLCxxkXexuq65tIzj/DDwm5bEgpwtr8lk97AAAc9klEQVRKy4tTe3LfiDBJZoUwVydPKnsLTmwD3BwYEeHNpmPFnK5uwMNJ/cGXjsrsEtpRo0bxt7/9DQcHBxITExkzZgwFBQXY29vz5ZdfkpKSwokTJ6ioqGDAgAGMGzeO7t27s3XrVr7++muSkpKwtrZmxIgRjBw5kkkWuGyeKfx2vJj5Xx7EgIEP7x7IuO7tWz8Y7uPMO7P789TEaJbvy+aHg3l8ujOTT3dmAkp9qYudNQagqr5JSV53JwBKUjMlJoC4KB+jRuQ6EgdbK26KCeCmmAByymr4eEcG3+7L4fHlh/jf9gxeuKkHQ8O9VIlt62OPAe2zsILBYGDZ3mxeX3uUmgYdcVE+/GFUOEPDPa+4gEdJZT07TpayNbWEbamlfLQ9g4+2Z1zx2GOjfXjuxh5E+7u09dMQQrTGjBnK3sIWVrjULQMC2X6ylDVJ+cwbFqZ2OB2W2SW0N9544/m/x8TEoNPpKC0tJSgoiG+++YZHHnkEKysrPD09mTVrFsuXL+fll1/mm2++4d5778XJyQmA+++/n6+//loS2itYcSCXZ1cm4WJvzcf3DmZAiIdqsQR7OvL0pO48NTGatJIq9mWe5nhhJUVn6iivaUSrVcoV6suLmDq8D/1DPAjzcuwUo2rBno68NK0XfxrXjf/8dpLPd2Uy+7+7mdzLn1dm9LKoRP5C9U06/vbjEb7dn0uguwNL58UwqpvPNR/j42LHzP6BzOwfiMFg4HhRJfszT5N1qpqaBh221lqi/VwYFuFFqJdTOz0TIYS4vsm9/fnrj0dYeTBPEtpWMLuE9kKffPIJERERBAUFAZCdnU3oBZc7w8LC2L9///nb4uLiLrptxYoV7RuwmTMYDHywJZ1//nKMQHcHPrs/lkhf81jRS6PREOnrQqTvlUfN4uPjmTQgqJ2jMg8eTrb8bWpP7h4WyqKfj/HzkUJ2pJXy1yk9WtWJwhzVNuh46Iv9bEstZXSUD+/O6Y+bg02zjqHRaOju70p3f9c2ilIIIUzH0daaG3sHsPJgLmklVbLSZgu1e0I7atQojh49esXbEhISCA5W6jg3btzIK6+8woYNGy66z4W/vA3niiiNuO1CixcvZvHixef/XV5eTnx8vPFPopXq6ura9XwAeoOB707q2ZRrINAJ/tSzgbSEHaS1axQtp8ZrZo5m+kDXGC3Ljjfx7MrDfLrpCPf00OJhd+Wk1pSvW1VVFUCbfR8adAbeTdJxohxGBmiY06WM3Vs3tcm5rkV+1lpGXrfmk9fs2kacfc/ZcclrZImvW4hemej61oqtzAw3/TwQS3zNLtXuCe22bduue58tW7Zw3333sXr1aqKjf2+nExISQmZmJoMHDwYgKyuLkJCQi24758LbLrVw4UIWLlx4/t9BQUHtWpoQHx/fruerb9Lx528T2ZRbQGxXTz68e1CzR73U1t6vmTmbBDxc28gba4/yzf4c/pGg5R+3xHBTTMBl9zXl67Z2yRLl/G3wfWjU6XnkiwOcKC/m3uFhvDStp2ojz/Kz1jLyujWfvGbX4ayMVF76Glni6zZBb+CbjE0klsOSCeNMPgHYEl+zS5ldF9+tW7cyb948fvrpJ/r27XvRbbfffjtLly5Fp9NRVlbGN998wx133HH+ts8++4zq6mrq6+v5+OOPmT17thpPwayUVtVz54d7WJNUwI29/fn8/tgOl8yKy7k52PDP2/rw33kDsdJqmP/VQZ76LpGq+ia1Q2s2g8HA898fZuOxYm4dEKRqMiuEEGrQajXcPCCQ/Io6dqSVqh1Oh2R2NbQPPPAA9fX13Hfffee/9sUXXxATE8O8efPYt28fUVFRADz99NP06NEDgDFjxjBr1ixiYmIAmD17NpMnT27/J2BGjhac4cHP9pNXXsuDI7vy/E09VF12VpjexF7+9Atx5+nvklhxIJe9GWW8fUc/BoaafqJfW3U3WPTLMVYcyGVcd18W3RojyawQQmHh3Q0uNWtQMEt+S2P53pzrToQVlzO7hDY1NfWqt1lZWbHk7GXPK3nxxRd58cUX2yKsDmdVYj7PrUyiUafnzdv6MGtQ2/eYFerwdbHn0/sG8/muLN5Yd5RZS3fxx7GRPDbO/Hs3Lt2SxtIt6QwIcWfJ3AHYXKEllxBCdAahXk6MiPRifUohpVX1shR6M8lvDwtTXFnHn75O4E9fJ+Boa8WyPwyVZLYT0Gg03DM8jNWPjSTKz4V3NqZy+9JdlNRefXJkc+Vv307+9u0mO97yvdn84+djRPu58PG9g3Gw7VgLYggh2lh8vLJ1InNiQ2jUGVh5IFftUDocsxuhtXSnquqpajBQ09CEvbWVyQq/K2oa+XRnJh9uS6eqvonxPXxZdGsf+YTXyUT5ufDjguG8FX+cD7dlcDQPHENyubl/YKsv5Se8+SYAXUaObHWc6w4X8MIPhwnxdOSLB2LNYmliIYSZOTd5uxOVHkzo6Yenky3L9+Xw0OhwKcFqBklo29k9n+zlSJ6OP+9QPnXaWmmxt9Hi5miDp6MtHk62BHk4EOblRLiPE129nQnycLjipdjKukb2ZpTx85FC1h0uoKZBR6iXI4tn9WVCTz/5j9BJ2Vlb8ZcpPYmL8mXBF3tY+G0ivx0v4e8ze5vFhMCVB3J57vskvJ3t+PKBIfha6AIRQgjRXHbWVtw2MIj/bk1nd3oZwyLUWRmyI5KEtp1NiemCp6ESb/8u1DfqqWvUUduoo6K2kZLKeo4VVlLfpL/oMdZaDSGejvi42GFrraW2QUfhmTryyms51263b7A79wwLZXrfLldcGlR0PiO7efPiYCt+Oe3N6sR8DmadZvGsvgxRaencRp2eN9Yd5ZMdmQS42fPpfbGEeDmqEosQQpirOwYH89+t6Szfly0JbTNIQtvOHh0TQXj9SSZN6nfF2w0GAyWV9aSXVpNxdksvqSajtIrU4ioamvTY21jh7WzLTTEBDAr1YHSUj6wsIq7I2VbDf+cN5Ou9Oby6Jpk5H+7m3uFd+fPEKJzs2u+/f/GZOh5ffohd6aeIDfNkyZ0D8HGRchghhLhUhI8zQ7p68vPhQl6e1oCHk5RkGUMSWjOj0WjwdbXH19WeoSqNpAnLotFomDskhNiunjy9IpGPd2QQn1zIazN7Ma67X5ue22AwsPJgHq+uTuZMXRPzhoby4rSe0s1ACCGuYe6QEPZklLF8Xw6PjolQO5wOQX6rCNFJRPo6s+KR4bw2oxcVtY3c/+l+7vl4L0fyKox6vHNwMM7BxnfMyD5Vw/2f7uOp7xKxtbZi6byBvDaztySzQgjjREQoWyd0Y+8AfF3s+GJXJk06/XXvL2SEVohOxUqrYd6wMCb09OcfPx9lVWI+W06UMCUmgIfjwokJdLvqZMK4a/SAvtDJ4ir+tz2D7/bn0KQ3cMuAQF6c2lM6GQghmqeNFnPpCGyttcwbGsq/N5wgPrmIKX0uX9pcXEwSWiE6IX83e96Z3Z9H4iL49/rjrD1cwNrDBXT3d2HWoGCm9AnAz8juAwaDgbSSaranlrDucCF7M8sAGBDizjOTu0vpjBBCtMDcISG8+9tJPtmRIQmtESShFaIT6xHgykf3DOZIXgXL92Xz06F8Xl2TwqtrUgj3caJ3Fze6ejvh42KHdtMaDBhoGjOVojN1FFTUkXe6lmOFZzhT1wSAvY2WqX0CuGd4GINCPaR1nBCi5c5dFVqwQN04VOLlbMeMvl347kAuSbnl9AlyVzsksyYJrRCC3oFu/D0whv/f3v1H1ZjvewB/736XUnRIpfRr51epdCKJGiM5Y2zMOTRzzGRyzDU6GHMnY2ZYM5mlMWtdF5cxZzodCXViDIaOMcy9Sxq/pohkkLRLW34lJopS+t4/uu0rlF2qZz/2+7XWXtvzPHt/n3ef9cXH9n32s2TCIOw/ex2ZBTeQU3ILu/OuaF/z3oHtAID/uuOu3WdraYr+fWww2MkWIzztEer1uy799gQieoF9/XXjs4E2tAAQM9Id205cxobDJVgV9fRvR6JG/JuHiLQsTI2h8nOCys8JQOPNOy7fvo+bVbW4esoSCgDp7wSjl405nOwsYGXGP0KIiDrLIKfuCPboiX+dvoKP/zCAN6JpBS83JqIW2ViYYqBjd4xS9kI3M2NYmRljhKc9vHpbs5klIuoCM0e6o+6hQPLhEqmj6DU2tERERER6auxAByh7WyP12CVU3quTOg6AxouBT5bextXK+1JH0WJDS0RERKSnjIwUiH3JE1W19Ug5UiJplt/uPcDa/ynE6P84gClfH0H6L6WS5nkU/8+QiHTSJyRE6ghEZEjGjZM6gd6YOMQJq34qRPLhYrw90g22lqZdev6auof4x89qfJ1ZhHsPHsKhuzlmhbpjvI/+fJ0YG1oi0kngRx9JHYGIDMmqVVIn0BsmxkaY/7IScdvykJSlRlxk/y4795mySsxPPwn1zWp4/K4b3hurxARfR5jo2V0f9SsNERERET1hSoAzlL2tsf5QMW7cren08zU0CCRlqTHl68PQ3L6Hj/4wAD8uGI1J/s5618wCbGiJSEcnvvwSJ778UuoYRGQo3n+/8UEAGm9dHhfZH/frHmLl/gudeq6q2nr8ZWMOEn44B5eeVtgZOxLvhnnCzER/20YuOSAinVw7ckTqCERkSPbvlzqB3hk3yAEjPOyx9bgGrw9zhb9Lx989rPxuLWam5CC/rBKvDXXGssk+sviaRv1ttYmIiIhIS6FQYOmkwTBSKPDprjOof9jQoeNfqqjGn745gvyySrz3shL/OdVPFs0swIaWiIiISDa8HWzwl1B3nL5cicQsdYeNm3+5En/82xFobt1DwhQfvB/hDYVC0WHjdzY2tEREREQy8u8R3lD2tsbq/76AM2WVzz3ewQvliPr7Udypqcff3gzE9OH9OiBl12JDS0RERCQjFqbGWBXlDwB4N/UEKqpq2z1W2i+XMDMlByZGCqTNGo7IwX06KmaXksfCCCKSnPL116WOQESGJDZW6gR6zcfZFssm+2DR9nzMSc3FxpnDYGlmrPP7GxoElu89h6Sfi+Ha0wobYoLg2cu6ExN3Lja0RKQT7z//WeoIRGRI/vpXqRPovaggV1y4XoX1h4oRk5KNf8wIgrX5s1u78ru1WPhdHjILyvH7fj3w9+jfo2c3sy5I3Hm45ICIiIhIppZMGIi3Q9xwTH0LU9Ydxtkrd1p8bUODwLfHNRi/OguZBeWYGtgXqbOGy76ZBfgJLRHp6OD/fVoStm6dxEmIyCCoVI3Pu3dLm0PPKRQKfDZxEPrYWmDFvgKovjoElb8T/hTYF77OtjAzMcKVaoHkQ8X4Z3YpLt6ogn03M6x9IwAT/Zykjt9h2NASkU6qNBqpIxCRISkqkjqBbCgUCrwb5okQT3t88cM57Mgtw47cssdedRbdLUwwO8wDseFesLU0lSRrZ2FDS0RERPQCGNLXDlv+bQTOlFXi0MWbuHDtLhqEQHXFVfwxLABh3r3bdOGYnLChJSIiInqB+DjbwsfZVru9b98NRPo4Spio8/GiMCIiIiKSNTa0RERERCRrXHJARDoJ+PBDqSMQkSFZuVLqBCQjbGiJSCdOoaFSRyAiQxIZKXUCkhEuOSAiIiIiWWNDS0Q62aNSYU/TF50TEXW2wYMbH0Q6YENLRERERLLGhpaIiIiIZI0NLRERERHJGhtaIiIiIpI1NrREREREJGv8Hloi0snotWuljkBEhmTXLqkTkIywoSUindj06yd1BCIyJF5eUicgGeGSAyLSyd1Ll3D30iWpYxCRobh4sfFBpAN+QktEOsmaNw8AMGH3bomTEJFBmDSp8fnXX6XNQbLAT2iJiIiISNbY0BIRERGRrLGhJSIiIiJZY0NLRERERLLGhpaIiIiIZE0hhBBSh5Caubk5evXq1WXnq6qqgrW1dZed70XAmrUP69Z2rFn7sG5tx5q1D+vWdi9KzcrLy1FbW/vUY2xoJdC3b19cvnxZ6hiywpq1D+vWdqxZ+7BubceatQ/r1naGUDMuOSAiIiIiWWNDS0RERESyZhwfHx8vdQhDNGLECKkjyA5r1j6sW9uxZu3DurUda9Y+rFvbveg14xpaIiIiIpI1LjkgIiIiIlljQ0tEREREssaGtgsVFhYiJCQE3t7eGDZsGM6ePSt1JL3k5uaGAQMGwN/fH/7+/ti6dSsA4MaNGxg/fjyUSiV8fHxw6NAhiZNKZ/78+XBzc4NCocCZM2e0+1ur0b179/DGG2/Ay8sL3t7e2LFjhxTRJdVS3cLDw+Hh4aGdc6tWrdIeM/S61dTUYPLkyfD29oa/vz/Gjx+PkpISAJxvLWmtZpxrrRs3bhyGDBkCf39/jBo1CqdOnQLAudaalmpmcHNNUJd56aWXxIYNG4QQQmzbtk0EBwdLG0hP9evXT+Tn5z+xPyYmRnz22WdCCCGys7OFq6urqKur6+J0+uHgwYNCo9E8UavWarR06VIxY8YMIYQQarVaODg4iFu3bnV1dEm1VLewsDCRkZHx1PcYet3u378v9uzZIxoaGoQQQqxdu1ZEREQIITjfWtJazTjXWnf79m3tr3fu3CkCAgKEEJxrrWmpZoY219jQdpHr168LW1tb7W/AhoYG4eDgIIqLi6UNpodaami7desmbty4od0OCgoSBw4c6MJk+ufxWrVWo0GDBons7GztsalTp2r/gWVo2tLQsm7N5eTkCE9PTyEE55uuHq0Z55ruUlJSRGBgoBCCc01Xj9bM0OYalxx0EY1GAycnJ5iYmAAAFAoFXF1dUVpaKnEy/TR9+nT4+vpi1qxZKC8vR0VFBRoaGprdotjNzY31e8SzalRaWop+/fo99RgBCxcuhK+vL6KioqBWq7X7Wbfm1qxZg4kTJ3K+tUFTzZpwrrUuOjoaLi4uWLJkCTZu3Mi5poPHa9bEkOYaG9oupFAomm0LfmPaU2VlZSEvLw+5ubmwt7fHjBkzALB+unhWjR49zvr9v82bN+PcuXM4ffo0Ro0ahVdffbXZcdat0RdffIHCwkIkJCQA4HzTxeM141x7tk2bNkGj0WDZsmVYuHAhAM61Z3lazQxtrrGh7SIuLi64fPky6uvrATROHo1GA1dXV4mT6Z+mmpiammLBggX4+eefYW9vDwAoLy/Xvu7SpUus3yOeVSNXV1fthSmPHzN0Li4uABr/gJ87dy7UajUqKioAsG5NVqxYgR07dmDv3r2wsrLifNPB4zUDONfaYsaMGThw4IB2m3Pt2ZpqVlFRYXBzjQ1tF+nduzcCAgKQmpoKANi+fTvc3Nzg5uYmbTA9U11djd9++027nZ6ejoCAAADA1KlTsW7dOgBATk4Orl27htDQUEly6qvWavToseLiYhw8eBAqlUqyrPqivr4e169f125v374dDg4O2oaNdQNWrlyJ9PR0/PTTT7Czs9Pu53xr2dNqxrnWujt37uDKlSva7Z07d8Le3h49e/bkXGtBSzXr3r274c01aZbuGqbz58+L4OBgoVQqRWBgoDhz5ozUkfROUVGR8Pf3F76+vsLHx0eoVCrthXPXrl0TERERwsvLSwwaNEhkZmZKG1ZCsbGxwtnZWRgbGwsHBwftBSet1aiqqkpMmzZNeHp6CqVSKbZt2yZVfMk8rW5VVVUiMDBQ+Pj4iCFDhogxY8aIU6dOad9j6HXTaDQCgPDw8BB+fn7Cz89PDBs2TAjB+daSlmrGuda60tJSERQUpK3Pyy+/LE6ePCmE4FxrSUs1M8S5xlvfEhEREZGscckBEREREckaG1oiIiIikjU2tEREREQka2xoiYiIiEjW2NASERERkayxoSUieoFt2LABLi4uMDIyQkpKitRxiIg6BRtaIqLnFB4eDoVCAYVCAUtLS3h6euLtt99GXl6epLkePHiA2NhYLFq0CGVlZYiKimp2PDU1FXZ2ds2+mB0A9u3bB0tLSxQUFHRlXCKidmNDS0TUARYsWICrV6+ioKAA69evR11dHYKCgpCRkSFZpitXrqCmpgYTJkyAo6MjLC0tmx1/8803ERoaijlz5mj3VVdXY/bs2YiPj0f//v07NE9tbW2HjkdE1IQNLRFRB+jWrRv69OkDV1dXhIeHIy0tDdHR0ZgzZw7q6uoAALt370ZwcDBsbGzg5OSE2NhYVFdXAwCOHTsGCwsL3L59u9m4o0aNQnx8fIvn3bRpE7y8vGBubg5fX1/s3bsXAJCZmQl3d3cAgIeHBxQKRbN7tzdJTEzEwYMHsXXrVgDAxx9/jF69eiEuLg5qtRoTJ06EtbU1nJycMHfuXNy7d0/73uXLl2PgwIGwsrKCUqnEmjVrmo0dHh6OuLg4vPPOO+jevTs++OCDthWViEhHbGiJiDrJvHnzUFZWhtzcXABATU0NFi9ejLy8PGzZsgUHDhzA0qVLAQDBwcFwd3fHt99+q31/cXExDh8+jLfeeuup4x85cgQzZ87E/Pnzcfr0aUyZMgWTJ09GSUkJQkJCcPToUQBAdnY2rl69ChcXlyfGcHZ2xooVKzBv3jxkZGQgKSkJycnJePjwISIjI6FUKnHixAns2rULOTk5zZpSc3NzJCUl4ddff0VCQgI++eQT/PDDD83GT0xMhKenJ3Jzc9nQElHnkfreu0REchcWFiYWL178xP6amhoBQGzZsuWp70tPTxfu7u7a7eXLl4uQkBDt9tKlS8XIkSNbPG9UVJSYOnVqs33Dhw8XcXFxQgghCgsLBQBRXFz8zJ9h7NixwsjISMTHxwshhNi4caMIDAxs9prDhw8LMzMzUV9f/9QxZs+eLWJiYrTbYWFhIjw8/JnnJiJ6XvyEloiokwghAAAKhQIAcPbsWUyZMgWurq6wsbFBTEwMNBqN9vXR0dH45ZdfUFRUBKDxoq3o6OgWxy8oKEBwcHCzfSNGjGjXxVyLFi1CQ0MDFi1aBADIz89HXl4erK2ttY+IiAg8ePAAZWVlAIA9e/YgNDQUDg4OsLa2RnJycrOfBwACAgLanIWIqK3Y0BIRdZLz588DANzc3AAAKpUKCoUCaWlpOH78ONasWYP6+nrt652cnDB27Fhs3rwZR48ehUajwbRp01ocv6lh7ggmJibNnquqqjB69GicOnVK+8jLy0NhYSEcHR2hVqvx2muvYcyYMdizZw9OnjyJ6Oho7XrhJlZWVh2WkYioJSZSByAielGtXbsWLi4uGDp0KG7evImioiJ899138Pf3B4Bm62WbxMTEYPHixbh+/TpUKhXs7OxaHH/AgAE4duxYs31Hjx7F6NGjnzu7n58fdu/ejb59+8LCwuKJ47m5ubC0tMTnn3+u3VdcXPzc5yUiag9+QktE1AGqq6tx7do1lJaWIjMzE9OnT0dqaiq++eYbmJiYoEePHujRoweSkpKgVquxdetWJCYmPjHO5MmTUVFRgfXr17e63AAA5s+fjx07duCrr77ChQsX8Omnn+LkyZOIjY197p9n+vTpMDMzQ1RUFHJycnDx4kVkZGQgLi4OAODp6Yk7d+4gJSUFFy9exLJly5CTk/Pc5yUiag82tEREHWD16tVwdHSEt7c3Zs6cCVNTU+Tk5OCVV14BABgbGyMtLQ379+/H4MGDkZiY2OzTzSbm5uaIiopCjx49EBkZ2eo5Q0JCkJycjNWrV8PHxwc7d+7E999/r13i8DxsbGyQmZkJMzMzREREwM/PD0uWLIGjoyOAxrWxCQkJ+PDDDzF06FCUlJRg9uzZz31eIqL2UIiOXIRFRETPTaVSwcvLCytXrpQ6ChGRLHANLRGRnqisrERWVhZ+/PFH5OfnSx2HiEg22NASEemJSZMm4fjx451y21kiohcZlxwQERERkazxojAiIiIikjU2tEREREQka2xoiYiIiEjW2NASERERkayxoSUiIiIiWWNDS0RERESyxoaWiIiIiGTtfwFan3N4fRsQkQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 800x640 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a=ndvi.columns[10]\n",
    "plt.figure(num=None, figsize=(10,8), dpi=80, facecolor='w', edgecolor='k')\n",
    "ndwi.loc[:,a].plot(label='NDWI')\n",
    "psri.loc[:,a].plot(label='PSRI')\n",
    "ndvi.loc[:,a].plot(label='NDVI')\n",
    "\n",
    "plt.plot(validation.loc[a,:].seeding, np.array(ndvi.loc[ndvi.index == validation.loc[a,:].seeding,a])[0], marker='o', color='red')\n",
    "plt.text(validation.loc[a,:].seeding+8,np.array(ndvi.loc[ndvi.index == validation.loc[a,:].seeding,a])[0]+30,'Seeding',fontsize=12)\n",
    "\n",
    "plt.plot(validation.loc[a,:].harvest, np.array(ndvi.loc[ndvi.index == validation.loc[a,:].harvest,a])[0], marker='o', color='red')\n",
    "plt.text(validation.loc[a,:].harvest+8,np.array(ndvi.loc[ndvi.index == validation.loc[a,:].harvest,a])[0]+30,'Harvest',fontsize=12)\n",
    "\n",
    "\n",
    "plt.axvline(validation.loc[a,:].seeding, color='brown', linestyle='--', label='Planting')\n",
    "plt.axvline(validation.loc[a,:].harvest, color='red', linestyle='--', label='Harvest')\n",
    "plt.xlabel('Day of Year',fontsize=12)\n",
    "plt.ylabel('EO Index Values',fontsize=12)\n",
    "plt.legend(fontsize = 12)\n",
    "plt.grid()\n",
    "plt.savefig('indices_validation.png', dpi=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "l=list()\n",
    "for i in range(1,ndwi.shape[0]+1):\n",
    "    l.append('NDWI_DOY%s'%i)\n",
    "ndwi.index = np.asarray(l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "l=list()\n",
    "for i in range(1,psri.shape[0]+1):\n",
    "    l.append('PSRI_DOY%s'%i)\n",
    "psri.index = np.asarray(l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_ndwi_psri = pd.concat([ndwi.T,psri.T],axis=1)\n",
    "data_ndwi_psri.to_csv('C:/Users/antonkout/Documents/Work/antonkoutpersolnal/data_ndwi_psri.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}