plots-appendix-heldout-1k.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Appendix - Repeated evaluation w/original setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T16:43:04.067095Z",
     "start_time": "2019-04-19T16:43:04.049687Z"
    }
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T16:43:04.753547Z",
     "start_time": "2019-04-19T16:43:04.068986Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cf.counterfactual as cf\n",
    "import cf.utils as utils\n",
    "import pandas as pd\n",
    "import pickle\n",
    "import itertools as it\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "from scipy.linalg import block_diag\n",
    "\n",
    "# Sepsis Simulator code\n",
    "from sepsisSimDiabetes.State import State\n",
    "from sepsisSimDiabetes.Action import Action\n",
    "from sepsisSimDiabetes.DataGenerator import DataGenerator\n",
    "import sepsisSimDiabetes.MDP as simulator \n",
    "\n",
    "import mdptoolboxSrc.mdp as mdptools\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)\n",
    "\n",
    "# Avoid Type 3 fonts\n",
    "import matplotlib\n",
    "matplotlib.rcParams['pdf.fonttype'] = 42\n",
    "matplotlib.rcParams['ps.fonttype'] = 42\n",
    "\n",
    "figpath = \"./figs\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T16:43:04.753547Z",
     "start_time": "2019-04-19T16:43:04.068986Z"
    }
   },
   "outputs": [],
   "source": [
    "fig_prefix = \"appendix-multiple-heldout-1k\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T16:43:04.769134Z",
     "start_time": "2019-04-19T16:43:04.755230Z"
    }
   },
   "outputs": [],
   "source": [
    "SEED = 1  # Note this is not the only random seed, see the loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T16:43:04.787255Z",
     "start_time": "2019-04-19T16:43:04.770642Z"
    }
   },
   "outputs": [],
   "source": [
    "np.random.seed(SEED)\n",
    "NSIMSAMPS = 1000  # Samples to draw from the simulator\n",
    "NSTEPS = 20  # Max length of each trajectory\n",
    "NCFSAMPS = 1  # Counterfactual Samples per observed sample\n",
    "DISCOUNT_Pol = 0.99 # Used for computing optimal policies\n",
    "DISCOUNT = 1 # Used for computing actual reward\n",
    "PHYS_EPSILON = 0.05 # Used for sampling using physician pol as eps greedy\n",
    "\n",
    "PROB_DIAB = 0.2\n",
    "\n",
    "# Number of iterations to get error bars\n",
    "N_REPEAT_SAMPLING = 100\n",
    "NHELDOUT = 1000 # Heldout samples for WIS\n",
    "\n",
    "# These are properties of the simulator, do not change\n",
    "n_actions = Action.NUM_ACTIONS_TOTAL\n",
    "n_components = 2\n",
    "\n",
    "# These are added as absorbing states\n",
    "n_states_abs = State.NUM_OBS_STATES + 2\n",
    "discStateIdx = n_states_abs - 1\n",
    "deadStateIdx = n_states_abs - 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T16:43:04.944514Z",
     "start_time": "2019-04-19T16:43:04.789044Z"
    }
   },
   "outputs": [],
   "source": [
    "# Get the transition and reward matrix from file\n",
    "with open(\"./data/diab_txr_mats-replication.pkl\", \"rb\") as f:\n",
    "    mdict = pickle.load(f)\n",
    "\n",
    "tx_mat = mdict[\"tx_mat\"]\n",
    "r_mat = mdict[\"r_mat\"]\n",
    "p_mixture = np.array([1 - PROB_DIAB, PROB_DIAB])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T16:43:05.064330Z",
     "start_time": "2019-04-19T16:43:04.946096Z"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.linalg import block_diag\n",
    "\n",
    "tx_mat_full = np.zeros((n_actions, State.NUM_FULL_STATES, State.NUM_FULL_STATES))\n",
    "r_mat_full = np.zeros((n_actions, State.NUM_FULL_STATES, State.NUM_FULL_STATES))\n",
    "\n",
    "for a in range(n_actions):\n",
    "    tx_mat_full[a, ...] = block_diag(tx_mat[0, a, ...], tx_mat[1, a,...])\n",
    "    r_mat_full[a, ...] = block_diag(r_mat[0, a, ...], r_mat[1, a, ...])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T16:43:14.717995Z",
     "start_time": "2019-04-19T16:43:05.066210Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 7s, sys: 220 ms, total: 1min 7s\n",
      "Wall time: 4.82 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "fullMDP = cf.MatrixMDP(tx_mat_full, r_mat_full)\n",
    "fullPol = fullMDP.policyIteration(discount=DISCOUNT_Pol, eval_type=1)\n",
    "\n",
    "physPolSoft = np.copy(fullPol)\n",
    "physPolSoft[physPolSoft == 1] = 1 - PHYS_EPSILON\n",
    "physPolSoft[physPolSoft == 0] = PHYS_EPSILON / (n_actions - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T16:43:14.761240Z",
     "start_time": "2019-04-19T16:43:14.740751Z"
    }
   },
   "outputs": [],
   "source": [
    "obs_reward = []\n",
    "offpol_opt_reward_WIS_hard_train = []                         \n",
    "offpol_opt_reward_WIS_hard_ho = []                         \n",
    "offpol_opt_reward_mb = []\n",
    "true_rl_reward = []\n",
    "\n",
    "# We will save the detailed samples from the first run\n",
    "saved_material = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_rl_policy(rl_policy, obs_samps, proj_lookup):\n",
    "\n",
    "    passes = True\n",
    "    # Check the observed actions for each state\n",
    "    obs_pol = np.zeros_like(rl_policy)\n",
    "    for eps_idx in range(NSIMSAMPS):\n",
    "        for time_idx in range(NSTEPS):\n",
    "            this_obs_action = int(obs_samps[eps_idx, time_idx, 1])\n",
    "            # Need to get projected state\n",
    "            if this_obs_action == -1:\n",
    "                continue\n",
    "            this_obs_state = proj_lookup[int(obs_samps[eps_idx, time_idx, 2])]\n",
    "            obs_pol[this_obs_state, this_obs_action] += 1\n",
    "\n",
    "    # Check if each RL action conforms to an observed action\n",
    "    for eps_idx in range(NSIMSAMPS):\n",
    "        for time_idx in range(NSTEPS):\n",
    "            this_full_state_unobserved = int(obs_samps[eps_idx, time_idx, 1])\n",
    "            this_obs_state = proj_lookup[this_full_state_unobserved]\n",
    "            this_obs_action = int(obs_samps[eps_idx, time_idx, 1])\n",
    "\n",
    "            if this_obs_action == -1:\n",
    "                continue\n",
    "            # This is key: In some of these trajectories, you die or get discharge.  \n",
    "            # In this case, no action is taken because the sequence has terminated, so there's nothing to compare the RL action to\n",
    "            true_death_states = r_mat[0, 0, 0, :] == -1\n",
    "            true_disch_states = r_mat[0, 0, 0, :] == 1\n",
    "            if np.logical_or(true_death_states, true_disch_states)[this_full_state_unobserved]:\n",
    "                continue\n",
    "\n",
    "            this_rl_action = rl_policy[proj_lookup[this_obs_state]].argmax()\n",
    "            if obs_pol[this_obs_state, this_rl_action] == 0:\n",
    "                print(\"Eps: {} \\t RL Action {} in State {} never observed\".format(\n",
    "                    int(time_idx / NSTEPS), this_rl_action, this_obs_state))\n",
    "                passes = False\n",
    "    return passes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the projection matrix for obs->proj states\n",
    "n_proj_states = int((n_states_abs - 2) / 5) + 2\n",
    "proj_matrix = np.zeros((n_states_abs, n_proj_states))\n",
    "for i in range(n_states_abs - 2):\n",
    "    this_state = State(state_idx = i, idx_type='obs', \n",
    "                       diabetic_idx = 1)  # Diab a req argument, no difference\n",
    "    # assert this_state == State(state_idx = i, idx_type = 'obs', diabetic_idx = 0)\n",
    "    j = this_state.get_state_idx('proj_obs')\n",
    "    proj_matrix[i, j] = 1\n",
    "\n",
    "# Add the projection to death and discharge\n",
    "proj_matrix[deadStateIdx, -2] = 1\n",
    "proj_matrix[discStateIdx, -1] = 1\n",
    "\n",
    "proj_matrix = proj_matrix.astype(int)\n",
    "\n",
    "proj_lookup = proj_matrix.argmax(axis=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T17:12:19.338614Z",
     "start_time": "2019-04-19T16:43:14.787770Z"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "81f56b497e114f22a84f94ee42e27d2d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, description='Outer Loop', style=ProgressStyle(description_width='initial')…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "for it in tqdm(range(N_REPEAT_SAMPLING), desc=\"Outer Loop\"):\n",
    "    np.random.seed(it)\n",
    "    dgen = DataGenerator()\n",
    "    states, actions, lengths, rewards, diab, emp_tx_totals, emp_r_totals = dgen.simulate(\n",
    "       NSIMSAMPS, NSTEPS, policy=physPolSoft, policy_idx_type='full', \n",
    "       p_diabetes=PROB_DIAB, use_tqdm=False) #True, tqdm_desc='Behaviour Policy Simulation')\n",
    "\n",
    "    obs_samps = utils.format_dgen_samps(\n",
    "       states, actions, rewards, diab, NSTEPS, NSIMSAMPS)\n",
    "    \n",
    "    emp_tx_mat = np.copy(emp_tx_totals)\n",
    "    emp_r_mat = np.copy(emp_r_totals)\n",
    "\n",
    "    ############## Construct the Transition Matrix w/proj states ##############\n",
    "    proj_tx_cts = np.zeros((n_actions, n_proj_states, n_proj_states))\n",
    "    proj_tx_mat = np.zeros_like(proj_tx_cts)\n",
    "\n",
    "    # (1) NOTE: Previous code marginalized here, but now we are just getting observed quantities out, no components\n",
    "    est_tx_cts = np.copy(emp_tx_totals)\n",
    "    assert est_tx_cts.ndim == 3\n",
    "\n",
    "    # (2) Add new aborbing states, and a new est_tx_mat with Absorbing states\n",
    "    death_states = (emp_r_mat.sum(axis=0).sum(axis=0) < 0)\n",
    "    disch_states = (emp_r_mat.sum(axis=0).sum(axis=0) > 0)\n",
    "\n",
    "    est_tx_cts_abs = np.zeros((n_actions, n_states_abs, n_states_abs))\n",
    "    est_tx_cts_abs[:, :-2, :-2] = np.copy(est_tx_cts)\n",
    "\n",
    "    death_states = np.concatenate([death_states, np.array([True, False])])\n",
    "    disch_states = np.concatenate([disch_states, np.array([False, True])])\n",
    "    assert est_tx_cts_abs[:, death_states, :].sum() == 0\n",
    "    assert est_tx_cts_abs[:, disch_states, :].sum() == 0\n",
    "\n",
    "    est_tx_cts_abs[:, death_states, deadStateIdx] = 1\n",
    "    est_tx_cts_abs[:, disch_states, discStateIdx] = 1\n",
    "\n",
    "    # (3) Project the new est_tx_cts_abs to the reduced state space\n",
    "    for a in range(n_actions):\n",
    "        proj_tx_cts[a] = proj_matrix.T.dot(est_tx_cts_abs[a]).dot(proj_matrix)\n",
    "\n",
    "    # Normalize\n",
    "    nonzero_idx = proj_tx_cts.sum(axis=-1) != 0\n",
    "    proj_tx_mat[nonzero_idx] = proj_tx_cts[nonzero_idx]\n",
    "\n",
    "    proj_tx_mat[nonzero_idx] /= proj_tx_mat[nonzero_idx].sum(axis=-1, keepdims=True)\n",
    "\n",
    "    ############ Construct the reward matrix, which is just based on absorbing states \n",
    "    proj_r_mat = np.zeros((n_actions, n_proj_states, n_proj_states))\n",
    "    proj_r_mat[..., -2] = -1\n",
    "    proj_r_mat[..., -1] = 1\n",
    "\n",
    "    proj_r_mat[..., -2, -2] = 0 # No reward once in aborbing state\n",
    "    proj_r_mat[..., -1, -1] = 0\n",
    "\n",
    "    ############ Construct the empirical prior on the initial state ##################\n",
    "    initial_state_arr = np.copy(states[:, 0, 0])\n",
    "    initial_state_counts = np.zeros((n_states_abs,1))\n",
    "    for i in range(initial_state_arr.shape[0]):\n",
    "        initial_state_counts[initial_state_arr[i]] += 1\n",
    "\n",
    "    # Project initial state counts to new states\n",
    "    proj_state_counts = proj_matrix.T.dot(initial_state_counts).T\n",
    "    proj_p_initial_state = proj_state_counts / proj_state_counts.sum()\n",
    "\n",
    "    # Because some SA pairs are never observed, assume they cause instant death\n",
    "    zero_sa_pairs = proj_tx_mat.sum(axis=-1) == 0\n",
    "    proj_tx_mat[zero_sa_pairs, -2] = 1  # Always insta-death if you take a never-taken action\n",
    "\n",
    "    # Construct an extra axis for the mixture component, of which there is only one\n",
    "    projMDP = cf.MatrixMDP(proj_tx_mat, proj_r_mat, \n",
    "                           p_initial_state=proj_p_initial_state)\n",
    "    try:\n",
    "        RlPol = projMDP.policyIteration(discount=DISCOUNT_Pol)\n",
    "    except:\n",
    "        assert np.allclose(proj_tx_mat.sum(axis=-1), 1)\n",
    "        RlPol = projMDP.policyIteration(discount=DISCOUNT_Pol, skip_check=True)\n",
    "\n",
    "    # Estimate the observed policy\n",
    "    obs_pol_proj = proj_tx_cts.sum(axis=-1)  # Sum over the \"to\" state\n",
    "    obs_pol_proj = obs_pol_proj.T # Switch from (a, s) to (s, a)\n",
    "    obs_states = obs_pol_proj.sum(axis=-1) > 0 # Observed \"from\" states\n",
    "\n",
    "    obs_pol_proj[obs_states] /= obs_pol_proj[obs_states].sum(axis=-1, keepdims=True)\n",
    "    \n",
    "    # Check if we always observe the RL policy in the non-absorbing states\n",
    "    prop_rl_obs = (obs_pol_proj[:-2, :][RlPol[:-2, :]==1] > 0).mean()\n",
    "    if prop_rl_obs < 1:\n",
    "        assert check_rl_policy(RlPol, obs_samps, proj_lookup), 'RL policy validation failed'\n",
    "\n",
    "    def projection_func(obs_state_idx):\n",
    "        if obs_state_idx == -1:\n",
    "            return -1\n",
    "        else:\n",
    "            return proj_lookup[obs_state_idx]\n",
    "\n",
    "    proj_f = np.vectorize(projection_func)\n",
    "    states_proj = proj_f(states)\n",
    "    assert states_proj.shape == states.shape\n",
    "\n",
    "    obs_samps_proj = utils.format_dgen_samps(\n",
    "       states_proj, actions, rewards, diab, NSTEPS, NSIMSAMPS)\n",
    "        \n",
    "    # Get the true RL reward as a sanity check\n",
    "    # Note that the RL policy includes actions for \"death\" and \"discharge\" absorbing states, which we ignore by taking [:-2, :]\n",
    "    NSIMSAMPS_RL = NSIMSAMPS\n",
    "    states_rl, actions_rl, lengths_rl, rewards_rl, diab_rl, _, _ = dgen.simulate(\n",
    "       NSIMSAMPS_RL, NSTEPS, policy=RlPol[:-2, :], policy_idx_type='proj_obs', \n",
    "       p_diabetes=PROB_DIAB, use_tqdm=False) #True, tqdm_desc='RL Policy Simulation')\n",
    "\n",
    "    obs_samps_rlpol = utils.format_dgen_samps(\n",
    "       states_rl, actions_rl, rewards_rl, diab_rl, NSTEPS, NSIMSAMPS_RL)\n",
    "\n",
    "    this_true_rl_reward = cf.eval_on_policy(\n",
    "        obs_samps_rlpol, discount=DISCOUNT, \n",
    "        bootstrap=False)  \n",
    "\n",
    "    # Get a soft version of the RL policy for WIS\n",
    "    RlPolSoft = np.copy(RlPol).astype(float)\n",
    "    RlPolSoft[RlPolSoft == 1] = 0.99\n",
    "    RlPolSoft[RlPolSoft == 0] = 0.01 / (n_actions - 1)\n",
    "\n",
    "    # This is the observed reward from the samples given\n",
    "    this_obs_reward = cf.eval_on_policy(\n",
    "        obs_samps_proj, discount=DISCOUNT, \n",
    "        bootstrap=False)\n",
    "\n",
    "    # This is the off-policy reward using WIS\n",
    "    this_offpol_opt_reward_WIS_hard_train, this_wis_samps, this_wis_ct = cf.eval_wis(\n",
    "        obs_samps_proj, discount=DISCOUNT,\n",
    "        bootstrap=False,\n",
    "        obs_policy=obs_pol_proj, new_policy=RlPol)\n",
    "\n",
    "    # Draw samples from the MDP under the new policy to get a model-based estimate of reward\n",
    "    BSampler = cf.BatchSampler(mdp=projMDP)\n",
    "    this_mb_samples_opt = BSampler.on_policy_sample(\n",
    "        policy=RlPol, n_steps=NSTEPS, n_samps=NSIMSAMPS_RL, \n",
    "        use_tqdm=False) #, tqdm_desc='Model-Based OPE') \n",
    "\n",
    "    this_offpol_opt_reward_mb = cf.eval_on_policy(\n",
    "        this_mb_samples_opt, discount=DISCOUNT,\n",
    "        bootstrap=False)\n",
    "\n",
    "    ###################################################\n",
    "    # Construct the held-out samples, freshly each time\n",
    "    ###################################################\n",
    "    ho_dgen = DataGenerator()\n",
    "    ho_states, ho_actions, ho_lengths, ho_rewards, ho_diab, ho_emp_tx_totals, ho_emp_r_totals = ho_dgen.simulate(\n",
    "       NHELDOUT, NSTEPS, policy=physPolSoft, policy_idx_type='full', \n",
    "       p_diabetes=PROB_DIAB, use_tqdm=False) #True, tqdm_desc='Behaviour Policy Simulation')\n",
    "\n",
    "    ho_obs_samps = utils.format_dgen_samps(\n",
    "       ho_states, ho_actions, ho_rewards, ho_diab, NSTEPS, NHELDOUT)\n",
    "\n",
    "    ho_emp_tx_mat = np.copy(ho_emp_tx_totals)\n",
    "    ho_emp_r_mat = np.copy(ho_emp_r_totals)\n",
    "\n",
    "    ############## Construct the Transition Matrix w/proj states ##############\n",
    "    ho_proj_tx_cts = np.zeros((n_actions, n_proj_states, n_proj_states))\n",
    "    ho_proj_tx_mat = np.zeros_like(ho_proj_tx_cts)\n",
    "\n",
    "    # (1) NOTE: Previous code marginalized here, but now we are just getting observed quantities out, no components\n",
    "    ho_est_tx_cts = np.copy(ho_emp_tx_mat)\n",
    "    assert ho_est_tx_cts.ndim == 3\n",
    "\n",
    "    # (2) Add new aborbing states, and a new est_tx_mat with Absorbing states\n",
    "    ho_death_states = (ho_emp_r_mat.sum(axis=0).sum(axis=0) < 0)\n",
    "    ho_disch_states = (ho_emp_r_mat.sum(axis=0).sum(axis=0) > 0)\n",
    "\n",
    "    ho_est_tx_cts_abs = np.zeros((n_actions, n_states_abs, n_states_abs))\n",
    "    ho_est_tx_cts_abs[:, :-2, :-2] = np.copy(ho_est_tx_cts)\n",
    "\n",
    "    ho_death_states = np.concatenate([ho_death_states, np.array([True, False])])\n",
    "    ho_disch_states = np.concatenate([ho_disch_states, np.array([False, True])])\n",
    "    assert ho_est_tx_cts_abs[:, ho_death_states, :].sum() == 0\n",
    "    assert ho_est_tx_cts_abs[:, ho_disch_states, :].sum() == 0\n",
    "\n",
    "    ho_est_tx_cts_abs[:, ho_death_states, deadStateIdx] = 1\n",
    "    ho_est_tx_cts_abs[:, ho_disch_states, discStateIdx] = 1\n",
    "\n",
    "    # (3) Project the new est_tx_cts_abs to the reduced state space\n",
    "    for a in range(n_actions):\n",
    "        ho_proj_tx_cts[a] = proj_matrix.T.dot(ho_est_tx_cts_abs[a]).dot(proj_matrix)\n",
    "\n",
    "    # Estimate the observed policy\n",
    "    ho_obs_pol_proj = ho_proj_tx_cts.sum(axis=-1)  # Sum over the \"to\" state\n",
    "    ho_obs_pol_proj = ho_obs_pol_proj.T # Switch from (a, s) to (s, a)\n",
    "    ho_obs_states = ho_obs_pol_proj.sum(axis=-1) > 0 # Observed \"from\" states\n",
    "\n",
    "    ho_obs_pol_proj[ho_obs_states] /= ho_obs_pol_proj[ho_obs_states].sum(axis=-1, keepdims=True)\n",
    "\n",
    "    def projection_func(obs_state_idx):\n",
    "        if obs_state_idx == -1:\n",
    "            return -1\n",
    "        else:\n",
    "            return proj_lookup[obs_state_idx]\n",
    "\n",
    "    proj_f = np.vectorize(projection_func)\n",
    "    ho_states_proj = proj_f(ho_states)\n",
    "    assert ho_states_proj.shape == ho_states.shape\n",
    "\n",
    "    ho_obs_samps_proj = utils.format_dgen_samps(\n",
    "       ho_states_proj, ho_actions, ho_rewards, ho_diab, NSTEPS, NHELDOUT)\n",
    "    \n",
    "    this_offpol_opt_reward_WIS_hard_ho, this_ho_wis_samps, this_ho_wis_ct = cf.eval_wis(\n",
    "        ho_obs_samps_proj, discount=DISCOUNT,\n",
    "        bootstrap=False,\n",
    "        obs_policy=ho_obs_pol_proj, new_policy=RlPol)\n",
    "\n",
    "    obs_reward.append(this_obs_reward)\n",
    "    offpol_opt_reward_WIS_hard_train.append(this_offpol_opt_reward_WIS_hard_train)                       \n",
    "    offpol_opt_reward_WIS_hard_ho.append(this_offpol_opt_reward_WIS_hard_ho)                       \n",
    "    offpol_opt_reward_mb.append(this_offpol_opt_reward_mb)\n",
    "    true_rl_reward.append(this_true_rl_reward)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T17:12:19.355539Z",
     "start_time": "2019-04-19T17:12:19.340377Z"
    }
   },
   "outputs": [],
   "source": [
    "# END OF LOOP\n",
    "def conv_to_np(this_list):\n",
    "    this_arr = np.array(this_list)[:, np.newaxis]\n",
    "    return this_arr\n",
    "\n",
    "obs_reward = conv_to_np(obs_reward)\n",
    "offpol_opt_reward_WIS_hard_train = conv_to_np(offpol_opt_reward_WIS_hard_train)\n",
    "offpol_opt_reward_WIS_hard_ho = conv_to_np(offpol_opt_reward_WIS_hard_ho)\n",
    "offpol_opt_reward_mb = conv_to_np(offpol_opt_reward_mb)\n",
    "true_rl_reward = conv_to_np(true_rl_reward)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T18:12:29.839468Z",
     "start_time": "2019-04-19T18:12:29.629901Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAE3CAYAAAA66vBzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xmc3HV9+PHXe5MgglYgC1hIIOAGz6Jt0yrVHyZCAvEA71pR1yOS2pq0UlsVLXeV1qLtUo9YryBYlLYKKUSTSDhUsAQ5FIVkLUGiINkkRM5kk33//vjOhsmwx+xmZmdn5/V8POYx+z3nPfvd3Xnv+/s5IjORJElS62hrdACSJEkaWyaAkiRJLcYEUJIkqcWYAEqSJLUYE0BJkqQWYwIoSZLUYkwAJUmSWowJoCRJUosxAZQkSWoxkxsdwHjX3t6eM2bMaHQYkiRJw7r55pt7MvPA4fYzARzGjBkzWLNmTaPDkCRJGlZE3FPNft4CliRJajEmgJIkSS3GBFCSJKnFmABKkiS1mKZIACNiWkRcGBE3RMSjEZERMaPKY9si4iMRsT4iHo+I2yLiDfWNWJIkafxqigQQ6ADeDGwBrh/hsecCZwH/BswHbgQui4hX1jJASZKkZtEsCeB1mXlwZr4SuKzagyLiIOCDwPmZ+c+ZuTozFwKrgfPrFKtUdz09PSxatIhNmzY1OhSNgtdPUqM1RQKYmX2jPPQEYC/g4or1FwO/FxFH7FFgUoN0dXVx22230dXV1ehQNApeP0mN1hQJ4B54PrAN6K5Yf0fp+XljG46053p6erjmmmsAWL16tVWkJuP1kzQeTPQE8ADgwczMivWby7ZLTaWyamQVqbl4/Zrf2rVrmT9/Pt3dlbUFqXlM9KngAqhM/vrXD35QxKnAqQCHHXZYHcJSq+jq6qr5h8Stt9662/Lq1avZsmVLTV+jo6ODxYsX1/SczagZr5/Xrv7OPvtsHnnkEc4880wuueSSRocjjcpErwBuBvaPiMqEb/+y7U+SmV/IzFmZOevAA4edT1mS1CLWrl3LvffeC8C9995rFVBNK558d3R8i4gFwL8DR2Tm+mH2fQewFJiZmd1l698JfAU4MjPvHuocs2bNyjVr1uxp2FLNnHLKKbs+gACmT59uFaKJeP2am9dP411E3JyZs4bdb4IngAcB9wIfz8yzy9avAg7OzN8b7vVMADXerF27lgULFuxa/vKXv0xHR0cDI9JIeP3GRj1u38OTb+EDvOhFL6rZ+b2Frz1VbQLYNG0AI+KNpS//sPQ8PyI2Ahsz89rSPjuApZn5HoDMfCAiPg18JCIeAn4M/CnwCuDkMX0DUo0cddRR7LXXXmzfvp3p06ebPDQZr5+k8aBpEkCePAD0Z0vP1wKzS19PKj3KfRR4GPgr4JnAXcCbM3NZfcKU6m/GjBl0d3dz9tlnD7+zxh2vX/3Vq4p2xhln7BrGB2DOnDleRzWlpkkAM3PInruD7ZOZO4HzSg9pQthnn304+uijrR41Ka9f81q8ePFuCaC3a9WsmiYBlOqpXu2F6mXdunVA83342L5Jza69vZ1nPOMZbN26lTlz5jB16tRGhySNigmgBHR3d7P2pz/msKftbHQoVdmrtxjB6fH1NzU4kur98uHK1hlSc5o2bRo7duzwnxk1NRNAqeSwp+3kY7MebnQYE9Z5a57W6BCkmpgyZQozZ860+qemZgIoSWqYZmt+Ac3ZBMPmF6pkAihJapju7m7uvPVWntnoQEagfwqtBwcYE3A8ur/RAWhcMgGUJDXUM4H3DD1Fu/bAl2iuCR80Nib6XMCSJEmqYAIoSZLUYkwAJUmSWowJoCRJUouxE4gkqWE2bNjAQ9hRoZ7uAx7esKHRYWicsQIoSZLUYqwASpIaZtq0aTzY0+MwMHX0JZL9pk1rdBgaZ6wASpIktRgTQEmSpBbjLWCJoiH6Iw9N4rw1T2t0KBPWPQ9NYl8bokvSuGAFUJIkqcVYAZQoGqI/vuM+Pjbr4UaHMmGdt+Zp7F2HhuhdXV10d3fX/Lz1tG7dOgAWL17c4EhGpqOjo+liljQwE0BJTa27u5tb7rgF9mt0JCPQVzzd8qtbGhvHSDzY6AAk1ZIJoKTmtx/0ze5rdBQTWts1thiSJhJ/oyVJklqMCaAkSVKLMQGUJElqMSaAkiRJLcYEUJIkqcWYAEqSJLUYh4GRJDXU/cCXyEaHUbVNpeepDY2ievfTXMNkamyYAEqSGqajo6PRIYzYxtJMLvvNnNngSKqzH835fVZ9mQBKkhqmGaeW64+5q6urwZFIo2cbQEmSpBZjAihJktRiTAAlSZJajAmgJElSizEBlCRJajEmgJIkSS3GBFCSJKnFmABKkiS1GBNASZKkFmMCKEmS1GJMACVJklqMCaAkSVKLMQGUJElqMSaAkiRJLcYEUJKkEejt7WXdunVs2rSp0aFIoza50QFIklRrXV1ddHd31+Xcd955Jzt37uTd7343hx9+eE3P3dHRweLFi2t6TmkgVgAlSapSb28vO3fuBGDLli309vY2OCJpdKwASpImnHpV0f7hH/6BO+64Y9fy9OnTOf300+vyWlI9mQBKamobNmyArdB2jTc06upB2JAbGh1Fw61atWq35ZUrV5oAqin5F1OSpCpFxJDLUrNomgpgREwHPg3MBQJYBfx1Zv6yimNzkE2/n5m31i5KSWNt2rRpbIyN9M3ua3QoE1rbNW1MO3Rao8NouOOOO47vfve7u5aPP/74BkYjjV5TVAAjYh/gauA5QCfwdmAmsDoi9q3yNF8Fjql4rK15sJKkCetNb3rTbstvfvObGxSJtGeaIgEE3gscCbw2M7+dmZcDJwGHAwurPMevMvPGisej9QpYkjTxLFu2bNdt34jgiiuuaHBE0ug0SwJ4EnBjZu4a1Ckz7wZ+AJzcsKgkSS1l5cqVZBatijKTFStWNDgiaXQGbQMYEVeP4DyZmcfVIJ7BPB+4fID1dwBvGmD9QN4XEX8L7ARuBM7MzOtrFJ8kqQXMnTuXq666it7eXqZMmcK8efMaHZI0KkNVANsoOlv0P54DzAZmAE8tPc8Gnl3aXk8HAFsGWL8Z2L+K4y8G/gI4HjgVmApcHRGzaxWgJGni6+zs3HULuK2tjc7OzgZHJI3OoAlgZs7OzDmZOQf4V6AXeElmHpmZx2TmkRQdKXpL2+ttoJ68VSWemfn2zPxGZl6fmRcDLwN+DZw30P4RcWpErImINRs3bhx9xJKkCaW9vZ358+cTEcyfP5+pU6c2OiRpVKptA3gu8PeZ+b/lKzPzR8BZDJJI1dAWiipgpf0ZuDI4pMx8CLgS+KNBtn8hM2dl5qwDDzxwpKeXJE1gnZ2dHH300Vb/1NSqHQdwJjBYKewBoKM24QzqDop2gJWeB/xslOcMBq4qSpI0qPb2di688MJGhyHtkWorgHcz+HArC4H1NYlmcFcAL4mII/tXRMQM4KWlbSMSEb8DvAr4UY3ikyRJahrVVgDPBi6JiJ8C/wn8BjgYeCNF55BT6hPeLv8OvB+4PCI+RlG5Oxe4F1jSv1NEHA78AjgnM88prfsgRUeV1RTt/g4HPgg8cwziliRJGneqSgAz89KI6KFIBD8CTKHo/HETcEJmfq9+IUJmPhIRr6CYCu5rFLdvv0cxFdzDZbsGMIndK5t3Aa8rPZ4B/JZi/MD3VLZplCRJagVVzwWcmauAVRHRBrQDPZk5ZpNvlub8fcMw+6ynomdwZi4DltUvMkmSpOYybBvAiNgrIjZHxEkAmdmXmQ+MZfInSZKk2hk2AczM7cAO4PH6hyNJkqR6q7YX8LcpOnxIkiSpyVXbBnA50BUR/0mRDN5HxRh6mTmSuYMlSZLUINVWAP8LOBR4PXARsBJYVfGsJtPT08OiRYvYtGlTo0ORJEljqNoK4Jy6RqGGWLJkCbfddhtLlizh9NNPb3Q4kiRpjFQ7DuC19Q5EY6unp4eVK4vC7YoVK1i4cKGTmkuS1CKqHgdQE8uSJUvo6ytG8unr67MKCPzy4Umct+ZpjQ6jKr95tGi9cfA+zTMa0y8fnsRRjQ5CkgSMIAGMiBcA76GYVm3vis2ZmcfVMjDV16pVq3ZbXrlyZUsngB0dHY0OYUS2r1sHwN4zZjY4kuodRfN9nyVpoqoqAYyIFwPXAuuBmcDtwP7AYcAGoLtO8alOImLI5VazePHiRocwIv3xdnV1NTgSSVIzqrYC+HHgv4G3U8wB/J7M/HFpft6vAefVKT5RfMh3d9c2x37605/Oli1bdluuZRLU0dHRdEmVJEmtotphYI4GLuaJsf8mwa6x/84DPlH70FRPhxxyyJDLkiRp4qq2AjgFeCQz+yJiM/C7ZdvuAl5Q88i0S70qaSeffDJbtmzhxBNPbOn2f5IktZpqK4C/oBgIGor2f++OiLaIaAPeBdxfj+BUX4cccgj77rsvCxcubHQokiRpDFVbAVwGzAa+TtEe8Ergt8BO4GlAyzf2qkc7vXpbv349AGeffXZjAxkh2xdKkrRnqh0I+qyyr1dFxEuANwD7AN/JzBX1Ca95XHPNNWzs2QSTmmhoxb6dANzykzsaHMgI7NzBhg0bTAAlSdoDo8pWMvMW4JYax9L8Jk2mbx9n06intkedt1iSpD1VVRvAiDg/IuZFxD71DqhZTZs2DWiusfTisa20PfzArkpgc4jS91qSJI1WtRXAtwF/B/RGxE3A6tLjB5m5rV7BNZNmnOHg9tsfoC/7OGDydo444ohGh1OlZzbl91qSpPGk2jaA0yLiKOAVwBzgVOCjwLaIuBG4OjPPrV+Y41+ztUnr6enh9a9/PQBbt27lzDPPZOpUb19LktQKqm4DmJlrgbXA5wEi4qXAWcBxwLFASyeA9VSPHsZ33333bsvvfOc7a1oFtKeuJEnjV9UJYEQ8FXgZT1QB/wB4FPgf4Oq6RKe62bp165DLkiRp4qoqAYyI64A/BrYDPwS+BSwCbs7MvvqFJ6jP7eVjjz32Seu6urpq/jqSJGn8qbYC+DLgMeBrwHeBazPTklETmz59Ovfee+9uy5IkqTVUOxXc0cBHgEOArwI9EXFTRPxjRJwYEfvWK0DVx5lnnrnbcrPNBiJJkkavqgQwM3+amV2Z+TpgKvBi4Bul5ysBR+dtMkcdddSuqt/06dMdWkWSpBZSbQUQgIiYQtHj96TS4xiK0Y+31D401duZZ57Jvvvua/VPkqQWU20nkNMpev8eAzyVouJ3LXAaxRiAP69bhKqbo446iuXLlzc6DEmSNMaq7QTyt8B1FIM/r87M2+oXkiRJkuqp2gRwqsO9SJIkTQzVTgXXBxAR7cBLKDqCLMvMzRGxN7DdBFGSJKk5VNUJJAqfBDYAVwBfBmaUNl9OcWtYkiRp3Orp6WHRokVs2uTgJdX2Av4I8H7gHIqhX6Js2zLg1TWOS5IkqaaWLl3K7bffztKlSxsdSsNVmwAuAM7JzI8DP67Y1g08q6ZRSZIk1VBPTw/Lly8nM1m+fHnLVwGrTQAPBW4cZNt2wJlAJEnSuLV06VIyE4C+vr6WrwJWmwD+CnjBINteCNxdm3AkSZJqb+XKlfT29gLQ29vLihUrGhxRY1WbAF4GnBERLy1blxFxFPA3wKU1j0ySJKlG5s6dy5QpUwCYMmUK8+bNa3BEjVVtAngWcCfFYNDrSusuA35SWj6/5pFJkiTVSGdnJxFFH9a2tjY6OzsbHFFjVZUAZuZjwGzgncAPgVXATcCpwNzM3F6n+CRJkvZYe3s78+fPJyKYP38+U6dObXRIDVXtTCBk5k7ga6WHJElSU+ns7GT9+vUtX/2D6m8BDyoifj8ivlWLYCRJkuqlvb2dCy+8sOWrfzBMBTAiJgF/CBwG/CIzbynbNgs4E3gl8FA9g5QkSVLtDFoBjIhpwI+AG4BvAmsi4hsRsVdEfLG07RXABcCRYxGsJEnSaDkV3BOGugV8PvAc4O8pqnzvB/4E+AHwbuAi4MjM/LvM3FzvQCVJkvbEkiVLuO2221iyZEmjQ2m4oRLA44CzMvPjmfmdzPwc0ElxS/jCzHxXZv5mTKKUJEnaAz09PaxcuRKAFStWtHwVcKgE8ECePP3bDaXny+oTjiRJUu0tWbKEvr4+oJgKrtWrgEMlgG0U8/yW619+tD7hSJIk1d6qVat2W+6vBraq4cYBfE1ElM8B3AYkcFJEvKh8x8z8cq2DkyRJqoX+WUAGW241wyWAHx1k/RkVywmYAEqSpD3W1dVFd3d3Tc/59Kc/nS1btuy2vHjx4pq+RkdHR83PWS9D3QI+YgSPug8DExHTI+I/I2JrRPw2Iv47Ig6r8ti9I+KTEXFfRDwWETdExLH1jlmSJI0PhxxyyJDLrWbQCmBm3jOWgQwlIvYBrga2UfRETuA8YHVEHJ2Zjwxzii8BrwL+Fvg/4C+B70bEMZl5a/0ilyRJI1WvKtrJJ5/Mli1bOPHEEzn99NPr8hrNouq5gBvsvRRVxmdnZjdARNwOrAMWAp8a7MCIeCHwVuDdmfmV0rprgTuAc4CT6hu6JEkaDw455BC2b9/OwoULGx1Kw+3xXMBj5CTgxv7kDyAz76YYlPrkKo7tBb5RduwO4FLghIh4Su3DlSRJ482UKVOYOXOmcwHTPAng84GfDrD+DuB5VRx7d2ZWDl1zB7AX0LHn4UmSJDWPZkkADwC2DLB+M7D/Hhzbv12SJKllNEsCCEXHj0rVDOITIz02Ik6NiDURsWbjxo3VxidJktQURpQARkRbRLwgIl4eEfvWK6gBbGHgSt3+DFzdK7d5iGP7t+8mM7+QmbMyc9aBBx44okAlSZLGu6oTwIj4S+B+4DaKIVmeXVr/7Yio96iHd1C05av0POBnVRx7RGkomcpjtwO1HWlSkiRpnKsqAYyI9wL/Cnwb+FN2v316PfCG2oe2myuAl0TErgGnI2IG8NLStuGOnQK8qezYyRTvY0Vmbqt1sJIkSeNZtRXA04ALMvNU4FsV2+6kVA2so38H1gOXR8TJEXEScDlwL7Ckf6eIODwidkTErqnqSgM9fwP4l4hYEBHHUQwBcwRwZp3jliRJGneqTQCPAL47yLZHgP1qE87ASjN9vAJYC3wNuAS4G3hFZj5ctmsAk3jy+3oX8BWK2UOuBKYDJ2bmj+sZtyRJ0nhU7UwgPcCMQbY9G/hVTaIZQmb+kmFuNWfmegbo3ZuZj1FUMU+rS3CSJElNpNoK4DLgjPI2eEBGRDvwAYq2gZIkSWoC1SaAHwO2UczGsYpiXL0u4OfAToo5dSVJktQEqkoAM3MTMAv4BEWP2l9Q3D7+N+CYzNxatwglSZJUU9W2ASQzHwLOLT0kSZLUpJppKjhJkiTVQFUVwIi4eojNfcBW4GbgS5n5m1oEJkmSpPqo9hZwAEcBv0sx/t5vgIMpxge8r7T8SuADEfHyzBxuejZJkiQ1SLW3gD8FPA78YWY+KzP/JDOfBfxRaf3ZwExgI/APdYlUkiRJNVFtAngecFZm3lK+MjNvpkj+zsvMDcAngWNrG6IkSZJqqdoE8CiK2UAGshHoKH39C2DfPQ1KkiRJ9VNtG8D1wAJg+QDbTi1tB2gHNu1xVJI0Eg9C2zVNNKhB/wzmT2toFCPzIHBoo4OQVCvVJoDnABdHxO3AfwEPAAdRzM37AuCtpf2OB35U6yAlaTAdHR3D7zTOrFu3DoCZh85scCQjcGhzfq8lDayqBDAz/yMieija+51OMRtIL7AGmJeZq0q7nkYxNZwkjYnFixc3OoQR64+5q6urwZFIalUjmQlkJbAyItoobvX2ZGZfxT6P1zg+SZI0TnV1ddHd3d3oMKrWX31vtn8cOzo6ah5z1Qlgv1LS90BNo5AkSU2nu7ubO37yc/bb56BGh1KVvu0BwK9+0TzdFR58tD4pV9UJYETsBcwHng3sXbE5M9M5giVJajH77XMQc57zlkaHMWGtvvPSupy32qngDgG+D8wAkmJmEEpf9zMBlCRJagLVjpvwSYrx/g6jSP5eDBxJMetHd+lrSZIkNYFqbwH/P+CDwK9Ly32ZuR44IyImAV3AybUPT5IkSbVWbQVwKvDrUgeQR4D9y7ZdDcyucVySJEmqk2oTwA0UQ79AMd3bvLJtfww4/IskSVKTqPYW8Grg5cC3gSXAZyLiRRSDQZ9QWidJkqQmUG0C+DHgAIDM/FxETAb+FNgH+CeKqeIkSZLUBKpNAHuBe/oXMvNC4MK6RCRJkqS6GjYBLFX7NgGvA5bVPSJJktQUNmzYwNZHH6rbYMUqZgLJDY/V/LzDdgLJzB3Ab4CdNX91SZIkjblqbwFfDCwArqpjLJIkqYlMmzaN2LbJqeDqaPWdl3LotKk1P2+1CeB64K0RcRNwOXAfu08DR2Z+ubahSZIkqR6qTQA/U3o+FPjDAbYnYAIoSZLUBKpNAI+oaxSSJEkaM1UlgJl5z/B7SZIkqRlUWwEEICKOBo6lmBt4SWbeHxEdwG8y86F6BChJkqTaqioBjIinUPQEfj0QFG3+lgH3U8wEshb4cJ1ilCRJUg0NOw5gyT8AxwNvBw6mSAL7LaeYD1iSJElNoNpbwH8GfCwzvx4Rkyq23Q3MqGlUkiRJqptqK4BTgZ8PcY6n1CYcSZIk1Vu1CeDdwDGDbPtj4K7ahCNJkqR6qzYBvAj4cEScAuxVWpcRMQf4AA4CLUmS1DSqbQP4T8ALga8BXyyt+z6wN3BpZl5Yh9gkSdI49+CjD7D6zksbHUZVHn58CwBP23v/BkdSvQcffYBDadBcwJm5E3hLRHyGosfvQcAm4DuZeW3No5IkSeNeR0dHo0MYkXXrNgNw6LNqn1DVy6FMrcv3eUQDQWfm9cD1NY9CkiQ1ncWLFzc6hBHpj7erq6vBkTReVW0AI+LHEfHXEXFwvQOSJElSfVXbCeQ3FO0A742IqyLiLRGxdx3jkjSE3t5e1q1bx6ZNmxodiiSpCVXbBnB+RBwEvBV4G/B14KGI+E/g4sxcXccYpabV1dVFd3d3zc971113sWPHDhYsWMD06dNrfv6Ojo6mu7UjSapetRVAMvOBzPyXzJwFPB/4DHAcsCoi7qlXgJJ219vby44dOwDYvHkzvb29DY5IktRsRtQJpF9m/jwizgHuAM4HptU0KmmCqEcV7YILLmDt2rX09vYyefJkZs6cyWmnnVbz15EkTVxVVwD7RcQrIuIrFO0CLwI2AItqHZikga1cuXJX1a+3t5cVK1Y0OCJJUrOpthfwCyLi/Ij4JbASeDnwr8BzMvOYzPxsPYOU9IS5c+cyZcoUAKZMmcK8efMaHJEkqdlUWwG8HVgIfAeYnZlHZuYZmbmufqFJGkhnZycRAUBbWxudnZ0NjkiS1GyqTQD/FHhmZp5aGgx6TEVEW0R8JCLWR8TjEXFbRLyhymO/GhE5wONf6h23VA/t7e3Mnz+fiGD+/PlMndo8I9pLksaHaoeBuWywbRHxcqAzM99ds6ie7Fzgg8BHgZuBtwCXRcSrM/OqKo7fCJxUse6+2oYojZ3Ozk7Wr19v9U+SNCqj6gUcER3AO4C3A4cDjwJ1SQBL4w9+EDg/M/+5tHp1KYbzgWoSwO2ZeWM94pMaob29nQsvvLDRYUiSmlTVvYAj4hkRcWpEfB+4i6IatwV4H3BIneIDOAHYC7i4Yv3FwO9FxBF1fG1pXOrp6WHRokXOBCJJGpUhE8BS27tXRsSlFLdMPw/MoBgEGuCvM3NJZv62jjE+H9gGVE6ncEfp+XlVnOOgiOiJiB0RsTYiPhQRk2oapTSGli5dyu23387SpUsbHYokqQkNmgBGxD8DvwKWAa8BvgWcCBwGnAHEWAQIHAA8mJlZsX5z2fah3Ar8DfBminaA1wKfAJbUMkhprPT09LB8+XIyk+XLl1sFlCSN2FAVwNOAgyja2B2Wmadk5orM7AMqk7GqRcTxg/TKrXxc03/IIK9XVQJamr7uwsy8OjOvysz3Uoxh+J6ImDlIjKdGxJqIWLNx48bRvE2pbpYuXUr//0N9fX1WASVJIzZUAvhl4CHgVcBdEfFvEfHHNXjNHwLPreLxjtL+m4H9o3/gsyfsX7Z9pP6j9DxroI2Z+YXMnJWZsw488MBRnF6qH2cCkSTtqUETwMxcADwTeBvF0Ct/DtwQET8HPsQoq4CZ+Whm3lnF45elQ+4AngI8q+JU/W3/fjaKMPqTyVFXMqVGcSYQSdKeGrITSGY+nplfz8wTgOnA6cBO4MMUSdT5EfG2iNi7jjF+B9gOnFKx/m3ATzPz7lGc860Uyd9NexibNOacCUSStKeqHgYmM+/LzH/MzBcALwY+C8wELqKOgypn5gPAp4GPRMRpETE7Ij4HvIIiId0lIr4XEd1ly4dHxHUR8RcRMS8iXhMRXwYWAUsy8xf1iluqF2cCkSTtqVENBJ2ZNwE3RcQHKHoIv2OYQ/bUR4GHgb+iuC19F/DmzFxWsd8kdn9PD1G0EfwQcDBF1e/nwGKKBFZqSs4EIknaE6NKAPtlZi/w36VH3WTmTuC80mOo/WZXLG8GXlu/yKTGcCYQSRq53t5e1q9fz6ZNm1r+7knVt4AlSZKa2f33388jjzzi8FmYAEqSpBbQ09PD5s3FyHEOor+Ht4AlSZJqrauri+7uyhlg98y99967axD9bdu2sWDBAqZPn17T1+jo6GDx4sU1PWe9WAGUJEkT3pYtW4ZcbjVWACVJ0rhSjyraBRdcwJVXXsmOHTuYPHkyr371qznttNNq/jrNwgqgJEma8Do7O+nr6wOKedRbfRgtE0BJkqQWYwIoSZImvKVLl9LWVqQ9bW1tLT8UjAmgJEma8FauXMmOHTsA2LFjBytWrGhwRI1lAihJkia8uXPnEhEARATz5s1rcESNZQIoSZImvNe85jW7xgHMTE466aQGR9RYJoCSJGnCW7Zs2W4VwCuuuKLBETWWCaAkSZrwVq5cuVsF0DaAkiRJE9zcuXOZMmUKAFOmTLENYKMDkCRJqrfOzs5dt4Db2tocCLrRAUiSJNVbe3s78+fPJyKYP38+U6dObXRIDeVcwJIkqSV0dnayfv36lq/+gQmgJElqEe3t7Vx44YWNDmNc8BawJElSizEBlCRJajEmgJItu4EQAAARl0lEQVQkSS3GBFCSJKnFmABKkiS1GBNASZKkFmMCKEmS1GJMACVJklqMCaAkSVKLMQGUJElqMSaAkiRJLcYEUJIkqcWYAEqSJLUYE0BJkqQWYwIoSZLUYkwAJUmSWowJoCRJUosxAZQkSWoxJoCSJEktxgRQkiSpxZgASpIktRgTQEmSpBZjAihJktRiTAAlSZJajAmgJI2x3t5e1q1bx6ZNmxodiqQWNbnRAUjSeNXV1UV3d3fNz3vXXXexY8cOFixYwPTp02t67o6ODhYvXlzTc0qaeKwAStIY6u3tZceOHQBs3ryZ3t7eBkckqRVZAZSkQdSjknbBBRewdu1aent7mTx5MjNnzuS0006r+etI0lCsAErSGFq5cuWuql9vby8rVqxocESSWpEJoCSNoblz5zJlyhQApkyZwrx58xockaRWZAIoSWOos7OTiACgra2Nzs7OBkckqRU1RQIYEadFxLKIuC8iMiLOGuHxL4uIH0bEYxFxf0R8KiKeWqdwJWlQ7e3tzJkzB4A5c+YwderUBkckqRU1RQIIvBc4CPj2SA+MiKOBlcADwKuBjwHvAr5aw/ikMdXT08OiRYscR65Jbdu2bbdnSRprzZIAPj8zXwwsGsWxZwMbgDdl5vcy84vAXwFvjog/qGWQ0lhZunQpt99+O0uXLm10KBqhnp4errvuOgCuvfZak3hJDdEUCWBm9o3muIiYApwIfDMzywfb+iawHTi5BuFJY6qnp4fly5eTmSxfvtwEosksWbKEvr7iT1pfXx9LlixpcESSWlFTJIB74FnA3sBPy1dm5uPAL4DnNSIoaU8sXbqUzASKBMIqYHNZtWrVbssrV65sUCSSWtlETwAPKD1vGWDb5rLtu4mIUyNiTUSs2bhxY92Ck0bDceSaW38P4MGWJWksjHkCGBHHl3ryDve4phYvV3rOIbY9SWZ+ITNnZeasAw88sAZhSLXjOHLN7bjjjttt+fjjj29QJJJaWSMqgD8EnlvF4x01eK3NpeeBKn37l22XmobjyDW3hQsX0tZW/Olta2tj4cKFDY5IUisa8wQwMx/NzDurePyyBi/3C2Ab8PzylRGxN3Ak8LMavIY0ptrb25k/fz4Rwfz58x1Hrsm0t7czd+5cAObNm+f1k9QQE7oNYGZuB75DMeTL5LJNbwSeAlzRkMCkPdTZ2cnRRx9t9a9JLVy4kBe+8IVW/yQ1TPT3JhzPImIWMIMiYf0GcBnFUC4AV2Xmo6X9vgR0ZubksmNfBNwALAc+UzrPJ4HvZeabhnvtWbNm5Zo1a2r2XiRJkuolIm7OzFnD7Td5uB3GifcD5aWON5UeAEcA60tfTyo9dsnMWyPiBOAfgSuBrcBFwOl1jFeSJGncaooKYCNZAZQkSc2i2grghG4DKEmSpCczAZQkSWoxJoCSJEktxgRQkiSpxZgASpIktRh7AQ8jIjYC9zQ6jjpqB3oaHYRGxWvX3Lx+zc3r17wm+rU7PDMPHG4nE8AWFxFrqukurvHHa9fcvH7NzevXvLx2BW8BS5IktRgTQEmSpBZjAqgvNDoAjZrXrrl5/Zqb1695ee2wDaAkSVLLsQIoSZLUYkwAJUmSWowJ4AQTEfMiYnlEbIqIxyNibUT8Y0TsX7FfRsR5jYpzvIqIPyt9b46tWH9waf1vBjjmL0vbXlBaPqu0PLlsn9+JiLMj4mcR8UhEbImIn0TEkog4qMrYlkXEhWXLs0uvVfPf4/73MMpjL4+Iz9Qojma6Hu8svU7H6N/xbuefXTrf7Cr2zYg4qxavO8j5XxsRpw2w/vcj4tGIOKxerz2Wyq5hRsRRA2yfXbb9+NK6s8rWZUTsiIh7IuJLEXHo2L+Liavi+zzYY32j42wWJoATSEScDnwXeBxYAJwAfB54J3BTRExvXHRN49rS87EV648FHgUOiojnDLBtE3DHQCeMiEnAKuB9wJeAk4BO4D+APwEOGS6oUgI0Fzi/bPVs4Ezq83v8ReCYUR57FvDegT5AR6GZrsdE9lrgSQlgZt4CrATOHfOI6ush4O0DrH9HadtAXkbxOzMH+DjwKuDKevyD1sKOqXjcT/GZV77udQ2LrslMHn4XNYOImAOcB/xLZn6gbNO1EfEt4GbgIoo/ThpEZv46Iv6PgROOq4Hnlr6+s2zb/wOuz8F7VL0c+CPgtZl5edn6K4CPV/kB8bfAssz8VRX7PklETAF2DBHjbjJzA7BhNK+VmbdExK3AXwN/MZpzlJ1rQl6PCWYJcHlEfCQzf93oYGrkv4G3RcQZ/T9HEfFU4A3Af1H8U13pR5m5o/T19RGxE/h34NnAz+sf8sSXmTeWL0fENqCncv1gIuIpmbmtLsE1If8zmTj+DtgMfKRyQ2beTVGpmB0RLy7bFBHx0YjYEBGPRcR1EfEidt/hhIj4YURsjYiHI+KuiDijru+k8a4Fjim/ZUiRZFwPfJ+yZCQiZgK/C1w3xPkOKD3fP9DGzOwbKpiIOASYD3y9bN1ZFNU/gN7+2x+lbTNKy38REf8UEb8GtgH7RcSBpduca0u37u6NiK9X3qoa6BZw6ZznRcTiiLg7Ih6KiGsj4vkDhH0pcErpQ3NPjfvrUaE9Ii6JiN9GxK8joisi9q44xz5RNM24OyK2l54/OlzyGRGTStfgvtL1u2aQ7z8RcWJE3FD63d4aEd+OiGdX7LM+Ir46wLG7bimXtncChw5ym20F8FsGToqa1deAwymqev1eB0yiSACr8dvS85QaxqUqRcSlEdEdEcdGxI0R8RhwTkTsXfoZ/nDF/s8prX9LxfrjS79nD5ceV0bEc8f0zdSJCeAEUPpgfDmwMjMfH2S3K0rPryhb9w7glcD7Kf54Hwx8LyIOKJ33yNJxdwN/SnGr7FPAvjV+C+PNdcDTgD8AiIj9gBdQJBzXs3s16tiyYwbzY2AHsCQiXhcV7TGrMJfig+f7Zeu+SHH7Ep649VR5y/ajwFHAqRQfXo9TJD+PU/yjcCJFJWsm8IPKJGUQb6O4tfVXwLuAwyiqP5V3E64DfmeAmEajGa5Hua8BvwBeD3wO+EvK/jErfa++S9FM418pkskvAn8PfHKY1z4LOB24hOK27Aqe+N3eJSJOBK4EHqb43X0fxffs+5XJfhXOBa4CNjLAbbZS1esGip+nieIeip+h8tvA7wC+RfE9HcikiJgcEU+NiD+kuE53AD+ta6QaSjvF7+NFFL9n/zmSgyPi9RS/qz3AWyl+Hg4ErouI361tqA2QmT6a/EGRuCXwiSH22bu0z2dLy0nxQ71v2T4zgF7g3NLyG0v7/U6j3+MYfz+PLL3vD5aWX0PR3mwvioQqgRmlbUuBrcCksuPPKu0zuWzdAooPjgT6KD4YPgkcUkU8nwN+NcD6J71O2XVMikQnhjn3JGB6af/XVZ67Yt8E1gFTytb1/4z8ScW+U4CdwOktdD3eWTrf2RXr/wdYW7b89tJ+x1bs91FgO3BQaXl2ab/ZpeX9SzF/vuK4D5X2O6ts3ZrStSp/z0dQ/H5/qmzdeuCrA7yXyvN9FdgwxPfkXIp/LNrq/ftZz0fZNewA3g1sofjb+bsU/zTMLbsux1f8fFU+fg48q9HvaSI/Sj+/Fw+y7dLSdTihYn3/Z+GHK9Y/p7T+LaXlNuBe4KqK/Q4AHgTOb/T739OHFcCJIUZ53FWZ+Uj/QmauB27kiarNrRQfGJdGxBujyt6RzS4z/4+i/Vt/NelYivY92zNzLfBAxbYfZObOYc75RYpE620Uo9C3AR8E7hjsFl6ZQyiqLyP17Sz9xSoXEe+LiNsi4mGKD7VfljY9u3LfAazMzN6y5Z+UnnfrBVraZytVdKgYThNejysrln/C7t+fEykqTD8sVYwml6qCKygS55cMct7fo6i+f7Ni/aXlCxGxL0W19Bv5RJs0smgK8gOKuwW1thF4Ck/cXp8ILqN4T68BTqFoMvC9IfZ/CUXb0hcDbwYeAVZExMF1jlODezQzvzvKY58PTAMurvg9/S1wE09ul9x0TAAnhh7gMYrKz2D6t91btu5JQ2iU1h0KkJndFD2J2yjK6PdHxI8ioh4fIOPNdcDLIiJ4or1Zv+8Dx0bENIrv61C3G3fJzC2ZeUlm/nlmPpfiFt7vAGcPc+jeFG34Ruq+yhURsQj4LEUv2NcDf8wTCUc1t4A3Vyz3xzXQsY8BtWgDCM11PQb6Hj2lbPkgivZlvRWP/y1tnzrIeftvOVX+3lYu70/xT+GTrj9FElOPJO2x0nOtrnfDZeZDwLcpKrbvAC7JoduH3pyZazLzfzPzMoqmEkcwQO9pjZkB2/lWqb/gcQlP/l09nsF/T5uGvYAngMzcERHXAXMjYu8cuB3gSaXnq8vWDfSf6cHArp6NmbkaWB0RTwFeCpxDMbTBjMzsqc07GJeuo2jz8RKKasrHyrZdT9G7tT8RvpZRyMzLI+I24HnD7LqJ4oNkxC8xwLq3AN/LzL/pXxERozl3NQ6g+OekFibC9Sg//m6KKtFA1g+yvj+hO5jdh7ip/D3eQnHtnznAOZ5Zev1+j1PcSt+lvw3wCPUfM9H+JlxEUdFtA/5sJAdm5m8iogc4uh6BqSoD/Q3spWieslfF+sqErv/35G8Y+J/KwdrbNw0rgBPHJyl+gD9euaH0Af8h4LrM/FHZpleWbhf17zeD4gP2hspzZOa2zLwa+CeK21D1ShrGi/4k4sMU1ZTy78n3KTpOvJmiLdqaoU4UEe0DdbAofe+nM3ClptydwPQBOlr0V6FGUnXZh+IPYLl3jeD4qkTEMykqZXfV6JTNcD2q9Z3S6zxcqhhVPgZLom6nuK1YmTju1mux1KzjZuBNUYx5CEBEHE4xzmF5gnwPReeQcq8e4LW3MfTP2RHAvZn52BD7NKOVFLfcP5+ZA44rOZhSJ4F2Rtd8Q3VSah7yK578c/+qiuWfAL8GnjvI72nTd+6xAjhBZOb3ohie5ZxSIncRRSXgDyg+NLfy5IFNH6Noo/JJiltUZ1O0b/g0QET8OcXttqsobh23U/Rm/DUTvGdbZt4ZEQ9QtP+5OTPLe/7dQtEY/zXA6oo2cQOZDXyuNJzG9RQNiA8HFlFUTj41zPHXUVyboyk6dvT7Wen5byJiObAzM4dMfiiSjw9FMWj4/1L0Cn/jMMeMRv9wQ1Xdjh1Ok1yPal1CkXR/LyIuAG6jqEY8i6JS/9rMfLTyoMx8MCI+DXw0Ih6iaDP4R8B7BniNv6eoXP1PRHyWohf12RR/By4o2+9S4Mul8/4P8EIGHs7lZ8ABEfE+igT78cz8Sdn2F1Ojaz2elJKFait/L45i7L82ip+nv6WoNH2+TuFp9C4FTouID1H8PM8B3lS+Q2bujIj3A5dFxD4Uw/9soqiiv5SiY9e/jW3YtWUCOIFk5rkRcRPwAeArFNWeX1Ikg5/IzMq2SRdRVBT+jSK5u4miB1T/frdRdJ3/BEV7iM0U1ZZTJuB/+gO5jiI5Km9v1v+H4QaKHoHVfOjdSDHMxysoxlPbnyLpuAmYW6qsDuV6iqT7NeyecPwPRXu+vwDOoKiMDdch6BxgP4qfkb0pqkEnAP9XxfsYiVdTJGrdNTzneL8eVcnM3og4geIfs1MpqmePUAwdcyVFT+DBnEVxjRdQDN/0o1Icu1WnMvM7EfEqirEiv1k65zXA3+XugzUvpahGvgdYWHpvrwMqr9sXKe4OfJzi5+ceSu2Ko5hh6IUUSWcr6x8WKCnant0M/Hlm/u/gh6hBzgaeTvF3cB9gGcU/PrsN7ZSZ34pikoXTKYbdeirFHYIbgIvHMN66iAE6CUoaZ0qD8p4CHDVQz97xpHR79T6KYVu+NNz+zaiZrke9laoo76MY8mTI3teSxg/bAErN4dMUlZc3NDqQKiykGJplaaMDqaNmuh51U0r2/wo4w+RPai4mgFITyMz+NpyVPdfGo23AO8vHoJtomux61NMMitlMvtbgOCSNkLeAJUmSWowVQEmSpBZjAihJktRiTAAlSZJajAmgJElSizEBlCRJajH/H/8iJtNLPrZ9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams.update({'font.size': 16})\n",
    "reward = np.concatenate([obs_reward, \n",
    "                         offpol_opt_reward_WIS_hard_train,\n",
    "                         offpol_opt_reward_WIS_hard_ho,                         \n",
    "                         offpol_opt_reward_mb,\n",
    "                         true_rl_reward,\n",
    "                         ], axis=1)\n",
    "reward_df = pd.DataFrame(reward, columns=['Obs', \n",
    "                                          'WIS (train)',\n",
    "                                          'WIS (heldout)',                                          \n",
    "                                          'MB',\n",
    "                                          'True'\n",
    "                                         ])\n",
    "\n",
    "plt.figure(figsize=(10,5))\n",
    "sns.boxplot(data=reward_df, whis=[2.5, 97.5])\n",
    "plt.ylabel(\"Average Reward\")\n",
    "plt.savefig(\"{}/{}-ope_wis_mb_cf_true.pdf\".format(\n",
    "    figpath, fig_prefix), bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-04-19T18:13:58.113270Z",
     "start_time": "2019-04-19T18:13:58.086468Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RESULTS:\n",
      "\tObserved Reward:  0.3097 \t 95% Range: 0.2694 to 0.3510\n",
      "\tTrue RL Reward:\t -0.2743 \t 95% Range: -0.5850 to 0.0549\n",
      "\tWIS (train) :\t 0.6074 \t 95% Range: -0.4217 to 0.9903\n",
      "\tWIS (heldout) :\t 0.3153 \t 95% Range: -0.9243 to 0.9879\n",
      "\tMB Estimate:\t 0.8128 \t 95% Range: 0.5700 to 0.9625\n"
     ]
    }
   ],
   "source": [
    "print((\"RESULTS:\"\n",
    "       \"\\n\\tObserved Reward:  {:.4f} \\t 95% Range: {:.4f} to {:.4f}\"\n",
    "       \"\\n\\tTrue RL Reward:\\t {:.4f} \\t 95% Range: {:.4f} to {:.4f}\"\n",
    "       \"\\n\\tWIS (train) :\\t {:.4f} \\t 95% Range: {:.4f} to {:.4f}\"\n",
    "       \"\\n\\tWIS (heldout) :\\t {:.4f} \\t 95% Range: {:.4f} to {:.4f}\"\n",
    "       \"\\n\\tMB Estimate:\\t {:.4f} \\t 95% Range: {:.4f} to {:.4f}\"\n",
    "      ).format(\n",
    "    obs_reward.mean(), \n",
    "    np.quantile(obs_reward, 0.025), \n",
    "    np.quantile(obs_reward, 0.975), \n",
    "    true_rl_reward.mean(), \n",
    "    np.quantile(true_rl_reward,0.025), \n",
    "    np.quantile(true_rl_reward, 0.975),\n",
    "    offpol_opt_reward_WIS_hard_train.mean(), \n",
    "    np.quantile(offpol_opt_reward_WIS_hard_train,0.025), \n",
    "    np.quantile(offpol_opt_reward_WIS_hard_train,0.975),\n",
    "    offpol_opt_reward_WIS_hard_ho.mean(), \n",
    "    np.quantile(offpol_opt_reward_WIS_hard_ho,0.025), \n",
    "    np.quantile(offpol_opt_reward_WIS_hard_ho,0.975),\n",
    "    offpol_opt_reward_mb.mean(), \n",
    "    np.quantile(offpol_opt_reward_mb,0.025), \n",
    "    np.quantile(offpol_opt_reward_mb,0.975)\n",
    "))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  },
  "toc-showcode": false,
  "toc-showmarkdowntxt": true,
  "toc-showtags": false
 },
 "nbformat": 4,
 "nbformat_minor": 2
}