Resampling_Paper.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0b64be41",
   "metadata": {},
   "source": [
    "# Code for reproducing plots in \"Did we personalize? Assessing personalization by an online reinforcement learning algorithm using resampling\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3df36b20",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Setup\n",
    "import subprocess\n",
    "import numpy as np\n",
    "import sys\n",
    "import os\n",
    "from datetime import date, datetime\n",
    "from itertools import product\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.pyplot import figure\n",
    "import math\n",
    "import pickle as pkl\n",
    "import seaborn as sns\n",
    "\n",
    "import matplotlib as mpl\n",
    "import pylab\n",
    "from matplotlib import rc\n",
    "\n",
    "lsize=80\n",
    "axSize=67\n",
    "\n",
    "def setPlotSettings(fix_plot_settings=True):\n",
    "    if fix_plot_settings:\n",
    "        plt.rc('font', family='serif')\n",
    "        plt.rc('text', usetex=True)\n",
    "        plt.rcParams['text.usetex'] = True\n",
    "        label_size = lsize\n",
    "        mpl.rcParams['axes.labelsize'] = label_size\n",
    "        mpl.rcParams['axes.titlesize'] = label_size\n",
    "        mpl.rcParams['xtick.labelsize'] = axSize \n",
    "        mpl.rcParams['ytick.labelsize'] = axSize \n",
    "        mpl.rcParams['figure.titlesize'] = label_size\n",
    "        mpl.rcParams['lines.markersize'] = label_size\n",
    "        mpl.rcParams['grid.linewidth'] = 2.5\n",
    "        mpl.rcParams['legend.fontsize'] = label_size\n",
    "    \n",
    "        pylab.rcParams['xtick.major.pad']=5\n",
    "        pylab.rcParams['ytick.major.pad']=5\n",
    "\n",
    "        lss = ['--',  ':', '-.', '-', '--', '-.', ':', '-', '--', '-.', ':', '-']\n",
    "        mss = ['>', 'o',  's', 'D', '>', 's', 'o', 'D', '>', 's', 'o', 'D']\n",
    "        ms_size = [25, 20, 20, 20, 20, 20, 20, 20, 20, 20]\n",
    "        colors = ['#e41a1c', '#0000cd', '#4daf4a',  'black' , 'magenta']\n",
    "    \n",
    "np.set_printoptions(edgeitems=30, linewidth=100000, formatter=dict(float=lambda x: \"%.7g\" % x))\n",
    "\n",
    "NUSERS = 91\n",
    "F_KEYS=[\"intercept\", \"dosage\", \"engagement\", \"other_location\", \"variation\"]\n",
    "\n",
    "## Int Score Computation, at day level sliding window\n",
    "def computeMetricSlidingDay(result, indexFS, x=2, delta1=.5, delta2=.5, IGNORE_ETA=False):\n",
    "    ndata=rindex(result['availability'][:T], 1)+1 # +1 for actual number of timepoints. \n",
    "    if ndata==0: #if no one is found\n",
    "        print(\"no availability :(\")\n",
    "        return {'isEqualAvail':False, 'isEqualEngAvail':False, 'r1':None, 'r2':None, 'stdEffectRatio': None}\n",
    "\n",
    "    last_day=math.floor((ndata-1)/NTIMES) #ndata is the ts index of the last available time\n",
    "    \n",
    "    # Generate list of standardized adv forecasts for this user.\n",
    "    varValues=[] #binary 0/1 for \\var (z)\n",
    "    stdEffects=[] #the standardized effects for the user (\\hat{\\Delta}(\\cdot))\n",
    "    etas=[] # eta values\n",
    "\n",
    "    effLastDay=min(last_day+1+1,90)\n",
    "    for day in range(effLastDay): #want it to iterate to last day so +1. One more in case we stop before day 90 and can still forecast next day.\n",
    "        for time in range(NTIMES):\n",
    "            ts = (day) * 5 + time\n",
    "            if result['availability'][ts]==1.:\n",
    "                ## get mean/std and eta\n",
    "                # the one at ts is really the posteriors from the last day, or (ts//5)-1 since we update at end of the day. \n",
    "                # Consequently, the values at day 0, or times 0-4, correspond to the prior params\n",
    "                beta=result['post_beta_mu'][day*5][-len(F_KEYS):] \n",
    "                mean =result['fs'][ts] @ beta\n",
    "\n",
    "                sigma=result['post_beta_sigma'][day*5][-len(F_KEYS):, -len(F_KEYS):]\n",
    "                std=math.sqrt((result['fs'][ts] @ sigma.T) @ result['fs'][ts])\n",
    "                \n",
    "                eta=0 if IGNORE_ETA else result['etas'][ts]\n",
    "                \n",
    "                ## compute stdEffect\n",
    "                etas.append(eta)\n",
    "                stdEffect=(mean-eta)/std\n",
    "\n",
    "                stdEffects.append(stdEffect)\n",
    "                varValues.append(result['fs'][ts, indexFS])\n",
    "            else: #if not available put none\n",
    "                stdEffects.append(\"NONE\")\n",
    "                varValues.append(\"NONE\")\n",
    "\n",
    "    varValues=np.array(varValues)\n",
    "\n",
    "    # for computing int scores, track a few variables\n",
    "    ## int score 2, \\var\n",
    "    nSlidingWindows=NDAYS\n",
    "    nInterestingDeterminedWindows=0\n",
    "    nDeterminedWindows=0\n",
    "    \n",
    "    ## int score 1\n",
    "    nSlidingWindows_intscore1=NDAYS\n",
    "    nInterestingDeterminedWindows_intscore1=0\n",
    "    nDeterminedWindows_intscore1=0\n",
    "\n",
    "    avgVarEffAll=[]\n",
    "    avgNonVarEffAll=[]\n",
    "    determinedTimes=[]\n",
    "    \n",
    "    # loop through each day, and (1) form sliding windows, (2) check for G_{d,1}, (3) compute int_score\n",
    "    for day in range(effLastDay):\n",
    "        avail_idx_pre2=np.array([]) # var used for determining if an update occurred 2 days before?\n",
    "        avail_idx_pre1=np.array([]) # --\"\"-- 1 day before\n",
    "        avail_idx_cur=np.array([]) # --\"\"-- current day\n",
    "        if day == 0: #length of sliding window is 2*NTIMES\n",
    "            startTime=0\n",
    "            endTime=NTIMES*2\n",
    "\n",
    "            # check day 0 has any updates\n",
    "            avail_idx_cur = np.logical_and(~np.isnan(result['reward'][0:NTIMES]), result['availability'][0:NTIMES] == 1)\n",
    "\n",
    "        elif day >= 1 and day < last_day: #length of sliding window is 3*NTIMES\n",
    "            startTime=day*NTIMES-NTIMES\n",
    "            endTime=day*NTIMES+NTIMES*2\n",
    "\n",
    "            if day>=2:\n",
    "                avail_idx_pre2 = np.logical_and(~np.isnan(result['reward'][(day*NTIMES-2*NTIMES):(day*NTIMES-NTIMES)]), result['availability'][(day*NTIMES-2*NTIMES):(day*NTIMES-NTIMES)] == 1)\n",
    "            avail_idx_pre1 = np.logical_and(~np.isnan(result['reward'][(day*NTIMES-NTIMES):(day*NTIMES)]), result['availability'][(day*NTIMES-NTIMES):(day*NTIMES)] == 1)\n",
    "            avail_idx_cur = np.logical_and(~np.isnan(result['reward'][(day*NTIMES):(day*NTIMES+NTIMES)]), result['availability'][(day*NTIMES):(day*NTIMES+NTIMES)] == 1)\n",
    "\n",
    "        else: #if last_day, length of sliding window is 2*NTIMES\n",
    "            startTime=day*NTIMES-NTIMES\n",
    "            endTime=day*NTIMES+NTIMES\n",
    "            \n",
    "            # check day lastday-1 has any updates\n",
    "            avail_idx_pre2 = np.logical_and(~np.isnan(result['reward'][(day*NTIMES-2*NTIMES):(day*NTIMES-NTIMES)]), result['availability'][(day*NTIMES-2*NTIMES):(day*NTIMES-NTIMES)] == 1)\n",
    "            avail_idx_pre1 = np.logical_and(~np.isnan(result['reward'][(day*NTIMES-NTIMES):(day*NTIMES)]), result['availability'][(day*NTIMES-NTIMES):(day*NTIMES)] == 1)\n",
    "            if day < 89:\n",
    "                endTime=day*NTIMES+2*NTIMES\n",
    "                avail_idx_cur = np.logical_and(~np.isnan(result['reward'][(day*NTIMES):(day*NTIMES+NTIMES)]), result['availability'][(day*NTIMES):(day*NTIMES+NTIMES)] == 1)\n",
    "\n",
    "        ## Subset above varValues and forecasts to sliding window durations\n",
    "        varWindow=varValues[startTime:endTime]\n",
    "        forecastsWindow=stdEffects[startTime:endTime]\n",
    "\n",
    "        # check that an update happened before any of the day windows in question. \n",
    "        enoughUpdates = (sum(avail_idx_pre2)>0) or (sum(avail_idx_pre1)>0) or (sum(avail_idx_cur)>0) #a function of non observed reward. \n",
    "        \n",
    "        # G_{d,1} condition\n",
    "        varIndices=np.where(varWindow==1)[0]\n",
    "        nonVarIndices=np.where(varWindow==0)[0]\n",
    "        nBlue=len(varIndices)\n",
    "        nRed=len(nonVarIndices)\n",
    "        isDetermined = (nBlue >=x) and (nRed >= x) and enoughUpdates\n",
    "\n",
    "        #if G_{d,1} for intscore 2\n",
    "        if isDetermined: \n",
    "            nDeterminedWindows=nDeterminedWindows+1\n",
    "            determinedTimes.append(day)\n",
    "\n",
    "            # calculate day-level avg forecasts\n",
    "            avgVarEffect=np.mean(forecastsWindow[varIndices])\n",
    "            avgNonVarEffect=np.mean(forecastsWindow[nonVarIndices])\n",
    "\n",
    "            avgVarEffAll.append(avgVarEffect)\n",
    "            avgNonVarEffAll.append(avgNonVarEffect)\n",
    "            \n",
    "            # compare to determine if an 'interesting' window\n",
    "            if avgVarEffect > avgNonVarEffect:\n",
    "                nInterestingDeterminedWindows=nInterestingDeterminedWindows+1\n",
    "\n",
    "        # if G_{d,1} for intscore 1\n",
    "        if sum(avail_idx_pre1) > 0 or day==0: # if an update occurred the day prior\n",
    "            effects_intscore1=stdEffects[(day*NTIMES):(day*NTIMES+NTIMES)][effects_intscore1!=\"None\"]\n",
    "            #effects_intscore1=effects_intscore1[effects_intscore1!=\"None\"]\n",
    "            if len(effects_intscore1)>=2: #if we have enough data in the day window\n",
    "                nDeterminedWindows_intscore1=nDeterminedWindows_intscore1+1\n",
    "                if np.mean(effects_intscore1)>0: # if the window is interesting enough\n",
    "                    nInterestingDeterminedWindows_intscore1=nInterestingDeterminedWindows_intscore1+1\n",
    "\n",
    "    nUndeterminedSlidingWindows=nSlidingWindows-nDeterminedWindows\n",
    "    nUndeterminedSlidingWindows_intscore1=nSlidingWindows_intscore1-nDeterminedWindows_intscore1\n",
    "\n",
    "    # output int scores (one and two sided) and G_{d,1} fractions\n",
    "    statistic={}\n",
    "    # int score 2\n",
    "    if nSlidingWindows >0 and nDeterminedWindows >0:\n",
    "        statistic[\"r1\"]=nUndeterminedSlidingWindows/nSlidingWindows\n",
    "        statistic[\"rawR2\"]=nInterestingDeterminedWindows/nDeterminedWindows\n",
    "        statistic[\"r2\"]=abs(nInterestingDeterminedWindows/nDeterminedWindows - 0.5)\n",
    "        statistic[\"isInteresting_2\"]=(statistic[\"r1\"]<=delta1) and (statistic[\"r2\"]>=delta2)\n",
    "    else: \n",
    "        statistic[\"r1\"]=None\n",
    "        statistic[\"r2\"]=None\n",
    "        statistic[\"rawR2\"]=None\n",
    "        statistic[\"isInteresting_2\"]=None\n",
    "        \n",
    "    # int score 1\n",
    "    if nSlidingWindows_intscore1>0 and nDeterminedWindows_intscore1 > 0:\n",
    "        statistic[\"r3\"]=nUndeterminedSlidingWindows_intscore1/nSlidingWindows_intscore1#modified to just be for if there are enough updates.\n",
    "        statistic[\"rawR4\"]=nInterestingDeterminedWindows_intscore1/nDeterminedWindows_intscore1 #\n",
    "        statistic[\"r4\"]=abs(nInterestingDeterminedWindows_intscore1/nDeterminedWindows_intscore1 - 0.5) #\n",
    "        statistic[\"isInteresting_1\"]=(statistic[\"r3\"]<=delta1) and (statistic[\"r4\"]>=delta2)\n",
    "    else:\n",
    "        statistic[\"r3\"]=None\n",
    "        statistic[\"r4\"]=None\n",
    "        statistic[\"rawR4\"]=None\n",
    "        statistic[\"isInteresting_1\"]=None\n",
    "    \n",
    "    # to reproduce twin curves of avg of engaged and not engaged effects at determiend times\n",
    "    statistic[\"determinedTimes\"]=determinedTimes\n",
    "    statistic[\"avgNonValAll\"]=avgNonEngageAll\n",
    "    statistic[\"avgValAll\"]=avgEngageAll\n",
    "    \n",
    "    # to reproduce plot of standardized posterior at engaged and not engaged states\n",
    "    statistic['varValues']=varValues\n",
    "    statistic['standardizedEffects']=stdEffects\n",
    "    statistic['etas']=etas\n",
    "    return statistic\n",
    "\n",
    "\n",
    "def plot_C1_histogram(r1s, baseline):\n",
    "    image=\"./plots_test\"+'/histogram_C1_'+baseline+'.pdf'\n",
    "    setPlotSettings(True)\n",
    "    fig, ax = plt.subplots(figsize=(15, 15))\n",
    "    barcol='gray'\n",
    "    bins=\"auto\"\n",
    "    df=pd.DataFrame(r1s, columns=['gd1'])\n",
    "    p = sns.histplot(data=df, x='gd1', bins=bins, stat='count', ax=ax, color=barcol, cbar_kws={\"linewidth\":0}, line_kws={\"linewidth\":0}, linewidth=0)\n",
    "\n",
    "    # cosmetics\n",
    "    for spine in ['top', 'right']:\n",
    "        ax.spines[spine].set_visible(False)\n",
    "    plt.xlabel(\"\")\n",
    "    plt.ylabel(\"\\# Users\")\n",
    "    plt.grid(axis='y', alpha=.5, zorder=0) \n",
    "    ax=plt.gca()\n",
    "    \n",
    "    from matplotlib.ticker import MaxNLocator\n",
    "    ax.yaxis.set_major_locator(MaxNLocator(integer=True))\n",
    "    yticks = ax.yaxis.get_major_ticks()\n",
    "    yticks[0].label1.set_visible(False)\n",
    "    \n",
    "    plt.xlim([0,1])\n",
    "    ax.set_xticklabels([0,.25,.5,.75,1])\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(image, format=\"pdf\", bbox_inches=\"tight\")\n",
    "    \n",
    "    print(image)\n",
    "    plt.clf()\n",
    "    return\n",
    "\n",
    "\n",
    "def plot_C2_histogram(r2s, baseline):\n",
    "    image=\"./plots_test\"+'/histogram_C2_'+baseline+'.pdf'\n",
    "    setPlotSettings(True)\n",
    "    plt.rcParams['text.usetex'] = True\n",
    "    fig, ax = plt.subplots(figsize=(15, 15))\n",
    "    barcol='gray'\n",
    "    df=pd.DataFrame(rawR2s, columns=['rawR2s'])\n",
    "    bins=10\n",
    "    p = sns.histplot(data=df, x='rawR2s', bins=bins, stat='count', ax=ax, color=barcol, cbar_kws={\"linewidth\":0}, line_kws={\"linewidth\":0}, linewidth=0)\n",
    "    \n",
    "    # cosmetics\n",
    "    for spine in ['top', 'right']:\n",
    "        ax.spines[spine].set_visible(False)\n",
    "    plt.xlabel(\"\")\n",
    "    plt.ylabel(\"\\# Users\")\n",
    "    plt.grid(axis='y', alpha=.5, zorder=0) \n",
    "    ax=plt.gca()\n",
    "    \n",
    "    from matplotlib.ticker import MaxNLocator\n",
    "    ax.yaxis.set_major_locator(MaxNLocator(integer=True))\n",
    "    yticks = ax.yaxis.get_major_ticks()\n",
    "    yticks[0].label1.set_visible(False)\n",
    "    \n",
    "    plt.xlim([0,1])\n",
    "    ax.set_xticklabels([0,.25,.5,.75,1])\n",
    "    plt.axvline(.1, color='k', ls='--', zorder=4, lw=6)\n",
    "    plt.axvline(.9, color='k', ls='--', zorder=4, lw=6)\n",
    "    \n",
    "    plt.tight_layout()\n",
    "    plt.savefig(image, format=\"pdf\", bbox_inches=\"tight\")\n",
    "    \n",
    "    print(image)\n",
    "    plt.clf()\n",
    "    return\n",
    "\n",
    "def plot_heatmap_histogram(ogR1, ogR2, baseline, B=500):\n",
    "    # range of \\gamma and \\delta to go over\n",
    "    gammas=[.75, .70,.65]\n",
    "    deltas=[.35,.40,.45]\n",
    "    \n",
    "    heatMap=[]\n",
    "    heatMapR1=[] # g_{d,1}\n",
    "    output_dir = os.path.join(\"./output\", \"Baseline-\"+ str(baseline)+\"_UserSpecific-\"+str(\"False\"))\n",
    "\n",
    "    for gamma in gammas:\n",
    "        heatMapRow=[]\n",
    "        heatMapRowR1=[]\n",
    "        for delta in deltas:\n",
    "            bootstrapNInteresting=[] \n",
    "            observed=sum(np.logical_and(ogR1 <= gamma, ogR2 >= delta))\n",
    "            \n",
    "            for bootstrap in range(B): #\n",
    "                output_dirB = os.path.join(output_dir, \"Bootstrap-\" + str(bootstrap))\n",
    "                results=os.path.join(output_dirB, \"results.csv\")\n",
    "                df=pd.read_csv(results)\n",
    "\n",
    "                bootstrapR1s=np.array(df[r1Key])\n",
    "                if 'None' in bootstrapR1s:\n",
    "                    bootstrapR1s=bootstrapR1s[bootstrapR1s.astype(str) != 'None']\n",
    "                bootstrapR1s=bootstrapR1s.astype(float)\n",
    "\n",
    "                bootstrapR2s=np.array(df[r2Key])\n",
    "                if 'None' in bootstrapR2s:\n",
    "                    bootstrapR2s=bootstrapR2s[bootstrapR2s.astype(str) != 'None']\n",
    "                bootstrapR2s=bootstrapR2s.astype(float)\n",
    "\n",
    "                bootstrapValue=sum(np.logical_and(bootstrapR1s <= delta1, bootstrapR2s >= delta2))\n",
    "                bootstrapNInteresting.append(bootstrapValue)\n",
    "\n",
    "            perc=stats.percentileofscore(bootstrapNInteresting, observed, 'weak')/100\n",
    "            heatMapRow.append(1-perc)\n",
    "            heatMapRowR1.append(sum(ogR1<=gamma))\n",
    "\n",
    "            #plot histogram too!\n",
    "            if delta1==.75 and delta2==.4: # can generate for other values too\n",
    "                image=\"./plots_test\"+'/histogram_Interesting_'+baseline+\"_delta1=\"+str(gamma)+\"_delta2=\"+str(delta)+\"_B=\"+str(B)+'.pdf'\n",
    "                plt.clf()\n",
    "                setPlotSettings(True)\n",
    "                fig, ax = plt.subplots(figsize=(15, 15))\n",
    "                barcol='gray'\n",
    "                bins=\"auto\"\n",
    "                df=pd.DataFrame(bootstrapNInteresting, columns=['nInt'])\n",
    "                p = sns.histplot(data=df, x='nInt', bins=bins, stat='probability', ax=ax, color=barcol, cbar_kws={\"linewidth\":0}, line_kws={\"linewidth\":0}, linewidth=0)\n",
    "                \n",
    "                # cosmetics\n",
    "                for spine in ['top', 'right']:\n",
    "                    ax.spines[spine].set_visible(False)\n",
    "                plt.xlabel(\"\")\n",
    "                plt.ylabel(\"Proportion\")\n",
    "                plt.grid(axis='y', alpha=.5, zorder=0)\n",
    "                plt.axvline(observed, color='b', ls='-.', zorder=4, lw=6)\n",
    "                plt.xlim(left=0)\n",
    "                if baseline==\"Zero\" and delta1==.75 and delta2==.4:\n",
    "                    labs=['0', '0.02','0.04','0.06','0.08','0.10','0.12','0.14']\n",
    "                    vals=[float(lab) for lab in labs]\n",
    "                    labs[0]=''\n",
    "                    plt.yticks(vals,labs)\n",
    "                elif baseline==\"other_location\" and delta1==.75 and delta2==.4:\n",
    "                    labs=['0', '0.05','0.10','0.15','0.20','0.25','0.30','0.35']\n",
    "                    vals=[float(lab) for lab in labs]\n",
    "                    labs[0]=''\n",
    "                    plt.yticks(vals,labs)\n",
    "                elif baseline==\"variation\" and delta1==.75 and delta2==.4:\n",
    "                    labs=['0', '0.025','0.050','0.075','0.100','0.125','0.150','0.175']\n",
    "                    vals=[float(lab) for lab in labs]\n",
    "                    labs[0]=''\n",
    "                    plt.yticks(vals,labs)\n",
    "                else:\n",
    "                    yticks = ax.yaxis.get_major_ticks() \n",
    "                    yticks[0].label1.set_visible(False)\n",
    "                plt.tight_layout()\n",
    "                plt.savefig(image, format=\"pdf\", bbox_inches=\"tight\")\n",
    "                print(image)\n",
    "        \n",
    "        heatMap.append(heatMapRow)\n",
    "        heatMapR1.append(heatMapRowR1)\n",
    "\n",
    "    # Now plot the heatmap!\n",
    "    sns.set(rc={'text.usetex': True})\n",
    "    plt.clf()\n",
    "    setPlotSettings(True)\n",
    "    fig, ax=plt.subplots(figsize=(15,15))\n",
    "    \n",
    "    heatMap=pd.DataFrame(heatMap)\n",
    "    heatMapR1=pd.DataFrame(heatMapR1)\n",
    "    \n",
    "    gammaLabels=['0.75', '0.70','0.65']\n",
    "    deltaLabels=['0.35','0.40','0.45']\n",
    "    s=sns.heatmap(heatMap,xticklabels=deltaLabels, ax=ax, yticklabels=gammaLabels, cmap=\"Blues_r\", fmt='.2f', annot=True, vmin=0, vmax=1, annot_kws={'fontsize': lsize})\n",
    "    plt.yticks(rotation=0)\n",
    "    s.set(xlabel=\"$\\delta$\", ylabel=\"$\\gamma$\")\n",
    "    f=s.get_figure()\n",
    "    \n",
    "    toWrite=baseline\n",
    "    if baseline==\"other_location\":\n",
    "        toWrite=\"other location\"\n",
    "    image=\"./plots_test\"+'/heatMap_Interesting_'+toWrite+\"_B=\"+str(B)+'.pdf'\n",
    "    print(image)\n",
    "    plt.tight_layout()\n",
    "    f.savefig(image, format=\"pdf\", bbox_inches=\"tight\")\n",
    "    plt.clf()\n",
    "\n",
    "# get original int scores and G_{d,1} ratios.\n",
    "def getOriginalResults(indexFS, delta=.4, gamma=.75):\n",
    "    baseline=\"Zero\"\n",
    "    if indexFS!=-1:\n",
    "        baseline=F_KEYS[indexFS]\n",
    "\n",
    "    original_result=\"./init/original_result_91.pkl\"\n",
    "    with open(original_result, 'rb') as handle:\n",
    "        original_result=pkl.load(handle)\n",
    "\n",
    "    r1s=[]\n",
    "    r2s=[]\n",
    "    rawR2s=[]\n",
    "    interesting=[]\n",
    "    \n",
    "    for result in original_result:\n",
    "        result=computeMetricSlidingDay(result, experiment, delta1=gamma, delta2=delta)\n",
    "\n",
    "        # get proper results\n",
    "        r1=result['r1']\n",
    "        r2=result['r2']\n",
    "        rawR2=result['rawR2']\n",
    "        interestingResult=result['isInteresting_2']\n",
    "        if baseline==\"Zero\":\n",
    "            r1=result['r3']\n",
    "            r2=result['r4']\n",
    "            rawR2=result['rawR4']\n",
    "            interestingResult=result['isInteresting_1']\n",
    "\n",
    "        # store relevant values\n",
    "        if r1 != None:\n",
    "            r1s.append(r1)\n",
    "            if r1 <= delta1: # if they are worth consideration (have enough data)\n",
    "                r2s.append(r2)\n",
    "                rawR2s.append(rawR2)\n",
    "        interesting.append(interestingResult)\n",
    "        \n",
    "    r1s=np.array(r1s)\n",
    "    r2s=np.array(r2s)\n",
    "    rawR2s=np.array(rawR2s)\n",
    "    \n",
    "    # now r2s by users\n",
    "    r2Users={}\n",
    "    r2RawUsers={}\n",
    "    for i in users:\n",
    "        r2Users[str(i)]=[]\n",
    "        r2RawUsers[str(i)]=[]\n",
    "        output_dir = os.path.join(\"./output\", \"Baseline-\"+ str(baseline)+\"_UserSpecific-\"+str(True)+\"_User-\" + str(i))\n",
    "        results=os.path.join(output_dir, \"results.csv\")\n",
    "        results=pd.read_csv(results)\n",
    "        r2RawUsers[str(i)]=results[r2RawKey]\n",
    "        r2Users[str(i)]=results[r2Key]\n",
    "\n",
    "    return r1s,r2s,rawR2s,interesting, r2Users, r2RawUsers, original_result\n",
    "\n",
    "def plot_r2_User(r2Users, r2RawUsers, rawR2, interesting, baseline, users):\n",
    "    percs=[]\n",
    "    for i in range(len(users)):\n",
    "        user=users[i]\n",
    "        if interesting[user]==True: #only added if its interesting\n",
    "            perc=stats.percentileofscore(r2Users[str(i)], float(holderR2s[i]), 'weak')/100\n",
    "            percs.append(1-perc)\n",
    "\n",
    "    #plot histogram of Frac\n",
    "    image=\"./plots_test\"+'/histogram_percs_'+baseline+'_x=2'+'.pdf'\n",
    "    plt.clf()\n",
    "    setPlotSettings(True)\n",
    "    fig, ax=plt.subplots(figsize=(15,15))\n",
    "    plt.grid(axis='y', alpha=.5, zorder=0)\n",
    "    barcol=\"gray\"\n",
    "    df=pd.DataFrame(percs, columns=['percs'])\n",
    "    p = sns.histplot(data=df, x='percs', stat='count', ax=ax, color=barcol, cbar_kws={\"linewidth\":0}, line_kws={\"linewidth\":0}, linewidth=0)\n",
    "    \n",
    "    # cosmetics\n",
    "    for spine in ['top', 'right']:\n",
    "        ax.spines[spine].set_visible(False)\n",
    "    ax=plt.gca()\n",
    "    from matplotlib.ticker import MaxNLocator\n",
    "    ax.yaxis.set_major_locator(MaxNLocator(integer=True))\n",
    "    yticks = ax.yaxis.get_major_ticks()\n",
    "    yticks[0].label1.set_visible(False)\n",
    "    \n",
    "    plt.xlabel(\"\")\n",
    "    plt.ylabel(\"\\# Users\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(image, format=\"pdf\",bbox_inches=\"tight\")\n",
    "    print(image)\n",
    "    plt.clf()\n",
    "    \n",
    "    # Now plot histogram of int scores!\n",
    "    for i in users:\n",
    "        if ((i == 4 and baseline==\"variation\") or (i==77 and baseline==\"Zero\")):\n",
    "            image=\"./plots_test\"+'/histogram_c2s_user_'+str(i)+'_'+baseline+'_x=2'+'.pdf'\n",
    "            plt.clf()\n",
    "            r2s=r2RawUsers[str(i)]\n",
    "            \n",
    "            setPlotSettings(True)\n",
    "            fig, ax=plt.subplots(figsize=(24,16))#, dpi=80)\n",
    "            plt.grid(axis='y', alpha=.5, zorder=0)\n",
    "            barcol=\"gray\"\n",
    "            df=pd.DataFrame(np.array(r2s), columns=['r2s'])\n",
    "            p = sns.histplot(data=df, x='r2s', stat='probability', ax=ax, color=barcol, cbar_kws={\"linewidth\":0}, line_kws={\"linewidth\":0}, linewidth=0)\n",
    "            \n",
    "            # cosmetics\n",
    "            for spine in ['top', 'right']:\n",
    "                ax.spines[spine].set_visible(False)\n",
    "            ax=plt.gca()\n",
    "            plt.xlabel(\"\")\n",
    "            plt.ylabel(\"Proportion\")\n",
    "            plt.axvline(allRawR2[i], ls='-.', color='b', zorder=4, lw=6)\n",
    "            plt.xlim(left=-.05)\n",
    "            plt.xticks([0, .25,.5,.75,1])\n",
    "            \n",
    "            yticks = ax.yaxis.get_major_ticks()\n",
    "            yticks[0].label1.set_visible(False)\n",
    "            plt.tight_layout()\n",
    "            plt.savefig(image, format=\"pdf\", bbox_inches=\"tight\")\n",
    "            print(image)\n",
    "            plt.clf()\n",
    "\n",
    "            \n",
    "def plotUserDayDateAndResims(result, resultRun, user, outputPath, stateName, bs_results, bs_resultsRuns, bs):\n",
    "    pathPeng=\"./init/all91_uid.pkl\"\n",
    "    resultPeng=pkl.load(open(pathPeng, \"rb\"))\n",
    "    uids=resultPeng[:,0,15]\n",
    "    otherDF=pd.DataFrame(uids, range(resultPeng.shape[0]), columns=[\"StudyID\"])\n",
    "    otherDF['indices']=range(len(uids))\n",
    "\n",
    "    pathBaselineInfo=\"./init/baseline_info.csv\"\n",
    "    baseline_info=pd.read_csv(pathBaselineInfo)\n",
    "    baseline_info=baseline_info[['start.date', 'StudyID']] # match by UIDS!\n",
    "\n",
    "    start_dates=pd.merge(otherDF, baseline_info, on=\"StudyID\")\n",
    "\n",
    "    # determined times, standardized effects\n",
    "    rc('mathtext', default='regular')\n",
    "    fig = plt.figure(figsize=(24, 16))\n",
    "    ax = fig.add_subplot()#111)#, aspect='equal')\n",
    "    spacing = .4\n",
    "\n",
    "    ndata=rindex(resultRun['availability'][:T], 1)+1 # +1 for actual number of timepoints. #np.where(data[:T,1]==0)[0] #first time the user is out\n",
    "    last_day=math.floor((ndata-1)/NTIMES) #ndata is the ts index of the last available time, reset the -1, then /NTIMES to get the day of the last available time.\n",
    "    start = dt.datetime.strptime(start_dates.iloc[user][2]+\" 00\", '%Y-%m-%d %H')#.date()\n",
    "    end = start + dt.timedelta(days=last_day)\n",
    "    plotEnd = start + dt.timedelta(days=last_day+1)\n",
    "\n",
    "    xs=np.array(range(len(result['varValues'])))\n",
    "    availTimes = np.logical_and(~np.isnan(resultRun['reward']), resultRun['availability'] == 1)[:len(xs)]\n",
    "    availTimes=np.where(resultRun['availability'] == 1)[0]\n",
    "    availsAndNonNanReward=xs[availTimes]\n",
    "    mapper=['red','blue']\n",
    "\n",
    "    markersize=650\n",
    "    opacity=.7\n",
    "    colors=[mapper[int(i)] for i in result['varValues'][availTimes]]\n",
    "    x = [ ]\n",
    "    hourInc=[0,5,10,15,20]\n",
    "    for i in range(NDAYS):\n",
    "        for j in range(NTIMES):\n",
    "            x.append((start + dt.timedelta(days=i)+dt.timedelta(hours=hourInc[j]))) \n",
    "    x=np.array(x)\n",
    "\n",
    "    xs=x[availsAndNonNanReward]\n",
    "    y=result['standardizedEffects'][availsAndNonNanReward]\n",
    "\n",
    "    ax.set_ylim([-2, 2.75])\n",
    "    vals=[-2,-1,0,1,2, 2.85]\n",
    "    labs=['-2','-1','0','1','2','']\n",
    "    ax.set_yticks(vals,labs)\n",
    "\n",
    "    ax.axhline(y=0, color='k', linestyle='--', linewidth=6, alpha=.75)\n",
    "    opacityB=.2\n",
    "    yb=bs_results[len(bs_results)-1]['standardizedEffects'][availsAndNonNanReward]\n",
    "    ax.scatter(xs[result['varValues'][availTimes] == 1], yb[result['varValues'][availTimes] == 1], marker='^', color='red',zorder=2, alpha=opacity,s=markersize)\n",
    "    ax.scatter(xs[result['varValues'][availTimes] == 0], yb[result['varValues'][availTimes] == 0], marker='o', color='blue',zorder=2, alpha=opacity,s=markersize)\n",
    "    \n",
    "    colors=[mapper[int(i)] for i in result['varValues'][availTimes]]\n",
    "\n",
    "    from matplotlib.lines import Line2D\n",
    "    markersize=35\n",
    "    lSize=40\n",
    "    labelBlue=stateName+\" = 0\"\n",
    "    labelRed=stateName+\" = 1\"\n",
    "    legend_elements = [\n",
    "                    Line2D([0], [0], marker='o', color='w', label=labelBlue,\n",
    "                            markerfacecolor='b', markersize=markersize, alpha=opacity),\n",
    "                    Line2D([0], [0], marker='^', color='w', label=labelBlue,\n",
    "                            markerfacecolor='r', markersize=markersize, alpha=opacity)\n",
    "                      ]\n",
    "   \n",
    "    ax.legend(handles=legend_elements, loc=\"upper left\",fancybox = True, shadow = True, fontsize=lSize, handletextpad=-.2)\n",
    "\n",
    "    for spine in ['top', 'right']:\n",
    "        ax.spines[spine].set_visible(False)\n",
    "    ax.grid(True, alpha=0.4)\n",
    "\n",
    "    lSize=80\n",
    "    ax.set_xlabel('Date (MM-DD)', fontsize=lSize)\n",
    "    ax.set_ylabel(\"Std. Advantage\", fontsize=lSize)\n",
    "\n",
    "    lSize=75\n",
    "    newLabels=[i.strftime('%m-%d') for i in xs]\n",
    "    newXTicks=[]\n",
    "    newLabels2=[]\n",
    "    for i in range(len(x)):\n",
    "        if i%(35*2)==0 and i!=0:\n",
    "            newXTicks.append(x[i])\n",
    "            newLabels2.append(x[i].strftime('%m-%d'))\n",
    "    ax.set_xticks(newXTicks, newLabels2)\n",
    "    ax.tick_params(axis='both', which='major', labelsize=lSize)\n",
    "    ax.tick_params(axis='both', which='minor', labelsize=lSize)\n",
    "\n",
    "    ax.set_xlim([start,plotEnd])\n",
    "    myFmt = mdates.DateFormatter('%m-%d')\n",
    "    ax.xaxis.set_major_formatter(myFmt)\n",
    "    xticks = ax.xaxis.get_major_ticks()\n",
    "    xticks[0].label1.set_visible(False)\n",
    "\n",
    "    print(outputPath+'/blueRedOverAll_user-'+str(user)+'-state-'+stateName + \"_resim-\"+str(bs)+'.pdf')\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(outputPath+'/blueRedOverAll_user-'+str(user)+'-state-'+stateName + \"_resim-\"+str(bs)+'.pdf')\n",
    "    plt.clf()\n",
    "    return \n",
    "\n",
    "def plotTrajectories(original_result, baseline, indexFS,gamma=.75,delta=.4):\n",
    "    subset=2\n",
    "    trajectories=random.sample(range(B), subset)\n",
    "    bs_results=[]\n",
    "    bs_results_sliding=[]\n",
    "    for bootstrap in trajectories:\n",
    "        print(\"processing \"+str(bootstrap))\n",
    "        output_dir = os.path.join(\"./output\", \"Baseline-\"+ str(baseline)+\"_UserSpecific-\"+str(True)+\"_User-\"+str(i), \"Bootstrap-\" + str(bootstrap))\n",
    "        onlyfiles = [f for f in listdir(output_dir) if isfile(join(output_dir, f))]\n",
    "        if len(onlyfiles)>0:\n",
    "            filepath=os.path.join(output_dir, onlyfiles[0])\n",
    "            bs_res=pkl.load(open(filepath, \"rb\"))\n",
    "            bs_results.append(bs_res)\n",
    "            bs_results_sliding.append(computeMetricSlidingDay(bs_res, indexFS, delta1=gamma, delta2=delta))\n",
    "            plotUserDayDateAndResims(resultSliding, original_result[i], i, \"./plots_test\", F_KEYS[indexFS], bs_results_sliding, bs_results, bootstrap)\n",
    "\n",
    "B=500\n",
    "for index in [2]:#[-1,2,3,4]: #for zero, engagement, location, and variation respectively\n",
    "    r1s,r2s, rawR2s,baseline,interesting, r2_BS_Users, r2Raw_BS_Users, original_result=getResults(index)\n",
    "    plot_C1_histogram(r1s, baseline)\n",
    "    plot_C2_histogram(rawR2s, baseline)\n",
    "    plot_heatmap_histogram(r1s, r2s, baseline, B=B)\n",
    "    \n",
    "    users=[4,77]\n",
    "    plot_C2_User(r2_BS_Users, r2Raw_BS_Users, interesting, baseline, users)\n",
    "    plot_Trajectories(original_result, baseline, index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e17bbc6b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}