HankelMatrix_and_DMDalg__model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Andernach</th>\n",
       "      <th>Bingen</th>\n",
       "      <th>Bonn</th>\n",
       "      <th>Frankfurt Osthafen</th>\n",
       "      <th>Kalkofen Neu</th>\n",
       "      <th>Kaub</th>\n",
       "      <th>Koblenz</th>\n",
       "      <th>Koblenz Up</th>\n",
       "      <th>Mainz</th>\n",
       "      <th>Oberwinter</th>\n",
       "      <th>Oestrich</th>\n",
       "      <th>Raunheim</th>\n",
       "      <th>Rockenau Ska</th>\n",
       "      <th>Speyer</th>\n",
       "      <th>Worms</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1980-12-31 23:15:00+00:00</th>\n",
       "      <td>307.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>199.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>289.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>153.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-12-31 23:30:00+00:00</th>\n",
       "      <td>307.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>199.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>290.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>153.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-12-31 23:45:00+00:00</th>\n",
       "      <td>307.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>198.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>289.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>153.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1981-01-01 00:00:00+00:00</th>\n",
       "      <td>307.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>198.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>289.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>153.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1981-01-01 00:15:00+00:00</th>\n",
       "      <td>307.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>198.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>289.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>153.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           Andernach  Bingen  Bonn  Frankfurt Osthafen  \\\n",
       "1980-12-31 23:15:00+00:00      307.0     NaN   NaN                 NaN   \n",
       "1980-12-31 23:30:00+00:00      307.0     NaN   NaN                 NaN   \n",
       "1980-12-31 23:45:00+00:00      307.0     NaN   NaN                 NaN   \n",
       "1981-01-01 00:00:00+00:00      307.0     NaN   NaN                 NaN   \n",
       "1981-01-01 00:15:00+00:00      307.0     NaN   NaN                 NaN   \n",
       "\n",
       "                           Kalkofen Neu   Kaub  Koblenz  Koblenz Up  Mainz  \\\n",
       "1980-12-31 23:15:00+00:00           NaN  199.0      NaN         NaN  289.0   \n",
       "1980-12-31 23:30:00+00:00           NaN  199.0      NaN         NaN  290.0   \n",
       "1980-12-31 23:45:00+00:00           NaN  198.0      NaN         NaN  289.0   \n",
       "1981-01-01 00:00:00+00:00           NaN  198.0      NaN         NaN  289.0   \n",
       "1981-01-01 00:15:00+00:00           NaN  198.0      NaN         NaN  289.0   \n",
       "\n",
       "                           Oberwinter  Oestrich  Raunheim  Rockenau Ska  \\\n",
       "1980-12-31 23:15:00+00:00         NaN       NaN       NaN           NaN   \n",
       "1980-12-31 23:30:00+00:00         NaN       NaN       NaN           NaN   \n",
       "1980-12-31 23:45:00+00:00         NaN       NaN       NaN           NaN   \n",
       "1981-01-01 00:00:00+00:00         NaN       NaN       NaN           NaN   \n",
       "1981-01-01 00:15:00+00:00         NaN       NaN       NaN           NaN   \n",
       "\n",
       "                           Speyer  Worms  \n",
       "1980-12-31 23:15:00+00:00     NaN  153.0  \n",
       "1980-12-31 23:30:00+00:00     NaN  153.0  \n",
       "1980-12-31 23:45:00+00:00     NaN  153.0  \n",
       "1981-01-01 00:00:00+00:00     NaN  153.0  \n",
       "1981-01-01 00:15:00+00:00     NaN  153.0  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_pickle(\"river.pickle\")\n",
    "\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that we have observations every 15 minutes. Thus in one day we have in total 4 * 24 observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Andernach</th>\n",
       "      <th>Bingen</th>\n",
       "      <th>Bonn</th>\n",
       "      <th>Frankfurt Osthafen</th>\n",
       "      <th>Kalkofen Neu</th>\n",
       "      <th>Kaub</th>\n",
       "      <th>Koblenz</th>\n",
       "      <th>Koblenz Up</th>\n",
       "      <th>Mainz</th>\n",
       "      <th>Oberwinter</th>\n",
       "      <th>Oestrich</th>\n",
       "      <th>Raunheim</th>\n",
       "      <th>Rockenau Ska</th>\n",
       "      <th>Speyer</th>\n",
       "      <th>Worms</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2014-12-31 22:00:00+00:00</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>217.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-12-31 22:15:00+00:00</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>217.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-12-31 22:30:00+00:00</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>216.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-12-31 22:45:00+00:00</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>216.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-12-31 23:00:00+00:00</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>216.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           Andernach  Bingen  Bonn  Frankfurt Osthafen  \\\n",
       "2014-12-31 22:00:00+00:00        NaN     NaN   NaN                 NaN   \n",
       "2014-12-31 22:15:00+00:00        NaN     NaN   NaN                 NaN   \n",
       "2014-12-31 22:30:00+00:00        NaN     NaN   NaN                 NaN   \n",
       "2014-12-31 22:45:00+00:00        NaN     NaN   NaN                 NaN   \n",
       "2014-12-31 23:00:00+00:00        NaN     NaN   NaN                 NaN   \n",
       "\n",
       "                           Kalkofen Neu  Kaub  Koblenz  Koblenz Up  Mainz  \\\n",
       "2014-12-31 22:00:00+00:00           NaN   NaN      NaN       217.0    NaN   \n",
       "2014-12-31 22:15:00+00:00           NaN   NaN      NaN       217.0    NaN   \n",
       "2014-12-31 22:30:00+00:00           NaN   NaN      NaN       216.0    NaN   \n",
       "2014-12-31 22:45:00+00:00           NaN   NaN      NaN       216.0    NaN   \n",
       "2014-12-31 23:00:00+00:00           NaN   NaN      NaN       216.0    NaN   \n",
       "\n",
       "                           Oberwinter  Oestrich  Raunheim  Rockenau Ska  \\\n",
       "2014-12-31 22:00:00+00:00         NaN       NaN       NaN           NaN   \n",
       "2014-12-31 22:15:00+00:00         NaN       NaN       NaN           NaN   \n",
       "2014-12-31 22:30:00+00:00         NaN       NaN       NaN           NaN   \n",
       "2014-12-31 22:45:00+00:00         NaN       NaN       NaN           NaN   \n",
       "2014-12-31 23:00:00+00:00         NaN       NaN       NaN           NaN   \n",
       "\n",
       "                           Speyer  Worms  \n",
       "2014-12-31 22:00:00+00:00     NaN    NaN  \n",
       "2014-12-31 22:15:00+00:00     NaN    NaN  \n",
       "2014-12-31 22:30:00+00:00     NaN    NaN  \n",
       "2014-12-31 22:45:00+00:00     NaN    NaN  \n",
       "2014-12-31 23:00:00+00:00     NaN    NaN  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1192128, 15)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The task "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The task is to predict the rivel level one step in the future."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The missing values / The dataset used for the prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to memory limitations, we use only the last 30 days ( December 2014 ) of the dataset in order to make our prediction.\n",
    "\n",
    "First we observe that in December 2014, there are a lot missing values in the dataset. We replace the NaN values by averages of  time, date and city in the last 20 years.\n",
    "\n",
    "For example, in the table above we see that the last observation for Andernach is a missing value. We replace it by the average of all non-missing values of Andernach at the same date and time in the last 20 years.\n",
    "\n",
    "After we fix the missing values problem, we define the Hankel matrix, based on this \"new\" dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "for row_index in range(1192128 - 4*24*30, 1192128) :  # the index runs over the observations of the last 30 days of the dataset\n",
    "    \n",
    "    # now we work for each row separately\n",
    "    v = df.iloc[row_index]\n",
    "    for j in range(15) :\n",
    "        if np.isnan(v[j]) == True :  # if it is a NaN value\n",
    "            v_past_av = 0 ; n = 0\n",
    "            for i in range(3):\n",
    "                if np.isnan(df.iloc[row_index - i*4*24*365][j]) == False :  # if it is a number\n",
    "                    v_past_av += df.iloc[row_index - i*4*24*365][j] \n",
    "                    n += 1\n",
    "            for i in range(4):\n",
    "                 if np.isnan(df.iloc[row_index - (i+3)*4*24*365 - 4*24 ][j]) == False :\n",
    "                    v_past_av += df.iloc[row_index - (i+3)*4*24*365 - 4*24 ][j]\n",
    "                    n += 1\n",
    "            for i in range(4) :\n",
    "                if np.isnan(df.iloc[row_index - (i+7)*4*24*365 - 2*(4*24) ][j]) == False :\n",
    "                    v_past_av += df.iloc[row_index - (i+7)*4*24*365 - 2*(4*24) ][j]\n",
    "                    n += 1\n",
    "            for i in range(4) :\n",
    "                if np.isnan( df.iloc[row_index - (i+11)*4*24*365 - 3*(4*24) ][j]) == False :\n",
    "                    v_past_av += df.iloc[row_index - (i+11)*4*24*365 - 3*(4*24) ][j]\n",
    "                    n += 1\n",
    "            for i in range(3) :\n",
    "                if np.isnan(df.iloc[row_index - (i+15)*4*24*365 - 4*(4*24) ][j]) == False  :\n",
    "                    v_past_av += df.iloc[row_index - (i+15)*4*24*365 - 4*(4*24) ][j]\n",
    "                    n += 1\n",
    "            for i in range(2) :   \n",
    "                if np.isnan(df.iloc[row_index - (i+19)*4*24*365 - 5*(4*24) ][j]) == False :\n",
    "                    v_past_av += df.iloc[row_index - (i+19)*4*24*365 - 5*(4*24) ][j]\n",
    "                    n +=1\n",
    "            v[j] = (v_past_av)/n\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Andernach</th>\n",
       "      <th>Bingen</th>\n",
       "      <th>Bonn</th>\n",
       "      <th>Frankfurt Osthafen</th>\n",
       "      <th>Kalkofen Neu</th>\n",
       "      <th>Kaub</th>\n",
       "      <th>Koblenz</th>\n",
       "      <th>Koblenz Up</th>\n",
       "      <th>Mainz</th>\n",
       "      <th>Oberwinter</th>\n",
       "      <th>Oestrich</th>\n",
       "      <th>Raunheim</th>\n",
       "      <th>Rockenau Ska</th>\n",
       "      <th>Speyer</th>\n",
       "      <th>Worms</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2014-12-01 23:15:00+00:00</th>\n",
       "      <td>241.055556</td>\n",
       "      <td>176.833333</td>\n",
       "      <td>273.055556</td>\n",
       "      <td>172.833333</td>\n",
       "      <td>226.444444</td>\n",
       "      <td>181.333333</td>\n",
       "      <td>198.444444</td>\n",
       "      <td>142.0</td>\n",
       "      <td>261.388889</td>\n",
       "      <td>204.222222</td>\n",
       "      <td>162.277778</td>\n",
       "      <td>142.277778</td>\n",
       "      <td>236.611111</td>\n",
       "      <td>319.555556</td>\n",
       "      <td>162.722222</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-12-01 23:30:00+00:00</th>\n",
       "      <td>240.888889</td>\n",
       "      <td>176.611111</td>\n",
       "      <td>273.000000</td>\n",
       "      <td>172.777778</td>\n",
       "      <td>226.666667</td>\n",
       "      <td>181.333333</td>\n",
       "      <td>198.111111</td>\n",
       "      <td>143.0</td>\n",
       "      <td>261.444444</td>\n",
       "      <td>204.388889</td>\n",
       "      <td>162.388889</td>\n",
       "      <td>142.611111</td>\n",
       "      <td>236.555556</td>\n",
       "      <td>319.722222</td>\n",
       "      <td>162.555556</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-12-01 23:45:00+00:00</th>\n",
       "      <td>240.944444</td>\n",
       "      <td>176.555556</td>\n",
       "      <td>272.888889</td>\n",
       "      <td>172.555556</td>\n",
       "      <td>226.666667</td>\n",
       "      <td>181.555556</td>\n",
       "      <td>198.055556</td>\n",
       "      <td>144.0</td>\n",
       "      <td>261.388889</td>\n",
       "      <td>204.444444</td>\n",
       "      <td>162.166667</td>\n",
       "      <td>142.166667</td>\n",
       "      <td>236.611111</td>\n",
       "      <td>319.833333</td>\n",
       "      <td>162.722222</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-12-02 00:00:00+00:00</th>\n",
       "      <td>240.944444</td>\n",
       "      <td>176.611111</td>\n",
       "      <td>273.111111</td>\n",
       "      <td>172.055556</td>\n",
       "      <td>226.833333</td>\n",
       "      <td>181.444444</td>\n",
       "      <td>198.166667</td>\n",
       "      <td>144.0</td>\n",
       "      <td>261.722222</td>\n",
       "      <td>204.500000</td>\n",
       "      <td>162.277778</td>\n",
       "      <td>142.444444</td>\n",
       "      <td>236.888889</td>\n",
       "      <td>319.611111</td>\n",
       "      <td>162.944444</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-12-02 00:15:00+00:00</th>\n",
       "      <td>241.055556</td>\n",
       "      <td>176.500000</td>\n",
       "      <td>273.166667</td>\n",
       "      <td>172.333333</td>\n",
       "      <td>227.111111</td>\n",
       "      <td>181.444444</td>\n",
       "      <td>197.833333</td>\n",
       "      <td>144.0</td>\n",
       "      <td>262.166667</td>\n",
       "      <td>204.444444</td>\n",
       "      <td>162.166667</td>\n",
       "      <td>142.000000</td>\n",
       "      <td>236.833333</td>\n",
       "      <td>319.722222</td>\n",
       "      <td>162.555556</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            Andernach      Bingen        Bonn  \\\n",
       "2014-12-01 23:15:00+00:00  241.055556  176.833333  273.055556   \n",
       "2014-12-01 23:30:00+00:00  240.888889  176.611111  273.000000   \n",
       "2014-12-01 23:45:00+00:00  240.944444  176.555556  272.888889   \n",
       "2014-12-02 00:00:00+00:00  240.944444  176.611111  273.111111   \n",
       "2014-12-02 00:15:00+00:00  241.055556  176.500000  273.166667   \n",
       "\n",
       "                           Frankfurt Osthafen  Kalkofen Neu        Kaub  \\\n",
       "2014-12-01 23:15:00+00:00          172.833333    226.444444  181.333333   \n",
       "2014-12-01 23:30:00+00:00          172.777778    226.666667  181.333333   \n",
       "2014-12-01 23:45:00+00:00          172.555556    226.666667  181.555556   \n",
       "2014-12-02 00:00:00+00:00          172.055556    226.833333  181.444444   \n",
       "2014-12-02 00:15:00+00:00          172.333333    227.111111  181.444444   \n",
       "\n",
       "                              Koblenz  Koblenz Up       Mainz  Oberwinter  \\\n",
       "2014-12-01 23:15:00+00:00  198.444444       142.0  261.388889  204.222222   \n",
       "2014-12-01 23:30:00+00:00  198.111111       143.0  261.444444  204.388889   \n",
       "2014-12-01 23:45:00+00:00  198.055556       144.0  261.388889  204.444444   \n",
       "2014-12-02 00:00:00+00:00  198.166667       144.0  261.722222  204.500000   \n",
       "2014-12-02 00:15:00+00:00  197.833333       144.0  262.166667  204.444444   \n",
       "\n",
       "                             Oestrich    Raunheim  Rockenau Ska      Speyer  \\\n",
       "2014-12-01 23:15:00+00:00  162.277778  142.277778    236.611111  319.555556   \n",
       "2014-12-01 23:30:00+00:00  162.388889  142.611111    236.555556  319.722222   \n",
       "2014-12-01 23:45:00+00:00  162.166667  142.166667    236.611111  319.833333   \n",
       "2014-12-02 00:00:00+00:00  162.277778  142.444444    236.888889  319.611111   \n",
       "2014-12-02 00:15:00+00:00  162.166667  142.000000    236.833333  319.722222   \n",
       "\n",
       "                                Worms  \n",
       "2014-12-01 23:15:00+00:00  162.722222  \n",
       "2014-12-01 23:30:00+00:00  162.555556  \n",
       "2014-12-01 23:45:00+00:00  162.722222  \n",
       "2014-12-02 00:00:00+00:00  162.944444  \n",
       "2014-12-02 00:15:00+00:00  162.555556  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df1 = df[1192128 - 4*24*30:1192128]\n",
    "df1.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Andernach</th>\n",
       "      <th>Bingen</th>\n",
       "      <th>Bonn</th>\n",
       "      <th>Frankfurt Osthafen</th>\n",
       "      <th>Kalkofen Neu</th>\n",
       "      <th>Kaub</th>\n",
       "      <th>Koblenz</th>\n",
       "      <th>Koblenz Up</th>\n",
       "      <th>Mainz</th>\n",
       "      <th>Oberwinter</th>\n",
       "      <th>Oestrich</th>\n",
       "      <th>Raunheim</th>\n",
       "      <th>Rockenau Ska</th>\n",
       "      <th>Speyer</th>\n",
       "      <th>Worms</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2014-12-31 22:00:00+00:00</th>\n",
       "      <td>368.888889</td>\n",
       "      <td>237.055556</td>\n",
       "      <td>400.277778</td>\n",
       "      <td>192.388889</td>\n",
       "      <td>301.333333</td>\n",
       "      <td>266.000000</td>\n",
       "      <td>307.888889</td>\n",
       "      <td>217.0</td>\n",
       "      <td>328.833333</td>\n",
       "      <td>321.444444</td>\n",
       "      <td>219.722222</td>\n",
       "      <td>161.666667</td>\n",
       "      <td>259.500000</td>\n",
       "      <td>389.944444</td>\n",
       "      <td>237.888889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-12-31 22:15:00+00:00</th>\n",
       "      <td>368.722222</td>\n",
       "      <td>237.111111</td>\n",
       "      <td>400.388889</td>\n",
       "      <td>192.055556</td>\n",
       "      <td>301.166667</td>\n",
       "      <td>265.944444</td>\n",
       "      <td>307.777778</td>\n",
       "      <td>217.0</td>\n",
       "      <td>328.722222</td>\n",
       "      <td>321.444444</td>\n",
       "      <td>219.666667</td>\n",
       "      <td>161.944444</td>\n",
       "      <td>258.944444</td>\n",
       "      <td>389.944444</td>\n",
       "      <td>237.888889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-12-31 22:30:00+00:00</th>\n",
       "      <td>368.666667</td>\n",
       "      <td>237.166667</td>\n",
       "      <td>400.555556</td>\n",
       "      <td>191.833333</td>\n",
       "      <td>301.055556</td>\n",
       "      <td>265.888889</td>\n",
       "      <td>307.833333</td>\n",
       "      <td>216.0</td>\n",
       "      <td>328.833333</td>\n",
       "      <td>321.444444</td>\n",
       "      <td>219.500000</td>\n",
       "      <td>162.333333</td>\n",
       "      <td>258.333333</td>\n",
       "      <td>389.833333</td>\n",
       "      <td>237.888889</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-12-31 22:45:00+00:00</th>\n",
       "      <td>368.555556</td>\n",
       "      <td>237.277778</td>\n",
       "      <td>400.555556</td>\n",
       "      <td>191.833333</td>\n",
       "      <td>301.000000</td>\n",
       "      <td>265.944444</td>\n",
       "      <td>307.777778</td>\n",
       "      <td>216.0</td>\n",
       "      <td>328.666667</td>\n",
       "      <td>321.277778</td>\n",
       "      <td>219.555556</td>\n",
       "      <td>162.833333</td>\n",
       "      <td>257.888889</td>\n",
       "      <td>389.611111</td>\n",
       "      <td>238.111111</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2014-12-31 23:00:00+00:00</th>\n",
       "      <td>368.555556</td>\n",
       "      <td>237.111111</td>\n",
       "      <td>400.277778</td>\n",
       "      <td>191.833333</td>\n",
       "      <td>301.166667</td>\n",
       "      <td>265.888889</td>\n",
       "      <td>307.500000</td>\n",
       "      <td>216.0</td>\n",
       "      <td>328.722222</td>\n",
       "      <td>321.222222</td>\n",
       "      <td>219.611111</td>\n",
       "      <td>163.111111</td>\n",
       "      <td>257.555556</td>\n",
       "      <td>389.500000</td>\n",
       "      <td>237.944444</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            Andernach      Bingen        Bonn  \\\n",
       "2014-12-31 22:00:00+00:00  368.888889  237.055556  400.277778   \n",
       "2014-12-31 22:15:00+00:00  368.722222  237.111111  400.388889   \n",
       "2014-12-31 22:30:00+00:00  368.666667  237.166667  400.555556   \n",
       "2014-12-31 22:45:00+00:00  368.555556  237.277778  400.555556   \n",
       "2014-12-31 23:00:00+00:00  368.555556  237.111111  400.277778   \n",
       "\n",
       "                           Frankfurt Osthafen  Kalkofen Neu        Kaub  \\\n",
       "2014-12-31 22:00:00+00:00          192.388889    301.333333  266.000000   \n",
       "2014-12-31 22:15:00+00:00          192.055556    301.166667  265.944444   \n",
       "2014-12-31 22:30:00+00:00          191.833333    301.055556  265.888889   \n",
       "2014-12-31 22:45:00+00:00          191.833333    301.000000  265.944444   \n",
       "2014-12-31 23:00:00+00:00          191.833333    301.166667  265.888889   \n",
       "\n",
       "                              Koblenz  Koblenz Up       Mainz  Oberwinter  \\\n",
       "2014-12-31 22:00:00+00:00  307.888889       217.0  328.833333  321.444444   \n",
       "2014-12-31 22:15:00+00:00  307.777778       217.0  328.722222  321.444444   \n",
       "2014-12-31 22:30:00+00:00  307.833333       216.0  328.833333  321.444444   \n",
       "2014-12-31 22:45:00+00:00  307.777778       216.0  328.666667  321.277778   \n",
       "2014-12-31 23:00:00+00:00  307.500000       216.0  328.722222  321.222222   \n",
       "\n",
       "                             Oestrich    Raunheim  Rockenau Ska      Speyer  \\\n",
       "2014-12-31 22:00:00+00:00  219.722222  161.666667    259.500000  389.944444   \n",
       "2014-12-31 22:15:00+00:00  219.666667  161.944444    258.944444  389.944444   \n",
       "2014-12-31 22:30:00+00:00  219.500000  162.333333    258.333333  389.833333   \n",
       "2014-12-31 22:45:00+00:00  219.555556  162.833333    257.888889  389.611111   \n",
       "2014-12-31 23:00:00+00:00  219.611111  163.111111    257.555556  389.500000   \n",
       "\n",
       "                                Worms  \n",
       "2014-12-31 22:00:00+00:00  237.888889  \n",
       "2014-12-31 22:15:00+00:00  237.888889  \n",
       "2014-12-31 22:30:00+00:00  237.888889  \n",
       "2014-12-31 22:45:00+00:00  238.111111  \n",
       "2014-12-31 23:00:00+00:00  237.944444  "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df1.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# are there any missing values in the dataset ?\n",
    "df1.isnull().values.any()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The \"time-delayed\" Hankel matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "val = df1.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = 2    # time steps\n",
    "n = 940\n",
    "l = 1000   # length of observations (rows) taken into account for each feature (column)\n",
    "H = np.zeros((val.shape[1]*n,l))\n",
    "for i in range(n):\n",
    "    for j in range(val.shape[1]):\n",
    "        H[i*val.shape[1]+j,:] = val[i*h:i*h+l,j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(14100, 1000)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "u,s,v = np.linalg.svd(H)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "/* Put everything inside the global mpl namespace */\n",
       "window.mpl = {};\n",
       "\n",
       "\n",
       "mpl.get_websocket_type = function() {\n",
       "    if (typeof(WebSocket) !== 'undefined') {\n",
       "        return WebSocket;\n",
       "    } else if (typeof(MozWebSocket) !== 'undefined') {\n",
       "        return MozWebSocket;\n",
       "    } else {\n",
       "        alert('Your browser does not have WebSocket support. ' +\n",
       "              'Please try Chrome, Safari or Firefox ≥ 6. ' +\n",
       "              'Firefox 4 and 5 are also supported but you ' +\n",
       "              'have to enable WebSockets in about:config.');\n",
       "    };\n",
       "}\n",
       "\n",
       "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n",
       "    this.id = figure_id;\n",
       "\n",
       "    this.ws = websocket;\n",
       "\n",
       "    this.supports_binary = (this.ws.binaryType != undefined);\n",
       "\n",
       "    if (!this.supports_binary) {\n",
       "        var warnings = document.getElementById(\"mpl-warnings\");\n",
       "        if (warnings) {\n",
       "            warnings.style.display = 'block';\n",
       "            warnings.textContent = (\n",
       "                \"This browser does not support binary websocket messages. \" +\n",
       "                    \"Performance may be slow.\");\n",
       "        }\n",
       "    }\n",
       "\n",
       "    this.imageObj = new Image();\n",
       "\n",
       "    this.context = undefined;\n",
       "    this.message = undefined;\n",
       "    this.canvas = undefined;\n",
       "    this.rubberband_canvas = undefined;\n",
       "    this.rubberband_context = undefined;\n",
       "    this.format_dropdown = undefined;\n",
       "\n",
       "    this.image_mode = 'full';\n",
       "\n",
       "    this.root = $('<div/>');\n",
       "    this._root_extra_style(this.root)\n",
       "    this.root.attr('style', 'display: inline-block');\n",
       "\n",
       "    $(parent_element).append(this.root);\n",
       "\n",
       "    this._init_header(this);\n",
       "    this._init_canvas(this);\n",
       "    this._init_toolbar(this);\n",
       "\n",
       "    var fig = this;\n",
       "\n",
       "    this.waiting = false;\n",
       "\n",
       "    this.ws.onopen =  function () {\n",
       "            fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n",
       "            fig.send_message(\"send_image_mode\", {});\n",
       "            if (mpl.ratio != 1) {\n",
       "                fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n",
       "            }\n",
       "            fig.send_message(\"refresh\", {});\n",
       "        }\n",
       "\n",
       "    this.imageObj.onload = function() {\n",
       "            if (fig.image_mode == 'full') {\n",
       "                // Full images could contain transparency (where diff images\n",
       "                // almost always do), so we need to clear the canvas so that\n",
       "                // there is no ghosting.\n",
       "                fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n",
       "            }\n",
       "            fig.context.drawImage(fig.imageObj, 0, 0);\n",
       "        };\n",
       "\n",
       "    this.imageObj.onunload = function() {\n",
       "        fig.ws.close();\n",
       "    }\n",
       "\n",
       "    this.ws.onmessage = this._make_on_message_function(this);\n",
       "\n",
       "    this.ondownload = ondownload;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_header = function() {\n",
       "    var titlebar = $(\n",
       "        '<div class=\"ui-dialog-titlebar ui-widget-header ui-corner-all ' +\n",
       "        'ui-helper-clearfix\"/>');\n",
       "    var titletext = $(\n",
       "        '<div class=\"ui-dialog-title\" style=\"width: 100%; ' +\n",
       "        'text-align: center; padding: 3px;\"/>');\n",
       "    titlebar.append(titletext)\n",
       "    this.root.append(titlebar);\n",
       "    this.header = titletext[0];\n",
       "}\n",
       "\n",
       "\n",
       "\n",
       "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n",
       "\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n",
       "\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_canvas = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var canvas_div = $('<div/>');\n",
       "\n",
       "    canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n",
       "\n",
       "    function canvas_keyboard_event(event) {\n",
       "        return fig.key_event(event, event['data']);\n",
       "    }\n",
       "\n",
       "    canvas_div.keydown('key_press', canvas_keyboard_event);\n",
       "    canvas_div.keyup('key_release', canvas_keyboard_event);\n",
       "    this.canvas_div = canvas_div\n",
       "    this._canvas_extra_style(canvas_div)\n",
       "    this.root.append(canvas_div);\n",
       "\n",
       "    var canvas = $('<canvas/>');\n",
       "    canvas.addClass('mpl-canvas');\n",
       "    canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n",
       "\n",
       "    this.canvas = canvas[0];\n",
       "    this.context = canvas[0].getContext(\"2d\");\n",
       "\n",
       "    var backingStore = this.context.backingStorePixelRatio ||\n",
       "\tthis.context.webkitBackingStorePixelRatio ||\n",
       "\tthis.context.mozBackingStorePixelRatio ||\n",
       "\tthis.context.msBackingStorePixelRatio ||\n",
       "\tthis.context.oBackingStorePixelRatio ||\n",
       "\tthis.context.backingStorePixelRatio || 1;\n",
       "\n",
       "    mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n",
       "\n",
       "    var rubberband = $('<canvas/>');\n",
       "    rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n",
       "\n",
       "    var pass_mouse_events = true;\n",
       "\n",
       "    canvas_div.resizable({\n",
       "        start: function(event, ui) {\n",
       "            pass_mouse_events = false;\n",
       "        },\n",
       "        resize: function(event, ui) {\n",
       "            fig.request_resize(ui.size.width, ui.size.height);\n",
       "        },\n",
       "        stop: function(event, ui) {\n",
       "            pass_mouse_events = true;\n",
       "            fig.request_resize(ui.size.width, ui.size.height);\n",
       "        },\n",
       "    });\n",
       "\n",
       "    function mouse_event_fn(event) {\n",
       "        if (pass_mouse_events)\n",
       "            return fig.mouse_event(event, event['data']);\n",
       "    }\n",
       "\n",
       "    rubberband.mousedown('button_press', mouse_event_fn);\n",
       "    rubberband.mouseup('button_release', mouse_event_fn);\n",
       "    // Throttle sequential mouse events to 1 every 20ms.\n",
       "    rubberband.mousemove('motion_notify', mouse_event_fn);\n",
       "\n",
       "    rubberband.mouseenter('figure_enter', mouse_event_fn);\n",
       "    rubberband.mouseleave('figure_leave', mouse_event_fn);\n",
       "\n",
       "    canvas_div.on(\"wheel\", function (event) {\n",
       "        event = event.originalEvent;\n",
       "        event['data'] = 'scroll'\n",
       "        if (event.deltaY < 0) {\n",
       "            event.step = 1;\n",
       "        } else {\n",
       "            event.step = -1;\n",
       "        }\n",
       "        mouse_event_fn(event);\n",
       "    });\n",
       "\n",
       "    canvas_div.append(canvas);\n",
       "    canvas_div.append(rubberband);\n",
       "\n",
       "    this.rubberband = rubberband;\n",
       "    this.rubberband_canvas = rubberband[0];\n",
       "    this.rubberband_context = rubberband[0].getContext(\"2d\");\n",
       "    this.rubberband_context.strokeStyle = \"#000000\";\n",
       "\n",
       "    this._resize_canvas = function(width, height) {\n",
       "        // Keep the size of the canvas, canvas container, and rubber band\n",
       "        // canvas in synch.\n",
       "        canvas_div.css('width', width)\n",
       "        canvas_div.css('height', height)\n",
       "\n",
       "        canvas.attr('width', width * mpl.ratio);\n",
       "        canvas.attr('height', height * mpl.ratio);\n",
       "        canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n",
       "\n",
       "        rubberband.attr('width', width);\n",
       "        rubberband.attr('height', height);\n",
       "    }\n",
       "\n",
       "    // Set the figure to an initial 600x600px, this will subsequently be updated\n",
       "    // upon first draw.\n",
       "    this._resize_canvas(600, 600);\n",
       "\n",
       "    // Disable right mouse context menu.\n",
       "    $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n",
       "        return false;\n",
       "    });\n",
       "\n",
       "    function set_focus () {\n",
       "        canvas.focus();\n",
       "        canvas_div.focus();\n",
       "    }\n",
       "\n",
       "    window.setTimeout(set_focus, 100);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_toolbar = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var nav_element = $('<div/>');\n",
       "    nav_element.attr('style', 'width: 100%');\n",
       "    this.root.append(nav_element);\n",
       "\n",
       "    // Define a callback function for later on.\n",
       "    function toolbar_event(event) {\n",
       "        return fig.toolbar_button_onclick(event['data']);\n",
       "    }\n",
       "    function toolbar_mouse_event(event) {\n",
       "        return fig.toolbar_button_onmouseover(event['data']);\n",
       "    }\n",
       "\n",
       "    for(var toolbar_ind in mpl.toolbar_items) {\n",
       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
       "\n",
       "        if (!name) {\n",
       "            // put a spacer in here.\n",
       "            continue;\n",
       "        }\n",
       "        var button = $('<button/>');\n",
       "        button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' +\n",
       "                        'ui-button-icon-only');\n",
       "        button.attr('role', 'button');\n",
       "        button.attr('aria-disabled', 'false');\n",
       "        button.click(method_name, toolbar_event);\n",
       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
       "\n",
       "        var icon_img = $('<span/>');\n",
       "        icon_img.addClass('ui-button-icon-primary ui-icon');\n",
       "        icon_img.addClass(image);\n",
       "        icon_img.addClass('ui-corner-all');\n",
       "\n",
       "        var tooltip_span = $('<span/>');\n",
       "        tooltip_span.addClass('ui-button-text');\n",
       "        tooltip_span.html(tooltip);\n",
       "\n",
       "        button.append(icon_img);\n",
       "        button.append(tooltip_span);\n",
       "\n",
       "        nav_element.append(button);\n",
       "    }\n",
       "\n",
       "    var fmt_picker_span = $('<span/>');\n",
       "\n",
       "    var fmt_picker = $('<select/>');\n",
       "    fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content');\n",
       "    fmt_picker_span.append(fmt_picker);\n",
       "    nav_element.append(fmt_picker_span);\n",
       "    this.format_dropdown = fmt_picker[0];\n",
       "\n",
       "    for (var ind in mpl.extensions) {\n",
       "        var fmt = mpl.extensions[ind];\n",
       "        var option = $(\n",
       "            '<option/>', {selected: fmt === mpl.default_extension}).html(fmt);\n",
       "        fmt_picker.append(option);\n",
       "    }\n",
       "\n",
       "    // Add hover states to the ui-buttons\n",
       "    $( \".ui-button\" ).hover(\n",
       "        function() { $(this).addClass(\"ui-state-hover\");},\n",
       "        function() { $(this).removeClass(\"ui-state-hover\");}\n",
       "    );\n",
       "\n",
       "    var status_bar = $('<span class=\"mpl-message\"/>');\n",
       "    nav_element.append(status_bar);\n",
       "    this.message = status_bar[0];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.request_resize = function(x_pixels, y_pixels) {\n",
       "    // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n",
       "    // which will in turn request a refresh of the image.\n",
       "    this.send_message('resize', {'width': x_pixels, 'height': y_pixels});\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.send_message = function(type, properties) {\n",
       "    properties['type'] = type;\n",
       "    properties['figure_id'] = this.id;\n",
       "    this.ws.send(JSON.stringify(properties));\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.send_draw_message = function() {\n",
       "    if (!this.waiting) {\n",
       "        this.waiting = true;\n",
       "        this.ws.send(JSON.stringify({type: \"draw\", figure_id: this.id}));\n",
       "    }\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
       "    var format_dropdown = fig.format_dropdown;\n",
       "    var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n",
       "    fig.ondownload(fig, format);\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype.handle_resize = function(fig, msg) {\n",
       "    var size = msg['size'];\n",
       "    if (size[0] != fig.canvas.width || size[1] != fig.canvas.height) {\n",
       "        fig._resize_canvas(size[0], size[1]);\n",
       "        fig.send_message(\"refresh\", {});\n",
       "    };\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_rubberband = function(fig, msg) {\n",
       "    var x0 = msg['x0'] / mpl.ratio;\n",
       "    var y0 = (fig.canvas.height - msg['y0']) / mpl.ratio;\n",
       "    var x1 = msg['x1'] / mpl.ratio;\n",
       "    var y1 = (fig.canvas.height - msg['y1']) / mpl.ratio;\n",
       "    x0 = Math.floor(x0) + 0.5;\n",
       "    y0 = Math.floor(y0) + 0.5;\n",
       "    x1 = Math.floor(x1) + 0.5;\n",
       "    y1 = Math.floor(y1) + 0.5;\n",
       "    var min_x = Math.min(x0, x1);\n",
       "    var min_y = Math.min(y0, y1);\n",
       "    var width = Math.abs(x1 - x0);\n",
       "    var height = Math.abs(y1 - y0);\n",
       "\n",
       "    fig.rubberband_context.clearRect(\n",
       "        0, 0, fig.canvas.width, fig.canvas.height);\n",
       "\n",
       "    fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_figure_label = function(fig, msg) {\n",
       "    // Updates the figure title.\n",
       "    fig.header.textContent = msg['label'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_cursor = function(fig, msg) {\n",
       "    var cursor = msg['cursor'];\n",
       "    switch(cursor)\n",
       "    {\n",
       "    case 0:\n",
       "        cursor = 'pointer';\n",
       "        break;\n",
       "    case 1:\n",
       "        cursor = 'default';\n",
       "        break;\n",
       "    case 2:\n",
       "        cursor = 'crosshair';\n",
       "        break;\n",
       "    case 3:\n",
       "        cursor = 'move';\n",
       "        break;\n",
       "    }\n",
       "    fig.rubberband_canvas.style.cursor = cursor;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_message = function(fig, msg) {\n",
       "    fig.message.textContent = msg['message'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_draw = function(fig, msg) {\n",
       "    // Request the server to send over a new figure.\n",
       "    fig.send_draw_message();\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_image_mode = function(fig, msg) {\n",
       "    fig.image_mode = msg['mode'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.updated_canvas_event = function() {\n",
       "    // Called whenever the canvas gets updated.\n",
       "    this.send_message(\"ack\", {});\n",
       "}\n",
       "\n",
       "// A function to construct a web socket function for onmessage handling.\n",
       "// Called in the figure constructor.\n",
       "mpl.figure.prototype._make_on_message_function = function(fig) {\n",
       "    return function socket_on_message(evt) {\n",
       "        if (evt.data instanceof Blob) {\n",
       "            /* FIXME: We get \"Resource interpreted as Image but\n",
       "             * transferred with MIME type text/plain:\" errors on\n",
       "             * Chrome.  But how to set the MIME type?  It doesn't seem\n",
       "             * to be part of the websocket stream */\n",
       "            evt.data.type = \"image/png\";\n",
       "\n",
       "            /* Free the memory for the previous frames */\n",
       "            if (fig.imageObj.src) {\n",
       "                (window.URL || window.webkitURL).revokeObjectURL(\n",
       "                    fig.imageObj.src);\n",
       "            }\n",
       "\n",
       "            fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n",
       "                evt.data);\n",
       "            fig.updated_canvas_event();\n",
       "            fig.waiting = false;\n",
       "            return;\n",
       "        }\n",
       "        else if (typeof evt.data === 'string' && evt.data.slice(0, 21) == \"data:image/png;base64\") {\n",
       "            fig.imageObj.src = evt.data;\n",
       "            fig.updated_canvas_event();\n",
       "            fig.waiting = false;\n",
       "            return;\n",
       "        }\n",
       "\n",
       "        var msg = JSON.parse(evt.data);\n",
       "        var msg_type = msg['type'];\n",
       "\n",
       "        // Call the  \"handle_{type}\" callback, which takes\n",
       "        // the figure and JSON message as its only arguments.\n",
       "        try {\n",
       "            var callback = fig[\"handle_\" + msg_type];\n",
       "        } catch (e) {\n",
       "            console.log(\"No handler for the '\" + msg_type + \"' message type: \", msg);\n",
       "            return;\n",
       "        }\n",
       "\n",
       "        if (callback) {\n",
       "            try {\n",
       "                // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n",
       "                callback(fig, msg);\n",
       "            } catch (e) {\n",
       "                console.log(\"Exception inside the 'handler_\" + msg_type + \"' callback:\", e, e.stack, msg);\n",
       "            }\n",
       "        }\n",
       "    };\n",
       "}\n",
       "\n",
       "// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n",
       "mpl.findpos = function(e) {\n",
       "    //this section is from http://www.quirksmode.org/js/events_properties.html\n",
       "    var targ;\n",
       "    if (!e)\n",
       "        e = window.event;\n",
       "    if (e.target)\n",
       "        targ = e.target;\n",
       "    else if (e.srcElement)\n",
       "        targ = e.srcElement;\n",
       "    if (targ.nodeType == 3) // defeat Safari bug\n",
       "        targ = targ.parentNode;\n",
       "\n",
       "    // jQuery normalizes the pageX and pageY\n",
       "    // pageX,Y are the mouse positions relative to the document\n",
       "    // offset() returns the position of the element relative to the document\n",
       "    var x = e.pageX - $(targ).offset().left;\n",
       "    var y = e.pageY - $(targ).offset().top;\n",
       "\n",
       "    return {\"x\": x, \"y\": y};\n",
       "};\n",
       "\n",
       "/*\n",
       " * return a copy of an object with only non-object keys\n",
       " * we need this to avoid circular references\n",
       " * http://stackoverflow.com/a/24161582/3208463\n",
       " */\n",
       "function simpleKeys (original) {\n",
       "  return Object.keys(original).reduce(function (obj, key) {\n",
       "    if (typeof original[key] !== 'object')\n",
       "        obj[key] = original[key]\n",
       "    return obj;\n",
       "  }, {});\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.mouse_event = function(event, name) {\n",
       "    var canvas_pos = mpl.findpos(event)\n",
       "\n",
       "    if (name === 'button_press')\n",
       "    {\n",
       "        this.canvas.focus();\n",
       "        this.canvas_div.focus();\n",
       "    }\n",
       "\n",
       "    var x = canvas_pos.x * mpl.ratio;\n",
       "    var y = canvas_pos.y * mpl.ratio;\n",
       "\n",
       "    this.send_message(name, {x: x, y: y, button: event.button,\n",
       "                             step: event.step,\n",
       "                             guiEvent: simpleKeys(event)});\n",
       "\n",
       "    /* This prevents the web browser from automatically changing to\n",
       "     * the text insertion cursor when the button is pressed.  We want\n",
       "     * to control all of the cursor setting manually through the\n",
       "     * 'cursor' event from matplotlib */\n",
       "    event.preventDefault();\n",
       "    return false;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
       "    // Handle any extra behaviour associated with a key event\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.key_event = function(event, name) {\n",
       "\n",
       "    // Prevent repeat events\n",
       "    if (name == 'key_press')\n",
       "    {\n",
       "        if (event.which === this._key)\n",
       "            return;\n",
       "        else\n",
       "            this._key = event.which;\n",
       "    }\n",
       "    if (name == 'key_release')\n",
       "        this._key = null;\n",
       "\n",
       "    var value = '';\n",
       "    if (event.ctrlKey && event.which != 17)\n",
       "        value += \"ctrl+\";\n",
       "    if (event.altKey && event.which != 18)\n",
       "        value += \"alt+\";\n",
       "    if (event.shiftKey && event.which != 16)\n",
       "        value += \"shift+\";\n",
       "\n",
       "    value += 'k';\n",
       "    value += event.which.toString();\n",
       "\n",
       "    this._key_event_extra(event, name);\n",
       "\n",
       "    this.send_message(name, {key: value,\n",
       "                             guiEvent: simpleKeys(event)});\n",
       "    return false;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.toolbar_button_onclick = function(name) {\n",
       "    if (name == 'download') {\n",
       "        this.handle_save(this, null);\n",
       "    } else {\n",
       "        this.send_message(\"toolbar_button\", {name: name});\n",
       "    }\n",
       "};\n",
       "\n",
       "mpl.figure.prototype.toolbar_button_onmouseover = function(tooltip) {\n",
       "    this.message.textContent = tooltip;\n",
       "};\n",
       "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Pan axes with left mouse, zoom with right\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n",
       "\n",
       "mpl.extensions = [\"eps\", \"jpeg\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\"];\n",
       "\n",
       "mpl.default_extension = \"png\";var comm_websocket_adapter = function(comm) {\n",
       "    // Create a \"websocket\"-like object which calls the given IPython comm\n",
       "    // object with the appropriate methods. Currently this is a non binary\n",
       "    // socket, so there is still some room for performance tuning.\n",
       "    var ws = {};\n",
       "\n",
       "    ws.close = function() {\n",
       "        comm.close()\n",
       "    };\n",
       "    ws.send = function(m) {\n",
       "        //console.log('sending', m);\n",
       "        comm.send(m);\n",
       "    };\n",
       "    // Register the callback with on_msg.\n",
       "    comm.on_msg(function(msg) {\n",
       "        //console.log('receiving', msg['content']['data'], msg);\n",
       "        // Pass the mpl event to the overridden (by mpl) onmessage function.\n",
       "        ws.onmessage(msg['content']['data'])\n",
       "    });\n",
       "    return ws;\n",
       "}\n",
       "\n",
       "mpl.mpl_figure_comm = function(comm, msg) {\n",
       "    // This is the function which gets called when the mpl process\n",
       "    // starts-up an IPython Comm through the \"matplotlib\" channel.\n",
       "\n",
       "    var id = msg.content.data.id;\n",
       "    // Get hold of the div created by the display call when the Comm\n",
       "    // socket was opened in Python.\n",
       "    var element = $(\"#\" + id);\n",
       "    var ws_proxy = comm_websocket_adapter(comm)\n",
       "\n",
       "    function ondownload(figure, format) {\n",
       "        window.open(figure.imageObj.src);\n",
       "    }\n",
       "\n",
       "    var fig = new mpl.figure(id, ws_proxy,\n",
       "                           ondownload,\n",
       "                           element.get(0));\n",
       "\n",
       "    // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n",
       "    // web socket which is closed, not our websocket->open comm proxy.\n",
       "    ws_proxy.onopen();\n",
       "\n",
       "    fig.parent_element = element.get(0);\n",
       "    fig.cell_info = mpl.find_output_cell(\"<div id='\" + id + \"'></div>\");\n",
       "    if (!fig.cell_info) {\n",
       "        console.error(\"Failed to find cell for figure\", id, fig);\n",
       "        return;\n",
       "    }\n",
       "\n",
       "    var output_index = fig.cell_info[2]\n",
       "    var cell = fig.cell_info[0];\n",
       "\n",
       "};\n",
       "\n",
       "mpl.figure.prototype.handle_close = function(fig, msg) {\n",
       "    var width = fig.canvas.width/mpl.ratio\n",
       "    fig.root.unbind('remove')\n",
       "\n",
       "    // Update the output cell to use the data from the current canvas.\n",
       "    fig.push_to_output();\n",
       "    var dataURL = fig.canvas.toDataURL();\n",
       "    // Re-enable the keyboard manager in IPython - without this line, in FF,\n",
       "    // the notebook keyboard shortcuts fail.\n",
       "    IPython.keyboard_manager.enable()\n",
       "    $(fig.parent_element).html('<img src=\"' + dataURL + '\" width=\"' + width + '\">');\n",
       "    fig.close_ws(fig, msg);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.close_ws = function(fig, msg){\n",
       "    fig.send_message('closing', msg);\n",
       "    // fig.ws.close()\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.push_to_output = function(remove_interactive) {\n",
       "    // Turn the data on the canvas into data in the output cell.\n",
       "    var width = this.canvas.width/mpl.ratio\n",
       "    var dataURL = this.canvas.toDataURL();\n",
       "    this.cell_info[1]['text/html'] = '<img src=\"' + dataURL + '\" width=\"' + width + '\">';\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.updated_canvas_event = function() {\n",
       "    // Tell IPython that the notebook contents must change.\n",
       "    IPython.notebook.set_dirty(true);\n",
       "    this.send_message(\"ack\", {});\n",
       "    var fig = this;\n",
       "    // Wait a second, then push the new image to the DOM so\n",
       "    // that it is saved nicely (might be nice to debounce this).\n",
       "    setTimeout(function () { fig.push_to_output() }, 1000);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_toolbar = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var nav_element = $('<div/>');\n",
       "    nav_element.attr('style', 'width: 100%');\n",
       "    this.root.append(nav_element);\n",
       "\n",
       "    // Define a callback function for later on.\n",
       "    function toolbar_event(event) {\n",
       "        return fig.toolbar_button_onclick(event['data']);\n",
       "    }\n",
       "    function toolbar_mouse_event(event) {\n",
       "        return fig.toolbar_button_onmouseover(event['data']);\n",
       "    }\n",
       "\n",
       "    for(var toolbar_ind in mpl.toolbar_items){\n",
       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
       "\n",
       "        if (!name) { continue; };\n",
       "\n",
       "        var button = $('<button class=\"btn btn-default\" href=\"#\" title=\"' + name + '\"><i class=\"fa ' + image + ' fa-lg\"></i></button>');\n",
       "        button.click(method_name, toolbar_event);\n",
       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
       "        nav_element.append(button);\n",
       "    }\n",
       "\n",
       "    // Add the status bar.\n",
       "    var status_bar = $('<span class=\"mpl-message\" style=\"text-align:right; float: right;\"/>');\n",
       "    nav_element.append(status_bar);\n",
       "    this.message = status_bar[0];\n",
       "\n",
       "    // Add the close button to the window.\n",
       "    var buttongrp = $('<div class=\"btn-group inline pull-right\"></div>');\n",
       "    var button = $('<button class=\"btn btn-mini btn-primary\" href=\"#\" title=\"Stop Interaction\"><i class=\"fa fa-power-off icon-remove icon-large\"></i></button>');\n",
       "    button.click(function (evt) { fig.handle_close(fig, {}); } );\n",
       "    button.mouseover('Stop Interaction', toolbar_mouse_event);\n",
       "    buttongrp.append(button);\n",
       "    var titlebar = this.root.find($('.ui-dialog-titlebar'));\n",
       "    titlebar.prepend(buttongrp);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._root_extra_style = function(el){\n",
       "    var fig = this\n",
       "    el.on(\"remove\", function(){\n",
       "\tfig.close_ws(fig, {});\n",
       "    });\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._canvas_extra_style = function(el){\n",
       "    // this is important to make the div 'focusable\n",
       "    el.attr('tabindex', 0)\n",
       "    // reach out to IPython and tell the keyboard manager to turn it's self\n",
       "    // off when our div gets focus\n",
       "\n",
       "    // location in version 3\n",
       "    if (IPython.notebook.keyboard_manager) {\n",
       "        IPython.notebook.keyboard_manager.register_events(el);\n",
       "    }\n",
       "    else {\n",
       "        // location in version 2\n",
       "        IPython.keyboard_manager.register_events(el);\n",
       "    }\n",
       "\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
       "    var manager = IPython.notebook.keyboard_manager;\n",
       "    if (!manager)\n",
       "        manager = IPython.keyboard_manager;\n",
       "\n",
       "    // Check for shift+enter\n",
       "    if (event.shiftKey && event.which == 13) {\n",
       "        this.canvas_div.blur();\n",
       "        event.shiftKey = false;\n",
       "        // Send a \"J\" for go to next cell\n",
       "        event.which = 74;\n",
       "        event.keyCode = 74;\n",
       "        manager.command_mode();\n",
       "        manager.handle_keydown(event);\n",
       "    }\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
       "    fig.ondownload(fig, null);\n",
       "}\n",
       "\n",
       "\n",
       "mpl.find_output_cell = function(html_output) {\n",
       "    // Return the cell and output element which can be found *uniquely* in the notebook.\n",
       "    // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n",
       "    // IPython event is triggered only after the cells have been serialised, which for\n",
       "    // our purposes (turning an active figure into a static one), is too late.\n",
       "    var cells = IPython.notebook.get_cells();\n",
       "    var ncells = cells.length;\n",
       "    for (var i=0; i<ncells; i++) {\n",
       "        var cell = cells[i];\n",
       "        if (cell.cell_type === 'code'){\n",
       "            for (var j=0; j<cell.output_area.outputs.length; j++) {\n",
       "                var data = cell.output_area.outputs[j];\n",
       "                if (data.data) {\n",
       "                    // IPython >= 3 moved mimebundle to data attribute of output\n",
       "                    data = data.data;\n",
       "                }\n",
       "                if (data['text/html'] == html_output) {\n",
       "                    return [cell, data, j];\n",
       "                }\n",
       "            }\n",
       "        }\n",
       "    }\n",
       "}\n",
       "\n",
       "// Register the function which deals with the matplotlib target/channel.\n",
       "// The kernel may be null if the page has been refreshed.\n",
       "if (IPython.notebook.kernel != null) {\n",
       "    IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n",
       "}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAyAAAAJYCAYAAACadoJwAAAgAElEQVR4XuzdB5QVRdrG8f9kco6SQXLOSBBREFEQFFQUjCAgSjChq9/uGtYcSYIKAgoKCmYEMxIl5ziSMwyZGYaJ3+nukaAgM9N9+/a99+lzPLpD91tVv6rVeqe7qsLS09PT0SUBCUhAAhKQgAQkIAEJSMAFgTAlIC4oqwgJSEACEpCABCQgAQlIwBRQAqKBIAEJSEACEpCABCQgAQm4JqAExDVqFSQBCUhAAhKQgAQkIAEJKAHRGJCABCQgAQlIQAISkIAEXBNQAuIatQqSgAQkIAEJSEACEpCABJSAaAxIQAISkIAEJCABCUhAAq4JKAFxjVoFSUACEpCABCQgAQlIQAJKQDQGJCABCUhAAhKQgAQkIAHXBJSAuEatgiQgAQlIQAISkIAEJCABJSAaAxKQgAQkIAEJSEACEpCAawJKQFyjVkESkIAEJCABCUhAAhKQgBIQjQEJSEACEpCABCQgAQlIwDUBJSCuUasgCUhAAhKQgAQkIAEJSEAJiMaABCQgAQlIQAISkIAEJOCagBIQ16hVkAQkIAEJSEACEpCABCSgBERjQAISkIAEJCABCUhAAhJwTUAJiGvUKkgCEpCABCQgAQlIQAISUAKiMSABCUhAAhKQgAQkIAEJuCagBMQ1ahUkAQlIQAISkIAEJCABCSgB0RiQgAQkIAEJSEACEpCABFwTUALiGrUKkoAEJCABCUhAAhKQgASUgGgMSEACEpCABCQgAQlIQAKuCSgBcY1aBUlAAhKQgAQkIAEJSEACSkA0BiQgAQlIQAISkIAEJCAB1wSUgLhGrYIkIAEJSEACEpCABCQgASUgGgMSkIAEJCABCUhAAhKQgGsCSkBco1ZBEpCABCQgAQlIQAISkIASEI0BCUhAAhKQgAQkIAEJSMA1ASUgrlGrIAlIQAISkIAEJCABCUhACYjGgAQkIAEJSEACEpCABCTgmoASENeoVZAEJCABCUhAAhKQgAQkoAREY0ACEpCABCQgAQlIQAIScE1ACYhr1CpIAhKQgAQkIAEJSEACElACojEgAQlIQAISkIAEJCABCbgmoATENWoVJAEJSEACEpCABCQgAQkoAdEYkIAEJCABCUhAAhKQgARcE1AC4hq1CpKABCQgAQlIQAISkIAElIBoDEhAAhKQgAQkIAEJSEACrgkoAXGNWgVJQAISkIAEJCABCUhAAkpANAYkIAEJSEACEpCABCQgAdcElIC4Rq2CJCABCUhAAhKQgAQkIAElIBoDEpCABCQgAQlIQAISkIBrAkpAXKNWQRKQgAQkIAEJSEACEpCAEhCNAQlIQAISkIAEJCABCUjANQElIK5RqyAJSEACEpCABCQgAQlIQAmIxoAEJCABCUhAAhKQgAQk4JqAEhDXqFWQBCQgAQlIQAISkIAEJKAERGNAAhKQgAQkIAEJSEACEnBNQAmIa9QqSAISkIAEJCABCUhAAhJQAqIxIAEJSEACEpCABCQgAQm4JqAExDVqFSQBCUhAAhKQgAQkIAEJKAHRGJCABCQgAQlIQAISkIAEXBNQAuIatQqSgAQkIAEJSEACEpCABJSAaAxIQAISkIAEJCABCUhAAq4JKAFxjVoFSUACEpCABCQgAQlIQAJKQDQGJCABCUhAAhKQgAQkIAHXBJSAuEatgiQgAQlIQAISkIAEJCABJSAaAxKQgAQkIAEJSEACEpCAawJKQFyjVkESkIAEJCABCUhAAhKQgBIQjQEJSEACEpCABCQgAQlIwDUBJSCuUasgCUhAAhKQgAQkIAEJSEAJiMaABCQgAQlIQAISkIAEJOCagBIQ16hVkAQkIAEJSEACEpCABCSgBERjQAISkIAEJCABCUhAAhJwTUAJiGvUKkgCEpCABCQgAQlIQAISUAKiMSABCUhAAhKQgAQkIAEJuCagBMQ1ahUkAQlIQAISkIAEJCABCSgB0RiQgAQkIAEJSEACEpCABFwTUALiGrUKkoAEJCABCUhAAhKQgASUgGgMSEACEpCABCQgAQlIQAKuCSgBcY360gXFxcXx/fffU758eXLmzHnpB3SHBCQgAQlIQAISkICrAqdOnWLbtm20b9+eIkWKuFp2sBSmBMRDPTlp0iR69uzpoRqpKhKQgAQkIAEJSEACFxKYOHEiPXr0EE42BJSAZAPNV4/MmzePli1bYgzo6tWr+6oYxZWABCQgAQlIQAISyKbA+vXrzV8Yz507lxYtWmQzSmg/pgTEQ/2/bNkyGjZsyNKlS2nQoIGHaqaqSEACEpCABCQgAQkYApqv2R8HSkDsGzoWQQPaMUoFkoAEJCABCUhAAj4R0HzNPqsSEPuGjkXQgHaMUoEkIAEJSEACEpCATwQ0X7PPqgTEvqFjETSgHaNUIAlIQAISkIAEJOATAc3X7LMqAbFv6FgEDWjHKBVIAhKQgAQkIAEJ+ERA8zX7rEpA7Bs6FkED2jFKBZKABCQgAQlIQAI+EdB8zT6rEhD7ho5F0IB2jFKBJCABCUhAAhKQgE8ENF+zz6oExL6hYxE0oB2jVCAJSEACEpCABCTgEwHN1+yzKgGxb+hYBA1oxygVSAISkIAEJCABCfhEQPM1+6xKQOwbOhZBA9oxSgWSgAQkIAEJSEACPhHQfM0+qxIQ+4aORdCAdoxSgSQgAQlIQAISkIBPBDRfs8+qBMS+oWMRNKAdo1QgCUhAAhKQgAQk4BMBzdfssyoBsW/oWAQNaMcoFUgCEpCABCQgAQn4REDzNfusQZeApKam8tprrzFmzBh27txJmTJl6N27N48//jgRERGXFFuzZg1Dhgxh7ty55r0tW7bk1VdfpVatWv/4bOvWrZk9ezY9evRg4sSJlyznQjdoQGeLTQ9JQAISkIAEJCAB1wQ0X7NPHXQJSP/+/Rk1ahT33nsvzZs3Z/78+YwbNw7j5yNHjvxHsdjYWBo3bkyhQoUYMGCAee+wYcM4evQoixYtonLlyhd8/sMPPzTjx8fHKwGxPyYVQQISkIAEJCABCXhWQAmI/a4JqgRk9erV1K1b10wehg4dekZn0KBBDB8+nJUrV1K7du2LqnXr1o2ZM2eyfv16882JcRlvUapXr06HDh347LPP/vaskZxUrVqVhx9+mH/9619KQOyPSUWQgAQkIAEJSEACnhVQAmK/a4IqAXn66ad58cUX2bJlCxUqVDijs3XrVipWrMhTTz3FCy+8cEG1kydPUrhwYbp3786ECRPOu+fuu+9mypQpxMXFkSdPnvP+7MEHH+THH3/E+HQrJiZGCYj9MakIEpCABCQgAQlIwLMCSkDsd01QJSDt27c333Ls27fvbzLFixenfv365huOC10LFiwwP9kyPt/q16/febcYPzM+sTLuadas2Zk/W7p0KU2aNOHbb78135CEhYUpAbE/Jj0XIT09ncTkNBKSUkhISiUxOdX8u/HXqeQUs99jIsPNv6IirL/y5ogkb44ockQZP7/02iPPNVoVkoAEJCABCUjgggJKQOwPjKBKQIzPq6KjozESg79eDRo0IDk5GeMzrQtd06ZNw/gE6+uvv6ZTp07n3WL8rHPnzkydOpWuXbuaf5aWlmYmIyVLluSrr74yf5aVBGTv3r0Yf517GZ9+9ezZ06y/UV9dvhFISklj77FT7D56ij1HE9lz9BT7jyey//hpDsWfJuF0KgnJKZw6k2Skkp6e/brkio6gWN4Y8ueMomDuaEoXzEnJ/DkpmieGcoVzUaFoborkjiE8PCz7hehJCUhAAhKQgARcEVACYp85qBKQSpUqYbzpMBae//Uy3m4cOHCAP/7444JqH330EXfddRfff/8911577Xn3/PDDDxhvV4x7jATBuEaPHm2u+1i3bt2Zz72ykoA888wzPPvssxesixIQ+wPbiBB38jQb9p5gw77jbD54ku2HEsy/jOQjzUZC4Uztzo8SHRFuJiPVS+ajbKFc1CqVn8rF81ChcG4lJr4AV0wJSEACEpBANgWUgGQT7pzHgioBcesNyMGDB82F5w899BDPPffcGc6sJCB6A2J/8J4b4eTpFJbvOMKSbUdYtuMI6/ceJ+5kUpYKMd5QFM0bQ+6YSHJFRWC8ucgZnfH3KOOfI82f/fnznBn3GG9HklLTzE+zUtLSMd6wHD+VQnyS9RbFeKsSdyKJ44nJHDhxmsPxma+XUUbNy/JRp3QB6pbJT4OyBSlTKFeW2qWbJSABCUhAAhJwTkAJiH3LoEpA3FoDYuyy9fHHHzNr1ixy5sx5pheMbXqNT7Vef/11ihQpQoECBbLUQxrQmeMy1mQYn08t3W4lHMbfjbccl3qrUTxfDOUK5aZs4VzmW4ZSBXJymflXDorny0GOKHfWasSfTuHgidPsO57I1rh4th2KNz8Fi91/gj8OnDSTmH+6jHpfUakwLS8vQsvKRSiSJyZzcLpLAhKQgAQkIAHbApqv2SYkqBIQY5erl156ySe7YE2ePJlDhw6Zu2B16dLlzLqPi3WBcRjiY489lqUe0oC+MFdqWjrr9hxnyfbDLNl+hKXbjpiT94tdeWMizU+YqpXMS7USxl/5zM+ZckVHZqk//HGz0dadhxNYvfsYm/afYO2e4+Y/GwnLxa4aJfPRumpRrqpSlPplCxIdGe6PqqtMCUhAAhKQQEgIaL5mv5uDKgExdsAydrq62DkgK1asoE6dOuZi9M2bN5M/f35zEfmfl7HA3FgDsmHDBkqXLm3++M9zQIy3K8ZCdeMydsPavXv33/RvueUWWrVqxcCBA81yqlSpkqUe0oA+y2V8pjR700F+3XjA/PuRhOSLWpYplJNG5QrRsFxB868qxfMSEUQLuo03PkbCtWz7URZtPcT8zYeIPXDygh7GZ2StqxSlXY3iXF2tmPk5mS4JSEACEpCABJwT0HzNvmVQJSAGh7GF7rvvvmuehN6iRQvmzZtnnoTet29fc+G4cW3bts1cOG6c7zF+/Pgzihs3bjS31TXOAzGSCOMyTkI33nwYJ6Eb6z7+6crKGpALxQnlAZ2Wlm7+pt9IOGZtPMjKXUcvuPNUZHgYNUvlp2HZgjQqX5BG5QpSLF8O+/9PCLAIxq5d8/6IY05sHHP/iLvgGxJjYbvxZqRzvctoW724a5+YBRilqisBCUhAAhcTSIqH/etg3yrYvxZOHYG0FOvu8Ehj+0/IWxJyFYIC5aBIZShcGaKDe61iKM/XnPo/S9AlICkpKbz66quMGTOGXbt2mW8yevfuzZAhQ4iMtH4bfLEExPizVatWmfcaiYtxtWzZkldeecV8o3GpSwnIpYTO//MjxluO2IP8tvEgv206yKELLM42/t1Wt3QBrqpalGYVC5v/bCwM13VWwEje1u87zi/rD/DzhgOs2Hn0bzzGuSQd61xG1walzLdExljVJQEJSEACEjhPID4O1n0F2+bAvtVwaDOQxW0jw6OgaDUo2xQqXAmVroaYvEEFrQTEfncGXQJin8R/EUJhQBtb4M5cs48Zq/eZazoutN66YK4orqxSlDZVi9GqchEKa5F1lgblgROJfL9mH9NX72XR1r8b1y6Vn1sblaZbwzJK5rIkq5slIAEJBKHAqaOwYTqsmQZbZkF66sUbGZU7480HkJpi3Zty8TWZZqCwcKh4FdS8Gap3hJwFAx4xFOZrvu4kJSC+Fs5C/GAd0LuOJJhJx3er97Jsx99/O28Q1S2dn9ZVi9GmalFzy9lgWsORhSHg+K3Gp1pfr9jDtGW72LDvxHnx88REcnuTMtzbooK5G5guCUhAAhIIIQHj7cb84bDiY0j9y0YnEdHWW4wSdaBEbeuv4jUh51929zT2oTc+yzp5AI5sg4PrYc8K2LUYjv99rSyROaF6J2jc23pDEqBXsM7X3OwOJSBual+irGAa0DsOJfDdmr3MWL2XlbuO/a3lxlkabaoV45pqxcy3HdpK1rcD0VjIvn7vCSbM38Y3q/aQkHT2N1zG11id6lxGnysrmruH6ZKABCQggSAWMNZ0zH7V+tQqPe1sQ41Ppyq3g1pdocp1EJPHHsLhrbD5Z9g4A7bOhtS/nIFVqiG0GAzVOkJ4YO3eGEzzNXudnP2nlYBk387xJwN9QG85eJIZGW86jO1j/3oZv3FvW70YHWqXNHdqcuvcDcc7KsADGgcifrZkF+Pnb2Xn4VPnteaKioV54KpK5qdvWicS4B2t6ktAAhI4VyDhMPz6IiwZe37iUbwWNOkDNW703edR8Ydg3Rew7EPYu/L8filaHdr8C6p1CphEJNDna174P4YSEC/0QkYdAnFAbz54kumr9pqfV/31Ex+jWflyRNKuRgmur13CPDQvJlILyL0y5IwzR35ct5/3Zm/+26dxtUrl46E2lbm2RnHCg2hLY6/Yqx4SkIAEXBMwPpNaORm+/5f1udSfV7mW0HIwXN7W2s3Kjcuoy64l8PtIWPvl+QvcL2sA7Z6DCq3cqImtMgJxvmarwT54WAmID1CzGzJQBnTcydN8s3IPXyzfzaoLfF5VIFcU7WuUoEPtEjSvVEQH42V3QLj43NLth3n3ty38tH7/eRsD1Cmdn4fbVTEPOdQbERc7REVJQAIScELAeOvxzSBY//XZaEWqwHUvw+XXOFFC9mMc3gKz34CVH5//Rsb4JOva/0GhCtmP7eMnA2W+5mMGW+GVgNjic/ZhLw/oE4nJ5m/LjcRjdmwcxm/Pz70K546mfa0SXF+rJE0rFiIqIrC+53S2JwM3Wuz+E4z+bQtfrth9Xh83Ll+QJ66rRqPyhQK3caq5BCQggVAS+ONn+LI/nNxntdpYAH7109C0H0REeUci7g/45TlrTcqfl1HXVo9Cy4chwnsH6np5vuadjv3nmigB8VBPeW1AJ6WkMWvjAXMy+tP6Axj/+9zLWNPRoVYJbqpfiqYVC2vnKg+NJbtV2X30FG/+sInPl+8670BI44T1Ie2rUrl4cO3pbtdLz0tAAhLwjEDyKfjxv7Do3bNVKlkPuo6xDgr06rXjd5j5JOxZfraGRarCTaOhVANP1dpr8zVP4WSyMkpAMgnlxm1eGNDGoXZLth8xkw5jbcexU8nnNd04idzYtapL/VK0q15c50i4MTD8WIaxscAbP2wyzxT584qKCOOOJmV5pF1V8ufy0G/R/OikoiUgAQl4QsDYWnfyHXBwg1Ud4wwO401C6ye89dbjYlhpabBsAvzyPCQcyrgrzGrDlY9DVA5PMHthvuYJCBuVUAJiA8/pR/05oI1tc6cu3cm0Zbsxfvt97mWsTWtWoTCd613GdbVKUCBXtNNNVzyPCxinq7/43XrzYMM/r/w5oxhyXVVub1xWC9U93n+qngQkEAICu5bCx7ecnbgXKAc3vx+Y522c2G8tmjcOR/zzuqw+dBkNxar5vTP9OV/ze+MdqoASEIcgnQjj9oA+mpBkbpv75fLdLDxnYvlnW6qVyMvNDUrRqe5llMyvg+qc6ONAj/HD2n28PGMDW+LizzTFWKj+Qpfa1C6tM0QCvX9VfwlIIEAFNv0An90NyQlWA2reBDcOh5gA/lzW2DFrw7fw9UA4lfHLr+i80OEVqN/Drx3l9nzNr431UeFKQHwEm52wbg7oXzccoM9HS0hOPX8xeYl8OczPq4x1HVVLBPC/uLLTAXomUwKJyanmgYbDf/mDk6dTzjxz9xXl9FlWpgR1kwQkIAEHBVZOgS8fgPSMA2abD4C2zwXMmRqXlIiPs9oX+8PZW+veAZ3ehsiYSz7uixvcnK/5ov5eiKkExAu9kFEHNwf04fgkmrzwEylp6cREhmMsLr6lURlaXl5Ei8k9NCa8XJV9xxLNz7K+XrnnTDVL5s/BfzvV4LpaJb1cddVNAhKQQHAIrPjEmpyT8cvE9i/CFQ8GR9vObYXxNuT3UfDjfyAtY22qcYBhj8+gQBnX2+vmfM31xrlUoBIQl6AzU4zbA/rNHzdRvnAurq1ZAmNHK10SyI7AvD/ieOqL1Ww/lPHqH8wDDJ/vUovi+byxYDA77dIzEpCABDwtcG7yYSw2v+ldqHOrp6tsu3LGOhfjU7NjO61QMfmh+yTXDy90e75m282DAZSAeKhTNKA91BmqSpYEjC2ah/8Sy+jfNp/5rC9fjkj+06kmXRuU0iGGWdLUzRKQgAQuIfDX5MNYbF67W2iwGYcrTr0PtvxqtddIvm54Exre49qJ7pqv2R9qSkDsGzoWQQPaMUoF8pPAhn3HeWLqKlbuOnamBtdUK8bLXetQNK9/vtX1E4WKlYAEJOAbgVBOPv4UTUuFn/4L84efNW79JLT5l2/M/xJV8zX7zEpA7Bs6FkED2jFKBfKjQGpaOuPmbcX4xC8hyVoUabwNefHm2txQu6Tehvixb1S0BCQQ4AJKPs7vwBUfw9cDIC1jQ5Ru46DWzT7vZM3X7BMrAbFv6FgEDWjHKBXIAwKbD54034YYB1v+eXWpdxnP3FhTZ8l4oH9UBQlIIMAE1n8DU+60Fpwbnx2F0mdX/9RVW+fARzdBidpw+yeQt4TPO1bzNfvESkDsGzoWQQPaMUoF8ohASmoaY+du5bXvN5o7rhlXmUI5ea1bXZpVLOyRWqoaEpCABDwusHcVfNDeOudDycffO2vrbCjVCKJzudKRmq/ZZ1YCYt/QsQga0I5RKpDHBDbuO8GAT5axaf9Js2aR4WH0b3M5g6+prFPUPdZXqo4EJOAxAeNU8DHXnN356frXocn9HqtkaFVH8zX7/a0ExL6hYxE0oB2jVCAPCpxOSTVPUR83b9uZ2jUqV5B3ejSgmLbr9WCPqUoSkIDfBZISYEJH2L3UqkqjXtDxTb9XK9QroPma/RGgBMS+oWMRNKAdo1QgDwv8uuEAj322kkPxSWYt8+eMYmj3elxVtZiHa62qSUACEnBZwDh8z9hudu3nVsEVWkPPaRAR5XJFVNxfBTRfsz8mlIDYN3Qsgga0Y5QK5HGB/ccTeXjKCuZvPnSmpg+2qcTgtlWIigj3eO1VPQlIQAIuCCx8F2YMsQoqUgV6/Qg5C7hQsIq4lIDma5cSuvSfKwG5tJFrd2hAu0atgjwgYGzXO/LXP3j7p01krE+nxeWFeevWevokywP9oypIQAJ+FNi+ACZ0grRkiM4LfWZBkcv9WCEVfa6A5mv2x4MSEPuGjkXQgLZBmXgc4mLh+C5IOASnjkLKaUhPs3bFyFnI2povf2nrN0l6hW0D29lH5/8Rx6ApKzh44rQZuFSBnAy7vR4NyxVytiBFk4AEJBAIAsd2w6jmkHjUqu0tE6Bml0CoecjUUfM1+12tBMS+oWMRNKAzSXl8L2yaCTsXwr7VcHCj9VuirFz5y0DxWlC6IZRvBaUbQ3hEViLoXgcF9h1LZOAny1m07bAZNSYynP92qskdTcs6WIpCSUACEvC4QHIijLsO9iy3KtrqMbjm3x6vdOhVT/M1+32uBMS+oWMRNKAvQmksxNs+H4xDmIzE48hWx8zPBIrOA5XbQfUboWoHiMrpfBmK+I8CSSlpvDRj/Xm7ZBkJyHM31iRS60I0eiQggVAQ+OHfMH+Y1dJKV0OPaRCudXFe63rN1+z3iBIQ+4aORdCA/gvl4a1gLMJb/RkkxF3cuWAF67OqguUgXynIVdhaqGckEeGRkJJkvco+sQ+O74bDW+DgBji648Ixo3JBrZuh4X3WGxJdrgp8sXwXT0xbjZGQGFeT8oUYc08j8uXQzi+udoQKk4AE3BXYMgs+ND61SgfjLf0D8yBHfnfroNIyJaD5WqaY/vEmJSD2DR2LoAENpKXCysmwcDTsW3VhW2MrwoqtocJVcFm97H86lRQP2+aB8S/92B/gUOzfy7usPrQYDDU6Q1iYY32tQP8ssGb3Me7/cAl7jyWaN15eLA8f3teEywrozZTGjgQkEIQCx/fAu1dC/EEgDO7+Biq0CsKGBkeTNF+z349KQOwbOhYhpAf06ZOwYATMHwFJJ/5uWq0j1LkVqnSAyGjHzM8LZCxiXz0VVk2GI2cPyzPvMd6wtHna+kRLr8N94/+XqMai9N4TFrNy1zGrC/LEMP7extQqpd8IutIBKkQCEnBHwPjM+OPbIPZ7q7wrh8DVT7tTtkrJlkBIz9eyJfb3h5SAOATpRJiQHNDG7lWzX4MFIyE99XzGkvWgwZ1Qr4e7azKM/xhsmwO/j4aN08+vU6lGcN1LUKaJE12uGJcQiD+dwqDJy/lp/QHzzpxREQy7vT7tahSXnQQkIIHgEDj3vI9yLay3H9oUxdN9G5LzNYd7RAmIw6B2woXUgDbWZcwbCrNfhVTrROwzV53boOUjUKyaHU5nnj20GX57FVZNsb7L/fMy6tj2WchX0plyFOWiAimpaTz/7TomLNhu3mN8Cfd/N9SgV8sKUpOABCQQ2ALGf2OMT6+STkJMfug7CwpVDOw2hUDtQ2q+5qP+VALiI9jshA2ZAb1mGsx4IuNb1wypsHBo0gdaPQp5imWHz7fPHNwEPz1z/huRnAXh2v9B/Z6+LVvRSU9PZ/z8bTz7zbozGve1qMC/O1YnTGtzNEIkIIFAFEhLs7bcNbaUN65OQ6HhPYHYkpCrc8jM13zYs0pAfIib1dBBP6D3roTpj8GuRefTNO4NVz0FuQtnlcz9+7f8Bt89BnGbzpZ9eVu46V3IXcT9+oRYib9s2M8DE5dxOmOHrA61SvB293rEROoMlxAbCmquBAJfwFjz+EPGWo+qN0D3SdrsJEB6Nejnay70gxIQF5AzW0TQDmhjgfnPz8LiMdbJ5H9e1TvBda9A/lKZJfLGfakpMO9ta+1KirVLE8bbkG7joFIbb9QxiGuxetcx7vxgIUcTrMMn21QtyqieDckRpSQkiLtdTZNAcAkY28CPaAIpp6ytdh+YD/lLB1cbg7g1QTtfc7HPlIC4iH2pooJyQG+YDl8PPP8cj2I1oMs7YGxxG8jXgfXwZX/Ys+xsK4ydsozPyLSA0Kc9u/NwAj3HLmT7oQSznPplC7zCUVcAACAASURBVDD+nibkz6WzQnwKr+ASkIB9AWOjk3HXw475VizjDXrd7vbjKoJrAkE5X3NNzypICYjL4P9UXFAN6Pg4K/E4dxepyBzWwm1jrUewbGVrLKb/9QXrjcifl7FlcOcR1lsRXT4T2HcskTve/50tcfFmGbVK5WNir6YUyOWjbZp91hIFloAEQkpg+ST4qr/VZOMT3h5T9elVgA2AoJqv+cleCYif4C9UbNAM6LVfwlcPnX+eR82b4frXA2OdR3bGxKYf4PPekGidWUHx2nDrBChcKTvR9EwmBY4mJNFrwhKWbj9iPlHZOLCwVxNK5teBhZkk1G0SkICbAif2w/CG1n8fjV/KPbQECpRxswYqywGBoJmvOWCR3RBKQLIr54PnAn5AJyXA1wNgzdSzOnmKQ+d3oHJbH4h5LOThLTC5JxxYa1UsTwm4bSKUaeyxigZXdU6eTuGeDxaxJCMJqVAkNxPubULZwrmCq6FqjQQkEPgC0+6H1Z9a7ejwGjTtE/htCsEWBPx8zQN9pgTEA53wZxUCekDvXAyf3gkn9p4VbXCXtcg8OoQmgqdPgPEfmE0zLIeo3HDrh6GRgPnx/0sJSSk8OGkZv248aNaiTKGcZhJSsWgeP9ZKRUtAAhI4R+CPn2BiV+sHJWpDn9+0XjBAB0hAz9c8Yq4ExCMdYVQjIAe0sZhu/jD46dmzJ5nnKGB9flTxKg/puliVtFSY+S9Y9K5VaHgkdB0DNW9ysRKhV1RiciqDJ69g5tp9ZuNLFcjJR72UhITeSFCLJeBBAeMLAePAwUOxYJx71XcOlKjlwYqqSpkRCMj5WmYa5uI9SkBcxL5UUQE3oBOPw+d9zv6232hgxTbQdWzwrvW4VCf++edGYjb3LWv7YfMKg5tGa6eTzPpl877TKak8+ulKvl1lvYkrmjeGqf2uoFzh3NmMqMckIAEJOCDw22vw6/+sQI16Qcc3HQiqEP4SCLj5mr+g/qFcJSAe6pSAGtAHNsDHt4Cxl/mfE+xr/gMtH9ZuHueOqSXjYPqjZ98OdRkN9W730KgLvqokp6aZScjXK/eYjStdMCeT+zSjdMEQ+hQw+LpVLZJA4Aoc2QbvNIfkeGtt4IClEKPPQwO3QwP0ixWPgSsB8VCHBEwCsu4ra51D6mlLL1dha7F1ueYe0vRQVVZOhq8ehLQUq1I3vw91bvVQBYOvKqlp6QyesoJvMpKQ4vmMNyHNKVNISUjw9bZaJAGPC0y9D9ZMsyrZeSTU7+nxCqt6lxIImPnapRrixz9XAuJH/L8W7fkBbXxW9PNzMPecV8dlr4DuH0OuQh6S9GBVVn0KX/Sz3oSERUC3D6BmFw9WNHiqZLwJGfDx8jNrQsoVzsVn/a6gWN4cwdNItUQCEvC2wNY5MKGjVccyzeC+mfpKwNs9lqnaeX6+lqlW+PcmJSD+9T+vdE8PaGN3p0/vgs2/nK3zFQ9Bu+eD51BBX4+FFR9bJ6eTbu3/3n2SdQiVLp8JGG9CBn6ynOmrrTUhFYvk5uP7m1Eiv5IQn6ErsAQkYAkYv7R7vw3sWW6tA+z1A5RpIp0gEPD0fC1AfJWAeKijPDugD26CiTfDsZ2WlrGr003vQu1uHtILkKosHgvTH7EqG5Mf7voSSjUIkMoHZjWNhem9JyxhTmyc2YCqxfPyad8ryJ8rKjAbpFpLQAKBIbB8ovX5rXHV6wldRgZGvVXLSwp4dr52yZp75wYlIN7pC29uw7thOnx279n1HsYCuh6fQsm6HpILsKrMfh1+ed6qdO5icP8vOgnXx11obNHbb+JSZmWcE1K/bAHG39uE/DmVhPiYXuElEJoCyYkwrJ51NpZxHlT/BVCwXGhaBGGrlYDY71QlIPYNHYvgqQFtvDo2Jsp/bhtotNJY73HbJG2xa7fHDdsZT5w9J6Tw5XD/r5Ajn93Iev4fBIzDCnuMWcjyHUfNuxqVK8jE3k3JERUhNwlIQALOCsx6BWa9aMW86im46gln4yuaXwU8NV/zq0T2C1cCkn07x5/0zIBOSYIv+sDaL8620di3vMOrEBHpeLtDMmBqCnx2N2z41mp++VbQ83OIjA5JDrcafTQhibvHLWblTisJaV+zOCPvaEBkRLhbVVA5EpBAsAvEH4K3a0FyAuQpDg8thhz5g73VIdU+z8zXAlhdCYiHOs8TAzrhsLXew1w0l7Hew0g8GvfykFSQVCX5FHzQHvautBrU+H644fUgaZx3m2EkIV1HzWfzwXizkjfUKcmw7vWJCA/zbqVVMwlIIHAEpj8Gi9+36nvjCGhwZ+DUXTXNlIAn5muZqql3b1IC4qG+8fuAjouFj246u9jc+I2N8clVhVYeUgqyqsTHwXtXnTXv/A7U7xFkjfRec3YdSeDW0QvYcyzRrNx9LSrwn041vFdR1UgCEggsgcNbYHhDSE8D4/PaBxbozXZg9WCmauv3+Vqmauntm5SAeKh//Dqgt82DT7rD6eOWSIFycNdXUKiCh4SCtCrG26ax7a2F/sYZIYa7kj6fd/bOwwnc9M584k5aB2oObluZwW2r+LxcFSABCQSxwNResGaq1cA7PoUq7YO4saHbNL/O14KEXQmIhzrSbwPaOCTPOJ8iLdnSKNMUekzVomg3x8bqqTAt4zM345vhfnMhTzE3axCSZa3edYzb3ltAQlKq2f5nb6zJ3c3Lh6SFGi0BCdgUMH6ZZLzRNq7LGlg7HIbp006bqp583G/zNU9qZK9SSkCy5+aTp1wf0MZuTHPeOLslrNGq2rdYZ3yEa2cgn3TyPwU1Tpk3+sO4KlxpLUqP0Daxvu6HRVsPc9cHC0lMTjPnCqN6NOS6WiV8XaziS0ACwSRg/Pf041sh9gerVb1+gjKNg6mFass5Aq7P14JQXwmIhzrV9QH9ywsw+9WzAlcOgTZP6Tc2/hoTxs5YEzrBjvlWDZr1h+te8ldtQqrcH9ftp+9HS0gzDqkPD2PsPY1pXaVoSBmosRKQgA2Bzb9YayiNq+JV1qe0uoJWwPX5WhBKKgHxUKe6OqC3zrEmu6RbAtqpwxsj4cQ+eP8aOL7Lqo/xKVzldt6oW5DXYtLC7Tz9xRqzlbmjI/isX3NqXKazWYK829U8CdgXMN5+jOsAOxYAYdD/dyhWzX5cRfCsgKvzNc8q2KuYEhB7fo4+7eqA3rcapt4Hxs5X5kK5ax1ti4LZENixEMZdZ+2ikrso9J0D+UraCKhHMyvw/uwtvPDdevP2Ynlj+KRPMyoVzZPZx3WfBCQQigJrv7TOdTKuGl3g1gmhqBBSbXZ1vhakskpAPNSxrg/opATY/DNUN96E6PKUwLmn6BqbAtw7Q+tyXOig9PR0/vPVWj76fbtZWplCOfmyfwsK54lxoXQVIQEJBJyA8fbDWHi+dwWER8GApVCwXMA1QxXOmoDr87WsVS8g7lYC4qFu0oD2UGf4uyrGepCPbwHju2LjuuY/0OpRf9cqJMpPTUvn0U9X8OWKPWZ7G5QtwKTezcgZrY0ZQmIAqJESyIrAuTsY1usBXd7JytO6N0AFNF+z33FKQOwbOhZBA9oxyuAIdPIgvNMMEuKs36wZWzqWrBMcbfN4K5JS0rhz7EIWbj1s1rRV5SJ8cE9joiLCPV5zVU8CEnBNwHj78W4rMD5pjswBDy2BAmVcK14F+U9A8zX79kpA7Bs6FkED2jHK4Am0/huY0tNqT9Hq0GcWROUInvZ5uCVHE5K49d0FbNp/0qzlHU3L8kKXWoRpX38P95qqJgEXBc59+9HwHug01MXCVZQ/BTRfs6+vBMS+oWMRNKAdowyuQF88ACs/ttrU9AHo8HJwtc/DrTlwIpEuI+ax51iiWctH2lVh4DWVPVxjVU0CEnBFwHj7Mao5HFgHETHw0GKt/XAF3huFaL5mvx+UgNg3dCyCBrRjlMEVKPEYjGoJx3ZY7brzC6h0dXC10cOtWb/3ON1GzSc+47T0V7rW5rbGZT1cY1VNAhLwucD6b2FKD6uYRr2g45s+L1IFeEdA8zX7faEExL6hYxE0oB2jDL5AOxfBB+2trXkLlocHFkB0ruBrp0dbtGDzIe4etwhjbYhxUOGHvZrQvFIRj9ZW1ZKABHwqYL79aAEH1lrr8x5cCIUr+bRIBfeWgOZr9vtDCYh9Q8ciaEA7RhmcgWY+Bb+PtNrWYjC0ezY42+nRVn29cg8DP1lu1s44qHBK3yuoVSq/R2uraklAAj4T+ONnmHizFb7RfdDxLZ8VpcDeFNB8zX6/KAGxb+hYBA1oxyiDM9Dpk9auWMd2QlgE9P0NStQOzrZ6tFXDf47ljR83mbUrkS8HXz3UguL5tCmAR7tL1ZKA8wLmuR+tYe9KCAuH/guhaBXny1FETwtovma/e5SA2Dd0LIIGtGOUwRso9keY1M1qX8l60PtniIgM3vZ6rGXGQYXPfL2WCQusgwprl8rPpPubki9HlMdqqupIQAI+Edg2D8Zfb4Wu1xO6ZLyV9klhCupVAc3X7PeMEhD7ho5F0IB2jDK4A03tBWumWm1s9zy0GBjc7fVY61JS07h3/GLmxMaZNWt5eRHG39uYSJ0R4rGeUnUk4LCA8fZj/A2wfZ4V+IH5ULymw4UoXCAIaL5mv5eUgNg3dCyCBrRjlMEdyDigcGQTOHUYInNC//lQqGJwt9ljrTuemEzPMQtZteuYWbN7mpfnmRs1EfFYN6k6EnBWYNcSGHONFbP6jXDbR87GV7SAEdB8zX5XKQGxb+hYBA1oxyiDP9DKKfBFH6udFa+CO78EHZDnar8bZ4R0HjGPvRlnhDxxXTUeuEo74bjaCSpMAm4KTO4BG761SuzzG1xWz83SVZaHBDRfs98ZSkDsGzoWQQPaMcrgD2R8CmCsBfnjJ6utt0yAml2Cv90ea+Ga3ce4ZfQCTiWnmjUb3bMB19Uq6bFaqjoSkIBtgX1rYHQLK0z5VnBPRiJiO7ACBKKA5mv2e00JiH1DxyJoQDtGGRqBDm22dsVKTYJ8peDBRRCTJzTa7qFWzo2NM88ISU1LJ2dUBFMfuIKal2l7Xg91kaoiAfsCXw+AZR9ace6dCeWusB9TEQJWQPM1+12nBMS+oWMRNKAdowydQL/8D2a/ZrW3xSBo91zotN1DLf1wwTb+89Vas0aFc0fz5YMtKFNIB0V6qItUFQlkX+DwVhjeENJTra3P+87RJ6/Z1wyKJzVfs9+NSkDsGzoWQQPaMcrQCZSUACObwrEdEB4J/eZCseqh036PtPSv2/PWKJmPaQ80J2d0hEdqqGpIQALZFvjh/2D+cOvx7p9AtYxteLMdUA8GuoDma/Z7UAmIfUPHImhAO0YZWoE2zoBPulttLtcC7pmu3875YQSkpaXTd+JSfly33yz9upolGHFHfW3P64e+UJEScEwg/hAMqwenj1u7DT60FMLDHQuvQIEpoPma/X5TAmLf0LEIGtCOUYZeoE/ugI3TrXZ3HQu1Mw4rDD0Jv7b4RGIyN70znz8OnDTrcX+rCjx9Qw2/1kmFS0ACNgRmvQKzXrQCdHwLGt1nI5geDRYBzdfs96QSEPuGjkXQgHaMMvQCHdkOIxpD6mnIVxoGLIGonKHn4IEW7zycQNdR8zlw4rRZm+G316dT3cs8UDNVQQISyJJAymkYVh+O74Y8JeDhNRARlaUQujk4BTRfs9+vSkDsGzoWQQPaMcrQDPTz8zDndavtbZ6G1kNC08EDrV66/Qjd31tAcmo6OaLCmdirKY3KF/JAzVQFCUgg0wKLx8L0R6zbWz8BbZ7K9KO6MbgFNF+z379KQOwbOhZBA9oxytAMdPqktVPLyX0QlQsGLIV8+s27vwbDxwt38NQXq83ii+SJ5tsBrSiRP4e/qqNyJSCBrAgYZy2NagEH1kJ0Xnh4NeQsmJUIujeIBTRfs9+5QZWApKam8tprrzFmzBh27txJmTJl6N27N48//jgREZfejWbNmjUMGTKEuXPnmrItW7bk1VdfpVatWmekExIS+PDDD/n6669ZvXo1hw4donz58nTs2JGnnnqKAgUKZLtXNKCzTacH/xRY8TF8+YD1v+rcBje/Jxs/Crw8YwOjf9ts1qBB2QJM7nMF0ZFawOrHLlHREsicQOxPMKmrdW+jXtDxzcw9p7tCQkDzNfvdHFQJSP/+/Rk1ahT33nsvzZs3Z/78+YwbNw7j5yNHjvxHrdjYWBo3bkyhQoUYMGCAee+wYcM4evQoixYtonLlyubPjCSlTp06tGrVivbt21OsWDGWLl1qJj1GImL8c758+bLVMxrQ2WLTQ+cKpKXBmGtgzzLrp71/htKNZOQnAeNwwnvGLWJObJxZA2MtyLDu9QgLC/NTjVSsBCSQKYEJnWDrbAgLtw55LWLNAXRJwBDQfM3+OAiaBMR4G1G3bl0zeRg6dOgZmUGDBjF8+HBWrlxJ7dq1LyrWrVs3Zs6cyfr16803J8ZlvEWpXr06HTp04LPPPjN/FhcXx549e8wk5Nzrgw8+oFevXrzxxhs88kjGN6NZ7B8N6CyC6fYLC+xYCB9ca/1Z6cbQ60dty+vHsXIkPonOI+ex43CCWYv/u6E6vVtV9GONVLQEJPCPAnGxMCLjFzc1OsOtGSegi00CGQKar9kfCkGTgDz99NO8+OKLbNmyhQoVKpyR2bp1KxUrVjQ/j3rhhRcuKHby5EkKFy5M9+7dmTBhwnn33H333UyZMsVMPPLkyXNR8ePHj5M/f37uu+8+xo4dm62e0YDOFpseupDA1PtgzTTrT24eA3VukZMfBbbFxdNpxFxOJKYQER5mLkq/olJhP9ZIRUtAAhcV+LwvrJps/fGdX0Clq4UlgfMENF+zPyCCJgExPocy3nLs27fvbyrFixenfv365huOC10LFiwwP9kyPt/q16/febcYPzM+4TLuadas2UXFN27cSLVq1XjyySd56aWXstUzGtDZYtNDFxI4usPaljclEQqWh4eWaPtIP4+UH9buo89HS81a5I2J5KuHWlCx6MV/qeHn6qp4CYSmQMJheKMqpCZBqUbQ+ye9QQ7NkfCPrdZ8zf6gCJoExPi8Kjo62lyD8derQYMGJCcnm4vGL3RNmzYN4xMsY2F5p06dzrvF+Fnnzp2ZOnUqXbtmLEi7QJC77rqLiRMnsnz5cvNTsEtde/fuxfjr3Mv4/Ktnz55mG4w665KALYEf/g3zh1khOr4Nje61FU4P2xcY+lMsb/20yQxUtXhevniwObmiI+0HVgQJSMAZgdmvwS//s2J1GQ31bncmrqIElYASEPvdGTQJSKVKlTDedBgLz/96GW83Dhw4wB9//HFBsY8++ggjgfj++++59tqMb+cz7vzhhx/MxebGPUZycKHrvffeo2/fvubaD2MNSGauZ555hmefffaCtyoByYyg7rmkgPGbvLfrQNIJyFcKBiyDKG0De0k3H96Qnp7OgE+W8+0q65cPxqL0obfVIzxci9J9yK7QEsicgHHw4JvVIeEQ5CkOA1dAdK7MPau7QkpACYj97g6aBMRfb0C+/PJL8+3J9ddfz+eff05kZOZ+m6k3IPYHryJkQuDXl+C3l60br3sZmmVs0ZuJR3WLbwTiT6dw44i5bD4YbxYwuG1lBret4pvCFFUCEsi8wJrPYWrGm+K2z0DLhzP/rO4MKQElIPa7O2gSEH+sATHejtx4443m+pHvvvuOHDns/XZZA9r+gFaEvwgkHrPegiQehdxFYdBKiM4tJj8LbI2Lp+uo+RyOT8LYkffD+5rQqnJRP9dKxUsghAWMgwffvRL2rYLwSHh4LeQtEcIgavo/CWi+Zn98BE0CYuxyZSz+9sUuWJMnTzYPHDx3F6zffvvN3J7XePPy888//+MOWZntJg3ozErpviwJzH0LfnrGeuSa/0Kr7G0TnaUydfMlBebGxnHnBwsx5j0Fc0Xxef8WVCii5PCScLpBAr4Q2LkIxrazIte/EzqP8EUpihkkApqv2e/IoElAjB2wjJ2uLnYOyIoVK8yzO4zF6Js3bza3zC1ZsuQZQWOBubEGZMOGDZQuXdr8+Z/ngBhvV4yF6n9eCxcupG3btuZ2v0YiUrBgQfs9oYNtHDFUkAsIJMXD0HoQfwByFLDeguQsICoPCJy7KP3yYnn46sEW5I7J3GecHqi+qiCB4BH4vA+smmK1p99cKHHxc8OCp9FqSXYFlIBkV+7sc0GTgJj/zujXj3fffdc8Cb1FixbMmzfPPAndWCA+evRos9Xbtm0zEwfjfI/x48efkTC20W3SpIl5HsjAgQPNnxsnoRtvPoyT0KtWrWr+bPv27Waic+LECV5++WVz4fu5l/G/27XL+C1KFvtHAzqLYLo98wK/j4aZT1j3t34C2jyV+Wd1p88E0tLSGTRlBd+s3GOWcXP9Urxxa12dlO4zcQWWwAUEjmyHYfUgPc06vNXYeleXBP5BQPM1+8MjqBKQlJQUXn31VcaMGcOuXbvMNxm9e/dmyJAhZxaHXywBMShXrVpl3mskLsbVsmVLXnnllfNOPZ81axZt2rS5qHzr1q0x7snOpQGdHTU9kymB5EQY3gCO74boPDBoFeTWQXiZsvPxTQlJKXQeMY/YAyfNknRSuo/BFV4CfxWY9TLMyji/6/bJULWDjCTwjwKar9kfIEGVgNjn8G8EDWj/+gd96UvHwzeDrGY2HwjXPh/0TQ6UBsbuP0HnkfNISErF2JF3Yu+mNK9UJFCqr3pKIHAFkhJgaF3rE9U8JazF5xH6DDJwO9Sdmmu+Zt9ZCYh9Q8ciaEA7RqlAFxJITbZORz+yFSJzwqAV2uXFQyPl1w0HuG/CYnNRepE8MXw3qCXF8trbWc9DzVNVJOBNgWUfwdcPWXW75j/Q6lFv1lO18pSA5mv2u0MJiH1DxyJoQDtGqUAXE1g5Bb7oY/1pkz5w/Wuy8pDA699vZMSv1oGp9csW4JP7m5EjKsJDNVRVJBBkAmOvhZ0LITIHPLpRG3QEWff6qjmar9mXVQJi39CxCBrQjlEq0MUE0lJhVHM4uAHCo2DgMihQVl4eEUhJTePOsYtYsOWQWaO7rijHc51reaR2qoYEgkxg+wIYd53VqNq3QNcxQdZANcdXApqv2ZdVAmLf0LEIGtCOUSrQPwms+wo+vcu6Q/vde26sHEtIpuOIOew8fMqs2zs9GnB97bNbhnuuwqqQBAJV4PO+sGqyVft+86CEkv1A7Uq36635mn1xJSD2DR2LoAHtGKUC/ZNAWhq819o68TcsAh5aDIUrycxDAqt2HTVPSk9OTSdXdASf9buCmpfl91ANVRUJBLhA/CF4uxYkJ0DpJtD7xwBvkKrvpoDma/a1lYDYN3Qsgga0Y5QKdCmBTT/Ax7dYd9W9A24adakn9OcuC0yYv43/fr3WLNU4If2bAS3Jo0MKXe4FFRe0Ar++BL+9bDWv8ztQv0fQNlUNc15A8zX7pkpA7Bs6FkED2jFKBbqUgLHV0phrYPdSCI+EgSugQJlLPaU/d1ngyWmrmLx4p1mqDil0GV/FBa+A8e+/oXXg6A7IWxIGr4aIqOBtr1rmuIDma/ZJlYDYN3Qsgga0Y5QKlBmB9d/AlJ7WnU0fgA4Zvw3MzLO6xxWBU0mp3Dhi7plDCv/dsQa9WlZwpWwVIoGgFdg4Az7pbjVPZyIFbTf7smGar9nXVQJi39CxCBrQjlEqUGYEjLUgI5vAoViIygWD1+h09My4uXzPxn0nuPmdecQnpRIZHsbUB5pTr0wBl2uh4iQQRALjO8K2OdYaOGMnwILlg6hxaoobApqv2VdWAmLf0LEIGtCOUSpQZgWWT4SvHrTubv0ktPlXZp/UfS4KzFyzj34Tl5olli2Ui28HtiRfDn0y4mIXqKhgETi6E96uDaRDjc5w64fB0jK1w0UBzdfsYysBsW/oWAQNaMcoFSizAilJMLQunNgDOQtab0Fi8mT2ad3nosC/v1zDR79vt3LFKkUZd09jwsPDXKyBipJAEAjM/Bf8/o7VkO6fQLXrg6BRaoLbApqv2RdXAmLf0LEIGtCOUSpQVgQWjITvn7KeaP8SXNE/K0/rXpcEEpNT6TZ6Pmt2HzdLfLJDNfq11vbJLvGrmGAQSE6E16vA6WNQuDI8uBDCI4KhZWqDywKar9kHVwJi39CxCBrQjlEqUFYETp+Et2pC4lHIV8raESsyOisRdK9LAnuPnaLD0DkcTUjWehCXzFVMEAmsmQZT77MadO0L0PyhIGqcmuKmgOZr9rWVgNg3dCyCBrRjlAqUVYFfX4TfXrGe0p74WdVz9f7v1+6j70fWepDi+WL4dkAriuaNcbUOKkwCASnw3lWwZzmER8HDayBviYBshirtfwHN1+z3gRIQ+4aORdCAdoxSgbIqYJwKbLwFSTkFRapAf+PThPCsRtH9Lgk8+81axs3bZpam9SAuoauYwBbYvxZGNbfaUK8ndBkZ2O1R7f0qoPmafX4lIPYNHYugAe0YpQJlR2DGE7BwtPXkbZOgesfsRNEzLgikpKZx+/u/s3jbEbM0nQ/iArqKCGyBGU/CwlFWG+6dCeWuCOz2qPZ+FdB8zT6/EhD7ho5F0IB2jFKBsiNgnAo8rD6kpUCphtD7ZwjTLkvZoXTjmV1HEsz1ICcSU8z1IBN7N6VZxcJuFK0yJBBYAonH4I1qkJyQsfh8kd7wBlYPeq62mq/Z7xIlIPYNHYugAe0YpQJlV+CLfrDyE+vpu7+FCq2yG0nPuSAwc81e+k1cZpZUIl8OZg5uRYFc2kDABXoVEUgCyyfBVxm7+3V8GxrdG0i1V109KKD5mv1OUQJi39CxCBrQjlEqUHYFDmyAd5paT1e6Bu78PLuR9JxLAi/NWM+7v20xS7u+dglG3tGAML25cklfxXheID0dxrSF3UsgIhoei4WcBTxfbVXQ2wKar9nvHyUg9g0di6AB7RilAtkR+OQO2DjditB3NpSsayeanvWxQFJKGjePmnfmfJCnMSfvtAAAIABJREFUrq9Gnyt1PoiP2RU+UAR2Loaxba3a1r0dbspY5xYo9Vc9PSmg+Zr9blECYt/QsQga0I5RKpAdgXP/g13zZrhlnJ1oetYFgc0HT9JlxDxOnLbWg3zevzl1Suu3vC7QqwivC3w3BBa9a9Wyz29wWT2v11j1CwABzdfsd5ISEPuGjkXQgHaMUoHsCoy7AbbPhbBweGgJFNZv1O2S+vr5H9ft5/4Pl5jFVCiSm+kDW5IrOtLXxSq+BLwrYCw+f6sWnD5uLT5/aLE21vBubwVUzTRfs99dSkDsGzoWQQPaMUoFsisQ+xNM6mpFaXgvdHrbbkQ974LAU1+s5uOFO8ySbqhdkhF31Nd6EBfcVYRHBRaPgemPWpW7/nVocr9HK6pqBZqA5mv2e0wJiH1DxyJoQDtGqUB2BYyFm6Nbwf7V1sLNwat1arBdUxeeT0hKofOIecQeOGmW9tLNtbm9SVkXSlYREvCgwAcdYMd8iMwJj/8BMXk8WElVKRAFNF+z32tKQOwbOhZBA9oxSgVyQmD1VJjWy4rUYjC0e9aJqIrhY4FtcfHcMGwO8Ump5IyK4LtBrcxPsnRJIKQE9qyA91pbTa7VFbp9EFLNV2N9K6D5mn1fJSD2DR2LoAHtGKUCOSGQmgIjGsKRbRCTDx5eAznyOxFZMXws8OninQyZtsospVqJvOaidK0H8TG6wntL4NuHYUlG0mEcqlq6kbfqp9oEtIDma/a7TwmIfUPHImhAO0apQE4JLB4L0x+xol3zH2iV8T21U/EVxycC6enp9J+0jBlr9pnxezYry/+61PZJWQoqAc8JpJyGN6tDwiEoVgP6L/BcFVWhwBbQfM1+/ykBsW/oWAQNaMcoFcgpgeREeLs2xB+AvCVh0CqI1EnbTvH6Mk786RQ6DZ/Llrh4s5gP7mnE1dWK+7JIxZaANwSWToBvBmb84uS/0CrjlyjeqJ1qEQQCmq/Z70QlIPYNHYugAe0YpQI5KTDrFZj1ohXx5jFQ5xYnoyuWDwVW7DxK11HzSU1Lp1DuaL56sAVlCuXyYYkKLQEPCIy9FnYuhKhc8PBayFXIA5VSFYJJQPM1+72pBMS+oWMRNKAdo1QgJwVOHoS3akBqElzWAO7/RXvpO+nr41jDfo7lzR83maU0KleQKX2vICI8zMelKrwE/CRweAsMq28VXvtW6Pq+nyqiYoNZQPM1+72rBMS+oWMRNKAdo1QgpwW+7A8rJllRe/0IZZo4XYLi+UggLS2dXhMW8+vGg2YJT1xXjQeu0sGSPuJWWH8LTH8MFmckHT2mQeW2/q6Ryg9CAc3X7HeqEhD7ho5F0IB2jFKBnBbYuwrebWVFrXkT3DLe6RIUz4cCB04k0v6t2RxJSCYqIozJfa6gYbmCPixRoSXgB4HUZHi9Mpw6AkWrQf/f9bbWD90QCkVqvma/l5WA2Dd0LIIGtGOUCuQLgfEdYdscCIuAwasgf2lflKKYPhKYuWYv/SYuM6OXKpCTGYNbkS9HlI9KU1gJ+EFg/bcwpYdVsHbt80MHhE6Rmq/Z72slIPYNHYugAe0YpQL5QuDc/7jrYEJfCPs85nPfrOODeVvNcro1LM3rt9T1eZkqQAKuCYy7AbbPtX5JMmglFCjjWtEqKLQENF+z399KQOwbOhZBA9oxSgXyhUBaqrW48+h2yFEAHlkP0dpRyRfUvoqZmJzKjSPmsmn/SbOIEXfUp2Ody3xVnOJKwD2BY7vgrZpWefpM1D33EC1J8zX7Ha8ExL6hYxE0oB2jVCBfCSwYCd8/ZUXv+BY0us9XJSmujwTW7D5Gl5HzSElLJ1d0BDMHXUnZwkokfcStsG4JzH4dfnneKq37J1DterdKVjkhKKD5mv1OVwJi39CxCBrQjlEqkK8EEo/BmzUg6SQUqQoPLtQiT19Z+zDuhwu28Z+v1polNClfiE/6NNPWvD70VmgfC6SmWG8/Tu6D3MXg4TUQGePjQhU+lAU0X7Pf+0pA7Bs6FkED2jFKBfKlwIwnYOFoq4Sen8Pl1/iyNMX2gUB6urE17xJ+2XDAjP5kh2r0a62teX1ArZBuCGz5DT680Sqp9ZPQ5l9ulKoyQlhA8zX7na8ExL6hYxE0oB2jVCBfChzaDMMbAulweTvoOdWXpSm2jwQOHE/k2rdnczQhGeNcQuOAwsbldWK0j7gV1pcCn/eFVZOtEh5cDEWr+LI0xZYAmq/ZHwRKQOwbOhZBA9oxSgXytcDH3WHTDP0H39fOPo7/8/r95psQ4ypbKBczBrUid0ykj0tVeAk4KBB/CN6oCmnJULox9P7JweAKJYELC2i+Zn9kKAGxb+hYBA1oxygVyNcCW2bBh52tUpr0hetf9XWJiu8jgWe+Xsv4+dvM6D2aluWFm2r7qCSFlYAPBJaOh28GWYFveg/q3uaDQhRSAucLaL5mf0QoAbFv6FgEDWjHKBXI1wLp6TCyKcRthJh81pa8MXl8Xari+0DgVFIq1w+bw9a4eDP68Nvr06mutub1AbVC+kJg7LWwcyFExMCQzRCT1xelKKYEzhPQfM3+gFACYt/QsQga0I5RKpAbAoveh+8es0q64U1o3MuNUlWGDwRW7DzKraMXkJSaRv6cUfzw8JUUz5fDByUppAQcFNi3Gka3tALWvhW6vu9gcIWSwMUFNF+zPzqUgNg3dCyCBrRjlArkhkDicXizurUlb9Hq0H+BtuR1w91HZYyZs4X/TV9vRr+6WjHG3t2IsLAwH5WmsBJwQOCnZ2Hum1ag+76Hss0cCKoQEri0gOZrlza61B1KQC4l5OKfa0C7iK2inBGY/igsHmPFumc6lM/4baQz0RXFRYHUtHS6v7eAxduOmKX+3w3V6d2qoos1UFESyIKAcfbHsHpwbCfkKwWD10B4eBYC6FYJZF9A87Xs2/35pBIQ+4aORdCAdoxSgdwSOLAB3mlqlVajC9w6wa2SVY4PBLYfiqfjsLmcOJ1CVEQY3w5oRdUS+qbeB9QKaVdgw3SYfIcVpdVjcM2/7UbU8xLItIDma5mmuuiNSkDsGzoWQQPaMUoFclNgfEfYNgfCI2HwasinBcxu8jtd1ner99J/0jIzbO1S+fm8f3OiIvSbZaedFc+mwNResCbjDKKH10H+UjYD6nEJZF5A87XMW13sTiUg9g0di6AB7RilArkpsO4r+PQuq8TWT0Cbp9wsXWX5QGDAJ8v5ZuUeM7JxQrpxUrouCXhG4NRReLMGJMdDmWbQ63vPVE0VCQ0Bzdfs97MSEPuGjkXQgHaMUoHcFDC+xX67NpzYA3mKW99iR0a7WQOV5bDAkfgk85T0gydOm5E/66dT0h0mVjg7Ar+PhplPWBE6DYWG99iJpmclkGUBzdeyTPa3B5SA2Dd0LIIGtGOUCuS2wG+vwq8vWKV2+wBqdXW7BirPYYHftxyi+3u/m1FL5s/BzEFXkj9XlMOlKJwEsiEwtj3s/B2i88JjmyA6VzaC6BEJZF9A87Xs2/35pBIQ+4aORdCAdoxSgdwWOLEf3qoJaclQtjncN8PtGqg8Hwg89806Ppi31YzctUFpXr+ljrbm9YGzQmZB4NBmGN7AeqBWN+g2NgsP61YJOCOg+Zp9RyUg9g0di6AB7RilAvlD4NxFof3mQYla/qiFynRQICkljRuGzSH2wEkz6og76tOxjjYZcJBYobIq8P3TsGCE9VSPaVC5bVYj6H4J2BbQfM02IUpA7Bs6FkED2jFKBfKHwI7f4YP2VsnGN9nGt9m6Al5g/d7jdB45DyMZKZQ7mm8GtKRUgZwB3y41IAAF0tLg7VpwfDcUrAADlunsjwDsxmCosuZr9ntRCYh9Q8ciaEA7RqlA/hBIT4fRrWD/aojKBY+sh5wF/FETlemwwMhf/+C17zeaUdtULcoH9zTWp1gOGytcJgS2zYXxN1g3tnwY2j6TiYd0iwScF9B8zb6pEhD7ho5F0IB2jFKB/CWwdAJ8M9Aq/bqXodkD/qqJynVQIC0tnR5jFrJgyyEzqk5JdxBXoTIv8OndsO5L6/7+v0Ox6pl/VndKwEEBzdfsYyoBsW/oWAQNaMcoFchfAkkJ8GY1SDwGRarAg4sgLMxftVG5DgrsOJRAx+FzOJ6YQkxkON8NakWlonkcLEGhJPAPAknx8GolSDkF5VrCvdPFJQG/CWi+Zp9eCYh9Q8ciaEA7RqlA/hT4bggseteqwT3fQfkW/qyNynZQ4KsVuxk0eYUZsWG5gnza9woiwpVgOkisUBcTWPUpfH6/9afXvw5NMv5ZYhLwg4Dma/bRlYDYN3Qsgga0Y5QK5E+BA+vhnWZWDWrfAl3H+LM2KttBgfT0dPp+tJQf1u03ow64+nIevbaqgyUolAQuIjCmHexaBBHR8PA6yFNUVBLwm4Dma/bplYDYN3Qsgga0Y5QK5G+BD66DHQusycIjGyB3YX/XSOU7JHDgRCLXD51D3MkkjJcf0x5oTv2yBR2KrjASuIDA0R3wdm3rD+rcBje/JyYJ+FVA8zX7/EpA7Bs6FkED2jFKBfK3wMop8EUfqxbtnocWGQvT/V0vle+IwJzYg9w5dpEZq3KxPHw7sCUxkRGOxFYQCfxNYO5b8FPGjle3T4aqHYQkAb8KaL5mn18JiH1DxyJoQDtGqUD+FkhOtBajnzoChSrBgKVajO7vPnG4/CenrWLy4p1m1HtblOe/nWo6XILCSQAwzv4YWgeO7YRcha3Pr6JyiEYCfhXQfM0+vxIQ+4aORdCAdoxSgbwgcO6JxXd9DRVbe6FWqoNDAscTk7nurdnsOZZoRvz4/qY0r1TEoegKI4EMgXMPONXZHxoWHhHQfM1+RygBsW/oWAQNaMcoFcgLAnGxMKKRVZMaXeDWCV6olergoMCyHUfoOmo+xhmUJfPnYMagVhTIFe1gCQoV8gLfPQ6LMtZ89J0DJeuEPIkA/C+g+Zr9PlACYt/QsQga0I5RKpBXBMZ3hG1zIDzSOhk9TzGv1Ez1cEjghenreH/OVjPaLQ1L89otdR2KrDAhL3D6BLxRDZJOQpGq8OBCfcoZ8oPCGwCar9nvByUg9g0di6AB7RilAnlFYPVUmNbLqs01/4VWj3ilZqqHQwLJqWnmrlixB06aEd/p0YDra5d0KLrChLTAmmkw9T6LQGd/hPRQ8FrjNV+z3yNKQOwbOhZBA9oxSgXyikDKaXizBiTEQYFyMHAFhId7pXaqh0MC6/Ycp8s780hKSaNInhimD2xJ8XxaKOwQb+iGmdwDNnwLhMFjm/QGNXRHgudarvma/S5RAmLf0LEIGtCOUSqQlwR+/A/MG2rVqOc0uLytl2qnujgkMOznWN78cZMZrX3N4ozu2ZCwMJ2S7hBv6IU5tgvergPpqVChNdz9degZqMWeFdB8zX7XKAGxb+hYBA1oxygVyEsChzbD8AZWjap1hO6TvFQ71cUhgdS0dG57dwFLth8xI77StTa3NS7rUHSFCTmB+SPgh6etZt82Cap3DDkCNdi7Apqv2e8bJSD2DR2LoAHtGKUCeU3gwy6w5VcIi4CH10I+rRHwWhc5UZ8/DpzgxhHzSEhKJV+OSGYOvpLLCuR0IrRihJrA+9fA7iUQnReGbIbImFATUHs9LKD5mv3OUQJi39CxCBrQjlEqkNcE1n0Fn95l1arN/0Hrx71WQ9XHIYHx87byzDfrzGjNKhbik/ub6VMsh2xDJsyeFfBexrlBdW6DmzO24Q0ZADXU6wKar9nvISUg9g0di6AB7RilAnlNIDUZ3qoJJ/dDvtIweBWER3itlqqPAwJpaencPW4Rc2LjzGjPdKrBPS0qOBBZIUJG4NxDTHv/DKUzzhMKGQA11OsCmq/Z7yElIPYNHYugAe0YpQJ5UeDn52DOG1bN7vgUqrT3Yi1VJwcEdh5O4Lq3ZxOflEqOqHB+fLg1ZQrlciCyQgS9QFoaDKsLR3dYO+cNWqmzP4K+0wOvgZqv2e8zJSD2DR2LoAHtGKUCeVHgyHYYahxSlw5VOsAdk71YS9XJIYHPl+3ikU9XmtHqlinA5w80JyJcu2I5xBu8YWJ/hEndrPY1HwjXPh+8bVXLAlZA8zX7XacExL6hYxE0oB2jVCCvCkzsBn/8CGHhMHg15C/t1ZqqXjYF0tPTuf/DJfy0/oAZ6enrq3P/lRVtRtXjQS8wrTes/sxq5sDlUEhjJuj7PAAbqPma/U5TAmLf0LEIGtCOUSqQVwU2TIfJd1i1a/0EtHnKqzVVvRwQiDt5mjavz+JEYgqR4WF8N6gVVYrndSCyQgSlgHFw6WuV4fQxKN0Eev8YlM1UowJfQPM1+32oBMS+oWMRNKAdo1QgrwqkpsDbteHEHshbEgavgYhIr9ZW9XJAYMbqvTwwaZkZqVapfEzt15wcUdqAwAHa4Auxcgp80cdqV7vnocXA4GujWhQUApqv2e9GJSD2DR2LoAHtGKUCeVng15fgt5etGuqAMS/3lGN1e2DiUmas2WfGe7htFQa1rexYbAUKIoEPO8OWWRARA4+sg9xFgqhxakowCWi+Zr83lYDYN3Qsgga0Y5QK5GWBY7vh7VqQngaV20OPT71cW9XNAYFjp5Jp++ZvHDxx2lyIPrXfFdQvW9CByAoRNAKnjsBrl0NaClTrCN0nBU3T1JDgE9B8zX6fKgGxb+hYBA1oxygVyOsCk26F2O+txejmyeiXeb3Gqp9NgR/X7TcXpRtXtRJ5+fLBFvoUy6ZpUD3++2iY+YTVpC6joF7GWrGgaqQaEywCmq/Z70klIPYNHYugAe0YpQJ5XWDd1/DpnVYtr/43XPmY12us+jkg8OinK5m2bJcZ6f5WFXj6hhoORFWIoBAY1QL2r4Ec+eGR9RCdOyiapUYEp4Dma/b7NagSkNTUVF577TXGjBnDzp07KVOmDL179+bxxx8nIuLSix7XrFnDkCFDmDt3rinbsmVLXn31VWrVqvU36azcm9lu0oDOrJTuC3iBlCR4szokxEHBCtZ2m2E6IyLg+/USDTiemMwNw+aw8/Ap81Osj3s3pWnFwsHebLXvUgKHt8Cw+tZdDe+BTkMv9YT+XAJ+FdB8zT5/UCUg/fv3Z9SoUdx77700b96c+fPnM27cOIyfjxw58h+1YmNjady4MYUKFWLAgAHmvcOGDePo0aMsWrSIypXPLprMyr1Z6SIN6Kxo6d6AF/j+aVgwwmrGPdOhfMuAb5IacGmBWRsPcM+4xeaNFYrkZsagVvoU69JswX3HLy/A7FetNvaYBpXbBnd71bqAF9B8zX4XBk0Csnr1aurWrWsmD0OHnv3tyaBBgxg+fDgrV66kdu3aFxXr1q0bM2fOZP369eabE+My3qJUr16dDh068NlnGQcjAVm5NytdpAGdFS3dG/ACBzfCyCZWM+p0h5vfDfgmqQGZE3jm67WMn7/NvPmOpmV58aaL/7s5cxF1V8AKpKVZb0NP7oN8pWHQSm3NHbCdGToV13zNfl8HTQLy9NNP8+KLL7JlyxYqVKhwRmbr1q1UrFiRp556ihdeeOGCYidPnqRw4cJ0796dCRMmnHfP3XffzZQpU4iLiyNPnjxk5d6sdo8GdFbFdH/AC4xpB7sWQWROeGyj9f23rqAXOJGYzI0j5rE1Lt788u7j3s24opI+xQr6jr9QA3cugrHtrD+58nG4+v9CkkGNDiwBzdfs91fQJCDt27c333Ls22ftNX/uVbx4cerXr2++4bjQtWDBAvOTLePzrX79+p13i/Ez4xMu455mzZqZf8/svVntHg3orIrp/v9n7zygrCqyNro70mQkRxWQnEEyCipBRcCAiooEQREkijIzOv+MTtARFclJUBADCqJiRBRBMojkHCVLUHLo+K97b4ughH6vqrtf+O5artHuqq/O2XUcz6FS0BNYNhE+TX1srOUgqN0l6F2SA2kjsGT7L9w3ZqHbuHCuOL7t35jsWfQoZdrohVCrz56EH8Z7Dj02G4qmngUJIRflSugRUL5mPqchU4A426tiY2NZtmzZn6jUrFmThIQEnG1aF/s+/PBDd1vV9OnTadWq1QVNnJ+1adOGqVOncs899+BL28tNz759+3D+Ov9ztn+1b9/e9cGxWZ8IhDyBs8fhlXKQcBKK1oTHvgt5l+Xg7wRe/HI9Y+Zsc39w3/XFGdi2mvCEE4H4U97bH86//wUrQvcFuowinOY/iH1VAWI+eSFTgJQuXRpnpcM5eP7Hz1mxOHDgAFu2bLkosUmTJtGhQwdmzJhB8+bNL2jz9ddf46yuOG2c4sCXtpebnueee47nn3/+ok1UgJgHthSCiMDHT8CKtz2DnQSkUKUgMl6mmhA4HZ/E7UPnuluxnG9SlzrcUKaAiaT6BhOB9Z/B+w95Ft82EOp2CybrZWsYE1ABYj75IVOAaAXEPBikIAKZQmDnInijhTd03e5w2/8yxQwNmjkE1uw5Suvh80hOgfw5svDNkzeSJ1ts5hijUTOWwJROsPYjb0znQdLcxTN2fI0mAn4SUAHiJ7jzuoVMAaIzIObBIAURyBQCKSkwvDYc3gxZr4L+GyE6S6aYokEzh8DArzYwcvZWd/C2tYrzyr3aipU5M5GBo5444N1+lZwI1zSEzl9k4OAaSgTMCKgAMePn9A6ZAsS55erFF19Ml1uwJk+ezOHDh9N0C9b5bX2dHgW0r8TUPmQIzBsM3/zTc+feCVDprpBxTY5cmUBCUjK3D5nL5gMn3MZjHq5Fi0qFr9xRLYKXwLIJ8Gkfz/57xkOVtsHriywPOwLK18ynPGQKEOcGLOemq0u9A7JixQqqVq3qHkbfunUruXPnpkiRIucIOgfMnTMgGzZsoHhxbxn4t3dAnNUV5/D5b58vbX2ZIgW0L7TUNqQIHP/Z+9PQlCQofQs8PC2k3JMzVyawft8xdytWQlIKBXNm4fPeN1Agp1bCrkwuSFu81Qa2zYaoWHh6K8TlClJHZHY4ElC+Zj7rIVOAOCicK3THjBnjvoTesGFD5s+f776E3q1bN0aPHu3S2rFjh/tOiPO+x4QJE84R3LhxI3Xq1HHfA+nd27sW1HkJ3Vn5cF5CL1eunF9tfZkiBbQvtNQ25Ai89yBs/NxZmIW+qyGP9yCovvAh8OrXGxk2y7ss5M7qRRncTleyhuTsH94Kw1Jveix/B7R7JyTdlFOhS0D5mvnchlQBkpiYyMCBAxk3bhy7d+92VzK6du3KgAEDiI727pe/VAHi/G7VqlVuW6dwcb5GjRrx0ksvuSsnf/x8aZvWaVJAp5WU2oUkgY1fwnvtPNeaPANN/hKSbsqpSxNITErmrpELWL3nqNtoSLvqtKleTMhCjcD3r8Csf3tetZ8G190Sah7KnxAnoHzNfIJDqgAxx5G5CgrozOWv0TOZQFIivFYRTvwMea6G3ishMjKTjdLwGU1g7d6j3D1yAWcTk90tWF/2ucG9HUtfCBEYcyPsW+ldOvHUFojSA5QhNLth4YryNfNpVgFiztCaggLaGkoJBSuBmf+E+YM96zt8AqWaBKsnstuAwLBvN/PqzE2uQvOKhRjb4XoDNXUNKAJ7lsHrN3sm1XgY2gwPKPNkjAikhYDytbRQunwbFSDmDK0pKKCtoZRQsBI4tAWG1/Ksr9wW2o4PVk9ktwEB51asdmMXseynX10VbcUygBloXb/8Kywe5Vn12GwoqnM+gTZFsufKBJSvXZnRlVqoALkSoQz8vQI6A2FrqMAl8MZtsHMBRGWBpzZ62zT0hR2BDfuP0XrYfOKTksmdNYbvnmpC3ux6oDCoAyE5GYZWgyM7IW8p6PUjREQEtUsyPjwJKF8zn3cVIOYMrSkooK2hlFAwE1jxLnzc3fPgtpeh7mPB7I1sNyAwbu42/vP5elfhpnIFeKNTbSKUsBoQzeSuW76Ft+/2jGjYB5r9K5MN0vAi4B8B5Wv+cTu/lwoQc4bWFBTQ1lBKKJgJxJ+EV8pB/HEoXBUenxvM3sh2AwLOVqz7xyzkx51HXBVtxTKAGQhdP3wUVn/gWeKsfuQrHQhWyQYR8JmA8jWfkf2pgwoQc4bWFBTQ1lBKKNgJTO8NP070vOj2PRSpFuweyX4/Cew4dJLmr33vbsXKHhvFrKeaUChXnJ9q6pZpBJIS4OXScOYoFK8NXb/JNFM0sAiYElC+ZkoQVICYM7SmoIC2hlJCwU5g9w8wLvVtgNqPQstXgt0j2W9AYML87Tz36TpXoYmzFatjbSIjdXbAAGnGd10zDaZ29sZt+hw06pfxNmhEEbBEQPmaOUgVIOYMrSkooK2hlFCwE0hJgZH14eB6iMsN/TdCTNZg90r2+0kgOTmFe0YvYHnqVqxX7q1G21rF/VRTt0whMOlu2PotRMZAv7WQs1CmmKFBRcAGAeVr5hRVgJgztKaggLaGUkKhQGDhCJjxjOfJ3eOg6r2h4JV88JPA7l9P0WzQ95xOSCJHlmj3gcISebP5qaZuGUrA2XY1sDQkJ0DZW+HB9zN0eA0mArYJKF8zJ6oCxJyhNQUFtDWUEgoFAicPwavlvaSlZGPoOD0UvJIPBgS0FcsAXmZ2/eEN+Cx1y1WroVCrY2Zao7FFwJiA8jVjhDoDYo7QnoIC2h5LKYUIgfcfhvWphUeflXDVtSHimNzwh0BKSgoPvL6IRdt+cbv/u00lHq6vmPCHZYb2GdMY9q2A2BzwpLOtMleGDq/BRMA2AeVr5kS1AmLO0JqCAtoaSgmFCoHNM+Gdtp43Nw6Am58NFc/kh58Edv1yituHzuX4mURyxkXzWa9GXJMvu59q6pbuBI7sgsGVvWGqPwR3jkz3ITWACKQ3AeVr5oRVgJgztKaggLaGUkKhQiA5CQZXgWN7IFcx6LsaIqNCxTv54SeBtxbu4B+frHV7Nyidj3e61tUDhX6yTPdu816Db57zhmn3LpRvme5DagARSG8CytfMCasAMWdoTUEBbQ2lhEKJwKz/wPcvex61/xCuaxpK3skXPwhdOoG5AAAgAElEQVQ4W7E6T1jK7I0H3d5/b1mBrjeU8kNJXdKVgHOb3ZBqcOQnyJYP+q2DGL3hkq7MJZ4hBJSvmWNWAWLO0JqCAtoaSgmFEoFftsPQ6p5HFe+E+1IfKAwlH+WLzwScrVgth87l2JlEssVGMaPvjboVy2eK6dxh30oYc6M3SINe0Pw/6Tyg5EUgYwgoXzPnrALEnKE1BQW0NZQSCjUCE+6AHXO9NwScN0Gy5ws1D+WPHwQ++GEXA6aucntWL5GHad0b6IFCPzimW5eZ/4T5gz35R76Gq+um21ASFoGMJKB8zZy2ChBzhtYUFNDWUEoo1Ais+gCmPep51eJFqN8j1DyUP34QcLZiPfrWMr5Z/7Pb+9nbK/DojdqK5QdK+10S4+HVcnD6F8hzNfReofNb9ilLMZMIKF8zB68CxJyhNQUFtDWUEgo1Agmn4ZVycPYoFKwI3RdARESoeSl//CCw/+gZmg2aw/GzicRGRTKj342UzK9bsfxAabfLlm/h7bs9zVv+ATf0t6svNRHIRALK18zhqwAxZ2hNQQFtDaWEQpHAZ0/CD+M9zx79DorVDEUv5ZMfBD5btZee7y53e1YqmouPn2hITFSkH0rqYo3Ap31g2QRPrucyyH+dNWkJiUBmE1C+Zj4DKkDMGVpTUEBbQymhUCSw50d4/SbPszqPwe2pN2OFoq/yyScCzlasbpOW8fU6byvWgFvL0aOJEl6fINpsfOYYvFoeEk5CocrQfb5NdWmJQKYTUL5mPgUqQMwZWlNQQFtDKaFQJOBc6TmyHhzcAFmv8g6jR2cJRU/lkx8EjpyK56ZXZvPrqQQiI2B6z0ZULpbbDyV1MSawagpM6+rJtBwEtbsYS0pABAKJgPI189lQAWLO0JqCAtoaSgmFKoH5Q2DmPzzv7psEFVuHqqfyyw8CM9f9zKNv/eD2rFAkl3srVtZYPVzpB0qzLh90gHWfQEQkPLUZsuc301NvEQgwAsrXzCdEBYg5Q2sKCmhrKCUUqgSO74dBFSAlGcreBg9ODlVP5ZefBHq/t5zpK/e6vZ+4qTRPtyjvp5K6+UXA+Xf0tUqQnAjX3gCdPvNLRp1EIJAJKF8znx0VIOYMrSkooK2hlFAoE3i7LWyZCRFR0H8D5CgYyt7KNx8JHD+TQIvXvmfv0TPuVqz3Hq1H3VJ6N8ZHjP43XzQKvvqr17/tG1D5Hv+11FMEApSA8jXziVEBYs7QmoIC2hpKCYUygTUfwtRHPA9bvAD1nwhlb+WbHwS+23iARyYsxTk2dF3BHEzv2ZBssdF+KKmLzwTebAk/zYOY7DBgK8Rk9VlCHUQg0AkoXzOfIRUg5gytKSigraGUUCgTSDgDr5T13gQpVAW6zwtlb+WbnwT+/vFq3l600+3dvt7V/OfOKn4qqVuaCRzaAsNrec0rtoH73kpzVzUUgWAioHzNfLZUgJgztKaggLaGUkKhTuDTvrDsTc/LbnOhSNVQ91j++UjgxNlEWg+fx7aDJ903Kyd2rsONZQv4qKLmPhGY9R/4PvV67Ic/gtI3+9RdjUUgWAgoXzOfKRUg5gytKSigraGUUKgT2LUUxjf1vKzXA259MdQ9ln9+EFi49TAPjlvkbsUqmjuOmU82JnsWbcXyA2XauoysDwfWQY7C8OR6iNRjkGkDp1bBRkD5mvmMqQAxZ2hNQQFtDaWEQp2Ak1EOrw2HN0O2/N5h9KiYUPda/vlB4KWvNjBq9la3553VizK4XQ0/VNTligT2LoexTbxmtTpBqyFX7KIGIhCsBJSvmc+cChBzhtYUFNDWUEooHAjMfRW+/Zfnabv3oPzt4eC1fPSRwKn4RO4asYCNPx93e77R6XpuLl/IRxU1vyKBz/vD0nFes66zoHjqWZArdlQDEQg+AsrXzOdMBYg5Q2sKCmhrKCUUDgSO7vHeGyAFyt8B7d4JB6/lox8E1uw56p4HSU6BPNli+K5/E67KHuuHkrpclICzIjmoIhzfC/nLQc8lAiUCIU1A+Zr59KoAMWdoTUEBbQ2lhMKFwFt3wrbvIDIG+m+E7HrvIVym3lc/h8/azCtfb3K73VqpMKPa1yTCOZ2uz5zAjnkwoaWn06gfNH3OXFMKIhDABJSvmU+OChBzhtYUFNDWUEooXAismgLTunre3jYQ6nYLF8/lp48EEpKSaTN8Puv2HXN7DnugBq2qFfVRRc0vSsB5l8d5n8f5eiyCghUESgRCmoDyNfPpVQFiztCaggLaGkoJhQuB+FPemyDxx6FINej2fbh4Lj/9ILD90En3lfT4pGRyZonmq343UiyPHsrzA+XvXZx3eV4uDfEnoEQ96DLDSE6dRSAYCChfM58lFSDmDK0pKKCtoZRQOBGY3gt+TH3wrPsCKOScC9EnAhcnMG7uNv7z+Xr3lzeVK8D4jrWJjNRWLL/jZe3HMKWj173FC1D/Cb+l1FEEgoWA8jXzmVIBYs7QmoIC2hpKCYUTgZ8Wwpu3eh7X7wkt/htO3stXHwmkpKRw35iFLN3xq9vzX20q0aH+tT6qqPk5AhNbw/Y5EBkNfVdDLm1rU3SEPgHla+ZzrALEnKE1BQW0NZQSCicCzg08Q2vAr9she0HvAbQoPTYXTiHgq6+7fz3FbUPmcvxMItljo/ikZ0OuK5jTVxm1P3kIXikDKclQoRXc/7aYiEBYEFC+Zj7NKkDMGVpTUEBbQymhcCMwZyB8l7ry8eAHULZFuBGQvz4SeGfxTzz70Rq3V+1rr2LyY/WJ0lYs3yg6734473843z3joUpb3/qrtQgEKQHla+YTpwLEnKE1BQW0NZQSCjcCR3bC4Cqe1xXvhPsmhhsB+esjAWcrVrdJy/h63c9uz35Ny9KnaRkfVcK4ubPyOKohHFgLsTm8lce4XGEMRK6HEwHla+azrQLEnKE1BQW0NZQSCkcCE+6AHXMhKtZ7EyRb3nCkIJ99ILD/6BnuGDaPQyfOEhsdybTuDahcLLcPCmHc9OBGGFHHA3D9I3DHa2EMQ66HGwHla+YzrgLEnKE1BQW0NZQSCkcCK96Djx/3PG/5KtROfR8kHFnI5zQT+HzVPp5490e3fZmCOfiizw3EREWmuX/YNpz9Esx+wXO/w3Qo1ThsUcjx8COgfM18zlWAmDO0pqCAtoZSQuFI4OwJ702QhJNQrBY8OiscKchnPwg8PWUlU5btdnt2aVSS/7ujoh8qYdQlKdHb8nh8b+rFD+sgKiaMAMjVcCegfM08AlSAmDO0pqCAtoZSQuFK4OMesOIdz/snlkCBcuFKQn77QODIqXj3Vqx9R88QEQFTH69PrWu0he+SCLd/DxNbeb++cQDc/KwPtNVUBIKfgPI18zlUAWLO0JqCAtoaSgmFK4Ed82BCS8/7hn2h2fPhSkJ++0hg0bbDtBu7yO1VOFcc3/RvTI4sus75ohg/fwqWvu79qsdiKFjeR9pqLgLBTUD5mvn8qQAxZ2hNQQFtDaWEwpVAcjIMrQbOrVg5i0C/tRAZFa405LePBP792TrGz9vu9rq3VnFevreajwph0NzZ6jioApw9BgXKQ49FuMtG+kQgjAgoXzOfbBUg5gytKSigraGUUDgT+O5FmPM/j0D7D+G6puFMQ777QOBMQpK7FWv7oZNur9c7XE+zioV8UAiDpqunwoddPEdvexnqPhYGTstFEbiQgPI184hQAWLO0JqCAtoaSgmFM4FftsPQ6h6Bym2h7fhwpiHffSSwYf8x7hg6j8TkFPJki+HLPjdQJHdWH1VCuPkHHWHdxxARCU9thuz5Q9hZuSYCFyegfM08MlSAmDO0pqCAtoZSQuFO4M3b4af5EB0HT22COL3tEO4h4Yv/w77dzKszN7ldbilf0F0JidQr6XDiAAyqCMkJcE0j6Py5L1jVVgRChoDyNfOpVAFiztCaggLaGkoJhTuBHyfB9J4ehVZDoFancCci/30gkJycQtvRC/hx5xG313OtKtKpYUkfFEK06cKRMONvnnP3jIcqbUPUUbklApcnoHzNPEJUgJgztKaggLaGUkLhTuDs8dQ3QU5BibrQ5etwJyL/fSSw65dTtBw6l2NnEskaE8W0Hg2oUCSXjyoh1vy3lcWY7DBgG8TEhZiDckcE0kZA+VraOF2ulQoQc4bWFBTQ1lBKSARg2mOw6n2PRK8fIV9pUREBnwh8sHQXAz5c5fapWjw3Ux9vQGx0mL6SfmgLDK/l8at4J9w30SeWaiwCoURA+Zr5bKoAMWdoTUEBbQ2lhEQAts2Gt9p4JG54Cm75P1ERAZ8J9H5vOdNX7nX7hfUr6d88B/Ne8/g9/DGUvslnluogAqFCQPma+UyqADFnaE1BAW0NpYREAJw3QQZXgWO7IVdx6LsaIsP0T68VD34TOHziLG1GzGf3r6dxzqFP6lKXhteF2c1PKSkw/Ho4vAVyFfPe19HbH37HlDoGPwHla+ZzqALEnKE1BQW0NZQSEgGPwLf/hrmveH/f4RMo1URkRMBnAvM2H6L9+MVuv0K5sjDzycbkiovxWSdoO/y0EN681TO/7uNw20tB64oMFwEbBJSvmVNUAWLO0JqCAtoaSgmJgEfg/H3rVdvB3WNERgT8IvDyjA2M+G6r27dphUK83qEWEeGyCjCtG6ya7HHrvhAKVfSLoTqJQKgQUL5mPpMqQMwZWlNQQFtDKSER+J3A+OawazFEZ019EyTMbzJSbPhFwHkl/Z5RC1i795jbf2Dbqtx3fQm/tIKqU1ICvFwazhyFYrXg0VlBZb6MFYH0IKB8zZyqChBzhtYUFNDWUEpIBH4n8MOb8Flf759bD4eaD4uOCPhFYPuhk7QY/D3xiclkiY7k6343ck2+7H5pBU2ntR/BlNR3dJo+B436BY3pMlQE0ouA8jVzsipAzBlaU1BAW0MpIRH4ncDpI/BqOUg8A9c0hM5fiI4I+E3g/Kt5KxXNxYfdGxAXE+W3XsB3nNgats+ByBjotwZyFg54k2WgCKQ3AeVr5oRVgJgztKaggLaGUkIicCGBqY/Amg+9n/VeAXn1qrVCxD8CKSkpPPrWD3yz/oAr0KNJaQbcWt4/sUDvdfIwvFIGUpKg/B3Q7p1At1j2iUCGEFC+Zo5ZBYg5Q2sKCmhrKCUkAhcS2PINvH2P97PGf4Wb/iZCIuA3gWNnEmj66hwOHD/rarz1SB1uLFvAb72A7bhoFHz1V8+8u8ZCtfsD1lQZJgIZSUD5mjltFSDmDK0pKKCtoZSQCFxIIDkJXqsEx/dBnquh90q9CaIYMSKwcOthHhq3iOQUKJwrjum9GlIwZ5yRZkB1dt7+GFEHDm2CuNzw5HqIDfHzLgE1ATImkAkoXzOfHRUg5gytKSigraGUkAj8mcDMf8L8wd7PO30O1zYSJREwIvDCF+sZ+/02V6NFpUKMbh9CV/P+vA5G1ff41HkMbn/ZiJU6i0AoEVC+Zj6bKkDMGVpTUEBbQykhEfgzgQMbYGRd7+fV28OdI0RJBIwIJCYlc++YhSzfecTV+ccdFXmkUYicL/r6/2DBUI9Px8+g5A1GrNRZBEKJgPI189lUAWLO0JqCAtoaSgmJwMUJjG0Ce5dDbA7vTRBtKVGkGBLYevAEd42Yz7EziWSNiWJq9/pUKprbUDWTuyfGe4fPzxyBq0pCr2UQGcI3fWUybg0ffASUr5nPmQoQc4bWFBTQ1lBKSAQuTmDxWPjyae93OlSrKLFEYMoPu3h66ipXrWyhHHzaqxFZooM4Yd/8DbyTemlD0+ehUeo7OpZ4SUYEgp2A8jXzGVQBYs7QmoIC2hpKCYnAxQmc+gVeKQvJCVDqJujwsUiJgDEB52pepwCZumy3q3X/9SV4qW1VY91ME/iwK6ye4g3feznkLZVppmhgEQhEAsrXzGdFBYg5Q2sKCmhrKCUkApcmMPkh2PAZEAH91kLuYqIlAsYEjp5KoM2Ieew4fMrVGt2+JrdWLmKsm+ECJw95D3cmJ0KJetBlRoaboAFFINAJKF8znyEVIOYMrSkooK2hlJAIXJrAhs9h8oPe72/5J9zwpGiJgBUCK3cd4a6R892reXNkieabJxtTOHeQXc27bCJ82tvjoW2KVuJCIqFHQPma+ZyqADFnaE1BAW0NpYRE4NIEnAO2g8rDqcOQvyw8sQQiIkRMBKwQGD9vO//+bJ2rVa14bqZ2b0BMVKQV7QwRmXAH7JgLUbHw9FaIy5Uhw2oQEQgmAsrXzGdLBYg5Q2sKCmhrKCUkApcn8OVfYPFor03XWVC8loiJgBUCyckpdHxzCXM3H3L1et50HU+1KGdFO91FDm70Hh90vkp3wb0T0n1IDSACwUhA+Zr5rKkAMWdoTUEBbQ2lhETg8gT2roCxjb02tbtCy1dFTASsEThyKp5bXp3D4ZPxruakLnW4oUwBa/rpJjT7fzD7RU++46dQ8sZ0G0rCIhDMBJSvmc+eChBzhtYUFNDWUEpIBC5PICUFRtaHg+shLg/03wgxQbZXX3Mc0AQWbTvMg68vcs+D5Mse617NWzRP1sC1OTnZW/04vBmyF/D+ndDbH4E7X7IsUwkoXzPHrwLEnKE1BQW0NZQSEoErE1gwDL7+u9eu7ZtQ+e4r91ELEfCBwMszNjDiu61ujxvK5OfNTrWJDtTzIFu+hbdT/x2o1wNuTV0J8cFfNRWBcCGgfM18plWAmDO0pqCAtoZSQiJwZQInDsCr5SElCUrfAg9Pu3IftRABHwg450Havb6IJdt/cXv1bVqGvk3L+qCQgU0/fgJWvO0N2OtHyFc6AwfXUCIQXASUr5nPlwoQc4bWFBTQ1lBKSATSRuC9B2DjF6lvgqyB3MXT1k+tRCCNBPYdPU3r4fM5ePwsUZER7irIjWUD7DxI/CkYVAHOHIEi1aDb92n0Ts1EIDwJKF8zn3cVIOYMrSkooK2hlJAIpI3A+W+C3Px3uPHptPVTKxHwgcA3637m0Uk/4Bw9KpI7jk96NqRgzgA6c3T+2x8tXoD6T/jgnZqKQPgRUL5mPucqQMwZWlNQQFtDKSERSBuBpATvT35PHoS8pbytJ3oTJG3s1MonAuefB6lfKh/vPlqXiECJtYmtYfsciM4KT22EuNw++abGIhBuBJSvmc+4ChBzhtYUFNDWUEpIBNJOYMazsHC4177zl3BNg7T3VUsRSCOB+MRkOr6xhIXbDrs9et18Hf2bB8D7IEf3wOAq3lmoCq3g/tRzIGn0S81EIBwJKF8zn3UVIOYMrSkooK2hlJAIpJ3AgfUwsp7XvvpDcOfItPdVSxHwgcDuX0/Ratg8fj2VQGQEvN21Lg1K5/dBIR2afvM8zBvkCd83CSq2TodBJCkCoUVA+Zr5fIZcAbJmzRoGDBjAvHnzXDqNGjVi4MCBVK5cOU20kpKSePnllxk3bhy7du2iRIkSdO3alaeffpqoqKhzGj/88ANvv/02s2bNYvv27WTPnp1KlSrxt7/9jaZNm6ZprD82UkD7hU2dRMCcwNibYO+PEJMdntoEWXKYa0pBBC5CYN7mQ7Qfv9j9TZ5sMczs15gCObNkDivnUMqQanDkJ8hdAvqs1NsfmTMTGjXICChfM5+wkCpANm/eTO3atcmbNy+9evVy6QwdOpQjR46wZMkSypQpc0ViPXr0YNSoUXTu3JkGDRqwYMEC3nzzTZyfjxgx4lz/tm3bMmfOHO655x5q1qzJiRMn3HZOATRy5Ei6d+9+xbFUgPiMSB1EIH0ILB0Pnz/pabcZATXap884UhUBYNDMTQz9drPLoubVeXi/W31iMuN9kJ2L4Y3m3pzU7wkt/qv5EQERSAMBFSBpgHSFJiFVgDhFwVdffcX69evdlQvnc1YxKlSowG233caUKVMui2P16tVUq1bNLV6GDBlyrm2fPn0YNmwYK1eupEqVKu7P58+fz/XXX0+WLL//ydXp06epXr06Bw8e5MCBA0RHR/s0Qwpon3CpsQjYI3D6CLxaDhLPwNX14ZGv7GlLSQT+QCAxKZkHX1/Mkh3e+yA9mpRmwK3lM57TlE6w9iNv3MfnQWHvv2/6REAELk9A+Zp5hIRMAeKsQOTLl4927doxceLEC8h07NiR999/n0OHDpEjx6W3Vjz77LO88MILbNu2jZIlS57TcLZYlSpVimeeeYb//vfyf0LUv39/Bg0axM6dO88VQWmdJgV0WkmpnQikA4GpXWDNVE+45zLIf106DCJJEfAI/HIynmaD5nD4ZLz7z290up6byxfKODxnT8DLpb2i+5qG0Nl5D0efCIhAWggoX0sLpcu3CZkCZOHChe6WKWf71OOPP36B187PnC1UTpt69VIPm16ES4sWLdxVjv379//pt4UKFaJGjRruCsvlvgceeICpU6e6276ccyG+fApoX2iprQhYJrD1O5h0pyfa6Elo+k/LA0hOBC4k4LyQ/uDri0hMTiFnXDTTezaiZH7f/rvhN9MV78HHqf+tbDkIanfxW0odRSDcCChfM5/xkClAPvzwQ5wtWNOnT6dVq1YXkHF+1qZNG7cwcM5sXOpztlfFxsaybNmyPzVxznkkJCTgbNO61Ods/XK2YN1xxx049lzu27dvH85f539O//bt27vjO+PpEwERyEACyckwpCoc3QU5i0C/tTqQm4H4w3Wo4bM288rXm1z3q5XIw3uP1iVbrG/bd/1iN6Yx7FsB0XHQbx1kz+eXjDqJQDgSUAFiPusBWYA4if7WrVvT5F22bNm4+uqrmTRpEh06dGDGjBk0b556qC5V4euvv8ZZ3XDaOAn+pb7SpUvjrHQ4B8//+DmrK865ji1btly0+9GjR6lfv75bVDirKI5Nl/uee+45nn/++Ys2UQGSpqlXIxGwT+C7F2DOS57uQ1OhTDP7Y0hRBM4jkJKSwqNv/cA36w+4P32w7tW8cFc6n8U4vBWGpf4hl3PhgnPxgj4REIE0E1ABkmZUl2wYkAXIjh07LjiDcTk3GzduzOzZs90Vh8xaAXEOnzsFztKlS90tWo5NV/q0AnIlQvq9CGQCgV93eNeSOl/FO+G+C8+TZYJFGjIMCBw5Fc/dIxew7dBJ19sX767CA3Uu/4dYRlhm/Re+H+hJPPwRlL7ZSE6dRSDcCKgAMZ/xgCxAnITeueI2LZ9z5W6dOnXc8x2ZcQYkPj6e1q1bu++BTJs2zd1+5e+ngPaXnPqJgEUCE+6AHXMhMgb6b9TWFItoJXVpAit2HaHd2IWcSUgmW2wU7z9WnyrFc9tHlhgPg8rDqcOQq5j39kdUjP1xpCgCIUxA+Zr55AZkAeKPW1e6BWvy5MkcPnz4srdgObdcvfjii2m+BSsxMdFddfn0009555133Bu4TD4FtAk99RUBSwRWvg8fPeaJ3foS1LvwUgtLo0hGBP5EYPKSnfx1mnfO8Jp82fisVyNyxlkuDjZ/A++knoW85R9wQ3/NhAiIgI8ElK/5COwizUOmAHF8cw6YO2dANmzYQPHixV13f3sHxNkidf7BcGcLlHNuwzn3ERPj/R+8c3bDuenqUu+ArFixgqpVq7ptk5OTefDBB93rfceOHcujjz5qPBsKaGOEEhABcwLxp7w3Qc4eg4KVoPt8iIgw15WCCFyBgHMe5NmP1/Du4p1uyxvK5OetR+oQYTP+pnSGtdM8S3qvgLy/XzmvCRIBEUgbAeVraeN0uVYhVYBs3LjR3Y7lvAfSu3dv12/nJXRn5cN5Cb1cuXLnWHTq1Ml9L8R54+Paa68993PnCt8xY8a4L6E3bNjQfXDQeeG8W7dujB49+ly7J598ktdee80979G1a9c/MW7WrJl7oN2XTwHtCy21FYF0JPBZP/jhDW+ArrOgeK10HEzSIvA7gVPxidw3ZiFr9hxzf9i3aRn6Ni1rB9GxffBaJUhJ0tsfdohKJUwJKF8zn/iQKkAcHKtWrWLAgAFu4eB8jRo14qWXXjq3cvEbsksVIM62qoEDBzJu3Dh2797trqQ4BYajef7L5k2aNLnsOZXvvvsOp40vnwLaF1pqKwLpSGDvChibeplEzQ7Qelg6DiZpEbiQwO5fT3Hb4LkcP5vo/uLtLnVpVCa/OaZFo+Grv3g6906ESqnv3pgrS0EEwoqA8jXz6Q65AsQcSeYpKKAzj71GFoE/ERh9A+xfBTHZ4amNkCWnIIlAhhFYsOUQD41fTEoKZI+N4qu+N1Iibzb/x3eExjgxvRpic8DTWyAmq/966ikCYUxA+Zr55KsAMWdoTUEBbQ2lhETAnMDScfB56gHdVkOhVkdzTSmIgA8Ezn+ksHzhnEzt3oAcWfx8pPCnhfDmrd7oWtXzYRbUVAT+TED5mnlUqAAxZ2hNQQFtDaWERMCcwJmj8Eo5SDwNxa6HR78115SCCPhAwDmU3mXiD8za4D1SeFeNYrx6bzUiI/24FOGLp2HJWG/0x+dD4co+WKKmIiAC5xNQvmYeDypAzBlaU1BAW0MpIRGwQ+Cj7rDyXU+r+wIoVMmOrlREII0EnEPpdwybx7aD3iOFf29Zga43lEpj79Rmzs1uzuHz079A/rLwxBLd7OYbQbUWgQsIKF8zDwgVIOYMrSkooK2hlJAI2CGwcxG80cLTqtMNbk99PdqOulREIE0ENv98nLajF3L0dALO4sf4TrW5qVzBNPV1Gy1/Gz55wmvf7F/QsE/a+6qlCIjAnwgoXzMPChUg5gytKSigraGUkAjYIeAc3B1RFw5thLjc3svoOrhrh61UfCLwxep9PPHuj+6h9Pw5YvmgW31KFciRNo237oRt30FULDy1CbJelbZ+aiUCInBRAsrXzANDBYg5Q2sKCmhrKCUkAvYILBgOXz/r6d39OlS9z562lETABwKDv9nE4G82uz0qFMnFh93rky32CofSf9kGw2pBSjKUux0eeM+HEdVUBETgYgSUr5nHhQoQc4bWFBTQ1lBKSATsETh5GAaVh6R4uKYRdP7cnraURMAHAknJKfR670e+WL3f7dWsYiHGPoGB/04AACAASURBVFzr8i+lz3gWFg73RnloKpRp5sOIaioCIqACJH1iQAVI+nD1S1UFiF/Y1EkE0p/AlM6wdpo3Ts8fIH+Z9B9TI4jARQgcPZXA/WMXsmH/cfe3vW6+jv7Ny12cVXKSd/j8+D64qiT0Xq7D54oqEbBAQPmaOUQVIOYMrSkooK2hlJAI2CWwbQ681drTbNALmv/Hrr7URMAHAjsOnaTVsHnnXkof9VBNbqtS5M8Kaz+CKZ28n984AG5O3Urow1hqKgIi8GcCytfMo0IFiDlDawoKaGsoJSQCdgkkJ8OwmvDrdsiWD55cD9FZ7I4hNRHwgcCibYd54PVF7qH02OhIPnmioXsu5IJvwh2wYy5ExkDfVZCrqA8jqKkIiMClCChfM48NFSDmDK0pKKCtoZSQCNgnMO81+OY5T7ftG1D5HvtjSFEEfCAwaeEO/u+TtW6P/Dmy8EXvRhTMFecpHNkFg6sAKVChNdw/yQdlNRUBEbgcAeVr5vGhAsScoTUFBbQ1lBISAfsEThyAQRUgORFKNoaO0+2PIUUR8IGA81L6Mx+t4b0lO91eVYvn5r1H65E9SzScf/i83btQvqUPymoqAiKgAiR9Y0AFSPry9UldBYhPuNRYBDKewPsPw/rUwqPXj5CvdMbboBFF4DwCzs1YzlasJdt/cX96R9UiDLm3MlGDynkvn+ct5V2cEBklbiIgApYIKF8zB6kCxJyhNQUFtDWUEhKB9CGw5Vt4+25P23lN2nlVWp8IZDKBY2cSuHPEfLYdPOlaMqTaLtps/Itn1S3/gBv6Z7KFGl4EQouA8jXz+VQBYs7QmoIC2hpKCYlA+hBwD6PXgF93QNa83mH0mNQ99+kzolRFIE0Eth48wf1jFnHoxFmmxD5H7chNEBHlHT7PXTxNGmokAiKQNgLK19LG6XKtVICYM7SmoIC2hlJCIpB+BOYPgZn/8PTvGgPV2qXfWFIWAR8IzN9yiH9N+IQZ0U+6vX6+5g4KdX7HBwU1FQERSAsB5WtpoXT5NipAzBlaU1BAW0MpIRFIPwLuy+gVIOksFK8DXWem31hSFgEfCWyc1JdyW990ez0e8X/07fYY5Qv/4XpeHzXVXARE4EICytfMI0IFiDlDawoKaGsoJSQC6Uvgo8dh5XveGN3mQpGq6Tue1EUgLQTiT8Ir5SD+OD8lF+Sm+EEUyZOd6T0bki+H3q1JC0K1EYG0EFC+lhZKl2+jAsScoTUFBbQ1lBISgfQlsGspjG/qjVGrE7Qakr7jSV0E0kJg1Qcw7VG35UeF+9BvR13378sXzsnHTzQkLkY3YaUFo9qIwJUIKF+7EqEr/14FyJUZZVgLBXSGodZAImBGwHl+esyNsH8VxGSD/hsgLreZpnqLgAkBJybHNYU9zpW70Zzps4HO729l4bbDruqNZQvwRsfriY6KNBlFfUVABADla+ZhoALEnKE1BQW0NZQSEoH0J7BsInza2xvntoFQt1v6j6kRROBSBHYvg3E3e7+t2g7uHoNzPW+b4fPZfsi7nrdtreIMvKcqkZER4igCImBAQPmaAbzUripAzBlaU1BAW0MpIRFIfwLOfvtXK8DZo5C/HDyxGCKU2KU/eI1wUQKf94el47xfPTYbitZw//bAsTPcPnQuh07Eu//ctVFJnm1ZgQjFqgJJBPwmoHzNb3TnOqoAMWdoTUEBbQ2lhEQgYwh8+RdYPNobq+OnUPLGjBlXo4jA+QScm9kGV4aEU1CgAvRYeEExvOnn49w7eiFHTye4vfo3K0uvW8qIoQiIgJ8ElK/5Ce68bipAzBlaU1BAW0MpIRHIGAIHN8GI2t5YFdvAfW9lzLgaRQTOJ7BwBMx4xvtJ6+FQ8+E/8Vm64xc6v7mUE2cT3d8937oSHRtcK44iIAJ+EFC+5ge0P3RRAWLO0JqCAtoaSgmJQMYRmNgKtn/vvTrdbw3kKppxY2skEXAOn49qCAfWQpbc8NRGiMl6US7fbTjA428v42xisvv7l9tW5d7rS4ihCIiAjwSUr/kI7CLNVYCYM7SmoIC2hlJCIpBxBNZNhw9S/8S58V/gptQ/ic44CzRSOBPYPBPeaesRuL4L3DHosjQ+XbmXPpOXk5wC0ZERDG5XnTuqqmgO5xCS774TUL7mO7M/9lABYs7QmoIC2hpKCYlAxhFISoQhVeHYHshRCPqugejYjBtfI4U3gffbw/pPgQjovRzylrwijw9+2MWAqavcdlmiIxnVviY3ly90xX5qIAIi4BFQvmYeCSpAzBlaU1BAW0MpIRHIWAJzXobv/uON2fYNqHxPxo6v0cKTwNHdMKQ6JCdAycbQcXqaOUyYv53nPl3nts8aE8UbnWpTv3S+NPdXQxEIZwLK18xnXwWIOUNrCgpoayglJAIZS+DEARhU0UsEr24Aj3yZseNrtPAkMONZWDjc8925AMG5CCGNX0pKCiNnb+XlGRvdHtljo3jn0XpUL5EnjQpqJgLhS0D5mvncqwAxZ2hNQQFtDaWERCDjCUztAmumeuM+Ph8KV854GzRi+BBISoBBFeDkQchbCnoug0jfXjl3ipCXvtrI6DlbzxUhUx5vQMWiucKHozwVAT8IKF/zA9ofuqgAMWdoTUEBbQ2lhEQg4wnsXAxvNPfGrfEwtEn9k+mMt0QjhgOB5W/DJ094nt70d2j8tF9eO0XIP6ev5a2FP7n9c2SJZmr3+pQvrCLEL6DqFBYElK+ZT7MKEHOG1hQU0NZQSkgEMp6Acx3q2CawbwVEx0G/dZBde+ozfiLCZMTRjWD/aojJBn1XQ/b8fjvuFCF/m7aayUt3uRo546KZ/Fg9KhXN7bemOopAKBNQvmY+uypAzBlaU1BAW0MpIRHIHAIrJ8NH3byxb/4/uPGpzLFDo4Y2gT0/wus3eT7W7Aithxr76xQhzs1YU5btPleETOpSV2dCjMlKIBQJKF8zn1UVIOYMrSkooK2hlJAIZA6BxLPwWmU4eQByFoW+qyAqJnNs0aihS+Dc1btA11lQvJYVX/+4EpIrLtq9Hev6a/Na0ZeICIQKAeVr5jOpAsScoTUFBbQ1lBISgcwj8N2LMOd/3vi6kjfz5iFURz7+s3f4PCUJSt4IHZ03QOx9ycnemZBJi7wzIU4RMrbD9dQrpe2E9ihLKdgJKF8zn0EVIOYMrSkooK2hlJAIZB4BJ0F8rZJ3JW/xOtB1ZubZopFDj8Cs/8D3L3t+3TsBKt1l3UenCPnXZ+uYsGCHq50zSzSjH65Fw+v8P2di3UgJikAmElC+Zg5fBYg5Q2sKCmhrKCUkAplLYFo3WDXZs+HRWVDMzhaZzHVKo2c6gYTT8Eo5OHsUchX3Xj6Pjk0Xs5wi5IUv1jNu3nZXP1tsFCMfqkmTcgXTZTyJikAwEVC+Zj5bKkDMGVpTUEBbQykhEchcAucfEq56P9w9NnPt0eihQeD8Sw5avAj1e6SrX04R8vLXGxk123snJDYqkmEP1qBFpcLpOq7ERSDQCShfM58hFSDmDK0pKKCtoZSQCGQ+gfHNYddiiIyBfmsgp5K2zJ+UILYgORlG1oNDGyEqCzy5PkOueXYOpg/+ZjNDvt3swouMgJfbVuOeWsWDGKZMFwEzAsrXzPg5vVWAmDO0pqCAtoZSQiKQ+QTWTIOpnT07bhwANz+b+TbJguAlsG02vNXGs7/OY3B76jmQDPLIeS39f19uODfa31tWoEujkkRERGSQBRpGBAKHgPI187lQAWLO0JqCAtoaSgmJQOYTSEqEIdXg2G7Ilg/6rYWYrJlvlywITgKTH4INn3m291wG+a/LcD+cm7Gem76WpOQUd+xujUvxlxbliXSWRfSJQBgRUL5mPtkqQMwZWlNQQFtDKSERCAwC8wbDN//0bGk9DGp2CAy7ZEVwETiw3tt+5XylboIOH2ea/Z+t2kv/D1ZyNjHZteGemsV54e7KZImOyjSbNLAIZDQB5WvmxFWAmDO0pqCAtoZSQiIQGARO/wqDKkLCKShQAXosBG1ZCYy5CSYrvvwrLB7lWdzpc7i2UaZaP2/zIbq/s4zjZxJdO5qUK8CwB2qQM06PbmbqxGjwDCOgfM0ctQoQc4bWFBTQ1lBKSAQCh8Dn/WHpOM+ehz+G0jcFjm2yJPAJnDgIQ6p6RWz+svDEkoAoYtfsOUqnN5dy6MRZl2GVYrkZ3/F6CuaKC3ymslAEDAkoXzMEqEPo5gBtKiigbdKUlggECIFDm2H49Z4xZZrDQ1MCxDCZERQEvn8FZv3bM7XNSKjxUMCYveuXU7Qfv5ifDp9ybSqaO45xHWtTsWiugLFRhohAehBQvmZOVSsg5gytKSigraGUkAgEFoF37oPNMzybui+EQhUDyz5ZE5gEEs/C0JreRQbZC8KT6yAqsLY5HTkVT8c3l7Jy1xGXofNg4aD7qnNrZV07HZhBJatsEFC+Zk5RBYg5Q2sKCmhrKCUkAoFFYPtcmHiHZ1PVdnD3mMCyT9YEJoEf3oTP+nq2BfBVzmcTk+j93nJmrP35HMenW5SjR5PSuqY3MCNLVhkSUL5mCFBbsMwB2lRQQNukKS0RCCACKSkwrins+QEio6H3cshzdQAZKFMCjoATM2NugP2rISY7PLkWsl4VcGb+ZtAfHyx0ft6qWlFeuqcK2WKjA9ZuGSYC/hBQvuYPtQv7aAXEnKE1BQW0NZQSEoHAI7D+M3g/df9+3cfhtpcCz0ZZFDgENs2Ad+/z7KnZEVoPDRzbLmPJpyv3MmDqKk4nJLmtKhbJxZiHa1Eib7agsF9GikBaCChfSwuly7dRAWLO0JqCAtoaSgmJQOARSE6GkXXh0CaIzuo9TJg9X+DZKYsCg8DEVrD9e4iIgp5LIV/pwLArDVY450G6v72MvUfPuK1zZ41h+IM1uKFMgTT0VhMRCHwCytfM50gFiDlDawoKaGsoJSQCgUlg+TvwSQ/PtsZ/gZueCUw7ZVXmEnC2XY1Ofeuj4p1w38TMtceP0Q8cO0PP95azZPsv53r3a1qWXjdfp5fT/eCpLoFFQPma+XyoADFnaE1BAW0NpYREIDAJJMbD0OpwbI+3n7/vGsiSIzBtlVWZR2DyQ7DhM2/8zl/CNQ0yzxaDkZ3D6S9+sYEJC3acU7mhTH6GtKtB3uyxBsrqKgKZS0D5mjl/FSDmDK0pKKCtoZSQCAQugYUjYEbqykeLF6F+6opI4FosyzKSwNHdMLgqpCRBycbQcXpGjp4uY32yYo97LuRsYrKrXzBnFvfl9LqltAUxXYBLNN0JKF8zR6wCxJyhNQUFtDWUEhKBwCVw9gS8VgnOHIFcxaD3CojWnwYH7oRlsGWf9oFlE7xB738HKqRe35zBZtgebvPPx3ls0jK2HzrpSkdEQK+by9DnljJERUbYHk56IpCuBJSvmeNVAWLO0JqCAtoaSgmJQGAT+O4FmJN6C9ado6D6g4Ftr6zLGAInD8OgCpB0FgpVhm5zITIyY8bOgFFOxyfxzEer+Wj5nnOj1bw6D8MerEmxPFkzwAINIQJ2CChfM+eoAsScoTUFBbQ1lBISgcAm4CSazipI4mnIXw56LAqpRDOw4QewdecXpneNgWrtAthY/037aPlunpm25txVvXExkfy7TWXa1iquhwv9x6qeGUhA+Zo5bBUg5gytKSigraGUkAgEPoEv/wKLR3t2tnsXyrcMfJtlYfoROHPMW/2IP+Ftzev1I8TEpd94mazsbMXqM3k5q3YfPWdJs4qFeOGuKhTImSWTrdPwInB5AsrXzCNEBYg5Q2sKCmhrKCUkAoFP4MhOGFLdO2xcvDZ0meltjNcXngQWjYKv/ur53vJVqN015DnEJyYz5NtNjJq9leQUz91ccdG8cHcV7qhaNOT9l4PBS0D5mvncqQAxZ2hNQQFtDaWERCA4CEzrBqsme7Z2+gKubRgcdstKuwTiT8HgKnDqEMTlhn7rwup65mU//ULv91aw58jpc1wbly3Ai3dXoajOhtiNNalZIaB8zRyjChBzhtYUFNDWUEpIBIKDwM/rYFR9z9brmkH7qcFht6y0S2DZRPi0t6fZ9Dlo1M+ufhCoXezNEOdsyF9vLU+H+tfq8cIgmMNwMlH5mvlsqwAxZ2hNQQFtDaWERCB4CLzbDjZ96dn7+DwoXCV4bJel5gQSz8Lw68HZkheTDZ5cD1nzmOsGqYKzGvLUlFXnrut13KhSLLe7GlK5WO4g9UpmhxoB5WvmM6oCxJyhNQUFtDWUEhKB4CGwczG80dyzt/I90PaN4LFdlpoTWPEufNzd02nYF5o9b64Z5Aq/nQ0Z8d3Wc544x6M61LuGp1qUI2dcTJB7KPODnYDyNfMZVAFiztCaggLaGkoJiUBwEXjjNti5wHmezbuSt2D54LJf1vpHICkRRjeCg+shOs5b/ciW1z+tEOy15cAJ+k9ZycpdR855lydbDANalKdd7RLalhWCcx4sLilfM58pFSDmDK0pKKCtoZSQCAQXgW2z4a02ns2V7oJ7U1/CDi4vZK2vBFa+Dx895vVybr1ybr/SdwGB5OQU3l2yk5e+3MDxs4nnfle2UA737ZC6pfKJmAhkOAHla+bIVYCYM7SmoIC2hlJCIhBcBFJSYEJL+Gm+Z/fj86Fw5eDyQdb6RiA5GV5vAvtWQlQW6LcGchT0TSOMWh89lcD/vtrAe0t2XuB1o+vy89fbyut8SBjFQiC4qnzNfBZUgJgztKaggLaGUkIiEHwEdsyHCbd7dpe/A9q9E3w+yOK0E1jzIUx9xGtf42FoMzztfcO45fp9x3hu+loWb//lAgrNKxaif/NylCucM4zpyPWMIqB8zZy0ChBzhtYUFNDWUEpIBIKTgLMNy9mO5XyPzYaiNYLTD1l9eQLOitfrN8He5d7qR+8fIXdxUfOBwNzNB/nn9LVsO3jygl63VynMsy0rUkzvh/hAU019JaB8zVdif26vAsScoTUFBbQ1lBISgeAksGsJjG/m2V6mBTz0QXD6IasvT2DddPjgYa9N9fZw5wgR84OAcz5k6o+7GTxzE3uPnrlA4c7qRfnrbRUonDvOD2V1EYHLE1C+Zh4hKkDMGVpTUEBbQykhEQheAm+3hS0zPfu7fAMlagevL7L8zwSc1Y8xN8L+VRAVC08shrylRMqAQGJSsns2ZNDMTfx6KuECpTbVi9KvaVmuzZ/dYAR1FYELCShfM48IFSDmDK0pKKCtoZSQCAQvgT0/ettznK/0zfDwR8Hriyz/M4ENX8DkB7yf1+wArYeJkiUCSckpTFywg9dmbrrgxixH/qZyBeje5DrqlNQ1x5Zwh7WM8jXz6VcBYs7QmoIC2hpKCYlAcBN47wHY+IXnQ+ev4Jr6we2PrPcIOKsfoxrAgXUQGQ09FkP+60THMoGziUmMn7ed4bO2cCo+6QL18oVz0qH+tdx7fXFioiItjyy5cCGgfM18plWAmDO0pqCAtoZSQiIQ3AT2r/YeqHO+a2+ATp8Ftz+y3iOw4XOY/KD399c/Ane8JjLpSMB5Uf2thTsY+/02Dhw/e8FIsdGR7mOGPZpcp3Mi6TgHoSqtfM18ZlWAmDO0pqCAtoZSQiIQ/AQ+6ADrPvH86PgplLwx+H0KZw+cdz9G1IbDW7zVjyeWQL7S4Uwkw3x3Dqt/umovb87fwYrzXlX/zYAaV+fh4XrX0KpaUa2KZNis2B/ImefIyAj7whdRVL5mjlkFiDlDawoKaGsoJSQCwU/gwHoY6Wy9SoES9eCRryAiY/7jGvzwAtCDtR/DlI6eYXW7w23/C0AjQ9+ktXuPMm7udj5avueizjrX+LapXoymFQoRlUHJbOhTTz8PE5KSWbj1MK/P3UaNEnl4snm59BvsPGXla+aYVYCYM7SmoIC2hlJCIhAaBKZ2gTVTPV/afwjXNQ0Nv8LNi6QEGFYTjuyEyBjovRzylAg3CgHl76n4RCYu+Mm9PWvnL6cuatst5Qtya+XCtK5elCzRUQFlfzga41wysO/oaT5ZsZcl239xV7OOnv791jPnz2fe6VKXBtflT3c8ytfMEasAMWdoTUEBbQ2lhEQgNAgc2gwj6kBKMhSuAo/NgUglQkE3uSvehY+7e2Y37APN/hV0LoSywat3H+XtRT+527T+eGj9N7+rFc/NzeUL0aJyIcoXzhXKOALGt4PHzzJvy0G+33SIlbuOsO3QhY9O/tHQMgVzMOKhmpQtlDPdfVC+Zo5YBYg5Q2sKCmhrKCUkAqFD4JOesHyS50+bkVDjodDxLRw8iT8FQ2vAif0QHQd910COAuHgedD56JwhmL/1ENN+3MPnq/YRn5R8UR+cnVkNr8tPg9L53dWRornjiND2SL/nOyUlxb0kYMHWQ8zZeJC1e4+x+cCJNOnFREVQ65qr3C1zD9W9hqyxGfMHNMrX0jQ9l22kAsScoTUFBbQ1lBISgdAhcHw/DK0JCSchZxHotQxi9aha0EzwotHw1V88c5v8DZr8NWhMD2dDnaR45e6jvL90JzPW/swvJ+Mvi8M5yO68MXJXjWKUKZhT50cuQsth6myZWrrjV37Y8Qvr9x9n0/7j7D924Sv2lwN9bb5s1Lj6KppVLETD0vnJnS0mU8JU+Zo5dhUg5gytKSigraGUkAiEFoHZL8HsF1KT2GegSWpCG1pehp43p37xzn6c/hXickO/tZAl/beHhB7IzPXISZy3HjzhHlyfs+kga/Ycu6JBRXLHUbFILioVzUW90vmoefVVZImODIuVkjMJSWw5cIJVu4/y485fXXZbD5zg2JnEK3I7v0GBnFlwtr45q02NyxagVIEcPvVPz8bK18zphlwBsmbNGgYMGMC8efNcOo0aNWLgwIFUrlw5TbSSkpJ4+eWXGTduHLt27aJEiRJ07dqVp59+mqioSy/tzZo1i1tuucUdY/PmzVx3ne+PSymg0zRFaiQC4Ucg/iQMqwXH90FMduj9I+QsHH4cgs3j716AOS95Vjf/LzToGWweyN6LEHBuXvphx69MX7mXxdsOX/FswvkSJfJmdVdInPMKta/NS7nCOSl+VdagKUycg+CHT55l35Ez7japNXuOsufIafb8etr93/MPhfsSPFfnzUaFIjmpViIPt1YqTIm82QL6SmTla77M7sXbhlQB4iT+tWvXJm/evPTq1cv1eOjQoRw5coQlS5ZQpkyZKxLr0aMHo0aNonPnzjRo0IAFCxbw5ptv4vx8xIgRF+0fHx9PtWrV3ILl5MmTKkCuSFkNREAEfCaw/B34pIfXrcbD0Ga4zxLqkIEEjuzyLhBIOAW5ikHvFRAdm4EGaKiMIpCYlOxu15qxdj9Ld3i3MzmP3vv65csei/On/s7qSaFccTivtpctnJMCObKQO2sMOeKi3du4frsd2OTcibOq49iYmJzC6fgkfj0V7xYWB46ddW8Fc4qJQyfO8vMx568z7hkN52FH069QrixcVzCHe1D8pnIFqVg0F/lzZDGVzfD+KkDMkYdUAdK2bVu++uor1q9f765cOJ9TFFSoUIHbbruNKVOmXJbY6tWr3ULCKV6GDBlyrm2fPn0YNmwYK1eupEqVKn/SePHFFxk8eDAPPvig+79aATEPTCmIgAj8gYDzkN3YG8F5JZ0IeHyudzOWvsAk8GlfWPamZ9udo6B66gvogWmtrLJM4PCJs+4WpFkbDrDx5+NsO3iCQycuf47EHxOcF91zZonG+V/nlaDk1KIiMTnZ/WenSHGKDafQOJuQfMmD9f6Mfak+ueKiuTZ/dkrmz07lormpeU0eKhbJnWEHxG36ciktFSDmlEOmADlx4gT58uWjXbt2TJw48QIyHTt25P333+fQoUPkyHHpPYTPPvssL7zwAtu2baNkyZLnNLZv306pUqV45pln+O9//3uB9k8//UTFihUZPnw4zt8///zzKkDM41IKIiACFyOwbQ681dr7Takm8PDHepwwECPFKRLHNIaUJChUxSsWdUtSIM5UhtvkrCY425acg9jO2YifDp90VxzOJJivLmSEM06hUzBnFgrniqNw7jicrVPli+SiXKGcbsHh3EplsjKTET7YGEMFiDnFkClAFi5c6G6ZcrZPPf744xeQcX7mbKFy2tSrV++S1Fq0aOGucuzfv/9PbQoVKkSNGjXcFZbzvzZt2nDgwAF3q5ZTfKgAMQ9KKYiACFyGwLv3w6bU/x96cAqUbS5cgUbgnXth89eeVQ9NhTLNAs1C2RNgBJxVirOJyez+9RS7fj3N+n3H2H7wJIdPxuOsphw/m+hulXLeKXEOeTttTT6nHs4eG01cTBQ5skS527tyZIkmT9ZYrsoe6xYYpQpkd89iFMuTlfw5YsOisEgrUxUgaSV16XYhU4B8+OGHOFuwpk+fTqtWrS7w2PmZUyhMnTqVe+6555I0nO1VsbGxLFu27E9tatasSUJCAs42rd++zz77zNV1zpfUqlWL5557Ls0FyL59+3D+Ov9zto61b9/eHd8ZT58IiIAI/InAwU0wsp73p+v5y0H3BRAVLVCBQmDT1/DuvZ411zWD9qkv2QeKfbIjpAg4hUtCUgrO4fCklBR321XUb4dEUj11/jkyIiJ1S5a3LUufGQEVIGb8nN4BWYA4if7WrVvT5F22bNm4+uqrmTRpEh06dGDGjBk0b37hnwh+/fXXOKsbThsnwb/UV7p0aZyVDmc144+fs7rirHRs2bLF/dXp06epVKkSzZo1Y8yYMe7PfClAfmt7MVtUgKRp6tVIBMKXwOf9Yek4z/+Wg6B2l/BlEUieO+d0RjWAg+shIgoe+w6KVAskC2WLCIiABQIqQMwhBmQBsmPHjgvOYFzOzcaNGzN79mwyegXk73//OyNHjmTTpk3kz5/f5wJEKyDmwSsFEQhbAicPea9rnz0G2fJBzx8gW96wxREwjp//6GDtR6HlKwFjmgwRARGwR0AFiDnLgCxAnNWFOXPmpMk758rdOnXquOc7MuoMyN69e91D6X379nXfi017jQAAHyBJREFUCPntc678dW7L+vbbb7n22mvdNr58CmhfaKmtCIQ5gflDYeb/eRCU7GZ+MJw5CoOrwpkj3qODTyzRWy2ZPyuyQATShYDyNXOsAVmA+OPWlW7Bmjx5MocPH77sLVjOLVfOlbpXugVrxYoV7oH0y33Zs2fHscmXTwHtCy21FYEwJ5AYD6MbwqFNEBEJj82BIlXDHEomuv9pH1g2wTPg1v9Bve6ZaIyGFgERSE8CytfM6YZMAeKgcA6YO2dANmzYQPHixV06v70D4pwBcbZp/fY5W6COHj2Kc+4jJibG/bFzA5ZTWFzqHRCn8Khatarbb+bMmX+i/8EHH7hvjTirIM74d955p08zpID2CZcai4AIbP0OJqX+/0yJuvDIDF33mhlR8fM67+wHKd7FAN2+h5i4zLBEY4qACGQAAeVr5pBDqgDZuHGjux3LeQ+kd+/eLh1nW5Sz8uHcVFWuXLlzxDp16uS+F+K88eFsl/rtc67wdQ6VOy+hN2zYkPnz57svoXfr1o3Ro0dflrgvh9AvJqSANg9oKYhA2BH4oAOs+8RzWw/eZfz0OwfPx90Ce3/0xu7wifdGiz4REIGQJaB8zXxqQ6oAcXCsWrWKAQMGuIWD8zVq1IiXXnrJXbk4/7tUAZKYmMjAgQMZN24cu3fvdlcynHMejmZ09OWvulQBYh6QUhABEfCRwJFdMKIOJJyCrHmh51LI7l2MoS8DCKycDB918waq0Brue0urUBmAXUOIQGYSUAFiTj/kChBzJJmnoIDOPPYaWQSCmsD5B9Irt4W244PanaAx/sRBGF4LnAPo0XHQeznkKho05stQERAB/wgoX/OP2/m9VICYM7SmoIC2hlJCIhBeBJISvW1A+1Z4fj/4AZRtEV4MMsPb8w+eN/8PNOiVGVZoTBEQgQwmoHzNHLgKEHOG1hQU0NZQSkgEwo/A/tUwtgkkJ0KuYtBjEcTlCj8OGeXxtjnwVmtvtPxloftCvUifUew1jghkMgHla+YToALEnKE1BQW0NZQSEoHwJPDtv2Duq57vtbtCy9S/D08a6ed1/Cl4/SY4uAGIgK7fQvFa6TeelEVABAKKgPI18+lQAWLO0JqCAtoaSgmJQHgSSDgDoxvB4c2e/52/gmvqhyeL9PT623/D3NRXzmt2gNbD0nM0aYuACAQYAeVr5hOiAsScoTUFBbQ1lBISgfAl8NMCePM2z/98ZeDxeXqTwmY07PkRxjeH5ATIVRyeWARZctocQVoiIAIBTkD5mvkEqQAxZ2hNQQFtDaWERCC8CXz2JPyQehNWoyeh6T/Dm4ct753D/uObwt7lnuJ9k6Bi6jkQW2NIRwREIOAJKF8znyIVIOYMrSkooK2hlJAIhDeBM8dgZD04tgciIr2tWFfXDW8mNryf9R/4/mVPqer9cPdYG6rSEAERCDICytfMJ0wFiDlDawoKaGsoJSQCIrB1Fky6y+NwVUlvK1aWHOLiL4F9q7yrjpPiIWcR6PY95Cjor5r6iYAIBDEB5Wvmk6cCxJyhNQUFtDWUEhIBEXAIfDEAlozxWNTqBK2GiIs/BJISYFRDOLTR6+28dl6xjT9K6iMCIhACBJSvmU+iChBzhtYUFNDWUEpIBETAIeBcFzu2MRza5PFo9x6Uv11sfCXw5V9h8SivV4320GaErwpqLwIiEEIElK+ZT6YKEHOG1hQU0NZQSkgEROA3Au6tTc28Bwrj8nhbsfKUEJ+0Etj+PUxs5bV2Hnh0+GXLm9beaicCIhCCBJSvmU+qChBzhtYUFNDWUEpIBETgfALzBsM3qTdhFa8Nnb6A6FgxuhKBk4dhRG04ddhr+cgMuLrelXrp9yIgAiFOQPma+QSrADFnaE1BAW0NpYREQATOJ5CcDO+1g80zvJ/W7wkt/itGlyOQnATv3gdbvvFa3fAU3PJ/YiYCIiACKF8zDwIVIOYMrSkooK2hlJAIiMAfCZz6BUbfAMd2e79p9y6UbylOlyIwfwjM/If326I1oMs3EBUtXiIgAiKgAsRCDKgAsQDRloQKEFskpSMCInBRAruWeK+ku+dBckO3uXDVNYL1RwLb56ae+0jxOPVYBLmKipMIiIAIuASUr5kHggoQc4bWFBTQ1lBKSARE4JJ/sj8UZqZuJSpaEzp/ATFZxes3Ar/ugPHN4cTPQETqSpFuDlOAiIAI/E5A+Zp5NKgAMWdoTUEBbQ2lhERABC5FICUF3nsANn3ptajcFu4ZBxERYhZ/Eia2hj0/eCwa9Ibm/xYXERABEbiAgPI184BQAWLO0JqCAtoaSgmJgAhcjsDpX2FcUzi8xWt107PQeEB4M3MKsw+7wJoPPQ5lmsMD70NkZHhzkfciIAJ/IqB8zTwoVICYM7SmoIC2hlJCIiACVyJweCu8fjOcOeK1bPsmVL77Sr1C9/ff/gvmvur5V6C8d1Vx9nyh6688EwER8JuA8jW/0Z3rqALEnKE1BQW0NZQSEgERSAsB55G9SXd5h9Kj47zzIMVqpaVnaLVZNgE+7eP5lPUq6PwVFCwfWj7KGxEQAWsElK+Zo1QBYs7QmoIC2hpKCYmACKSVwPnJd/aC0GUG5C2V1t7B327TDJj8oFeERURBx+lwbaPg90seiIAIpBsB5WvmaFWAmDO0pqCAtoZSQiIgAr4Q+OoZWDTC65HnGu/F71xFfFEIzrbOtcQTW0HiGc/+cN+GFpyzKKtFIMMJKF8zR64CxJyhNQUFtDWUEhIBEfCFgPNS+rSuvx/ALlgROn0O2fL6ohJcbfevhjduhfgTnt0tXoD6TwSXD7JWBEQgUwgoXzPHrgLEnKE1BQW0NZQSEgER8JVAYjxMfgC2fOP1dN4IeXiadyYi1L6f13nFx9mjnmeNnoSm/ww1L+WPCIhAOhFQvmYOVgWIOUNrCgpoayglJAIi4A8B5x0M51D6rsVe7yLV4OGPQ2slxFn5mNASzqQWH7W7wu2v6B0Uf+JFfUQgTAkoXzOfeBUg5gytKSigraGUkAiIgL8EnMT87bawe4mnUKgKdPgYsuf3VzFw+jlnPhzfflv5qNkRWg1R8RE4MyRLRCAoCChfM58mFSDmDK0pKKCtoZSQCIiACYGzx+Gde2HnQk/lqpLw0BTIX8ZENXP7OrddTekECac8O7TykbnzodFFIIgJKF8znzwVIOYMrSkooK2hlJAIiIApgbMn4L12sGOupxSXB+5/G0reYKqc8f2dq4Y/7+9dtet8DftA0+e18pHxM6ERRSAkCChfM59GFSDmDK0pKKCtoZSQCIiADQKJZ2F6L1j1vqcWGQOth0H1B2yop79GchLM/AcsHP77WM3+DQ16qfhIf/oaQQRCloDyNfOpVQFiztCaggLaGkoJiYAI2CKQkgJzXoLZL/6uWL8nNH0OomJsjWJf59QvMKUjOK+9u8VTNNw9FirfY38sKYqACIQVAeVr5tOtAsScoTUFBbQ1lBISARGwTWDl+zC9JyTFe8rFa3sP9+UpYXskc72di7zzHsf3eVpZ83pnWIpfb64tBREQgbAnoHzNPARUgJgztKaggLaGUkIiIALpQWDnYpjaGY7t8dSz5IaWr0KVtoGxpSkpAeYMhLmvQEqyZ2PRGtDuvfB42T095lyaIiACfyKgfM08KFSAmDO0pqCAtoZSQiIgAulF4ORh+KgbbJn5+wilb4YWL0LB8uk16pV19/wI03vDz6t/b1unG7T4b2BvFbuyZ2ohAiIQYASUr5lPiAoQc4bWFBTQ1lBKSAREID0JJCfDkrHeAe+ks95IEVFQuws0+VvGPlx4bK93PuXHt373OGcRuGMwlLs1PSlIWwREIEwJKF8zn3gVIOYMrSkooK2hlJAIiEBGEDi0Gb7624WrIc51vY36Qc0O6VuIHN0Ni0bBD2/8/raH43PV+6H5fyFHgYwgoDFEQATCkIDyNfNJVwFiztCaggLaGkoJiYAIZCSBzTNhxjNwaNPvo0bHeTdOOasixWrZsca5kct5l2TZRFj/6e+rL456/nJw+8tQqrGdsaQiAiIgApcgoHzNPDRUgJgztKaggLaGUkIiIAIZTcA5AL50nHdl7+lfLxy9cFUo1QRK1IUSdSBHwbRb5zyI+NMCb5XFKTr+v737jrmiSuM4/kiRRdBVROyAqBhABUtsIHYxsREriokYC9ZgiVgTS6xI7F2I/vGisWKNYjeKuhQrK4INBVSKrqwiVdj8Drnv3vcCL3M8Z96Zufc7yWYT38Mz53yeuffOM3PmTGllq1KE9Tua9R1q1mugWbNmyePSEgEEEPibApyv/U24sn9GARJuGC0CB3Q0SgIhgEBWAksWmP179IpiZObEVfdi/U4ripH2Xc3+8c8V/2vW3ExFzKL/rigy/jPNbPZkszlTzGz5ynG0DLAeMu/Rn4fMs8o1+0WgRgU4XwtPPAVIuGG0CBzQ0SgJhAACeRDQylQTRppNHWM2f054j3S3o3t/s54DzDp0z8fSv+GjIgICCBRMgPO18IRRgIQbRovAAR2NkkAIIJAnAT27oTsaM8abTf+X2fRxZrMm/f9dHavqq95c3q6LmaZv6W6Hnu3o0C1Po6IvCCBQowKcr4UnngIk3DBaBA7oaJQEQgCBvAssnr/iroieF1k4z2zZUrNmLc1atTVr08FMS+k2b5H3UdA/BBCoQQHO18KTTgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtAgd0NEoCIYAAAggggAACqQhwvhbOSgESbhgtwtixY61Pnz5WV1dn3bp1ixaXQAgggAACCCCAAAJxBCZPnmwnnXSSvffee9a7d+84QWssCgVIjhI+atQod0CzIYAAAggggAACCORbQBeMBw4cmO9O5rR3FCA5SszcuXNtzJgx1rlzZ2vdunXqPStV8NxxSZ06tR2Qw9RomywwOWwy6tR2RA5To22ywOSwyahT2VFT52/BggU2bdo069evn7Vv3z6VMVV7UAqQas9wI+NjDmPxk08OyWHxBYo/Aj6H5LD4AsUeAZ/B4uWPAqR4OYvWYz6w0SgzC0QOM6OPtmNyGI0ys0DkMDP6aDsmh9EoMwlE/jJhD9opBUgQX7H/MR/YYudPvSeH5LD4AsUfAZ9Dclh8gWKPgM9g8fJHAVK8nEXrMR/YaJSZBSKHmdFH2zE5jEaZWSBymBl9tB2Tw2iUmQQif5mwB+2UAiSIr9j/+KeffrIHHnjABg8ebJtuummxB1OjvSeHxU88OSSHxRco/gj4HBY7h+SvePmjAClezugxAggggAACCCCAAAKFFaAAKWzq6DgCCCCAAAIIIIAAAsUToAApXs7oMQIIIIAAAggggAAChRWgACls6ug4AggggAACCCCAAALFE6AAKV7O6DECCCCAAAIIIIAAAoUVoAApbOroOAIIIIAAAggggAACxROgAClezoJ7/Ndff9ktt9xiI0aMsOnTp9uWW25pp512ml188cXWvHnz4PgE8BeYMGGC1dXV2ZtvvmnfffedtWnTxnr06GGXXXaZHXjggQ0C+uTPp61/r/kXjQkolwcccIBr8tVXX9k222xT33zBggV29dVX26OPPmpz5sxxfxsyZIidfvrpK4X0aUtGwgW0nOe1115rL730ks2aNcs23HBD22233dyS5RtvvHH9Dh566CG744477Ouvv7aNNtrITjzxRJfT1q1br9QJn7bhI6jtCPpNU/7eeOMNUy6Vs3322ceuuOIK69q1az2Oz3ejT9va1vcb/R9//GHDhw+3iRMnmn4Df/75Zzv55JPtkUceWSmQTw7Saus3OlqvSYACZE1CVfj3s88+2+677z475ZRTbK+99rL333/fHn74YdN/v+eee6pwxPkf0jHHHGPvvPOOHX300bbzzjubvpiVk0mTJtm9995rZ511Vv0gfPLn0zb/SsXp4eLFi61nz56uwJ8/f/5KBcihhx5qr776qp177rnWvXt3d7L73HPP2bBhw9yFgPLNp21xhPLZUxWKffv2dUXEoEGDbIsttnAF4gcffOAu2my77bau48rTJZdcYkceeaQpP1988YXddddddsghh9iLL77YYHA+bfOpUpxe/fLLL7b99tubPn/6ztxqq61cgajfu7XWWss+//xzl1NtPt+NPm2Lo5V9T6dNm+ZypPeQ7bLLLu6zs7oCxCcHabXNXqy6ekABUl35XONo9AWsE6PzzjvPXb0rbbr6qh/QTz/91HbYYYc1xqFBXIGxY8farrvuaq1ataoPrCvfvXr1cidAs2fPthYtWrgf0KT582kbdzREu/HGG+322293V8X1/+V3QPQje/jhh9utt95qF1xwQT2WTmZfe+01+/77790VdW0+bVEPE1i+fLntvvvupqunuhjQtm3bVQbU57FTp0520EEHuaKxtN1222124YUXupypKNHm0zas9/xrCehizTnnnGPPP/+8+4yVtqefftp0kUc5Ov/88/kezcnhsmjRIps7d65tvvnmtnTpUmvZsuUqCxCf37K02uaErKq6QQFSVelc82B0G/qGG26wb7/91l15KG2a9tOlSxe7/PLL7frrr19zIFo0icBFF13kTlR/+OEHN1XOJ38+bZtkMDWyExUQuqtx9913u2LimmuuaVCADBw40J555hn79ddfG0zXeeutt2z//fe3Bx98sH4qlk/bGuFNbZilKXMvvPCCHXbYYbZw4UJr1qyZrb322g32qelUZ5xxhpsuud9++9X/7c8//3TTtY466igbNWqU++8+bVMbWA0Fvummm9y01fHjx7sLOqVNd7B0t1/T6JQ7n+9Gn7Y1RB19qI0VID45SKtt9AET0ChAauwg6Nevn7vLobmWlZvmyu600072yiuv1JhKfod7wgkn2FNPPWW//fabey7EJ38+bfMrULye6U6G7lhpaqOKj8oCZLvttrMNNtjAPvzwwwaD0x2vddZZxwYPHmz333+/+5tP2+JJ5avHQ4cOddOs3n77bXchRvnTtB3dFdFFgD333NN1WPlRkaiCo/J5D7WdN2+effnll95t86VRzN6o8NDzOnvssYd7tqBz585uCpbuTGlaq/6+3nrr8T2aw/Q2VoD4/Jal1TaHZIXvEgVI4VPoNwBNr9IVPT30Vbnp2YMlS5a429Ns2QtMnjzZTcHS1VhNIdDmkz+fttmPtjp6oOk3KkDGjRvn5jTroeTKAmTddde1gw8+uD6n5SNv166d9e7d23QVXptP2+oQzG4U/fv3d1Oq2rdvb3vvvbep+P/xxx/dA80qNpRTfaY0tUfFiZ43qNx090PT6H7//Xf3J5+22Y28uvas4lAFZHl+9FyP7jrqDhXfo/nMd2MFiM9vWVpt86lW7F5RgBQ7f96933rrrd2qIPoBrdx0i1pXbnXFiC1bAV1F1RVXreKiO1YdO3Z0HfLJn0/bbEdbHXvXHQytXKZnAzTVQ9uqChCtNHf88ce7FbAqt80228xN33r99dfdn3zaVodidqPQanNaOUkrl5X81Zt3333XPZh+7LHH2hNPPOH+rosDKk4qNz3zozY6mdLm0za7kVfXnkePHu3uIOpKuFaX0wU13dnSw+ljxoxxd5J9vht92laXZNOOprECxCcHabVtWo3a2BsFSG3kuX6UPlcHaowmN8PViax+PDVdQNPhtIRkafPJn0/b3Ay+wB258sor3UOwU6dOdVfRV1eA+NzV8GlbYLpcdF13K3QHS0uAaiWe8k1TefS51LK8Pnc1fNrmAqHgndBdDhX3H330UYPFVLTinL5Tb775ZtNUO5/vRp+2BefLtPvcAcmUP5OdU4Bkwp7dTn3mR2bXy9rds5aPPOKII9wDrvox1fSr8s0nfz5ta1c8zsh1NVyLOGiFHb1Tp7TdeeedbnU5XVnXSaza+DzX4dM2zkhqN8qZZ57p7ly9/PLLbjnd8k3PFOikVp9PngHJ7zGiizVaeUzLIlduevZDd7JUZPp8N/q0za9M/nvGMyD5z1HsHlKAxBbNeTzNjdUSoayClb9E6QtYS0Vq/r9W0RkwYMBKnfTJn0/b/GkUq0effPKJW8ChsU1TP/QgrKbpaJpIklWwfNoWSyx/vR05cqQrHrVyVXkRqZ6W3h0xY8YM9wC6ipDKVbB0h0TP8JSvguXTNn8ixeuRCnZtU6ZMadB5LbGsu4l9+vRxd5V9vht92hZPLD89bqwA8clBWm3zI1U9PaEAqZ5cJhqJnifQidLq3gOiE6kdd9wxUSwaxRNYtmyZOzF9/PHHGyzDWrkHn/z5tI03ktqMpGd29PBx5abnAZ588kl3F0QnsXrQWQWm7nKt7j0gejlXhw4dXCiftrUpH2/Ueh+B3u+hKTd6L4+ev9FWehfLqaeeaiNGjHDPyamdrow/++yz9R0ovQek/B0UPm3jjaR2I2kBCH1mlL/SqmXS0EqCeoZHS7Red9117rm6pL+DPm1rVz585I0VID45SKtt+AiJUClAAVKDx0RpqoHehK4Vd/Rlrbduly//WYMsmQ5Zy0TqBEZTCCqvvqpjerBZiwdo88mfT9tMAap056t6CF1D1RQfTcvShQA9dK6TXK3ApLuTl156aQMNn7ZVythkw9LLWTWNTqtgHXfccTZz5kzTNDrdvdLKgXoXj7bS+yZUUJa/CV0PsmsKV/nm07bJBlqlO9LiKvvuu697oavehq0HkvUQuu5E6e7Uxx9/bJtssgnfoznKv96XpGXmdRHuqquucoWh7iJq04Wa0gVRn9+ytNrmiK0qukIBUhVp9BuErjQMGzbMXc3TlAJdmdVJrx7O09u22ZpeQD+aevvy6ja9pE5ttPnkz6dt04+6+ve4ugJEy7rqx/axxx5zc9Z1ojRkyBB3EaBy82lb/aLpj7Curs7dndJzBHovi4p/FYZ6fqd80/MiKli++eYb9+Z6LdurJZf1byo3n7bpj7C69/DZZ5+5pZMnTJjgVipT4aEc6s6H7lyVNp/vRp+21a0bf3R6Nk4vbF3VpgujgwYNSvV3j9zGz2nSiBQgSaVohwACCCCAAAIIIIAAAsECFCDBhARAAAEEEEAAAQQQQACBpAIUIEmlaIcAAggggAACCCCAAALBAhQgwYQEQAABBBBAAAEEEEAAgaQCFCBJpWiHAAIIIIAAAggggAACwQIUIMGEBEAAAQQQQAABBBBAAIGkAhQgSaVohwACCCCAAAIIIIAAAsECFCDBhARAAAEEEEAAAQQQQACBpAIUIEmlaIcAAggggAACCCCAAALBAhQgwYQEQAABBBBAAAEEEEAAgaQCFCBJpWiHAAIIIIAAAggggAACwQIUIMGEBEAAAQQQQAABBBBAAIGkAhQgSaVohwACCCCAAAIIIIAAAsECFCDBhARAAAEEEEAAAQQQQACBpAIUIEmlaIcAAggggAACCCCAAALBAhQgwYQEQAABBBBAAAEEEEAAgaQCFCBJpWiHAAIIIIAAAggggAACwQIUIMGEBEAAAQQQQAABBBBAAIGkAhQgSaVohwACCCCAAAIIIIAAAsECFCDBhARAAAEEEEAAAQQQQACBpAIUIEmlaIcAAggggAACCCCAAALBAhQgwYQEQAABBBBAAAEEEEAAgaQCFCBJpWiHAAIIIIAAAggggAACwQIUIMGEBEAAAQQQQAABBBBAAIGkAhQgSaVohwACCCCAAAIIIIAAAsECFCDBhARAAAEEEEAAAQQQQACBpAIUIEmlaIcAAggggAACCCCAAALBAhQgwYQEQAABBBBAAAEEEEAAgaQCFCBJpWiHAAIIIIAAAggggAACwQIUIMGEBEAAAQQQQAABBBBAAIGkAhQgSaVohwACCCCAAAIIIIAAAsECFCDBhARAAAEEEEAAAQQQQACBpAIUIEmlaIcAAggggAACCCCAAALBAhQgwYQEQAABBBBAAAEEEEAAgaQCFCBJpWiHAAIIIIAAAggggAACwQIUIMGEBEAAAQQQQAABBBBAAIGkAhQgSaVohwACCCCAAAIIIIAAAsECFCDBhARAAAEEEEAAAQQQQACBpAIUIEmlaIcAAggggAACCCCAAALBAhQgwYQEQAABBBBAAAEEEEAAgaQC/wPk+6a49cOGnwAAAABJRU5ErkJggg==\" width=\"640\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1efc19ab048>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib notebook\n",
    "plt.plot(np.arange(v.shape[0]),v[1,:])\n",
    "plt.plot(np.arange(v.shape[0]),v[2,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remark :\n",
    "\n",
    "We observe that the graph of the principal components resembles a sinusoidal graph. This is an indication that the underlying dynamical system is \"close\" to become a linear dynamical system.\n",
    "\n",
    "Consequently, we can now use the DMD algorithm to determine the dmd modes of the system. Then we use them to make predictions for the future behaviour of the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The \"dmd\" and the \"time evolution\" functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dmd (A,r):\n",
    "    \"\"\"\n",
    "    A : data-matrix, snapshots of a \"linear\" dynamical system\n",
    "    Phi : the dmd modes\n",
    "    \"\"\"\n",
    "    \n",
    "    X = A[:,:-1]\n",
    "    Y = A[:,1:]\n",
    "    \n",
    "    u,s,v = np.linalg.svd(X)\n",
    "    U = u[:,:r]\n",
    "    S = np.diag(s)[:r,:r]\n",
    "    V = np.conjugate(v).T[:,:r]\n",
    "    \n",
    "    Atil = (((np.conjugate(U).T@Y)@V)@np.linalg.inv(S))\n",
    "    mu,W = np.linalg.eig(Atil)\n",
    "\n",
    "    Phi = (((Y@V)@np.linalg.inv(S))@ W) \n",
    "    \n",
    "    b = ((np.linalg.pinv(Phi))@X[:,0])\n",
    "    \n",
    "    return mu, Phi, b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_evolution (A,r,f):\n",
    "    \"\"\"\n",
    "    f : future time steps \n",
    "    \"\"\"\n",
    "    t = np.arange(A.shape[1]+f)\n",
    "    \n",
    "    Psi = np.zeros((r, t.size), dtype='complex')\n",
    "    j = dmd(A,r)\n",
    "    for i,_t in enumerate(t):\n",
    "        Psi[:,i] = (j[0]**_t)*j[2]\n",
    "    E = j[1]@Psi\n",
    "    \n",
    "    return E"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9458511828148355"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# How many principal components should we take into account ?\n",
    "\n",
    "np.sum(s[:90])/np.sum(s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = time_evolution(H,90,1) # We predict one step ahead in the future (f = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The prediction is :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "prediction for Andernach is 368.99214864096757\n",
      "prediction for Bingen is 237.18560332555586\n",
      "prediction for Bonn is 400.7814870555772\n",
      "prediction for Frankfurt Osthafen is 191.73546525223122\n",
      "prediction for Kalkofen Neu is 300.84628450330416\n",
      "prediction for Kaub is 266.087503738341\n",
      "prediction for Koblenz is 308.186699873881\n",
      "prediction for Koblenz Up is 215.04861597981798\n",
      "prediction for Mainz is 329.31826503519284\n",
      "prediction for Oberwinter is 321.7599448638986\n",
      "prediction for Oestrich is 219.56848271744104\n",
      "prediction for Raunheim is 163.65692096580145\n",
      "prediction for Rockenau Ska is 259.02734366078545\n",
      "prediction for Speyer is 390.26129566612656\n",
      "prediction for Worms is 237.82413634358193\n"
     ]
    }
   ],
   "source": [
    "for i,j in enumerate(df.columns):\n",
    "    print('prediction for',j, 'is' ,D.real[i-15,:][-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare with the last \"real\" (averaged) values : "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Andernach             368.555556\n",
       "Bingen                237.111111\n",
       "Bonn                  400.277778\n",
       "Frankfurt Osthafen    191.833333\n",
       "Kalkofen Neu          301.166667\n",
       "Kaub                  265.888889\n",
       "Koblenz               307.500000\n",
       "Koblenz Up            216.000000\n",
       "Mainz                 328.722222\n",
       "Oberwinter            321.222222\n",
       "Oestrich              219.611111\n",
       "Raunheim              163.111111\n",
       "Rockenau Ska          257.555556\n",
       "Speyer                389.500000\n",
       "Worms                 237.944444\n",
       "Name: 2014-12-31 23:00:00+00:00, dtype: float64"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = df.iloc[-1]\n",
    "d"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}