<jupyter_start><jupyter_text>Vaccine Bandits We consider the problem of identifying optimal strategies for allocating influenza vaccine.Following~\citep{Libin:2017}, we format this problem as an instance of the BAI where each vaccine allocation strategy is an arm,and the reward of a strategy is defined as the proportion of individuals that did not experience symptomatic infection.  We generate 1000 rewards for each strategy by simulating the epidemic for 180 days using (FluTE)[https://github.com/dlchao/FluTE]  (with basic reproduction number $R0 = 1.3$).  The $0.5-$quantile (i.e. the median) is more robust to those outliers compared with the mean.  We consider the task of identifying $m=1$ (optimal arm: No. 9) and $m=3$ (optimal arms: No. 9, 25 and 13) arms with highest $0.5-$quantile values,   out of totally $K = 32$ arms with fixed budget.<jupyter_code># direct to proper path
import os
import sys
module_path = os.path.abspath(os.path.join('..'))
if module_path not in sys.path:
    sys.path.append(module_path)
    
import numpy as np
from ipywidgets import IntProgress
from IPython.display import display
import matplotlib.pyplot as plt
from collections import defaultdict
import seaborn as sns
import pickle
from prettytable import PrettyTable

from codes.BAI_FixedBudget import *
from codes.Environment import Mixture_AbsGau, setup_env, Exp
from codes.SimulatedGames import simulate
from codes.plots import *

%matplotlib inline

# import matplotlib.pylab as pylab
# params = {'legend.fontsize': 'x-large',
#           'figure.figsize': (15, 5),
#          'axes.labelsize': 'x-large',
#          'axes.titlesize':'x-large',
#          'xtick.labelsize':'x-large',
#          'ytick.labelsize':'x-large'}
# pylab.rcParams.update(params)

# font = {'family' : 'normal',
#         'weight' : 'bold',
#         'size'   : 15
#         }
# import matplotlib
# matplotlib.rc('font', **font)
# matplotlib.rcParams['pdf.fonttype'] = 42
# matplotlib.rcParams['ps.fonttype'] = 42
rewards_array = np.load('rewards.npy')
data_type = 'R0_1_dot_3'

# rewards_array = np.load('rewards-R0-1.5.npy')
# data_type = 'RO_1.5'

# rewards_array = np.load('rewards-R0-5.npy')
# data_type = 'RO_5'

# setting
results = defaultdict(dict)
num_exper = 5000
epsilon = 0 # no error
tau = 0.5 # quantile level

# Single best arm identification 
m= 3
budget_list = [500, 750, 1000, 1250, 1500, 1750, 2000, 2250, 2500]

# multiple best arms identification
# m = 3
# budget_list = [1000, 1250, 1500, 1750, 2000, 2250, 2500, 2750, 3000, 3250, 3500, 3750, 4000]

# output/plot setting

output_file = '../saved_results/Vaccine_' + data_type + '_' + 'm' + str(m) + str(tau)[-1] +'quant_' + str(num_exper) + 'exper_' + str(budget_list[0]) + '_' + str(budget_list[-1]) + '_budget'
# '../Group_Two_Vaccine_R0_5_m3_05quant_budget_500_2500_exper_2000.pickle'
budget_plot_title = 'Vaccine with '+ str(tau) + '-quantile (m=' + str(m) + ')'
# draw violin plot

plot_array = []
# rewards_array = 1 - rewards_array
for i in rewards_array:
    plot_array.append(list(i))

fig = plt.figure(1, figsize=(15, 4))
ax = fig.add_subplot(111)
b= ax.boxplot(plot_array, meanline=True, sym = 'k+', notch = 1, showfliers=False)
# code from https://matplotlib.org/gallery/statistics/customized_violin.html#sphx-glr-gallery-statistics-customized-violin-py

def adjacent_values(vals, q1, q3):
    upper_adjacent_value = q3 + (q3 - q1) * 1.5
    upper_adjacent_value = np.clip(upper_adjacent_value, q3, vals[-1])

    lower_adjacent_value = q1 - (q3 - q1) * 1.5
    lower_adjacent_value = np.clip(lower_adjacent_value, vals[0], q1)
    return lower_adjacent_value, upper_adjacent_value


def set_axis_style(ax, labels):
    ax.get_xaxis().set_tick_params(direction='out')
    ax.xaxis.set_ticks_position('bottom')
    ax.set_xticks(np.arange(0, len(labels)))
    ax.set_xticklabels(labels)
    ax.set_xlim(0.25, len(labels) + 0.75)
    ax.set_xlabel('Sample name')

quartile1, medians, quartile3 = np.percentile(plot_array, [25, 50, 75], axis=1)
whiskers = np.array([
    adjacent_values(sorted_array, q1, q3)
    for sorted_array, q1, q3 in zip(plot_array, quartile1, quartile3)])
whiskers_min, whiskers_max = whiskers[:, 0], whiskers[:, 1]

inds = np.asarray(list(np.arange(0, len(medians))))

fig,  ax2 = plt.subplots(nrows=1, ncols=1, figsize=(15, 4), sharey=True)

ax2.set_title('Vaccine Reward Violin Plot')
parts = ax2.violinplot(
        plot_array, positions=inds, showmeans=False, showmedians=False,
        showextrema=False)

for pc in parts['bodies']:
    pc.set_facecolor('skyblue')
    pc.set_edgecolor('black')
    pc.set_alpha(0.9)

sorted_idx = np.argsort(medians)[::-1]

ax2.scatter(inds, medians, marker='o', color='white', s=30, zorder=3)
ax2.scatter(inds[sorted_idx[0]], medians[sorted_idx[0]], marker='o', color='red', s=30, zorder=3)
ax2.scatter(inds[np.asarray(sorted_idx[1:3])], medians[np.asarray(sorted_idx[1:3])], marker='o', color='orange', s=30, zorder=3)
ax2.vlines(inds, quartile1, quartile3, color='grey', linestyle='-', lw=3)
ax2.vlines(inds, whiskers_min, whiskers_max, color='grey', linestyle='-', lw=1)
# ax2.hlines(np.max(medians),0,len(medians), colors='grey', lw=1)
# ax2.hlines(np.sort(medians)[::-1][2],0,len(medians), colors='grey', lw=1)

# set style for the axes
for i in plot_array:
    set_axis_style(ax, list(range(0,32)))

plt.subplots_adjust(bottom=0.15, wspace=0.05)
plt.xlabel('Arm (Vaccine Allocation Strategy)')
plt.ylabel('Reward')
# plt.show()
plt.savefig('Vaccine_Reward_Violin_Plot.pdf', bbox_inches='tight')
# for i in plot_array:
#     plt.hist(i)
class Vaccine_Env():
    """Env for Exponential Distribution.
    """
    def __init__(self, reward):
        self.reward = reward

    def sample(self, size = None):
        return np.random.choice(self.reward)

    def est_hazard_rate(self,bandwidth = 0.002):
        sorted_reward = np.sort(self.reward)
        x_list = np.arange(sorted_reward[0], sorted_reward[-1], 0.004)
        est_h_list = []
        for i, x in enumerate(x_list):
            count1 = 0
            count2 = 0
            for reward in sorted_reward:
                if x <= reward <= x + bandwidth:
                    count1 +=1
                if x<= reward:
                    count2+=1
            est_h = count1/count2/bandwidth
            est_h_list.append(est_h)

        plt.scatter(x_list, est_h_list)
        plt.title('Hazard Rate Estimation.')
        plt.xlabel('Reward')
        plt.ylabel('Estimate Hazard Rate')
Vaccine_Env(rewards_array[0]).est_hazard_rate()
def setup_env(rewards_array, ss_list = ['quantile_0.5'], random_set = None):
    """Setup environment for Vaccine simulations.
    """
    
    rewards_env = defaultdict(list)
    true_ss_dict = {}

    name = 'vaccine'
    true_ss_dict[name] = defaultdict(list)

    for arm_no, reward in enumerate(rewards_array):
        
        current_env = Vaccine_Env(reward)
                
        rewards_env[name].append(current_env)
                
        for ss in ss_list:
            ss_name = ss.split('_')[0]
            if len(ss.split('_'))> 1:
                ss_para = float(ss.split('_')[-1])
            if ss_name == 'quantile':
                true_ss_dict[name][ss].append(np.quantile(reward, ss_para))
            elif ss_name == 'mean':
                true_ss_dict[name][ss_name].append(np.mean(reward))
            else:
                assert True, 'Unknown summary statistics!'

    return rewards_env, true_ss_dict
rewards_env, true_ss_dict = setup_env(rewards_array, ['quantile_' + str(tau)])
rewards_env_mean, true_ss_dict_mean = setup_env(rewards_array, ['mean'])
np.sort(true_ss_dict['vaccine']['quantile_' + str(tau)])[::-1]
np.sort(true_ss_dict_mean['vaccine']['mean'])[::-1]
def run_exper(rewards_env, summary_stats, policy, epsilon, m, budget, num_exper, 
              fixed_samples_list = None):
    for key in rewards_env.keys():
        name = key + '_' + str(num_exper)
        policy_name = str(policy).split('.')[-1].split('\'')[0] + '-'
        subname = policy_name  + str(budget) 
        p = IntProgress(max = num_exper)
        p.description = 'Running'
        display(p)
        results[name][subname] = []

        result = simulate(rewards_env[key],summary_stats[key], policy, 
                epsilon, m, budget, num_exper, p, fixed_samples_list)
        results[name][subname].append(result)
if not os.path.isfile(output_file):
    for budget in budget_list:
        run_exper(rewards_env, true_ss_dict, Q_SR, epsilon, m, budget, num_exper)
        run_exper(rewards_env, true_ss_dict, Q_SAR_Simplified, epsilon, m, budget, num_exper)
        run_exper(rewards_env, true_ss_dict, batch_elimination, epsilon, m, budget, num_exper)
        run_exper(rewards_env, true_ss_dict, uniform_sampling, epsilon, m, budget, num_exper)
# save results

if os.path.isfile(output_file):
    with open(output_file, 'rb') as handle:
        results = pickle.load(handle)
else:
    with open(output_file, 'wb') as handle:
         pickle.dump(results, handle, protocol=pickle.HIGHEST_PROTOCOL)

plot_eva_budget(results, 'pe', paper_flag = True, with_para = True, title = budget_plot_title, log_scale = False)<jupyter_output>Budget  Probability of Error       std     Policy
0    500                0.5936  0.049116       Q-SR
1    500                0.6372  0.048081      Q-SAR
2    500                0.6544  0.047556       Q-BS
3    500                0.7728  0.041902  Q-Uniform
4    750                0.4972  0.049999       Q-SR
{'fraction_non_alphanumeric': 0.134134767478646, 'fraction_numerical': 0.03943899609828114, 'mean_word_length': 3.4990512333965844, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '636484', 'license': 'MIT', 'n_tokens_mistral': 3649, 'n_tokens_neox': 3386, 'n_words': 750, 'path': 'notebooks/Vaccine.ipynb', 'repo_name': 'Mengyanz/Quantile_Bandit', 'size': '9483'}