swissre_economic.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "universal-percentage",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import yaml\n",
    "import pycountry\n",
    "import pycountry_convert as pc\n",
    "from datetime import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "heard-spray",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/francescoposa/opt/anaconda3/envs/personal/lib/python3.8/site-packages/IPython/core/interactiveshell.py:3165: DtypeWarning: Columns (70,71,72,73,74,75,76,77,78) have mixed types.Specify dtype option on import or set low_memory=False.\n",
      "  has_raised = await self.run_ast_nodes(code_ast.body, cell_name,\n"
     ]
    }
   ],
   "source": [
    "data_swissre = pd.read_csv('./data/economic_weighted.csv', parse_dates=['date']).set_index('date')\n",
    "\n",
    "iso_codes = ['AUT','BEL', 'BGR', 'HRV', 'CYP', 'CZE', 'DNK', 'EST', 'FIN', 'FRA', 'DEU', 'GRC', 'HUN', 'IRL', 'ITA', 'LVA', 'LTU', 'LUX', 'MLT', 'NLD', 'POL', 'PRT', 'ROU', 'SVK', 'SVN', 'ESP', 'SWE', 'GBR', 'CHE', 'NOR']\n",
    "\n",
    "data_swissre = data_swissre[data_swissre['iso_code'].isin(iso_codes)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "swiss-reservoir",
   "metadata": {},
   "source": [
    "# Selecting columns for the first model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "given-championship",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('config/economic_features.yaml', 'r', encoding='utf-8') as file:\n",
    "    swissre_features = yaml.load(file, Loader=yaml.FullLoader)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cosmetic-appreciation",
   "metadata": {},
   "outputs": [],
   "source": [
    "columns_to_use= swissre_features['main'] +\\\n",
    "            swissre_features['demography'] +\\\n",
    "            swissre_features['sanitary'] +\\\n",
    "            swissre_features['weather'] +\\\n",
    "            swissre_features['economic'] +\\\n",
    "            swissre_features['policies']+\\\n",
    "            swissre_features['outcome']\n",
    "            #swissre_features['vaccines'] +\\\n",
    "\n",
    "rep_rate = 'epiforecasts_effective_reproduction_number_median'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "moved-louisiana",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['name', 'iso_code', 'cases', 'weight',\n",
       "       'epiforecasts_effective_reproduction_number_median',\n",
       "       'total_population_all', 'total_population_below_1',\n",
       "       'total_population_1_to_4', 'total_population_5_to_14',\n",
       "       'total_population_15_to_24', 'total_population_25_to_34',\n",
       "       'total_population_35_to_54', 'total_population_55_to_74',\n",
       "       'total_population_above_75', 'population_percentage_over_60',\n",
       "       'hospital_beds_per_1000', 'physicians_per_1000', 'nurses_per_1000',\n",
       "       'maxtempC', 'mintempC', 'tempC', 'FeelsLikeC', 'humidity', 'pressure',\n",
       "       'HeatIndexC', 'uvIndex', 'totalSnow_cm', 'sunHour', 'windspeedKmph',\n",
       "       'precipMM', 'unemployment_rate_2021_09_30',\n",
       "       'unemployment_rate_2021_12_31', 'unemployment_rate_2019_03_31',\n",
       "       'unemployment_rate_2019_06_30', 'unemployment_rate_2019_09_30',\n",
       "       'unemployment_rate_2019_12_31', 'gdp_2019_03_31', 'gdp_2019_06_30',\n",
       "       'gdp_2019_09_30', 'gdp_2019_12_31', 'gdp_2020_03_31', 'gdp_2020_06_30',\n",
       "       'gdp_2020_09_30', 'gdp_2020_12_31', 'gdp_2021_03_31', 'gdp_2021_06_30',\n",
       "       'gdp_2021_09_30', 'gdp_2021_12_31', 'school_closing',\n",
       "       'workplace_closing', 'cancel_public_events',\n",
       "       'restrictions_on_gatherings', 'close_public_transport',\n",
       "       'stay_at_home_requirements', 'public_information_campaigns',\n",
       "       'restrictions_on_internal_movement', 'international_travel_controls',\n",
       "       'unemployment_rate_idx'],\n",
       "      dtype='object')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = data_swissre[columns_to_use].copy()\n",
    "df.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "accomplished-avenue",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Taking the data starting from the 1st of april. ALTERNATIVE: use days since outbreak\n",
    "df = df[df.index >= '2020-04-01'].sort_index() #"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "existing-pride",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Filtering for valid countries: \n",
    "def filter_valid_countries(data, min_ratio):\n",
    "    iso_codes = []\n",
    "    for iso_code in data.iso_code.unique():\n",
    "        data_country = data[data.iso_code == iso_code]['epiforecasts_effective_reproduction_number_median'].dropna()\n",
    "        n_values = data_country.shape[0]\n",
    "\n",
    "    # If rows exist for that country with a valid r_estim value\n",
    "        if n_values > 0:\n",
    "            n_days = (data_country.index[-1] - data_country.index[0]).days + 1\n",
    "            #print(n_days)\n",
    "            ratio = n_values / n_days\n",
    "\n",
    "            if ratio > min_ratio:\n",
    "                iso_codes.append(iso_code)\n",
    "            \n",
    "    return data[data['iso_code'].isin(iso_codes)]\n",
    "\n",
    "df = filter_valid_countries(df, 0.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "administrative-boston",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of countries: 30\n"
     ]
    }
   ],
   "source": [
    "print(f\"Number of countries: {len(df.iso_code.unique())}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "theoretical-banner",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Putting NaN where the estimated reproduction number is larger than 4\n",
    "df[rep_rate] = np.where(df[rep_rate] >= 4, np.nan, df[rep_rate])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "challenging-evans",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Shifting the reproduction rate by 11 days for each country\n",
    "df['date_index'] = df.index\n",
    "df = df.sort_values(by = ['iso_code', 'date_index'])\n",
    "df['shifted_r_estim'] = df[rep_rate]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "growing-flour",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of countries: 30\n"
     ]
    }
   ],
   "source": [
    "# retaining only the countries that exceeded 2000 casese\n",
    "iso_codes = list(df[df['cases'] > 2000]['iso_code'].unique())\n",
    "df = df[df['iso_code'].isin(iso_codes)]\n",
    "print(f\"Number of countries: {len(df.iso_code.unique())}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "remarkable-museum",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The dataset does not include this data, while it is available\n",
    "df.loc[df['iso_code'] == 'CYP', 'population_percentage_over_60'] = 19.79"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sunset-questionnaire",
   "metadata": {},
   "source": [
    "# Training the first model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "palestinian-bloom",
   "metadata": {},
   "outputs": [],
   "source": [
    "from classes.context import Context\n",
    "from classes.pipeline import Pipeline\n",
    "from classes.plot import Plot\n",
    "import matplotlib.patches as mpatches"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "vital-newspaper",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Global seed set to 42\n",
      "Global seed set to 42\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training the model...\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Global seed set to 42\n",
      "GPU available: False, used: False\n",
      "TPU available: False, using: 0 TPU cores\n",
      "\n",
      "  | Name     | Type   | Params\n",
      "------------------------------------\n",
      "0 | lstm_1   | LSTM   | 2.6 K \n",
      "1 | linear_1 | Linear | 1.6 K \n",
      "2 | linear_2 | ReLU   | 0     \n",
      "3 | mixed_1  | Linear | 1.1 K \n",
      "4 | mixed_2  | ReLU   | 0     \n",
      "5 | mixed_3  | Linear | 16    \n",
      "6 | mixed_4  | ReLU   | 0     \n",
      "------------------------------------\n",
      "5.2 K     Trainable params\n",
      "0         Non-trainable params\n",
      "5.2 K     Total params\n",
      "0.021     Total estimated model params size (MB)\n",
      "The dataloader, val dataloader 0, does not have many workers which may be a bottleneck. Consider increasing the value of the `num_workers` argument` (try 12 which is the number of cpus on this machine) in the `DataLoader` init to improve performance.\n",
      "The dataloader, train dataloader, does not have many workers which may be a bottleneck. Consider increasing the value of the `num_workers` argument` (try 12 which is the number of cpus on this machine) in the `DataLoader` init to improve performance.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "..........................."
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7sAAADrCAYAAACo/IGcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7IUlEQVR4nO3deXxU1fnH8c/Dnhiwyg4RoiIgICICCiokisiPtcWFqm2NSgEVtCiKFJW1ilWx1oWKoiCKWgTbigIqomC1KpsC4kYJaowsIlICSoDz++PexJnJJJkkk2Um3/frdV/JnHvuuWeemSzPnHPPNeccIiIiIiIiIvGkWkV3QERERERERCTalOyKiIiIiIhI3FGyKyIiIiIiInFHya6IiIiIiIjEHSW7IiIiIiIiEneU7IqIiIiIiEjcUbIrIiIiIiIicUfJroiIxBwz+6WZrTSzHWZ2wMy2mdk/zKxvKdtNNTNnZqkBZX8wsyGl7XMJ+zPJzFwFnDfFj0N6hPUfNLOXwpS3MLOHzOxzM/vRzPaZ2QdmNsHMjg6o58xsWgFtP21mGSFlGf4x4bY/+HXMzNaZ2c3FeOoiIhJHalR0B0RERIrDzK4HHgCeAO4BsoETgf7AucDSUjS/FugOfBxQ9gfgbWBRKdqNW2Z2IjAC6BFS3hP4F7AD+CuwEagJnAlcBzQAxpTi1MuASWHKMwCcc87MpgBPmNls59zuUpxLRERikJJdERGJNWOBfzjnrg4oewN4zMxKNWPJObcX+E9p2ogGM6vtnPupovsRoT8AHzrnVucWmNkxwAvAZqC3cy47oP6rZnYfIclxCexyzhX1Wv0L+BEYBvy5lOcTEZEYo2nMIiISa44Fvg23wzl3BMDMGprZETP7Te4+MxvoT3N9OqAs0cwOmtm1/uOgacz+9NmWwOUB02TnBNQLt80Jaf9uM9vqn2erP4W3WkCd3LaGmNljZrYT2F7QkzezUWb2rpntNrM9ZvYfM+sfUid3GvIIM5tiZll+3ZfMLDmkbqKZPWJm3/nTjP8FBNUppC+1gd8A80N2/R5oCIwOSXQBcM5lO+dei+QcpeGcOwwswEt2RUSkitHIroiIxJr3gSvM7L/AP51zn4VWcM7tNLONeNOac5Pbc4EDQFpA1XPwptauKOBcvwJeAT7k5ymzO/2te0jdNOBOvNFMzKwG3lTbdsBUYAPeFN7b8RL2m0KOfxBYAvwWqFNAfwBSgMfxpuvWAAYCi82sn3NuSUjd8cA7wFVAI+A+4BmgV0CdR4GhwGTgA+B88ievBTkT+AWwKqS8N/Bt4GhvBMyPWb7y4tR3zh0KKVoJjDazE5xz/y1Gf0REJMYp2RURkVgzEm+K7J+BP5vZd8BrwJPOuVcD6q0ABgU8TgNmAjeaWRvn3Kd+2bfOuc3hTuScW2dmPxF+ymzeYzM7CbgZWMjP02UvBc4GejnnVvply80MYKKZ3e2c2xHQ3vvOuSJHIJ1zYwPOWw1YDrTGi0tosrvNOXdZQP2GwD1m1sw5942ZtQEuAyY456b71V41syS/vaKcCTjgo5Dy4/CvnS2GP/pbONvClF3mb0HMrGtIkr3O/3omoGRXRKQK0TRmERGJKf5I7ml4o5N/AtbjjcAuM7PbAqquAFLM7Hgzqw90BOYBn+GN8uJ/LWhUNyL+9amLgS+A3zrncldP7ouXpL1jZjVyN+BVfl6oKdCLEZ7vdDNbbGbbgUNADt5obJsw1V8OebzB/9rC/3oG3v8Cfw+p91wkfQGaAXudcwcjrF+YJ4CuYbbQBD7XkgLqfxxSb2dAX0VEpArRyK6IiMQc/1rMlf6GmTXDW4V5opk97Jz7HngTOII3evsD8D3edOQVQJqZPQN0Bh4raT/85PUFvGnHvZxzBwJ2N8K73jengMPrhzzOiuB8x+GN5H4MjAa+xEt4pwInhzkkdAXi3EWvcqdJN/W/hl4jXOA1wyHqBLQZ6CvglAjbyJUVbtqzmRW0ivLuCKdJ574mCcXsj4iIxDgluyIiEvP8KbmP492S6CS8KcF7zGw93ujtD8Cb/u1o3gAeAlKB6pRuZPdhvNHEs5xzoYtmfQdsBS4p4NiM0KcRwfn6AkcDlzjnvs4tNLPEiHqbX26C3ZjgKb6NIzz+O+CYMOWvA+eb2enOuTUl7Fu0HOt/3VWhvRARkXKnacwiIhJT/NHNcNr6XwOTzhV4I7tpeLcnyi1rAFwPfOWc+6KIU/5EmFFBMxsDXA382jm3Id9R3kjzccA+59zqMFtJkq/cpDZvtNjMWgNnlaAtgPfwRr9DE/JfR3j8J0DN0BWe8RbQ2gU8ZGZHhR7krwDdu7idLaHj/a+fltP5RESkktDIroiIxJqNZrYC7xrXrUA9oB/egkp/d859GVD3DbxVj5vhj+D6KzVvAs4DnorgfB8D55jZALxEepff3r3+8bvNLPD6253OuS14qx5fibco1X14U6hrASfiLZz1S+fc/mI+99fxpi0/5bfZFG8V5S8pwQfYzrlPzWw+MMVf7Cp3NeZ+ETaRu/BWNyBvpNk5t9vMLsS7z+1aM3sQ2Ih3rXI3fl5k7PXi9jlAg5C45/rWOZcR8PgMvA8HKvz+ySIiUr6U7IqISKwZh5eMTcGbbnsYb9GpW4G/hNRdhZccfuecC1y46A2gA5FNYR6Pd13v3/FGeOfiXQ9cDUj3t0BzgXTnXI6ZXeD3azjeCGM2sAVv4ahiL+rknNtkZpfjPfd/+W3dije9ObW47flGAPuAsXjJ+Bt4qxy/HUF/MszsfbzbHy0K2bfSzE7FW6V6DN69e3Pwbs30MPBICfub6wJ/C/UwMCrg8QDgXyX4YEFERGKc/bxopIiIiEjxmFk63rXSTStbQukvXPYlcIFzbnlF90dERMqXrtkVERGR0pgHZALXVnRHwrgZWKlEV0SkalKyKyIiIiXm3wbqKqBSjer6soDrKroTIiJSMTSNWUREREREROKORnZFREREREQk7sT1aswNGjRwKSkpFd0NERERERERKQNr1qzZ5ZxrGG5fXCe7KSkprF69uqK7ISIiIiIiImXAzLYVtE/TmEVERERERCTuKNkVERERERGRuKNkV0REREREROKOkl0RERERERGJO0p2RUREREREJO7E9WrMRdm7dy87duwgJyenorsiInGgZs2aNGrUiHr16lV0V2JCyq0vl7qNjOn9o9ATERERiUdVNtndu3cv27dvp3nz5iQkJGBmFd0lEYlhzjkOHDhAZmYmgBJeERERkQpWZZPdHTt20Lx5cxITEyu6KyISB8yMxMREmjdvzjfffJMv2dUopoiIiEj5qrLJbk5ODgkJCRXdDRGJMwkJCbo0QkpMH4qIiIhET5VNdgFNXRYpQx99vafUbXRM/kWp2yhv+r0iIiIiUjloNWYRERERERGJO1V6ZFcqvwEDBtCgQQPmzJlTaL309HROOOEE7rjjjjLtzwsvvMDFF1+Mcy6oPHAU8/+6d+TXV/yeK0aOLlbblW0UMzU1lQ4dOvDQQw+V+7m/3/0dqae2YsWKFaSmprJhwwb69u3LZ599xlFHHVXu/RGRihGNad2gqd0iIlVVuSa7ZjYeGAK0AX4C/gOMd85tLOK4U4CHgG7AbuBRYKoLzTiioMu019i172C0my1Qg6RarL7t/Ijrp6enM3fuXACqV69Os2bN6N+/P3feeSfHHHNMWXWzWFJSUhg1ahRjx44tl/Nt2LCBf/7zn2zbtq1czleUZxa/QUIxFj7L/OpL+vU4lQ8++IAuXbqUYc/K1ptvvklaWho7d+6kQYMGUW37lFNO4cwzz2TGjBncfvvtBdabM2cOo0aNYt++ffn2mRkLFizgoosuIiMjg+OPPz5sG0uWLKFv377MmTOHK6+8Mq+8UaNGdOvWjenTp9O+ffvSPykRERERKVPlPY05FXgE6AGcCxwCXjezYws6wMzqAa8B24GuwPXAzcCNZdHB8kx0S3q+3r17k5WVRUZGBo8//jgvvfQS1157bRn0rngOHizf2OV68MEHufDCCyvNrV6Ord+AhISKWeW7ol6D8nDllVcyc+ZMDh06FLU2ly5dSlZWVtB27rnn5u1PTEwkKyuLb775hpdffpns7Gz69+8f13EWERERiRflmuw65y5wzj3pnNvonNsA/BZoCJxVyGGXA4nAFf5xC4G7gRutiq4EU7t2bZo0aUJycjJ9+vRh6NChvPrqq0F1nnzySdq1a0edOnVo3bo1999/P0eOHMnbb2Y89NBD9O/fn8TERFq2bMnTTz8d1MaGDRvo3bs3CQkJHHvssaSnp/PDDz/k7U9PT2fAgAHcfffdJCcnk5ycTGpqKtu2bePmm2/GzIIW63nnnXfo1atX3u1ZrrnmGvbu3Zu3f//+/aSnp5OUlETjxo258847i4zF4cOH+fvf/86gQYOCyp9++mm6du1K3bp1adSoERdffHHe/U/BG4U0M5YvX84ZZ5xBYmIiXbp0Ye3atUHtPPXUU7Rs2ZLExEQGDBjA9u3bi+zT/3XvyNy/PZj3+NTjjuGFZ+YwdmQ6Z7RuTr+zOrF40fN5+/v1OBWArl27Ymakpqbm7YvkdXz44YcZMmQIRx11FH/84x/zntvixYvp1KkTderU4fTTT2fNmjVB/Vy0aBGnnHIKtWvX5rjjjuNPf/pTvunZkcY0IyODtLQ0ABo2bIiZcfsY7wMY5xxPznyA/medRrdWTbmwd4+g5w+wcf1aft0vla6tmnBJ355sWLc63/n79OnD7t27efPNNwsLf7HUr1+fJk2aBG21atXK229mNGnShKZNm9KlSxfGjBnDtm3b+PTTT6PWBxGR4ki59eVSbyIiVUVFL1BV1+/D94XU6Q6scs4dCChbBjQDUsqua7Hhv//9L0uXLqVmzZp5ZY899hh//OMfmTJlCps3b+a+++7j7rvv5pFHHgk6duLEiQwaNIj169czfPhwfve737F6tZdk7N+/n759+5KUlMT777/Piy++yDvvvMNVV10V1MZbb73FRx99xNKlS1m+fDmLFi0iOTmZO+64I2+kDLzEuU+fPgwaNIgPP/yQRYsWsX79+qD2xo4dy2uvvcbChQtZvnw569atY+XKlYU+/48++ogffvgh3/TfgwcPMnnyZD788EMWL17Mrl27uPTSS/MdP378eKZPn87atWupX78+l19+eV7C995775Gens7w4cNZv349AwcOLPE1wY/+5R5S+/wfC5at4oKBv2Li2NF88/WXADzz0nLg51HGRYsWAZG/jpMnT6Zfv35s2LCB6667Lq987Nix3H333axevZoTTjiB/v37s3//fgDWrFnDxRdfzJAhQ9iwYQPTp0/nrrvuKvT63MJietxxx7Fw4UIANm3aRFZWFrdMvguAh/48jRefe5o/TruHRW/8h6uvG8PUW29k5fJlAOzfn82o9KEkt0hh/uI3uOHWicyYlj/OtWrVolOnTrz11lvFfwGiYM+ePcyfPx8g6OdNRERERCqnil6g6gFgPfBuIXWaAF+HlG0P2Lc1cIeZDQeGA7Ro0SIqnaxsli5dSlJSEocPH+bHH38EYMaMGXn7p06dyp///GcuuugiAI4//nhuvfVWHnnkEUaNGpVXb8iQIYwYMQKACRMmsGLFCv7yl7/w9NNP88wzz7Bv3z7mzZtH3bp1AZg1axZpaWl88cUXtGrVCoA6derwxBNPULt27bx2q1evTt26dWnSpEle2T333MPQoUO56aab8spmzpzJaaedxo4dO0hMTGT27Nk88cQTXHDBBYA3qpmcnFxoLLZt25Y3+hYoMIk+4YQTmDlzJieffDJff/11UJtTp07NG5G84447OPvss8nMzCQ5OZkHHniA8847jwkTJgDQunVrPvjgA2bPnl1on8IZcOElDBgyFIDrxk5g/uxHWfv+uzRLbsEx9b3rW3NHGQP7FsnrOHToUIYNGxYUE4Dbb789Xyznz5/PsGHDmDFjBr169WLy5Ml5z+3zzz/n7rvvZvTo8AtrFRXTY4/1rkZo1KgRDRo0YMehPezfn828xx7hb88spPMZPQBIbtGSjR+u5fm5j9PzvAt45cUFHMo5yJT7HiLxqCROatuO73ffyIQbRubrQ7NmzcjIyChG5AvXs2dPqlUL/swvMzOTo48+GoDs7GySkpJwzuV9UDBo0CDatm0btT6IiIiISNmosGTXzGYAZwNnO+cOF1E9dG6lFVCOc24WMAugS5cuUV/AqjLo2bMns2bN4sCBAzz22GNs2bKF66+/HoCdO3fy1VdfMWLECK655pq8Yw4dOpRvimr37t3zPX75ZW960+bNm+nYsWNeogvQo0cPqlWrxscff5yX7Hbo0CEo0S3ImjVr+OKLL3j++Z+nr+b2Z8uWLSQmJnLw4MGgPiUlJXHKKacU2u6BAweoWbNmvoRl7dq1TJ48mfXr17N79+68c3355ZdByW7Hjh3zvm/WrBkAO3bsIDk5mc2bNzNw4MCgdrt3716iZLd1258XNKpRowbH1K/P7l27CqxfnNexoEWtwsXy448/BrzXt3//4NVJzz77bCZPnszevXvDXv8caUwD/fezT/nppx+55rcXB01pP3Qoh2bJ3odRWz//jJPatifxqKS8/ad27ha2vYSEBA4cOBB2X0nMnz+fDh06BJUFvucTExNZv349hw4dYuXKldx77708+uijUTu/iIiIiJSdCkl2zex+4NdAmnPuv0VU/xZvBDdQI/9r0RdQxqHExMS8ZPOvf/0raWlpTJ06lUmTJuVdz/m3v/2NHj16lPgczjkKuiQ6sDzS28AcOXKEYcOGMWbMmHz7mjdvXuJrIBs0aMDBgwfZv38/if4KyNnZ2VxwwQX07t2befPm0ahRI3bt2sU555yTb2GhwOmouc8rN4bRXOy7Rsi0VzPDBVx7G6o4r2NJbsUT6eubqzgxDXTEec/jr08+S9NmwQlxjZrerx+X/zOrAu3evZuUlJQC99erV48DBw6Qk5MT9Nru2bMHIG/ENldycnLez1I4Zpa3v23btmRlZXHppZeyYsWKiPssIiIiIhWj3K/ZNbMHgMuAc51zn0RwyLvAOWZWJ6DsfOAbICP6PYw9EydO5O677+abb76hcePGNG/enC1bttCqVat8W6D//Oc/+R6ffPLJALRr144PP/yQ//3vf3n733nnHY4cOZJXpyC1atXi8OHgwfrOnTuzadOmsH1KSEigVatW1KxZM6hP2dnZbNxY6F2p6NSpE0DeiCXAJ598wq5du7jzzjvp2bMnbdu2ZceOHYW2E067du3CxijacpOywJgV53UsSLhYBr6+b7/9dlD9t99+m+Tk5KCRzVyRxDR3YafA53HiSW2oVbs2WV9/RYvjTwjackd2TzipDZ9/+jH792fnHffRug/CPqeNGzfSuXPnAp9zmzZtOHLkCOvWrQsqz114rE2bNgUeG4kxY8awdu3avOuqRURERKTyKu/77D6MtwLzL4HvzSx3xHafc26fX+cuoJtz7jx/33xgIjDHzKYBrYFbgcllcZ/dWJSamkr79u2ZNm0ajzzyCJMmTWL06NH84he/oF+/fuTk5LB27VoyMzMZP3583nGLFi2ia9eupKam8sILL7B8+XLee+89AC6//HImTpzI7373O6ZMmcL333/PiBEjGDJkSJHJVkpKCqtWreI3v/kNtWvXpkGDBowbN44zzzyTkSNHMmLECOrWrcsnn3zCSy+9xKOPPkpSUhJXX30148aNo2HDhjRr1owpU6bkS5pDNWzYkM6dO/P222/nTedt0aIFtWvX5qGHHuK6665j8+bNhd6btSDXX389PXr04K677uKiiy7izTff5MUXXyx2O0U5tkFD6tRJYNmyZaSkpFCnTh2OPvroiF/HgkybNi0olrVq1eKyyy4D4KabbqJr165MmjSJyy67jA8++ID77ruvwBWwI4lpy5YtMTNefvllBg4cyP7snzgqqS5XDB/FjGm345zj9DN6sH9/Nh+t/YBq1apx0eXp9PvlRTz052lMvGk0I/5wMzu3f8vjD87I14eMjAwyMzPp06dPgc+5ffv29OnTJ++65BNPPJHPPvuMG264gUsuuSTfdfzfffcd3377bVDZ0UcfTUJCQtj269Wrx7Bhw5g4cSK//OUv802fFxGR8heNFaYzpvcvupKIxJzy/k/tWrwVmJcDWQHb2IA6TYETcx84537AG8ltBqwGHgbuA/L/NxwFDZJqFV2pEp7vxhtvZPbs2Wzbto1hw4bxxBNPMG/ePE499VTOOeccZs2axfHHHx90zKRJk1i4cCEdO3Zk5syZPPnkk3Tt2hXwpkovW7aMvXv30q1bNwYPHkz37t154okniuzLlClT+OqrrzjxxBNp2LAh4F0bu3LlSjIyMujVqxennnoq48ePp3HjxnnH3XvvvaSlpfGrX/2KtLQ0OnToQM+ePYs83/Dhw3nmmWfyHjds2JC5c+fyj3/8g3bt2jF58uSgBbwideaZZzJ79mxmzpxJx44dWbRoEZMmTSp2O0WpUaMG46ZM5/HHH6dZs2YMHjwYIOLXsSDTp0/npptuonPnznz++ecsXrw4b8pz586dWbBgAQsXLqRDhw7ceuut3HrrrUELXwWKJKbNmzdn8uTJTJgwgcaNG3PXbbcAcN3NExh54608NeshhvTuzsjLfsXyJS/R/LiWACQelcRf5zzHlxlb+HW/VGZMu50bxk/M14dnn32WPn360LJly0Kf9/PPP09aWhojR46kffv2jB49msGDB4e91rpv3740bdo0aAt8L4Vzww038Mknn/Dcc88VWk9EREREKpbF8+Boly5dXO6tdEJt3ry5yOm48czMWLBgQd5Kv7Hsp59+om3btjz11FOcc845FdKHj77eU+o2Oib/otRtgHcP4bS0NHbu3EmDBg2i0mZJRDMmP/30EyeddBLPPvssZ51V2G25K4dwv1808pCfYpKfYhIsWveEVUyCxVM8QDERqerMbI1zLuyKrZqDJzGvdu3azJ07l927d1d0V6SMbNu2jQkTJsREoisiIiIilUNF32dXJCoime4ssat169a0bt26orshIiIiIjFEyW4VFc/T16u61NRUvb4iIiJVlKb/i/xM05hFREREREQk7ijZFRERERERkbijZFdERERERETijpJdERERERERiTtKdkVERERERCTuaDVmERERERGJW9FYoVqrU8cmjexKgV544QXMLO/xnDlzSEpKqpC+DBgwgPT09Ao5t4iIiIiIxB6N7Ia65yTI3lF+5zuqEdz8ecTV09PTmTt3LgA1atTguOOOY8iQIUyePJmjjjqqrHoJwNChQ+nXr1/E9VNSUhg1ahRjx44tw1555syZw6hRo9i3b1++fWbGggULuOiii8jIyOD4448P28aSJUvo27cvc+bM4corr8wrb9SoEd26dWP69Om0b9++zJ6DiIiIiIhEj5LdUOWZ6JbwfL1792bevHnk5OSwatUqhg0bRnZ2NjNnzsxX99ChQ1SvXj1ohLakEhISSEhIKHU7lcHSpUs59dRTg8qOPfbYvO8TExPZsmULzjkyMzO55ZZb6N+/P5999hm1atUq7+6KiIiIiERNVZnarWnMMah27do0adKE4447jssuu4zLL7+cf/zjHwBMmjSJDh06MGfOHE488URq165NdnY2P/zwA8OHD6dRo0bUrVuXXr16sXr16qB2n3rqKVq2bEliYiIDBgxg+/btQfvDTWN++eWXOeOMM0hISKB+/foMHDiQH3/8kdTUVLZt28bNN9+MmQUl2++88w69evUiMTGR5s2bc80117B37968/fv37yc9PZ2kpCQaN27MnXfeGeUIQv369WnSpEnQFpjEmhlNmjShadOmdOnShTFjxrBt2zY+/fTTqPdFRERERESiT8luHEhISCAnJyfv8datW5k/fz4LFizgww8/pHbt2vTv35/MzEwWL17MunXr6NmzJ+eeey5ZWVkAvPfee6SnpzN8+HDWr1/PwIEDueOOOwo979KlSxk8eDDnn38+a9asYcWKFfTq1YsjR46waNEikpOTueOOO8jKyso7z4YNG+jTpw+DBg3iww8/ZNGiRaxfv56rrroqr92xY8fy2muvsXDhQpYvX866detYuXJlGUQuMnv27GH+/PkA1KxZs8L6ISIiIiIikdM05hj3/vvvM3/+fM4777y8soMHDzJv3jwaN24MwBtvvMH69evZuXNn3jTkqVOn8tJLLzFv3jxuueUWHnjgAc477zwmTJgAQOvWrfnggw+YPXt2geeeOnUqF110EdOmTcsr69ixI+BNA65evTp169alSZMmefvvuecehg4dyk033ZRXNnPmTE477TR27NhBYmIis2fP5oknnuCCCy4A4MknnyQ5Obm0oQrSs2dPqlUL/qwnMzOTo48+GoDs7GySkpJwzrF//34ABg0aRNu2baPaDxERERERKRtKdmPQ0qVLSUpK4tChQ+Tk5DB48GAefPDBvP3Jycl5iS7AmjVr2L9/Pw0bNgxq58cff2TLli0AbN68mYEDBwbt7969e6HJ7rp164q9QvKaNWv44osveP755/PKnHMAbNmyhcTERA4ePEj37t3z9iclJXHKKacU6zxFmT9/Ph06dAgqq1u3bt73iYmJrF+/nkOHDrFy5UruvfdeHn300aj2QUREREREyo6S3RjUs2dPZs2aRc2aNWnWrFm+qbWhqzIfOXKExo0bs2rVqnxt1atXD/g54SxrR44cYdiwYYwZMybfvubNm5f4mth69epx4MABcnJyguKxZ88egLwR21zJycm0atWqwPbMLG9/27ZtycrK4tJLL2XFihUl6p+IiIiIiJQvXbMbgxITE2nVqhUtW7aM6BrSzp07s337dqpVq0arVq2CtkaNGgHQrl07/vOf/wQdF/o41Gmnncby5csL3F+rVi0OHz6cry+bNm3K149WrVqRkJBAq1atqFmzZtC5s7Oz2bhxY6F9adOmDUeOHGHdunVB5WvXrs3bXxpjxoxh7dq1LFq0qFTtiIiIiIhI+VCyWwX07t2bs846i8GDB7NkyRK2bt3Ku+++y8SJE/NGe6+//npef/117rrrLj7//HMee+wxXnzxxULbnTBhAgsWLOC2227j448/ZtOmTdx///1517impKSwatUqMjMz2bVrFwDjxo3j/fffZ+TIkaxbt44vvviCxYsXM2LECMCbsnz11Vczbtw4XnvtNTZt2sRVV12VL2kO1b59e/r06cOwYcN4/fXX2bp1K8uWLePaa6/lkksuoUWLFkH1v/vuO7799tug7cCBAwW2X69ePYYNG8bEiRM5cuRI4QEXEREREZEKp2Q31FGN4u58ZsYrr7zCueeey+9//3vatGnDJZdcwqeffkqzZs0AOPPMM5k9ezYzZ86kY8eOLFq0iEmTJhXabr9+/XjxxRdZsmQJp512Gr169WLFihV5Cz9NmTKFr776ihNPPDHveuGOHTuycuVKMjIy6NWrF6eeeirjx48Pusb43nvvJS0tjV/96lekpaXRoUMHevbsWeTzfP7550lLS2PkyJG0b9+e0aNHM3jw4LDXHfft25emTZsGbc8880yh7d9www188sknPPfcc0X2RUREREREKpaV17WaFaFLly4u9F6yuTZv3szJJ59czj2SePXR13tK3UbH5F+Uuo3KpCrHJNzvl6py8/biUEzyU0yCRSMeoJiEiqd4gGISSj83+ek9kl88xcTM1jjnuoTbp5FdERERERERiTtajVmKrSqP2ImIiIiISGzQyK6IiIiIiIjEHSW7IiIiIiIiEnfKPdk1s55m9i8zyzQzZ2bpRdRP8euFbn1L25d4XpxLRCqGfq+IiIiIVA4VMbKbBGwEbgAKvrFpfn2BpgHbG6XpRM2aNQu9r6qISEkcOHCAmjVrVnQ3RERERKq8cl+gyjn3CvAKgJnNKcah3znnvo1WPxo1akRmZibNmzcnISEBM4tW0yJSBTnnOHDgAJmZmUH3jRYRERGRihFLqzEvMrM6wOfA/c65F0rTWL169QD45ptvyMnJiUL3qo7t35d+RHzz/xKi0JPKQzHJryrGpGbNmjRu3Djv94uIiIiIVJxYSHb3AWOBfwOHgEHA82Z2hXPu6dDKZjYcGA7QokWLQhuuV6+e/iktgf+Lo5tQR4tikp9iIiIiIiIVqdInu865XcB9AUWrzawBcAuQL9l1zs0CZgF06dJFK8WIiIiIiIhUQbF666H3gJMquhMiIiIiIiJSOcVqstsJyKroToiIiIiIiEjlVO7TmM0sCWjlP6wGtDCzTsBu59yXZnYX0M05d55f/wogB1gHHAEGAtcB48q77yIiIiIiIhIbKuKa3S7AioDHk/1tLpCOdw/dE0OOuQ1oCRwGPgOuCrc4lYiIiIiIiAhUzH123wQKvKmtcy495PFcvERYREREpPgmHR2FNn4ofRsiIlKuYvWaXREREREREZECKdkVERERERGRuFPp77MrIiJSIE1PFRERkQIo2S0v+odMRERERGKF/neVOKBkV0REREREpDBK/mOSkl0RERGRqkT/tOenmIjEJS1QJSIiIiIiInEnqsmumZ0azfZERERERERESiIq05jN7HTgDmAAUD0abYqISIhoTLMDTbUTERGRKiGiZNfMhgDDgBbAF8A059xqM2sN3Af0A/YD08uqoyIiIiIiIlJJxMC17kUmu2Z2JTAb2AN8BnQFVpjZVcCTfht/AaY753aWWU9FREREREREIhTJyO4o4A1gsHMu28yqAX8FngW2Av2dc5+VYR9FREREREREiiWSBaraAg8457IBnHNHgDv9Y29ToisiIiIiIiKVTSQjuwnA9pCy3MdbotsdqTK00I6IiIiIiJShSG895AooPxytjoiIiIiIiIhES6S3HnrJzA6GKX/FzHICHjvnXMso9EtEJCZW+RMRERGRyimSZHdumfdCREREREREJIqKTHadc1eWR0dEREREREREoiXSacwiIiISCzT9X0REBIhggSozm2Fmx4WUDTGzo0PKWpvZv6LdQREREREREZHiimRk9wbgOeArADOrDiwAugJrA+odA/SPdgdFqgyNxoiIiIiIRE0ktx6yCMtEREREREREKoVI77MrIiIiIiIiEjO0QFUEUm59udRtZNSJQkdEREREREQkIpGO7LoIy0REREREREQqXKQju7PM7H8hZbPNbF/A47qRNGRmPYGxwOlAM+BK59ycIo45BXgI6AbsBh4FpjrnlHCLiIiIiIhIPpEkuyvJP4r7Vph6e/y6RUkCNgJP+VuhzKwe8JrfdlegDTAHyAbui+B8IiIiIiIiUsUUmew651KjeULn3CvAKwBmNieCQy4HEoErnHMHgI1mdjJwo5nN0OiuiIiIiIiIhCoy2TWzcwvZfQjY7pz7NHpdyqc7sMpPdHMtA6YCKcDWMjy3iIiIiIiIxKBIpjG/jjeNuaB76zozywLGO+fmRa1nP2sCfB1Stj1gn5JdERERERERCRJJsptWyL7qeItMXQTMMbPvnXOLo9KzYKFTla2AcsxsODAcoEWLFmXQFREREREREansIrlmN9xiVKGeNrO/AzcD0U52v8UbwQ3UyP+6PaQc59wsYBZAly5ddD2viIiIiIhIFRTpfXYj8SzQKYrt5XoXOMfM6gSUnQ98A2SUwflEREREREQkxkUz2f0fUKuoSmaWZGadzKyTf/4W/uMW/v67zGx5wCHzgf1406Q7mNkQ4FZAKzGLiIiIiIhIWNFMds8AtkVQrwuwzt8SgMn+91P8/U2BE3MrO+d+wBvJbQasBh7Gu7/ujGh1XEREREREROJLJAtUFcrMagCDgHHAX4qq75x7k4JXdsY5lx6mbAPQs6R9FBERERERkaolkvvsfkWYVY991YEGfjuvA9Oi1zURERERERGRkolkZHc5BSe7h4AdwArn3OtR65WIiIiIiIhIKURy66H0cuiHiIiIiIiISNREc4EqERERERERkUpBya6IiIiIiIjEHSW7IiIiIiIiEneU7IqIiIiIiEjcUbIrIiIiIiIicUfJroiIiIiIiMQdJbsiIiIiIiISd5TsioiIiIiISNxRsisiIiIiIiJxR8muiIiIiIiIxB0luyIiIiIiIhJ3lOyKiIiIiIhI3FGyKyIiIiIiInFHya6IiIiIiIjEHSW7IiIiIiIiEneU7IqIiIiIiEjcUbIrIiIiIiIicUfJroiIiIiIiMQdJbsiIiIiIiISd5TsioiIiIiISNxRsisiIiIiIiJxR8muiIiIiIiIxB0luyIiIiIiIhJ3KiTZNbNrzWyrmf1oZmvM7JxC6qaYmQuz9S3PPouIiIiIiEjsKPdk18yGAg8AdwKnAe8AS8ysRRGH9gWaBmxvlGU/RUREREREJHZVxMjujcAc59xjzrnNzrnRQBZwTRHHfeec+zZgO1j2XRUREREREZFYVK7JrpnVAk4HXg3Z9SrQo4jDF5nZDjP7t5ldVCYdFBERERERkbhQ3iO7DYDqwPaQ8u1AkwKO2QeMBS4B+gHLgefN7DfhKpvZcDNbbWard+7cGZ1ei4iIiIiISEypUUHndSGPLUyZV9G5XcB9AUWrzawBcAvwdJj6s4BZAF26dAnbpoiIiIiIiMS38h7Z3QUcJv8obiPyj/YW5j3gpGh1SkREREREROJLuSa7/qJSa4DzQ3adj7cqc6Q64S1qJSIiIiIiIpJPRUxjngHMM7P3gX8DI4FmwN8AzOwuoJtz7jz/8RVADrAOOAIMBK4DxpV/10VERERERCQWlHuy65x73szqA7fh3S93I9DPObfNr9IUODHksNuAlnhToD8DrnLO5bteV0RERERERAQqaIEq59wjwCMF7EsPeTwXmFsO3RIREREREZE4Ud4LVImIiIiIiIiUOSW7IiIiIiIiEneU7IqIiIiIiEjcUbIrIiIiIiIicUfJroiIiIiIiMQdJbsiIiIiIiISd5TsioiIiIiISNxRsisiIiIiIiJxR8muiIiIiIiIxB0luyIiIiIiIhJ3lOyKiIiIiIhI3FGyKyIiIiIiInFHya6IiIiIiIjEHSW7IiIiIiIiEneU7IqIiIiIiEjcUbIrIiIiIiIicUfJroiIiIiIiMQdJbsiIiIiIiISd5TsioiIiIiISNxRsisiIiIiIiJxR8muiIiIiIiIxB0luyIiIiIiIhJ3lOyKiIiIiIhI3FGyKyIiIiIiInFHya6IiIiIiIjEnQpJds3sWjPbamY/mtkaMzuniPqnmNlbZnbAzDLN7A4zs/Lqr4iIiIiIiMSWck92zWwo8ABwJ3Aa8A6wxMxaFFC/HvAasB3oClwP3AzcWC4dFhERERERkZhTESO7NwJznHOPOec2O+dGA1nANQXUvxxIBK5wzm10zi0E7gZu1OiuiIiIiIiIhFOuya6Z1QJOB14N2fUq0KOAw7oDq5xzBwLKlgHNgJRo91FERERERERinznnyu9kZs2ATKCXc25lQPkdwOXOuTZhjnkV+No5d1VAWQtgG9DDOfduSP3hwHD/YRvg06g/kZJpAOyq6E5UIopHfopJfopJMMUjP8UkP8UkP8UkmOKRn2KSn2ISTPHIr7LEpKVzrmG4HTXKuye+0AzbwpQVVT9cOc65WcCsknetbJjZaudcl4ruR2WheOSnmOSnmARTPPJTTPJTTPJTTIIpHvkpJvkpJsEUj/xiISblfc3uLuAw0CSkvBHeAlThfFtAfQo5RkRERERERKqwck12nXMHgTXA+SG7zsdblTmcd4FzzKxOSP1vgIxo91FERERERERiX0WsxjwDSDezYWZ2spk9gLfY1N8AzOwuM1seUH8+sB+YY2YdzGwIcCsww5XnBcelV+mmVlcwxSM/xSQ/xSSY4pGfYpKfYpKfYhJM8chPMclPMQmmeORX6WNSrgtU5Z3U7FrgFqApsBEYk7tglZnNAVKdcykB9U8BHga6Ad/jJcZTYizZFRERERERkXJSIcmuiIiIiIiISFmqiGnMIiIiIiIiImWqyia7ZjbezD4ws71mttPMXjKzDiF1zMwmmdk3ZnbAzN40s/YB+481swfN7BN//1dmNtPM6oe0c4yZzTOzH/xtnpn9IoI+nmJmb/ltZ5rZHWZmAftTzcyF2dpWxXj4dWqZ2RQz22pmP5nZl2Z2fXHjES8xMbM5BbxHsqtqTPw6l5nZejPbb2bfmtnTZha66ntVisd1ZrbZr/Opmf2uuLGIlZiYWR3/5+IjM8sxszcLqNfLzNaY2Y9m9l8zG1mVY2JmTc1svn/+w+ZdclSV4zHEzF71+/8/M3vPzAbFSEwmmNm/zSzbzCKe3mdF/72J2nvEby8eYhK190mcxKOXmb1jZt/5dT4xs7EliUe8xCSk7tlmdsjMNhY3FgFtxHxMLIr5DQDOuSq5AcuAK4EOwCnAi3i3OTo2oM444H/AhX69v+OtAl3X398BWAQMAloBvYBNwKsh51ril/cAuvvfv1RE/+r5/fm7f54L/b7cFFAnFe9ew+3wbs+Uu1WvivHw6y0E3sdbsTsFOAPvGvCq+h45OuS90QTYAjxZhWNyFt4t0MYAxwNnAmuB5VU0HtcA+4BLgROAX/t1Bsbpe+QovHUfhgP/AN4MU+d4IBt4EDgZ+D2QA1xYhWOSAvwVSMe7e8KcksQijuLxAN5imd3880/E+71yTgzEZApwE/AnwEXYv0h+l0TtPRJHMYna+yRO4nE63t+Y9ni/Z3+D97v22qr6HgmoewzwX/85baziPzepRCm/cc5V3WQ3TPCT8H4BDfQfG5AFTAiok+C/ICMKaacfcASo5z8+2X/Bzgqoc7Zf1qaQdq4B9gIJAWW3AZn8fK117puhgeLhAPoAP5RFPGI1JmGOOctvt0dVjQkwFtgWctyVwL4qGo93gPtDjrsPeDse3yMhbT5E+ETmbuDzkLLHgXerakxC6iymlIlMPMUjoO77wH2VOSYh+y4i8n9Qi/X3JtrvkXiISbTfJ3EUj0XAs1X9PeLHYSIwiVIku/EQE6Kc31TZacxh1MWb1v29//h4vE8RXs2t4Jw7AKzE+4S4IPWAn/BulwTeJ8n7CL6P8L/xPskqrJ3uwCr/nLmW4d2mKSWk7mozyzKz5WaWVkibxRGL8fgl8AFwo5l9bWafm9lfzSypkHaLIxZjEur3wCbnXEH3tS6uWIzJv4GmZjbQn8rTAO+T5lcKaTdSsRiP2sCPIccdALqZWc1C2o5UZYtJJLoH9s+3DOhShWNSluIlHnX5+TlEo62yiElJleTvTbTFS0yi9T6J+XiY2Wl+394q5blzxWRMzLtLTRNgWinPF05MxsQXlfxGye7PHgDWA+/6j3Ov39seUm97wL4g5l0HNBV4zDl3KKCdnc7/qAL8jz5gR0HtBBwX7tyBfcvC+4TkQmAI8Cmw3Mx6FtJupGIxHifgfWp/Kl5MRgF9gTmFtFscsRiTwHMfDVwMPFZIm8UVczFxzr2LN2X3GeAgsBPvk84rCmk3UjEXD7w/MleZWVc/+e8CDANqAg0KaTtSlS0mkSgobjWoujEpSzEfDzO7DkgG5kWpybKKSUkV6+9NGYn5mET5fRKz8fAHJH4CVgOPOOf+Vspz54q5mJh3e9WJwOXOucOlPF84MRcTopzfKNkFzGwGXpJ0YZg3mgutHqYMMzsKeAlvGP6WItoIasfMNpnZPn9bUsS588qdc5865/7mnFvjnHvXOXctsBRvmmaJxWo88N7PDrjMOfeec24ZXsJ7oZk1DnPOiMVwTAL9BqhOlP4Zi9WYmFk7vOvKpuJdP9QX7xfso2HOF7FYjQdeHF7GG/3KAf4JzPX3leoPbyWOSSSK87MVsRiPSdTFQzzM7ELgHrx/VreVpI2Q9so6JkWdPxp/b6IqHmISzfdJHMTjHKALMBL4g5n9tjjnL6BPMRcTM6sNPAeMdc5tLc75IuxTzMUEop/f1CjJQfHEzO7Hm8KY5pz7b8Cub/2vTYCvAsobEfKJhD9NNncK5ADnXOCUwG+BRmZmuZ8um5kBDQPa6Yc3igLe9MHc40I/YWnkfw39RCTQe/7zKZEYj0cWkOmc+yGgzmb/a4vQfkYqxmMS6PfAQufc7nDPszhiPCbjgfedc/f4jz8yb3XqVWY2wTn3FcUUy/HwpxJdZWYjgMZ4P0fD8a7f2VX4My9YJY5JJAqK2yHgu2K0EyTGYxJ18RAPP4GZB/zOOfev4h4fpr2yjkkkovk/SanFQ0yi+T6Jh3gEJHYbzBuMmEQpPoiP4Zg0xVuE6Ukze9Ivr+Z1xw4B/ZxzoZfURCSGY1KQkuc3LkoXQMfihje0vx04Ocy+3Au4/xhQVgfvouoRAWV1gbfxrgOqG6adkwlZEAhvTrwjsoVl6gSU/ZGiL/R/EXijKsYD7x/0/UBSQJ3z/LYbVcWYBJR389tLreo/N3grdr8Qclx3v+0WVS0eBRz3FjA/Ht8jIW0UtkDVpyFlsyjFAlWxHpOQOqVefCge4gFcgvcP3CWliUV5xiSkzZIsKhPR75JovEfiJSbRfJ/EQzzCHHcH8HVVjAleItghZHsE+Nz/PimS88RTTAo5ruT5TUnfXLG+AQ/7wT6X4GWtAxOlcX6dIf6b7jmCl+auizcHfhNwUkg7tQLaWQJswLvFSXf/+6JufXA03qcfz/nnHuL3JXBp7j/gLcp0Et4y7nfh/REfUkXjkYT3KdUCPx5nARuBBVX1PRJQ93HgMwr5RVJVYoJ3W4wcvF+4J/jvkw+ANVU0Hq2B3/rn7ubX/Q5Iicf3iH9cO6CTf97V/vedAvYfj7eI0V/wEqZheNd3l/TWQzEfE79OJ39bCfzL/75dVYwH3ghDDnBDyLmPLW48KiAmLfznMxbvf4bc17XAf6yJ8O9NtN4j8RKTaL5P4iQeo4EB/rlPAq7260yvqu+RMMdMonS3Hor5mBDF/Ma5qp3sugK2SQF1zH/TZeGtVvoW0CFgf2oh7aQG1DsWeNp/Mff63/8igj6egvcH40e/DxMJ/rTwFuALvE8MdwOr8KY8VMl4+HXa4K0wtx/vU6KHKeITqSoQk7p4K5Leop+bvDqj8X6J7/frzAeSq2I88JK5dX4sfsC7r2hEo14xHJOMcG2H1OmFd//ln4CtwEjFJOy5M6piPIA3Czj3mzHwHplTVJ2S/C6J5nskXmISzfdJnMTjD3h/e7Px/t6sBa4FqlXV90iY+pMoXbIb8zEhivmNcy5vWp+IiIiIiIhI3KhW0R0QERERERERiTYluyIiIiIiIhJ3lOyKiIiIiIhI3FGyKyIiIiIiInFHya6IiIiIiIjEHSW7IiIiIiIiEneU7IqIiFQCZpZuZi5gyzazDDN70cwuMbNi/802s05mNsnMji2LPouIiFRmSnZFREQql4uB7kA/4HbgJ+BZ4FUzSyhmW52AiYCSXRERqXJqVHQHREREJMh659wXAY/nmdkCYAHwZ2B0xXRLREQktmhkV0REpJJzzi0E/gn83swSAcxsspmtNbMfzGyXmb1hZmfmHmNm6cCT/sPPA6ZHp/j7a5jZeDP7xMx+MrNvzOw+M6tTvs9ORESkbCjZFRERiQ2vALWBLv7j5sD9wC+BdGAHsNLMOvr7Xwam+d/nTo3uDmT5ZU8DtwHzgf7AXcDVwDNl+BxERETKjaYxi4iIxIYv/a9NAZxzw3J3mFl1YCmwCS9hvcE5t9PMtvhVgqZGm9k5wFDgCufcU37x62a2G3jazDo559aX6bMREREpYxrZFRERiQ3mf3UAZtbbzFaY2XfAISAHaA20iaCtvsBBYKE/nbmGmdUAXvX394xu10VERMqfRnZFRERiw3H+1ywz64w3rXkZ3khuFnAYeByI5JrbRkAtYF8B++uXrqsiIiIVT8muiIhIbOgP/AisASbgjeYOcc7l5FYws2OAPRG09Z3f1jkF7P+mVD0VERGpBJTsioiIVHJmNgQYBDzgnNvvr8h8GH9Ks1/nXKAFsDXg0J/8r6H3510KjAOOds4tL7OOi4iIVCAluyIiIpVLJzNrgDfNuAUwAG815deA8X6dpcAfgDlm9iTetbq3A5khbX3sf73OzObiXdf7kXPuTTN7FnjBzGYA7wNHgBSgHzDOOfdZ2Tw9ERGR8mHOuaJriYiISJkKuS8ueNOMdwBr8W4P9IIL+KNtZqOBG4EmwEa8RPg2AOdcakC9icBwv1414HjnXIaZVQNGA1fhLWr1E5CBdx3wn5xzP5TB0xQRESk3SnZFREREREQk7ujWQyIiIiIiIhJ3lOyKiIiIiIhI3FGyKyIiIiIiInFHya6IiIiIiIjEHSW7IiIiIiIiEneU7IqIiIiIiEjcUbIrIiIiIiIicUfJroiIiIiIiMQdJbsiIiIiIiISd/4fqALsBWmf9u0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "final_iso_code = 'CHE'\n",
    "\n",
    "# Final training colums\n",
    "train_cols= swissre_features['demography'] +\\\n",
    "            swissre_features['sanitary'] +\\\n",
    "            swissre_features['economic'] +\\\n",
    "            swissre_features['policies']+ ['shifted_r_estim']# +\\\n",
    "              #swissre_features['weather'] # swissre_features['vaccines']+\\\n",
    "\n",
    "final_train_cols = [col for col in train_cols if not df[df['iso_code'] == final_iso_code][col].isnull().all()]\n",
    "final_train_cols = [col for col in final_train_cols if not df[col].std() == 0]\n",
    "\n",
    "context = Context()\n",
    "pipeline = Pipeline(df, final_train_cols, 'unemployment_rate_idx', final_iso_code, 28, context, 'weight') #shifted_r_estim\n",
    "\n",
    "pipeline.fit_pipeline(save_model=False)\n",
    "\n",
    "predictions = pipeline.predict()\n",
    "\n",
    "test_indices = df['iso_code'] == final_iso_code\n",
    "index = df.loc[test_indices].index\n",
    "ground = df.loc[test_indices]['unemployment_rate_idx']\n",
    "error_curve = np.abs(ground - predictions)\n",
    "\n",
    "df2 = pd.DataFrame({'index': index, 'ground': ground, 'pred': predictions, 'error': error_curve}).sort_values(by = 'index')\n",
    "df2 = df2[(df2['index'] >= '2020-05-01') & (df2['index'].dt.day == 1)]\n",
    "\n",
    "df2['index'] = df2['index'].dt.strftime('%Y-%m')\n",
    "\n",
    "df2.plot.bar(x = 'index', y = ['ground', 'pred'], rot = 0, fontsize=14, legend=False, figsize = (16, 3))\n",
    "\n",
    "pre_patch = mpatches.Patch(color=(1.0, 0.4980392156862745, 0.054901960784313725, 1.0), label='Predicted UER')\n",
    "ground_patch = mpatches.Patch(color=(0.12156862745098039, 0.4666666666666667, 0.7058823529411765, 1.0), label='Reported (and interpolated) UER')\n",
    "plt.legend(handles=[ground_patch, pre_patch], fontsize=14)\n",
    "\n",
    "plt.title(\"Switzerland (CHE)\", fontsize=16)\n",
    "plt.xlabel(\"Date\", fontsize=16);\n",
    "plt.ylabel(\"UER\", fontsize=16);"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}