stan_sirvc_train_sequential.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using backend: pytorch\n"
     ]
    }
   ],
   "source": [
    "import csv\n",
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "#from owid_downloader import GenerateTrainingData\n",
    "#from utils import date_today, gravity_law_commute_dist\n",
    "\n",
    "os.environ['NUMEXPR_MAX_THREADS'] = '16'\n",
    "os.environ['NUMEXPR_NUM_THREADS'] = '8'\n",
    "\n",
    "import pickle\n",
    "import matplotlib.pyplot as plt\n",
    "import dgl\n",
    "import torch\n",
    "from torch import nn\n",
    "import torch.nn.functional as F\n",
    "from model import STAN\n",
    "\n",
    "import sklearn\n",
    "from sklearn.metrics import mean_absolute_error\n",
    "from sklearn.metrics import mean_squared_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#GenerateTrainingData().download_jhu_data('2020-08-01', '2020-12-01')\n",
    "start_date = '2021-01-01'\n",
    "end_date = '2021-05-31'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "163\n"
     ]
    }
   ],
   "source": [
    "# data processing\n",
    "raw_data = pd.read_csv('https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/owid-covid-data.csv', usecols=[\"location\", \"date\", \"total_cases\", \"new_cases_smoothed\", \"total_deaths\",\n",
    "                    \"new_deaths\", \"total_vaccinations\", \"people_fully_vaccinated\", \"new_vaccinations\", \"population\"])\n",
    "raw_data['date'] = pd.to_datetime(raw_data['date'])\n",
    "mask = (raw_data['date'] >= start_date) & (raw_data['date'] <= end_date) # & (raw_data['location'].isin(countries))\n",
    "raw_data = raw_data.loc[mask]\n",
    "#print(raw_data[raw_data['location'] == 'United States']['total_cases'].values[0])\n",
    "countries = []\n",
    "loc_list = list(raw_data['location'].unique())\n",
    "# only include countries that have more than 1000 total cases on start date and at least 1 death\n",
    "for loc in loc_list:\n",
    "    if raw_data[raw_data['location'] == loc][\"total_cases\"].values[0] > 1000 and \\\n",
    "        raw_data[raw_data['location'] == loc][\"total_deaths\"].values[0] > 0:\n",
    "        countries.append(loc)\n",
    "# hard-coded; these are problematic locations (non-countries) that need to be removed\n",
    "to_remove = [\"European Union\", \"Europe\", \"Africa\", \"Asia\", \"North America\", \"Oceania\", \"South America\", \"World\",\n",
    "            \"Tajikistan\", \"High income\", \"Low income\", \"Lower middle income\", \"Upper middle income\"]\n",
    "for r in to_remove:\n",
    "    countries.remove(r)\n",
    "mask = raw_data['location'].isin(countries)\n",
    "raw_data = raw_data.loc[mask]\n",
    "print(len(raw_data['location'].unique()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Generate Graph\n",
    "# add flight neighbors\n",
    "# for now, add a connection if there is any flight between the two countries between start and end date\n",
    "loc_list = list(raw_data['location'].unique())\n",
    "flight_counts = pd.read_csv('processed_flights/flight_counts_2021_all_to_05.csv')\n",
    "adj_map = {}\n",
    "for each_loc in loc_list:\n",
    "    df = flight_counts.loc[flight_counts[\"origin_country\"] == each_loc]\n",
    "    adj_map[each_loc] = set(df[\"destination_country\"].unique())\n",
    "flight_counts['day'] = pd.to_datetime(flight_counts['day'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# add land neighbors\n",
    "import csv\n",
    "neighbor_reader = csv.reader(open('neighbors.csv', 'r'))\n",
    "neighbors = {}\n",
    "for row in neighbor_reader:\n",
    "   neighbors[row[0]] = row[1].split(',')\n",
    "for each_loc,connected in adj_map.items():\n",
    "    for neighbor in neighbors[each_loc]:\n",
    "        if neighbor in loc_list:\n",
    "            connected.add(neighbor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<bound method DGLHeteroGraph.number_of_nodes of Graph(num_nodes=163, num_edges=3765,\n",
      "      ndata_schemes={}\n",
      "      edata_schemes={})>\n"
     ]
    }
   ],
   "source": [
    "# create graph\n",
    "rows = []\n",
    "cols = []\n",
    "for each_loc in adj_map:\n",
    "    for each_loc2 in adj_map[each_loc]:\n",
    "        if each_loc in loc_list and each_loc2 in loc_list:\n",
    "            rows.append(loc_list.index(each_loc))\n",
    "            cols.append(loc_list.index(each_loc2))\n",
    "#print(rows)\n",
    "#print(cols)\n",
    "g = dgl.graph((rows, cols))\n",
    "print(g.number_of_nodes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import networkx as nx\n",
    "# nx_g = g.to_networkx()\n",
    "# pos = nx.kamada_kawai_layout(nx_g)\n",
    "# plt.figure(1,figsize=(8,8)) \n",
    "# nx.draw(nx_g, pos, with_labels=True, node_color=[[.7, .7, .7]])\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "done\n"
     ]
    }
   ],
   "source": [
    "#Preprocess features\n",
    "\n",
    "#active_cases = []\n",
    "confirmed_cases = []\n",
    "new_cases = []\n",
    "new_vaccinations = []\n",
    "fully_vaccinated = []\n",
    "death_cases = []\n",
    "static_feat = []\n",
    "\n",
    "for i, each_loc in enumerate(loc_list):\n",
    "    confirmed_cases.append(raw_data[raw_data['location'] == each_loc]['total_cases'])\n",
    "    new_cases.append(raw_data[raw_data['location'] == each_loc]['new_cases_smoothed'])\n",
    "    new_vaccinations.append(raw_data[raw_data['location'] == each_loc]['new_vaccinations'])\n",
    "    fully_vaccinated.append(raw_data[raw_data['location'] == each_loc]['people_fully_vaccinated'])\n",
    "    death_cases.append(raw_data[raw_data['location'] == each_loc]['total_deaths'])\n",
    "    static_feat.append(np.array(raw_data[raw_data['location'] == each_loc][['population']]))\n",
    "confirmed_cases = np.nan_to_num(np.array(confirmed_cases))\n",
    "death_cases = np.nan_to_num(np.array(death_cases)[:, 14:])\n",
    "new_cases = np.nan_to_num(np.array(new_cases)[:, 14:])\n",
    "new_vaccinations = np.nan_to_num(np.array(new_vaccinations)[:, 14:])\n",
    "fully_vaccinated = np.nan_to_num(np.array(fully_vaccinated))\n",
    "static_feat = np.nan_to_num(np.array(static_feat)[:, 0, :])\n",
    "\n",
    "import copy\n",
    "# active = confirmed(today) - confirmed(14 days ago)\n",
    "cases_copy = copy.deepcopy(confirmed_cases)\n",
    "active = []\n",
    "for loc in confirmed_cases:\n",
    "    active_loc = []\n",
    "    for i in range(14, len(loc)):\n",
    "        active_loc.append(loc[i] - loc[i-14])\n",
    "    active.append(active_loc)\n",
    "active_cases = np.array(active)\n",
    "\n",
    "confirmed_cases = confirmed_cases[:, 14:]\n",
    "fully_vaccinated = fully_vaccinated[:, 14:]\n",
    "\n",
    "recovered_cases = confirmed_cases - active_cases - death_cases + 0.94*fully_vaccinated\n",
    "susceptible_cases = np.expand_dims(static_feat[:, 0], -1) - active_cases - recovered_cases\n",
    "\n",
    "# Batch_feat: new_cases(dI), dR, dS\n",
    "#dI = np.array(new_cases)\n",
    "dI = np.concatenate((np.zeros((active_cases.shape[0],1), dtype=np.float32), np.diff(active_cases)), axis=-1)\n",
    "dR = np.concatenate((np.zeros((recovered_cases.shape[0],1), dtype=np.float32), np.diff(recovered_cases)), axis=-1)\n",
    "dS = np.concatenate((np.zeros((susceptible_cases.shape[0],1), dtype=np.float32), np.diff(susceptible_cases)), axis=-1)\n",
    "# number of new fully vaccinated each day\n",
    "Vt = np.concatenate((np.zeros((fully_vaccinated.shape[0],1), dtype=np.float32), np.diff(fully_vaccinated)), axis=-1)\n",
    "print(\"done\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#Build normalizer\n",
    "normalizer = {'S':{}, 'I':{}, 'R':{}, 'dS':{}, 'dI':{}, 'dR':{}, 'Vt':{}}\n",
    "\n",
    "for i, each_loc in enumerate(loc_list):\n",
    "    normalizer['S'][each_loc] = (np.mean(susceptible_cases[i]), np.std(susceptible_cases[i]))\n",
    "    normalizer['I'][each_loc] = (np.mean(active_cases[i]), np.std(active_cases[i]))\n",
    "    normalizer['R'][each_loc] = (np.mean(recovered_cases[i]), np.std(recovered_cases[i]))\n",
    "    normalizer['dI'][each_loc] = (np.mean(dI[i]), np.std(dI[i]))\n",
    "    normalizer['dR'][each_loc] = (np.mean(dR[i]), np.std(dR[i]))\n",
    "    normalizer['dS'][each_loc] = (np.mean(dS[i]), np.std(dS[i]))\n",
    "    normalizer['Vt'][each_loc] = (np.mean(Vt[i]), np.std(Vt[i]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prepare_data(data, sum_I, sum_R, Vt, edges_df, start, history_window=8, pred_window=14, slide_step=7):\n",
    "    # Data shape n_loc, timestep, n_feat\n",
    "    # Reshape to n_loc, t, history_window*n_feat\n",
    "    n_loc = data.shape[0]\n",
    "    timestep = data.shape[1]\n",
    "    n_feat = data.shape[2]\n",
    "    \n",
    "    x = []\n",
    "    y_I = []\n",
    "    y_R = []\n",
    "    last_I = []\n",
    "    last_R = []\n",
    "    concat_I = []\n",
    "    concat_R = []\n",
    "    concat_Vt = []\n",
    "    edges = []\n",
    "    for i in range(0, timestep, slide_step):\n",
    "        if i+history_window+pred_window-1 >= timestep or i+history_window >= timestep:\n",
    "            break\n",
    "        x.append(data[:, i:i+history_window, :].reshape((n_loc, history_window*n_feat)))\n",
    "        \n",
    "        concat_I.append(data[:, i+history_window-1, 0])\n",
    "        concat_R.append(data[:, i+history_window-1, 1])\n",
    "        last_I.append(sum_I[:, i+history_window-1])\n",
    "        last_R.append(sum_R[:, i+history_window-1])\n",
    "\n",
    "        y_I.append(data[:, i+history_window:i+history_window+pred_window, 0])\n",
    "        y_R.append(data[:, i+history_window:i+history_window+pred_window, 1])\n",
    "\n",
    "        concat_Vt.append(Vt[:, i+history_window:i+history_window+pred_window])\n",
    "        \n",
    "        e_matrix = np.zeros((n_loc, n_loc)) # edge weight matrix for every time period \n",
    "        df = edges_df.groupby([\"origin_country\"])\n",
    "        for loc in range(n_loc):\n",
    "            try:\n",
    "                src_df = df.get_group(loc_list[loc])\n",
    "                src_df = src_df.groupby([\"destination_country\"])\n",
    "                for loc2 in range(n_loc):\n",
    "                    try:\n",
    "                        dst_df = src_df.get_group(loc_list[loc2])\n",
    "                        dst_df = dst_df.loc[(dst_df['day'] >= (pd.to_datetime(start) + pd.DateOffset(days=i)))]\n",
    "                        dst_df = dst_df.loc[(dst_df['day'] < (pd.to_datetime(start) + pd.DateOffset(days=i+history_window-1)))]\n",
    "                        e_matrix[loc, loc2] = dst_df['flight_count'].sum()\n",
    "                    except:\n",
    "                        continue\n",
    "            except:\n",
    "                continue\n",
    "        # normalize matrix (doubly stochastic, see https://arxiv.org/pdf/1809.02709.pdf)\n",
    "        # step 1: row normalize\n",
    "        norm = np.sum(e_matrix, axis=1, keepdims=True)\n",
    "        norm[norm==0] = 1e-10\n",
    "        norm = 1.0 / norm\n",
    "        P = e_matrix * norm\n",
    "\n",
    "        # step 2: P @ P^T / column_norm\n",
    "        norm = np.sum(P, axis=0, keepdims=True)\n",
    "        norm[norm==0] = 1e-10\n",
    "        norm = 1.0 / norm\n",
    "\n",
    "        PT = np.transpose(P, (1, 0))\n",
    "        P = np.multiply(P, norm)\n",
    "        T = np.matmul(P, PT)\n",
    "        edges.append(T) # n_rows = # countries, n_cols = # countries\n",
    "        \n",
    "    \n",
    "    x = np.array(x, dtype=np.float32).transpose((1, 0, 2))\n",
    "    last_I = np.array(last_I, dtype=np.float32).transpose((1, 0))\n",
    "    last_R = np.array(last_R, dtype=np.float32).transpose((1, 0))\n",
    "    concat_I = np.array(concat_I, dtype=np.float32).transpose((1, 0))\n",
    "    concat_R = np.array(concat_R, dtype=np.float32).transpose((1, 0))\n",
    "    y_I = np.array(y_I, dtype=np.float32).transpose((1, 0, 2))\n",
    "    y_R = np.array(y_R, dtype=np.float32).transpose((1, 0, 2))\n",
    "    concat_Vt = np.array(concat_Vt, dtype=np.float32).transpose((1, 0, 2))\n",
    "    return x, last_I, last_R, concat_I, concat_R, y_I, y_R, concat_Vt, edges"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def squish_edges(E):\n",
    "    # use src and dst (rows and cols) to make E matrix\n",
    "    edges = []\n",
    "    for M in E:\n",
    "        edges_flat = []\n",
    "        for i in range(len(rows)):\n",
    "            edges_flat.append(M[rows[i]][cols[i]])\n",
    "        edges.append(edges_flat)\n",
    "    return np.array(edges)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "valid_window = 29\n",
    "test_window = 29\n",
    "\n",
    "history_window=14 # days of information\n",
    "pred_window=14 # predicts future # of days\n",
    "slide_step=3 # increment\n",
    "\n",
    "dynamic_feat = np.concatenate((np.expand_dims(dI, axis=-1), np.expand_dims(dR, axis=-1), np.expand_dims(dS, axis=-1)), axis=-1)\n",
    "    \n",
    "#Normalize\n",
    "for i, each_loc in enumerate(loc_list):\n",
    "    dynamic_feat[i, :, 0] = (dynamic_feat[i, :, 0] - normalizer['dI'][each_loc][0]) / normalizer['dI'][each_loc][1]\n",
    "    dynamic_feat[i, :, 1] = (dynamic_feat[i, :, 1] - normalizer['dR'][each_loc][0]) / normalizer['dR'][each_loc][1]\n",
    "    dynamic_feat[i, :, 2] = (dynamic_feat[i, :, 2] - normalizer['dS'][each_loc][0]) / normalizer['dS'][each_loc][1]\n",
    "    # vaccinations don't need to be normalized\n",
    "    #mean_vax = normalizer['Vt'][each_loc][0]\n",
    "    #if mean_vax != 0:\n",
    "    #   Vt[i] = (Vt[i] - mean_vax) / normalizer['Vt'][each_loc][1]\n",
    "dI_mean = []\n",
    "dI_std = []\n",
    "dR_mean = []\n",
    "dR_std = []\n",
    "\n",
    "for i, each_loc in enumerate(loc_list):\n",
    "    dI_mean.append(normalizer['dI'][each_loc][0])\n",
    "    dR_mean.append(normalizer['dR'][each_loc][0])\n",
    "    dI_std.append(normalizer['dI'][each_loc][1])\n",
    "    dR_std.append(normalizer['dR'][each_loc][1])\n",
    "\n",
    "dI_mean = np.array(dI_mean)\n",
    "dI_std = np.array(dI_std)\n",
    "dR_mean = np.array(dR_mean)\n",
    "dR_std = np.array(dR_std)\n",
    "\n",
    "#Split train-test\n",
    "train_feat = dynamic_feat[:, :-valid_window-test_window, :]\n",
    "val_feat = dynamic_feat[:, -valid_window-test_window:-test_window, :]\n",
    "test_feat = dynamic_feat[:, -test_window:, :]\n",
    "\n",
    "valid_start_date = pd.to_datetime(end_date) + pd.DateOffset(days=-valid_window) + pd.DateOffset(days=-test_window)\n",
    "test_start_date = pd.to_datetime(end_date) + pd.DateOffset(days=-test_window)\n",
    "\n",
    "train_edges = flight_counts[(flight_counts[\"day\"] >= start_date) & (flight_counts[\"day\"] < valid_start_date)]\n",
    "val_edges = flight_counts[(flight_counts[\"day\"] >= valid_start_date) & (flight_counts[\"day\"] < test_start_date)]\n",
    "test_edges = flight_counts[(flight_counts[\"day\"] >= test_start_date) & (flight_counts[\"day\"] < end_date)]\n",
    "\n",
    "train_x, train_I, train_R, train_cI, train_cR, train_yI, train_yR, train_Vt, train_edges = prepare_data(train_feat, active_cases[:, :-valid_window-test_window], recovered_cases[:, :-valid_window-test_window], Vt[:, :-valid_window-test_window], train_edges, start_date, history_window, pred_window, slide_step)\n",
    "\n",
    "val_x, val_I, val_R, val_cI, val_cR, val_yI, val_yR, val_Vt, val_edges = prepare_data(val_feat, active_cases[:, -valid_window-test_window:-test_window], recovered_cases[:, -valid_window-test_window:-test_window], Vt[:, -valid_window-test_window:-test_window], val_edges, valid_start_date, history_window, pred_window, slide_step)\n",
    "\n",
    "test_x, test_I, test_R, test_cI, test_cR, test_yI, test_yR, test_Vt, test_edges = prepare_data(test_feat, active_cases[:, -test_window:], recovered_cases[:, -test_window:], Vt[:, -test_window:], test_edges, test_start_date, history_window, pred_window, slide_step)\n",
    "\n",
    "train_edges = squish_edges(train_edges)\n",
    "val_edges = squish_edges(val_edges)\n",
    "test_edges = squish_edges(test_edges)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(18, 3765)\n",
      "(163, 18, 42)\n",
      "(1, 3765)\n",
      "(163, 1, 42)\n",
      "(1, 3765)\n",
      "(163, 1, 42)\n"
     ]
    }
   ],
   "source": [
    "print(train_edges.shape) # one edge array (len = # edges) for each timestep\n",
    "print(train_x.shape) # one array of features for each timestep for each location\n",
    "print(val_edges.shape)\n",
    "print(val_x.shape)\n",
    "print(test_edges.shape)\n",
    "print(test_x.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Build STAN model\n",
    "\n",
    "in_dim = 3*history_window\n",
    "hidden_dim1 = 32\n",
    "hidden_dim2 = 32\n",
    "gru_dim = 32\n",
    "num_heads = 1\n",
    "device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n",
    "\n",
    "g = g.to(device)\n",
    "model = STAN(g, in_dim, hidden_dim1, hidden_dim2, gru_dim, num_heads, pred_window, device).to(device)\n",
    "optimizer = torch.optim.Adam(model.parameters(), lr=1e-2)\n",
    "criterion = nn.MSELoss()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "STAN(\n",
       "  (layer1): MultiHeadGATLayer(\n",
       "    (heads): ModuleList(\n",
       "      (0): GATLayer(\n",
       "        (fc): Linear(in_features=42, out_features=32, bias=True)\n",
       "        (attn_fc): Linear(in_features=64, out_features=1, bias=True)\n",
       "      )\n",
       "    )\n",
       "  )\n",
       "  (layer2): MultiHeadGATLayer(\n",
       "    (heads): ModuleList(\n",
       "      (0): GATLayer(\n",
       "        (fc): Linear(in_features=32, out_features=32, bias=True)\n",
       "        (attn_fc): Linear(in_features=64, out_features=1, bias=True)\n",
       "      )\n",
       "    )\n",
       "  )\n",
       "  (gru): GRUCell(32, 32)\n",
       "  (nn_res_I): Linear(in_features=34, out_features=14, bias=True)\n",
       "  (nn_res_R): Linear(in_features=34, out_features=14, bias=True)\n",
       "  (nn_res_sir): Linear(in_features=34, out_features=3, bias=True)\n",
       ")"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_x = torch.tensor(train_x).to(device)\n",
    "train_I = torch.tensor(train_I).to(device)\n",
    "train_R = torch.tensor(train_R).to(device)\n",
    "train_cI = torch.tensor(train_cI).to(device)\n",
    "train_cR = torch.tensor(train_cR).to(device)\n",
    "train_yI = torch.tensor(train_yI).to(device)\n",
    "train_yR = torch.tensor(train_yR).to(device)\n",
    "train_Vt = torch.tensor(train_Vt).to(device)\n",
    "train_edges = torch.tensor(train_edges).to(device)\n",
    "\n",
    "val_x = torch.tensor(val_x).to(device)\n",
    "val_I = torch.tensor(val_I).to(device)\n",
    "val_R = torch.tensor(val_R).to(device)\n",
    "val_cI = torch.tensor(val_cI).to(device)\n",
    "val_cR = torch.tensor(val_cR).to(device)\n",
    "val_yI = torch.tensor(val_yI).to(device)\n",
    "val_yR = torch.tensor(val_yR).to(device)\n",
    "val_Vt = torch.tensor(val_Vt).to(device)\n",
    "val_edges = torch.tensor(val_edges).to(device)\n",
    "\n",
    "test_x = torch.tensor(test_x).to(device)\n",
    "test_I = torch.tensor(test_I).to(device)\n",
    "test_R = torch.tensor(test_R).to(device)\n",
    "test_cI = torch.tensor(test_cI).to(device)\n",
    "test_cR = torch.tensor(test_cR).to(device)\n",
    "test_yI = torch.tensor(test_yI).to(device)\n",
    "test_yR = torch.tensor(test_yR).to(device)\n",
    "test_Vt = torch.tensor(test_Vt).to(device)\n",
    "test_edges = torch.tensor(test_edges).to(device)\n",
    "\n",
    "dI_mean = torch.tensor(dI_mean, dtype=torch.float32).to(device).reshape((dI_mean.shape[0], 1, 1))\n",
    "dI_std = torch.tensor(dI_std, dtype=torch.float32).to(device).reshape((dI_mean.shape[0], 1, 1))\n",
    "dR_mean = torch.tensor(dR_mean, dtype=torch.float32).to(device).reshape((dI_mean.shape[0], 1, 1))\n",
    "dR_std = torch.tensor(dR_std, dtype=torch.float32).to(device).reshape((dI_mean.shape[0], 1, 1))\n",
    "\n",
    "N = torch.tensor(static_feat[:, 0], dtype=torch.float32).to(device).unsqueeze(-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_adaptive_loss(pred_I, pred_I_sir, pred_R, pred_R_sir, true_I, true_R, pred_window=14):\n",
    "    total_loss = 0\n",
    "    for timestep in range(len(pred_I)):\n",
    "        for day in range(1, pred_window+1):\n",
    "            sir_weight = day/(pred_window+1)\n",
    "            pred_weight = 1-sir_weight\n",
    "            sir_loss = sir_weight*((pred_I_sir[timestep][day-1] - true_I[timestep][day-1])**2) + sir_weight*((pred_R_sir[timestep][day-1] - true_R[timestep][day-1])**2)\n",
    "            pred_loss = pred_weight*((pred_I[timestep][day-1] - true_I[timestep][day-1])**2) + pred_weight*((pred_R[timestep][day-1] - true_R[timestep][day-1])**2)\n",
    "            total_loss += sir_loss + pred_loss\n",
    "    return total_loss / pred_window\n",
    "\n",
    "def val_adaptive_loss(pred_I, pred_I_sir, true_I, pred_window=14):\n",
    "    total_loss = 0\n",
    "    for day in range(1, pred_window+1):\n",
    "        sir_weight = day/(pred_window+1)\n",
    "        pred_weight = 1-sir_weight\n",
    "        sir_loss = sir_weight*((pred_I_sir[day-1] - true_I[day-1])**2)\n",
    "        pred_loss = pred_weight*((pred_I[day-1] - true_I[day-1])**2)\n",
    "        total_loss += sir_loss + pred_loss\n",
    "    return total_loss / pred_window"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----Save best model-----\n",
      "Epoch 0, Loss 24.06, Val loss 0.33\n",
      "Epoch 1, Loss 24.13, Val loss 0.34\n",
      "-----Save best model-----\n",
      "Epoch 2, Loss 23.92, Val loss 0.34\n",
      "Epoch 3, Loss 23.97, Val loss 0.35\n",
      "Epoch 4, Loss 24.21, Val loss 0.35\n",
      "Epoch 5, Loss 24.25, Val loss 0.35\n",
      "Epoch 6, Loss 24.29, Val loss 0.35\n",
      "Epoch 7, Loss 24.07, Val loss 0.34\n",
      "Epoch 8, Loss 24.36, Val loss 0.34\n",
      "Epoch 9, Loss 24.15, Val loss 0.35\n",
      "-----Save best model-----\n",
      "Epoch 10, Loss 23.91, Val loss 0.35\n",
      "-----Save best model-----\n",
      "Epoch 11, Loss 23.86, Val loss 0.35\n",
      "Epoch 12, Loss 24.05, Val loss 0.36\n",
      "-----Save best model-----\n",
      "Epoch 13, Loss 23.78, Val loss 0.35\n",
      "Epoch 14, Loss 23.95, Val loss 0.36\n",
      "-----Save best model-----\n",
      "Epoch 15, Loss 23.51, Val loss 0.37\n",
      "Epoch 16, Loss 24.66, Val loss 0.35\n",
      "Epoch 17, Loss 23.54, Val loss 0.35\n",
      "Epoch 18, Loss 23.77, Val loss 0.37\n",
      "Epoch 19, Loss 24.36, Val loss 0.37\n",
      "Epoch 20, Loss 24.27, Val loss 0.36\n",
      "Epoch 21, Loss 24.47, Val loss 0.36\n",
      "Epoch 22, Loss 23.58, Val loss 0.35\n",
      "Epoch 23, Loss 23.87, Val loss 0.36\n",
      "Epoch 24, Loss 24.04, Val loss 0.36\n",
      "-----Save best model-----\n",
      "Epoch 25, Loss 23.43, Val loss 0.37\n",
      "Epoch 26, Loss 24.18, Val loss 0.36\n",
      "Epoch 27, Loss 24.42, Val loss 0.36\n",
      "Epoch 28, Loss 23.71, Val loss 0.35\n",
      "Epoch 29, Loss 23.68, Val loss 0.36\n",
      "Epoch 30, Loss 23.65, Val loss 0.36\n",
      "Epoch 31, Loss 24.03, Val loss 0.36\n",
      "Epoch 32, Loss 23.79, Val loss 0.37\n",
      "Epoch 33, Loss 24.21, Val loss 0.37\n",
      "Epoch 34, Loss 23.84, Val loss 0.37\n",
      "Epoch 35, Loss 23.83, Val loss 0.37\n",
      "Epoch 36, Loss 23.85, Val loss 0.37\n",
      "Epoch 37, Loss 24.13, Val loss 0.37\n",
      "Epoch 38, Loss 23.85, Val loss 0.37\n",
      "Epoch 39, Loss 23.73, Val loss 0.37\n",
      "Epoch 40, Loss 23.50, Val loss 0.37\n",
      "Epoch 41, Loss 23.62, Val loss 0.37\n",
      "Epoch 42, Loss 24.14, Val loss 0.36\n",
      "Epoch 43, Loss 23.53, Val loss 0.37\n",
      "Epoch 44, Loss 23.53, Val loss 0.37\n",
      "Epoch 45, Loss 23.89, Val loss 0.37\n",
      "Epoch 46, Loss 23.96, Val loss 0.37\n",
      "Epoch 47, Loss 24.00, Val loss 0.37\n",
      "Epoch 48, Loss 24.25, Val loss 0.37\n",
      "Epoch 49, Loss 23.96, Val loss 0.37\n"
     ]
    }
   ],
   "source": [
    "#Train STAN\n",
    "all_loss = []\n",
    "file_name = './save/stan_sequential'\n",
    "min_loss = 1e20\n",
    "\n",
    "train_locs = [\"United States\", \"United Kingdom\"]\n",
    "\n",
    "for epoch in range(50):\n",
    "    global_train_loss = 0\n",
    "    global_val_loss = 0\n",
    "    for loc in train_locs:\n",
    "        cur_loc = loc_list.index(loc)\n",
    "        model.train()\n",
    "        optimizer.zero_grad()\n",
    "        \n",
    "        active_pred, recovered_pred, phy_active, phy_recover, _ = model(train_x, train_cI[cur_loc], train_cR[cur_loc], N[cur_loc], train_I[cur_loc], train_R[cur_loc], V=train_Vt[cur_loc], e_weights=train_edges)\n",
    "        phy_active = (phy_active - dI_mean[cur_loc]) / dI_std[cur_loc]\n",
    "        phy_recover = (phy_recover - dR_mean[cur_loc]) / dR_std[cur_loc]\n",
    "        # SIR loss = (day) / (pred_window + 1); dynamics loss = 1 - SIR loss \n",
    "        #loss = criterion(active_pred.squeeze(), train_yI[cur_loc].squeeze())+criterion(recovered_pred.squeeze(), train_yR[cur_loc].squeeze()) \\\n",
    "        #    + 0.1*criterion(phy_active.squeeze(), train_yI[cur_loc].squeeze())+0.1*criterion(phy_recover.squeeze(), train_yR[cur_loc].squeeze())\n",
    "        loss = train_adaptive_loss(active_pred.squeeze(), phy_active.squeeze(), recovered_pred.squeeze(), phy_recover.squeeze(), train_yI[cur_loc].squeeze(), train_yR[cur_loc].squeeze(), pred_window)\n",
    "\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        all_loss.append(loss.item())\n",
    "        \n",
    "        model.eval()\n",
    "        _, _, _, _, prev_h = model(train_x, train_cI[cur_loc], train_cR[cur_loc], N[cur_loc], train_I[cur_loc], train_R[cur_loc], V=train_Vt[cur_loc], e_weights=train_edges)\n",
    "        val_active_pred, val_recovered_pred, val_phy_active, val_phy_recover, _ = model(val_x, val_cI[cur_loc], val_cR[cur_loc], N[cur_loc], val_I[cur_loc], val_R[cur_loc], prev_h, V=val_Vt[cur_loc], e_weights=val_edges)\n",
    "        \n",
    "        val_phy_active = (val_phy_active - dI_mean[cur_loc]) / dI_std[cur_loc]\n",
    "        # SIR loss = (day) / (pred_window + 1); dynamics loss = 1 - SIR loss \n",
    "        # change loss here \n",
    "        #val_loss = criterion(val_active_pred.squeeze(), val_yI[cur_loc].squeeze()) + 0.1*criterion(val_phy_active.squeeze(), val_yI[cur_loc].squeeze())\n",
    "        val_loss = val_adaptive_loss(val_active_pred.squeeze(), val_phy_active.squeeze(), val_yI[cur_loc].squeeze(), pred_window)\n",
    "        global_train_loss += loss\n",
    "        global_val_loss += val_loss\n",
    "    if (global_val_loss + global_train_loss) / 2 < min_loss:   \n",
    "        state = {\n",
    "            'state': model.state_dict(),\n",
    "            'optimizer': optimizer.state_dict(),\n",
    "        }\n",
    "        torch.save(state, file_name)\n",
    "        min_loss = (global_val_loss + global_train_loss) / 2\n",
    "        print('-----Save best model-----')\n",
    "    print('Epoch %d, Loss %.2f, Val loss %.2f'%(epoch, global_train_loss.item(), global_val_loss.item()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated beta in SIR model is 0.00\n",
      "Estimated gamma in SIR model is 0.01\n",
      "Estimated theta in SIR model is 0.00\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEOCAYAAACjJpHCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA84UlEQVR4nO3debzU8/7A8de7fdFeKCcVxW1RUbK1EBFFZe26soTsl2vNHl3XLltKthJJonIjft0IkVJatKAkSqGitJ10znn//nh/pzMdp3Omzpn5zsx5Px+PeZyZ78x35n3mNL3ns70/oqo455xze6pU2AE455xLbZ5InHPOFYknEuecc0XiicQ551yReCJxzjlXJJ5InHPOFYknEueiiMhQEbmzGJ9vuYicsIfnLhSRY4srljzPrSLSOB7P7UoeTyQureT3H6SIDBCRV2I5X1UvV9WBwXnHisjKeMQZPP9wEfl31O3mIrJaRG4IYmmuqlPj9frOFRdPJM4lARFpDXwI3Keqj4YcjnO7xROJK1EirQwRuUFEfg1aABdF3T9cRP4tIpWBSUA9EdkUXOqJSCkR6S8i34nIOhEZIyI1o87vIyI/BPfdHmNM7YD/Abep6tNRx3d0iwWtqjEi8rKIbAy6vdpGPfYwEZkT3PeGiLyep7VzU/C7rhKRvnlev1rwvGuC2O8QkVLBfReKyKciMkhE1ovIMhE5Oji+IngPL9jdv4NLL55IXEm0L1AN2A+4GBgsIjWiH6Cqm4GTgVWquldwWQX8E+gJdALqAb8DgwFEpBkwBOgT3FcLyCgklnbAe8C/VPX5Qh57GjAaqA68DTwdvG45YBwwHKgJvAb0ipwkIl2BG4EuQBMg75jNU8H7cUDwe50PXBR1/xHA/OD3GRXEcDjQGDgPeFpE9iokdpfGPJG4kmg7cK+qblfVd4FNwMExnnsZcLuqrlTVbcAA4EwRKQOcCUxU1Y+D++4Ecgp5viOBDVjrpzDTVPVdVc0GRgKtop6jDPBk8Du9BcyMOu9s4CVVXRAkyAGRO0SkNHAOcKuqblTV5cCjWDKM+F5VXwpe93WgPvb+bVPV/wP+xJKKK6E8kbh0kw2UzXOsLJY8ItapalbU7S1ArN+oGwDjgm6e9cDi4DX3wVohKyIPDP7TXlfI8w0GvgAm520V5ePnPDFXCBJYPeAn3bkC64qo6/Xy3P4h6nptoFyeYz9grbWIX6KubwVQ1bzHvEVSgnkicenmR6BhnmON2Pk/yljlVxp7BXCyqlaPulRQ1Z+A1di3dQBEpBLWHVSQbOAfQdzvi0jVPYhzNbCfiEjUsfp57o++vX/U9bVYkm2Q5/6f9iAOV0J5InHp5nXgDhHJCAbGTwBOBcbuwXP9AtQSkWpRx4YC94lIAwARqSMiPYL7xgLdRaR9MG5xLzF8xlR1O3AW9p/6u8FA/+6YjiWkq0WkTBBPu6j7xwAXikizILndHfXa2cH994lIleD3uh6Iabq0c+CJxKWfe4HPgGnYQPhDwD9UdcHuPpGqfo0NXC8LurLqAU9gA93/JyIbgc+xwWhUdSFwFTYgvTp4/ZjWoajqn8DpQCbwXxGpuBtxRs69GFiPDYBPBLYF908CHgc+AJYGP6NdA2wGlmHv2yjgxVhf3znxja2cSz8iMgMYqqovhR2LS3/eInEuDYhIJxHZN+jaugBoiU0rdi7uyoQdgHOuWByMjXXsBXwHnKmqq8MNyZUU3rXlnHOuSLxryznnXJGUuK6t2rVra8OGDcMOwznnUsrs2bPXqmqd/O4rcYmkYcOGzJo1K+wwnHMupYjILhf1eteWc865IvFE4pxzrkg8kTjnnCuSEjdGkp/t27ezcuVKMjMzww7FARUqVCAjI4OyZfMW8XXOJSNPJMDKlSupUqUKDRs2ZOcCqi7RVJV169axcuVKGjVqFHY4zrkYeNcWkJmZSa1atTyJJAERoVatWt46dC6FeCIJeBJJHv63cK54rVwJzzwD334bn+f3ri3nnEszqrBgAUyYAOPHw+zZdnzQIDjooOJ/PW+RJIlffvmFc889lwMOOIA2bdpw1FFHMW7cuITGsHz5clq0aJHv8VGjRu3Rcz7++ONs2bJlx+299vIdWZ2Lh6ws+OgjuP56aNwYWraEO++EsmXh/vth8WK47rr4vLYnkiSgqvTs2ZOOHTuybNkyZs+ezejRo1m58q97ImVlZeXzDPFVUCIpLJ68icQ5V3w2b7YWx4UXwr77wrHHwuDB8Le/wbPPwqpVMH069O9vx+LFu7aSwAcffEC5cuW4/PLLdxxr0KAB11xzDQDDhw/nnXfeITMzk82bNzN27Fj69u3LsmXLqFSpEsOGDaNly5YMGDCAvfbaixtvvBGAFi1aMHHiRABOPvlk2rdvz2effcZ+++3HhAkTqFixIrNnz6Zv375UqlSJ9u3b5xtf//79Wbx4Ma1bt+aCCy6gRo0aO8Vz11138cgjj+x4rauvvpq2bdvyxx9/sGrVKo477jhq167Nhx9+CMDtt9/OxIkTqVixIhMmTGCfffaJ23vrXLr59VeYONESyOTJkJkJ1atD9+7QowecdBJUqZLYmDyR5HXddTB3bvE+Z+vW8Pjju7x74cKFHHbYYQU+xfTp05k/fz41a9bkmmuu4dBDD2X8+PF88MEHnH/++cwtJOYlS5bw2muv8dxzz3H22Wfz5ptvct5553HRRRfx1FNP0alTJ2666aZ8z33ggQd2ShTDhw/fKZ6pU6fme94///lPHnvsMT788ENq164NwObNmznyyCO57777uPnmm3nuuee44447CozduZJuyZLc8Y7PPrMxkP33h379LHl06GBdWGHxRJKErrrqKqZNm0a5cuX44osvAOjSpQs1a9YEYNq0abz55psAdO7cmXXr1rFhw4YCn7NRo0a0bt0agDZt2rB8+XI2bNjA+vXr6dSpEwB9+vRh0qRJMcUYHc/uKFeuHN27d98Rx+TJk3f7OZxLdzk58MUXuclj8WI73ro13HUX9OwJrVpBskxw9ESSVwEth3hp3rz5jsQAMHjwYNauXUvbtm13HKtcufKO6/ltRiYilClThpycnB3HotdilC9ffsf10qVLs3XrVlR1j6faRsdT0OvmVbZs2R2vWbp06VDGfJxLRllZ1lU1YQK8/TasXg2lS0OnTnD55dbyaNAg7Cjz54PtSaBz585kZmYyZMiQHccKGqDu2LEjr776KgBTp06ldu3aVK1alYYNG/Lll18C8OWXX/L9998X+LrVq1enWrVqTJs2DWDHc+ZVpUoVNm7cuMvnadCgAYsWLWLbtm1s2LCBKVOmxHyucw4WLoQjjoBTToFXX4VjjoGRI208ZMoU+Oc/kzeJgLdIkoKIMH78eP71r3/x0EMPUadOHSpXrsyDDz6Y7+MHDBjARRddRMuWLalUqRIjRowA4IwzzuDll1+mdevWHH744RwUw4Txl156acdg+0knnZTvY1q2bEmZMmVo1aoVF154ITVq1Njp/vr163P22WfTsmVLmjRpwqGHHrrjvn79+nHyySdTt27dHYPtzjmTnW2dILffDlWrWhI5/XSoUCHsyHZPiduzvW3btpp3Y6vFixfTtGnTkCJy+fG/iUt3339v03Y//ti6rYYNg733DjuqXROR2araNr/7vGvLOecSSBWef94WDM6dC8OHw7hxyZ1ECuNdW845lyCrV8Oll8I770DnzvDSSzaNN9V5i8Q55xLgjTegRQsbPH/iCZuhlQ5JBDyROOdcXP32G5x7Lpx9Nhx4IMyZY7OwSqXR/75p9Ks451xyef99OOQQa43ce6+tSo9nzauweCJxzrlitmkTXHEFdO1qdbBmzLBKvGXSdFTaE0mSKF26NK1bt6ZFixacddZZSVsxd1el5p1z5tNPrZTJs8/CjTfaXiCFlNJLeZ5IkkTFihWZO3cuCxYsoFy5cgwdOnSn+7Ozs0OKzDkXi23b4JZbrIBiTg5MnQoPP5x6iwv3hCeSJNShQweWLl3K1KlTOe644zj33HM55JBDyMzM5KKLLuKQQw7h0EMP3bFSfPjw4fTo0YOuXbty8MEHc8899+x4rscee4wWLVrQokULHg/qiG3evJlu3brRqlUrWrRoweuvvw7A7Nmz6dSpE23atOGkk05i9erVO463atWKo446isGDByf2zXAuBcydC4cfDg89BJdcAvPmQceOYUeVOGnaY7fnQqgiv5OsrCwmTZpE165dAZg5cyYLFiygUaNGPProowB89dVXfP3115x44ol8G2zCHHlcpUqVOPzww+nWrRsiwksvvcSMGTNQVY444gg6derEsmXLqFevHu+88w4AGzZsYPv27VxzzTVMmDCBOnXq8Prrr3P77bfz4osvxlRq3rmSKCvLkseAAVCrlu0T0q1b2FElnieSJLF169YdZd47dOjAxRdfzGeffUa7du1o1KgRYOXjI5td/e1vf6NBgwY7EkmXLl2oVasWAKeffjrTpk1DROjVq9eOSr2nn346n3zyCV27duXGG2/klltuoXv37nTo0IEFCxawYMECunTpAlhXWt26dYtUat65dLZkCZx/Pnz+uU3tfeYZSyYlkSeSPEKoIg/kjpHkVVj5+Ii85eBFZJePP+igg5g9ezbvvvsut956KyeeeCK9evWiefPmTJ8+fafHrl+/fo9LzTuXjnJyYMgQuOkmG/947TXo3TvsqMLlYyQpJLp8/LfffsuPP/7IwQcfDMDkyZP57bff2Lp1K+PHj+eYY46hY8eOjB8/ni1btrB582bGjRtHhw4dWLVqFZUqVeK8887jxhtv5Msvv+Tggw9mzZo1OxLJ9u3bWbhwYcyl5p0rCVassK1sr77a9glZsMCTCHiLJKVceeWVXH755RxyyCGUKVOG4cOH79iwqn379vTp04elS5dy7rnn7tgU68ILL6Rdu3YAXHLJJRx66KG8//773HTTTZQqVYqyZcsyZMgQypUrx9ixY/nnP//Jhg0byMrK4rrrrqN58+YxlZp3Lp1lZcHgwbYWJCcHhg61bW69sW68jDypX7J8+PDhzJo1i6effjrsUIpNqv9NXPqYPt0WF86bZ62RwYOt1ElJ42XknXNuN61da1N5jz7aro8dC5MmlcwkUpi4JxIRKS0ic0RkYnB7gIj8JCJzg8spUY+9VUSWisg3InJS1PE2IvJVcN+TEoz+ikh5EXk9OD5DRBrG+/dJRhdeeGFatUacC1NODjz3HBx8MIwYYYPqX38NZ5zhXVm7kogWybXA4jzHBqlq6+DyLoCINAN6A82BrsAzIlI6ePwQoB/QJLh0DY5fDPyuqo2BQUD+e9PGoKR18SUz/1u4sMyZYy2Qfv2s5PvcubZOZK+9wo4sucU1kYhIBtANeD6Gh/cARqvqNlX9HlgKtBORukBVVZ2u9j/My0DPqHNGBNfHAsfLHsxVrVChAuvWrfP/wJKAqrJu3ToqlIS6Ei5pbNhgpd3btrUtcF9+2UqcNG8edmSpId6zth4Hbgaq5Dl+tYicD8wCblDV34H9gM+jHrMyOLY9uJ73OMHPFQCqmiUiG4BawNroFxORfliLhv3z2UkmIyODlStXsmbNmt3/DV2xq1ChAhkZGWGH4UoAVRg1Cm64AdassUH1f//bKva62MUtkYhId+BXVZ0tIsdG3TUEGAho8PNRoC+QX0tCCzhOIfflHlAdBgwDm7WV9/6yZcvuWD3unCsZFi2Cq66ylke7drb9bZs2YUeVmnara0tESolI1RgffgxwmogsB0YDnUXkFVX9RVWzVTUHeA5oFzx+JVA/6vwMYFVwPCOf4zudIyJlgGrAb7vzOznnSpbNm6F/f2jVyqb0PvusTfH1JLLnCk0kIjJKRKqKSGVgEfCNiBRauU9Vb1XVDFVtiA2if6Cq5wVjHhG9gAXB9beB3sFMrEbYoPpMVV0NbBSRI4Pxj/OBCVHnXBBcPzN4DR/ocM79hSqMGwdNm8KDD0KfPvDNNzawnk7b3oYhlrevmar+gQ1wvwvsD/Qpwms+FEzlnQ8cB/wLQFUXAmOwZPUecJWqRjbhuAIbsF8KfAdEqga+ANQSkaXA9UD/IsTlnEtTy5ZB9+5w+ulQowZMmwYvvgh16oQdWXqIZYykrIiUxRLJ06q6XUR261u/qk4FpgbXd5mEVPU+4L58js8C/rItn6pmAmftTizOuZIjM9Om7/7nP1C2LAwaZHWy0nXL27DE8nY+CywH5gEfi0gD4I94BuWcc0X13nuWNL77zgorPvoo1KsXdlTpqdCuLVV9UlX3U9VT1PyAdUk551zSWbECzjwTTj4ZSpeGyZOt1LsnkfiJZbB9HxF5QUQmBbebkTvA7ZxzSeGPP2ynwqZN4d134b77YP58OOGEsCNLf7EMtg8H3gci+fxb4Lo4xeOcc7slMxMeewwOOADuuQe6doWFC+G22yDYZcHFWSyJpLaqjgFywFaQA9kFn+Kcc/GVlWXFFZs0sZXpbdrAF19YlV5fX5xYsSSSzSJSi2DFuIgcCWyIa1TOObcLOTnw+uvQrJmtAcnIgA8/hPfft1pZLvFimbV1Pbbw70AR+RSogy3+c865hFG1/UBuv92q8rZoARMmwKmnenn3sBWaSFT1SxHpBByM1bb6RlW3xz0y55wLTJsGt95qPxs1gpEj4e9/t1lZLnyxzNo6C6gYrDzvCbwuIofFOzDnnJs7F7p1gw4dYOlSeOYZ22TqvPM8iSSTWMZI7lTVjSLSHjgJ2/9jSHzDcs6VZEuWWIvj0EOtoOIDD9jCwiuugHLlwo7O5RVLIonM0OoGDFHVCYD/KZ1zxW7lShtAb9oU3n7bxkOWLYNbboFKlcKOzu1KLIPtP4nIs8AJwIMiUp7EbNHrnCsh1q6F+++HwYNtVtaVV1oS2WefsCNzsYglkZyN7ZH+iKquD8rAF1pG3jnnCrNxoy0mfPRR2yekTx9bnd6wYdiRud0Ry6ytLcBbIrK3iET2qf06vmE559JZZiYMGWJVedeutfLuAwfa2hCXemKZtXWaiCwBvgc+Cn5OKvgs55z7q+3b4fnn4aCD4PrroXVrmDkT3nzTk0gqi2WsYyBwJPCtqjbCxko+jWtUzrm0kpUFI0bA3/4Gl14KdevClClWmffww8OOzhVVLIlku6quA0qJSClV/RBoHd+wnHPpIDsbXn3VWhsXXgjVq8PEifD559C5c9jRueISy2D7ehHZC/gYeFVEfgWy4huWcy6V5eRY8cQBA2DxYmjZ0vZL79HDy5mko1haJD2ALdje6u9he6afGs+gnHOpSRXeegtatYJzzrGkMWYMzJkDPXt6EklXu0wkItJYRI5R1c2qmqOqWao6ApgLVE9UgM655KcK//2vlXI/4wz4808YNco2ljrrLCjlK8/SWkF/3seBjfkc3xLc55wr4VRtb/QjjoDTToMNG2xQfeFCL6pYkhSUSBqq6vy8B1V1FtAwbhE555Keqs26at/e9kb/9Vd44QUrqHj++VAmltFXlzYKSiQVCrivYnEH4pxLDR99BMcea3uh//gjDB0K334LfftC2bJhR+fCUFAi+UJELs17UEQuBmbHLyTnXDL67DNLHscea9V5n3rKSrtfdplX5C3pCmqAXgeME5F/kJs42mKVf3vFOS7nXJKYORPuusu2st17bxg0yJJHRe+XcIFdJhJV/QU4WkSOA1oEh99R1Q8SEplzLlRffgl3320LCGvVgocesqq8lSuHHZlLNrEUbfwQ+DABsTjnksCCBdYCGTcOatSwwopXXw1VqoQdmUtWPrfCOQfYeMfdd8Nrr1nSuOceuPZaqFYt7MhcsvNE4lwJt2KFlXB/8UUoX952I7zpJqhZM+zIXKqIKZGISAOgiar+T0QqAmVUNb/Fis65FPHrr9ZtNWSI3b7ySrjtNth333Djcqmn0EQSTAHuB9QEDgQygKHA8fENzTkXD7//Do88Ak88YRtMXXCBjYk0aBB2ZC5VxdIiuQpoB8wAUNUlIrJ3XKNyzhW7TZsseTzyCKxfD7172zjIQQeFHZlLdbGUUtumqn9GbohIGUBjfQERKS0ic0RkYnC7pohMFpElwc8aUY+9VUSWisg3InJS1PE2IvJVcN+TIlZDVETKi8jrwfEZItIw1ricKykyM+Hxx+GAA+COO6BDB5g71wbVPYm44hBLIvlIRG4DKopIF+AN4L+78RrXAoujbvcHpqhqE2BKcBsRaQb0BpoDXYFnRCRS8m0I1r3WJLh0DY5fDPyuqo2BQcCDuxGXc2lt+3Z47jlo0gT+9S845BCYPh3eftvKvDtXXGJJJP2BNcBXwGXAu8AdsTy5iGQA3YDnow73AEYE10cAPaOOj1bVbar6PbAUaCcidYGqqjpdVRV4Oc85kecaCxwfaa04V1JFdiVs2hT69YOMDCuwOGUKHHlk2NG5dBTLGEkP4GVVfW4Pnv9x4GYgeinTPqq6GkBVV0eNt+wHfB71uJXBse3B9bzHI+esCJ4rS0Q2ALWAtdFBiEg/rEXD/vvvvwe/hnPJTxUmTIA777RFhS1b2h4h3br5hlIuvmJpkZwGfCsiI0WkWzBGUigR6Q78qqqxFnjM75+6FnC8oHN2PqA6TFXbqmrbOnXqxBiOc6lBFf7v/2xPkF69bFOp0aNtV8Lu3T2JuPgrNJGo6kVAY2xs5FzgOxF5vuCzADgGOE1ElgOjgc4i8grwS9BdRfDz1+DxK4H6UednAKuC4xn5HN/pnCDBVQN+iyE259LCp5/CccfBSSfBL7/YniALF9o2t74roUuUmP6pqep2YBKWEGZj3V2FnXOrqmaoakNsEP0DVT0PeBu4IHjYBcCE4PrbQO9gJlYjbFB9ZtANtlFEjgzGP87Pc07kuc4MXiPmGWXOpao5c6zLqn1720zqqady9wTxTaVcosWyILErlgiOA6ZiA+dnF+E1HwDGBPua/AicBaCqC0VkDLAIyAKuUtXs4JwrgOHYhlqTggvAC8BIEVmKtUR6FyEu55Le11/b4sE33rCCig88YAUVvSKvC5MU9gVeREZjLZFJqrotIVHFUdu2bXXWrFlhh+HcbvnhB1s8OGIEVKpk03lvuMELKrrEEZHZqto2v/tiKSPv3/KdC8nPP8N998Gzz9qYx3XXQf/+4HNGXDLZZSIRkWmq2l5ENrLzTCgBVFWrxj0650qo336Dhx+2kiZ//gkXX2zTejMyCj/XuUQraIfE9sFP387GuQTZuNGSx8MP2/W//926tBo3Djsy53at0FlbIjIylmPOuT0XqYd14IHW8jjuOJg3z1aoexJxyS6W6b/No28E6zXaxCcc50qW7dvh+edz62G1bGn1sMaPt9pYzqWCXSaSoBLvRqCliPwRXDYCv5C7jsM5twdycqz6brNmcOmlsN9+8L//2cXrYblUs8tEoqr3B+MjD6tq1eBSRVVrqeqtCYzRubShavWvDj0Uzj0XKlSw+ljTp8PxvlWcS1GxdG3NFJEds9VFpLqI9IxfSM6lpw8+gKOOgtNOgy1bYNQoGwc57TSvh+VSWyyJ5G5V3RC5oarrgbvjFpFzaWbGDDjhBGtxrFwJw4bBokU2I8vrYbl0EMs/4/we49V8nCvEvHnQs6eNecybB4MGwdKlNiZStmzY0TlXfGJJCLNE5DFgMLYw8RqscKNzLo9ffrFB9JEj4csvoWpVGDgQrr0WqviKLJemYkkk1wB3Aq9jq9rfB26PZ1DOpZItW2y67siRMHmy7VB42GHWAjn/fKhZM+wInYuvWGptbSbYVx1ARPYHrgIejmNcziW17GyYOtWSx5tvwqZNUL8+3HwznHeeTet1rqSIdbfD2li5979j29uOi2dQziWrr76y5DFqFPz0k3VdnXOOJY+OHX3w3JVMBRVtrAL0wnZFPAhLHgeoqpeNcyXKqlWWOEaOhPnzbeOok0+Gxx6DU0+FihXDjtC5cBXUIvkVmAncAUxTVRWRXokJy7lwbdoEb71lyWPKFFtIeMQRthPhOed4GXfnohWUSG7DdhwcAowSkdcTE5Jz4cjKshIlI0fa4PmWLdCoEdxxh3VdHXRQ2BE6l5wKKiM/CBgkIgdgYyPjgXoicgswTlW/TUyIzsWPqu1//sor1n31yy+2hW2fPnY5+mhfde5cYWKZtbUMuA+4T0QOwZLKJODAOMfmXNxs327b1j7+OCxcaAsEu3e35HHKKVC+fNgROpc6dmuFuqp+BXyFdXs5l3Kysqz1MXAgLFsGbdrAkCFw9tm+3sO5PeWlTlyJkJ1tK87vucfKlBx2mFXh7dbNu66cKyqf9e7SWk4OjB4NLVpYt1WlSjaQPmuWdWV5EnGu6Ara2GpK8PPBxIXjXPHIyYGxY23Hwb//HUqXtttz5kCPHp5AnCtOBXVt1RWRTsBpIjIaq7O1g6p+GdfInNsDqrZR1N132+LBv/3NWiRnneWrzp2Ll4ISyV1Yja0M4LE89ynQOV5BObe7VOGddyyBfPml7YH+yivQu7e1Rpxz8VPQOpKxwFgRuVNVByYwJudipgrvvw933QVffAEHHADDh8M//mGlTJxz8RfLOpKBInIa0DE4NFVVJ8Y3LOcKpmqlS+66y/Y7b9AAnn/eyrb7plHOJVahvcYicj9wLbAouFwbHHMuFFOnQqdO0KULrFgBQ4fCt9/CxRd7EnEuDLE0/rsBrVU1B0BERgBzgFvjGZhzeU2bZi2QDz+EevXg6afhkkt8FbpzYYt1Hkv1qOvV4hCHc7s0fTqceCJ06ACLFllZk6VL4aqrPIk4lwxiaZHcD8wRkQ+xKcAd8daIS4BNm6BvX3jjDSvb/sgjcMUVtqjQOZc8Yhlsf01EpgKHY4nkFlX9Od6BuZJt5UrbNGr+fLj3Xrj+eqhcOeyonHP5ialrS1VXq+rbqjoh1iQiIhVEZKaIzBORhSJyT3B8gIj8JCJzg8spUefcKiJLReQbETkp6ngbEfkquO9JEVuXLCLlReT14PgMEWm4W7+9S0pz5tgmUkuXwsSJcOednkScS2bxXOu7Deisqq2A1kBXETkyuG+QqrYOLu8CiEgzbCOt5kBX4BkRiSwlGwL0A5oEl67B8YuB31W1MTAI8HIuKe7tt6F9e1tE+OmntqWtcy65xS2RqNkU3CwbXLSAU3oAo1V1m6p+DywF2olIXaCqqk5XVQVeBnpGnTMiuD4WOD7SWnGpRRUGDYKePaFZM5gxw+pkOeeSX4GJRERKiciCPX1yESktInOx/d8nq+qM4K6rRWS+iLwoIjWCY/sBK6JOXxkc2y+4nvf4TueoahawAaiVTxz9RGSWiMxas2bNnv46Lk6ysmwG1vXXQ69e8NFHULdu2FE552JVYCIJ1o7ME5H99+TJVTVbVVtj9braiUgLrJvqQKy7azXwaPDw/FoSWsDxgs7JG8cwVW2rqm3r1KmzW7+Di68//rBy7kOGwM032wwtn5XlXGqJZfpvXWChiMwENkcOquppsb6Iqq4PZn51VdVHIsdF5DkgUm5lJVA/6rQMYFVwPCOf49HnrBSRMtgal99ijcuF64cfLIl8/TUMGwaXXhp2RM65PRFLIrlnT55YROoA24MkUhE4AXhQROqq6urgYb2ASNfZ28AoEXkMqIcNqs9U1WwR2RgM1M8AzgeeijrnAmA6cCbwQTCO4pLczJlw2mmQmQmTJsEJJ4QdkXNuT8WyjuQjEWkANFHV/4lIJSCWwtx1gRHBzKtSwBhVnSgiI0WkNdYFtRy4LHidhSIyBqvnlQVcparZwXNdAQwHKgKTggvAC8BIEVmKtUR6xxCXC9nYsbZbYd26Vu6kadOwI3LOFYUU9gVeRC7Fpt7WVNUDRaQJMFRVj09EgMWtbdu2OmvWrLDDKJFU4aGHoH9/OOoo2/J2773Djso5FwsRma2qbfO7L5bpv1cBxwB/AKjqEsA//m63/PmnjYH072+bTX3wgScR59JFLIlkm6r+GbkRDGr7OISL2e+/28LCF16wVeqvvgoVKoQdlXOuuMQy2P6RiNwGVBSRLsCVwH/jG5ZLF999B926wbJlMGKEbTzlnEsvsbRI+gNrgK+wgfF3gTviGZRLD59+CkceCWvWwP/+50nEuXQVy6ytnGAzqxlYl9Y3PsXWFWbUKLjoItsC9513oEmTsCNyzsVLLFvtdgO+A54EngaWioiX0nP5UrWy7//4h7VGpk/3JOJcuoula+tR4DhVPVZVOwHHYZV2ndvJtm3WfXX33fbz//4Pav2l8plzLuG2bIGXXoIff4zL08eSSH5V1aVRt5dhRRid22HtWlud/sor8O9/w/Dhvg2uc6H79lurhpqRYduNjhkTl5fZ5RiJiJweXF0oIu8CY7AxkrOAL+ISjUtJ33xjM7NWroTRo+Gcc8KOyLkSLCvLdoR75hmYPBnKlIHTT4crr4SOHePykgUNtp8adf0XoFNwfQ1Q468Pd+kqOxt++cUSRX6XOXOgYkUrd3LUUWFH61wJ9fPP8Pzz8Oyz9sHMyICBA+Hii+O+L8MuE4mqXhTXV3ZJISsLVq+GFSt2nShWrbJkEq18eft3mpEBZ51lCw0bNQrnd3CuxFKFjz+2fRjefNM+0F26wFNPWWntMrEsFSy6Ql9FRBoB1wANox+/O2XkXXjWrYMlS+D77/NPEj//DDk5O59TqRLUr29JonPn3IQRfalVC3wvSudC8scfMHKkdV8tWgTVq8M118Dll8NBByU8nFjS1Xisyu5/gZyCH+rCsH69JYv8Lr//vvNjq1bNTQYtWuSfJKpX9yThXFKaP99aHyNHwubN0KaN1R7q3TvUHeFiSSSZqvpk3CNxBdq0adfJInr3YBFrTTRpYoPeTZrYF5QDDrAkUbVqeL+Dc24PbNsGb71lrY9p06xQXe/eNnh++OFhRwfElkieEJG7gf8DtkUOquqXcYuqhNq6FZYuzT9ZrF6982Pr1bMk0aOH/YxOGBUrhhO/c64Y/fCDDZw//7x9W2zcGB55BC68MOkWaMWSSA4B+gCdye3a0uC2K6KsLHj5ZXj4YdtyNtree1uCOOmk3GTRpIn9e9prr3Didc7FUU6OreR95hmrLQRw6qnW+jjhBCgVy9K/xIslkfQCDoguJe+KLjsbXnsN7rnHWiFt2tj16IRRrVrYUTrnEmLtWlt5PnSolcree2+49Vbo1w/23z/s6AoVSyKZB1THV7MXi5wc22p2wABYvBhatrSdAk87zQe4nStRVOHzz6318cYbNhbSoQPcd58tICxXLuwIYxZLItkH+FpEvmDnMRKf/rsbVGHCBKtDNX++7VM+ZgyccUbStladc/GwaZPt7jZkCMybB1WqwCWX2NTdFi3Cjm6PxJJI7o57FGlMFSZNgrvugtmzbXzjlVds0kXp0mFH55xLmIULLXm8/DJs3AitWllX1j/+kfKDnrHsR/JRIgJJN6owZYolkOnToWFDePFF6NMnYYtNnXNhi0zdHTIEPvnEuqvOPtsGz488Mm36s2NZ2b6R3D3aywFlgc2q6isSduHjj61kyMcf29qNoUNtk6cU6vJ0zhXF8uUwbFju1N0DDoCHHrL/CGrXDju6YhdLi6RK9G0R6Qm0i1dAqezzzy2B/O9/sO++8OSTcOmltn7IOZfmsrPhvfes9fHuu9ba6N7dWh9duqT1YOhud7Ko6ngR6R+PYFLV7NnWhfXuu/Zl45FH4IorQq1Y4JxLlF9/tX7rZ5+1lsi++8Ltt9u3yBSYulscYunaOj3qZimgLbldXSVLZuZOzYuvvrIEMn481KgB//mP1U1L8XEz51xhVK1cyZAhNp9/+3Y49ljrvurZE8qWDTvChIqlRRK9L0kWsBzoEZdoktmzz8KDD8KcOXy9uhoDBtj03SpVbE3Iddf5AkLn0t4ff9i0yyFDYMEC+9BfcYVN3W3aNOzoQhPLGInvSwJw6KEsXV6Ge9su4tVlR1Gxoi08veEGqFkz7OCcc3GlahV3r7/e9mY47DAbSO/dGypXDju60BW01e5dBZynqjowDvEkradmtONfLKbc0m1cf+Zybn6mIXXqhB2Vcy7uli2Dyy6zWTRHHw2PPQbt2qXN1N3iUNA0gs35XAAuBm6Jc1xJ54gj4KorlO8OPImHZ3aiToWNYYfknIunrCyrptqiBcyYAYMH21qQI47wJJJHQVvtPhq5LiJVgGuBi4DRwKO7Oi9dtWsH7dqVgc8ehPbt4ZZbrEaOcy79fPmllS2ZM8f2anj6aVsU5vJV4MRmEakpIv8G5mNJ5zBVvUVVS24Bx6OPhn/9ywbbPvgg7Gicc8Vpyxa46Sb75rh6tc3IGjfOk0ghdplIRORh4AtgI3CIqg5Q1d939fh8zq8gIjNFZJ6ILBSRe4LjNUVksogsCX7WiDrnVhFZKiLfiMhJUcfbiMhXwX1Pili7UkTKi8jrwfEZItJw99+CPTBwoBXNuuQSK8DmnEt9kydbN9Yjj0Dfvlae+4wzvBsrBgW1SG4A6gF3AKtE5I/gslFE/ojhubcBnVW1FdAa6CoiRwL9gSmq2gSYEtxGRJoBvYHmQFfgGRGJlDUcAvQDmgSXrsHxi4HfVbUxMAh4MLZfu4gqVbIFSMuXw223JeQlnXNxsnYtXHABnHiirf+YOtXKm1SvHnZkKWOXiURVS6lqRVWtoqpVoy5VYqmzpSbydb1scFFsDcqI4PgIoGdwvQcwWlW3qer3wFKgnYjUBaqq6nRVVeDlPOdEnmsscHyktRJ3HTrA1VfDU09ZUS3nXGpRtXLuTZvCqFFwxx1W1r1Tp7AjSzlxLf4iIqVFZC62KdZkVZ0B7KOqqwGCn3sHD98PWBF1+srg2H7B9bzHdzpHVbOADUDiNjO+/34rxta3r/WtOudSw/LlcMopcN55cOCBNrg+cKAXxttDcU0kqpqtqq2BDKx1UdCuLfm1JLSA4wWds/MTi/QTkVkiMmvNmjWFRL0bKleGF16A776zbzPOueSWlWXrQJo3txInTz4Jn34KhxwSdmQpLSHlKFV1PTAVG9v4JeiuIvgZmQG2EqgfdVoGsCo4npHP8Z3OEZEyQDXgt3xef5iqtlXVtnWKexXhscdaiYTHH4fPPive53bOFZ958+Coo6wcxXHH2UZT11zjO8wVg7glEhGpIyLVg+sVgROAr4G3gQuCh10ATAiuvw30DmZiNcIG1WcG3V8bReTIYPzj/DznRJ7rTOCDYBwlsR580Kp89u0LW7cm/OWdcwXYuhX694c2beDHH2H0aPjvf0tMZd5EiGeLpC7woYjMx6YRT1bVicADQBcRWQJ0CW6jqguBMcAi4D3gKlXNDp7rCuB5bAD+O2BScPwFoJaILAWuJ5gBlnBVqljdnW++sU3ZnXPJ4YMPrNvqwQdtZtbixXDOOT6lt5hJGF/gw9S2bVudNWtWfJ68Xz8bM/nsMyuj4JwLx2+/wY03wksv2WD6sGHQuXPYUaU0EZmtqm3zuy99t+wKw8MPQ7161sW1bVvY0ThX8mRnW9dV06bw8svWpfXVV55E4swTSXGqVg2eew4WLYJ77w07GufSX04OzJ0LgwZZTazateHvf7fxj9mzbYp+xYphR5n2dnurXVeIrl3hoousT/b0022AzzlXPHJyrIXx4Ye2Av3jj+H3oHJT48Zw5plwwgn202djJYyPkcTD+vU2T71WLZg1C8qVi+/rOZeuIolj6lS7fPRRbuI48ECbfn/ssbYavX79XT+PK7KCxki8RRIP1avb1rynngr33Qf33BN2RM6lhpwc28I2OnH8FiwNO+AA6NUrN3H49N2k4YkkXrp3hz594D//sX/8rVuHHZFzyScnxxYGRieOdevsvkaNbNwj0urwxJG0PJHE0+OPW2nqiy6CmTOtsqhzJVlOjk1GiU4ca9fafQ0bWis+kjgaNAgtTLd7PJHEU82atgFWr17wwANw551hR+RcYmVl2S6Dn3xiA+OffJLbVbX//tCtW27iaNgwxEBdUXgiibeePaF3b6ss2rOnF4dz6S0zE774wpLGxx/b4tzI5m+NG9tnoGNHuzRqFGqorvj4rK1EWLsWmjWzb2Cffw5lPH+7NLFpE0yfnps4ZszIXYzbokVu0ujQwRbrupTls7bCVrs2PPMMnHWWrX6/9dawI3Juz/z+u5VfjySO2bNtNXnp0nDYYbbZW8eOcMwxNv3dlQjeIkmks86Ct9+2PuNmzcKJwcXftm2wciX89FPut3OR3Ev07YLuK+yxpUpB+fK2Tim/n2XKFL044c8/545vfPyxrelQtdc44ojcFsdRR1nxUpe2vEWSLAYPthW5ffvaZjq+8jb1bN9uCWLlSlixYudL5Nivvxb+PIkgYv/h7yrRFJSEcnJsrGPJEnuuypWtlXHWWZY42rXz3QTdDp5IEmnvveHpp60W0KBBVp3UJY/sbFi9etcJYsUK+4aetxVfrRpkZNjK6sMOs58ZGXapVMkeH7lA/tcLu53ffdnZ8Oef1uop6s9Nm3a+nZ0NrVrBZZdZ4mjd2qevu13yRJJo55wDr79uW/OeeiocfHDYEZUsmZm2b8yiRbY3xbff2mZHK1ZYEsnO3vnxlStbYqhf3waPI0kicqx+fe/ScSWeJ5JEE7G1Jc2aWRfXxx97F1c8bNwIX3+dmzAWLbLL999btw3Y+96okc2mO/74/JNEtWq+CZJzhfBEEoZ994UnnoDzz4ennoLrrgs7otT122+5iSI6YaxYkfuYsmWt5demDZx3niXxZs2gSRMbD3DOFYknkrCcdx6MGQO33WZ1uRo3Djui5KVqA9iRJBGdMH75JfdxFSvahkYdO+Ymi2bNrNifr91xLm780xUWERg61MrNX3yxzeYqVYL2GduyxYrzrVtnCzbz+xm5vmxZbulwgKpVLUGcckpusmja1GozlaT30Lkk4YkkTPvtZ7O3+va1BYtXXx12RLtH1abDZmZaF1NhCSH659atu37eatVsMVutWjbTrU0bS7hNm1rSqFfPxy2cSyKeSMJ24YU2i6t/fxsgLl1650Vn+S1E29Wxgo5nZ9u0zsglM7Pg27E+pqAFrSJQo4YlhNq1bSC7devc25FkEbleu7YVuvRpps6lFE8kYROxfd6POsrGSxL5uhUq5C5MK1/+r7crVrRNugp6TORYzZo7J4RatSyJ+Iw059KeJ5JkUL8+LF+eW04jeuFZfovRdnWsoMdGymlELsVRPsM55/BEkjzKlPGZRc65lORTXJxzzhWJJxLnnHNF4onEOedckXgicc45VySeSJxzzhWJJxLnnHNFUuK22hWRNcAPYcexC7WBtWEHsYdSNfZUjRs89rCU1NgbqGqd/O4ocYkkmYnIrF3tiZzsUjX2VI0bPPaweOx/5V1bzjnnisQTiXPOuSLxRJJchoUdQBGkauypGjd47GHx2PPwMRLnnHNF4i0S55xzReKJxDnnXJF4InHOOVcknkicc84ViSeSBBBzs4g0DTuW3ZWqsadq3OCxh8Vj33OeSOJMROoD7wAPAOeLSOWQQ4pZqsaeqnGDxx4Wj71oPJHEXyngIaAZ0AZIpdIKqRp7qsYNHntYPPYi8E3Ci5mI7A30A1YAY1T1BxFZo6pbROQj7BvDAlVdF26kf5Wqsadq3OCxh8VjL17eIilGInIV8CFQGTgPGCgie6vqluAhjwH1gZNEpHTUeZLwYPNI1dhTNe4gBo89BB578fNEUkxEpCxQC7hcVW8FLgOOAMoF95dR1a3Ay8BpQCUROUFE9tKQywukauypGncQm8ceAo89TlTVL3t4AVoDvYEDgtv7Bj/LBj/fATrlc95HwE/AF8CBHnv6x+2xe+zpHHvC35h0uQBXAMuAJ4A5wLF57q8CzAUyoo6VAc4HNgBXeuwlI26P3WNP99hDeYPS4QKMBLoE1y8FJgONgtuCzZx4M+rx+wQ/WwKVoo6X8tjTO26P3WNP99h9jCRGItJaROoF1wX4FaghIqKqzwE/AmeJSFm1v96+wBIR6SAi04FuwWPnq82uKAWgqjkee3rF7bF77CUt9oRm2FS8AI2At4HPgXHA1cHx+4EbgKrB7TZYf2TN4PZzwJ/ARKCXx57+cXvsHntJiz1y8RbJLohIKRGpg/VPfqaqR2JNzWZiK0nHAp2Bv4lIeVWdDawGLgieYjZwr6p2V9VxHnv6xu2xe+wlLfa/CDOLJfMFm0LXAWgKVAmO1QBmAAcHt/+NrSiN9GM+STAgRrBpWHA90WMJKRl7qsbtsXvsJS32vBdvkUQRkUpRN/8PuFhVF6vqRhEpB/wBrAEqBo95FPgGuEVEPgUOB74G0OCvG1yPfx+lxR9ZgJRSsadq3EHskc9QSsUe9Z6nYuz+OQ0p9l0KM4slywUbtBoLTIk6tj9wH7B/1LGDsSZo3vObAR1Cir0ucDpQLepYg2SPHfvmVSu4HtnyOZXe895A/VR6z4PXrgcMjH79FHrf/XMaQuyxXEp8i0RE7gTeBbYDv4pIg+Cu8sBe2DeEiBbAlOC8W0WkH4CqLlLVT4LjCXtPReQ64D2gK/CUiJwZ3FWJJI5dRAYAq4CzRaSSBp8UUuM9vwf793I8MExETg/uqpAisU8EsiKvH0iF990/pyHEHqukCyiRROSf2De044Ebsfo1GwFUdQmwH9At6pS6QG8R+QTr13wz73NqYpvHrYBLVLUf8ArwoohUVdXFQAZJGLuINME+/E8BjYNYIq+f7O95N6A6cJKqXgpMA/YOYviGJH3PAUSkI9ADeExV78kTwxKsPlOyxn4tqfs5LUMKfk53W9hNokRf2HmhTsU89y0CTo+63RkYSm7XyzvAPKBN1GMkXrHuKnZsMVJjYCpQJ3IfNqPjgeD2ickSO9AQ2C+4Xh6oCZQO4ruOYDpjkr7nDQlWDgPlo44fDizFVh8fkWzveT7ve83gvb4NOBYYFsTeLbj/pCSM/S/ve3A72T+nDYF6wfVGWLmSpP+cFuVSYsrIi0g14AXgIGzlJ0BmcF9ZVd0OvBrcH/Eb8ANQGxsA66OqvwXnCOw84JWo2IPXXCoiPwP/FpG3gE7AGOBCERmElVZYEWbsIlIDq0baHPhdRF4E/hsVxxtYOYevRGSqqmYD60mO9zy/2CcC20SkMbbSeASwBXhNRDoA3wIrkzD2l7CxhWnAf7AS5IOwAd2XRORQbEA3GWN/EfvPdZuIlFPVP0nez+lf3ndVHS0iPwD3icibJOHntDiUpK6ta7F+eRGRW4JjkZWf24Pbe2HfGCJ+w755/hnc3oA9QSkNxD1qk1/sYN8wZwN9gAqqei0wHiudsAVbwBRK7MEMlBeAP1S1HTAK6Ih1QwCgqlOwb2edsBYK2OrdI8KKO8bYv8MqsA5U1RHYN85rsb7uw5Iw9k5YF8pX2CK3lqr6hKo+hi2EuwrrKkrG2Dti3T0ESQSS8HO6i9iPE1sncjMwiyT8nBabsJtE8b5gf9zeQI3gdmvgd3K7iUoRzMHGmsjf5Dl/ONA7GWOPelx0l8urQOPg+ogwYg/iPhdr4kdW4ZYFFgCHBrdLBz/3xrpZnsC+FUdaX2G+5wXGns85/wHOCPM9jyH2NpHbec55EOiZ5LHn/TeTjJ/TAt/34FhSfU6L9T0IO4A4/WGjF+p0AZ7Nc/8Y4OXgepmo4w2wmSGto46VS+LYI6WkywOHYkXdxpFbajphsecT9zByE3SF4OcEoEWe8+pi3yhnE/R7k6dPPBljx74VHxb8PT4kdwFZ2P9eYok9MgD8RqrFHhxLxs9pLO976J/TeF3SqmtLRCqJSH+smRuxGvhBRGpHTZvrB/QQkWaqmhU0SwEUeFhV50ZO1tzmdDLGvj2YvVUaOAX4P1Xtpao/Jyr2AuJejg3woqqZQRO/FvB9nqc4BXvP26jqW8GxsN/z5RQee13gdmCaqh6nNmsrGf69xBL7flh3yyepFHtkzADIIfk+p8sp/H0vS0if03hLm8F2EbkCuAQrfDY/6q5NWJfQdlXNEatZs15E7gaeFZHngLoi8qiq/oj10adU7FgJhf9o5KuS9a8mYipvoXFHHesAzFTVzSJyGbAZeFVVX4h6vlKqmhP5PZI49k2q+qqI9NZgfC1R73kxxn6BqmalWOybReQ1VV2BDVAnVDG97/dH3utEvu/xlhYtEhE5BptF87yqXqWqmZH7VHU5Nq5waXAokv3nAcdg/zDGRz5UiVZMsauqatQsj0Qkkd2JG6AOcLSIvAf0BD7NmzAS+J9ZUWP/LHjs9kS+58Uce1YKxj5NbWZfwhXj+56T6Pc9EVK2RSK2sO1C7D/VKVh/5FYRaYeVIlgCLFXVj7BFQF1FpKKqbg3+UVyPDXCNCZ5PEvFNOJ6xxzv+PYy7sqpuxhJfA+AOVf1fPONMZOwJaj157P5vJqmlZItERJpjBc+2AlcCF2Hf1k8A3sKm1DUAHhaR1thgbja5Tc/pqnpqSEkkJWMvQtzbgqe4UVXrRj5UUf3dceexe+wlKfYwSComRxHpCzRT1RtF5CBsAKsRMBN4T1XXBY+7G5tvfhf2jeIfqvpDpG8ykQkk1WMvhrgl6H4rnejuCY/dYy9JsYchJVsk2GrQI8RWpH+L7RqWhU2/Wxf1uErAPFXdBnwMnAy5fZMhNTFTNfaixq3BzzA+VB67x767Ujn2hEvVRLIYW+xzdnB7Plb3aD8RKSsiB4rIUGzmxOLgMf9W1aGJD/UvUjX2VI0bPPaweOwlRKomkjVYE7OziOynqhuxPssWQFVslfRqVT1aVecAqOqW0KLdWarGnqpxg8ceFo+9hEjJRBJ070zE/tgPB4dLA5uDZmdvDUplS5LV7k/V2FM1bvDYw+KxlxwpOdgeISLlsaqmpbCy6n9X1S+D+xI+kL47UjX2VI0bPPaweOzpL6UTCez4Q9dR1ZVhx7K7UjX2VI0bPPaweOzpLeUTSTRJ4ZIDqRp7qsYNHntYPPb0k1aJxDnnXOKV+EEi55xzReOJxDnnXJF4InHOOVcknkicc84ViScS55xzReKJxDnnXJH8P3vEuMNWJ4cFAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE of Models\n",
      "Proposed: 127802974.90992846\n",
      "[[30837.148 32307.781 32724.695 33039.703 35491.53  36949.934 37990.27\n",
      "  39492.18  40988.383 42451.09  43964.33  45710.516 47681.    48571.383]]\n",
      "[31077. 26731. 26722. 26803. 27021. 27489. 27812. 27804. 27731. 28426.\n",
      " 29147. 30979. 32160. 33348.]\n",
      "Estimated beta in SIR model is 0.00\n",
      "Estimated gamma in SIR model is 0.00\n",
      "Estimated theta in SIR model is 0.00\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEOCAYAAAC0BAELAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABLa0lEQVR4nO2debzV4/bH36vTKEol3ZQGkqG503ipjBVFAxIicSVDXMMl0zVfswhXyJCxkOJSiAqRcqKJ0MxRSCmp9GtYvz/Wc2p3nHPa1dn7u/c56/16fV/7+32+3+f5fs53n73Xfp5nPWuJquI4juM4hU2JqAU4juM4RRM3MI7jOE5CcAPjOI7jJAQ3MI7jOE5CcAPjOI7jJAQ3MI7jOE5CcAPjOIWAiAwVkRsLsb3FInJsYbXnOFHgBsZxABFREamXq+xmEXkhnvqqOkBVbwv1jhSR7EToDO3XFJFRIvKriKwWkdkick44Vyf8LSV3oj03Zk5CiPuf0HGclOF5YCZQG9gANAL+Fqkix8kD78E4Thzk9EpE5EoR+UVElolIv5jzz4rI7SJSHhgH7Ccif4RtPxEpISKDRGSBiKwQkVdEpHJM/bNEZEk4d/0O5LQEnlXVtaq6SVW/VNVx4dxH4XVVuHdbETlQRCaEtn8VkRdFZO9w3+eBWsD/wvVXh/I2IvKpiKwSkZkicmSM1nNEZKGIrBGRRSJy5m49XKfI4gbGceLnb0BFoAZwHvCoiFSKvUBV1wLHA0tVdc+wLQUuBboDHYD9gN+ARwFE5DDgMeCscK4KULMAHZ+Fe/cWkVq5zrUPr3uHe08BBLgztH0osD9wc9B7FvA9cGK4/h4RqQG8DdwOVAauAkaJSNVgQIcAx6vqXsDfgRlxPDunGOIGxnHiZyNwq6puVNWxwB/AwXHWvQC4XlWzVXUD9gV/SpgrOQV4S1U/CuduBLYU0NapwMfhukUiMkNEWuZ3sarOV9XxqrpBVZcDD2CGLj/6AGNVdayqblHV8UAWcEI4vwVoKCLlVHWZqn4VzwNwih9uYBzH2AyUylVWCjMqOaxQ1U0xx+uAPeNsvzYwOgw5rQLmhntWw3oWP+RcGHpBK/JrSFV/U9VBqtog1J8BjBERyet6EdlXREaIyI8i8jvwArDPDrSemqM16D0CqB60nQYMAJaJyNsickh8j8ApbriBcRzje6BOrrK6wJJdaCuvEOU/YMNKe8dsZVX1R2AZNmwFgIjsgQ2T7fhGqr8C92FGqnI+974zlDdW1QpYDyXWGOWu8wPwfC6t5VX1rnDPd1X1OKA68A3wZDxaneKHGxjHMUYCNwQX4BLBbfdE4LVdaOtnoIqIVIwpGwrcISK1AcJ8Rrdw7jWgq4gcISKlgVsp4LMpIneLSEMRKSkiewEXAvNVdQWwHBvCOiCmyl7YcN6qML/yrzz0xl7/AnCiiHQSkQwRKRucHGqKSDUROSnMxWwI7W7eiWfjFCPcwDiOcSvwKTAZm4C/BzhTVefsbEOq+g3wMrAwDDHtBzwEvAm8JyJrsIn61uH6r4CLgZew3sxvQEHraPYARgOrgIXYkNZJoa11wB3AJ+HebYBbgObAamzy/vVc7d2JGddVInKVqv4AdAOuwwzWD5hRKhG2K4GlwEpsLueinX1GTvFAPOGY4ziOkwi8B+M4juMkBDcwjuM4TkJwA+M4juMkBDcwjuM4TkLwYJeBffbZR+vUqRO1DMdxnLRi+vTpv6pq1bzOuYEJ1KlTh6ysrKhlOI7jpBUiku9iZB8icxzHcRKCGxjHcRwnIbiBcRzHcRKCz8E4jpOybNy4kezsbP7888+opRR7ypYtS82aNSlVKnfQ8fxxA+M4TsqSnZ3NXnvtRZ06dcgnG4GTBFSVFStWkJ2dTd26deOu50NkjuOkLH/++SdVqlRx4xIxIkKVKlV2uifpBmY3UYV77oHFi6NW4jhFEzcuqcGuvA9uYHaT+fPhxhuhXj3o2xe++SZqRY7jOKmBG5jd5KCDYOFCGDgQXn0VDjsMTjkFvvgiamWO4xQGP//8M2eccQYHHHAAmZmZtG3bltGjRydVw+LFi2nYsGGe5S+99NIutfnggw+ybt26rcd77hlv9u/4cQNTCNSoAYMHw5IlcN118P77kJkJxx8PkydHrc5xnF1FVenevTvt27dn4cKFTJ8+nREjRpCd/dd8cJs2bUq6voIMzI705DYwicANTCFStSrcfrsZmv/8B6ZPh3btoH17ePddm69xHCd9mDBhAqVLl2bAgAFby2rXrs3AgQMBePbZZzn11FM58cQT6dixIytXrqR79+40btyYNm3aMGvWLABuvvlm7rvvvq1tNGzYkMWLF7N48WIOPfRQzj//fBo0aEDHjh1Zv349ANOnT6dJkya0bduWRx99NE99gwYN4uOPP6Zp06YMHjz4L3omTZpE165dt15/ySWX8OyzzzJkyBCWLl3KUUcdxVFHHbX1/PXXX0+TJk1o06YNP//8824/P3dTTgAVK8K118Jll8GwYXDvvdC5s/VqrrsOuneHEm7aHWfn+Oc/YcaMwm2zaVN48MF8T3/11Vc0b968wCamTJnCrFmzqFy5MgMHDqRZs2aMGTOGCRMmcPbZZzNjB5rnzZvHyy+/zJNPPkmvXr0YNWoUffr0oV+/fjz88MN06NCBf/3rX3nWveuuu7jvvvt46623ADN4sXomTZqUZ71LL72UBx54gIkTJ7LPPvsAsHbtWtq0acMdd9zB1VdfzZNPPskNN9xQoPYdkdCvORFZLCKzRWSGiGSFssoiMl5E5oXXSjHXXysi80XkWxHpFFOeGdqZLyJDJLgziEgZERkZyqeKSJ2YOn3DPeaJSN9E/p35sccecOmlsGCBGZrVq+Hkk6FhQ3j+eYigR+04zm5w8cUX06RJE1q2bLm17LjjjqNy5coATJ48mbPOOguAo48+mhUrVrB69eoC26xbty5NmzYFIDMzk8WLF7N69WpWrVpFhw4dALa2GQ+xenaG0qVLb+3t5OjYXZLRgzlKVX+NOR4EfKCqd4nIoHB8jYgcBvQGGgD7Ae+LSH1V3Qw8BvQHPgPGAp2BccB5wG+qWk9EegN3A6eJSGXgJqAFoMB0EXlTVX9Lwt/7F0qXhvPOMy+z116z4bOzz4abboJrrrHysmWjUOY4aUQBPY1E0aBBA0aNGrX1+NFHH+XXX3+lRYsWW8vKly+/dV/zGAcXEUqWLMmWLVu2lsWuJylTpszW/YyMDNavX4+q7rJ7dqyegu6bm1KlSm29Z0ZGRqHMKUUxUNMNGB72hwPdY8pHqOoGVV0EzAdaiUh1oIKqTlF7957LVSenrdeAY0LvphMwXlVXBqMyHjNKkVKyJPTubb38N9+EffeFAQPggAPggQfgjz+iVug4TixHH300f/75J4899tjWsoImxtu3b8+LL74IwKRJk9hnn32oUKECderU4YvgWvrFF1+waNGiAu+79957U7FiRSYHL6GcNnOz1157sWbNmnzbqV27Nl9//TUbNmxg9erVfPDBB3HXLQwSbWAUeE9EpotI/1BWTVWXAYTXfUN5DeCHmLrZoaxG2M9dvl0dVd0ErAaqFNDWdohIfxHJEpGs5cuX7/IfubOUKAEnnghTppjH2aGHwpVXQp06cNtt8Fsk/SzHcXIjIowZM4YPP/yQunXr0qpVK/r27cvdd9+d5/U333wzWVlZNG7cmEGDBjF8uP3+Pfnkk1m5ciVNmzblscceo379+ju89zPPPMPFF19M27ZtKVeuXJ7XNG7cmJIlS9KkSRMGDx78l/P7778/vXr1onHjxpx55pk0a9Zs67n+/ftz/PHHbzfJX+ioasI2YL/wui8wE2gPrMp1zW/h9VGgT0z5U8DJQEvg/ZjydsD/wv5XQM2YcwswA/Mv4IaY8huBKwvSmpmZqbvEpk2q/fqpTpy4a/UDU6aonniiKqjutZfqNdeo/vTTbjXpOGnP119/HbUEJ4a83g8gS/P5Xk1oD0ZVl4bXX4DRQCvg5zDsRXj9JVyeDewfU70msDSU18yjfLs6IlISqAisLKCtwmfxYvNBPuoo2z76aJeaadPGhs1mzoQuXczzrE4dW8D5/feFqthxHCcpJMzAiEh5EdkrZx/oCMwB3gRyvLr6Am+E/TeB3sEzrC5wEDBNbRhtjYi0CfMrZ+eqk9PWKcCEYFHfBTqKSKXgpdYxlBU+Bx5obmIPPWRxYjp0gGOPhU8+2aXmGjeGl1+2ps48Ex5/3OZoTj3VFm36WhrHcdKFRPZgqgGTRWQmMA14W1XfAe4CjhORecBx4RhV/Qp4BfgaeAe4WM2DDOBCYBg28b8A8yADG0arIiLzgSswjzRUdSVwG/B52G4NZYmhbFnzR1640GbrZ8+GI46ATp3gs892qcmDDjLX5gULbH7mgw9s0WbLlubivGFDIf8NjuM4hYyo/yQGoEWLFpqVlVU4ja1dC489ZmGWly+3VZa33AKtWu1Wk88/D0OGwNy5UK0aXHiheaFVq1Y4sh0n1Zg7dy6HHnpo1DKcQF7vh4hMV9UWeV3v68kTQfnycNVV1qO56y74/HNo3Rq6drX4MbvY5IAB8NVXNuWTmQk33wy1atk6Gg+u6ThOquEGJpHsuaetpFy0yFZXfvoptGgB3brBl1/uUpMi0LEjvP02fPst9O8Po0aZwWnf3vY9QoDjOKmAG5hksNdeFpxs8WJb6PLRR9C8OfToYW5ju0j9+vDww/Djjzb1k51tqQIOPNC80Hw9jePsPhkZGTRt2pSGDRty6qmnJjwC8a6SX0j/KHEDk0wqVIAbbjBDc/PNMHGiBds75RSYM2eXm61YES6/HObNgzFjzMBcfTXUrGnzNHPnFpJ+xymGlCtXjhkzZjBnzhxKly7N0KFDtzu/efPmfGo6bmCioGJFC0S2aJGlw3zvPfNPPu00+PrrXW42I8NG3yZMsI5R797wzDOWBK1TJxg3DmLCEjmOs5O0a9eO+fPnM2nSJI466ijOOOMMGjVqxJ9//km/fv1o1KgRzZo1Y+LEiYBFN+7WrRudO3fm4IMP5pZbbtna1gMPPEDDhg1p2LAhD4Y4a2vXrqVLly40adKEhg0bMnLkSMBC93fo0IHMzEw6derEsmXLtpbvKKR/lHi4/iipVAluvdXCkD/wgK2lefVVswz//jcccsguN924MTz1lPkYPPEE/Pe/cMIJNqx26aXmGJCABHaOkzAiiNa/HZs2bWLcuHF07mxhDadNm8acOXOoW7cu999/PwCzZ8/mm2++oWPHjnz33XfbXbfHHnvQsmVLunTpgojwzDPPMHXqVFSV1q1b06FDBxYuXMh+++3H22+/DcDq1avZuHEjAwcO5I033qBq1aqMHDmS66+/nqeffjqukP5R4j2YVKByZctUtmiROQW8+SY0aAB9+kD4J91VqlaF66+3pl96yWzaJZfY8NmVV1q54zj5s379epo2bUqLFi2oVasW5513HgCtWrWibt26wPZh+g855BBq16691cAcd9xxVKlShXLlytGzZ08mT57M5MmT6dGjB+XLl2fPPfekZ8+efPzxxzRq1Ij333+fa665ho8//piKFSvy7bffMmfOHI477jiaNm3K7bffTnZ29m6F9E8W3oNJJfbZB+68E664wmbpH33UlvX36wd33w1Vquxy06VLw+mn2/bZZ9ZZGjLEfr2dfroZIV9u4KQyEUTrB7bNweRmR2H6c8gddl9E8r2+fv36TJ8+nbFjx3LttdfSsWNHevToQYMGDZgyZcp2165atWqXQ/onC+/BpCJVq9oizYULbTxr+HD79n/55UKJFdOmjTW1eLH1YsaMsQ7T6afbOhvHcXaO2DD93333Hd9//z0HH3wwAOPHj2flypWsX7+eMWPGcPjhh9O+fXvGjBnDunXrWLt2LaNHj6Zdu3YsXbqUPfbYgz59+nDVVVfxxRdfcPDBB7N8+fKtBmbjxo189dVXcYf0jxI3MKlMtWoweLAtzqxTB844wyJhLllSKM3XqGF2bPFiGDQI3nrLsm2eeiqEVOKO48TBRRddxObNm2nUqBGnnXYazz777NZEYkcccQRnnXUWTZs25eSTT6ZFixY0b96cc845h1atWtG6dWv+8Y9/0KxZM2bPnk2rVq1o2rQpd9xxBzfccAOlS5fmtdde45prrqFJkyY0bdqUTz/9FIgvpH+UeKiYQKGGikkEmzfDI4/YWBbYnM3AgeY6VkisWGHDEEOGwO+/2zKdG2+EmBQSjpNU0j1UzLPPPktWVhaPPPJI1FIKBQ8VU1TJyIDLLrMxrA4dbOFL27a7tVAzN1Wq2DrQxYvNi3rCBFsP2q0bpLLtdRwnNXEDk27Urm1jWSNG2FBZZqZFCVi/vtBuUamSrQNdssS8qD/+2KI4d+kCU6cW2m0cp8hzzjnnFJney67gBiYdEbFFmXPn2oKWu+6CRo0spn8hUrGiDZEtXmyh1KZONQeBzp0trJrjJAMfxk8NduV9SLiBEZEMEflSRN4KxyNFZEbYFovIjFBeR0TWx5wbGtNGpojMFpH5IjIkJB4jJCcbGcqnikidmDp9RWRe2PpSFKlc2VZTTphgRufYY82lecWKQr1NhQrWSVq0yLylp0+Hww+32338caHeynG2o2zZsqxYscKNTMSoKitWrKBs2bI7VS/hk/wicgXQAqigql1znbsfWK2qtwbj8Jaq/iVam4hMAy4DPgPGAkNUdZyIXAQ0VtUBItIb6KGqp4lIZSAr3FeB6UCmquYb/jHlJ/l3xPr1NoFy7702xvXQQxYRIAF+8mvXwtCh5oH2yy9w5JE2Z9OhQ0Ju5xRjNm7cSHZ2Nn/++WfUUoo9ZcuWpWbNmpQqVWq78oIm+VHVhG1ATeAD4GjMeMSeE+AH4KBwXAeYk0cb1YFvYo5PBx4P++8CbcN+SeDX0O7Wa8K5x4HTC9KamZmpRYKZM1VbtlQF1eOPV128OGG3WrtWdfBg1b/9zW7Xrp3q+PGqW7Yk7JaO46QYQJbm8726U0NkIlJCRCrsRJUHgauBvEIstgN+VtV5MWV1w3DahyLSLpTVALJjrskOZTnnfgBQ1U3AaqBKbHkedYo2jRvDlCnmb/zRR7aC8sEHzc25kNljD4sPtXChuTYvWADHHWfZot99t1DWhDqOk8bs0MCIyEsiUkFEygNfA9+KyA6jqolIV+AXVc0vhePpwMsxx8uAWqraDLgCeCkYs7wGXXK+uvI7V1CdWI39RSRLRLKWL1+ej8w0JAkuzbGUK2dLchYssOg2339vjgBt2sAbbyTEtjmOkwbE04M5TFV/B7pj8x+1gHiiqh0OnCQii4ERwNEi8gKAiJQEegIjcy5W1Q2quiLsTwcWAPWx3kfNmHZrAkvDfjawf0ybFYGVseV51NmKqj6hqi1UtUXVqlXj+JPSjCS4NMdStixcdBHMnw+PPw4//wzdu1sE54cessWbjuMUH+IxMKVEpBRmYN5Q1Y3k0RvIjapeq6o1VbUO0BuYoKp9wuljsXmVrUNfIlJVRDLC/gHAQcBCVV0GrBGRNsF77GzgjVDtTSDHQ+yUcA/F5mY6ikglEakEdAxlxY+8XJobNzbPswRRpoylcp43D0aOhL/9zYbSata0ztTChQm7teM4KUQ8BuZxYDFQHvhIRGoDu/tbtDfbD48BtAdmichM4DVggKquDOcuBIYB87GezbhQ/hRQRUTmY8NqgwBCvduAz8N2a0xbxZNYl2aAY46Bc8+FlYl7LKVKQa9e8MkntobmxBMt2k29etazmTTJ52kcpyizS27KIlIyTKoXGdLeTXlniHVprlzZMpJ165aUWy9dasnPhg615TpNmljvpndvG2JzHCe92K1YZCJSTUSeEpFx4fgwtg1LOelIuXK2NH/6dAup3L07nH8+/PFHwm+9334Wp/OHH+DJJ80BoF8/qFXL1tL89FPCJTiOkyTiGSJ7Fpu/2C8cfwf8M0F6nGTSuLFlHxs0yIbPmjWz4yRQrhz84x+WFuD9983j7LbbzND07QtffJEUGY7jJJB4DMw+qvoKYS1LGBpzx9OiQunSlkVz0iTYuNEWsdx8M2xKzgioiE0HvfkmfPstDBgAo0aZw1v79vD66+7m7DjpSjwGZq2IVCF4jolIG2xBo1OUaN/e1smccQbccosZmvnzkyrhoINswWZ2Ntx/vw2jnXyyOQXcfz+sWpVUOY7j7CbxGJgrMHfgA0XkE+A5YGBCVTnRULEiPPec+RZ/9x00bQrDhiXd1WvvveGKK8y+vf66Lee56ipzcx440NyfHcdJfXZoYFT1C6AD8HfgAqCBqnpC3aJMr142OdKmjU3+9+gBEUQ6yMiwW0+aZHMyp5xiDm8HHwxdu8L48bAlryBEjuOkBPF4kZ0KlFPVr7DFliNFpHmihTkRU7MmvPcePPAAjBtn+WbGjo1MTrNm8OyzFobmppvg88+hY0c45BCTWMgZChzHKQTiGSK7UVXXiMgRQCdgOPBYYmU5KUGJErb0PisL9t3XUlpefDGsWxeZpGrVzMB8/z08/7wdX3mleVuffbbF+fTFm46TGsRjYHJ8eLoAj6nqG0DpxElyUo5GjWDaNPsm/+9/oXlzW0MTIWXKQJ8+lvBs1ixzeR4zBv7+d5s6euwxj33mOFETj4H5UUQeB3oBY0WkTJz1nKJE2bJw332WlvmPP2x+5j//SQkf4kaNLATN0qU2R5ORYUE3a9Qwt+cZM6JW6DjFk3gMRS9soWVnVV0FVAZ2GK7fKaIcfTTMnm3+w9dfb+kAFi2KWhUAe+5pPgnTp1uH69RTzSmuWTOzh8OHJyyQtOM4eRCPF9k6VX0dWC0itYBSwDcJV+akLpUqwcsvwwsvmLFp0sS+vVNk8kMEWraEp5+GH3/clirgnHMsVM3ll8M3/h/sOAknHi+yk0RkHrAI+DC8jiu4llPkEYEzz7QJkGbN7Nu7V6+Uc+eqVAkuvdRyr02aZInQHn0UDj3UOmOvvAL/939Rq3Scokk8Q2S3AW2A71S1LpbL5ZOEqnLSh9q1LQXA3Xdb+spGjcy9OcUQsdG8l1+2CAF33mkje6edBvvvD9ddB4sXR63ScYoW8RiYjSHTZAkRKaGqE4GmiZXlpBUZGXD11Zb0Ze+9oVMni8GfohMe1apZfM8FC2yJT5s2Zh8POMA8sf/3v5TwXXCctCceA7NKRPYEPgJeFJGHgLgjIYpIhoh8KSJvheObReRHEZkRthNirr1WROaLyLci0immPFNEZodzQ0JmS0SkjIiMDOVTRaROTJ2+IjIvbJ5eIBk0a2Yz7AMH2sRHy5Y2hJailChhQ2ZvvGG9lxtvhC+/hJNOgrp1zQhNm5YyU0uOk3bEY2C6AeuAy4F3sIySJ+7EPS4D5uYqG6yqTcM2FrbmmekNNAA6A//NSaGMLezsj6VRPiicBzgP+E1V6wGDgbtDW5WBm4DWQCvgppA62Uk05cpZxMp33rFsma1bmwNAirP//hbjc8kSi+Z82GEWYLN1a0shcOmlNoeTpCDTjlMkyNfAiEg9ETlcVdeq6hZV3aSqw4EZwN7xNC4iNbEFmsPiuLwbMEJVN6jqIiw9cisRqQ5UUNUpauk3n8NC1uTUyfn2eg04JvRuOgHjVXWlqv4GjGebUXKSQadOtgClbVtzALjgAtiwIWpVO6RUKejZ0+zjzz+bbczMtORoRx0F1avbos6xY9Piz3GcSCmoB/MgsCaP8nXhXDw8CFxNyCUTwyUiMktEno7pWdQAfoi5JjuU1Qj7ucu3qxPy1KwGqhTQ1naISH8RyRKRrOURBHMs8uy7r034X3ONrYA84gjrIqQJlStb+JkxY+DXX+G11yz+2auv2lxN1aqW3eC115KSDNRx0o6CDEydvKImq2oWUGdHDYtIV+AXVc0dU+Qx4EDMUWAZcH9OlTya0QLKd7XOtgLVJ1S1haq2qFq1ah5VnN2mZEm46y4YPdpSADRvDu++G7WqnaZ8eVtb+uKL8Msv1oM57TSL6HzqqWZsune3hZ0rV0at1nFSg4IMTNkCzpWLo+3DgZNEZDEwAjhaRF5Q1Z9VdbOqbgGexOZIwHoZ+8fUrwksDeU18yjfro6IlAQqAisLaMuJiu7dLWhmjRpw/PFw661pG2u/TBn7E558EpYts7mZ/v3Nv6FvX/NS69jR4qEtWxa1WseJjoIMzOcicn7uQhE5D9hhpENVvVZVa6pqHWzyfoKq9glzKjn0AOaE/TeB3sEzrC42mT9NVZcBa0SkTZhfORt4I6ZOjofYKeEeioW26SgilcIQXMdQ5kTJQQfBZ59ZlMqbboITT0z7n/slS9r6mocesgjP06ZZcrQlS7bFQzv8cHMYWLgwarWOk1xE8/HBFJFqwGjg/9hmUFpgkZR7qOpPcd9E5EjgKlXtKiLPY8NjCiwGLghGBBG5HjgXc4P+p6qOC+UtgGexntM4YKCqqoiUBZ4HmmE9l96qujDUORe4Lki4Q1WfKUhjixYtNCsrK94/ydkdVGHoULjsMvsGHjXKhs6KEKowd65l5Hz9dXN/Bov03LOnJU879NBIJTpOoSAi01W1RZ7n8jMwMZWPAhqGw69UdUIh60sJ3MBEwNSp9k27fLnFbznvvKgVJYxFi2wa6vXX4dNPzQC1aGFOBKefDvvsE7VCx9k1dsvAFBfcwETE8uXmivX++2ZgHnnEUgMUYZYtg5EjzQV6xgxzje7SxeZvTjgBSnu2JSeNKMjAeF4XJ1qqVrVFJ9dfD089ZRMWKRL+P1FUr26RdL78EmbOtEWcU6ZAjx4W7XngQPOH8N9+TrrjBsaJnowMuP12CwK2cKGtbBw7NmpVSaFxY8vjlp0Nb78Nxxxj3mktW0LDhnDPPZZywHHSkbgMjIjUFpFjw345EdkrsbKcYknXrubrW7u27d90U7GJOlmypA2PjRxpQ2iPP25xQ6+5xkLVdOoEL70E69ZFrdRx4ieefDDnY2FYHg9FNYExCdTkFGcOOMBmwc85x9bKdOmScjlmEk2lSrau5pNPbG3qdddZgrQzz4S//c2mqj76KG2XETnFiHh6MBdjiyZ/B1DVecC+iRTlFHPKlbP5mCeegIkTzYX588+jVhUJBx0Et91m01ITJ1o0gZEjbe1NvXpw882WdsBxUpF4DMwGVd2a8y+smPfpRyexiMD559vPeBGLY/bEE8V25rtECTjySHjmGQvC+dxzcOCB1smrVw/atYNhw2D16qiVOs424jEwH4rIdUA5ETkOeBX4X2JlOU6gRQublznqKIvI3K9fyiYySxbly8NZZ1kctCVL4D//MW/v88+3IbTTT7ccN3/+GbVSp7gTz0LLEljelY5YEMl3gWFaxBbQ+DqYFGfzZhsruvVWc70aNcp+wjuAdeymTbOezYgRFoFnzz0tGs8pp1jstHLxRBB0nJ1kd1fy9wDGqmqRzn7hBiZNGDfOZru3bLFv05NOilpRyrFxI0yYYGkERo82H4ny5c1f4tRTzdiULx+1SqeosLsLLU8CvhOR50WkS5iDcZxoOP54+OILm3jo1s22b76JWlVKUaqUuTU/+ST89JMNpfXpY04Cp55qaXpOPdWcBTyPjZNIdmhgVLUfUA+bezkDWCAi8WSodJzEUKcOTJ5skw8TJ9qKxAsvtNlvZztKloRjj7XYokuXWs+mb1/4+GPo3dsCKfTsaWtsfv89arVOUSPuWGQiUgpLO9wPaKeqRSpDlw+RpSnLl9u8zNChFsPs6qvhiit8DGgHbN5sDnqvvmrTWcuWWZ6bTp2sd3PiiVCxYtQqnXRgt4bIRKSziDwLzMdyrgwDqhdYyXGSRdWq8PDD8NVXluXr3/+G+vVtHU0xiQKwK2RkQPv29uiys61HM2CAOeyddZYNo3XtagE5f/starVOuhLPHMw52Mr9+qraV1XHquqmeG8gIhki8qWIvBWO7xWRb0RkloiMFpG9Q3kdEVkvIjPCNjSmjUwRmS0i80VkSEg8RkhONjKUTxWROjF1+orIvLD1xSna1K9vP8UnT7bYKv/4hyVfGTeu2K6diZcSJWyZ0YMPWtK0Tz+Fiy+GWbMsoEK1ahbG5umn0z4/nJNkEh6uX0SuwBKVVQgJxzpimSc3icjdAKp6TTAOb6lqwzzamAZcBnwGjAWGqOo4EbkIaKyqA0SkN5YI7TQRqQxkhfsqljAtU1Xz/S3mQ2RFCFUzNoMG2TL3Y46Be++FZs2iVpZWqFoAhVdfNY+0xYttTqd5czjsMEuYlrPVrWu9Iqf4sUtDZCIyObyuEZHfY7Y1IhLXdKCI1AS6YMNqAKjqezE9oM+w2GYFtVEdM05Twtqb54Du4XQ3YHjYfw04JvRuOgHjVXVlMCrjsfkjpzggYos/vv7achnPmGERms8+236iO3EhAq1amW1euNBSCFx1lU1vvfOOBeI86SQLZ1O+PDRpYo4Dt9wCr7wCs2fDhiK9uMHZEfm6HKvqEeF1dyInPwhcDeTXxrnAyJjjuiLyJRb37AZV/RioAWTHXJMdygivPwSdm0RkNVAltjyPOlsRkf5Af4BatWrtzN/lpAOlS1uylbPPhrvusjGgV16Byy+33o3PYseNiNnozMxtZb/9Zh7ic+fa9vXXttjzlVe2jUqWKGHxS3P3eA45BCpUiOZvcZLHDte0iMjzqnrWjsryqNcV+EVVp4vIkXmcvx7YBLwYipYBtVR1hYhkAmNEpAEWPSA3OeN6+Z0rqM62AtUngCfAhsgK+nucNGbvvc3AXHQR3HCD7Q8bZg4BF1zgKSR3kUqVoG1b22JZt86iQMcanrlzbTps48Zt19WoYcYm1vg0bAhVqiT373ASRzyLJhvEHoSFlpn5XBvL4cBJInICUBaoICIvqGqfMOneFTgmJ+RMiBSwIexPF5EFQH2s9xE7jFYTWBr2s4H9geygqyKwMpQfmavOpDg0O0WZWrVs9f/ll8O//mW9myFD4M47LUyx5PW7xNlZ9tjD/CuaNt2+fONGG2rLbXieegrWrrVrSpSwYbeBAy38nL8laY6q5rkB1wJrsF7G72FbA6wA7syvXj5tHYlN4IPNhXwNVM11TVUgI+wfAPwIVA7HnwNtsJ7JOOCEUH4xMDTs9wZeCfuVgUVApbAtymkrvy0zM1OdYsSWLapjx6o2bKgKqm3bqn7ySdSqiiWbN6suWaL6zjuq11yjus8+9pYceqjqo4+qrlkTtUKnIIAsze+7P78Tuu2Lf6eMST5txBqY+dj8yIyw5RiIk4GvgJnAF8CJMfVbAHOABcAjbPN+K4tFGJgPTAMOiKlzbiifD/TbkUY3MMWUTZtUhw1TrV7dPg4nn6z63XdRqyrWrF+v+uyzqpmZ9pZUqKB66aWq334btTInLwoyMPEGu5ygqqvD8d7Akao6psCKaYa7KRdz1q6FBx6Au+8216cLL4Qrr7T0zU4kqMLUqfDII+Y4sHGjRRq45BILSedu0anB7kZTnqGqTXOVfamqRWpRgRsYB7B4ZjffbJEiN2+2GezTTrP4KfvtF7W6YstPP9lbkhNT7YADzGejXz+oXDlqdcWb3Y2mnNc1HlHZKZpUqwaPPQbz55u32fr18M9/Qs2allLyscfgl1+iVlns+Nvf4MYbbbHnyJHmgXbVVfa2nH8+zJwZtUInL+LpwTwNrAIexVx9BwKVVPWcRItLJt6DcfLl22/tW23ECHN7ysiAo4+2nk3Pnuav6ySdmTNt+OzFF+13QLt2NnzWo4elLHCSw+4OkZUHbgSOZVtGy9tVdV1hC40SNzDODlGFOXO2GZsFC+ybrGNHW8J+0km+ejACVq6EZ56BRx+FRYtsJHPAAOjf3zqkTmLZLQOTR2O1gNNU9d7CEJcquIFxdgpVS3w2YoQZnB9+sHj3XbpYz6ZrV1sQ4iSNzZttMecjj8C775rt79XLejWtW/uamkSxu3MwiMg+InKhiHwETAT8d4FTvMmJnXLvvTYx8MknFhXg00/NwOy7L5x+OowZA3/+GbXaYkFGhtn1d96xUc0LL4Q33zQ/jVatLPWAvxXJJd8ejIjsBfTAsljWB0ZjPZcCg1OmK96DcQqFzZstucqIERaCeMUKGzbr3t0Mz3HH+QRBElmzBl54wXo1X39tHmdnn23DZ4ceGrW6osEuDZGJyHps8eINwGRVVRFZqKoHJE5qdLiBcQqdjRstR/HIkfD667B6tX3D9ewJffpYxi8ft0kKqpZd+4kn7K3YuNGcAi64wKIElS0btcL0ZVeHyK7DVso/BlwrIgcmQpzjFFlKlbKVgU8/betr3nzTVgiOGGEuz02aWNDN9eujVlrkETHHvxEjLIPn3Xfbepo+fczl+fLLzUHQKVzyNTCqOlhVWwMnYd5jY4D9ROQaEamfJH2OUzQoU8YS3b/wghmbp56yb73zz7fFHNdea44CTsLZd1+4+mqL+Pz++3DsseaBdthh1ql88UWfqyksdjjJr6oLVfUOVW0EtMQiFo9LuDLHKarssQece64lQps0CTp0gHvusbSQvXqZw4CneU44JUpYstORI/Pu1VxxhfdqdpeEp0xOF3wOxomUxYvtZ/SwYbBqleUlvuwycwwoUyZqdcWGLVu2zdWMHm1zNe3bm1OAz9XkzW67KTuOk2Dq1DGX5+xsC0ezfj307Ws5bG66CZYti1phsSCvXs2PP3qvZldxA+M4qUT58rYM/auv4L33bAHHbbdZVOc+feDzz6NWWGzIPVdzzDHw8MM2V9Ohg8/VxEO+BkZEPgivd+/ODUQkQ0S+FJG3wnFlERkvIvPCa6WYa68Vkfki8q2IdIopzxSR2eHcEBHz7RSRMiIyMpRPFZE6MXX6hnvMCxk0HSd9ELE1M//7n33DXXSReaG1amUrB0eM2D7/sJMwcno1r7xivZq77rLX2F7NN99ErTI1KagHU11EOmBpj5uJSPPYbSfucRkQ26kcBHygqgcBH4RjROQwLCtlAyzr5X9FJCfjw2NAf+CgsHUO5ecBv6lqPWAwcHdoqzJwE9AaaAXcFGvIHCetqFcPHnzQxmqGDLHFm6efbsNqd9wBy5dHrbDYUK0aXHMNzJsH48dv69Uceigcfjjcf78F4nYC+WUiA07BvMXWYOFhYrcJ+dXL1UZNzIgczbaMlt8C1cN+deBb3Zai+dqYuu8CbcM138SUnw48HntN2C8J/Iq5VG+9Jpx7HDi9IK2e0dJJGzZvVn37bdWOHS3lY5kyqv36qX75ZdTKiiU//aR6552qjRvb2wGqDRqoXned6tSp9nYVZSggo2VB62BeU9XjgXtU9ahc29Fx2q8HgauBLTFl1VR1WbjHMmDfUF4DS6WcQ3YoqxH2c5dvV0dVNwGrgSoFtOU46U+JEnDCCRbR8euvzeV55Eho1swmB0aNgk2bolZZbKhWDQYNsvQBCxfC4MFQtao5CLRuDfvvb3HR3n3XkqUWJ+JZB3ObiJwkIveFrWs8DYfrflHV6XFqyStmhhZQvqt1YjX2F5EsEcla7sMMTjpy6KHw3//a8Nl998H338Mpp9iw2v33m8uzkzTq1rX8dBMn2nra4cOhTRt4/nno3NkMz2mnwUsvFY+3ZocGRkTuxOZRvg7bZaFsRxyOzd8sBkYAR4vIC8DPIlI9tF0dyEkPmA3sH1O/JrA0lNfMo3y7OiJSElsEurKAtrZDVZ9Q1Raq2qJq1apx/EmOk6LsvTdceaVNALz+us3P5KR8HDjQJg2cpFKligXWHDUKfv0V3nrLjMukSXDmmWZsjjvOAnEW2SAO+Y2d6bb5i1lAiZjjDGDWjurlauNIts3B3AsMCvuDsCE4sMn9mUAZoC6wEMgI5z4H2mA9k3HACaH8YmBo2O8NvBL2KwOLgEphWwRULkijz8E4RY4vvlDt21e1VClVEdWuXVXff191y5aolRVrNm9W/fRT1auvVj344G3zNs2bq95yi+rMmen1FlHAHEy8BqZyzHHl3TQwVbCJ/3nhNbbt64EFmCPA8THlLYA54dwjbItAUBZ4FZiPRX4+IKbOuaF8PtBvRxrdwDhFlmXLVP/9b9WqVe0j37Ch6rBhquvWRa3MUdW5c1Xvuku1bVv7HQCqdeqoXnaZ6sSJqhs3Rq2wYAoyMPGkTD4duAvzHhOgPebtNaLAimmGh4pxijx//gkvv2wuz7NmwT772OzzhRdC9epRq3OAn36ypU9vvGGLOzdssAwPXbva+tu2baNW+Fd2O2VymCtpiRmYqar6U+FKjB43ME6xQdUmAh580L7NSpaE3r1tdrr5zixxcxLJH3+Y59kbb9jbtGqVOQkOGmRZIFIlldBuG5jigBsYp1gyf74t3nz6aVi71rJw/fOf0K2b5SB2UoK1a+HJJ80xMDsbmjY1Q3PKKdG/TR7s0nGcvKlXzwxMdrZ9e33/vYUNrlcPHnjAsnA6kVO+vNn9BQvgmWdstLN3bzj4YHj88dSNieYGxnEcc3O+4grr0YwaZasDr7zS3JwvvdTjn6QIpUvDOedYLNTXX7f5mQEDbP3NPffA779HrXB7CjQwIlJCROYkS4zjOBFTsiT07AkffQRZWdCjBwwdCvXr27DZxImeDC0FKFHC3pqpU+GDD6BRI4uRVqsWXH89/PLLjttIBgUaGFXdAswUkVpJ0uM4TqqQmQnPPQdLlti31qefWmL7Zs1sKbqHo4kcEXtL3nvPMjkcdxzceadld7jkEstjFyXxDJFVB74SkQ9E5M2cLdHCHMdJEapXt5w0339vM80bNthS9Jz5m7Vro1boAC1awKuvWuqAM8+0rJz16llagdmzo9EUzzqYDnmVq+qHCVEUEe5F5jhxsmWLxT255x745BObCLj4YvvJvO++O67vJIUff7TAm0OH2m+Arl3N8+zwwwv3PrvlRRYMyWKgVNj/HPiiUBU6jpM+lCgBJ50Ekyfb1q7dtqybF11krk5O5NSosS3+6a23wpQpcMQR0L49jB2bnKm0eIJdng+8huVUAQt7PyaBmhzHSRcOPxzGjLG0AWeeCcOGmUPAaaeZk4ATOZUrw4032lTaQw/ZvEyXLraW5uWXEzuVFs8czMVYZOTfAVR1HttyuDiO41jagGHD7NvrX/+Cd96Bli0t5eO777rnWQpQvrx5nC9YAM8+axm3zzjD1tIMHZqYtygeA7NBVf8v5yCExff/Fsdx/sp++1nS+h9+gHvvtRnnzp3t5/KLL9q3mhMppUpB374wZ451PqtWhdGjExN6Jh4D86GIXAeUE5HjsOjF/yt8KY7jFBkqVLB8NIsW2dLzjRvNnalePRun+eOPqBUWe0qUsKVNU6aY91lC7hHHNYOA5cBs4AJgLHBDYuQ4jlOkyFl6PmcOvPmmrQT85z/t9cYbU2dFYDFGxH4PJKTtOKMplwYOwYbGvo0dMisquJuy4ySJTz+14bM33oAyZaBfPwtTU69e1MqcXWC33JRFpAuW6GsIluxrvogcH0e9siIyTURmishXInJLKB8pIjPCtlhEZoTyOiKyPubc0Ji2MkVktojMF5EhIjZaKCJlQnvzRWSqiNSJqdNXROaFre+O9DqOkyT+/ncb9P/6axs2e+opm2nu1cuWoztFh/wykeVswDdAvZjjA4Fv4qgnwJ5hvxQwFWiT65r7gX+H/TrAnHzamga0ZVvK5OND+UVsnzJ5pG7LurkwvFYK+5UK0usZLR0nIpYuVR00SLViRUvneNRRquPGpVfe4GIMBWS0jGcO5hdVjQ2luhDY4cBpuHfOTF6psG0djwu9kF7AywW1E5KdVVDVKeGPeQ7oHk53A4aH/deAY0K7nYDxqrpSVX8DxgOdd6TZcZwIqF7dAmh9/70NnX33HRx//DbPM495lrbka2BEpKeI9MTikI0VkXPCUNP/sNX8O0REMsIQ2C/YF/7UmNPtgJ/V1tXkUFdEvhSRD0WkXSirAWTHXJMdynLO/QCgqpuA1UCV2PI86sTq6y8iWSKStXz58nj+JMdxEkWO59nChZYALdbzzGOepSUF9WBODFtZ4GegA3Ak5lFWKZ7GVXWzqjYFagKtRKRhzOnT2b73sgyoparNgCuAl0SkAjYs9pemw2t+5wqqE6vvCVVtoaotqlatusO/x3GcJFC6tE38z5ljjgA1a8Jll5nn2U03wa+/Rq3QiZOS+Z1Q1X6FdRNVXSUik7BhqjlhsWZPIDPmmg3AhrA/XUQWAPWx3kfNmOZqAkvDfjawP5Ad2qwIrAzlR+aqM6mw/h7HcZJATsyzk06yoJr33GNBte69F8491xKi1a0btUqnAOLxIqsrIg+IyOs7E65fRKqKyN5hvxxwLOYwQM6+qmbnuj4j7B8AHAQsVNVlwBoRaRPmV84G3gjV3gRyPMROASaEeZp3gY4iUklEKgEdQ5njOOnI4Ydbb+arryxX8BNPwEEHWayTL7+MWp2TD/FM8o/Boik/jHl95Ww7ojowUURmYXM241X1rXCuN3+d3G8PzBKRmdiE/QBVXRnOXQgMA+ZjLtPjQvlTQBURmY8Nqw0CCPVuC/f9HLg1pi3HcdKVww6z+ZmFC+Hyyy1tQPPm0KmTpXb0mGcpRTz5YKaqausk6YkMX2jpOGnIqlUWqfHBB+Hnny0L59VXw8knQ0ZG1OqKBbu10BJ4SERuEpG2ItI8ZytkjY7jODvP3ntbFq3Fi23Y7PffLVVA/frw2GOwfn3UCos18RiYRsD5wF1sGx67L5GiHMdxdoqyZeH882HuXBg1CvbZx5Kf1a4Nt98OK32EPAriMTA9gANUtYOqHhW2oxMtzHEcZ6fJyICePeGzz2DSJEtUf+ON5uJ8+eU2jOYkjXgMzExg7wTrcBzHKTxEoEMHyw08cyb06AEPPwzNmll8eicpxGNgqgHfiMi7O+Om7DiOkxI0bgzPPw9ffAHlypnhSVQKR2c78l1oGcNNCVfhOI6TaBo3hqwsOPNMuPBCi9z86KM2f+MkhB0aGFX9MBlCHMdxEk6lSvC//8HNN9vk/+zZ5hSw//5RKyuSxLOSf42I/B62P0Vks4j8ngxxjuM4hU5GBtx2m+Wk+eYbWzszaVLUqookOzQwqrqXqlYIW1ngZCzxmOM4TvrSvTtMmwZVqsCxx8LgwT4vU8jEM8m/Hao6BnA3Zcdx0p9DDoGpUy2g5hVXWHqAdeuiVlVk2OEcTMgJk0MJoAV5hL53HMdJSypUgNdeg7vughtusDQBo0fDAQdErSztiacHc2LM1glYg2WSdBzHKRqUKAHXXWfrZn74wRZovvNO1KrSnni8yAotL4zjOE5K07mzuTL36AEnnGCeZtdeaws3nZ0mXwMjIv8uoJ6q6m0J0OM4jhMtBxwAn35qsc2uv94MzvDhsNdeUStLOwoaIlubxwZwHnDNjhoWkbIiMk1EZorIVyJySyi/WUR+FJEZYTshps61IjJfRL4VkU4x5ZkiMjucGxISjyEiZURkZCifKiJ1Yur0FZF5YeuL4zhOvJQvDy++CA88AG++Ca1awbffRq0q7cjXwKjq/Tkb8ARQDugHjADimf3aABytqk2ApkBnEWkTzg1W1aZhGwsgIodhicgaYKmV/5uT4RJ4DOiPZbk8KJwHM3a/qWo9YDBwd2irMhaBoDXQCrgpZLZ0HMeJDxELkDl+PKxYAS1bWlZNJ24KnOQXkcoicjswCxtOa66q16jqLztqWI0/wmGpsBXkfdYNGKGqG1R1EZa9spWIVAcqqOqUkA75OaB7TJ3hYf814JjQu+mEZdBcqaq/AePZZpQcx3Hi56ijYPp0OPhgWztz442weXPUqtKCfA2MiNyLpRteAzRS1ZvDl3XciEiGiMwAfsG+8KeGU5eIyCwReTqmZ1ED+CGmenYoqxH2c5dvV0dVNwGrgSoFtJVbX38RyRKRrOXLl+/Mn+Y4TnFi//3h44/h3HNt4v/EE+G3nfo6LJYU1IO5EtgPuAFYGhMuZk28oWJUdbOqNgVqYr2Rhthw14HYsNkyLIEZQF5uGlpA+a7WidX3hKq2UNUWVatWLeAvcRyn2FO2LAwbZpky33/fhsxmz45aVUpT0BxMCVUtlytUTIWc4525iaquAiYBnVX152B4tgBPYnMkYL2M2IhzNYGlobxmHuXb1RGRkkBFYGUBbTmO4+w6IjBggMUuW7cO2rSBkSOjVpWy7HSomHgRkaoisnfYLwcci+WVqR5zWQ9gTth/E+gdPMPqYpP501R1GbBGRNqE+ZWzgTdi6uR4iJ0CTAjzNO8CHUWkUhiC6xjKHMdxdp+//93mZZo1g9694aqrYNOmqFWlHPHkg9lVqgPDgydYCeAVVX1LRJ4XkabYkNVi4AIAVf1KRF4BvgY2AReras5M2oXAs5gn27iwATwFPC8i87GeS+/Q1koRuQ2bQwK4VVU9KbfjOIVH9eowYYLFMLv/fjM4w4dbemYHAFGPHgpAixYtNCsrK2oZjuOkI8OHw8UXW8iZwYPNGaCYrP4Xkemq2iKvcwkbInMcxyk29O1rE/6ZmfCPf1iYmezsHdcr4riBcRzHKQzq1oUPPoCHH4aPPoKGDeG554p1jhk3MI7jOIVFiRJwySUwcyY0amQ9m27dYNmyqJVFghsYx3GcwqZePXNlfuABCzXToAG89FKx6824gXEcx0kEGRkWy2zGDAszc+aZcMop8MsOI20VGdzAOI7jJJKDD4bJk+Huu+Gtt6w38+qrUatKCm5gHMdxEk1GBlx9NXz5JdSpA7162QLNX3+NWllCcQPjOI6TLA47DKZMsYCZr79uvZkxY6JWlTDcwDiO4ySTkiW3Zcrcbz9Lz9ynD6wsesFG3MA4juNEQePGMG0a3HSTBcxs2NDmaIoQbmAcx3GiolQpuPlmmDoVqlSxPDP9+sGqVVErKxTcwDiO40RN8+Y2ZHbddbb6v1EjeDf9A8C7gXEcx0kFypSBO+4wJ4C99oLOnaF/f/g9rvyOKYkbGMdxnFSiVSv44gv4178sg2ajRhbjLA1JZMKxsiIyTURmishXInJLKL9XRL4RkVkiMjomKVkdEVkvIjPCNjSmrUwRmS0i80VkSEg8RkhONjKUTxWROjF1+orIvLD1xXEcJ10oWxbuuccWaJYpA8ceC506wRtvpFVis0T2YDYAR6tqE6Ap0FlE2gDjgYaq2hj4Drg2ps4CVW0atgEx5Y8B/bEslwcBnUP5ecBvqloPGAzcDSAilYGbgNZYSuabQmZLx3Gc9OHvf7dQM7feCnPmQPfuFrX5ttvSIoBmwgyMGn+Ew1JhU1V9T1VzTPBnQM2C2gkpliuo6pSQDvk5oHs43Q0YHvZfA44JvZtOwHhVXamqv2FGrTOO4zjpxh57wI03wpIltjjzkEPg3/+2zJm9esHEiSkbRDOhczAikiEiM4BfsC/8qbkuOZdt6Y8B6orIlyLyoYi0C2U1gNjMPdmhLOfcDwDBaK0GqsSW51HHcRwn/ShZ0hZljh8P334Ll14K778PRx9tEQGGDEk59+aEGhhV3ayqTbFeSisRaZhzTkSuBzYBL4aiZUAtVW0GXAG8JCIVgLzyjuaY6/zOFVRnKyLSX0SyRCRr+fLlcf5VjuM4EVO/Ptx/P/z4IzzzjHmdXXYZ1KhhGTW/+CJqhUCSvMhUdRUwiTBMFSbduwJnhmEvVHWDqq4I+9OBBUB9rPcRO4xWE1ga9rOB/UObJYGKwMrY8jzqxOp6QlVbqGqLqlWrFsaf6jiOkzzKlYNzzrGFmllZcPrplncmMxNat4bhw2H9+sjkJdKLrGqMh1g54FjgGxHpDFwDnKSq63JdnxH2D8Am8xeq6jJgjYi0CfMrZwNvhGpvAjkeYqcAE4LBehfoKCKVwuR+x1DmOI5TNMnMNLfmpUvhwQdh9WozPjVqwJVXwrx5SZeUyB5MdWCiiMwCPsfmYN4CHgH2AsbnckduD8wSkZnYhP0AVc2J/nYhMAyYj/VscuZtngKqiMh8bFhtEECod1u47+fArTFtOY7jFF323tuGy+bOhQkTzMV5yBAbVuvYEUaPTpqrs2iKeh8kmxYtWmhWVlbUMhzHcQqfZcusd/PEE5Cdbb2a/v1tvma//XaraRGZrqot8jrnK/kdx3GKOtWrm6vzokXWg2nQwKI416plaZw/+CAhrs4lC71Fx3EcJzUpWdIWa3bvbnMyjz9uXmjff2+pAwoZ78E4juMURw46CO67z4bMXn45IbdwA+M4jlOcKVcODjwwIU27gXEcx3ESghsYx3EcJyG4gXEcx3ESghsYx3EcJyG4gXEcx3ESghsYx3EcJyF4qJiAiCwHlkStIx/2AX6NWsQu4tqjIV21p6tuKL7aa6tqnuHo3cCkASKSlV+sn1THtUdDumpPV93g2vPCh8gcx3GchOAGxnEcx0kIbmDSgyeiFrAbuPZoSFft6aobXPtf8DkYx3EcJyF4D8ZxHMdJCG5gHMdxnITgBsZxHMdJCG5gHMdxnITgBiZCxLhaRA6NWsvO4tqjIV21p6tucO27gxuYiBCR/YG3gbuAs0WkfMSS4sa1R0O6ak9X3eDadxc3MNFRArgHOAzIBNIpxIRrj4Z01Z6uusG17xYlk33D4oqI7Av0B34AXlHVJSKyXFXXiciH2C+MOaq6Ilqlf8W1R0O6ak9X3eDaCxvvwSQBEbkYmAiUB/oAt4nIvqq6LlzyALA/0ElEMmLqSdLF5sK1R0O6ak9X3UGDay9k3MAkGBEpBVQBBqjqtcAFQGugdDhfUlXXA88BJwF7iMixIrKnRhxmwbVHQ7pqT1fdQZtrTwSq6lshb0BToDdwQDj+W3gtFV7fBjrkUe9D4Efgc+BA1+7aU117uup27cnRnvQHU9Q34EJgIfAQ8CVwZK7zewEzgJoxZSWBs4HVwEWu3bWng/Z01e3ak6c9kgdUlDfgeeC4sH8+MB6oG44F8+QYFXN9tfDaGNgjpryEa3ftqaw9XXW79uRp9zmY3UREmorIfmFfgF+ASiIiqvok8D1wqoiUUntX/wbME5F2IjIF6BKunaXm7VECQFW3uHbXnkra01W3a49Oe1Itb1HagLrAm8BnwGjgklB+J3AlUCEcZ2LjnZXD8ZPA/wFvAT1cu2tPde3pqtu1R6c9Z/MezE4iIiVEpCo2/vmpqrbBuqyHia2cfQ04GjhERMqo6nRgGdA3NDEduFVVu6rqaNfu2lNVe7rqdu3Raf8LUVq3dNwwV792wKHAXqGsEjAVODgc346toM0ZJx1CmIgjJHkL+8keM3ftrr3I63bt0WnPvXkPJg5EZI+Yw/eA81R1rqquEZHSwO/AcqBcuOZ+4FvgGhH5BGgJfAOg4V0P+4kfAzX9OQurXHtyted8vtJGe8zzTivdQbt/TiPSni9RWrdU37DJsteAD2LKagF3ALViyg7GurK56x8GtItIe3WgJ1Axpqx2mmivBFQJ+zlpvdPpufcG9k+n5w7sB9wWe+80eub+OY1Aezyb92DyQURuBMYCG4FfRKR2OFUG2BP7RZFDQ+CDUO9aEekPoKpfq+rHoTxpz1pE/gm8A3QGHhaRU8KpPUh97TcDS4FeIrKHhk8R6fHcb8H+Z44BnhCRnuFU2VTWHnS/BWzKuXcgHZ65f04j0B4vKScoFRCRS7FfdMcAV2HxfdYAqOo8oAbQJaZKdaC3iHyMjZuOyt2mJreb3QT4h6r2B14AnhaRCqo6F6hJ6mo/CPtieBioF/TkaEj1594F2BvopKrnA5OBfYOGb0nR5y4i7YFuwAOqekuu+8/D4lelnG4AEbmM9P2cliRNP6c7RdRdqFTZ2H4BUrlc574GesYcHw0MZdvwzdvATCAz5hpJlNb8tGOLrOoBk4CqOecwD5O7wnHHFNNeB6gR9ssAlYGMoPGfBNfLFH3udQirpYEyMeUtgfnYiuvWqfbccz3zyuE5XwccCTwRdHcJ5zuliu6Cnnk4TvXPaR1gv7BfFwvbkhaf013din24fhGpCDwF1MdWugL8Gc6VUtWNwIvhfA4rgSXAPtjE21mqujLUEdh+oi1Z2sM954vIT8DtIvI60AF4BThHRAZjISZ+SAHtlbAIrw2A30TkaeB/MVpexUJbzBaRSaq6GVhFajz3vLS/BWwQkXrY6urhwDrgZRFpB3wHZEepPQ/dz2BzF5OB/2Ch3gdjE8nPiEgzbCI5Ut35aH8a+9LdICKlVfX/SN3P6V+eu6qOEJElwB0iMooU/ZzuLj5EBpdhY/4iIteEspyVrhvD8Z7YL4wcVmK/Uv8vHK/GGiihgYSrNvLSDvaLdDpwFlBWVS8DxmAhJNZhC7Mi0x68Yp4CflfVVsBLQHtsSAMAVf0A+0XXAevRgK1Ybp3i2hdgUW1vU9Xh2K/Uy7Dx9OZRac9HdwdsKGY2tnivsao+pKoPYAv8LsaGnCLTXYD29tiwEcG4QAp+TvPRfpTYOpergSxS9HNaKETdhYpqw9703kClcNwU+I1tw00lCD7kWFf721z1nwV6p6L2mOtih21eBOqF/eERaz8DGy7IWXlcCpgDNAvHGeF1X2zI5iHsl3ROjy1ltedR5z/AyVE+9x3ozsw5zlXnbqB7mv2/pOLntMDnHspS7nNaaM8gagFJfsNjFyAdBzye6/wrwHNhv2RMeW3MU6VpTFnpFNaeE7K7DNAMC4Y3mm0hvaPW/gTbjHfZ8PoG0DBXverYr9DphLF1co27p6J27Jd08/CeTGTb4rikPfdd1J0z8fxqVLp38/8lFT+n8Tz3lPicJmIrFkNkIrKHiAzCuss5LAOWiMg+Me59/YFuInKYqm4K3VsABe5V1Rk5lXVbtzwVtW8M3mQZwAnAe6raQ1V/ShHti7HJZVT1zzBcUAVYlKuJE7Dnnqmqr4eydNBeHbgemKyqR6l5kSXlue+m7hrYsM3Hyda9O9pz5iSALaTe53QxO37upYjwc5pIivwkv4hcCPwDCxg3K+bUH9jQ0kZV3SIW02eViNwEPC4iTwLVReR+Vf0eG/9PK+1YKIn/aM5PKxu/TZYb5g61x5S1A6ap6loRuQBYC7yoqk/FtFdCVbfk/C0prP0PVX1RRHprmMNL1nMvJN19VXVTMnUXgva1IvKyqv6ATYwnlUJ67nfmPOtkPvdEU6R7MCJyOObRM0xVL1bVP3POqepibN7i/FCU82thJnA49g8zJufDlmwKSbuqqsZ4nSTry2JntANUBf4uIu8A3YFPchuSNNL+abh2YzKfeyHq3pSG/y+T1bwMk04hPvctyX7uyaDI9WDEFuudg33ZfoCNd64XkVZYSIZ5wHxV/RBb3NRZRMqp6vrwz3IFNrH2SmhPkvGrOZHak/Srf1e0l1fVtZhRrA3coKrvJ1prsrQn+rmnq+5Eak8G6fzck02R6sGISAMsUNx64CKgH/br/ljgdcz1rzZwr4g0xSaQN7OtCztFVU+MyLgUR+0bQhNXqWr1nA9czJi6ay9iul17dNqjQIqS0RSRc4HDVPUqEamPTZzVBaYB76jqinDdTZi//L+xXyBnquqSnLHPZH45u3ZdkqNZRDKSPdSRrtrTVbdrj057FBSpHgy2+rW12Ar877Asb5swN8EVMdftAcxU1Q3AR8DxsG3sM6KuanHWruE1ig9cumpPV93g2qPSnnSKmoGZiy1i6hWOZ2ExoWqISCkROVBEhmKeHHPDNber6tDkS/0Lrj0a0lV7uuoG115sKGoGZjnWVT1aRGqo6hpsTLQhUAFbEb5MVf+uql8CqOq6yNRuj2uPhnTVnq66wbUXG4qUgQnDRG9h/wT3huIMYG3ovvbWEJJcUix3gmuPhnTVnq66wbUXJ4rUJH8OIlIGixJbAgtff7qqfhHOJX0SfGdw7dGQrtrTVTe49uJAkTQwsPUfoKqqZketZWdx7dGQrtrTVTe49qJOkTUwsUgah15w7dGQrtrTVTe49qJIsTAwjuM4TvIp9pNQjuM4TmJwA+M4juMkBDcwjuM4TkJwA+M4juMkBDcwjuM4TkJwA+M4juMkhP8Hp3CnII64C9gAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE of Models\n",
      "Proposed: 1638765842.9084125\n",
      "[[508427.6  496424.88 485543.9  473169.25 462228.3  452013.62 443993.22\n",
      "  434710.6  426443.03 422412.6  411802.12 401770.03 394003.12 382215.94]]\n",
      "[497888. 485219. 469264. 451661. 432205. 418394. 409830. 399007. 388030.\n",
      " 376471. 364864. 344987. 328111. 317800.]\n"
     ]
    }
   ],
   "source": [
    "#Pred with STAN\n",
    "file_name = './save/stan'\n",
    "checkpoint = torch.load(file_name)\n",
    "model.load_state_dict(checkpoint['state'])\n",
    "optimizer.load_state_dict(checkpoint['optimizer'])\n",
    "model.eval()\n",
    "\n",
    "prev_x = torch.cat((train_x, val_x), dim=1)\n",
    "prev_I = torch.cat((train_I, val_I), dim=1)\n",
    "prev_R = torch.cat((train_R, val_R), dim=1)\n",
    "prev_cI = torch.cat((train_cI, val_cI), dim=1)\n",
    "prev_cR = torch.cat((train_cR, val_cR), dim=1)\n",
    "prev_Vt = torch.cat((train_Vt, val_Vt), dim=1)\n",
    "prev_edges = torch.cat((train_edges, val_edges), dim=0)\n",
    "\n",
    "train_locs.reverse()\n",
    "\n",
    "for loc_name in train_locs:\n",
    "    cur_loc = loc_list.index(loc_name)\n",
    "    prev_active_pred, _, prev_phy_active_pred, _, h = model(prev_x, prev_cI[cur_loc], prev_cR[cur_loc], N[cur_loc], prev_I[cur_loc], prev_R[cur_loc], V=prev_Vt[cur_loc], e_weights=prev_edges)\n",
    "\n",
    "    test_pred_active, test_pred_recovered, test_pred_phy_active, test_pred_phy_recover, _ = model(test_x, test_cI[cur_loc], test_cR[cur_loc], N[cur_loc], test_I[cur_loc], test_R[cur_loc], h, V=test_Vt[cur_loc], e_weights=test_edges)\n",
    "\n",
    "    print('Estimated beta in SIR model is %.2f'%model.alpha_scaled)\n",
    "    print('Estimated gamma in SIR model is %.2f'%model.beta_scaled)\n",
    "    print('Estimated theta in SIR model is %.2f'%model.theta_scaled)\n",
    "\n",
    "    #Cumulate predicted dI\n",
    "    pred_I = []\n",
    "    sir_I = []\n",
    "\n",
    "    for i in range(test_pred_active.size(1)):\n",
    "        # below is regular prediction\n",
    "        cur_pred = (test_pred_active[0, i, :].detach().cpu().numpy() * dI_std[cur_loc].reshape(1, 1).detach().cpu().numpy()) + dI_mean[cur_loc].reshape(1, 1).detach().cpu().numpy()\n",
    "        # below is SIR model prediction\n",
    "        sir_pred = test_pred_phy_active[0, i, :].detach().cpu().numpy()\n",
    "        # below is average of the two predictions\n",
    "        #cur_pred = (cur_pred + test_pred_phy_active[0, i, :].detach().cpu().numpy()) / 2\n",
    "        cur_pred = np.cumsum(cur_pred)\n",
    "        cur_pred = cur_pred + test_I[cur_loc, i].detach().cpu().item()\n",
    "        pred_I.append(cur_pred)\n",
    "\n",
    "        sir_pred = np.cumsum(sir_pred)\n",
    "        sir_pred = sir_pred + test_I[cur_loc, i].detach().cpu().item()\n",
    "        sir_I.append(sir_pred)\n",
    "    pred_I = np.array(pred_I)\n",
    "    sir_I = np.array(sir_I)\n",
    "\n",
    "    #I_true = get_real_y(active_cases[:], history_window, pred_window, slide_step)\n",
    "    I_true = active_cases[cur_loc][-test_window:]\n",
    "    I_true = I_true[history_window:history_window+pred_window]\n",
    "\n",
    "\n",
    "    import matplotlib.dates as mdates\n",
    "    dates = [test_start_date + pd.DateOffset(days=i+history_window) for i in range(pred_window)]\n",
    "    plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%m/%d/%Y'))\n",
    "    #plt.plot(I_true[cur_loc, -1, :],c='r', label='Ground truth')\n",
    "    #stan_res = pickle.load(open(f\"stan_results/{loc_name}_stan_results_new\", \"rb\"))\n",
    "    #stan_sirv_res = pickle.load(open(f\"stan_results/{loc_name}_stan_results_sirv\", \"rb\"))\n",
    "    plt.plot(dates, I_true, c='r', label='Ground truth')\n",
    "    #plt.plot(dates, stan_res[-1, :], c='g', label=\"STAN (SIRV)\")\n",
    "    #plt.plot(dates, sir_I[-1, :], c='purple', label=\"SIRVIC\")\n",
    "    plt.plot(dates, pred_I[-1, :],c='b', label='Proposed')\n",
    "    plt.gcf().autofmt_xdate()\n",
    "    plt.legend()\n",
    "    plt.ylabel(\"Number of Active Cases\")\n",
    "    plt.title(loc_name)\n",
    "    plt.show()\n",
    "\n",
    "    prop_mse = ((pred_I[0] - I_true)**2).mean()\n",
    "    sirvic_mse = ((sir_I[0] - I_true)**2).mean()\n",
    "    #stan_mse = ((stan_res[0] - I_true)**2).mean()\n",
    "    print(\"MSE of Models\")\n",
    "    print(f\"Proposed: {prop_mse}\")\n",
    "    print(pred_I)\n",
    "    print(I_true)\n",
    "    #print(f\"STAN: {stan_mse}\")\n",
    "    #print(f\"SIRVIC: {sirvic_mse}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# prop_mse = ((pred_I[0] - I_true)**2).mean()\n",
    "# sirvic_mse = ((sir_I[0] - I_true)**2).mean()\n",
    "# stan_mse = ((stan_res[0] - I_true)**2).mean()\n",
    "# print(\"MSE of Models\")\n",
    "# print(f\"Proposed: {prop_mse}\")\n",
    "# print(f\"STAN: {stan_mse}\")\n",
    "# print(f\"SIRVIC: {sirvic_mse}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# final graph visualization\n",
    "import networkx as nx\n",
    "import seaborn as sns\n",
    "\n",
    "def plot(g, attention, ax, nodes_to_plot=None, nodes_labels=None,\n",
    "         edges_to_plot=None, nodes_pos=None, nodes_colors=None,\n",
    "         edge_colormap=plt.cm.viridis, labels=[], edge_layer_2=None):\n",
    "    \"\"\"\n",
    "    Visualize edge attentions by coloring edges on the graph.\n",
    "    g: nx.DiGraph\n",
    "        Directed networkx graph\n",
    "    attention: list\n",
    "        Attention values corresponding to the order of sorted(g.edges())\n",
    "    ax: matplotlib.axes._subplots.AxesSubplot\n",
    "        ax to be used for plot\n",
    "    nodes_to_plot: list\n",
    "        List of node ids specifying which nodes to plot. Default to\n",
    "        be None. If None, all nodes will be plot.\n",
    "    nodes_labels: list, numpy.array\n",
    "        nodes_labels[i] specifies the label of the ith node, which will\n",
    "        decide the node color on the plot. Default to be None. If None,\n",
    "        all nodes will have the same canonical label. The nodes_labels\n",
    "        should contain labels for all nodes to be plot.\n",
    "    edges_to_plot: list of 2-tuples (i, j)\n",
    "        List of edges represented as (source, destination). Default to\n",
    "        be None. If None, all edges will be plot.\n",
    "    nodes_pos: dictionary mapping int to numpy.array of size 2\n",
    "        Default to be None. Specifies the layout of nodes on the plot.\n",
    "    nodes_colors: list\n",
    "        Specifies node color for each node class. Its length should be\n",
    "        bigger than number of node classes in nodes_labels.\n",
    "    edge_colormap: plt.cm\n",
    "        Specifies the colormap to be used for coloring edges.\n",
    "    \"\"\"\n",
    "    if nodes_to_plot is None:\n",
    "        nodes_to_plot = sorted(g.nodes())\n",
    "    if edges_to_plot is None:\n",
    "        edges_to_plot = sorted(g.edges())\n",
    "    nx.draw_networkx_edges(g, nodes_pos, edgelist=edges_to_plot,\n",
    "                           edge_color=attention, edge_cmap=edge_colormap,\n",
    "                           width=2, alpha=0.5, ax=ax, edge_vmin=0,\n",
    "                           edge_vmax=1)\n",
    "    if edge_layer_2 is not None:\n",
    "        nx.draw_networkx_edges(g, nodes_pos, edgelist=edge_layer_2[0],\n",
    "                           edge_color=edge_layer_2[1], edge_cmap=edge_colormap,\n",
    "                           width=2, alpha=0.5, ax=ax, edge_vmin=0,\n",
    "                           edge_vmax=1)\n",
    "\n",
    "    if nodes_colors is None:\n",
    "        nodes_colors = sns.color_palette(\"deep\", max(nodes_labels) + 1)\n",
    "\n",
    "    nx.draw_networkx_nodes(g, nodes_pos, nodelist=nodes_to_plot, ax=ax, node_size=20,\n",
    "                           node_color=[nodes_colors[nodes_labels[v - 1]] for v in nodes_to_plot],\n",
    "                            alpha=0.9)\n",
    "    loc_labels = {}\n",
    "    for node in g.nodes():\n",
    "        if node in labels:\n",
    "            loc_labels[node] = loc_list[node]\n",
    "    nx.draw_networkx_labels(g, nodes_pos, loc_labels, font_size=16)\n",
    "\n",
    "# nx_g = model.g.cpu().to_networkx(edge_attrs=['e'])\n",
    "# #print(nx_g.nodes())\n",
    "# attention = []\n",
    "# mid_attention = []\n",
    "# mid_edges = []\n",
    "# for e in sorted(nx_g.edges()):\n",
    "#     a = nx_g.get_edge_data(e[0], e[1])[0]['e'].item()\n",
    "#     attention.append(a)\n",
    "#     if a > 0.1 and a < 0.9:\n",
    "#         mid_attention.append(a)\n",
    "#         mid_edges.append(e)\n",
    "# labels = [0 if i not in train_locs else 1 for i in range(len(nx_g.nodes()))]\n",
    "# loc_i_list = [loc_list.index(loc) for loc in train_locs]\n",
    "# pos = nx.spring_layout(nx_g, k=0.25, iterations=20)  # positions for all nodes\n",
    "\n",
    "# fig, ax = plt.subplots(figsize=(16, 12))\n",
    "# plot(nx_g, attention, ax, nodes_pos=pos, nodes_labels=labels, nodes_colors=[\"red\" if i not in train_locs else \"blue\" for i in range(len(nx_g.nodes()))], edge_layer_2=(mid_edges, mid_attention))\n",
    "# ax.set_axis_off()\n",
    "# sm = plt.cm.ScalarMappable(cmap=plt.cm.viridis, norm=plt.Normalize(vmin=0, vmax=1))\n",
    "# sm.set_array([])\n",
    "# plt.colorbar(sm, fraction=0.046, pad=0.01)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "cdfa9ace6613655d28fa2370fe498fbea8375f41a5fc4e643d7f8a581612fdbd"
  },
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit ('tf-gpu': conda)",
   "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.12"
  },
  "metadata": {
   "interpreter": {
    "hash": "cdfa9ace6613655d28fa2370fe498fbea8375f41a5fc4e643d7f8a581612fdbd"
   }
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}