stan_sirvc_world_predict.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "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",
    "from stan_sirvc_layer import run_layer, get_features\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": 73,
   "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'\n",
    "#end_date = '2021-06-30'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "163\n",
      "1\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",
    "countries.remove(\"European Union\")\n",
    "countries.remove(\"Europe\")\n",
    "countries.remove(\"Africa\")\n",
    "countries.remove(\"Asia\")\n",
    "countries.remove(\"North America\")\n",
    "countries.remove(\"Oceania\")\n",
    "countries.remove(\"South America\")\n",
    "countries.remove(\"World\")\n",
    "countries.remove(\"Tajikistan\")\n",
    "mask = raw_data['location'].isin(countries)\n",
    "world = [\"World\"]\n",
    "world_mask = raw_data['location'].isin(world)\n",
    "world_raw_data = raw_data.loc[world_mask]\n",
    "raw_data = raw_data.loc[mask]\n",
    "print(len(raw_data['location'].unique()))\n",
    "print(len(world_raw_data['location'].unique()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "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",
    "flight_counts = pd.read_csv('processed_flights/flight_counts_2021_all.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": 76,
   "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": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<bound method DGLHeteroGraph.number_of_nodes of Graph(num_nodes=163, num_edges=3939,\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": 78,
   "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": 79,
   "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_cases = np.nan_to_num(np.array(new_cases))\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",
    "# total_cases = []\n",
    "# for i,loc in enumerate(confirmed_cases):\n",
    "#     total_loc = [loc[0]]\n",
    "#     for j in range(1, len(loc)):\n",
    "#         total_loc.append(total_loc[-1] + new_cases[i][j])\n",
    "#     total_cases.append(total_loc)\n",
    "# total_cases = np.array(total_cases)\n",
    "\n",
    "# confirmed_cases = total_cases\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": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "done\n"
     ]
    }
   ],
   "source": [
    "# load world data for prediction\n",
    "_, world_features, world_active_cases, world_static_feat, world_norms = get_features(world_raw_data, start_date, end_date, world, edges=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "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": 82,
   "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": 83,
   "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": 84,
   "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": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "world_train_features = world_features[0]\n",
    "world_val_features = world_features[1]\n",
    "world_test_features = world_features[2]\n",
    "\n",
    "# x, I, R, cI, cR, yI, yR, Vt, edges\n",
    "# want everything but outer \n",
    "\n",
    "train_I = world_train_features[1]\n",
    "train_R = world_train_features[2]\n",
    "train_cI = world_train_features[3]\n",
    "train_cR = world_train_features[4]\n",
    "train_yI = world_train_features[5]\n",
    "train_yR = world_train_features[6]\n",
    "train_Vt = world_train_features[7]\n",
    "\n",
    "val_I = world_val_features[1]\n",
    "val_R = world_val_features[2]\n",
    "val_cI = world_val_features[3]\n",
    "val_cR = world_val_features[4]\n",
    "val_yI = world_val_features[5]\n",
    "val_yR = world_val_features[6]\n",
    "val_Vt = world_val_features[7]\n",
    "\n",
    "test_I = world_test_features[1]\n",
    "test_R = world_test_features[2]\n",
    "test_cI = world_test_features[3]\n",
    "test_cR = world_test_features[4]\n",
    "test_yI = world_test_features[5]\n",
    "test_yR = world_test_features[6]\n",
    "test_Vt = world_test_features[7]\n",
    "\n",
    "# norms = [dI_mean, dI_std, dR_mean, dR_std]\n",
    "dI_mean = world_norms[0]\n",
    "dI_std = world_norms[1]\n",
    "dR_mean = world_norms[2]\n",
    "dR_std = world_norms[3]\n",
    "\n",
    "#true_dI_train = world_features[0][5] # 0=train, 5=yI; locations, timesteps, values\n",
    "#true_dI_val = world_features[1][5]\n",
    "\n",
    "# target for train and test\n",
    "#true_dI_prev = np.concatenate((true_dI_train, true_dI_val), axis=1)\n",
    "#true_dI_test = world_features[2][5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(18, 3939)\n",
      "(163, 18, 42)\n",
      "(1, 3939)\n",
      "(163, 1, 42)\n",
      "(1, 3939)\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 (now just world)\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": 87,
   "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": 88,
   "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": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "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(world_static_feat[:, 0], dtype=torch.float32).to(device).unsqueeze(-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "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": 91,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----Save best model-----\n",
      "Epoch 0, Loss 58285764.00, Val loss 2014279.62\n",
      "Epoch 1, Loss 65589160.00, Val loss 2988135.75\n",
      "Epoch 2, Loss 77341848.00, Val loss 2157064.00\n",
      "-----Save best model-----\n",
      "Epoch 3, Loss 55480648.00, Val loss 1263402.25\n",
      "-----Save best model-----\n",
      "Epoch 4, Loss 30606026.00, Val loss 679229.38\n",
      "-----Save best model-----\n",
      "Epoch 5, Loss 20227968.00, Val loss 230861.67\n",
      "-----Save best model-----\n",
      "Epoch 6, Loss 8186369.00, Val loss 99342.07\n",
      "-----Save best model-----\n",
      "Epoch 7, Loss 1765047.12, Val loss 35034.74\n",
      "-----Save best model-----\n",
      "Epoch 8, Loss 1095962.38, Val loss 15693.32\n",
      "-----Save best model-----\n",
      "Epoch 9, Loss 323821.75, Val loss 12071.44\n",
      "-----Save best model-----\n",
      "Epoch 10, Loss 198822.16, Val loss 5456.21\n",
      "-----Save best model-----\n",
      "Epoch 11, Loss 98984.95, Val loss 5567.65\n",
      "-----Save best model-----\n",
      "Epoch 12, Loss 41276.39, Val loss 2606.51\n",
      "-----Save best model-----\n",
      "Epoch 13, Loss 32739.57, Val loss 2246.84\n",
      "-----Save best model-----\n",
      "Epoch 14, Loss 24184.75, Val loss 1284.85\n",
      "-----Save best model-----\n",
      "Epoch 15, Loss 12495.68, Val loss 1338.57\n",
      "-----Save best model-----\n",
      "Epoch 16, Loss 8215.79, Val loss 482.17\n",
      "Epoch 17, Loss 9041.10, Val loss 432.42\n",
      "-----Save best model-----\n",
      "Epoch 18, Loss 5348.86, Val loss 575.39\n",
      "-----Save best model-----\n",
      "Epoch 19, Loss 4843.04, Val loss 322.68\n",
      "-----Save best model-----\n",
      "Epoch 20, Loss 2482.50, Val loss 349.20\n",
      "-----Save best model-----\n",
      "Epoch 21, Loss 2065.90, Val loss 194.99\n",
      "-----Save best model-----\n",
      "Epoch 22, Loss 1616.36, Val loss 213.29\n",
      "-----Save best model-----\n",
      "Epoch 23, Loss 1438.91, Val loss 160.36\n",
      "Epoch 24, Loss 1721.14, Val loss 188.97\n",
      "-----Save best model-----\n",
      "Epoch 25, Loss 1204.32, Val loss 188.06\n",
      "-----Save best model-----\n",
      "Epoch 26, Loss 955.84, Val loss 99.90\n",
      "-----Save best model-----\n",
      "Epoch 27, Loss 862.61, Val loss 113.48\n",
      "-----Save best model-----\n",
      "Epoch 28, Loss 727.11, Val loss 93.14\n",
      "-----Save best model-----\n",
      "Epoch 29, Loss 672.12, Val loss 84.81\n",
      "-----Save best model-----\n",
      "Epoch 30, Loss 529.71, Val loss 53.79\n",
      "Epoch 31, Loss 676.17, Val loss 62.37\n",
      "Epoch 32, Loss 666.25, Val loss 67.96\n",
      "-----Save best model-----\n",
      "Epoch 33, Loss 325.12, Val loss 43.57\n",
      "Epoch 34, Loss 512.59, Val loss 46.89\n",
      "Epoch 35, Loss 311.15, Val loss 66.98\n",
      "-----Save best model-----\n",
      "Epoch 36, Loss 259.60, Val loss 42.53\n",
      "Epoch 37, Loss 273.61, Val loss 54.53\n",
      "Epoch 38, Loss 288.85, Val loss 36.89\n",
      "Epoch 39, Loss 359.40, Val loss 34.23\n",
      "Epoch 40, Loss 325.01, Val loss 35.06\n",
      "-----Save best model-----\n",
      "Epoch 41, Loss 173.24, Val loss 40.55\n",
      "Epoch 42, Loss 195.04, Val loss 35.75\n",
      "Epoch 43, Loss 339.25, Val loss 28.94\n",
      "Epoch 44, Loss 275.48, Val loss 36.87\n",
      "-----Save best model-----\n",
      "Epoch 45, Loss 139.50, Val loss 22.78\n",
      "Epoch 46, Loss 186.51, Val loss 47.85\n",
      "Epoch 47, Loss 162.32, Val loss 29.24\n",
      "Epoch 48, Loss 139.46, Val loss 29.45\n",
      "-----Save best model-----\n",
      "Epoch 49, Loss 131.68, Val loss 24.81\n",
      "Epoch 50, Loss 205.81, Val loss 16.56\n",
      "-----Save best model-----\n",
      "Epoch 51, Loss 118.92, Val loss 33.88\n",
      "Epoch 52, Loss 164.03, Val loss 44.13\n",
      "Epoch 53, Loss 152.57, Val loss 43.81\n",
      "Epoch 54, Loss 144.15, Val loss 32.87\n",
      "-----Save best model-----\n",
      "Epoch 55, Loss 118.82, Val loss 32.50\n",
      "Epoch 56, Loss 171.70, Val loss 37.78\n",
      "Epoch 57, Loss 147.41, Val loss 39.13\n",
      "Epoch 58, Loss 132.07, Val loss 27.80\n",
      "Epoch 59, Loss 180.19, Val loss 26.78\n",
      "Epoch 60, Loss 182.00, Val loss 36.84\n",
      "-----Save best model-----\n",
      "Epoch 61, Loss 125.44, Val loss 23.56\n",
      "Epoch 62, Loss 194.44, Val loss 22.33\n",
      "-----Save best model-----\n",
      "Epoch 63, Loss 118.97, Val loss 27.67\n",
      "Epoch 64, Loss 153.63, Val loss 18.47\n",
      "-----Save best model-----\n",
      "Epoch 65, Loss 97.76, Val loss 31.43\n",
      "Epoch 66, Loss 145.91, Val loss 27.90\n",
      "Epoch 67, Loss 106.64, Val loss 32.22\n",
      "Epoch 68, Loss 132.34, Val loss 26.26\n",
      "Epoch 69, Loss 124.72, Val loss 21.07\n",
      "Epoch 70, Loss 162.64, Val loss 18.77\n",
      "Epoch 71, Loss 185.71, Val loss 26.51\n",
      "Epoch 72, Loss 148.40, Val loss 31.91\n",
      "Epoch 73, Loss 117.13, Val loss 30.83\n",
      "Epoch 74, Loss 120.50, Val loss 23.67\n",
      "Epoch 75, Loss 128.37, Val loss 26.83\n",
      "Epoch 76, Loss 117.90, Val loss 27.78\n",
      "-----Save best model-----\n",
      "Epoch 77, Loss 95.73, Val loss 18.51\n",
      "Epoch 78, Loss 129.59, Val loss 13.17\n",
      "Epoch 79, Loss 91.86, Val loss 23.55\n"
     ]
    }
   ],
   "source": [
    "#Train STAN\n",
    "\n",
    "loc_name = 'World'\n",
    "cur_loc = world.index(loc_name)\n",
    "\n",
    "all_loss = []\n",
    "file_name = './save/stan_' + loc_name\n",
    "min_loss = 1e20\n",
    "\n",
    "for epoch in range(80):\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",
    "    #if val_loss < min_loss: \n",
    "    if (val_loss + 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 = val_loss\n",
    "        min_loss = (val_loss + loss) / 2\n",
    "        print('-----Save best model-----')\n",
    "    \n",
    "    print('Epoch %d, Loss %.2f, Val loss %.2f'%(epoch, all_loss[-1], val_loss.item()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Pred with STAN\n",
    "file_name = './save/stan_' + loc_name\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",
    "\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",
    "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",
    "#_, _, _, _, 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",
    "\n",
    "#test_pred_active, test_pred_recovered, test_pred_phy_active, test_pred_phy_recover, _ = model(val_x, val_cI[cur_loc], val_cR[cur_loc], N[cur_loc], val_I[cur_loc], val_R[cur_loc], h, V=val_Vt[cur_loc], e_weights=val_edges)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimated beta in SIR model is 0.01\n",
      "Estimated gamma in SIR model is 0.16\n",
      "Estimated theta in SIR model is 0.05\n"
     ]
    }
   ],
   "source": [
    "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "#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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[10059281.  9874620.  9704572.  9113609.  8907296.  8699681.  8534607.\n",
      "  8368268.  8161630.  7970445.  7799422.  7582181.  7435430.  7278143.]\n"
     ]
    }
   ],
   "source": [
    "#I_true = get_real_y(active_cases[:], history_window, pred_window, slide_step)\n",
    "I_true = world_active_cases[cur_loc][-test_window:]\n",
    "#I_true = world_active_cases[cur_loc][-test_window-valid_window:]\n",
    "I_true = I_true[history_window:history_window+pred_window]\n",
    "print(I_true)\n",
    "\n",
    "\n",
    "# test_pred_phy_active = (test_pred_phy_active - dI_mean[cur_loc]) / dI_std[cur_loc]\n",
    "# change loss here \n",
    "# loss = criterion(test_pred_active.squeeze(), val_yI[cur_loc].squeeze()) + 0.1*criterion(test_pred_phy_active.squeeze(), val_yI[cur_loc].squeeze())\n",
    "# #print(loss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEOCAYAAACTqoDjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5eUlEQVR4nO3dd3hUZfbA8e9JIPQiEhEBSRTpJVRRBASliYpiWYqgCCKKiu5PFxDrYsG1oCiLYgFdCyoKWLAiIqgrvXcQEXEVQXqX8/vjvUOGMEkmJDN3Jjmf57nPzNwyOZlk5szbRVUxxhhjMkrwOwBjjDGxyRKEMcaYkCxBGGOMCckShDHGmJAsQRhjjAnJEoQxxpiQLEEY4yMReUBEXs/i+AYRuTCaMRkTYAnCmGyIyFARmZph35pM9nWLbnTGRI4lCGOy9w3QQkQSAUTkVKAw0CjDvmreuWERkUIRiNWYPGMJwpjszcElhDTvcStgOrAqw751ACLygYhsE5G1InJD4Em86qSJIvK6iOwErsv4g0Skl4j8JCJbRWRYpH4hY8JhCcKYbKjqQeAHXBLAu50JzMqw7xvgLWATcBpwJfCIiFwQ9HRdgIlAWeCN4J8jIrWBMUAv7/qTgcp5/gsZE6Z8lyBE5BUR+V1EloZx7kgRWehtq0VkexRCNPFpBunJoCUuQczMsG8GcB4wWFX3q+pC4CXcB37A96o6WVWPqOq+DD/jSuAjVf1GVQ8A9wJHIvLbGBOGfJcggPFAx3BOVNU7VDVNVdOAZ4H3IxiXiW/fAOeJyElAsqquAb4DzvX21QVWAttUdVfQdT8BlYIe/5zFzzgt+Liq7gG25lH8xuRYvksQqvoNsC14n4icKSKfisg8EZkpIjVDXNodVz1gTCjfA2WA/sC3AKq6E9js7dvsbeVEpFTQdacDvwQ9zmr65F+BKoEHIlIcV81kjC/yXYLIxFjgVlVtDNwJ/Dv4oIhUBVKBr3yIzcQBrzpoLvB3XNVSwCxv3zeq+jOuVPGoiBQVkfpAXzK0NWRhInCxiJwnIknAPyk471ETg/L9P5+IlATOBd4VkYXAC0DFDKd1Ayaq6l9RDs/ElxnAKbikEDDT2xfo3todSMGVJiYB96vqF+E8uaouAwYCb+JKE3/iGryN8YXkxwWDRCQF19hXV0RKA6tUNWNSCD5/ATBQVb+LVozGGBPr8n0Jwqsn/lFErgIQp0HguIjUAE7C1TEbY4zx5LsEISJv4T7sa4jIJhHpC/QE+orIImAZri96QHdggubHopQxxuRCvqxiMsYYk3v5rgRhjDEmb+SrycLKly+vKSkpfodhjDFxY968eX+oanKoY/kqQaSkpDB37ly/wzDGmLghIj9ldsyqmIwxxoRkCcIYY0xIliCMMcaEZAnCGGNMSJYgjDHGhGQJwhhjTEj5qpurMcYUFPv3w5o1sHIl7NgB/frl/c+wBGGMMTFs2zaXBFasOPb2xx/hiLcgbZky0LcviOTtz7YEYYwxPjtyBDZudB/8GZPBli3p5xUpAjVqQOPG0LMn1KoFNWtC9ep5nxzAEoQxxkRNoFookAACSWDVKti3L/28cuXch3+XLi4B1KzpHletComJ0YvXEoQxxuTQkSOwe7er+9++PfvbP/+EdetctVBgAm0R94Ffqxa0aZOeBGrWhOSQMyNFnyUIY0yBcuCA+3Dftev4LfChnt0H/s6d6fX/mSlSBMqWde0DZctCkybQq1d6iaB6dShePMK/bC5ZgjDGxDxV+O03+OOP4z/UM/uwz+zYoUPZ/zwRKF362A/4008/9nFWt2XKQNGiEXxBosQShDHGd6qwdaurgtmw4fjbDRtc/X1WEhOhVKnjtwoVQu8vWfL4fYEP+FKlIMFGiVmCMMZEx/btoT/8A7e7dx97/kknQWoq1K4NnTtDSkr6h32oD/ciRSLTk6cgswRhjMkTu3enf9sP/vAP3N++/djzS5VyCeCMM+CCC1wCSE11tykp7pu88ZclCODSS+HgQVdETUx0RcvA/cz2Zfc4eF+JEq63QuCfv2xZv39jY3Ju3z746afMSwF//HHs+cWKpX/on3tu+v9/4LZcOfvGH+ssQeAarvbuhb/+ctuRI+n3Qz0O55xAV7ZQypY99o2S8bZkyaj82sYc4+BBN1grsyqg//3v2POTktwXn5QU6No1/Zt/aqrbTjnFEkC8swQBTJ+e98+pmp40du5037wyvuFWroRPPz12gAxA+fKZJ5CqVd03M2My+usv10Mn1Hb4cPr9bdtCJ4Bffjn2i01iouu5k5oKF110/P9ixYrWkJvfWYKIEJH0qqby5d3WuPHx56nC77+HfsMuWgRTprhvdsFOPfXYb2s1a0KdOm6QTaz3qy6oVF0pdds2t/35Z+b3//zT9djJ7MM+4wd+YMuq1BqKCFSu7P6P2rY9PgFUqgSF7BOiQLM/v89EXM+MChXg7LOPP37kiCvahyr2z54NEye6D4vAc6WkuGRRp47r/WGJI/dUXZI+cCB927s3/YM9qw/7wP1t27Luf1+4sOu1U66cq4IsVsz9zQoXTt8KFTr2cXZbZueXLu0SQJUqrprImMxYgohxCQlw2mlua9Hi+OOHDrkh/MuWuW35cnf72WfpH0h+Jo59+1wJKbD99lv6/W3b0uNLSHBb4H64t1kdO3TIfZhn/HDP6ZaxBJeV0qXTP+jLlXOvc7lyx+4Ldb9ECauvN7HHEgS4hoBmzdw7Nc4ULpw+dP+KK9L3ByeOQNIIlTgC/cyDk0dWiePIEffBntmHfsbHu3aFfp4SJdzLnZDgnjPQZhO4DbUvu9vA/WBFimS/lSkT3nnBW/Hix3/Qly3r/h7G5BeiOa24jGFNmjTRuXPn5uyiHTvSW9v69YM77nAtwfnU4cOwdu3xiWPVqtCJo0yZYz/wt2xxjaEZJSS4CcZOOeXYrUKF4/edcopLEJESSBiB0oQxJnMiMk9Vm4Q8FqkEISKvABcDv6tq3RDHBXgGuAjYC1ynqvO9Yx29Y4nAS6o6IpyfeUIJAmDxYnjiCXjrLffp8re/wV13QVpazp8rTmWWOPbuDf0hn3FfuXLRnYbYGJM3/EoQrYDdwGuZJIiLgFtxCeJs4BlVPVtEEoHVQDtgEzAH6K6qy7P7mSecIAJ+/hmefhrGjnXDQtu3d4niggvsq6gxJl/KKkFErBezqn4DbMvilC645KGq+l+grIhUBJoBa1V1vaoeBCZ450ZelSrw5JMuUTz6qCtZtGvn+qe+9VZ6dyFjjCkA/BzmUgn4OejxJm9fZvtDEpH+IjJXROZuCV6bLzfKloUhQ1x/0pdecvUsPXrAWWfBqFGwZ0/e/BxjjIlhfiaIUHU2msX+kFR1rKo2UdUmyXm9DFORIm4l8OXL3Yi1SpVg0CA3vPTee12rrTHG5FN+JohNQJWgx5WBzVns909CgpvRb9Ys+PZbaNUKHn7Y9Xa66SbXumuMMfmMnwniA6C3OM2BHar6K65R+iwRSRWRJKCbd25sOPdcmDTJrTTeqxeMG+fWDrzySvjhB7+jM8aYPBOxBCEibwHfAzVEZJOI9BWRASIywDtlKrAeWAu8CNwMoKqHgVuAz4AVwDuquixScZ6wGjVcb6cNG2DoUJg2DZo3h9at4aOPsl+w1hhjYpwNlMsru3bByy/DU0+5XlC1a8Odd7rG7SJF/InJGGOy4Us31wKnVCm4/XY3v8Xrr7uZ0q6/3i2X9dhjbsFdY4yJI5Yg8lrhwtCzJyxc6CY+qlXLdZmtXNn1iFqwwO8IjTEmLJYgIkXEjcT+8ks34O7aa2HCBGjUyE3LOmFCzqYJNcaYKLMEEQ316sHzz7slu0aOdOMnund33WQfeAA2+9uL1xhjQrEEEU1ly7p2ilWrYOpUV5r45z9douje3Y2xyEedBowx8S1HCUJEEkSkdKSCKTASEqBTJ/j4Y1i9Gm67DT75BM47zyWNl19203sYY4yPsk0QIvKmiJQWkRLAcmCViNwV+dAKiGrV3ASBv/wCL7zgFlvo1881at91l1tb1BhjfBBOCaK2qu4ELsMNbjsd6BXJoAqkEiWgf39YtAhmzHBTjI8cCWee6ab5+PxzG3xnjImqcBJEYREpjEsQU1T1EFlMnmdyScTN9fTuu26U9rBhbgqPDh1cl9lRo2DnTr+jNMYUAOEkiBeADUAJ4BsRqQrYJ1Q0VK4Mw4fDxo1u8N1JJ7nZZCtVgoED3SyzxhgTISc01YaIFPLmTIopvk61ES1z5sDo0W4cxYED0Lata9ROSfE7MmNMHMrVVBsiUkFEXhaRT7zHtYFr8zhGE66mTWH8+PRV7+bOdSO3//rL78iMMflMOFVM43Ezq57mPV4N3B6heEy4kpPdFB6jR8N337lJAo0xJg+FkyDKq+o7wBE4Oh23fV2NFT17wuWXwz33wLLYmxXdGBO/wkkQe0TkZLyeS4HFfSIalQmfiJvGo0wZ6N0bDh3yOyJjTD4RToL4O25FtzNF5FvgNeDWiEZlcuaUU1ySmD/fLYVqjDF5oFB2J6jqfBFpDdQABFjljYUwsaRrV7jmGnjoIbjkEmjc2O+IjDFxLpxeTFcBxbxlPy8D3haRRpEOzJyAUaOgQgVX1bR/v9/RGGPiXDhVTPeq6i4ROQ/oALwKjIlsWOaEnHSSGxOxfDnce6/f0Rhj4lw4CSLQY6kzMEZVpwBJkQvJ5ErHjnDjjW4CwFmz/I7GGBPHwkkQv4jIC8DVwFQRKRLmdYhIRxFZJSJrRWRIiOMnicgkEVksIrNFpG7QsQ0iskREFopIPh8enccef9yNrL72Wti92+9ojDFxKpwP+qtxA+U6qup2oByQ7XTfIpIIjAY6AbWB7t4o7GB3AwtVtT7QG3gmw/E2qpqW2TBwk4lSpWDcODdV+ODBfkdjjIlT2SYIVd2rqu8DO0TkdKAwsDKM524GrFXV9ap6EJgAdMlwTm1gmvdzVgIpIlIhJ7+AyUTr1m71un//G774wu9ojDFxKJxeTJeKyBrgR2CGd/tJGM9dCfg56PEmb1+wRUBX7+c0A6oClb1jCnwuIvNEpH8W8fUXkbkiMnfLli1hhFWAPPww1KwJ118P27f7HY0xJs6EU8U0HGgOrFbVVOBC4NswrpMQ+zJOHTsCOElEFuIG3y0AArPEtlDVRrgqqoEi0irUD1HVsaraRFWbJCcnhxFWAVKsGLz6Kvz6qytNGGNMDoSTIA6p6lYgQUQSVHU6kBbGdZuAKkGPKwObg09Q1Z2q2kdV03BtEMm4Egqqutm7/R2YhKuyMjnVrBkMHeoSxZQpfkdjjIkj4SSI7SJSEvgGeENEniH9W35W5gBniUiqiCQB3XBTdhwlImW9YwD9gG9UdaeIlBCRUt45JYD2wNLwfiVznHvvhbQ0t6SpVcMZY8IUToLoAuwF7gA+BdYBl2R3kTfr6y24HlArgHdUdZmIDBCRAd5ptYBlIrISV5U0yNtfAZglIouA2cDHqvpp+L+WOUZSErz2Gvz5J9x0E5zAIlHGmIIn0xXlRKQaUEFVv82wvxXwi6qui0J8OVIgVpTLjREjXHXTm29C9+5+R2OMiQEnuqLc08CuEPv3esdMvLnzTmje3K1nvXlz9ucbYwq0rBJEiqouzrhTVecCKRGLyEROoUKusXr/fujXz6qajDFZyipBFM3iWLG8DsRESfXq8Nhj8MknbmI/Y4zJRFYJYo6I3JBxp4j0BeZFLiQTcQMHQps2cMcdsGGD39EYY2JUVgsG3Q5MEpGepCeEJriZXC+PcFwmkhIS4JVXoH596NMHpk1z+4wxJkimnwqq+puqngs8CGzwtgdV9RxV/V90wjMRk5ICI0fC11/Ds8/6HY0xJgZl2s01Hlk31xxSdcuTTpsGCxdCjRp+R2SMibIT7eZq8jsRePFFN2fTtdfC4XAGyBtjCgpLEAVdxYpuSvAffnALDRljjCfcleGqisiF3v1igXmSTD7xt7/BVVfB/ffDokV+R2OMiRHhrAdxAzAReMHbVRmYHMGYTLSJuFJEuXLQuzccPOh3RMaYGBBOCWIg0ALYCaCqa4BTIhmU8UH58jB2LCxeDP/8p9/RGGNiQDgJ4oC3ZCgAIlKI4xf+MfnBpZfCddfBo4+6NgljTIEWToKYISJ3A8VEpB3wLvBhZMMyvnn6aahUyfVq2rfP72iMMT4KJ0EMAbYAS4AbganAPZEMyvioTBk3ynrVKrj7br+jMcb4KKupNgK6AK+p6ouRDsbEiAsvdPM1Pf00dOkC55/vd0TGGB+EU4K4FFgtIv8Rkc5eG4TJ7x57DKpVg44dYdAg+J/NrmJMQZNtglDVPkA1XNtDD2CdiLwU6cCMz0qUgOnToWdPGD0azjgD/vEP+OMPvyMzxkRJWAPlVPUQ8AkwATeza5dIBmViROXKbs2IFSvgiivgiScgNRXuucetb22MydfCGSjXUUTGA2uBK4GXgIoRjsvEkrPOgv/8B5YuhU6d4OGHXaIYPhx27vQ7OmNMhIRTgrgON3K6uqpeq6pTVTWsWd285LJKRNaKyJAQx08SkUkislhEZotI3XCvNT6oXRveecfN/Hr++XDffS5R/OtfsGeP39EZY/JYOG0Q3VR1sqoeyMkTi0giMBroBNQGuotI7Qyn3Q0sVNX6QG/gmRxca/zSoAFMngyzZ0OzZjB4sGujePppt961MSZfyDRBiMgs73aXiOwM2naJSDj1Cs2Ataq63huJPYHj2y5qA9MAVHUlkCIiFcK81vitaVO3tvWsWVCnjlvCtFo1GDPG5nMyJh/IakW587zbUqpaOmgrpaqlw3juSsDPQY83efuCLQK6AohIM6AqbjLAcK7Fu66/iMwVkblbtmwJIyyT51q0gK++cltKCtx8M1Sv7gbc2RoTxsStcBqp/xPOvlCXhtiXcQ6nEcBJIrIQuBVYABwO81q3U3WsqjZR1SbJyclhhGUipk0bmDnTlSqSk6FvX6hVC15/Hf76y+/ojDE5FE4jdZ3gB95AucZhXLcJqBL0uDKwOfgEVd2pqn1UNQ3XBpEM/BjOtSZGibjBdbNnw5QpbjxFr15Qrx68+y4cOeJ3hMaYMGXVBjFURHYB9YPbH4DfgClhPPcc4CwRSRWRJKAb8EGGn1HWOwbQD/hGVXeGc62JcSJudtj5811iALj6amjY0CWOfLQWujH5VVZtEI+qaing8QztDyer6tDsntjrCnsL8BmwAnhHVZeJyAARGeCdVgtYJiIrcT2WBmV1bS5+T+OXhAS48kpYssRVNe3dC5dd5no/ffqpJQpjYphoNm9QEbkc+EpVd3iPywLnq+rkiEeXQ02aNNG5c+f6HYbJyuHD8NprblGin35yDdwPPWQTAhrjExGZp6pNQh0Lpw3i/kByAFDV7cD9eRSbKWgKFYLrr4fVq90ypz/+6Bq327WD//7X7+iMMUHCSRChzrEZXU3uJCXBTTfB2rXw1FOwaBGccw5ccgksWOB3dMYYwksQc0XkKRE5U0TOEJGRuAn7jMm9YsXcALv16+GRR9ygu0aN4KqrYPlyv6MzpkALJ0HcChwE3sZN+b0PuDmSQZkCqGRJGDrUVTndd59rwK5b13WRXbfO7+iMKZDCmYtpj6oO8QajNQaeBwZGPjRTIJUtCw8+6BLFXXfBe+9BjRrQvz9s3Oh3dMYUKGGtByEi5UXkJhH5BpgOVIhsWKbAK1/erWq3bp2buuPVV92047fdBr/+6nd0xhQIWQ2UKyUivUXkU2A2blW5M1T1TFW9M2oRmoKtYkUYNQrWrIFrr3U9n84801a3MyYKsipB/A70BR4GzlTV/8O1RRgTfaefDmPHwsqVbuBdYHW7++6D7dv9js6YfCmrBHE3UBQYAwwVkTOjE5IxWahWzQ20C6xuN3y4W4vi0Udh926/ozMmX8lqqo2Rqno2cCludtXJwGkiMlhEqkcpPmNCC6xuN3++G419990uUYwcCfv2+R2dMflCOL2Y1qvqw6paD2gKlAE+iXhkxoSjYUP48EP4/nu30t3f/26LFhmTR8LqxRSgqktU9W5VteomE1uaN4cvvoDp013bxM03u3aLYcPcnE/GmBzLUYIwJuadf75btOjzz92MsSNGuKqnSy6BqVNt4SJjcsAShMl/RNzkfx984AbcDR0Kc+ZA586u+mnECPj9d7+jNCbmZTUOYpp3+1j0wjEmj51+uptOfONGePttt2b20KFQuTL06OHmfrI1KYwJKasSREURaQ1cKiINRaRR8BatAI3JE0lJbkW76dPdJIA33eSqnFq2hPr13QC8nTv9jtKYmJLpgkEiciVuoNx5QMZVeFRV20Y4thyzBYNMjuzZAxMmuB5P8+a59bOvucYljwYN/I7OmKjIasGgcFaUu1dVh0cksjxmCcKcsDlzXCliwgTYv9+tTXHTTW7a8aJF/Y7OmIjJ1YpyqjpcRC4VkSe87eK8D9EYnzVtCuPGwS+/uAWMtm6F3r1dW8Vdd9mU46ZAyjZBiMijwCBgubcN8vYZk/+UK+cWMFq5Er780nWbHTnS9X7q0AEmT3brahtTAITTzbUz0E5VX1HVV4CO3r5siUhHEVklImtFZEiI42VE5EMRWSQiy0SkT9CxDSKyREQWiojVG5noEoELLoCJE10PqAcfhGXL4PLL3UC84cNt2nGT74U7DqJs0P0y4VwgIonAaKATUBvoLiK1M5w2EFiuqg2A84EnRSQp6HgbVU3LrH7MmKg47TQ3a+yGDTBpkpsH6r77XBfav/0NZsywrrImXwonQTwKLBCR8SLyKm496kfCuK4ZsNaby+kgMAHokuEcBUqJiAAlgW2Ald9NbCpUCC67DD77DFavhkGD3PQe55/vlkd97jnYscPvKI3JM+E0Ur8FNAfe97ZzVHVCGM9dCfg56PEmb1+w54BawGZgCTBIVY8EfjTwuYjME5H+mf0QEekvInNFZO6WLVvCCMuYPHDWWW5Nil9+cY3bJUrArbdCpUowYAAsWuR3hMbkWlhVTKr6q6p+oKpTVPV/YT63hHqqDI87AAuB04A04DkRKe0da6GqjXBVVANFpFUmsY311stukpycHGZoxuSRYsXguutg9my3XX21Wx41LQ3OOw/efBMOHPA7SmNOSCTnYtoEVAl6XBlXUgjWB3hfnbXAj0BNAFXd7N3+DkzCVVkZE7uaNoVXXknvKvvbb9CzJ1Sp4qb32LDB7wiNyZFIJog5wFkikuo1PHcDPshwzkbgAgARqQDUANaLSAkRKeXtLwG0B5ZGMFZj8k6gq+yqVW5W2RYt4F//Sp9V9pNP4MiR7J/HGJ9lmSBEJEFETuiDWVUPA7cAnwErgHdUdZmIDBCRAd5pw4FzRWQJMA0YrKp/ABWAWSKyCJgNfKyqn55IHMb4JiHBzSo7aZIrPQwb5kZsX3SRa8N4/HH44w+/ozQmU+FMtfEGMFRVN0YnpBNnU22YmHfwoEsY//43fPMNFCniusredBOcfbYbf2FMFOVqqg2gIrBMRKaJyAeBLW9DNKaASEpKHzuxZAn06+cSxjnnQOPG8NJLbhJBY2JAOCWI1qH2q+qMiESUC1aCMHFp1y544w1XqliyBMqUcV1m/+//oGxZv6Mz+VxuJ+ubAWwACnv35wDz8zRCYwqyUqXSx07MnAkXXugWOUpJcbe7dvkdoSmgwpms7wZgIvCCt6sSMDmCMRlTMIm4sRMTJ8LChdC6Ndx7r5v76fHHYe9evyM0BUw4bRADgRbATgBVXQOcEsmgjCnwGjSAKVPc4LumTeEf/3DdZEeNcutVGBMF4SSIA95cSgCISCGOHxFtjImEpk3duIlZs6BWLTf/01lnwQsvuB5RxkRQOAlihojcDRQTkXbAu8CHkQ3LGHOMFi3cetrTprlZZAcMgBo1YPx4W5/CREw4CWIIsAU3md6NwFTgnkgGZYzJRNu2rjQxdSqcfDL06QN16sBbb9nobJPnwunFdAR4FTfq+UHgVc2ub6wxJnJEoFMnNyp70iQ32K5HD6hfH95/39amMHkmnF5MnYF1wCjc9NxrRaRTpAMzxmRDxK1PsXAhTJjgqpquuMINuPv4Y0sUJtfCqWJ6Erey2/mq2hpoA4yMbFjGmLAlJLjR2UuXuqnGd+yAiy92o7O//NIShTlh4SSI372puAPWA79HKB5jzIkqVAh694aVK2HsWNi82U0WeP75bt4nY3Io0wQhIl1FpCtuHqapInKdiFyL68E0J2oRGmNypnBhuOEGWLMGnn3WLY/aujW0bw8//OB3dCaOZFWCuMTbigK/Aa2B83E9mk6KeGTGmNwpUgRuuQXWrXPLoy5YAM2bQ+fO8N13fkdn4kC2k/XFE5usz5gs7N7tShRPPglbt7pSxdChrmRh04wXWLmarM9bEe4pEXnfpvs2Jo6VLOkSwk8/wdNPw9q10LEjNGkC771n4yjMccJppJ6Mm831WVyPpsBmjIlHJUq4KTvWrXPrT+zcCVde6QbcjR8Phw75HaGJEeEkiP2qOkpVp6vqjMAW8ciMMZFVpAj07et6Pb39tnvcpw9UqwbPPQf79vkdofFZOAniGRG5X0TOEZFGgS3ikRljoiMxEa6+2jViT50KVaq4BYtSUuDRR924ClMghZMg6gE3ACNIr156IpJBGWN8EJjCY9YsN26iUSO4+243OeCwYfC7DX8qaMJJEJcDZ6hqa1Vt421tw3lyEekoIqtEZK2IDAlxvIyIfCgii0RkmYj0CfdaY0wEtWzpphmfPx86dHAliZQUuO022LjR7+hMlISTIBYBZXP6xCKSCIwGOgG1ge4iUjvDaQOB5araADfG4kkRSQrzWmNMpDVsCO+8AytWQPfuMGYMnHmma6tYudLv6EyEhZMgKgArReSzHHZzbQasVdX13oJDE4AuGc5RoJSICFAS2AYcDvNaY0y01KgBL7/sej7dfLNr1K5dG666ypUyTL4UToK4H1fN9Ag56+ZaCfg56PEmb1+w54BawGbcehODvOnFw7kWABHpLyJzRWTuli1bwgjLGHPCTj8dnnnGjaW4+2744gs3e2zHjq7dIh8NvDXhrQcxI9QWxnOHGpqZ8b+nA7AQOA1IA54TkdJhXhuIb6yqNlHVJsnJyWGEZYzJteRkeOghlyhGjHA9oFq3hvPOg4kTbZW7fCKckdS7RGSnt+0Xkb9EZGcYz70JqBL0uDKupBCsD/C+OmuBH4GaYV5rjPFbmTIweDBs2ODGTmze7KqdUlPhkUfASvVxLZwSRClVLe1tRYErcFVD2ZkDnOVN1ZEEdAMytl1sBC4AEJEKQA3cdOLhXGuMiRXFisHAgW76jilToGZN1zW2ShW47jqwOdLiUjhtEMdQ1clAtt1cVfUwcAvwGbACeEdVl4nIABEZ4J02HDhXRJYA04DBqvpHZtfmNFZjTJQlJsKll7q2ieXLoV8/N89T06ZuAaM334SDB/2O0oQp29lcvTUhAhKAJkBrVT0nkoGdCJvN1ZgYtGOHW+nuuefcGhWnngo33ui2ihX9jq7Ay9VsrqSvC3EJrlF5F9bl1BgTrjJl3AC7lSvd4LtGjeCf/3Q9onr0gO+/t95PMcrWgzDGRN/atfDvf8Mrr7gSRqNGbv6nbt2gaFG/oytQsipBZJogROS+LJ5TVXV4XgSXlyxBGBNndu+G1193CxktXw7ly7vlUm+6yTVwm4g70SqmPSE2gL7A4DyN0BhTMJUsCQMGwNKlMG2aG0fx2GOum+yVV8LXX1v1k48yTRCq+mRgA8YCxXDjFiYAZ0QpPmNMQSACbdvCpEmwfj3ceSdMnw5t2kCDBjB2LOzZk/3zmDyVZSO1iJQTkYeAxUAhoJGqDlZVm/fXGBMZVau60dmbNrn5nxITXY+nypXd2ApbnyJqMk0QIvI4bsDaLqCeqj6gqn9GLTJjTMFWrBhcf72bDHDmTLjgAjc6+4wzYORIOHDA7wjzvawaqY8AB3CzqwafJLhG6tKRDy9nQjVSHzp0iE2bNrF//36fojIBRYsWpXLlyhQuXNjvUEy8mj/fTe3x5ZeupPHQQ66rbEKOx/wazwn1YopHoRLEjz/+SKlSpTj55JNxs4obP6gqW7duZdeuXaSmpvodjol3X3zhEsWCBa6NYsQIt7CRvcdzLLcD5eLa/v37LTnEABHh5JNPtpKcyRvt2rn5nd58E3budEulXnABzJnjd2T5Sr5PEIAlhxhhfweTpxIS3Cp3K1fCqFGwZAk0awZ/+5sbiGdyrUAkCGNMPpaU5EZhr1sH994LH30EtWq52WV/+83v6OKaJYgo+O233+jRowdnnHEGjRs35pxzzmHSpElRjWHDhg3UrVs35P4333zzhJ7z6aefZu/evUcflyxZ8oTjMybXSpd2czytW+dGY7/wgls/+4EHYNcuv6OLS5YgIkxVueyyy2jVqhXr169n3rx5TJgwgU2bNh137mEfVuHKKkFkF0/GBGFMTDj1VDfP0/Llrm3iwQehWjUYPdqmGs8hSxAR9tVXX5GUlMSAAQOO7qtatSq33norAOPHj+eqq67ikksuoX379mzbto3LLruM+vXr07x5cxYvXgzAAw88wBNPPHH0OerWrcuGDRvYsGEDtWrV4oYbbqBOnTq0b9+effv2ATBv3jwaNGjAOeecw+jRo0PGN2TIEGbOnElaWhojR448Lp6vv/6aiy+++Oj5t9xyC+PHj2fUqFFs3ryZNm3a0KZNm6PHhw0bRoMGDWjevDm/WfHe+Kl6dXj3Xfjvf12V0y23QO3a8PbbcOSI39HFhUJ+BxBVt98OCxfm7XOmpcHTT2d6eNmyZTRq1CjLp/j+++9ZvHgx5cqV49Zbb6Vhw4ZMnjyZr776it69e7Mwm5jXrFnDW2+9xYsvvsjVV1/Ne++9xzXXXEOfPn149tlnad26NXfddVfIa0eMGMETTzzBRx99BLiEFRzP119/HfK62267jaeeeorp06dTvnx5APbs2UPz5s15+OGH+cc//sGLL77IPffck2XsxkTc2We7aTs++QSGDHEzxj7xhJvzqW22a58VaFaCiLKBAwfSoEEDmjZtenRfu3btKFeuHACzZs2iV69eALRt25atW7eyI5upBVJTU0lLSwOgcePGbNiwgR07drB9+3Zat24NcPQ5wxEcT04kJSUdLW0E4jAmJojARRe5cROvvgq//+66xXbsCIsW+R1dzCpYJYgsvulHSp06dXjvvfeOPh49ejR//PEHTZqkj0spUaLE0fuhBi6KCIUKFeJIULE4eDxBkSJFjt5PTExk3759qOoJdysNjiern5tR4cKFj/7MxMREX9pUjMlSYiL07g1XX+3aKR5+GBo2hJ49YfhwSEnxO8KYYiWICGvbti379+9nzJgxR/dl1bDbqlUr3njjDQC+/vprypcvT+nSpUlJSWH+/PkAzJ8/nx9//DHLn1u2bFnKlCnDrFmzAI4+Z0alSpViVxY9PKpWrcry5cs5cOAAO3bsYNq0aWFfa0zMKloU/v531+Np8GCYONE1ZHfuDBMmgNeOV9BZgogwEWHy5MnMmDGD1NRUmjVrxrXXXstjjz0W8vwHHniAuXPnUr9+fYYMGcKrr74KwBVXXMG2bdtIS0tjzJgxVK9ePdufPW7cOAYOHMg555xDsWLFQp5Tv359ChUqRIMGDRg5cuRxx6tUqcLVV19N/fr16dmzJw0bNjx6rH///nTq1OmYRmpj4krZsvDoo25g3V13weLFbvBdhQrQt69bj6IAN2hHdC4mEekIPAMkAi+p6ogMx+8CenoPCwG1gGRV3SYiG3Azyf4FHM5srpBgoeZiWrFiBbVq1crtr2LyiP09TEw7cgRmzIDXXnOlit273drZ11wDvXpBzZp+R5jnfJmLSUQSgdFAJ6A20F1Eagefo6qPq2qaqqYBQ4EZqrot6JQ23vFsk4MxxuRaQoJbpGjcODcK+403XNfYESNcV9lmzeC55+CPP/yONCoiWcXUDFirqutV9SBuJbouWZzfHXgrgvEYY0z4ihd3U4l/8olbvOjJJ+HQITetR8WKcOmlrpSRjyegjGSCqAT8HPR4k7fvOCJSHOgIvBe0W4HPRWSeiPTP7IeISH8RmSsic7ds2ZIHYRtjTAYVK7pG7QULXDvFHXfAvHlw1VVu5PaNN8KsWflu/exIJohQfSwze/UuAb7NUL3UQlUb4aqoBopIq1AXqupYVW2iqk2Sk5NzF7ExxmSnXj34179g40b4/HNXknj9dWjZ0s39dP/9+WY22UgmiE1AlaDHlYHNmZzbjQzVS6q62bv9HZiEq7IyxpjYkJjo1qV47TXXXvHaay5BDB8OZ50F554LY8bAtm3ZP1eMimSCmAOcJSKpIpKESwIfZDxJRMoArYEpQftKiEipwH2gPbA0grEaY8yJK1nS9XL64gv4+WdXwti1C26+2VVBde0KH34IcTZ4NGIJQlUPA7cAnwErgHdUdZmIDBCRAUGnXg58rqp7gvZVAGaJyCJgNvCxqn4aqVgjLTExkbS0NOrWrctVV10VszOgZjYluDEmBypVSh9TsWCBa9T+7jtXFXX66TBsmBugFwciOlBOVaeqanVVPVNVH/b2Pa+qzwedM15Vu2W4br2qNvC2OoFr41WxYsVYuHAhS5cuJSkpieeff/6Y43/99ZdPkRljIkbETeb55JOuVDF5MjRu7LrMVqvm5oJ6662Y7gVlI6mjrGXLlqxdu5avv/6aNm3a0KNHD+rVq8f+/fvp06cP9erVo2HDhkyfPh1ws6t26dKFjh07UqNGDR588MGjz/XUU09Rt25d6taty9PePFN79uyhc+fONGjQgLp16/L2228Dburv1q1b07hxYzp06MCvv/56dH92U4IbY3KpcGHo0sVVM23cCA89BD/+6LrRnnYaDBrklkyNMQVqsj4fZvs+xuHDh/nkk0/o2LEjALNnz2bp0qWkpqby5JNPArBkyRJWrlxJ+/btWb169THnFS9enKZNm9K5c2dEhHHjxvHDDz+gqpx99tm0bt2a9evXc9ppp/Hxxx8DsGPHDg4dOsStt97KlClTSE5O5u2332bYsGG88sorYU0JbozJQ5UquWqmoUPhq6/gpZfg+efdutpnnw39+rl1tUuV8jtSK0FEw759+0hLS6NJkyacfvrp9O3bF4BmzZqRmpoKHDvNd82aNalaterRBNGuXTtOPvlkihUrRteuXZk1axazZs3i8ssvp0SJEpQsWZKuXbsyc+ZM6tWrx5dffsngwYOZOXMmZcqUYdWqVSxdupR27dqRlpbGQw89xKZNm3I1JbgxJpcSEuDCC93kgL/8AiNHuqk9brjBjbvo188tduTj2IoCVYLwYbZvIL0NIqPspvkOyDhtt4hken716tWZN28eU6dOZejQobRv357LL7+cOnXq8P333x9z7vbt2094SnBjTB4qX95VcQwaBD/84EoVEybAyy9DnTouWVxzjTsviqwEESOCp/levXo1GzdupEaNGgB88cUXbNu2jX379jF58mRatGhBq1atmDx5Mnv37mXPnj1MmjSJli1bsnnzZooXL84111zDnXfeyfz586lRowZbtmw5miAOHTrEsmXLwp4S3BgTJSLQvLlLEL/+Ci++6LrQ3nGHq5rq1g2+/DJqM8wWqBJELLv55psZMGAA9erVo1ChQowfP/7oQkDnnXcevXr1Yu3atfTo0ePoYkPXXXcdzZq58YP9+vWjYcOGfPbZZ9x1110kJCRQuHBhxowZQ1JSEhMnTuS2225jx44dHD58mNtvv506deowbtw4rr/+eooXL06HDh18+/2NMRmUKuVKDv36uS6zL78M//mPW1M7JcVNR37ddVC5csRCiOh039GWH6f7Hj9+PHPnzuW5557zO5Q8Ee9/D2N8tX8/TJrkShhffeXaMTp1cu0WF1/sRnfnkC/TfRtjjMljRYu6BY2mTXPzPQ0ZAvPnuxHbEfiybyUIE1X29zAmjx0+7EZme22WOVXgSxD5KQnGM/s7GBMBhQqdcHLITr5PEEWLFmXr1q324eQzVWXr1q0ULVrU71CMMWHK972YKleuzKZNm7DFhPxXtGhRKkewx4UxJm/l+wRRuHDho6OVjTHGhC/fVzEZY4w5MZYgjDHGhGQJwhhjTEj5ahyEiGwBfvI7jhDKA3/4HcQJstj9YbFHX7zGDbmLvaqqJoc6kK8SRKwSkbmZDUSJdRa7Pyz26IvXuCFysVsVkzHGmJAsQRhjjAnJEkR0jPU7gFyw2P1hsUdfvMYNEYrd2iCMMcaEZCUIY4wxIVmCMMYYE5IlCGOMMSFZgjDGGBOSJYhcEOcfIhJ3S6RZ7P6w2KMvXuMG/2O3BHGCRKQK8DEwAugtIiV8DilsFrs/LPboi9e4ITZitwRx4hKAfwG1gcZAPA3Rt9j9YbFHX7zGDTEQe75fMCiviMgpQH/gZ+AdVf1JRLao6l4RmYHL8EtVdau/kR7PYveHxR598Ro3xGbsVoIIg4gMBKYDJYBrgOEicoqq7vVOeQqoAnQQkcSg6yTqwWZgsfvDYo++eI3biyEmY7cEkQ0RKQycDAxQ1aHAjcDZQJJ3vJCq7gNeAy4FiovIhSJSUn0epm6x+8Nit7hzIqZjV1XbMmxAGtANOMN7fKp3W9i7/RhoHeK6GcAvwBzgTIvdYrfYLe54jj3qL0ysb8BNwHrgGWABcH6G46WAhUDloH2FgN7ADuBmi91it9gt7vwQuy8vUCxvwH+Adt79G4AvgFTvseB6ErwXdH4F77Y+UDxof4LFbrFb7BZ3PMde4NsgRCRNRE7z7gvwO3CSiIiqvghsBK4SkcLq/iqnAmtEpKWIfA909s5drK63QQKAqh6x2C12i71gxx3vsUc1c8bSBqQCHwD/BSYBt3j7HwX+DyjtPW6Mq+8r5z1+ETgIfARcbrFb7Ba7xZ3fYg9sBa4EISIJIpKMq//7TlWb44p8tcWNXJwItAVqikgRVZ0H/Apc6z3FPOCfqnqxqk6y2C12i93izi+xH8fP7ORTVn8NaAnUAkp5+04CfgBqeI8fwo1gDNQTjsJrSMJbZMm7H+16V4vdYi8Qscdr3PEee8atQJQgRKR40MPPgb6qukJVd4lIErAT2AIU8855ElgFDBaRb4GmwEoA9f5q3v3I1wG6+AMDYyz26MYeeH/EVexBr3lcxW7vU/9iz5Sf2SkKmfxUXHFuWtC+04GHgdOD9tXAFQUzXl8baOlT7BWBrkCZoH1V4yT2k4CTvfuBZW3j6XXvBlSJw9f9NGB48M+Ph9fd3qf+xB7Olm9LECJyLzAVOAT8LiJVvUNFgJK4jB5QF5jmXTdURPoDqOpyVZ3p7Y/aayUitwOfAh2BZ0XkSu9QcWI/9geAzcDVIlJcvXcB8fG6P4j7n7kAGCsiXb1DReMk9o+Aw4Gf74np193ep/7EHq6YCygviMhtuG9TFwB34uY32QWgqmuASkDnoEsqAt1EZCau3vC9jM+p0S2mNgD6qWp/4HXgFREpraorgMrEbuxn4d7YzwLVvHgCMcT6694ZKAt0UNUbgFnAKV4Mq4jt170V0AV4SlUfzBDDGtwcPjEXu4gMIn7fp4WI0/dpjvhdhMmrjWMHkBTLcGw50DXocVvgedKrPz4GFgGNg86RSMWaWey4QTLVgK+B5MAxXA+HEd7j9jEWewpQybtfBCgHJHox3o7XdS9GX/cUvNGqQJGg/U2BtbgRr2fHweteznut7wbOB8Z6sXf2jneIldgze829x7H+Pk0BTvPup+KmvYiL9+mJbnE/3beIlAFeBqrjRhoC7PeOFVbVQ8Ab3vGAbcBPQHlcw1EvVd3mXSNwbENRtGL3fuZaEfkf8JCIvA+0Bt4BrhORkbgh+j/HQOwn4WaYrAP8KSKvAB8GxfIubmqAJSLytar+BWwnNl73ULF/BBwQkWq40a2vAnuBt0SkJbAa2BSDsY/D1d/PAh7BTRc9EtcYOk5EGuIaQ32NPZPX/GPca56kqgeJ3ffpca+5qk4QkZ+Ah0XkPWL0fZpb+aGKaRCuzltEZLC3LzDS8JD3uCQuwwdsw31LPOg93oF7ggT1RDxqJ1Ts4L4NzgN6AUVVdRAwGTcEfy9uYI1vsXu9Ml4GdqpqM+BNoBWuSgAAVZ2G+0bVGleiADdi9OwYj30dblbN4ar6Ku5b4iBcfXKjGIy9Na46YwluAFZ9VX1GVZ/CDdIaiKu28S32LF7zigBecoAYfJ9mEnsbceMc/gHMJUbfp3nC7yLMiW64P1o34CTvcRrwJ+nVNQl4fYhxRdVVGa4fD3SLxdiDzguu9ngDqObdf9Xn2HvgituBkZ+FgaVAQ+9xond7Cq664xnct9hAiSlmYw9xzSPAFTH+ujcOPM5wzWPAZX7GnsP/l1h8n2b5mnv7Yu59mmevgd8B5PAPFjyApB3wQobj7wCvefcLBe2viuspkRa0LymGYw9M+VsEaIibzGsS6VMC+x37WNKTb1HvdgpQN8N1FXHfAufh1S2Tod45FmPHfZNt5P1NppM+uCnmX3fcrJ8NgHf9ij0X/y+x+D4N5zWPifdpJLa4qGISkeIiMgRX3Az4FfhJRMoHdQ/rD3QRkdqqetgrHgIo8LiqLgxcrOnF2liM/ZDXmykRuAj4XFUvV9X/xUjsG3ANo6jqfq+4fTLwY4anuAj3ujdW1fe9ffEQe0VgGDBLVduo68UUL697JVzVx8xox36icQfq5IEjxN77dAPZv+aF8fF9Gkkx30gtIjcB/XATXi0OOrQbVzVzSFWPiJvTZLuI3A+8ICIvAhVF5ElV3Yir/46r2HFD8R/RwFcbV38ZrW582cYetK8lMFtV94jIjcAe4A1VfTno+RJU9Ujgd4nh2Her6hsi0k29Nqw4et0DsV+rqoejGXtu/19E5C1V/RnXsBtVefSaPxp4naP5/xJpMV2CEJEWuB4lL6nqQFXdHzimqhtw9fY3eLsC2XoR0AL3B58ceKNEWx7FrqqqQb0eovUhlZPYAZKBc0XkU+Ay4NuMiSCOYv/OO/dQHL7ugdgPRzP2PIh7lrpeblGXh6/5kWj/v0RDzJUgxA22ug73YTkNV9+3T0Sa4Ya0rwHWquoM3OCUjiJSTFX3eX/sv+Maht7xnk+i8a01krFH6Vv3icReQlX34JJaVeAeVf0y0rFGK/Z4ft0jHbv9v/jz/xJtMVWCEJE6uImu9gE3A31w364vBN7HdR2rCjwuImm4BtC/SC8Cfq+ql/iUHApi7Ae8p7hTVSsG3jBBdcoWez6MPV7jjvfY/SCxlPRE5HqgtqreKSLVcQ0/qcBs4FNV3eqddz+uv/R9uG8APVX1p0DdXzQ/XC12/SkQs4gkRruqwGKPfuzxGne8x+6HmCpB4EYfni1uBPRq3CpLh3HdzLYGnVccWKSqB4BvgE6QXvfnU1GvIMeu3q0fbxiLPfqxx2vcEN+xR12sJYgVuEEoV3uPF+PmxKkkIoVF5EwReR7Xk2CFd85Dqvp89EM9jsXuD4s9+uI1bojv2KMu1hLEFlxRr62IVFLVXbg6wbpAadyI3F9V9VxVXQCgqnt9i/ZYFrs/LPboi9e4Ib5jj7qYShBeNctHuD/i497uRGCPV/zrpt50xhJjc6db7P6w2KMvXuOG+I7dDzHVSB0gIkVwM1Qm4Ka/7q6q871jUW/EzQmL3R8We/TFa9wQ37FHU0wmCDj6B0xW1U1+x5JTFrs/LPboi9e4Ib5jj5aYTRDBJI6Hrlvs/rDYoy9e44b4jj2S4iJBGGOMib4C3whjjDEmNEsQxhhjQrIEYYwxJiRLEMYYY0KyBGGMMSYkSxDGGGNC+n9VItgsL3Tu0AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.dates as mdates\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",
    "dates = [test_start_date + pd.DateOffset(days=i+history_window) for i in range(pred_window)]\n",
    "#dates = [valid_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(dates, I_true, c='r', label='Ground truth')\n",
    "plt.plot(dates, pred_I[-1, :],c='b', label='Proposed')\n",
    "# stan_res = pickle.load(open(f\"stan_results/{loc_name}_stan_results_new\", \"rb\"))\n",
    "# world_pred = pickle.load(open(f\"results/{loc_name}_world_input\", 'rb'))\n",
    "# plt.plot(dates, world_pred[-1, :], c='purple', label=\"After Calibration\")\n",
    "# stan_sirv_res = pickle.load(open(f\"stan_results/{loc_name}_stan_results_sirv\", \"rb\"))\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.gcf().autofmt_xdate()\n",
    "plt.legend()\n",
    "plt.ylabel(\"Number of Active Cases\")\n",
    "plt.title(loc_name)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Proposed: 3505110628934.5713\n"
     ]
    }
   ],
   "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}\")\n",
    "# loc_list.index(\"Germany\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "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):\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",
    "\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",
    "\n",
    "# nx_g = model.g.cpu().to_networkx(edge_attrs=['e'])\n",
    "# #print(nx_g.nodes())\n",
    "# attention = []\n",
    "# for e in sorted(nx_g.edges()):\n",
    "#     attention.append(nx_g.get_edge_data(e[0], e[1])[0]['e'].item())\n",
    "# attention = np.array(attention)\n",
    "# attention = (attention - np.min(attention)) / np.ptp(attention)\n",
    "# labels = [0 for i in range(len(nx_g.nodes()))]\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\" for i in range(len(nx_g.nodes()))])\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.12 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
}