mimic_ope.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from tqdm.notebook import tqdm\n",
    "import pickle\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from src.offline_pg import ISPG\n",
    "from src.offline_ql import discrete_PQL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute IS weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/s45da_mimic_train_episodes\", 'rb') as f:\n",
    "    train_dataset = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/s45da_mimic_val_episodes\", 'rb') as f:\n",
    "    val_dataset = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/s45da_mimic_test_episodes\", 'rb') as f:\n",
    "    test_dataset = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# computer IS weights: select pibs, compute pies\n",
    "def compute_is_weights(dataset, q_policy):\n",
    "    n = dataset[\"observations\"].shape[0]\n",
    "    horizon = dataset[\"observations\"].shape[1]\n",
    "    \n",
    "    rewards = dataset[\"rewards\"][:,:,0]\n",
    "    \n",
    "    pies = np.zeros((n, horizon))\n",
    "    pibs = np.zeros((n, horizon))\n",
    "    length = np.zeros((n))\n",
    "    \n",
    "    for idx in range(n):\n",
    "        for t in range(horizon):\n",
    "            pie = q_policy.get_prob(dataset[\"observations\"][idx,t,:], dataset[\"estm_pibs\"][idx,t,:], \n",
    "                                    dataset[\"nn_action_dist\"][idx,t,:])\n",
    "            a = dataset[\"actions\"][idx,t,0]\n",
    "            \n",
    "            pies[idx,t] = pie[a]\n",
    "            pibs[idx,t] = dataset[\"pibs\"][idx,t,a]\n",
    "    \n",
    "            if not dataset[\"not_done\"][idx,t,0]:\n",
    "                length[idx] = t+1\n",
    "                break\n",
    "    return pibs, pies, rewards, length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wis_ope(pibs, pies, rewards, length, weighted=True):\n",
    "    n = pibs.shape[0]\n",
    "    weights = np.ones((n,20))\n",
    "    \n",
    "    for i in range(n):\n",
    "        last = 1\n",
    "        for t in range(int(length[i])):\n",
    "            weights[i,t] = last*(pies[i,t]/pibs[i,t])    \n",
    "            last = weights[i,t]\n",
    "        weights[i,length[i]:] = weights[i,length[i]-1]\n",
    "    weights = np.clip(weights, 0, 1e3)\n",
    "    if weighted:\n",
    "        weights_norm = weights.sum(axis=0)\n",
    "    else:\n",
    "        weights_norm = weights.shape[0]\n",
    "    weights /= weights_norm\n",
    "    \n",
    "    return (weights[:,-1] * rewards.sum(axis=-1)).sum(axis=0), weights[:,-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import norm, t\n",
    "def bca_bootstrap(pibs, pies, rewards, length, alpha=0.05, B = 1000):\n",
    "    wis_list = []\n",
    "    n_patients = pibs.shape[0]\n",
    "    n_subsample = n_patients\n",
    "    for b in range(B):\n",
    "        ids = np.random.choice(n_patients, n_subsample)\n",
    "        sam_rewards = rewards[ids,:]\n",
    "        sam_pibs = pibs[ids,:]\n",
    "        sam_pies = pies[ids,:]\n",
    "        sam_length = length[ids]\n",
    "        wis_pie, _ = wis_ope(sam_pibs, sam_pies, sam_rewards, sam_length)\n",
    "        wis_list.append(wis_pie)\n",
    "    y = []\n",
    "    for i in range(n_patients):\n",
    "        sam_rewards = np.delete(rewards,i,axis=0)\n",
    "        sam_pibs = np.delete(pibs,i,axis=0)\n",
    "        sam_pies = np.delete(pies,i,axis=0)\n",
    "        sam_length = np.delete(length,i,axis=0)\n",
    "        wis_pie, _ = wis_ope(sam_pibs, sam_pies, sam_rewards, sam_length)\n",
    "        y.append(wis_pie)\n",
    "        \n",
    "    wis_list = np.array(wis_list)\n",
    "    wis_list = np.sort(wis_list)\n",
    "    y = np.array(y)\n",
    "    avg,_ = wis_ope(pibs, pies, rewards, length)\n",
    "    \n",
    "    ql, qu = norm.ppf(alpha), norm.ppf(1-alpha)\n",
    "\n",
    "    # Acceleration factor\n",
    "    num = np.sum((y.mean() - y)**3)\n",
    "    den = 6*np.sum((y.mean() - y)**2)**1.5\n",
    "    ahat = num / den\n",
    "\n",
    "    # Bias correction factor\n",
    "    zhat = norm.ppf(np.mean(wis_list < avg))\n",
    "    a1 = norm.cdf(zhat + (zhat + ql)/(1-ahat*(zhat+ql)))\n",
    "    a2 = norm.cdf(zhat + (zhat + qu)/(1-ahat*(zhat+qu)))\n",
    "\n",
    "    print('Accel: %0.3f, bz: %0.3f, a1: %0.3f, a2: %0.3f' % (ahat, zhat, a1, a2))\n",
    "    return np.quantile(wis_list, [a1, a2]), wis_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bca_bootstrap_mc(x, alpha=0.05, B=1000):\n",
    "    values = []\n",
    "    n_patients = x.shape[0]\n",
    "    n_subsample = n_patients\n",
    "    for b in tqdm(range(B)):\n",
    "        ids = np.random.choice(n_patients, n_subsample)\n",
    "        value = x[ids,:].mean()\n",
    "        values.append(value)\n",
    "    y = []\n",
    "    for i in range(n_patients):\n",
    "        x_sub = np.delete(x,i,axis=0)\n",
    "        value = x_sub.mean()\n",
    "        y.append(value)\n",
    "        \n",
    "    values = np.array(values)\n",
    "    values = np.sort(values)\n",
    "    y = np.array(y)\n",
    "    avg = x.mean()\n",
    "    \n",
    "    ql, qu = norm.ppf(alpha), norm.ppf(1-alpha)\n",
    "\n",
    "    # Acceleration factor\n",
    "    num = np.sum((y.mean() - y)**3)\n",
    "    den = 6*np.sum((y.mean() - y)**2)**1.5\n",
    "    ahat = num / den\n",
    "\n",
    "    # Bias correction factor\n",
    "    zhat = norm.ppf(np.mean(values < avg))\n",
    "    a1 = norm.cdf(zhat + (zhat + ql)/(1-ahat*(zhat+ql)))\n",
    "    a2 = norm.cdf(zhat + (zhat + qu)/(1-ahat*(zhat+qu)))\n",
    "\n",
    "    print('Accel: %0.3f, bz: %0.3f, a1: %0.3f, a2: %0.3f' % (ahat, zhat, a1, a2))\n",
    "    return np.quantile(values, [a1, a2]), values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot Result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pick the Best Model on Validation Set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "ESS_lb = 200"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_model_from_param(param, ess_constraint=True):\n",
    "    \n",
    "    FEATURES_UB = [20, 0.5000, 0.5000, 5.1672, 2.7543, 7.9158, 2.7949, 5.2763, 7.1967,\n",
    "                      9.0866, 11.8737, 7.8319, 236.7704, 3.0286, 9.9485, 8.7237, 6.2627, 15.6281,\n",
    "                      8.5191, 13.5491, 6.1870, 40.5525, 10.0892, 6.7538, 20.4172, 8.5691, 6.5705,\n",
    "                      12.9162, 6.6171, 17.8618, 29.6399, 5.0306, 11.4043, 40.8348, 5.0392, 2.3167,\n",
    "                      23.0958, 7.7701, 7.5767, 5.6300, 4.9695, 4.2235, 1.7962, 2.8099, 2.0322]\n",
    "    \n",
    "    if param[\"name\"] in {\"POELA\", \"POELA-IS\"}:\n",
    "        settings = f\"PG_{param['env']}_masknn_action_dist_weightedTrueam0.0\"\\\n",
    "                   f\"thresh{param['thresh']}var{param['var_coeff']}clip0.01seed{param['seed']}\"\n",
    "        model = ISPG(state_dim=45, num_actions=25, device=\"cpu\", \n",
    "                     threshold=10000, action_mask_type=\"nn_action_dist\") #threshold=param['thresh']     \n",
    "    \n",
    "    elif param[\"name\"] in {\"norm-POEM\", \"norm-IS\", \"POELA-kNN\", \"POELA-kNN-IS\", \"norm-POEM*\", \"norm-IS*\"}:\n",
    "        settings = f\"PG_{param['env']}_maskstep_weightedTrueam0.0\"\\\n",
    "                   f\"thresh{param['thresh']}var{param['var_coeff']}clip0.01seed{param['seed']}\"\n",
    "        model = ISPG(state_dim=45, num_actions=25, device=\"cpu\", \n",
    "                     threshold=param['thresh'], action_mask_type=\"step\")\n",
    "    \n",
    "    elif param[\"name\"] == \"BCQ\":\n",
    "        settings = f\"BCQ_{param['env']}_thresh{param['thresh']}_seed{param['seed']}\"\n",
    "        model = discrete_PQL(num_actions=25, state_dim=45, device=\"cpu\", state_clipping=False, \n",
    "                             max_state=FEATURES_UB, action_threshold=param['thresh'], log_pibs=False)\n",
    "    \n",
    "    elif param[\"name\"] == \"PQL\":\n",
    "        settings = f\"PQL_{param['env']}_thresh{param['thresh']}_seed{param['seed']}\"\n",
    "        model = discrete_PQL(num_actions=25, state_dim=45, device=\"cpu\", state_clipping=True, \n",
    "                             max_state=FEATURES_UB, action_threshold=param['thresh'], log_pibs=False)\n",
    "        \n",
    "    if \"POELA\" in param[\"name\"] or \"norm-\" in param[\"name\"]:\n",
    "        val_score = np.load(f\"./results/{settings}_val_wis_scores.npy\")\n",
    "        train_score = np.load(f\"./results/{settings}_train_wis_scores.npy\")\n",
    "        #print(param['thresh'], param['var_coeff'], param['seed'], train_score.shape, val_score.shape)\n",
    "        #print(np.max(val_score[:,0]), np.max((val_score[:,-1]>200)*val_score[:,0]), val_score[-1,0])\n",
    "        if ess_constraint and np.sum(val_score[:,-1]>200):\n",
    "            model.load(f\"models/{settings}_best_\")\n",
    "            train_score = train_score[np.argmax((val_score[:,-1]>200)*val_score[:,0]),:]\n",
    "            val_score = val_score[np.argmax((val_score[:,-1]>200)*val_score[:,0]),:]\n",
    "        else:\n",
    "            model.load(f\"models/{settings}_best_noess_\")\n",
    "            train_score = train_score[np.argmax(val_score[:,0]),:]\n",
    "            val_score = val_score[np.argmax(val_score[:,0]),:]\n",
    "    else:\n",
    "        val_score = np.load(f\"./results/{settings}_val_scores.npy\")\n",
    "        train_score = np.load(f\"./results/{settings}_q_scores.npy\")\n",
    "        #print(np.argmax(val_score[:,0]), np.argmax((val_score[:,-1]>200)*val_score[:,0]))\n",
    "        if ess_constraint and np.sum(val_score[:,-1]>200):\n",
    "            model.load(f\"models/{settings}_best_\")\n",
    "            train_score = train_score[np.argmax((val_score[:,-1]>200)*val_score[:,0])]\n",
    "            val_score = val_score[np.argmax((val_score[:,-1]>200)*val_score[:,0]),:]\n",
    "        else:\n",
    "            model.load(f\"models/{settings}_best_noess_\")\n",
    "            train_score = train_score[np.argmax(val_score[:,0])]\n",
    "            val_score = val_score[np.argmax(val_score[:,0]),:]\n",
    "    \n",
    "    return model, val_score, train_score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "seeds = [1,2,3,4,5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best wis: 97.38  ESS: 91.02 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.4, 'var_coeff': 0.0, 'seed': 1}\n",
      "best wis: 97.87  ESS: 77.87 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.4, 'var_coeff': 0.0, 'seed': 3}\n",
      "best wis: 98.5  ESS: 85.41 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.4, 'var_coeff': 0.1, 'seed': 3}\n",
      "best wis: 99.25  ESS: 88.19 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.4, 'var_coeff': 1.0, 'seed': 3}\n",
      "reliable wis: 93.6  ESS: 376.5 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.6, 'var_coeff': 0.0, 'seed': 1}\n",
      "reliable wis: 93.87  ESS: 374.5 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.6, 'var_coeff': 0.0, 'seed': 4}\n",
      "reliable wis: 94.01  ESS: 407.28 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.6, 'var_coeff': 0.1, 'seed': 2}\n"
     ]
    }
   ],
   "source": [
    "# Our method\n",
    "var_coeff_list = [0.0, 0.1, 1.0, 10.0]\n",
    "threshold_list = [0.4, 0.6, 0.8, 1.0]\n",
    "\n",
    "best_wis = -100\n",
    "best_reliable_wis = -100\n",
    "for threshold in threshold_list:\n",
    "    for var_coeff in var_coeff_list:\n",
    "        for seed in seeds:\n",
    "            param = {\"name\": \"POELA\", \"env\": \"mimic45\", \"thresh\": threshold, \"var_coeff\": var_coeff, \"seed\": seed}\n",
    "            model, score, _ = get_model_from_param(param)\n",
    "        \n",
    "            wis_pie = score[0]\n",
    "            ESS = score[-1]\n",
    "        \n",
    "            if ESS >= ESS_lb and wis_pie > best_reliable_wis:\n",
    "                best_reliable_wis = wis_pie\n",
    "                ours_rel_param = param\n",
    "                print(\"reliable wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)\n",
    "                \n",
    "            model, score, _ = get_model_from_param(param, False)\n",
    "            wis_pie = score[0]\n",
    "            ESS = score[-1]\n",
    "            if wis_pie > best_wis:\n",
    "                best_wis = wis_pie\n",
    "                ours_param = param\n",
    "                print(\"best wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New wis: 97.38  ESS: 91.02 at: {'name': 'POELA-IS', 'env': 'mimic45', 'thresh': 0.4, 'var_coeff': 0.0, 'seed': 1}\n",
      "New wis: 97.87  ESS: 77.87 at: {'name': 'POELA-IS', 'env': 'mimic45', 'thresh': 0.4, 'var_coeff': 0.0, 'seed': 3}\n",
      "New best reliable wis: 93.6  ESS: 376.5 at: {'name': 'POELA-IS', 'env': 'mimic45', 'thresh': 0.6, 'var_coeff': 0.0, 'seed': 1}\n",
      "New best reliable wis: 93.87  ESS: 374.5 at: {'name': 'POELA-IS', 'env': 'mimic45', 'thresh': 0.6, 'var_coeff': 0.0, 'seed': 4}\n"
     ]
    }
   ],
   "source": [
    "# Our method no var\n",
    "var_coeff_list = [0.0]\n",
    "threshold_list = [0.4, 0.6, 0.8, 1.0]\n",
    "\n",
    "best_wis = -100\n",
    "best_reliable_wis = -100\n",
    "for threshold in threshold_list:\n",
    "    for var_coeff in var_coeff_list:\n",
    "        for seed in seeds:\n",
    "            param = {\"name\": \"POELA-IS\", \"env\": \"mimic45\", \"thresh\": threshold, \"var_coeff\": var_coeff, \"seed\": seed}\n",
    "            model, score, _ = get_model_from_param(param)\n",
    "        \n",
    "            wis_pie = score[0]\n",
    "            ESS = score[-1]\n",
    "        \n",
    "            if ESS >= ESS_lb and wis_pie > best_reliable_wis:\n",
    "                best_reliable_wis = wis_pie\n",
    "                ours_novar_rel_param = param\n",
    "                print(\"New best reliable wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)\n",
    "            model, score, _ = get_model_from_param(param, False)\n",
    "            wis_pie = score[0]\n",
    "            ESS = score[-1]\n",
    "            if wis_pie > best_wis:\n",
    "                best_wis = wis_pie\n",
    "                ours_novar_param = param\n",
    "                print(\"New wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New wis: 99.37  ESS: 14.23 at: {'name': 'POELA-kNN', 'env': 'mimic45', 'thresh': 0.01, 'var_coeff': 0.0, 'seed': 1}\n",
      "New wis: 99.98  ESS: 11.81 at: {'name': 'POELA-kNN', 'env': 'mimic45', 'thresh': 0.01, 'var_coeff': 0.0, 'seed': 3}\n",
      "New best reliable wis: 88.81  ESS: 744.85 at: {'name': 'POELA-kNN', 'env': 'mimic45', 'thresh': 0.01, 'var_coeff': 10.0, 'seed': 1}\n",
      "New best reliable wis: 88.82  ESS: 771.94 at: {'name': 'POELA-kNN', 'env': 'mimic45', 'thresh': 0.01, 'var_coeff': 10.0, 'seed': 4}\n",
      "New best reliable wis: 89.94  ESS: 402.48 at: {'name': 'POELA-kNN', 'env': 'mimic45', 'thresh': 0.01, 'var_coeff': 10.0, 'seed': 5}\n",
      "New best reliable wis: 90.31  ESS: 227.4 at: {'name': 'POELA-kNN', 'env': 'mimic45', 'thresh': 0.02, 'var_coeff': 0.1, 'seed': 1}\n",
      "New best reliable wis: 90.82  ESS: 272.41 at: {'name': 'POELA-kNN', 'env': 'mimic45', 'thresh': 0.02, 'var_coeff': 10.0, 'seed': 4}\n",
      "New best reliable wis: 92.21  ESS: 378.09 at: {'name': 'POELA-kNN', 'env': 'mimic45', 'thresh': 0.02, 'var_coeff': 10.0, 'seed': 5}\n"
     ]
    }
   ],
   "source": [
    "# kNN\n",
    "var_coeff_list = [0.0, 0.1, 1.0, 10.0]\n",
    "threshold_list = [0.01, 0.02, 0.05, 0.1]\n",
    "\n",
    "best_wis = -100\n",
    "best_reliable_wis = -100\n",
    "for threshold in threshold_list:\n",
    "    for var_coeff in var_coeff_list:\n",
    "        for seed in seeds:\n",
    "            try:\n",
    "                param = {\"name\": \"POELA-kNN\", \"env\": \"mimic45\", \"thresh\": threshold, \"var_coeff\": var_coeff, \"seed\": seed}\n",
    "                model, score, _ = get_model_from_param(param)\n",
    "            except:\n",
    "                continue\n",
    "            wis_pie = score[0]\n",
    "            ESS = score[-1]\n",
    "        \n",
    "            if ESS >= ESS_lb and wis_pie > best_reliable_wis:\n",
    "                best_reliable_wis = wis_pie\n",
    "                knn_rel_param = param\n",
    "                print(\"New best reliable wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)\n",
    "            \n",
    "            model, score, _ = get_model_from_param(param, False)\n",
    "            wis_pie = score[0]\n",
    "            #print(wis_pie)\n",
    "            ESS = score[-1]\n",
    "            if wis_pie > best_wis:\n",
    "                best_wis = wis_pie\n",
    "                knn_param = param\n",
    "                print(\"New wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New wis: 99.37  ESS: 14.23 at: {'name': 'norm-POEM', 'env': 'mimic45', 'thresh': 0.01, 'var_coeff': 0.0, 'seed': 1}\n",
      "New wis: 99.98  ESS: 11.81 at: {'name': 'norm-POEM', 'env': 'mimic45', 'thresh': 0.01, 'var_coeff': 0.0, 'seed': 3}\n",
      "New best reliable wis: 88.81  ESS: 744.85 at: {'name': 'norm-POEM', 'env': 'mimic45', 'thresh': 0.01, 'var_coeff': 10.0, 'seed': 1}\n",
      "New best reliable wis: 88.82  ESS: 771.94 at: {'name': 'norm-POEM', 'env': 'mimic45', 'thresh': 0.01, 'var_coeff': 10.0, 'seed': 4}\n",
      "New best reliable wis: 89.94  ESS: 402.48 at: {'name': 'norm-POEM', 'env': 'mimic45', 'thresh': 0.01, 'var_coeff': 10.0, 'seed': 5}\n"
     ]
    }
   ],
   "source": [
    "# pg\n",
    "var_coeff_list = [0.0, 0.1, 1.0, 10.0]\n",
    "threshold_list = [0.01]\n",
    "\n",
    "best_wis = -100\n",
    "best_reliable_wis = -100\n",
    "for threshold in threshold_list:\n",
    "    for var_coeff in var_coeff_list:\n",
    "        for seed in seeds:\n",
    "            param = {\"name\": \"norm-POEM\", \"env\": \"mimic45\", \"thresh\": threshold, \"var_coeff\": var_coeff, \"seed\": seed}\n",
    "            model, score, _ = get_model_from_param(param)\n",
    "            wis_pie = score[0]\n",
    "            ESS = score[-1]\n",
    "        \n",
    "            if ESS >= ESS_lb and wis_pie > best_reliable_wis:\n",
    "                best_reliable_wis = wis_pie\n",
    "                npoem_rel_param = param\n",
    "                print(\"New best reliable wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)\n",
    "            model, score, _ = get_model_from_param(param, False)\n",
    "            wis_pie = score[0]\n",
    "            #print(wis_pie)\n",
    "            ESS = score[-1]\n",
    "            if wis_pie > best_wis:\n",
    "                best_wis = wis_pie\n",
    "                npoem_param = param\n",
    "                print(\"New wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New wis: 100.0  ESS: 1.0 at: {'name': 'BCQ', 'env': 'mimic45', 'thresh': 0.0, 'seed': 1}\n",
      "New wis: 100.0  ESS: 1.0 at: {'name': 'BCQ', 'env': 'mimic45', 'thresh': 0.0, 'seed': 2}\n",
      "New best reliable wis: 90.1  ESS: 235.56 at: {'name': 'BCQ', 'env': 'mimic45', 'thresh': 0.01, 'seed': 1}\n",
      "New best reliable wis: 90.55  ESS: 266.1 at: {'name': 'BCQ', 'env': 'mimic45', 'thresh': 0.01, 'seed': 3}\n",
      "New best reliable wis: 91.24  ESS: 226.76 at: {'name': 'BCQ', 'env': 'mimic45', 'thresh': 0.01, 'seed': 4}\n"
     ]
    }
   ],
   "source": [
    "# BCQ\n",
    "threshold_list = [0.0, 0.01, 0.05, 0.1, 0.3]\n",
    "\n",
    "best_wis = -100\n",
    "best_reliable_wis = -100\n",
    "for threshold in threshold_list:\n",
    "    for seed in seeds:\n",
    "        param = {\"name\": \"BCQ\", \"env\": \"mimic45\", \"thresh\": threshold, \"seed\": seed}\n",
    "        model, score, _ = get_model_from_param(param)\n",
    "        wis_pie = score[0]\n",
    "        ESS = score[-1]\n",
    "        \n",
    "        if ESS >= ESS_lb and wis_pie > best_reliable_wis:\n",
    "            best_reliable_wis = wis_pie\n",
    "            bcq_rel_param = param\n",
    "            print(\"New best reliable wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)\n",
    "        \n",
    "        model, score, _ = get_model_from_param(param, False)\n",
    "        wis_pie = score[0]\n",
    "        ESS = score[-1]\n",
    "        if wis_pie > best_wis:\n",
    "            best_wis = wis_pie\n",
    "            bcq_param = param\n",
    "            print(\"New wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New wis: 100.0  ESS: 1.0 at: {'name': 'PQL', 'env': 'mimic45', 'thresh': 0.0, 'seed': 1}\n",
      "New best reliable wis: 89.77  ESS: 243.64 at: {'name': 'PQL', 'env': 'mimic45', 'thresh': 0.01, 'seed': 1}\n",
      "New best reliable wis: 91.18  ESS: 234.14 at: {'name': 'PQL', 'env': 'mimic45', 'thresh': 0.01, 'seed': 4}\n"
     ]
    }
   ],
   "source": [
    "# PQL\n",
    "threshold_list = [0.0, 0.01, 0.05, 0.1, 0.3]\n",
    "\n",
    "best_wis = -100\n",
    "best_reliable_wis = -100\n",
    "for threshold in threshold_list:\n",
    "    for seed in seeds:\n",
    "        param = {\"name\": \"PQL\", \"env\": \"mimic45\", \"thresh\": threshold, \"seed\": seed}\n",
    "        model, score, _ = get_model_from_param(param)\n",
    "        wis_pie = score[0]\n",
    "        ESS = score[-1]\n",
    "        \n",
    "        if ESS >= ESS_lb and wis_pie > best_reliable_wis:\n",
    "            best_reliable_wis = wis_pie\n",
    "            pql_rel_param = param\n",
    "            print(\"New best reliable wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)\n",
    "        \n",
    "        model, score, _ = get_model_from_param(param, False)\n",
    "        wis_pie = score[0]\n",
    "        ESS = score[-1]\n",
    "        if wis_pie > best_wis:\n",
    "            best_wis = wis_pie\n",
    "            pql_param = param\n",
    "            print(\"New wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute WIS and BCa confidence interval on test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accel: -0.018, bz: -0.033, a1: 0.039, a2: 0.938\n",
      "WIS: 87.63949777272308\n",
      "ESS: 528.1810651990509\n",
      "[85.06 90.  ]\n",
      "Accel: -0.017, bz: -0.048, a1: 0.037, a2: 0.934\n",
      "WIS: 86.2944708138815\n",
      "ESS: 491.71875054933196\n",
      "[83.51 88.59]\n",
      "Accel: -0.070, bz: -0.116, a1: 0.017, a2: 0.897\n",
      "WIS: 82.36211922281583\n",
      "ESS: 21.239846409705695\n",
      "[62.65 92.78]\n",
      "Accel: -0.070, bz: -0.050, a1: 0.024, a2: 0.917\n",
      "WIS: 82.36211922281583\n",
      "ESS: 21.239846409705695\n",
      "[64.34 93.6 ]\n",
      "Accel: -0.078, bz: -0.068, a1: 0.020, a2: 0.909\n",
      "WIS: 83.78235170786718\n",
      "ESS: 9.505121394305865\n",
      "[52.23 99.99]\n",
      "Accel: 0.164, bz: 0.108, a1: 0.132, a2: 0.995\n",
      "WIS: 96.32784699205546\n",
      "ESS: 1.2731855860578365\n",
      "[ 62.41 100.  ]\n"
     ]
    }
   ],
   "source": [
    "CIs = []\n",
    "ESSs = []\n",
    "wis_tests = []\n",
    "\n",
    "for param in [ours_param, ours_novar_param, knn_param, npoem_param, bcq_param, pql_param]:\n",
    "    model, val_score, _ = get_model_from_param(param, False)\n",
    "    pibs, pies, rewards, length = compute_is_weights(test_dataset, model)\n",
    "    length = length.astype(int)\n",
    "    wis_pie, weights = wis_ope(pibs, pies, rewards, length)\n",
    "\n",
    "    v_pib = rewards.sum(axis=-1).mean(axis=0)\n",
    "    ESS = 1/(np.sum(weights**2,axis=0))\n",
    "    lb_ub, wis_list = bca_bootstrap(pibs, pies, rewards, length, alpha=0.05)\n",
    "\n",
    "    wis_tests.append(wis_pie)\n",
    "    ESSs.append(ESS)\n",
    "    CIs.append(lb_ub)\n",
    "    print(\"WIS:\", wis_pie)\n",
    "    print(\"ESS:\", ESS)\n",
    "    print(lb_ub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "ESSs_noess = np.array(ESSs)\n",
    "CIs_noess = np.array(CIs)\n",
    "wis_tests_noess = np.array(wis_tests)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accel: -0.021, bz: 0.048, a1: 0.054, a2: 0.954\n",
      "WIS: 91.46953279566598\n",
      "ESS: 624.9277837733208\n",
      "[89.5919446  93.24882146]\n",
      "Accel: -0.024, bz: -0.005, a1: 0.043, a2: 0.942\n",
      "WIS: 91.89766854852357\n",
      "ESS: 541.6944688343083\n",
      "[89.76281071 93.77737454]\n",
      "Accel: -0.023, bz: -0.048, a1: 0.035, a2: 0.932\n",
      "WIS: 87.95621154928499\n",
      "ESS: 372.00541193367064\n",
      "[84.77195178 90.58275047]\n",
      "Accel: -0.021, bz: -0.005, a1: 0.043, a2: 0.943\n",
      "WIS: 87.71563821179501\n",
      "ESS: 399.56605713789367\n",
      "[84.90612873 90.04370031]\n",
      "Accel: -0.020, bz: 0.070, a1: 0.060, a2: 0.958\n",
      "WIS: 82.67837938055808\n",
      "ESS: 228.8200254461097\n",
      "[78.2510896 86.8322061]\n",
      "Accel: -0.022, bz: 0.035, a1: 0.051, a2: 0.951\n",
      "WIS: 84.40907192144543\n",
      "ESS: 231.93008071238776\n",
      "[80.13960247 88.29341149]\n"
     ]
    }
   ],
   "source": [
    "CIs = []\n",
    "ESSs = []\n",
    "wis_tests = []\n",
    "\n",
    "# for param in [ours_rel_param, ours_novar_rel_param, knn_rel_param, knn_novar_rel_param, \n",
    "#               npoem_rel_param, npoem_novar_param, nomask_rel_param, nomask_novar_param, \n",
    "#               bcq_rel_param, pql_rel_param]:\n",
    "for param in [ours_rel_param, ours_novar_rel_param, knn_rel_param, npoem_rel_param, bcq_rel_param, pql_rel_param]:\n",
    "    model, _, _ = get_model_from_param(param)\n",
    "    pibs, pies, rewards, length = compute_is_weights(test_dataset, model)\n",
    "    length = length.astype(int)\n",
    "    wis_pie, weights = wis_ope(pibs, pies, rewards, length)\n",
    "\n",
    "    v_pib = rewards.sum(axis=-1).mean(axis=0)\n",
    "    ESS = 1/(np.sum(weights**2,axis=0))\n",
    "    lb_ub, wis_list = bca_bootstrap(pibs, pies, rewards, length, alpha=0.05)\n",
    "\n",
    "    wis_tests.append(wis_pie)\n",
    "    ESSs.append(ESS)\n",
    "    CIs.append(lb_ub)\n",
    "    print(\"WIS:\", wis_pie)\n",
    "    print(\"ESS:\", ESS)\n",
    "    print(lb_ub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "444853adc6674f658e595c08a4a13e2b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/html": [
       "<p>Failed to display Jupyter Widget of type <code>HBox</code>.</p>\n",
       "<p>\n",
       "  If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n",
       "  that the widgets JavaScript is still loading. If this message persists, it\n",
       "  likely means that the widgets JavaScript library is either not installed or\n",
       "  not enabled. See the <a href=\"https://ipywidgets.readthedocs.io/en/stable/user_install.html\">Jupyter\n",
       "  Widgets Documentation</a> for setup instructions.\n",
       "</p>\n",
       "<p>\n",
       "  If you're reading this message in another frontend (for example, a static\n",
       "  rendering on GitHub or <a href=\"https://nbviewer.jupyter.org/\">NBViewer</a>),\n",
       "  it may mean that your frontend doesn't currently support widgets.\n",
       "</p>\n"
      ],
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=1000.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Accel: -0.005, bz: -0.020, a1: 0.045, a2: 0.944\n"
     ]
    }
   ],
   "source": [
    "ci_pib, _ = bca_bootstrap_mc(test_dataset[\"rewards\"].sum(axis=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "81.10183639398998 [79.83305509 82.18267816]\n"
     ]
    }
   ],
   "source": [
    "print(v_pib, ci_pib)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "ESSs_3seed = np.array(ESSs)\n",
    "CIs_3seed = np.array(CIs)\n",
    "wis_tests_3seed = np.array(wis_tests)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "ESSs = ESSs_3seed\n",
    "CIs = CIs_3seed\n",
    "wis_tests = wis_tests_3seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.set_printoptions(precision=2, linewidth=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Weights on low reward samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WIS: 97.18793922867032\n",
      "ESS: 2936.7939155941785\n",
      "WIS: 97.87511531583822\n",
      "ESS: 2722.5502164135955\n",
      "WIS: 99.81041850227159\n",
      "ESS: 1054.8303740271786\n",
      "WIS: 99.65193736473071\n",
      "ESS: 1102.3197481723705\n",
      "WIS: 85.16059005417785\n",
      "ESS: 589.0570844000422\n",
      "WIS: 85.08755549123165\n",
      "ESS: 610.4690433047814\n"
     ]
    }
   ],
   "source": [
    "CIs = []\n",
    "ESSs = []\n",
    "wis_tests = []\n",
    "\n",
    "# for param in [ours_rel_param, ours_novar_rel_param, knn_rel_param, knn_novar_rel_param, \n",
    "#               npoem_rel_param, npoem_novar_param, nomask_rel_param, nomask_novar_param, \n",
    "#               bcq_rel_param, pql_rel_param]:\n",
    "for param in [ours_rel_param, ours_novar_rel_param, knn_rel_param, npoem_rel_param, bcq_rel_param, pql_rel_param]:\n",
    "    model, _, _ = get_model_from_param(param)\n",
    "    pibs, pies, rewards, length = compute_is_weights(train_dataset, model)\n",
    "    length = length.astype(int)\n",
    "    wis_pie, weights = wis_ope(pibs, pies, rewards, length)\n",
    "\n",
    "    v_pib = rewards.sum(axis=-1).mean(axis=0)\n",
    "    ESS = 1/(np.sum(weights**2,axis=0))\n",
    "    #lb_ub, wis_list = bca_bootstrap(pibs, pies, rewards, length, alpha=0.05)\n",
    "\n",
    "    wis_tests.append(wis_pie)\n",
    "    ESSs.append(ESS)\n",
    "    #CIs.append(lb_ub)\n",
    "    print(\"WIS:\", wis_pie)\n",
    "    print(\"ESS:\", ESS)\n",
    "    #print(lb_ub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.02812060771329683,\n",
       " 0.021248846841617874,\n",
       " 0.0018958149772840827,\n",
       " 0.0034806263526928616,\n",
       " 0.14839409945822157,\n",
       " 0.1491244450876834]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# because low reward means reward = 0\n",
    "# reward in dataset is either 0 or 100, and weighted IS is a normalized weighted average.\n",
    "[1.0 - x/100.0 for x in wis_tests]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Show the effect of different hyperparameter \\delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 2, 3, 4, 5]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seeds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "best wis: 97.38  ESS: 91.02 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.4, 'var_coeff': 0.0, 'seed': 1}\n",
      "best wis: 97.87  ESS: 77.87 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.4, 'var_coeff': 0.0, 'seed': 3}\n",
      "best wis: 98.5  ESS: 85.41 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.4, 'var_coeff': 0.1, 'seed': 3}\n",
      "best wis: 99.25  ESS: 88.19 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.4, 'var_coeff': 1.0, 'seed': 3}\n",
      "reliable wis: 93.6  ESS: 376.5 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.6, 'var_coeff': 0.0, 'seed': 1}\n",
      "best wis: 93.6  ESS: 376.5 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.6, 'var_coeff': 0.0, 'seed': 1}\n",
      "reliable wis: 93.87  ESS: 374.5 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.6, 'var_coeff': 0.0, 'seed': 4}\n",
      "best wis: 93.87  ESS: 374.5 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.6, 'var_coeff': 0.0, 'seed': 4}\n",
      "reliable wis: 94.01  ESS: 407.28 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.6, 'var_coeff': 0.1, 'seed': 2}\n",
      "best wis: 94.01  ESS: 407.28 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.6, 'var_coeff': 0.1, 'seed': 2}\n",
      "reliable wis: 91.9  ESS: 521.49 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.8, 'var_coeff': 0.0, 'seed': 1}\n",
      "best wis: 91.9  ESS: 521.49 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.8, 'var_coeff': 0.0, 'seed': 1}\n",
      "reliable wis: 92.48  ESS: 456.04 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.8, 'var_coeff': 0.0, 'seed': 2}\n",
      "best wis: 92.48  ESS: 456.04 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.8, 'var_coeff': 0.0, 'seed': 2}\n",
      "reliable wis: 93.0  ESS: 560.36 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.8, 'var_coeff': 0.1, 'seed': 2}\n",
      "best wis: 93.0  ESS: 560.36 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.8, 'var_coeff': 0.1, 'seed': 2}\n",
      "reliable wis: 93.06  ESS: 558.67 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.8, 'var_coeff': 1.0, 'seed': 2}\n",
      "best wis: 93.06  ESS: 558.67 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 0.8, 'var_coeff': 1.0, 'seed': 2}\n",
      "reliable wis: 89.37  ESS: 365.75 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 1.0, 'var_coeff': 0.0, 'seed': 1}\n",
      "best wis: 89.37  ESS: 365.75 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 1.0, 'var_coeff': 0.0, 'seed': 1}\n",
      "reliable wis: 91.82  ESS: 432.74 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 1.0, 'var_coeff': 0.0, 'seed': 2}\n",
      "best wis: 91.82  ESS: 432.74 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 1.0, 'var_coeff': 0.0, 'seed': 2}\n",
      "reliable wis: 92.66  ESS: 568.62 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 1.0, 'var_coeff': 0.0, 'seed': 3}\n",
      "best wis: 92.66  ESS: 568.62 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 1.0, 'var_coeff': 0.0, 'seed': 3}\n",
      "reliable wis: 92.81  ESS: 545.87 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 1.0, 'var_coeff': 0.1, 'seed': 3}\n",
      "best wis: 92.81  ESS: 545.87 at: {'name': 'POELA', 'env': 'mimic45', 'thresh': 1.0, 'var_coeff': 0.1, 'seed': 3}\n"
     ]
    }
   ],
   "source": [
    "# Our method\n",
    "var_coeff_list = [0.0, 0.1, 1.0, 10.0]\n",
    "threshold_list = [0.4, 0.6, 0.8, 1.0]\n",
    "\n",
    "ablation_models = []\n",
    "for threshold in threshold_list:\n",
    "    best_wis = -100\n",
    "    best_reliable_wis = -100\n",
    "    ours_rel_param = None\n",
    "    ours_param = None\n",
    "    for var_coeff in var_coeff_list:\n",
    "        for seed in seeds:\n",
    "            param = {\"name\": \"POELA\", \"env\": \"mimic45\", \"thresh\": threshold, \"var_coeff\": var_coeff, \"seed\": seed}\n",
    "            model, score, _ = get_model_from_param(param)\n",
    "        \n",
    "            wis_pie = score[0]\n",
    "            ESS = score[-1]\n",
    "        \n",
    "            if ESS >= ESS_lb and wis_pie > best_reliable_wis:\n",
    "                best_reliable_wis = wis_pie\n",
    "                ours_rel_param = param\n",
    "                print(\"reliable wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)\n",
    "                \n",
    "            model, score, _ = get_model_from_param(param, False)\n",
    "            wis_pie = score[0]\n",
    "            ESS = score[-1]\n",
    "            if wis_pie > best_wis:\n",
    "                best_wis = wis_pie\n",
    "                ours_param = param\n",
    "                print(\"best wis:\", round(wis_pie, 2), \" ESS:\", round(ESS, 2), \"at:\", param)\n",
    "    if ours_rel_param:\n",
    "        ablation_models.append(ours_rel_param)\n",
    "    else:\n",
    "        ablation_models.append(ours_param)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WIS: 86.62125725983307\n",
      "ESS: 1278.0850529000188\n",
      "WIS: 90.07865228412467\n",
      "ESS: 819.646050816067\n",
      "WIS: 90.87345341488546\n",
      "ESS: 706.1999061273318\n",
      "WIS: 91.23976246994216\n",
      "ESS: 642.5375063111641\n"
     ]
    }
   ],
   "source": [
    "CIs = []\n",
    "ESSs = []\n",
    "wis_tests = []\n",
    "\n",
    "# for param in [ours_rel_param, ours_novar_rel_param, knn_rel_param, knn_novar_rel_param, \n",
    "#               npoem_rel_param, npoem_novar_param, nomask_rel_param, nomask_novar_param, \n",
    "#               bcq_rel_param, pql_rel_param]:\n",
    "for param in ablation_models:\n",
    "    model, _, _ = get_model_from_param(param)\n",
    "    pibs, pies, rewards, length = compute_is_weights(test_dataset, model)\n",
    "    length = length.astype(int)\n",
    "    wis_pie, weights = wis_ope(pibs, pies, rewards, length)\n",
    "\n",
    "    v_pib = rewards.sum(axis=-1).mean(axis=0)\n",
    "    ESS = 1/(np.sum(weights**2,axis=0))\n",
    "    #lb_ub, wis_list = bca_bootstrap(pibs, pies, rewards, length, alpha=0.05)\n",
    "\n",
    "    wis_tests.append(wis_pie)\n",
    "    ESSs.append(ESS)\n",
    "    #CIs.append(lb_ub)\n",
    "    print(\"WIS:\", wis_pie)\n",
    "    print(\"ESS:\", ESS)\n",
    "    #print(lb_ub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1278.0850529000188, 819.646050816067, 706.1999061273318, 642.5375063111641]"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ESSs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[86.62125725983307, 90.07865228412467, 90.87345341488546, 91.23976246994216]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "wis_tests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "86.62125725983307, 90.07865228412467, 91.46345341488546, 91.23976246994216"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WIS: 91.61999140054891\n",
      "ESS: 3601.118917637893\n",
      "WIS: 98.41025151771905\n",
      "ESS: 2242.0739140638734\n",
      "WIS: 98.89883450479996\n",
      "ESS: 1993.081707432032\n",
      "WIS: 99.11940231493642\n",
      "ESS: 1769.4637894867988\n"
     ]
    }
   ],
   "source": [
    "CIs = []\n",
    "ESSs = []\n",
    "wis_tests = []\n",
    "\n",
    "# for param in [ours_rel_param, ours_novar_rel_param, knn_rel_param, knn_novar_rel_param, \n",
    "#               npoem_rel_param, npoem_novar_param, nomask_rel_param, nomask_novar_param, \n",
    "#               bcq_rel_param, pql_rel_param]:\n",
    "for param in ablation_models:\n",
    "    model, _, _ = get_model_from_param(param)\n",
    "    pibs, pies, rewards, length = compute_is_weights(train_dataset, model)\n",
    "    length = length.astype(int)\n",
    "    wis_pie, weights = wis_ope(pibs, pies, rewards, length)\n",
    "\n",
    "    v_pib = rewards.sum(axis=-1).mean(axis=0)\n",
    "    ESS = 1/(np.sum(weights**2,axis=0))\n",
    "    #lb_ub, wis_list = bca_bootstrap(pibs, pies, rewards, length, alpha=0.05)\n",
    "\n",
    "    wis_tests.append(wis_pie)\n",
    "    ESSs.append(ESS)\n",
    "    #CIs.append(lb_ub)\n",
    "    print(\"WIS:\", wis_pie)\n",
    "    print(\"ESS:\", ESS)\n",
    "    #print(lb_ub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[91.62, 98.41, 98.9, 99.12]\n"
     ]
    }
   ],
   "source": [
    "print([round(x, 2) for x in wis_tests])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3601.12, 2242.07, 1993.08, 1769.46]\n"
     ]
    }
   ],
   "source": [
    "print([round(x, 2) for x in ESSs])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "def toprint(name):\n",
    "    if name == \"POELA-kNN\":\n",
    "        return \"PO-mu\"\n",
    "    elif name == \"POELA-kNN-IS\":\n",
    "        return \"PO-mu-NVar\"\n",
    "    if name == \"norm-POEM\":\n",
    "        return \"PO-CRM\"\n",
    "    elif name == \"norm-POEM-IS\":\n",
    "        return \"PO-CRM-NVar\"\n",
    "    else:\n",
    "        return name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "BIGGER_SIZE = 20\n",
    "plt.rc('font', size=BIGGER_SIZE)  \n",
    "plt.rc('axes', titlesize=BIGGER_SIZE)  \n",
    "plt.rc('axes', labelsize=BIGGER_SIZE)  \n",
    "plt.rc('xtick', labelsize=BIGGER_SIZE)  \n",
    "plt.rc('ytick', labelsize=BIGGER_SIZE)  \n",
    "plt.rc('legend', fontsize=BIGGER_SIZE)  \n",
    "plt.rc('figure', titlesize=BIGGER_SIZE)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hlines(v_pib, 0, 5, color=\"grey\",label=\"Behavior Policy\")\n",
    "plt.hlines(ci_pib[0], 0, 5, linestyle='--', color=\"grey\")\n",
    "plt.hlines(ci_pib[1], 0, 5, linestyle='--', color=\"grey\")\n",
    "# plt.fill_between([0, len(results)-1], ci_pib[1], ci_pib[0],\n",
    "#                  color='C0', alpha=0.2, rasterized=True)\n",
    "plt.legend([\"Behavior Policy\"])\n",
    "\n",
    "errors = np.abs(CIs.transpose() - wis_tests)\n",
    "for i, alg_idx in enumerate([0,1,2,3,4,5]):\n",
    "    plt.errorbar(i, wis_tests[alg_idx], yerr=errors[:,alg_idx:alg_idx+1], fmt = 'o', capsize=4, color = 'k') \n",
    "\n",
    "\n",
    "plt.xticks(range(6),[\"POELA\", \"no CRM\", \"PO-mu\", \"PO-CRM\", \"PQL\", \"BCQ\"],rotation=20)\n",
    "#plt.xticks(range(len(ESSs)),[str(round(float(e), 1)) for e in ESSs])\n",
    "#plt.xlabel(\"Effective Sample Size\")\n",
    "\n",
    "#plt.ylim([80,100])\n",
    "plt.ylabel(\"95% BCa CI\")\n",
    "#plt.xlabel(\"Different algorithms\")\n",
    "#plt.title(\"Confidence Interval of WIS Estimates on Test Set\")\n",
    "#plt.savefig(\"plots/mimic_ci.png\", bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trade-offs in hyperparameter and random seeds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdoAAAFeCAYAAAAixd5hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOydaZgU5dWw7zMbDAyyKAISFRKVIUaEOJiI36tgEo1LlBgTSMAteWNckrwaHcTExCUmKmTRJEJiFo1CBFeiuCeAC25g2FRATdToIDuDDAzM0uf7UVVNdU91d1VPd8927uvqq6mqp6rPPN3UqXOes4iqYhiGYRhGfihqawEMwzAMozNjitYwDMMw8ogpWsMwDMPII6ZoDcMwDCOPmKI1DMMwjDxS0tYCdDT2228/HTJkSKuvs3PnTnr27Nl6gboINl/RsTmLjs1ZdMLO2auvvrpZVfsXQKR2hynaiAwZMoSlS5e2+jqLFi1i7NixrReoi2DzFR2bs+jYnEUn7JyJyHv5l6Z9Yq5jwzAMw8gjpmgNwzAMI4+YojUMwzCMPBJa0YrIT0TkuAxj/kdEftJ6sQzDMAyjcxDFor0WGJthzHHANdkKYxiGYRidjVy7jkuAWI6vaRiGYRgdllwr2qOAzVFOEJGzROS3IvKciHwkIioiszKcM0ZEHhORrSKyS0RWisilIlKc5pzTRGSRiGwXkToReVlEzo0iq2EYhmFEJW0erYgsSNp1noiMDRhaDBwIHAzcE1GGq4EjgTrgA6Ayg0xnAA8Au4G5wFbgS8CvgWOBrwac813gt8AWYBbQAJwF3CkiR6jqFRFlzgmqioik3DYMwzA6PpkKVoz1/VuBIe4rmRiOEpsLXBZRhstwFOzbwPHAwlQDRWQf4I9AMzBWVZe6+38MLADOEpGJqjrHd84Q4Bc4CrlKVd91918PLAEuF5EHVPXFiHJH5rYLneeWgR88S8+T+/PgL59h/eD/AZRLBp6JAA0l3Si7eqMpXcMwjE5CWtexqhZ5L0CAa/37fK8SVR2gqt9Q1U1RBFDVhar6lobrQH8W0B+Y4ylZ9xq7cSxjgIuSzvkm0A34nadk3XO2AT93Ny+MInPrUCrfngtA5VtzcJ5fmuNHy5r2wLW9kev60HDD/oUTyzAMw8gLUdZozwfm5UuQkJzgvj8RcOxZYBcwRkS6hTzn8aQxeeV/vn4IAz94Fs9OFWBgzfMc12tm4Piypj2Ee/4wDMMw2ivSnm7k7vrvQmC2qk4OOL4EqMJxAb8acPw14HDgk6q62t23CdgP2E9VtwScUwf0BHqq6q4Ucl0AXAAwYMCAo+bMmRM0LBS7X38dgIb996ds40YAuvdrTH3CoJFZf1Znoq6ujoqKirYWo0NhcxYdm7PohJ2zcePGvaqqVQUQqd0RqamAiBQBlwCTgOE4yqnEPTYK+DZwi6q+mWtBXXq779tTHPf294l4Tk93XKCiVdXbgdsBqqqqNNui41v+dg8bfvs7BHjve9/l4N/+DgUGHrWVfofuDjxHJ9baWi1W7D0bbM6iY3MWHZuzzESpDFUGPA3cAnwC2AH4NcA7OOuhk3IpYEQ8eaKY6dmckxUbf/pTJOljBGXDq31TnmNK1jAMo2MTZY22GhgHXAcMAP7kP6iqtTjrpCflTLqWeFZp7xTH90kaF+Wcj1ohVzhUSXw2IWA7+ZT249o3DMMwohPFdTwJWKyq1wOISJAGeAcnpzVfrMVZoz0MSFijFZESYCjQBPwn6Zz93HNeTDpnEI7b+INU67O5pHL1G8h1jld7Q79Ghk9cl/GcmStmcvHIi7P6vHnLapj+5FrW1dZzQJ9yqk8axvhRg7O6lmEYhpEdUSzaocBLGcZsBfplL05GvAIaXww4dhzQA3hBVfeEPOfkpDF5RURoKCmLdM5Hez7Kyqqdt6yGqx5cRU1tPQrU1NZz1YOrmLesJvK1DMMwjOyJomjrSQwyCuIgoDZ7cTJyP06Jx4kiEo9eE5HuwA3uZnKuzB3AHuC7bvEK75y+wA/dzd/nSd4WlP5oY+ixO4uK6FXWK6t12ulPrqW+sTlhX31jM9OfXBv5WoZhGEb2RHEdLwdOFJEyVW1IPigivXHWZ1+IIoCIjAfGu5sD3fdjRORO99+bvRKJqvqRiHwbR+EuEpE5OFb06cAwd/9c//VV9R0RqQZ+AywVkbnsLcH4MeCXhagK5eFYtd0yjjti6EEATG6si1QlynMX19TWBx5fl2K/YRiGkR+iKNo/ArOB2SLyLf8BEemDYzn2Jbp1OBJILu7/cfcF8B4Qr0WsqvNE5HjgR8BXgO445Rt/APwmqMKUqv5WRN51r3MOjiX/BnC1qv41orytpuzqjbBokbN6nIbKvpVUV1VHUrJXPbiqhSXr54A+5REkNQzDMFpLaEWrqveIyOdxKkSdDmwDEJGlOEUiugG3qepjUQRQ1Wtxet1GOWcxcErEcx4BHolyTr5pKOnmlFxMos5VrGu2rWH60ulMGT0llLINchf7KS8tpvqkYdkLbBiGYUQmUsEKVf2WiDwH/B8wAic35dPA68CvVPWO3IvYefE3D7ht2W0sen8Ra7atiR+v7FtJRWlFaIs2nVu4b49SrvnS4RZ1bBiGUWAi96NV1TtVdRRQgbPG2UtVjzAlmx0igqpS11jHmm1rmDx8MivPWcnk4ZNZs20Nde4abRjSuYVrdzVy6dzlHHvTAos8NgzDKCCRLFo/qlqPE4lstBIRoVdZLyYPnxx3E08ZPQUgUtRx9UnDUq7ReqraS/MBzLo1DMMoAKEVrZsOMwj4tz9PVUTOx4ka3olT5/iVnEvZBbh45MUJ0cWeso2S2uMpzkvnLk87zkvzMUVrGIaRf6K4jn8OvOw/R0S+h1OK8UvARJyUm0/mVMIuRLJSzSZ/dvyowQwOEVlsaT6GYRiFIYqiPRb4p+sy9rgCqMGpyvQ1d98PciSbkSXVJw2jvLQ47RhL8zEMwygMUdZoBwP/9DZcy/VA4EpVfd7d91UcpWu0EVfPW8U9L79Pc5oAKgFL8zEMwygQUSzacsDfNPVYnBibf/j2/RtHIRttwNXzVjHrpf+mVbLgfGm2PmsYhlEYoijaGqDSt30STmu5Fb59fbFI5DbjnpffDzUuzBquYRiGkRuiuI4XAueKyHdxLNvTgQdUNeYbcwgQ7m5v5JxMlqyHuY0NwzAKRxSL9kagDrgVuB1H2V7rHRSR/YHjidhUwMgdxSGilHuUFpnb2DAMo4CEVrSq+g5OTeP/A74PfEpV/WXxDwZuA+7MpYBGeL7+mQMzjmlsVqsMZRiGUUCi1jpeD/wuxbElwJJcCGWEZ96yGq575HW27WoEoLQIGmOpxzfG1IpVGIZhFJCsSzAabc+8ZTVU37+Cxua9a7PplKyHFaswDMMoHJGbChjth+seeT1ByYbFilUYhmEUDlO0HZR5y2ri7uJU9O1R2mKf9aQ1DMMoLKZoOyjTn1yb9nixCMt+ciK3TBjJ4D7lCE7+7I1nHmHrs4ZhGAXE1mg7KJnWWZtVOfamBVSfNIzFU08okFSGYRhGMqZoOxDzltUw/cm1rKutp0gkY4EK6z1rGIbR9pjruIMwb1kNVz24ipraepTwVaC83rOGYRhG2xDZohWRYmAYTl3jwF5sqvpsK+Uykpj+5FrqG5uzOtfSeQzDMNqOSIpWRH4MXAb0zjA0fTNUIzKtUZaWzmMYhtF2hFa0IjIFuA7YDtyN0zygKU9yGUkc0KecmiyUraXzGIZhtC1RLNpv47TK+7SqbsqTPEYKqk8axlUPrkrrPh7cp5xxlf1ZuGYT62rrOaBPOdUnDbNAKMMwjDYkiqI9EPijKdm2wVOWl9+7IjAQanCfckvjMQzDaIdEiTregKUDtSnjRw3ml187kvLSxCVwcw8bhmG0X6IoznuBL4tIN1Xdky+BuhKqivh6yCZvB+FZtl4+rbmHDcMw2jdRFO1PgM8C94vI993+tEaWbKrfxLQl05gyegoigqoybck0epX14uKRF6c9d/yowQmKdd6yGo69aYEpXsMwjHZIFEX7OlAKHACcIiLbgdqAcaqqn8iFcJ0VVaVZm5m1ehYAU0ZPYdqSacxaPYvJwyeHsmw9vEIWXpCUVYMyDMNoX0RZoy3CSef5r/vaDkjAK+/VpsThmyLykojsEJFdIrJMRL7vFtTwjx0iIprmNSff8gbIz8AeA5k8fDKzVs9ixF0j4krWs3DDElTIwqpBGYZhtB9CW7SqOiSPckTlr8DZwEZgLrAT+DxwK3CciHxVtUVo7gpgXsC1XsunoOmYMnpK3Kr1tqMoWUhdyMKqQRmGYbQPOlwUsYiMx1Gy7wBHq+pmd38pTsDWV4BzgTuTTl2uqtcWTtLMTFsyrcV2VGWbqpCFVYMyDMNoH2Tt5hWRfUTkQBHZJ5cCheBM9/2XnpIFUNVG4Mfu5vcKLFMkVJX1u9bH3cUrz1kZdyNPWzKNlsZ4aqpPGmbpPoZhGO2YqLWOi4Fq4H+Bob797wB/An6hqvkuyzjQff9PwDFv36dFpI+q+oO1DhCR7wD7AluAF1V1ZR7lTImIUCzFCWuyU0ZPAaBXWa9IFq2l+xiGYbRvJKz1JCJlwBPA8YACHwAfAoOAj+EEQj0HnKiqDXmR1pHjb8DXgUtUdUbSsU8Bq9zNY1T1JREZguNmDmIRcK6q/jfDZ14AXAAwYMCAo+bMaX38VF1dHRUVFa2+TlfB5is6NmfRsTmLTtg5Gzdu3KuqWlUAkdofqhrqBUwFYsDDwKFJxz6BE2jUDEwNe81sXsA3cBT920A/3/4S4AH3mAInu/v3B64HPg30cV/HAQvccW8BPcN+/lFHHaW5YOHChTm5TiGJxWJpt/NJR5yvtsbmLDo2Z9EJO2fAUs2jbmjPryhrtN/AidAdr6pvJSnrf+Osnb4OTIpwzWyYAzyOo9zfEJHbReQWYDlwCo7iBEfpo6obVfUnqvovVa11X88CJwIvA4fguMKNNMxYPiNh/VjdAhszls/IcKZhGEbXJoqiPQR4XFVjQQfd/Z4CzBvu55wOXAGsx4lA/iaOK/v/4ay/gpP6k+46TTjryuBYuEYKVJUdDTsSgrW8Ahs7GnZECt4yDMPoakQJhmoAMjniewKN2YsTDldJ/tJ9xRGRcmAkUI9jXWfC60TUM6cCdjL8wVqzVs+K5/5mU2DDMAyjqxHFol0JnCUi/YMOish+wFk4hSHairOB7sC96qT7ZOKz7ntQBLPhw69sPUzJGoZhZCaKov0d0B94RUS+JSIfF5FyERkqIufjrHf2d8fllaDcXREZDdwE1OEEP3n7P+NGTCePPwG4zN2clXy80LR396vnLvYTNee3q+I1fRg69VGOvWkB85bVtLVIhmEUkCglGO8VkZE40ce3BwwRYJqq3psr4dLwtIjU4wRn7QAOxwmE2gOcqap+C/Vm4HARWYSzjgswAvC6pP9YVV8ogMxxZiyfQb9d/eLNAzwl5nXu8fZ7JG8XGv+arOcu9rah61q285bVZMxfvnreKma/9F+8xxFr+mAYXY9IBStU9Yci8jDwLWAU0BunucAy4C+q+mLuRQzkfmAiMBkoB9bhBDbdpKrvJo29G/gyMBo4GacD0Qacco2/U9XnCiQzsDewiPq9JRf9Suy2ZbdR11iXVfu8fCEi9CrrlZMCG52BectquO6R19m2a+/qRJACnbesJkHJenhNH0zRGkbXIHKtY1V9CXgpD7JEkWE6MD3k2D8Df86vROHxlNTczXOZuXpmQmBRdVU105dOz0n7vFRkay0nW9re39EVlay/LaGfZAU6/cm1LZSshzV9MIyuQ4drKtAZ8Nrk+Um2FPMR3Ttj+Qx2NOzI2lpO/vyupmQhuC2hH78CTadMi0SYt6ymhVUbxh1tGEbHImUwlIgc5L6Kk7YzvgonfsdE1Wkq4McLLMpXdK/lwuaGTJZonx6lgf9OplmVS+cuZ9T1T8WDozxruaa2HmWvO9qCpwyjY5POon0Xp0ThcOBN33YmNMN1uzRxK7K+V2Bgkec+9pNN+7xkPAWuaIK1PGn4pC7pAs6WVG0JPfzPK2GeXbbtaoyv7QZZy7aeaxgdn3QK8S4cpbk9adtoBV5gUb/yflw0+qIEC7aitCK+RpuP6N7Gnw1gatMepvp3vnMjDU/9irKr0xbSMlyqTxqWco0WYHt9Y+C/01Hf2Mx1j7xO7a7g8baeaxgdm5SKVlXPS7dtZM/FIy9mUe2iwMCiGctn5CW6V1Upa9oTeKysaU9CoFN7Sy9qT3iW5eX3rqA5wGQ9oE95fJ01ylPptl2N9O1RmhDJ7L9mW2DrxYaRG8zF207wFFm66N5CKMDWBkx1BTxlk2zZlpcWM66yf1qLNx2qzjWSr1l90rDWCx2R5Ojqmtp6Lpu7nKXvbeWG8UcUXB7D6MiEVrQi0gxcq6o/TTPmR8B1qmoKvBUERfc23LB/C4tUgIaSbqHcvpkUsqdUvYApyE96UWfBn8Ljt/gyRSWnY3t9I7+eMLLNrEi/BVsk0sJiV2D2S/+l6uB+ZtkaRgSiKERxX2HGGTkkitu3NVjzgGiMHzW4hcK5bO7ylONLi4TGWGqH8gF9ygOvWQiSLdggtzg4ytaCswwjGlFqHYehL7A7x9fs8oSxRjMRiwV2N4zjpfdY84DWkW49taJ7CYPd48mz2VYuYo8olnhNbb2lHBlGBNJatCKS3Kd1SMA+gGLgIJym72tzJJuRI7x118tKugVaxg0l3SjzBUIFNQ8wZRuO6pOGcWkKq7Z2VyPLfnIiED7QKNW4XAcqRY1stnrNhhGeTK7jRexN6VHgXPcVhAAx4PKcSGZkRSwWo6hor6Oiubl577rriZcnlHn0yj6WueOteUA0Uim75DrIHn5rN4yLOCgg6aoHV7H0va088GpNi/0enkxTR8aoDag+FUSm/OBkLL83GhbB3bXJpGivx1GwAvwER/E+EzCuGdgCLFTVNbkU0HCjizMdF+H8J86nrqGOuafNpaioiFgsxsRHJ9KztCeTh0/OuO5qzQPCs1cJNgESV3aqyjVfOjwwIjmqazhVAYt7Xn6/xRqql4u7uzEWP6ehOUb1fSviObrpbvDVJw2j+v4VNDaHT0qy/N5wpHpgAvMIdBXSKlpVvdb7t4icC8xT1d/kWygjERGJP+0ko+7xWCxGXUMda7atYcL8Ccw9bS4T5k9gzbY1VPat5IqjrogrWUhtnVrzgHBMf3Itzb2foFvxbvZsOA0Q6hubuOGlm/jWmE9y45lnMP3JtdTU1lMsErcAIfzNNZUiSxWoFGRFN8Y0vr+mtp7q+1dw7cOvs70+QPFGLEfTVvm9HQ2r+GVE6Uc7NJ+CGKlJZ9GKe7yoqChBuR5595EAVPatZM6pc/jFq79IOC/duqs1D8jMutpdlA3YTVm/xQDs2XAa3QbMp6HnYnY0HMgZow8AaJUlk8qdWxyQehOWxmaltn6v4vWXfwyKiPY+S0jUw+ks9DBu0q7kSk31wGQega5DrqOOjTwQNurYU7Z+PCXrrbuuPGdl3I3sNRcwonNAnx7s2XAaDVuPpazfYnoNv4qyfosp23l8/AEmnSUThuqThlFeWpywr7y0mK9/5sDA/X3KUzcxSIUnTzrr+d2bTuXXE0YyuE85AgzuU86NZx6RMngrU2OErtY8IZXlbx6BrkOkwhLi3NHPAk4CBgPdAoapqn4uB7IZEYnFYkyYPyFh38RHJ3L8x463ddccE695vOG0uFULcPVnp8bnNKwlk8q6S1UUY/yowVQd3K/FfmhZrSoMNbX1DE5hPYsrX9j83jBu0q7mSg2qj93W6VxGYYlSGaob8BgwFuKeJP9dWn37jQLjKVlvTdbvRgbHsrV119wxftRgVJUbXrqJBt/+N5tmo+rMbSrXr9+SyRQok0rBpVN8ngIuLhJKiyVjgFOxCNUnDeOyuctb/OeNWqAizMNFV3OlpntgMroGUVzHVwLjgBuA/jhK9VrgAOAbwPvAHKAstyIamdy73hptRVlFXMl6buTKvpVUlFVQXJzoajQl2zpUlTebZtPQ85mULvlUrl+/JXPtw6+3yr3sJ9kyPqBPOdPPOjLu8k1Fs6rz4JDieBQFGMZNGsWVOm9ZDcfetIChUx/l2JsWdFj38vhRg1k89QTeuelUFk89wZRsFyOKov0q8C9VvUZVt3g7VXW9qs4BTgBOAy7NsYxdHhGhoSTIS+8Um/CU5h1fvCOuZGHvmu0dX7wj0uclK3Zbx21JqlSoycMnx13y40cN5sYzj6CvrwF8t5K9/+XmLauJByYlE9W6C1r3rNnmXGPx1BP49YSRFKd4uPKqVQ3OwVpimIeLMGNS/U2deS3X6LxEWaP9BPBH37YC8TuIqv5HRB4FzgN+mRPpjDhlV28M7N5TltTVp6ioqMV2FKx7T3jCpkLtbtxb/rK2PrHReyqiBsoErXvGVOOfcdWDqwIjlf0KLhdriWHcpGFdqV1tLdfovERRtI0k1jHegeNC9vMecHprhTKCCUq7yaVibMvuPY0bN/LepMkM+dtsSvon/6zaL5lSodIpi3RWa9RAmXTrnqnqGBeLJEQPR1lLTJeeEyZwKsyYrraWa3ReoijaD3AijT3eBI5JGjMK2NpaoYxw5FoxtmX3ns0zZtLw/gcsnPAj3jxsYsKxnn26MeSLefvovJJOWaQKlurbozSyxZYu8CqVDDF3bdZPa0pDeufnijDBZIbREYjiV1wMjPFtzwOOEJE/i8ipIjId+DxOmUajAPjXBWetnsWIu0Yk1CnORjEWonvPnVMXc9uFC+KvP57/EFvufQBBGbT+Jcr2bI+PLSoWhh65X84+u9CkC/xJtVZ5zZcOj/w5QdcqcqOJe6fIr021PxOtzQ8OQ8tgMsftbWkxRkckiqL9G7BWRIa427cAS4DzgYdxmgn8G5iaQ/mMDORaMabq3pPLgKihI/ZNSAwb8u7joO46psYY8t7j8WNSJFSdMiSlrOm22wPpAn+8YKkwhSAgfQRu0LUG93X626b6KWT77JTKQq6prc9JdPCM5TOYtmQaZ4w8gBvPPIJBfbrRbcB89v3YIm488wjOGHlA1tc2jLYgSgnGRfisVVXdJSLHAmcAhwDvAo+o6q7cimikI5dt7QrVvafq1KGsfmE9zU0xyvZsZ9CGlyhWx0Iq1mYGrX+Jdw8+maYefRg+ZhA9e7eMuO4oQVuZ1j3DFoII465NvtaiRYsApz1fEKn2ZyJdpx9/dLBftrAkL4dUlFbwsU8uom7bGr42fDKnHzkoq+85KJAwyrZhtIZIlaGSUdUm4IEcyWJEJNeKsVDde3r27sbwMQN57bl1idasJ4dr1f778G8EWrNtGbSVDWGVaTpaE4Gb67XOoOjkZLKNDk6OE/Co7FvZosVj2O8500OZ/zjQbh/ajI5LqxSt0bbkQzEWqnuPZ9Xut2Vl3Jr1KNJm9tu8krIxlwdasyLCDw4+n89dfj+XTry7oEFbbUVrInBTKcade5ri5RX9eBHFXuehZlUGpykNmYtCF36835xf0fobZUT5njM9lMVisYTjw3Z9mqqZ36Oh5zOU7TyeQbEP+PKnP5bV32EYHikVrYick+1FVfWubM81opEPxViI7j2eVftC088BOPy4Axj7jUoAdm7fw4PTX+XMFGuzAJtnzqTn5p18ZTH85SRnDbSzKllonVXqKcbkhvT+nF5vTLKL2su9TVca8tibFuTUYg5aDvET5XsOE0nvP35RRS8aej5Dw9Zj2bHhi/zww9fixUeC5Gwrd3NX6n7UGUgXDHUncIfvlbwd9PLG5BVx+KaIvCQiO0Rkl4gsE5Hvi0hxinPGiMhjIrLVHb9SRC5NNb4j0VHb2lWdOpSKft2gooFXDng0HszUY58y1p3+HH9958+B5zVu3Mj2hx4CVcatVHrXOed15m5EYasppWL8qMH0KGv5XJ0cLZwq59Ybe/m9K1oEOrVWNj/+5ZDKvpUsn7ycyr6VCWOifs+ZAgaDju/tMRwcTe0FbHlyeHLPWD4jtFyZSBXsZxWzOh7pFO35wDd9r/OB+Tjxos8C1wMXu+/Pufsfdsfmm78CfwaGAnNxKlaVAbcCcyVJ04jIGa7MxwEPAbe543+NU5/ZaAN69u7GOT8bQ93XVnD3e3fGb1zejXZHww5UtcUNZ9OMGTQ1O5ZZCUXct/mrnb71X9QI5SBaU/Dfo1m1xU09F7J5iAgVpU7N7jXb1jDx0Yms2baGvt36As5abdTvOVMkfdDxbgPm46UUJc+J3x2d6jfbWtIp8kKkVxm5JaXrWFX/6t8WkVOALwJnqOojScOvc5XZvcDvcy5lohzjgbOBd4CjVXWzu7/U/fyvAOfiWNeIyD44irgZGKuqS939PwYWAGeJyES3XrNRYIJce33qlDvmlHPkA+czc8VMdjTs4DN8BoDGDRvZ+sD9FDc5AVTFTTG2P/gQP7joKaBzt/5rbVBVGPdzuohij/rGpoRAJ3ULX+TKdXnJqEu46MiL4oFPANv2bGPy8MnxgKiw33OmgEF/gFXZzuNpLhkc7zEMjmV7QJ8eCdfMd2GXoHXl7zz2E17cPI+Grceyp/YgCGgTYRWz2i9R8mh/BDwUoGQBUNW/4xSx+HEuBEvDme77Lz0l635+o++zv+cbfxZOqcg5npJ1x+8GrnY3L8qfuEYmkl13X3k+Rs/NO9k8c2b8hrN+13pUlaevu4BYLPFpXmMxNs/8PVNGT2lXUaLtLc8324L/fsr2e5puA+azrtbJ4suHyxScGt1B7l5vf9jvOVPzh6Kiovjxqz87lSIR9mw4jYatx6LN3SkvLQl0gadyRyeT/J2H6UYUVIgmrmRdl3YQVjGr/RIl6vhIYGGGMW8Dp2QvTigGuu//CTjm7fu0iPRR1VqcrkIATwSMfxbYBYwRkW6quie3ohoQLmfRc931qVPGrVJQEizVre9vZcRdI5j5ShOlyUuIjY3ULViAXPOTvP8tYWm4YX/KmhJ/ToLTbans6o1tIlPUgv8tLVtFindT1m8xFd1KUD01b2lVQe7c7zz2E15bNZYPa3eHDgCat6yG2U9Wsq52Fw8vWBg/x295+gMK521YzeA+xaxzLdnqM49y18QAACAASURBVIM/I0i+SY9NYkT/EVw5+srANKIopSuDIq/TKVmrmNW+iaJoG3CUbTqOxGk+kE88K3ZowLGP+/5dCbwEeL++N5MHq2qTiLwDHO6euzqHchpkzmFMdu2d89hutnAvEKOpoYlFE6+m4rCJ9B+zkwtfvJXFY2qZcuuX27V7WFVbKFmPsqY9bZrnG7Xgf6JycKy9kiKhoc8zjLhrBJDZZRo1QjbI3Rt3nXbbhHJaqKIY85bVUH3fChpjCgg1tfVU37ci8BxP9j7lpSyeOjbt/ATJd/MrNzN7zWxWbXZkunL0lS0eQhLXVh2ZvLXVM0YekPLh06PbgPktlK1A4JxaAY72hYR1Z4nIfThu2/8DblPfiW7w0XdxyjI+oKpfy4Os3md9A5iNU+7xaFXd6u4vwQmM8lzLp6jq4yLyJnAocKiqvh1wPa+G8xhVfTHFZ14AXAAwYMCAo+bMaf1ybl1dHRUVFa2+Tntn/a71sKU7xdrSHVksMXru3siOQfvQVKQMKNuPPW++CQm/SaGu5wEU9xIadwq7S3ZS1kcY2GNgi+u1Kz5cnvrYoJEFESFXv7Ha+kY2bN9NQ3OMsuIiBuzTnXX1b8WPf3LfT6Y9t2ZbPTHfd1okwuC+5fRJU2t5U/0mmrU5/j2vWb+D5qJtqAratE98XEmRUCSyV7be3ePXfePDj2iOtby/FRcJnxy0T4v9EH7OkuUDWL9zPfVN9dQ37fUE9CvvFx+zqsap4S0lHyGixBp7x8f171tPsRTTv9zpXLV+13q21m+lX3k/aj/qSXPRNqS4Dm2uiJ9XVlzEsIG9wsm2a33C9XNJ2DkbN27cq6palXMBOgBRLNqpwDicyN5LReR5YAMwAPh/OBbmVvJf63gOMBk4GXhDRB7Gcf9+Hqdn7ls4ijV12ZpEvMe8lE8cqno7cDtAVVWVjh07NivB/SxatIhcXKe9o6rccuscitb2o0QTb6yHrZ3D4A+fZ8OQ4xn7+AzWX3c9tQ88AI17nSLNUsyHg8ZQe96X2LqkFx995S3ufu/OrAJP7py6mJ21LS3Nnn26cd5Nx2b/RwZx7Rmpj319e+pjOSQfvzG/Necx+cDU34WTY9vyIWtwn+JQlqN3zfOnPorSj1SuUy/cpLy0mRvP/CTjRw3mvKmPprz2u18P/uwocxZkNQJxSx9g5Zkr42N+dNMCamp30W3AYsr6LY6vue570BM0ND3D5OGTOWv0WQntLy8afRF/X76Oqx5cSazfC2hzdxo2f4Hy0mJuPPMIjvdZwjOWz+CjPR9BOcxeM5tJlZNAYOWmlazavCrh+rmkq9zLWkPoYChV/TfwWeCfOG7Wc4Bq9/3jwD+AY1Q1aO00Z6hqDKfn7RXAepwI5G/itPH7f8AWd6i3EObd1XoTzD5J44wcIiJ8+7wvk/wc49U4FpQBH7xA8+bN1C1cmKBkwal9vN/mlQAMHzOI6uN/EA9kiXrDGDpiX4qKE8/p6N2BCkmyy3TlOSszplW1pqKV//sNDvRp+XmFTHMJ+v2lSyNyAs1K4sFWZf0W02v4VTT0fKbFg+PFIy+ObzvpUyPYd89Xadz8hXj61Dr5e/z6qspHDR8xe81sVm5eyaTKScxeM5vZqx139qTKSZ26oEt7J1IJRtf1eqKIDMbpPdsbR0EtU9WCZUu7NZZ/6b7iiEg5MBKoB153d68FqoDDgFeTxpfgWOJNBAdXGa1EVbntzVt5d/89VG74DCU4Vq2/xnFMm1h4/SV84ZlFCefu3L6Hu69+kebGGPvLLqpOGdKqyldVpw5l9YvroXnvDTpdd6Agwqw3qmpKuyt+vAPe8LIp+ZmrOstHHfkKtW++x84PT8WxapVuA+bHLTw/nhLv26M0oRKWR98e2bUHTEeYuuMJwWgbTounEHnH0/0mxo8anLCO63xeYgqQ99yxavOq+Fqxx5VHX9khf3OdhSjpPXFUtUZV56vqbPe9vZQkORvoDtzrpvuAkysLTg5wMscBPYAXLOI49/hvPkPHVdCt1KlbXLqnNrFjT3OMQYveoP6N1bz9hRNp2rQJcMs0HjMQBMp7lsbrHmd7w/Cu51m1RcWSsjtQEGEr8ogIDSXB12wo6dahb3h+Swv2RsemSrfJRdUoVWXwvkJRn+fZ96AnEJR9D3qCsn6LkeLdJFu2nhK/5kuHU5rkwSgtlqz6/WYiUxqRN1/jRw3m+SvHceGZryecn+wRyFR5KigFaPaa2Ryx7xGB8k16bFKbp5d1ZbJStG2NW4Qied9o4CagDqdalcf9OJHKE0Wkyje+O3CDuzkzf9J2Xfw3n+rjL3OUJjA0oGNPMUWsmzKFxg8+YNOMvV9H1alD2Wff7vToXZYTmapOHYoUuUoiojUbpSJP2dUb0Wtq4drt8ZdeU9tmqT25JErJz1xUjfIrlYaez1DhuluP2W88RVvPwL9u61fi40cNZvpZRyZ89vSzjsxbTeAwDyFh3O9hK08F5fI2xBoStof1deZi1eZV3PzKzaZs24h0TQX+gvOo+ENV3eBuh0FV9Vs5kS41T4tIPfAasAMnPecUYA9wpn+dWFU/EpFv4yjcRSIyBydo63Sc1J/7caKVjTzgz1GsOnUo/31jK4OWvtGiYw+NjTS87QSFb3/wQfpffBEl/fvTs3c3zr5hTLy3amvp2bsbnl5oboxx55V73XeZgqKirjd21BrUuSYXVaOC8kr/cMr1/P2Ade2quH6m7zys+z1M5amgFKC129ayb/d92bJ7S3x7WN9hHLX/UezTbZ8u+xtsa9Kt0Z6Ho2hvxokuPi/kNRXIt6K9H5iIE31cDqwD/gTcpKrvthBIdZ6IHI9T3eorOO7lt4EfAL9Re8zLK95/7p69uzH5p8cw/fNfoemhjzHvU7dQX7Yjnj/rRRxrLMamGTMZlKcCFJ8Y1Z+1L29I2BcmKGp4WQOXz/8FVxz3XbZ13+tU6WwVedpjZ5hU9Yr9a5/JhC0QEaagSi5zUsN03Ap6sAhSsv414UmPTWLV5lWcdPBJ/G3t3+LnjTtwHBePvNiUbBuSTtF6BSFqkrbbHFWdDkyPeM5i8l+1ykhD/Obw3iwmT5rMy6MXM23JNOYvuZtTH4Bizy3b2Jhg1eaaY848hLf/tYnmxr3u6zBu5Ku3vkjvXdv4+pp/MGOkk67d2SryRKleVCjCBBoFKZF0rn7vb2lRwWvYdch1Z8QreOWrwlcmyzfVg0V1VTVFRUVxy9gfTTz7lNnc9PJNrNqSGAhV11iXtZxGbkjXVOC9dNuGEZVUbrMj7nwR9K2Esfm0ar2gqDcWf0isWUMFRTVu3Ei/555CUU56fwlzKj9P+cAB7cLaaw3J1uuuhqaMyimXhLEWs4l2hsyu/kwVvGKxWJtU+Er1YHHJ4zdS9NjP4+O8ld+Gp38VV/pSJPGc2bAPJEb+iZTeYxitJchtdujr22hqSgyO8uoXkyf3sT/VJ4w1u3nGTDTmyFhWBI/2fYtBU7+eF9kKRZD1mgq/csqVGzVTeU4/YdytyWRKLcokd1FR+ljRfCmtoAeL6qrqBCXrx6/09ynbJ97lyP9AUlFaEVleK+OYO1L+kkTkoGxfhfwDjNxSiI4zyf9ZD33mGYavWd3idWhSXm0u8acOhbFmtz/00N5iGq5r20tD6qika/KezAF9ynPa7DxsZK2fqMFluWxIX2iSI5jDKv2LR15MRWkF05dOjyvG6qpq6hrrIn1PhWhs35VIZ9G+S5qyhGnQDNc18kRrn0Dz3XGm/o3VvPvVrzL0gfvpXlkZ+fxcl1D0oqCjWLMe+Q7YKgRh+5eWlxZzxYmH8WbDkgQ3ZGs69/itrXz0dIVw3YraM9nMgapS11iX8D15/XbDfk/+hyDvGvnq0tRVSKcQ7yI7RWsUGFWl8WcDWqUkC9FxZt2UKdDcTM0V1XxifmBb47QMHbFvfF3VozUlFL3UoUwElYbMt2u7EKRyrfYpL6Vnt5IWykk1t4oxU2RtLkiXWpSpglcsFktbaKCQCmfG8hmE6cCbiweYQjwEdTXSBUOdV0A5jCzxrNBU5RzCKsnWHs9E/Rur43myDW+/ze41ayJbtbkooZgN+XRhF5LGjRt5b9JkhvxtNiX9+1N90rCENVpwrNdrTz88UDnlWjFmiqzNN14Fr6AHzIaSbpQVFaU8Xl9UzB0rZnKAnpF3i9mrY5xpTPL6dWu+p0I8BHUlOmRlKMMhnRXqpz3851g3JbGCTc0V1ZGv0doSiu2Nxo0bE0pO5pvNM2YmVN6KWrUplWLMZh3fvyZb2beSFWeviFdImjB/Arctuy3yNbOhRQWvQSMTKniV/mgDN518FUcMPYibT/khek0tN518FUcfPJjlNR9y1YMrM5bkbC0iwpWjr6S+qGUXJGhZ1jMX31Muv2vD1lI7NO1BgYbBb816pLJqk62uZKJGC3sUtEVeSPyKL99rvfGALtWEHOWwVZuyzWVNhYhQUVpBZd9K1mxbw/Sl06muqmbp+qWs2baGqoFVBXPNpguyEpF4Qf5kN+rfF1RR37g74dx8pUKJCN1/vIUj/C34znG6WpUFKNnWfE+5/q6NLCxaERktIleLyEwR+UvA68/5ENTIL5meVFvzJJtszXoEWbXJVlcyUaKF/bS3FnnJii+TVZttNLg3zh/Q5QVyRSFVLmvYloVB8l8y6hLmnjY3bskeefeRrNm2pt2tBfrXLD2mjJ7Ch7W7A8eHDTKLQioLM5nWfk+5uoaRSGiLVpzZvROn7KHTpyqxC7P69ue7BKMRgbBrtGnXq1rxn2vPv/8dGHSy59//TthOZXUlEzZaOPmctljfTUWQ4ktl1UbJN/WzqX4T05ZM4wcHnx+YnhS18lY2uaxh5G/va4GplNygPlWsC1C22ZTkTJcxENXCzPZ78pOLaxh7iWLRfhenDd3dOP1dBbgFGAP8EKe4/xycJvBGAYglpZwEEaUtWz46zqgqj/z1XL52VQkP33Uelavf4OG7zuNrV5XwyF/PTbB0wlpdXrRwlLXZ9rS+G5SXu2Hu/Zz8k4darO9lk2/qndeszcxaPYunr7sgZXpSVKLmsmaSPxaLpVwLLPQadhDpuu186ohFlJcm3kKzydMN0xKvLSzMqN+1kZooa7TnAmu9aGR30mtV9SXgJRF5EngJeBq4I8dyGkl4VsJlaazQ0h9tiGyJ5vo/V9hUgVRFIXJZ7zjb9d1cE5SXKxpj3CvzuUoda8hb40s3f171H49kC2Rgj4FMHj6ZfX57JyT3Py9QelIm+f05nsmW2jmP7S7YGnY6+dOVfzzlzBGtijoOm7MaxcLM1gNi5I8oinYYTm5t4PmqukxE5uOU4DRFm0OS3UqxWGzvf84TL29xw6quqqasAOkRYQmTKhC1KEQ2wU2eVfvac+tyZs1mI0dQXm6ZNnPM+teY0Xhmi2CaoPnzqv9kuplOGT2FEd/be97Kc1YW3DJJ9/2nUmL77iwKtYxQCDIpudYEPkXJWQ3zEGzFJtonURStANt92zuBfklj3gJObK1Qxl6Cnk6nL51ORWlF3IXV3hPK07U482QNUxQilVLzCBPclM36bjqyKaLh5eUOnfpoYEWY5GCaoPlb9P4i1mxbA6S/mWaa90KQ7vtPpcTWX3c9tSHXsAtBPt2oucxZtWIT7ZMoirYG8D+6/Qc4KmnMoTgK2MgBmZ5Oq6uq23UQCSSucV2wYjpFu9zSGi/CjL8sBFzrL0RRiCCl5ieMOzhsNaiwtCbIKlPRe0gfCFPZtzLlzVRVWb9rvdOSsA1TNLJJFWnatCnvywjtiTAPolGwYhPtjyiK9hUSFevjQLWI/Bh4EBgLnAHMz5l0XZww61t+2sJayYTfPXhE6RBWt6KEYpBS818nF+7gTHm8yaRquXffjUszupRTVWbyB9Okc69WlFbErVpIvJmKCMVSHLm1XK7JtMYZJEdnrS0dRD5yVnOtuI3WE2Uh7wGgWESGutvTgPeA64CVwG+BWmBqTiXs4gTl8F1x1BUJa7L+qjrtsXqLl8Ix+tShSFGSCy5CUFJy5LAUOa+o10lHpjzeoEjYKt/f5ckRJm83bGWm5E4uIns7svhJ/u77l/dvcZ7nri0kQfKnkyPtMkInI9cRxemipNvjvaGrENqiVdV5wDzf9lYRGQV8G/gETrefu1T1w1wL2ZUJejqd+OhExh44NsGy9dZs22tCuYhk1XA9Gb9VW1RcxCFH7c/al9fnzJrNFIATVM0pKMgqrEs5bGWm5OjidNG6yUot1XWiWu+tIcoa5+JjfpbSG3BoziVre3KZs5qNB8HIP60qwaiq24Ff5EgWI4nkp9PqqmomzJ8QdxfOOXVOi2jjQhRjj4r/JtLaFJtkpVZ1yhA+fLs2Z9ZsuiIS6RRxcpBV2IeKbKKWc3UzLWQJyLCkC3hrqypehSCXwVZWbKL9EaUy1MnAk6qauUqCkROCbqhzT5sbV7YjZ40E2ndEYXKP257A8JILeK3xJIaPGZyVFepXarkKbkqVx7tgWxW1u7sDcNjaOQxqaKIYaG5oamHVJssR5qEi29Z/2dxM/ePDVuEqNCkD3oQ2y3vuiFixifZFFPPnUeB9EZkmIp/Kl0BGIsnrW0VFRcw9bW7CmPaqZFN1F6qquI99ijdSdfLBWV03m8pQmUgVgFO5+WmKioWyPdsZtOElitUJXCrS5ow1isPUZa5qxbp1lJtpcvWhzTNm0tTcGP87N82YmXU95VwSNB8Awz47sMN2aTKMKIr2dqA7cAWwQkSWiMh3RWTf/IhmeAStz/lpr0EOqW78PYu3cXb/i+jZp3uBJUpNqgCcXu/9CykShrz7OGj0MoZVpw5ln327p1Sc2ZSGvHPqYm67cEGL151TFweOTy6D2LhhI1seuI/iplj879zywL3c+vS1KcsAFork+QBA4JjxnyioHIaRS0IrWlW9EBgETMBJ7TkSuBWoEZEHReR0EbG2e3mkK0YUelG+NDXl9XMOfWYRw9esbvE67LlnGH7MQPpvXRW3ZvcKlzkSNoz1HRS1nI50Ec1BVqk/knXW6lncdulYYrGkvyWm6J33RaqnnC+SrVqzZo2OTiTFqKoNwH3AfSKyP04nn3OB8Tg5tFtE5G+qemnOJTW6ZEShF7DT2IaF5atOHcrdL95Ic2OM4tIizr7hmJze+KOWhkwV0fzawYt4acm8hKWEoC45VW8ppUl6trgpxvHv9ODcdlBNyD8fpWXFZs0aHZ6sQ1RVdaOq/kpVjwRGAb8BegPfy5VwRkui5iS2Ja3tcesP2Gnetq3Nurhk2wM3CplczEHy+N3NlccMpK6kNsG7sX7X+sAuORd9ryTeTemQZ56h9MADOfS5Zxn1wtKEz2nLtX9vPiZd/1mzZo0OT6tzQUTkMOBrwJlAaaslMjLSUSIKRZwet0GEad+XHKCUTVu3XBFFEWZD1ACvZHfz6FOHJriHR9w1gq31WwO75PiXHJ6+7gI3xWdGynZ1bUE+At4Mo63Iak1VRPoAE3HcxkfjNBz4CPgzTnN4wwDcHrcBTa0zte9rkW7Txmkoyek7hSz2kEqeIHdzlC453Wp3cdBzD4AqW+6/j0f2g8mjz26zusiG0VkJbdGKSJGInCoi9wIfArcBo4F/4qzVDlTVb6tqcOij0WXJxgJPV++2PZCpVGMhSLayU9W4VdXAJYcJLxRR7N0CVPnhG4cUtLG4YXQVoriO1wEPA2fh1Di+GjhYVU9U1b+p6u58CGh0TdpzvdvkYg9tuXbsuVeDItL7lfdLWLNNLr/o9xgUN8X4+PPv0rx5M9C+1/4No6MRRdF2B/4IjFHVSlW9UVVr8iRXKFwL+ykR+UBE6kXkPyJyn4gckzRuiIhomtectvobjGCS0226H344w9esjvdybUuCSjW2NUHu4QFl+3HHn7qz766iFlZpGI+BWbKGkRuirNEOUNXUXbcLjIjcDEwBtuA0O9gMHIKTZvQVETlHVWclnbYCX2MEH6/lU1aj85CqVGN7KGGYXJaxadMmem7eyanP1DtNLH2k9Ri0k7rHhtFZiNK9pz0p2YE4Fao2ACNUdaPv2DhgAXA9kKxol6vqtYWS08g9bR2E1N57pfprGTdv25YyiKw9eAYMo6vQ/lq9hONgHNlf9itZAFVdCOwA2r5CupFz2joIqT2vHfvZ7Juf9uLeNoyuSkctmfgW0AAcLSL7qepm74CIHAf0IthFfICIfAfYF8fl/KKqriyEwEbraQ8dZzqCJRifpwu+7e5oP+5tw+iKdEiLVlW3AlcCA4A3ROR2EbnRTT16Cnga+E7AqV8Afg/8zH1fISILReSgAolutIL2GITUHmnvqVGG0dWQjlyIXkTGA38B+vp2vw1co6p/843bH/gujpX7H3f3COBaYJx7zkhV3Znicy4ALgAYMGDAUXPmtD5Iua6ujoqKilZfp6uwY8cOSt9/H/y/VxG6H3YYlHRUx0x+2LP2TbSpkYb996ds496VFSkppduww9pQsvaP/b+MTtg5Gzdu3KuqWlUAkdodHVbRisgU4Oc4NZZ/B6wHKoEbgROB6ao6JcM1SoDngc8Al6rqrZk+t6qqSpcuXZppWEYWLVrE2LFjW32drsLT99zDx35+Y+L6aGkpfc46q10EIbVH7DcWHZuz6ISdMxHpsoq2Q7qORWQscDPwsKr+QFX/o6q7VPVfwJeBGuByEfl4uuuoahPwJ3fzuHzKbLSO2Ec7OkQQkmEYRjLZ1jruCfQBioOOq+p/WyNUCE5z3xcGfPYuEXkFR+GOYq+rOBVeWZ+euRPPyDXdhh3G8DWr21oMwzCMyERStCJyNk4Q0vA0wzTqdbPAa+mRKoTS298Q4lqfdd8zKWTDMAzDiExohSgi5+EEHjUDzwHvA035ESsjz+EEN10gIn/wl4IUkZOBY4HdwAvuvs8Ay9zG9fjGngBc5m4mF7cwDMMwjFYTxfK8AtgG/D9VbWsf3v3AP4DPA6tF5CGcYKjhOG5lAaaq6hZ3/M3A4SKyCPjA3TcCOMH9949V9YUCyW4YhmF0IaIo2kOAv7YDJYuqxkTkFOASnL64XwZ6AFuBx4DfqOpTvlPudseMBk7GaVC/AbgX+J2qPldA8Q3DMIwuRBRFuxXHHdsuUNVG4Bb3lWnsn3Ga0huGYRhGQYmS3jMfGCvWO8swDMMwQhNF0V6FE+37exGx0imGYRiGEYIoruP7gF3A/wLfEJG3gNqAcaqqn8uFcIZhGIbR0YmiaMf6/t0TGJliXMes6WgYhmEYeSBK4/cOWa7RMAzDMNoSU56GYRiGkUdM0RqGYRhGHomsaEVkooj8Q0S2iEiTiGwVkadFZGI+BDQMwzCMjkyUWscC3AV8A6fEYTNO55v9gM8BJ4jI6ar6jXwIahiGYRgdkSgW7XeAScC/cGoMd1fVQUB3d/tVYIKIXJhzKQ3DMAyjgxJF0X4TeBc4TlUXqGozgKo2q+oC4Hj3+LdyLaRhGIZhdFSiKNpPAg+pan3QQXf/PNL3qjUMwzCMLkUURas4a7PpsDrIhmEYhuEjiqJdDZwpIuVBB93944E3ciGYYRiGYXQGoijavwAHAc+KyOdEpARARIpFZBywEDjYHWcYhmEYBtFqHf8B+B/g68BTQExEtgL9cBS2APeq6u9zLqVhGIZhdFBCW7TqMAknxWcBsB1HyW53tyepqhWtMAzDMAwfUSxaAFT1HuCePMhidHJUFafuSfC2YRhGZ8RqHRsFYcbyGUxbMg1Vp4uiqjJtyTRmLJ/RxpIZhmHkF1O0Rt5RVXY07GDW6llxZTttyTRmrZ7FjoYdceVrGIbRGYnsOjaMqIgIU0ZPAWDW6lnMWj0LgMnDJzNl9BRzHxuG0akxi9YoCH5l62FK1jCMroApWqMgeO5iP/41W8MwjM6KKVoj7/jXZCcPn8zKc1YyefjkhDVbwzCMzoqt0Rp5R0ToVdYrYU3WcyP3Kutl7mPDMDo1rVK0ItIfOBanKtSLqro+J1IZnY6LR16ckDfrKVtTsoZhdHbSuo5FZKSIXC8iIwOOfQt4D3gAuB94T0T+Lz9iGp2BZKVqStYwjK5ApjXaycBVwIf+nSIyAvg90B1YDDwBNAO/EpExeZDTMAzDMDokmRTtscCrqrohaf93gWLgp6p6nKqeCpzoHrswxzIanZDkACgLiDIMo7OSSdEeCKwJ2P95YDdwk7dDVZ/HaS5wTM6kS4OInCoiT4nIByJSLyL/EZH7RCTw80VkjIg8JiJbRWSXiKwUkUtFpLgQ8hp7mbF8BjcvuTmhHOPNS262coyGYXRKMina/YAa/w4R6QcMAV5W1fqk8a8Dg3MmXQpE5GZgPvBpHLf1rcC/gDOAxSIyOWn8GcCzwHHAQ8BtQBnwa2BOvuU19qKqPF/zPLNXz44r25uX3Mzs1bN5vuZ5s2wNw+h0ZIo6bsRphednlPv+r4DxO4G83ilFZCBwBbABGKGqG33HxuFY1dcDs9x9+wB/xFlDHquqS939P3bHniUiE1XVFG6BGLHfCFZtXsXs1bOZvXp2wn7DMIzORiaL9m0cN7GfE3GU6UsB4wcB+U7xORhH7pf9ShZAVRcCO4D+vt1nudtzPCXrjt0NXO1uXpRXiY04IsKVR1/JpMpJCfsnVU7iyqOvtEhkwzA6HZkU7SPAJ0TkdhEZISJn4QQ77QaeDBj/WeCdHMuYzFtAA3C0iOznPyAixwG9gH/4dp/gvj8RcK1ngV3AGBHplgdZjVQk61PTr4ZhdFIyKdpf4uTKfgtYBszFUWS/VtWP/ANF5JNAJY47Nm+o6lbgSmAA8Ib7EHCjiNwLPAU8DXzHd8ow9/3NgGs14TwYlAAfz6fchoN/TdaPf83WMAyjMyGZbmwiMgj4KU408RZgrqreFjDuCmAScJ6qrsiDrMmfNx74C9DXt/tt4BpV/Ztv3JvAocChqvp2wHUWA2OAMar6YorPugC4AGDAgAFHzZnT+uXcuro6frg+cgAAIABJREFUKioqWn2djsg729+hvqmeft37MbDnQNbvXM/W3VspLylnaO+hged05fnKFpuz6NicRSfsnI0bN+5VVa0qgEjtjowlGFX1Q+B/Q4z7BfCLXAiVCRGZAvwc+A3wO5x14UrgRmC2iIxU1SlpLpFwOfc95ROHqt4O3A5QVVWlY8eOzVLyvSxatIhcXKcj8sbyN/hoz0dcdPRFiIhj5b5yM3SDsSPHBp7TlecrW2zOomNzFh2bs8x0uKYCIjIWuBl4SFV/4Dv0LxH5Mo6L+HIR+b2q/gfY7h7vneKS+7jv21McN3JMUN1jC4QyDKOz0hHb5J3mvi9MPqCqu4BXcP4uLw1prft+WPJ4ESkBhgJNwH9yLqmREqt7bBhGVyGtRSsi2QQ2qap+Lkt5wuBFB/dPcdzb3+C+L8BZO/4icE/S2OOAHsCzqronl0IahmEYBmR2HY/N4pr5Dht9DqfW8gUi8gdVjVeuEpGTceoz7wZecHffj+Nqnigiv/UVrOgO3OCOmZlnmQ3DMIwuSiZFO64gUkTjfpw82c8Dq0XkIZxgqOE4bmUBpqrqFgBV/UhEvu2et0hE5gBbgdNxUn/ux0lbMgzDMIyck1bRquozhRIkLKoaE5FTgEuAicCXcdy/W4HHgN+o6lNJ58wTkeOBHwFfwWnv9zbwA3e8JW8ahmEYeSHTGm2pqjYWSpiwuDLd4r7CnrMYOCVvQhmGYRhGAJmijreJyOMiUi0iny6IRIZhGIbRici0RlsCnITb1F1EanHSav4J/ENV38qveIZhGIbRscmkaPvgRPGe4L6qgDNx1kURkRocpftP4J9uFSnDMAzDMFwyBUPtZq8iRUR6AcezV/EeAZwLnOMeX4tj6X4/jzIbhmEYRochUglGVd0BzHdfiMi+7FW6J+LUGx4GmKI14sRiMYqKilJuG4ZhdGZae7c7yPfav/XiGJ2N8584nwnzJxCLxQBHyU6YP4Hznzi/jSUzDMMoDJEUrYhUisjFInK/iGwGlgLTgZHAwzit5FrUFDa6JrFYjLqGOtZsWxNXthPmT2DNtjXUNdTFla9hGEZnJlMe7cHsdQ2PAwbhVF7aiFOdaRGwUFXXprqG0XUpKipi7mlz48r1yLuPBKCybyVzT5tr7mPDMLoEmdZovY42W4BncVJ7FqrqG3mVyug0eMrWU7KAKVnDMLoUme52Xu+ylcCrOK7iNXmVyOhUeO5iP/41W8MwjM5OJkU7CfgLTs/Wn+F0xNkmIvNF5HIRGZX2bKNL41+TrexbyYqzV1DZtzJhzdYwDKOzkymP9h7cHq7ueu3ncdZrx+LUDVa3WtQzOG7lBar6ej4FNjoORUVFVJRVJKzJemu2FWUV5j42DKNLEDqPVlXfA/7svhCR4cDncBTvccAZ7v5Nqjow96IaHZE7vnhHQt6sp2xNyRqG0VWIVLDCj6quZm8/2C8Al+FUiuqfI9mMTkKyUjUlaxhGVyKyonWrQY3DsWQ/BxziHQIUJ3DKMAzDMAxCKFoRqWBvfePPAZ/CUapeRPK/cWohL8BZo92cH1ENwzAMo+ORqWDFCzgde4rZq1jXsbdV3j9V9f28SmgYhmEYHZhMFu1ngW04FaD+iWOxWh6tYRiGYYQkk6KtApapqhZCGMMwDMPobGTKo/1XoQQxDMMwjM6ImLEaDRHZBLyXg0vtB1jgWHhsvqJjcxYdm7PohJ2zg1W1S6Z/mqJtI0RkqapWtbUcHQWbr+jYnEXH5iw6NmeZscoBhmEYhpFHTNEahmEYRh4xRdt23N7WAnQwbL6iY3MWHZuz6NicZcDWaA3DMAwjj5hFaxiGYRh5xBStYRiGYeQRU7SGYRiGkUdM0RYQEfmYiPxFRNaJyB4ReVdEbhGRvm0tWz4RkX1F5H9F5CEReVtE6kVku4g8LyLfEpGipPFDRETTvOak+axzReQVEalzP2ORiJyW/78y97i/j1RzsD7FOWNE5DER2Soiu0RkpYhcKiLFaT7nNHeetrvz9rKInJu/vyw/iMh5GX43KiLNvvFd5ncmImeJyG9F5DkR+cj9+2ZlOKcgv6WONpfZkHXjdyMaIvIJ4AVgf+DvwBrgaOD/gC+KyLGquqUNRcwnXwVmAh/idH76LzAAOBP4E3CyiHw1oKb2CmBewPVeC/oQEfkFcDnwAfBHoAyYCDwiIt9T1d/l4G8pNNuBWwL21yXvEJEzgAeA3cBcYCvwJeDXwLE430PyOd8FfgtsAWYBDcBZwJ0icoSqXpGbP6MgLAeuS3Hsf3BafT4ecKwr/M6uBo7E+d18AFSmG1yo31IHncvoqKq9CvACngQU+F7S/l+5+3/f1jLm8W8/Aec/aVHS/oE4SleBr/j2D3H33RnhM8a457wN9E261hacG8aQtp6LiPP2LvBuyLH7ABuBPUCVb393nAc8BSYmnTPEnZct/rkB+rrzqMAxbT0POZrLF92/5/Su+DsDxgGH4rQ7Hev+DbPa8rfUUecym5e5jguAiHwcOBHnxnlb0uFrgJ3A2SLSs8CiFQRVXaCqj6hqLGn/euD37ubYVn7Mhe77z1R1m+8z3sWZ827A+a38jPbMWUB/YI6qLvV2qupuHGsG4KKkc76JMy+/c+fJO2cb8HN380I6OCLyKZyWnzXAo628XIf8nanqQlV9S11NloFC/ZY65FxmgynawnCC+/5UgLLZASwGeuDcDLoaje57U8CxA0TkOyLyQ/d9RJrreHP8RMCxx5PGdCS6ichkdw7+T0TGpVgjS/f3PwvsAsaISLeQ53TkOUvmO+77n1W1OeC4/c4SKdRvqSvMJWBrtIVimPv+Zorjb+FYvIcB/yyIRO0AESkBznE3g/6zfcF9+c9ZBJyrqv/17esJDAbqVPXDgOu85b4f1lqZ24CBwN1J+94RkfNV9RnfvpS/MVVtEpF3gMOBjwOrQ5zzoYjsBD4mIj1UdVdr/oi2QkTKgclADCceIAj7nSWS999SF5pLwCzaQtHbfd+e4ri3v08BZGlP3AR8CnhMVZ/07d8F/BQ4CmeNpy9wPE4g1Vjgn0lu9s46v3cAn8NRtj2BI4A/4KxhPS4iR/rGZjMHYc/pneJ4R+BrOH/z46r6ftIx+50FU4jfUleZS8AUbXtB3PcuUw9TRL6PE224Bjjbf0xVN6rqT1T1X6pa676exbH6XwYOAf43i4/tUPOrqte569sbVHWXqr6mqhfiBNCVA9dGuFw2v7HO8Lu8wH3/Q/IB+51lTSF/S51iLk3RFoZMlsE+SeM6NSJyCXAr8AYwTlW3hjlPVZvY6/47znco0/xmenruaHgBZFHmIOg3FvacjyJJ104Qkf/f3pnHRXVke/xX3TSg0AjIjgiIisuoQVxwm4jOqDHiEjV5JmY08QWN4Wn0JWYZk2iWz2SSGJfEiLgRzLxMnIx7jHvUuCNuSSZBouCGAiIIItJNd70/6l5tm9sLS0Mj5/v59OdCbbeq7rl1bp3aOkHMbL0CYJu98UjO6kWWmkpdAiBFW19kSldL4w3tpKulMdyHBsbYywA+h1ijGC/NPK4OBdL1nkmPc14GMaPUkzEWrBDnYavffOlqata0KGPSWHgkxISzC3bGCZbSv9JYx2dhexKUNZqynDlclppQXQIgRVtf/CBdh7CquyBpIRaAlwM4Wt8Zq08YY69BLHg/DaFk821EUUKemX3BzH2vdB2mEOcxszCNnT7S1bQOrJX/jxCz2g9zzivsjNOo64wx5g4xJGEEsKoGSTRlOasvWWoKdSlo6IW8TeWHJrxhhVTOt6RyngDgayNsbwCuCu6DIBaxcwB9zfweqsXvELM6q9QTgHCIGZkcwJsm7l4QvbDqbDIQiYd0wwoIJcsBbCE5q1K+gbC9YYXDZelhqEt7f3QebT2hsAXjrxAvejyEeaQvf0i3YJT2Ok0FYIDYok1p3CWHc54qhd8HoWj2QYyvAUBX3F9T9xbn/H2F+ywAMFuK8y3Edm5PAWgJ8YHTaLZzY4zNA/A6hDUkG0ApgCgAj0M0eNsAjOGc60zijIYo910A/4TYNm8kxNKLbwE8yc1eeMbY/wBYAtGwfYP72+a1ArCAN64tGO/BGPsRQH+InaC2WAizD01EziTZGC39GwRgKERv/UfJ7Ybps64vWWqMdVkjGlrTN6UfgDCIJRvXIITwIsSkIKs9vMb+g5gdy2389pmEnwJgK8ROWrchvqwvQby8A2zcaxKAdIjdtkoB7AcwoqHroAZ19iiAryFmZRdDbOxRAGAXxNpjZiFePwglXAQxHPETgFkA1FbulSDVU6lUb+kQa0gbvB5qWHcdJZm6bKPcTUbO7HgHcxpKlhpbXdbkRz1agiAIgnAgNBmKIAiCIBwIKVqCIAiCcCCkaAmCIAjCgZCiJQiCIAgHQoqWIAiCIBwIKVqCIAiCcCCkaAmCIAjCgTRpRcsYC2KMfckYu8IYMzDGOGPMW/LzYowtYYzlMMYqJb9HnCDPXNrRhqgnGGOTpXqf3NB5kXHGPD2sKL1zjLF5kvvAaqSzjzHm8I0LqI1wPlwaOgN1hZ0CHM8532fyfyrE2ZNf4/5+nHclv48gTv/YCmAtxPaB1T1pptowxnIAgHMe4eh71RWMsT8DmA6xpaQfxA4vBQDOQGzx9hmnnVGswhhTA3gewESIA961EDvyXAdwHMBmzvnmhssh4Sw0xjbCUUgfmmsAPMelLVydkYdG0Zow34pfjvwHY8wVwJ8B7OacP6MQdgSAc5zzhLrNXq3pCMBpji1jjL0J4AOIY7O2QxyXpYHYZPxRiL1Ov5D8CQUkJbsV4hSTYgDfQez96guxv/HTADoAMFW0GyBOe7pWr5klZD6H2AP4UkNnRAGnaiOIh1DRcs7n2Rk0CMJ0nmvBPwTAgbrIU13COf+tofMgwxgLB/AuxGHO/TnnP5n5qyA+Zqp7FmhTYwKEkj0D4FHO+QOHLjDGmkNYC+4hhXkoDsVujHDObwC40dD5UMKZ2ghCoqE3W66rH6TNse0MmwPljbVTIU7ysLrpvZTGUIgNt29AbEZ+HsDHALwt3LMVxMkWWRDm6ZsQJsG3JP+BFu7LAaSalXOfyf/LJbeRFu4bJ/n/y8y9OYA3IM6GLYPYVP0IgAnVqPMnpbQ3VvNZTQbwb4jTQ8ohFPUhABMthJefiQbA21Jd34XYcP8Fk3DTIDY+L4foEc4HoDJLK8LkWXcAsFF6FmUADgIYYiG/HMBkC8/1c6ksFRAnl2wG0LMa9fGFlP7L1azDB/Iklam6G8dPgDghqEiq018BzAXgVs1nGgxgKcS7pYMYOlgPINZa3iFOr9oHsZl8CURvvqOd95wgpfOpBX833De/u0huLQC8CnHO6RWTvG4GEGchHaX3f57kPlAh/H8ByJDkMB9i6ClElmOzsK4AkiDakouSDN0EsBvAY2ZhB1p5thbbCBP3FgD+BmF1uivVzQ4Af1IIK99rHoBHpOdSDNFT3g+z4wPteFYjAeyBsMBUQHRw9gOYrhDWV8rnr1Id3pLiDjELt89KfURUJ3+O/j10PVo7WQTR4M6E6EVslNxPA/CGeIDvQAh+quSXI0dmjL0N0YjfhDD55UMcr/UKgOGMsT6c8xKT8D0gBNoXope8HkLRdYIQ5Pek9OcDeNkkjzKnrZQlFUAixAkYSmN4f5GuX5rkxxuioYkBcBLAaoje/VAA/8cY68w5n2vlnjLysX5tGGNqzrm9PddlAP4DURfXII7EGg5gLWMsmnP+loV4/4To2W2DOM1mHIAUxpgeov4nQTyPPRAv9tsQDcPfFdKKhPiw+BniYyUY4niu7xljT3POv7FVCMZYdwA7IZ7rDojn6gdxHNlBxtgYzvk2W+ngfj22tyOsNTbCRE5N6ALgCZiZExljqyDGha9A5L0Y4sPsPQCDGWN/5pzbNPkzxiIhPlJCIOTqa4iTqsYDeJwxNpZzvlUh6ggAowB8DyAZ4n0YDqAnY6wTF71Ga2yAaISfYYzNUcjrKIj3eYGJX0eIoY4DEMqjCEBrCHl5jDGWwDnfbqvMlmCMzYI4Y7oYQJp0HQpxRKaSBcIX4gSvwxAnMxVAyGICgG2MsRc45yulsDmoWRshv/OHIOo4XYrrB/GxvJMx9iLnfLlC1B4A5kC8Kysh6mosgD2MsUc455nW7ivdOxHiHbsOYAtE5yQA4p19DuJDUw4bDtH+RkDM79gOwANCVrYzxqZyzldIwVMh6ncUxNGjpnVQbCtf9UpDa/q6+uH+l8w8C7/XzcJHwOxLUCE9pa/CeMnvMMx6r7j/pb7QxM0V4jxRDuBphfTCzP7PgULPw1q+IL5QKwC0NHN3g/gYyIP0RS+5p0rpzDEL7w4h2EYAj9hR5x64bx04ANFod4aVY7SkeFEKbq4QClIPINTMb590j3TTOgfQBqJHUiTVcaiJnzfEC11gVnb5uXMAH5vdp4d0/yIAXgrPdbKJmwvEBLq7EOZe03RCAFyF+Iiw2TOE+ODRSfW+FkIphtuIUyVPFsK1glCk5TDpsZnEXw+gmVmceZLfTDvfvR1S+L+aufeFGJsvBOCpcO9KAIPN4vxNSTat3Fu26FQ5Vg1CkXIAXUzcWgDws1BPuQB+tfOdk+tooJlsyT3SCBN3FYQFh6Nqj9YNQCuFe7aA+Ai8qfB8clD9NkKup+UwOWIRQDuID4AKszwPxP33ZLJZWlMl9y/sfEYZUvoBCn5+Zv/vk94D84PlvSEUaTmAwOq+Bw39a/AM1FlBbJ93WmwWPgI1U7QbJL/OFuKdApBv8v9YKfwmO8tRk5foTcn9JTP3cTAzrUH0HisBpFtIv5sU5yM789tVKrNpXcvmpemohgkSQsFwAH8xc98nuQ9WiLNX8ntewW+N5Beu8NyLAWgV4qRK/pNM3Kq8zBBf0VWUtYn/TMl/uJ1lfxJCMZvWY6EkbwkK4W02MBAzl09DNFzjFORUD4WhDgBqiI+U43bku5WUj4sANAr+a82fqUnev1IIHyn5fWtnvfWF8tBIkCTnJ6shf0uktFqbuduraP8quc1XSLsNxFwFXo38zJbS+6OZew6q0UZADLnIZ71WOfsawoLBAbxt4jZQcjuoEF4jyc4JO8uRId3fx0Y4ue35lwV/+Z2bbuJm8z1wht9DZzrmnDMH36IPhJCNZ4yNV/B3BeDPGGvJOS+EMMUBwjzmKNIgXpZJEONkMpOk65cmbj0hGlLOGJunkJZGuna058ac87MAYiTzeDyAWIgy/1H6JTLG4jnnRXIcxlhrAK8BGAxhimpmlmyohdudUHCTJ7NlKPhdla6tIBSBKSc556UKcfZB1FsMHqw3c/pI13AL9dhOunaEMHVbhXO+jjG2AaIO+0v37w9hhh7NGEuDaEy4rbSAezOZ10E0XnM459+a+DWX3G8AeJkxxVemAvbJQIx0/ZFzrlfw3wuxZCkGQk5NUXqel6Wrjx33Buf8MGPsHIAExpiPiZw9AyHnqeZxGGP9ID6E+kCYMF3NgoSiZrOJu0vX/Qr5vMAYuwwgXCE/nSHGjf8IYTZ2V8hPbegAMVR1iHN+U8F/L8S4fIyCX5VnxDnXM8byYOczAvAPAAsA/MIY+waifg5xzgvMwsnvVAsL75S/dLWrbXImHjpFWw+0hKi3d2yE84TokXhL/1+1ErZWcM6vMMb2APgzY6wj5/xXxlgAxEzW05zzMybBW0rXntLPEp7VzMMJmLyUjLFeEIqqG0RdvSy5t4GYBOYDMQazE8J0ZYDobU6CMKcp3UNpjEsee7Pmp1Hwy7NQFHmtdAsL/jJyPSp9bJlidz1Kimqn9JOV5ViIMfS/QPRuN1pM4EGWQjz/5Zzzj838fAAwiIbLlhzbQq4nS8uMZHdvBb8q42ic80pJ8aurkYcvIcZd/wti/B8QcqSHGC++B2NsDIBvIUz+uyAm1pVB9PoHQixJU5Q/O5DrwppsPaBoGWNxEIrOBWLoZDPEpDAjxCSkUbXIj3m+6uQZSVTCzmfEOf+UMXYDwsI1A6It4Iyx/QBeldoO4P479WfpZ4lqtU3OACna6nMLYiarr53hZUGt7VepLb6EEM5JAF6H+KJ3QdVemayQFnLOZzsqM5zz44yxJIjZk4NMvGZDvFBVFpgzxibgfi/c0QRacA+SrraWzsj+o7iDNpLgYnLZOsZYF4gexyDYoWgZY3MgxtG2A3hJIYic91Oc8+4K/tVBTivIgn+wWThHsBb3LTrLGGMxEBPANin0mt6DGA/vwTn/1dSDMbYcQtHWFLmMgQB+UfBXqqO5EBYd8810wBh7A0LR1pYGf0ac8zQAadKkrL4AxkDM6dghdQ7yTe4/k3O+xFF5aQia9BaMNeQoAB/J3GNveAB4zM7wBlTva15mPcSX8ERp/eokiK/O/zMLdxzia3lADe5RXWTTrKltsq10/bdC+No0ctWlO2NMq+A+ULqeshFffq4NVY+KMMbGAfgQYjb9k1xhJjjn/DaEIujMGLP3g9EScj31Z4wpfbjHS9eTtbyPRTjnlyF6hb0ZY9FQHjKRaQvgPwpKVgVhqq8NchmryLFkyQmzkJ+b5krWUjoS1W0jMiHmTTzCGFMy9zr8Gclwzos559s45y9AmPV9cf8dqsk7Jct3TdrMeoMUbfVZKF1XMMZCzD0ZYx6SOUhmC8TkhZFSj808vHlPtxBijNd83NIqnPNyiDG5UACzIEy226QvRdNw+RBjJj0YY28pNY6MsShpyYZVGGO9pD13q+SVMaaBGIcFHtz4I0e6DjQLPxTAf9u6Zx3SAmL5j2keekBYAm5BmGmtsQnC7PgSY2y4UgDGWB9pPNQqjLEJjLE/S429uV8QgBekf61uoCLJ3VqIcesRFsagZT6FGJtcLfUyzNPykZYvWYVzfgXCBBuB+8tO5DR6Q+xqVQTb9VlbUqXrFIj1tYUQS73MyQHQzvTdZcJW/Q7E0pfa8A8Ic/X/MMYiTNJXQayxV2pvcwD4Msa6mjoyxqZALAtSolptBOdcJ+XNE2KDGdP7REGYc/UQslPnMMaGWfgIC5Cud6R8noAYTnqCMfa8hbS6SMNiMvLSuNZ1lV9H8NCZji0Mosts5JxbXW9mC875HsbY6xDLELIYY9sglpZ4Qoy/PAqxpnCYFF4nTZraCbFGdSrEl5s7xKD+YDz4HPZAjJ1uZ4wdgJiUcoZzvsWO7H0Joaz+ZvK/EkkQk3XeBfAsY+wgxLhSiJSnnhCNVbaN+4VAzOz9XErjPxBjX8FS+YMglsCYvtxfQKyd+xdj7N8QY9d/kMKvg1jLWh8cAPDfkjI4hPvraFUApnKTddBKSBNCnoBY2vIdY+wwxAzfOxA9l54QM02DYXs7vN4Qk3OuS/Uo13skgMchTIubIMYWrbEaQq6OSWUz9y/mnC+S8r+aMRYLMW52njG2A2ICkK903z9CPNtpNu4JKcwhAB8zxoZAjNXL62iNEMME1pR+XSBbdF6GGJP/zMLkrIUQa3ZPSfKnB9APQslugVi/WiM45zlS27BASv8biI+2oRDjn2chZumbskjyP8gYWyeF7wHRu/4WYuWAOTVpI16H6CkmMcZ6QmxSIq+j1QJI4pzbet9ryj8B3JVkOwfCMjNAKkMGxPCSzNMQ1olVjLEZELJcDDGhsStEW9EHYu8CQKzvvQMxqc8X98fHP7Mwp6NhaOhpz3X1g+3lPebLMyJQg+U9Jv79IRRDLu7vLnMaoqfQQyF8awglky2FL4QQIvO1hx4QEzquQJh+H8ijHfnKwv2lIa5Wwsk70sgL6SsgGto9EI1VS0txTdLQQijkNRCNyA0pzzeldF+H8hKavhAvUxGEWfQgxOzagVLe55mF3wcLyyJwfzlOhILfPCivdeRSvI4QCqwI4mU9BGCoQjqTzeXHxC8AwlT7s5TGbekZfAsx29ZFKd9maYRBjKVugDDzlUgycg1ixvJEVN3hqkqeYHnHM/mXo3DvEbi/6YoO9w8xeB9Ah2q8f6GS3F6U0rkBMZ5cZYcsa/Vpj4xbycNKk7JW2ZHK7P7yjmg3pHrvoiQvlvJjKazkNwHCDHsXol34ChZ2hjJ5Bkch3oViiI/yP1qqJ9SwjYBQ9n+X5LNCutcuKO+GNhAK76KZrFWRJwthp0l1fAHiHbkJMeQwB8rtgxZiyWIGxPtUDtFufgexOY+HWfhhEAr3tsnzr9IeNOSPSRkliCaBZNLLBvAl53xyg2aGIIgmAY3REgRBEIQDIUVLEARBEA6EFC1BEARBOBAaoyUIgiAIB0I9WoIgCIJwIKRoCYIgCMKBkKIlCIIgCAdCipYgCIIgHEiNtmDMyMiIUKvViSqV6jHOub1nEhIEQRDEQwNjrMhoNH5vMBhSYmNjcyyGq+6s44yMjAiNRrM+MDDQ29vbu9TV1VVv4eBogiAIgngo4ZxDp9NpiouLtXl5ecV6vf4JS8q22qZjtVqdGBgY6B0YGHjTzc2NlCxBEATR5GCMwc3NTR8YGHgzMDDQW61WJ1oKW21Fq1KpHvP29nb0SRwEQRAE0Sjw9vYuValUFs8cr7ai5Zz7uLq6Kh0/RRAEQRBNDldXV721+Uo1mnVM5mKCIAiCENjSibS8hyAIgiAcCClagiAIgnAgpGgJgiAIwoGQoiUIgiAIB0KKtoYwxmJNf2q1OtbHx6dbXFxc+2XLlvlainfgwIHm48aNi2jVqlUXd3f37p6enjHt27fvNHXq1FbZ2dkapThbt27Vmt9P6WcaZ+zYsRGMsdglS5a0rE65SktLVVqt9hHGWGxCQkJkdeISDUd9yqO9nDp1yn3SpElh7dq166zVah/RaDTdAwICug4cOLDtwoUL/e7cufPADBKlMrRo0eKRXr16RS9ZsqSl0Wisco/MzExXOXzz5s1jioqKFNs0o9GIsLCwP8hht27dqq1N2Yi6Rak9c3V17R4aGtrliSeeiDh58qS7tfgbNmzwGjlyZGRoaGiXZs2axbi7u3dv3bpLVoN3AAAY80lEQVT1H0aPHh25bt06L6U4RqMRqamp3oMHD44KCAjoqtFount7ez/So0eP6Pfeey/AXD5rQ422YCTuM2vWrGsAoNfrWVZWltvu3bt9jh07ps3IyGi+cuXKK3I4o9GIl156KTQ5OTlIrVbzfv36lSQkJNzU6XSq9PR0j5SUlMC0tDT/L774Iue5554rUrpXSEiI7qmnnip0ZHlWr17tc/v2bTVjDDt37vS5fv36paCgIIMj70nUHfUpj9Z45ZVXghcuXBhiNBrRrVu3srFjx5Z4enoa8/PzXY4cOaKdPXt2+KpVq/x//vnnX62V4cKFC247d+70Tk9P9zxx4oRHWlraJaX7qdVqXl5erlq1apXvK6+8csPcf/PmzdorV664qdVqbjAYaNmEkyI/ewC4deuW+tSpUx4bNmxouX37dp/du3f/1rdv33LT8EVFRaonn3wycvfu3d5ubm48Li6upF27dnc1Gg2/ePGi2759+1ps2rTJd/fu3XkpKSn35P/GjRvqUaNGtTl8+LCXp6enYeDAgbfCw8N1N2/eVO/fv7/F22+/HbZq1aqALVu2ZHXp0qWituUiRVtLPv3001zT/zdt2qQdM2ZM+9WrVwe++uqr+dHR0ToAmDNnTnBycnJQSEiIbtOmTVk9evS4axovNTXVe9q0aW1eeOGFNn5+fucSEhKqbAoSGhqqM79fXbNmzRp/lUqFF154IW/58uWBycnJfvPmzctz5D2JuqM+5dESr7/+etCCBQtCgoKCdP/4xz8uDBo0qMw8zNdff91i0aJFgfaUYefOnR6PPfZYh6+++sr/zTffvN6hQwedeZzOnTvfyc3NdU1NTfVXUrQrVqzwd3V15XFxcSUHDhxoYW9ZiPpFqX2bNGlSWFpaWsCCBQsC+/btmyO7GwwGjBw5MurgwYNevXv3Lv3nP/+ZHRER8cAeD+Xl5eyTTz7xP3funLtZvDZHjhzx6t+/f8m//vWvC6adCb1ej1mzZoUuXbo0aNiwYe1Pnz79H39//1p1Nsh0XMeMGjWqNDIy8i7nHIcOHfIAhHlr0aJFwS4uLnzDhg2/mzdqADB58uTid99997LBYMCMGTNaGwz134lMT093P3PmjEdcXFzJvHnzrmk0Gp6WluZX7xkh6oz6lsfMzEzXBQsWhLi4uPDNmzdnKSlZAJgwYcKtffv2ZdmT5pAhQ8rkMhw5csRDKYyLiwt/6qmnCn/55Zfmhw8fbmbqd+3aNZddu3Z5Dx06tMjb25usM42MoUOHlgBAYWHhAx3DlJQU34MHD3q1bt26YufOnb+bK1kAaNasGX/rrbfyU1JSLstuy5cv9z1y5IhXq1atKrZt23be3GKn0Wjw+eefX3388ceLcnNzXV977bWQ2paBFK0DkA9qkBcxJycn+xkMBjZkyJDiXr16lVuKN2vWrAJ/f399Tk6O+7Zt2+p9DGnp0qX+APDss88WBgUFGeLj429lZ2e7b9++3bO+80LUHfUpj8nJyX6VlZVs2LBhRT179qyiwE1p1qyZ3SeayGXQaDQW40yfPr2AMYbk5GR/szy11Ov1LDExsUpPl3B+du3a5QUAMTExd0zdV69e7Q8ASUlJ1728vKoO4JtgKmtr1qyR4+VptVqL8ebPn58LAN98841fRUVFrYYbnNZ0/NXRi75L9mSFFpRWuPpr3XQzBre7OjEu/GZD58sWGzdu1Obk5LgzxtCvX78yADh69KgnAAwaNKjEWlyNRoO4uLjSLVu2+P7444+e5ua6q1evus6ePVvx66pDhw7liYmJ1R5Lk7lz5w7bsGFDS09PT8PEiROLAGDSpEk3du7c6Z2cnOw/bNiw2zVNu9HzQXAM9HeqfpRqmhvx12unGiBHduNIeVRCTjs+Pr7O9kP//vvvPXNyctw1Gg0fMGCAYg8ZADp16qSLi4sr3bhxo+/t27cve3p6cgBYu3atX3h4eMXw4cNLV6xY0WQtNPrcXE3O089ER3z9f79pgoMrGzo/Spi2byUlJapTp055nDp1yjM+Pv7W22+/fV320+v1OHPmjAcAPPbYY3bLmmm84cOHW5X/2NjYu/7+/vqCggLN4cOHm8XHx9+xFt4aTqlovzp60fe9rf8Jr6g0qgAgv7TC9b2t/wkHAGdTtrJgmE4+4ZxjypQpee3bt9cBQH5+vgYAwsPDq4wtmdOqVSsdAOTm5laZ8Zmbm+u6cOHCYKV4gwcPLq6Nol29erVvSUmJesKECTfkBmr8+PG3ZsyYUbl9+3afgoKCS7Udp2i0KClZa+4NSH3KoxJy2q1bt7aZtj1lkCdDcc7xzjvvXAkPD7e6z/pzzz1XkJiY2CY1NdU3KSmpcPv27Z7Z2dnuf/3rX6+oVE73uOqV/IWLgivz8tzyFy4MCf3oI8VJZQ2NUvsWFRV198knn7zp4+Nzr/eZn5/votfrGQC0adPGblkzjRcVFWUzXnBwsK6goEBz6dIlVwAPl6JdsicrVFayMhWVRtWSPVmhzqZoZcFgjEGr1RpiY2NLJ02adGP69On38mluurOGtbA9e/a8ffz48cy6yfmDrFmzxg8ApkyZcs+8ptFoMGbMmMKVK1cGLl++vOXcuXPzHXFvou5wtDxmZma6Ll++vEqvUJ7EUp20bZVBhjGGhQsX5sycOdPmjPtnn322eM6cOZVr1qzxS0pKKkxOTvZ3cXHhL774okNn6zs7+txcTemOHX7gHKXbd/jpZ83KdcZeLec8Q/67pKRElZGR4f7GG2+0evHFFyN/+eWXZp999tlVQMyarwnVjSfLc3l5ea2+0pxS0RaUVrhWx70hMRUMSwQEBOizs7Pdc3JybOb/6tWrrgAQHBxcbycknTx50v3kyZOekZGRdwcPHvyAaS4xMfHGypUrA9PS0vxI0To/jpbHrKwsN6Veh6xoAwMD9dnZ2e5SD6BGyGUoKSlR7d2712P69OkRr776anhkZKRu5MiRVs2E7u7ufNy4cYUrV64M3L17t8f27dt9Bg8eXBwaGup0SqU+yV+4KJhLSoYbjXDmXq2Ml5eXMT4+/s7WrVvPt27duuuyZcsCZ86cmd+2bVt9UFBQpUaj4Xq9nmVnZ7t27tzZriU4pvHOnz/vamvpzvXr110BIDAwsFby45S2FH+tm2KX3pK7sxMXF3cbAPbu3au4cFqmsrISR48e1QLAgAED6m1MVJ4ElZ2d7W6+aLxXr16dASArK6vZrl27FGd8Eo2L2sjjiBEjSjnnGeY/hbRrPZnPy8vLOHr06NINGzb8bjQaWWJiYmRpaanNNmv69Ok3AGDixIlRFRUVTX4S1L3ebGWlMDNUVrLS7Tv89NeuOWVHyxw/Pz9DZGTkXYPBwI4ePeoBCGtbt27dygDg+++/t1vWNBoNunbtWgYA27Ztsyr/J0+edC8oKNCoVCrExcXV2GwMOKminTG43VU3F9UDfXw3F5VxxuB2VxsqT7UhMTHxhlqtxs6dO71PnDhhcYeTxYsX+xUUFGgiIiLuDh8+vM4mk1ijvLycrV+/3lelUmH8+PE3nnzyySq//v37lwDA8uXL/W2lRzg/jpTHadOm3XBxceE7duzwycjIsLqbT3l5uV325d69e5c/9dRTBXl5eZr3338/wFb4mJiYu7Gxsbfz8vI0ISEhulGjRlmd9PKwY9qblZF7tQ2UpWpz69YtF+BB0+/zzz9fAACff/55kK0PMFNZk+N98cUXgbdv37Yog/PmzQsGgL59+5YE19LM7pSKdmJc+M23RnS6GKB10zEAAVo33VsjOl10tvFZe+nUqZMuKSnpWmVlJRszZkxbpQZo7dq13nPnzg1Tq9VYvHjxJbVaXS95+/LLL32Ki4td+vfvf2vdunUXv/nmmyq/zZs3n2/WrJnxu+++8yksLKyfjDkTmubKAzuW3J0cR8pjdHS07n//939z9Xo9GzlyZLsDBw40Vwr37bffesXHx7ezN8/vv//+NTc3N75s2bKggoICm5lZsWLFxbS0tPPr1q37vb7eJWel7NAh73u9WZnKSlZ28JB3A2WpWqxdu9b76tWrri4uLjw+Pv6epS8xMfFm//79Sy5evOg2dOjQqIsXL1aZsHf37l32t7/9zX/q1KlhpvH69OlTeunSJbcRI0ZEmctTZWUlXn755ZAtW7b4uru7Gz/66KMr5ulWF6c1HUyMC7/ZWBWrEgsWLMgtKytTrVy5MrB3796dBgwYUNKhQ4dyvV7P0tPTPc+ePevh7u5uTElJuWBpHMra8h4AePPNN/P8/PwemBn85Zdf+u3fv1/RtPLMM8/clNeiPf/88xbNaz4+Psbhw4cX/fvf/26ZkpLi+8YbbxTYV+qHBCdfwlMT6kIeLfHhhx9er6ysZAsXLgx59NFHO8bExJR169atTN6C8dixY9qLFy+6de7c2W5zXGRkpP7pp58uWLNmTcC8efOCli5datW6FRMTczcmJsbqOt6mQvvDh842dB7sxbR9KysrU2VmZrrLO3m98cYbV8PCwu71LNVqNTZv3nx+/PjxkXv27PGOjo7u0qdPn5L27dvfVavV/PLly66HDh3yKioqcklMTLy3u52Liws2b958fuTIkW3279/fIioqqkt8fPyt1q1b39uC8erVq66urq48JSUlu0+fPhbXmtsLk2dV2cuZM2dyunXr1qTHPACxCTZg3+QTU3744Yfmn332WcCxY8e0hYWFGpVKxUNDQ3Xx8fG3XnvttbyoqKgqk6C2bt2qTUhIaG8r7d9+++0neYu9sWPHRqxfv97qgQJz58698v7777fy9fWtzM3NPevm5mZRGHbt2uUxZMiQDtHR0eW//fbbf+wpK1F/1Kc82svJkyfdFy9e7H/48GGva9euuVZUVDBvb+/Kjh07lo8ZM6Zo2rRphaYbCdgqw+XLl13at2/fBQDOnTv3U1hYWGVmZqZrhw4dunTv3v12RkaGzRn5o0aNity8ebPvli1bzo0YMaJehmcI25gfigIIRerj46Pv2rVrWVJSUsGYMWMsDgGsX7/ea82aNS1PnjzpWVhYqOGcw9/fXx8TE3N70qRJhePGjasSVzpUwOerr75qefbsWY/i4mIXeQe0Rx55pCwtLS27Ovscnzlzxq9bt24RiuUjRUsQBEEQwPHjx5sNHjw42t3d3bh3797MulK0TjlGSxAEQRD1Ta9evcpTU1MvFBUVuQwZMqT9uXPn6mRJKSlagiAIgpAYO3ZsSWpq6oUJEyYU7tq1q072eXfayVAEQRAE0RBMnDixeOLEicV1lR71aAmCIAjCgZCiJQiCIAgHQoqWIAiCIBwIKVqCIAiCcCCkaAmCIAjCgZCiJQiCIAgHQoqWIAiCIBwIKVqCIAiCcCCkaAmCIAjCgZCiJQiCIAgHQoqWIAiCIBwIKdoawhiLNf2p1epYHx+fbnFxce2XLVvmaynegQMHmo8bNy6iVatWXdzd3bt7enrGtG/fvtPUqVNbZWdna+qzDMTDQ13KY3R0dKcXX3wx9NKlS4p7oW/dulXLGIvt1atXtONKRBAPD3SoQC2ZNWvWNQDQ6/UsKyvLbffu3T7Hjh3TZmRkNF+5cuUVOZzRaMRLL70UmpycHKRWq3m/fv1KEhISbup0OlV6erpHSkpKYFpamv8XX3yR89xzzxU1XImIxkxdyWNycnJQWlpawMqVKy9MmDDhVsOViCAaP3Twew1hjMUCAOc8w9R906ZN2jFjxrQHgF9//fWn6OhoHQC88sorwQsWLAgJCQnRbdq0KatHjx53TeOlpqZ6T5s2rU1lZSXbsGHDuYSEhNL6KgvR+HGUPBoMBuzYsSNz0KBBZbLf1q1btQkJCe179ux5+/jx45mOLx1BOD908Hs9MmrUqNLIyMi7nHMcOnTIAwAyMzNdFy1aFOzi4sI3bNjwu3mjBgCTJ08ufvfddy8bDAbMmDGjtcFgsOt+mZmZroyx2LFjx0b88ssvbsOGDWvj7e39iIeHR0y/fv3apaenuwNAbm6uy4QJE8L9/f27urm5df/DH/7QccuWLVrz9MaOHRvBGIvNzMyscuCxbDKcPXt2SLUrhmgQaiuPlZWVbPr06eH1n3OCeHhwWkVr5Ear/zszspWAMQYASE5O9jMYDGzIkCHFvXr1KrcUb9asWQX+/v76nJwc923btlVRgta4fPmyW//+/TsUFBRoxo8ff6N///4lR44c8RoyZEj0Tz/95Na7d+8Op0+f9khISCgaPnz4zczMzGbjxo1rl5WVVUWhEg8XtZXHzMzMZnv27PGon9wSxMOHU47RfpT+UUiprlQ9v+/8yyqmgpEb8c7hd8K0rlrDnJ5zchs6f9bYuHGjNicnx50xhn79+pUBwNGjRz0BYNCgQSXW4mo0GsTFxZVu2bLF98cff/Ssjvk4PT3dc86cOVf//ve/X5fdXn311eBPPvkkpH///h0ff/zxm2vXrr2kVqsBAEuXLi1JSkqK/PDDDwNXrVp1uUaFbSKsfvXHruWl+ioT1ZppNfrnPx5wtiHyZC91JY8//PCD5+DBg8ushScIQhmnU7RGbkSprlS98feNAQAwv+/8y+8cfids4+8bA0a3HZ1v5EaomPN0xGUzqunkE845pkyZkte+fXsdAOTn52sAIDw8XGcrvVatWukAIDc3t1ozkENCQnQffPDBdVO3xMTEwk8++SREp9OxpUuXXpGVLABMnTr15syZMyN+/vnnZtW5T1MkrKNv8e8Z+X5GA2eym0rNeFgn3+KGzJcSjpLHK1eukOWDIGqI0ylaFVNhft/5lwFg4+8bA2SFO7rt6Hy5h+tMLFy4MBgQZjmtVmuIjY0tnTRp0o3p06fflMOYm+6sYR42MzPTdfny5X7m4T799NMHevadOnW64+Ly4OOMiIjQSdcKHx+fB2zvLi4uaNmyZeX169epAbVB3Kioa7+fzH/gGTAGxI2KcjrriqPksaKiwrlePIJoRDidogXuK1tZyQKiZ+tsShaoOstTiYCAAH12drZ7Tk6OTaV29epVVwAIDg7WA0BWVpab3HiaYq5ovby8qsye0mhEp1ir1SrOrHJxceGVlZW2W9smjralu75t94Abcq9WpWY8KjbghtbXvbKh82aOo+TRz8/P6cpKEI0F59NcwL0xWVO3dw6/E9aYJkSZEhcXdxsA9u7d62UtXGVlJY4ePaoFgAEDBtwGgBEjRpRyzjPMf47Kq0olREKv11dRwEVFReoqEZoIcaOirkGqEWftzdpLTeSxR48eND5LEDXE6RStrGTlMdkzfzmTMbrt6PyNv28MaKzKNjEx8YZarcbOnTu9T5w44W4p3OLFi/0KCgo0ERERd4cPH94g62i9vb0rASA7O7tKbyc9Pb15/efIOZB7tWCAs/Zm7eXFF1+8oVarub3y2KJFC8MTTzxhdeIUQRCWcTpFq2IqaF21BtMx2fl9518e3XZ0vtZVa3BG87EtOnXqpEtKSrpWWVnJxowZ0zYjI6NK47Z27VrvuXPnhqnVaixevPiS6cSl+qRXr15lAJCSkvLAmOTx48ebrVq1KrBBMuUkxI2KuubRwq2iMfdmASA6Olo3c+ZMu+QRAN5+++0rWq228X3hEoST4JRjtHN6zsk1nV0sK9vGqGRlFixYkFtWVqZauXJlYO/evTsNGDCgpEOHDuV6vZ6lp6d7nj171sPd3d2YkpJyYeTIkQ22K9QzzzxT/MEHH1Rs3brVt0ePHq6xsbFlV65ccd29e7f3n/70p+Jt27b5NFTeGhptS3f95A/7/dzQ+agLPv7442tlZWXq5cuXW5RHAJg2bdr12bNnK+4Ed+HCBfexY8dGKPmFhYXpFi1a1Kg/SAiirnBKRQugyhKexqxkAUCtVmPFihVXnn766ZufffZZwLFjx7RHjhzxUqlUPDQ0VPfCCy/kvfbaa3lRUVH6hsxn8+bN+e7duzNnzpwZdvDgQa+ff/7Zo23btuUpKSkXWrZsaWjKivZhQqVSITk5+cqECRNuLlmyJODo0aPaw4cPe+l0OgYA/v7++pSUlOzRo0db/OgrLCx0Wb9+fUslv+jo6HIApGgJArTXMUEQJhQVFani4uI6nD9/3n3NmjUXnn32WadbK0wQzgjtdUwQhF34+PgYv/vuuyxfX9/KKVOmtPn222+tzkwmCMI2pGgJgniAtm3b6jdt2pSVlJR0/fTp083u3r1La60JohY47RgtQRANR58+fcr79Olj8cABgiDsh3q0BEEQBOFASNESBEEQhAMhRUsQBEEQDqRGira6S4IIgiAI4mHFlk6stqJljBXpdLpqnZVKEARBEA8rOp1OwxgrsuRfbUVrNBq/Ly4u1tYuWwRBEATxcFBcXKw1Go3fW/KvtqI1GAwpeXl5xXl5eb4VFRUaMiMTBEEQTQ3OOSoqKjR5eXm+eXl5xQaDIcVS2GpvwQgAGRkZEWq1OlGlUj3GOae9bwmCIIgmB2OsyGg0fm8wGFJiY2NzLIajHilBEARBOA5a3kMQBEEQDoQULUEQBEE4EFK0BEEQBOFASNESBEEQhAMhRUsQBEEQDuT/AfMJ2fOc38XMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "idx=0\n",
    "seeds = [1,2,3,4,5]\n",
    "ess = []\n",
    "wis = []\n",
    "for seed in seeds:\n",
    "    for var_coeff in [0.0, 0.1, 1.0, 10.0]:\n",
    "        for threshold in [0.4, 0.6, 0.8, 1.0]:\n",
    "            param = {\"name\": \"POELA\", \"env\": \"mimic45\", \"thresh\": threshold, \"var_coeff\": var_coeff, \"seed\": seed}\n",
    "            model, val_score, train_score = get_model_from_param(param)\n",
    "            wis.append(val_score[0])\n",
    "            ess.append(val_score[-1])\n",
    "plt.scatter(ess,wis,marker='o',color='C0',label=\"POELA\")\n",
    "ess = []\n",
    "wis = []\n",
    "for seed in seeds:\n",
    "    for var_coeff in [0.0, 0.1, 1.0, 10.0]:\n",
    "        for threshold in [0.01, 0.02, 0.05, 0.1]:\n",
    "            param = {\"name\": \"POELA-kNN\", \"env\": \"mimic45\", \"thresh\": threshold, \"var_coeff\": var_coeff, \"seed\": seed}\n",
    "            model, val_score, train_score = get_model_from_param(param)\n",
    "            wis.append(val_score[0])\n",
    "            ess.append(val_score[-1])\n",
    "plt.scatter(ess,wis,marker='x',color='C2',label=\"PO-mu\")\n",
    "ess = []\n",
    "wis = []\n",
    "for seed in seeds:\n",
    "    for var_coeff in [0.0, 0.1, 1.0, 10.0]:\n",
    "        param = {\"name\": \"norm-POEM\", \"env\": \"mimic45\", \"thresh\": 0.01, \"var_coeff\": var_coeff, \"seed\": seed}\n",
    "        model, val_score, train_score = get_model_from_param(param)\n",
    "        wis.append(val_score[0])\n",
    "        ess.append(val_score[-1])\n",
    "plt.scatter(ess,wis,marker='s',color='C1',label=\"PO-CRM\") \n",
    "ess = []\n",
    "wis = []\n",
    "for seed in seeds:\n",
    "    for threshold in [0.0, 0.01, 0.05, 0.1, 0.3]:\n",
    "        param = {\"name\": \"PQL\", \"env\": \"mimic45\", \"thresh\": threshold, \"seed\": seed}\n",
    "        model, val_score, train_score = get_model_from_param(param)\n",
    "        wis.append(val_score[0])\n",
    "        ess.append(val_score[-1])\n",
    "plt.scatter(ess,wis,marker='v',color='C4',label=\"PQL\") \n",
    "ess = []\n",
    "wis = []\n",
    "for seed in seeds:\n",
    "    for threshold in [0.0, 0.01, 0.05, 0.1, 0.3]:\n",
    "        param = {\"name\": \"BCQ\", \"env\": \"mimic45\", \"thresh\": threshold, \"seed\": seed}\n",
    "        model, val_score, train_score = get_model_from_param(param)\n",
    "        wis.append(val_score[0])\n",
    "        ess.append(val_score[-1])\n",
    "plt.scatter(ess,wis,marker='^',color='C3',label=\"BCQ\")   \n",
    "    \n",
    "# plt.hlines(v_pib, 0, 3000, color=\"grey\", label=\"Behavior Policy\")\n",
    "# plt.hlines(ci_pib[0], 0, len(results)-1, linestyle='--', color=\"grey\")\n",
    "# plt.hlines(ci_pib[1], 0, len(results)-1, linestyle='--', color=\"grey\")\n",
    "# plt.xscale(\"log\")\n",
    "# plt.ylim([75, 101])\n",
    "plt.xlabel(\"Effective Sample Size on validation set\")\n",
    "plt.ylabel(\"WIS on validation set\")\n",
    "plt.grid(True)\n",
    "#plt.title(\"WIS-ESS trade-off for different algorithms\")\n",
    "plt.legend(bbox_to_anchor=[1.2, -0.4], loc='right', ncol=3)\n",
    "plt.savefig(\"new_plots/mimic_ess_tradeoff.png\", bbox_inches='tight', pad_inches=0.1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MIMIC Visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize heatmap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_actions(dataset, model, mask=None):\n",
    "    n = dataset[\"observations\"].shape[0]\n",
    "    horizon = dataset[\"observations\"].shape[1]\n",
    "    na = dataset[\"pibs\"].shape[-1]\n",
    "    \n",
    "    action_counts = np.zeros(na)\n",
    "    \n",
    "    for idx in range(n):\n",
    "        for t in range(horizon):\n",
    "            if (mask is None) or (mask[idx, t] == True):\n",
    "                if model is None:\n",
    "                    a = dataset[\"actions\"][idx,t,0]\n",
    "                    action_counts[a] += 1\n",
    "                else:\n",
    "                    pie = model.get_prob(dataset[\"observations\"][idx,t,:], dataset[\"estm_pibs\"][idx,t,:], \n",
    "                                         dataset[\"nn_action_dist\"][idx,t,:])\n",
    "                    action_counts += pie\n",
    "            \n",
    "            if not dataset[\"not_done\"][idx, t]:\n",
    "                break\n",
    "    return action_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.palplot(sns.color_palette(\"Greens\",12))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(font_scale=1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cmap = sns.color_palette(\"Purples\", n_colors=100) #“Purples” rocket_r\n",
    "action_counts = count_actions(test_dataset, None).reshape((5,5))\n",
    "ax = sns.heatmap(action_counts, cmap=cmap, xticklabels='auto', yticklabels='auto', vmin=0, vmax=12000)\n",
    "ax.invert_yaxis()\n",
    "ax.tick_params(left=False, bottom=False)\n",
    "plt.xlabel(\"Vasopressors\") \n",
    "plt.ylabel(\"IV fluid\") \n",
    "plt.title(\"Physician Policy\")\n",
    "plt.savefig(\"plots/physician_heatmap.png\", bbox_inches='tight')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "for param in [ours_rel_param, ours_novar_rel_param, knn_rel_param, knn_novar_rel_param, \n",
    "              npoem_rel_param, npoem_novar_param, nomask_rel_param, nomask_novar_param, \n",
    "              bcq_rel_param, pql_rel_param]:\n",
    "    model, _, _ = get_model_from_param(param)\n",
    "    action_counts = count_actions(test_dataset, model).reshape((5,5))\n",
    "    #print(action_counts.round(2))\n",
    "    \n",
    "    ax = sns.heatmap(action_counts, cmap=cmap, xticklabels='auto', yticklabels='auto', vmin=0, vmax=12000)\n",
    "    ax.invert_yaxis()\n",
    "    ax.tick_params(left=False, bottom=False)\n",
    "    plt.xlabel(\"Vasopressors\") \n",
    "    plt.ylabel(\"IV fluid\") \n",
    "    plt.title(toprint(param['name']))\n",
    "    plt.savefig(f\"plots/{param['name']}_heatmap.png\", bbox_inches='tight')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'test_dataset' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-4950174db2b3>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mlow_sofa\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mtest_dataset\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'observations'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0msofa_id\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0msofa_std\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0msofa_mean\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0;36m5\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'test_dataset' is not defined"
     ]
    }
   ],
   "source": [
    "sofa_id=34\n",
    "sofa_mean = 6.0695\n",
    "sofa_std = 3.3598\n",
    "\n",
    "low_sofa = (test_dataset['observations'][:,:,sofa_id]*sofa_std+sofa_mean) < 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### For different SOFA value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sofa_id=34\n",
    "sofa_mean = 6.0695\n",
    "sofa_std = 3.3598"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "high_sofa = (test_dataset['observations'][:,:,sofa_id]*sofa_std+sofa_mean) >= 13"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mid_sofa = ((test_dataset['observations'][:,:,sofa_id]*sofa_std)+sofa_mean < 15) & ((test_dataset['observations'][:,:,sofa_id]*sofa_std)+sofa_mean >= 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "low_sofa = (test_dataset['observations'][:,:,sofa_id]*sofa_std+sofa_mean) < 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "action_counts = count_actions(test_dataset, None, high_sofa).reshape((5,5))\n",
    "cmap = sns.color_palette(\"rocket_r\",n_colors=100) #“Purples” rocket_r\n",
    "ax = sns.heatmap(action_counts, cmap=cmap, xticklabels='auto', yticklabels='auto', vmin=0, vmax=800)\n",
    "ax.invert_yaxis()\n",
    "ax.tick_params(left=False, bottom=False)\n",
    "plt.xlabel(\"Vasopressors\") \n",
    "plt.ylabel(\"IV fluid\") \n",
    "plt.title(\"Clinician Policy\")\n",
    "plt.savefig(\"plots/physician_highsofa_heatmap.png\", bbox_inches='tight')\n",
    "plt.show()\n",
    "\n",
    "for param in [ours_rel_param, ours_novar_rel_param, knn_rel_param, knn_novar_rel_param, \n",
    "              npoem_rel_param, npoem_novar_param, nomask_rel_param, nomask_novar_param, \n",
    "              bcq_rel_param, pql_rel_param]:\n",
    "    model, _, _ = get_model_from_param(param)\n",
    "    action_counts = count_actions(test_dataset, model, high_sofa).reshape((5,5))\n",
    "    #print(action_counts.round(2))\n",
    "    \n",
    "    ax = sns.heatmap(action_counts, cmap=cmap, xticklabels='auto', yticklabels='auto', vmin=0, vmax=800)\n",
    "    ax.invert_yaxis()\n",
    "    ax.tick_params(left=False, bottom=False)\n",
    "    plt.xlabel(\"Vasopressors\") \n",
    "    plt.ylabel(\"IV fluid\") \n",
    "    plt.title(toprint(param['name']))\n",
    "    plt.savefig(f\"plots/{param['name']}_highsofa_heatmap.png\", bbox_inches='tight')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "action_counts = count_actions(test_dataset, None, mid_sofa).reshape((5,5))\n",
    "ax = sns.heatmap(action_counts, cmap=\"Blues\", xticklabels='auto', yticklabels='auto')\n",
    "ax.invert_yaxis()\n",
    "ax.tick_params(left=False, bottom=False)\n",
    "plt.xlabel(\"Vasopressors\") \n",
    "plt.ylabel(\"IV fluid\") \n",
    "plt.title(\"Physician Policy\")\n",
    "plt.savefig(\"plots/physician_midsofa_heatmap.png\", bbox_inches='tight')\n",
    "plt.show()\n",
    "\n",
    "for param in [ours_rel_param, ours_novar_rel_param, knn_rel_param, knn_novar_rel_param, \n",
    "              npoem_rel_param, npoem_novar_param, nomask_rel_param, nomask_novar_param, \n",
    "              bcq_rel_param, pql_rel_param]:\n",
    "    model, _, _ = get_model_from_param(param)\n",
    "    action_counts = count_actions(test_dataset, model, mid_sofa).reshape((5,5))\n",
    "    #print(action_counts.round(2))\n",
    "    \n",
    "    ax = sns.heatmap(action_counts, cmap=\"Blues\", xticklabels='auto', yticklabels='auto')\n",
    "    ax.invert_yaxis()\n",
    "    ax.tick_params(left=False, bottom=False)\n",
    "    plt.xlabel(\"Vasopressors\") \n",
    "    plt.ylabel(\"IV fluid\") \n",
    "    plt.title(toprint(param['name']))\n",
    "    plt.savefig(f\"plots/{param['name']}_midsofa_heatmap.png\", bbox_inches='tight')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "action_counts = count_actions(test_dataset, None, low_sofa).reshape((5,5))\n",
    "ax = sns.heatmap(action_counts, cmap=\"Greens\", xticklabels='auto', yticklabels='auto')\n",
    "ax.invert_yaxis()\n",
    "ax.tick_params(left=False, bottom=False)\n",
    "plt.xlabel(\"Vasopressors\") \n",
    "plt.ylabel(\"IV fluid\") \n",
    "plt.title(\"Physician Policy\")\n",
    "plt.savefig(\"plots/physician_lowsofa_heatmap.png\", bbox_inches='tight')\n",
    "plt.show()\n",
    "\n",
    "for param in [ours_rel_param, ours_novar_rel_param, knn_rel_param, knn_novar_rel_param, \n",
    "              npoem_rel_param, npoem_novar_param, nomask_rel_param, nomask_novar_param, \n",
    "              bcq_rel_param, pql_rel_param]:\n",
    "    model, _, _ = get_model_from_param(param)\n",
    "    action_counts = count_actions(test_dataset, model, low_sofa).reshape((5,5))\n",
    "    #print(action_counts.round(2))\n",
    "    \n",
    "    ax = sns.heatmap(action_counts, cmap=\"Greens\", xticklabels='auto', yticklabels='auto')\n",
    "    ax.invert_yaxis()\n",
    "    ax.tick_params(left=False, bottom=False)\n",
    "    plt.xlabel(\"Vasopressors\") \n",
    "    plt.ylabel(\"IV fluid\") \n",
    "    plt.title(toprint(param['name']))\n",
    "    plt.savefig(f\"plots/{param['name']}_lowsofa_heatmap.png\", bbox_inches='tight')\n",
    "    plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}