Figure 4.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# imports\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy import signal\n",
    "import pandas as pd\n",
    "from glob import glob as glob\n",
    "import os\n",
    "from neo import io\n",
    "import seaborn as sb\n",
    "from sklearn import preprocessing\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['pdf.fonttype'] = 42\n",
    "mpl.rcParams['ps.fonttype'] = 42\n",
    "plt.rcParams[\"font.size\"]= 12\n",
    "plt.rcParams[\"font.family\"] = \"Arial\"\n",
    "from sklearn.metrics import roc_auc_score\n",
    "\n",
    "# functions \n",
    "from skimage.measure import label\n",
    "\n",
    "def clustertest(true,permutations,thresh):\n",
    "    nperm = permutations.shape[-1]\n",
    "    t_val = abs((true - np.mean(permutations,-1)) / np.std(permutations,-1)) # 2-sided\n",
    "    above_threshold = t_val>thresh\n",
    "    sig_clust = label(above_threshold)\n",
    "    clust_sum = np.zeros(np.max(sig_clust))\n",
    "    for c in range(1,max(sig_clust)+1):\n",
    "        \n",
    "        clust_sum[c-1] = np.sum(t_val[sig_clust==c])\n",
    "    # get random clusters\n",
    "    clust_sum_rand = np.zeros(nperm)\n",
    "    for perm in range(nperm):\n",
    "        compare_to = np.delete(permutations,perm,axis=-1)\n",
    "        t_val = abs((permutations[:,perm] - np.mean(compare_to,-1)) / np.std(compare_to,-1))\n",
    "        above_threshold = t_val>thresh\n",
    "        sig_clust_rand = label(above_threshold)\n",
    "        if np.max(sig_clust_rand)>0:\n",
    "            clust_sum_r = np.zeros(np.max(sig_clust_rand))\n",
    "            for c in range(1,max(sig_clust_rand)+1):\n",
    "                clust_sum_r[c-1] = np.sum(t_val[sig_clust_rand==c])\n",
    "            clust_sum_rand[perm] = np.max(clust_sum_r)\n",
    "    # compard random and true\n",
    "    P = np.zeros(true.shape)*np.nan\n",
    "    for c in range(len(clust_sum)):\n",
    "        P[sig_clust==c+1] = np.mean(clust_sum[c]<=clust_sum_rand)\n",
    "    return P\n",
    "\n",
    "def getAUROC(true,pred):\n",
    "    ''' compute AUC of classifier for a temporal generalization matrix \n",
    "    pred is a matrix of shape ntrials x nbins x nbins'''\n",
    "    if len(pred.shape)==3:\n",
    "        # AUC time x time matrix\n",
    "        nbins = pred.shape[1]\n",
    "        AUC = np.zeros((nbins,nbins))\n",
    "        for b in range(nbins):\n",
    "            for b2 in range(nbins):\n",
    "                AUC[b,b2] = roc_auc_score(true,pred[:,b,b2])\n",
    "    elif len(pred.shape)==2:\n",
    "        nbins = pred.shape[1]\n",
    "        AUC = np.zeros(nbins)\n",
    "        for b in range(nbins):\n",
    "            if len(true.shape)>1:\n",
    "                AUC[b] = roc_auc_score(true[:,b],pred[:,b])\n",
    "            else:\n",
    "                AUC[b] = roc_auc_score(true,pred[:,b])\n",
    "    elif len(pred.shape)==1:\n",
    "        AUC = roc_auc_score(true,pred)\n",
    "    return AUC\n",
    "\n",
    "\n",
    "def designMatrix(globdev,local,stimulus=None,include_stim=False):\n",
    "    '''\n",
    "    design matrix of nsamples x nvariables\n",
    "    Variables are: Stimulus, Local deviance, Global deviance, Blocktype (xx or xY)'''\n",
    "    \n",
    "    ntr = len(local)\n",
    "    loc = local.copy();\n",
    "    glob = globdev.copy()\n",
    "    if include_stim:\n",
    "        stim = stimulus.copy()\n",
    "        X = np.vstack((loc,glob,stim,np.ones(ntr))).astype(float).T\n",
    "    else:\n",
    "        X = np.vstack((loc,glob,np.ones(ntr))).astype(float).T\n",
    "        \n",
    "    return X\n",
    "\n",
    "def designMatrix_first3(blocks,stimulus):\n",
    "    '''\n",
    "    design matrix of nsamples x nvariables\n",
    "    Variables are: Stimulus, Local deviance, Global deviance, Blocktype (xx or xY)'''\n",
    "    \n",
    "    ntr = len(blocks)\n",
    "    context = blocks.copy()\n",
    "    stim = stimulus.copy()\n",
    "    X = np.vstack((context,stim,np.ones(ntr))).astype(float).T\n",
    "    return X\n",
    "\n",
    "def projectMUA(Coef,Y,Sess,Testtrials,compute_matrix = False):\n",
    "    '''\n",
    "    Projection of binned population MUA onto coefficient vector from \n",
    "    multivariate linear regression with cross-validation\n",
    "    Projection onto single coefficient vector\n",
    "    '''\n",
    "    ntr,nch,nbins = Y.shape\n",
    "    nsess = len(np.unique(Sess))\n",
    "    nvar = Coef.shape[1]   \n",
    "    \n",
    "    Proj = np.zeros((nsess,nvar,ntr,nbins))# pre-allocate\n",
    "    # project individual trials\n",
    "    # use hold-out data per fold\n",
    "    K = Coef.shape[-1]\n",
    "    for sess in range(nsess):\n",
    "        print('Session %s'%(sess+1))\n",
    "        sessind = np.where(Sess==sess)[0] # trial indices from session\n",
    "        Ysess = Y[sessind,:] # MUA from session\n",
    "        \n",
    "        proj = np.zeros((nvar,len(sessind),nbins))\n",
    "        for k in range(K): \n",
    "            test_ind = Testtrials[sess][k]\n",
    "            # projection onto coefficients from each bin\n",
    "            coef = preprocessing.normalize(Coef[sess,:,:,k],norm='l2',axis=1)\n",
    "            for b in range(nbins):\n",
    "                \n",
    "                # project all bins onto this axis \n",
    "                proj[:,test_ind,b] = coef@Ysess[test_ind,:,b].T\n",
    "                        \n",
    "                # get inter-session projection, only for one fold\n",
    "                if k==0:\n",
    "                    Proj[sess,:][:,Sess!=sess,b] = coef@Y[Sess!=sess,:,b].T\n",
    "            \n",
    "        Proj[sess,:][:,sessind,:] = proj\n",
    "        \n",
    "    return Proj\n",
    "\n",
    "def projectMUA_cv(Coef,Y,Sess,Testtrials,compute_matrix = False):\n",
    "    '''\n",
    "    Projection of binned population MUA onto coefficient vector from \n",
    "    multivariate linear regression with cross-validation\n",
    "    compute_matrix: \n",
    "    '''\n",
    "    nbins_train = Coef.shape[-2]\n",
    "    ntr,nch,nbins = Y.shape\n",
    "    nsess = len(np.unique(Sess))\n",
    "    nvar = Coef.shape[1]   \n",
    "    if compute_matrix == True:\n",
    "        Proj = np.zeros((nsess,nvar,ntr,nbins_train,nbins))# pre-allocate\n",
    "    else:\n",
    "        Proj = np.zeros((nsess,nvar,ntr,nbins_train))# pre-allocate\n",
    "    # project individual trials\n",
    "    # use hold-out data per fold\n",
    "    K = Coef.shape[-1]\n",
    "    for sess in range(nsess):\n",
    "        print('Session %s'%(sess+1))\n",
    "        sessind = np.where(Sess==sess)[0] # trial indices from session\n",
    "        Ysess = Y[sessind,:] # MUA from session\n",
    "        \n",
    "        if compute_matrix == True:\n",
    "            proj = np.zeros((nvar,len(sessind),nbins_train,nbins))\n",
    "        else:\n",
    "            proj = np.zeros((nvar,len(sessind),nbins_train))\n",
    "        for k in range(K): \n",
    "            test_ind = Testtrials[sess][k]\n",
    "            # projection onto coefficients from each bin\n",
    "            for b in range(nbins_train):\n",
    "                coef = preprocessing.normalize(Coef[sess,:,:,b,k],norm='l2',axis=1)\n",
    "\n",
    "                # project all bins onto these axes (temporal generalization matrix)\n",
    "                if compute_matrix == True:\n",
    "                    for btest in range(nbins):\n",
    "                        proj[:,test_ind,b,btest] = coef@Ysess[test_ind,:,btest].T\n",
    "                        \n",
    "                        # get inter-session projection, only for one fold\n",
    "                        if k==0:\n",
    "                            Proj[sess,:][:,Sess!=sess,b,btest] = coef@Y[Sess!=sess,:,btest].T\n",
    "                # project only same bin as training\n",
    "                else:\n",
    "                    proj[:,test_ind,b] = coef@Ysess[test_ind,:,b].T\n",
    "                        \n",
    "                    # get inter-session projection, only for one fold\n",
    "                    if k==0:\n",
    "                        Proj[sess,:][:,Sess!=sess,b] = coef@Y[Sess!=sess,:,b].T\n",
    "            \n",
    "        Proj[sess,:][:,sessind,:] = proj\n",
    "        \n",
    "    return Proj\n",
    "\n",
    "\n",
    "def balance_data(Xtrain, train_ind,include_stim):\n",
    "    '''balance classes in data'''\n",
    "    N = []\n",
    "    if include_stim:\n",
    "        for i,v1 in enumerate(np.unique(Xtrain)):\n",
    "            for j,v2 in enumerate(np.unique(Xtrain)):\n",
    "                for l,v3 in enumerate(np.unique(Xtrain)):\n",
    "                # count trials of condition\n",
    "                    n = sum((Xtrain[:,0]==v1) & (Xtrain[:,1]==v2) & (Xtrain[:,3]==v3))\n",
    "                    N.append(n)\n",
    "        nmin = np.min(N)\n",
    "        # subsample\n",
    "        keep_ind = []\n",
    "        for i,v1 in enumerate(np.unique(Xtrain)):\n",
    "            for j,v2 in enumerate(np.unique(Xtrain)):\n",
    "                for l,v3 in enumerate(np.unique(Xtrain)):\n",
    "                    # count trials of condition\n",
    "                    condind = np.where((Xtrain[:,0]==v1) & (Xtrain[:,1]==v2) & (Xtrain[:,3]==v3))[0] #balance local, global, stimulus\n",
    "                    keep_ind.extend(condind[np.random.permutation(len(condind))[:nmin]])\n",
    "    else:\n",
    "        for i,v1 in enumerate(np.unique(Xtrain)):\n",
    "            for j,v2 in enumerate(np.unique(Xtrain)):\n",
    "                # count trials of condition\n",
    "                n = sum((Xtrain[:,0]==v1) & (Xtrain[:,1]==v2))\n",
    "                N.append(n)\n",
    "        nmin = np.min(N)\n",
    "        # subsample\n",
    "        keep_ind = []\n",
    "        for i,v1 in enumerate(np.unique(Xtrain)):\n",
    "            for j,v2 in enumerate(np.unique(Xtrain)):\n",
    "                # count trials of condition\n",
    "                condind = np.where((Xtrain[:,0]==v1) & (Xtrain[:,1]==v2))[0] #balance local, global, stimulus\n",
    "                keep_ind.extend(condind[np.random.permutation(len(condind))[:nmin]])\n",
    "    train_ind = train_ind[keep_ind]\n",
    "    return train_ind\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load preprocessed data and extract sequence variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['PFC_MU', 'PPC_MU', 'PFC_SU', 'PPC_SU', 'TrialID', 'ItemID', 'StimID', 'StimName', 'StimOn', 'blockID', 'blockType', 'date', 'StimDur', 'ISIDur', 'RewardOn']\n",
      "Recording dates: ['20200226' '20200228' '20200305' '20200306' '20200311' '20200312']\n",
      "Day: 20200226\n",
      "Day: 20200228\n",
      "Day: 20200305\n",
      "Day: 20200306\n",
      "Day: 20200311\n",
      "Day: 20200312\n"
     ]
    }
   ],
   "source": [
    "# parameters\n",
    "animal_id = 'A11' # A11, H07\n",
    "exp = 'classical' # 'wang or classical\n",
    "out_path = '/Volumes/Bellet/Local_Global/processedData/%s/'%animal_id # path to the data\n",
    "\n",
    "# load preprocessed data\n",
    "df = pd.read_pickle(os.path.join(out_path,'%s_stims_spikes_dataframe.pkl'%(animal_id)))\n",
    "\n",
    "# extract basic variables\n",
    "dates = np.unique(df.date)\n",
    "ndates = len(dates)\n",
    "nitems = np.max(df.ItemID)\n",
    "stimID = np.unique(df.StimID)\n",
    "nstim = len(stimID)\n",
    "print(list(df.keys()))\n",
    "print('Recording dates:',dates)\n",
    "\n",
    "soa = 600\n",
    "Local = []\n",
    "Global = []\n",
    "BlockID = []\n",
    "Blocktype = []\n",
    "Block = [] # which block\n",
    "Include = []\n",
    "Stimon = []\n",
    "SessionDate = []\n",
    "StimID = []\n",
    "sess = 0\n",
    "blocks = ['aa','aB','bb','bA'] # coded like this in the experiment\n",
    "for d,date in enumerate(dates):\n",
    "\n",
    "\n",
    "    \n",
    "    # get conditions\n",
    "    last_item = np.where((df.ItemID==3) & (df.date==date) )[0]#& ((df.StimDur+df.ISIDur)==soa))[0]\n",
    "    Ntrials = len(last_item)\n",
    "    if Ntrials>0:\n",
    "        \n",
    "        print('Day:',date)\n",
    "        SessionDate.append(date)\n",
    "        Local.append(np.zeros(Ntrials).astype(int))\n",
    "        for i,ind in enumerate(last_item):\n",
    "            if df.StimID.iloc[ind]!=df.StimID.iloc[ind-1]: # local deviant = transition of two stim\n",
    "                Local[sess][i] = 1\n",
    "        \n",
    "        BlockID.append(np.array(df.blockID[last_item])) # blockID\n",
    "        Block.append(np.array(df.blockType[last_item])) # which block (classical experiment: 0,1,2 or 3)\n",
    "        Blocktype.append(np.zeros(Ntrials).astype(int)) # xx=0 xY=1\n",
    "        for bl in np.unique(BlockID[sess]):\n",
    "            ind = np.where(BlockID[sess]==bl)[0]\n",
    "            Blocktype[sess][ind] = Local[sess][ind[0]] # blocktype defined by the first habituation trial in block\n",
    "        \n",
    "        Global.append((Local[sess]!=Blocktype[sess]).astype(int))\n",
    "        \n",
    "        # filter habituation trials\n",
    "        nhab = 50\n",
    "        Include.append(np.ones(len(Global[sess])).astype(int))\n",
    "        for bl in np.unique(BlockID[sess]):\n",
    "            ind = np.where(BlockID[sess]==bl)[0]\n",
    "            Include[sess][ind[:nhab]] = 0\n",
    "            \n",
    "        Stimon.append(np.array(df.StimOn.iloc[last_item])) # last stim onset times\n",
    "        StimID.append(np.array(df.StimID.iloc[last_item]))\n",
    "\n",
    "        \n",
    "        sess += 1\n",
    "    \n",
    "Session = np.concatenate([np.ones(len(Local[i]))*i for i in range(len(Local))])\n",
    "Local = np.concatenate(Local)\n",
    "Global = np.concatenate(Global)\n",
    "Blocktype = np.concatenate(Blocktype)\n",
    "Block = np.concatenate(Block)\n",
    "Include = np.concatenate(Include)\n",
    "StimID = np.concatenate(StimID)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### AUC context (blocktype) on global axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## first project entire sequence on local or global axis\n",
    "bins_around = 5 # how many bins to include pre and post max. decoding performance bin\n",
    "nperm = 100\n",
    "K = 10\n",
    "\n",
    "# get MUA\n",
    "MUA = []\n",
    "sess = 0\n",
    "for d,date in enumerate(dates):\n",
    "    # get conditions\n",
    "    last_item = np.where((df.ItemID==3) & (df.date==date) )[0]#& ((df.StimDur+df.ISIDur)==soa))[0]\n",
    "    Ntrials = len(last_item)\n",
    "    if Ntrials>0:\n",
    "        \n",
    "        print('Day:',date)\n",
    "        \n",
    "        # PSTH\n",
    "        Rb_last = np.load(os.path.join(out_path,'Rb_last_gauss_%s.npy'%date))\n",
    "        Rb = np.load(os.path.join(out_path,'Rb_seq_gauss_%s.npy'%date))\n",
    "        time_bins = np.load(os.path.join(out_path,'time_bins_seq_gauss_%s.npy'%date))\n",
    "        \n",
    "        _,nch,nbins = Rb.shape\n",
    "        # centre each channel\n",
    "        Rc = Rb.copy()\n",
    "        for ch in range(nch):\n",
    "            Rc[:,ch,:] = (Rb[:,ch,:] - np.mean(Rb_last[:,ch,:])) / (np.std(Rb_last[:,ch,:]) + 10e-6) # centre each channel\n",
    "        MUA.append(Rc)    \n",
    "        \n",
    "        \n",
    "        sess += 1\n",
    "        \n",
    "del Rb, Rc\n",
    "\n",
    "nbins = len(time_bins)\n",
    "\n",
    "# load coefficients\n",
    "Coef = np.load(os.path.join(out_path,'Coef_last_%s_%s.npy'%(exp,soa)))\n",
    "# load AUC of last item and get bin of maximum performance on the diagonal\n",
    "AUC_glob = np.load(os.path.join(out_path,'AUC_Global_last_True_classical_600.npy'))\n",
    "maxbin_glob = np.argmax(np.mean(np.diagonal(np.diagonal(AUC_glob,axis1=0,axis2=1),axis1=1,axis2=2)[0,:],0))\n",
    "GlobCoef = np.mean(Coef[:,1,:,maxbin_glob-bins_around:maxbin_glob+bins_around+1,:],2)\n",
    "\n",
    "# load test folds:\n",
    "Folds = []\n",
    "for i,d in enumerate(dates):\n",
    "    Folds_dict = np.load(os.path.join(out_path,'Testfolds_%s_last_classical_600.npy'%(d)),allow_pickle=True).item()\n",
    "    Folds.append([])\n",
    "    for j in range(K):\n",
    "        Folds[i].append(Folds_dict[j])\n",
    "nsess = len(MUA)        \n",
    "\n",
    "# Projections\n",
    "AUC_glob_block = np.zeros((nbins))\n",
    "RandAUC_glob_block = np.zeros((nbins,nperm))\n",
    "Proj_glob_all = []\n",
    "for sess,date in enumerate(dates):\n",
    "    # project MUA on axes within each recording session\n",
    "    Rc = MUA[sess][Include[Session==sess]==1,:]\n",
    "    ntr = Rc.shape[0]\n",
    "    Proj_glob = np.zeros((ntr,nbins))\n",
    "    for k in range(K):\n",
    "        test_ind = Folds[sess][k]\n",
    "        \n",
    "        # normalize coefficient vector:\n",
    "        coef_glob = preprocessing.normalize(GlobCoef[sess,:,k].reshape(-1, 1),norm='l2')\n",
    "        for b in range(nbins):\n",
    "            Proj_glob[test_ind,b] = coef_glob[:,0]@Rc[test_ind,:,b].T\n",
    "    Proj_glob_all.extend(Proj_glob)        \n",
    "# compute AUC from projections\n",
    "Proj_glob_all_arr = np.array(Proj_glob_all)    \n",
    "ind = Global[:-1][Include[1:]]==0\n",
    "    \n",
    "AUC_glob_block = getAUROC(Blocktype[1:][Include[1:]==1][ind].astype(float),Proj_glob_all_arr[1:,:][ind,:])\n",
    "for perm in range(nperm): # permutations for diagonal of temp gen matrix only\n",
    "    randind = np.random.permutation(sum(ind))\n",
    "    \n",
    "    RandAUC_glob_block[:,perm] = getAUROC(Blocktype[1:][Include[1:]==1][ind][randind].astype(float),Proj_glob_all_arr[1:,:][ind,:])              \n",
    "    \n",
    "AUC_glob_block[np.argmin(abs(time_bins - 3*soa/1000)):] = np.nan\n",
    "RandAUC_glob_block[np.argmin(abs(time_bins - 3*soa/1000)):,:] = np.nan\n",
    "\n",
    "threshold = 3\n",
    "alpha = .01\n",
    "\n",
    "time_bins = np.load(os.path.join(out_path,'time_bins_seq_gauss_%s.npy'%date))\n",
    "time_bins_last = np.load(os.path.join(out_path,'time_bins_gauss_%s.npy'%date))\n",
    "binsize = np.mean(np.diff(time_bins))\n",
    "\n",
    "t_sub = time_bins_last[maxbin_glob-bins_around:maxbin_glob+bins_around+1] + 3*soa/1000\n",
    "\n",
    "# plot trajectories\n",
    "colors = [['#229eba','#48e5ca'], # local standards\n",
    "         ['#f4884d','#f2dd20']] # local deviants\n",
    "plt.figure(figsize=(3,3))\n",
    "for c1 in range(2):\n",
    "    for c2 in range(2):\n",
    "        ind = (Global[Include==1][1:]==c2) & (Local[Include==1][1:]==c1) & (Global[Include==1][:-1]==0)\n",
    "        av = np.mean(Proj_glob_all_arr[1:,:][ind,:],0)\n",
    "        sem = np.std(Proj_glob_all_arr[1:,:][ind,:],0)/np.sqrt(sum(ind))\n",
    "        plt.plot(time_bins+binsize,av,label='Local %s, Global %s'%(c1,c2),color=colors[c1][c2])\n",
    "        plt.fill_between(time_bins+binsize,av-sem,av+sem,alpha=.2,color=colors[c1][c2],linewidth=0)\n",
    "# plot decoding window\n",
    "ymin,ymax = plt.ylim()\n",
    "plt.fill_between(t_sub,ymin*np.ones(len(t_sub)),ymax*np.ones(len(t_sub)),color='k',alpha=.1,linewidth=0)\n",
    "sb.despine(trim=True)\n",
    "plt.tight_layout()\n",
    "plt.savefig(os.path.join(out_path,'plots','Projections_Late_Global_%s_%s.pdf'%(exp,soa)),bbox_inches='tight')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "fig,ax = plt.subplots(2,1,figsize=(3,4),sharex=True)\n",
    "\n",
    "# plot AUC\n",
    "m = AUC_glob_block\n",
    "ax[1].plot(time_bins+binsize,m,'k',linewidth=1)\n",
    "ymin,ymax = ax[1].get_ylim()\n",
    "for i in range(4):\n",
    "    ax[1].fill_between(np.linspace(i*soa/1000,i*soa/1000+soa/1000/2,2),\n",
    "                        np.ones(2)*ymin,np.ones(2)*ymax,alpha=.05,color='k',linewidth=0)\n",
    "\n",
    "# stats\n",
    "p = clustertest(AUC_glob_block,RandAUC_glob_block,threshold)\n",
    "\n",
    "# stats:\n",
    "l = label(~np.isnan(p))\n",
    "ymin,ymax = ax[1].get_ylim()\n",
    "for j in np.unique(l):\n",
    "    if p[l==j][0]<=alpha:\n",
    "        print(p[l==j][0])\n",
    "        ax[1].plot(time_bins[l==j]+binsize,np.ones(sum(l==j))*ymax,color='k',linewidth=1.5)\n",
    "\n",
    "\n",
    "for c1 in range(2):\n",
    "    for c2 in range(2):\n",
    "        ind = (Global[Include==1][1:]==c2) & (Local[Include==1][1:]==c1) & (Global[Include==1][:-1]==0)\n",
    "        av = np.mean(Proj_glob_all_arr[1:,:][ind,:],0)[:int((soa/1000*3-binsize-.05)/binsize)]\n",
    "        sem = np.std(Proj_glob_all_arr[1:,:][ind,:],0)[:int((soa/1000*3-binsize-.05)/binsize)]/np.sqrt(sum(ind))\n",
    "        ax[0].plot(time_bins[:int((soa/1000*3-binsize-.05)/binsize)]+binsize,av,\n",
    "                   label='Local %s, Global %s'%(c1,c2),color=colors[c1][c2])\n",
    "        ax[0].fill_between(time_bins[:int((soa/1000*3-binsize-.05)/binsize)]+binsize,av-sem,av+sem,alpha=.2,color=colors[c1][c2],linewidth=0)\n",
    "ax[1].set_xlim(-soa/2/1000+binsize,soa/1000*3-binsize-.05)\n",
    "plt.xlabel('Time from sequence onset')\n",
    "ax[1].set_ylabel('Blocktype AUROC')\n",
    "ax[1].set_ylim(None,.6)\n",
    "ax[1].set_yticks((.5,.6))\n",
    "ax[1].plot(time_bins+binsize,np.ones(nbins)*.5,'k--',linewidth=1)\n",
    "sb.despine()\n",
    "plt.tight_layout()\n",
    "sb.despine(trim=True)\n",
    "plt.show()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learning during habituation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Coef = np.load(os.path.join(out_path,'Coef_first3_classical_600.npy')) # load coefficients of glob axis\n",
    "\n",
    "from scipy.stats import ranksums,ttest_rel\n",
    "# get global deviance response on every trial during habituation\n",
    "block_ax = 0\n",
    "nhab = 50\n",
    "\n",
    "# time vectors\n",
    "time_bins = np.load(os.path.join(out_path,'time_bins_gauss_%s.npy'%dates[0])) # time vector\n",
    "time_bins_seq = np.load(os.path.join(out_path,'time_bins_seq_gauss_%s.npy'%dates[0])) # time vector\n",
    "binsize = np.round(np.mean(np.diff(time_bins)),3)\n",
    "time_bins = time_bins + binsize\n",
    "\n",
    "# load PSTH per session\n",
    "MUA = []\n",
    "sess = 0\n",
    "for d,date in enumerate(dates):\n",
    "    # get conditions\n",
    "    last_item = np.where((df.ItemID==3) & (df.date==date) )[0]#& ((df.StimDur+df.ISIDur)==soa))[0]\n",
    "    Ntrials = len(last_item)\n",
    "    if Ntrials>0:\n",
    "        \n",
    "        print('Day:',date)\n",
    "        \n",
    "        # PSTH\n",
    "        Rb = np.load(os.path.join(out_path,'Rb_seq_gauss_%s.npy'%date))\n",
    "        \n",
    "        ntr,nch,nbins = Rb.shape\n",
    "        # centre each channel\n",
    "        Rc = Rb.copy()\n",
    "        for ch in range(nch):\n",
    "            Rc[:,ch,:] = (Rb[:,ch,:] - np.mean(Rb[:,ch,:])) / (np.std(Rb[:,ch,:]) + 10e-6) # centre each channel\n",
    "        MUA.append(Rc)\n",
    "        \n",
    "        sess += 1\n",
    "        \n",
    "del Rb, Rc\n",
    "MUA = np.concatenate(MUA)\n",
    "nbins = len(time_bins_seq)\n",
    "nsess = len(dates)\n",
    "K = Coef.shape[-1]\n",
    "Proj = np.zeros((nsess,4,nhab,nbins,K)) # shape is sessions x blocks x trials x bins\n",
    "for sess in range(nsess):\n",
    "    coef = Coef[sess,block_ax,:]\n",
    "    for bl in range(4): # separately for each block in session\n",
    "        ind = np.where((Session==sess) & (Block==bl))[0][:nhab]\n",
    "        # project those trials onto global axis\n",
    "        for b in range(nbins):\n",
    "            for k in range(K):\n",
    "                Proj[sess,bl,:,b,k] = coef[:,k]@MUA[ind,:,b].T\n",
    "          \n",
    "\n",
    "    # PLOT\n",
    "# which bins to average:\n",
    "soa = .6 # sec\n",
    "binsize = np.round(np.mean(np.diff(time_bins_seq)),3)\n",
    "b_start_glob = int(np.round((abs(time_bins_seq[0]) + 0)/binsize))\n",
    "b_end_glob = int(np.round((abs(time_bins_seq[0]) + 3*soa - .2)/binsize)) # until 3d stim + soa - .2\n",
    "\n",
    "Proj_cv_av = np.mean(Proj,-1) # average across cv folds\n",
    "\n",
    "window = 1\n",
    "proj_xx = np.reshape(Proj_cv_av[:,[0,2],:],(nsess*2,nhab,nbins))\n",
    "proj_xx = np.mean(proj_xx[:,:,b_start_glob:b_end_glob],-1)\n",
    "proj_xx_smooth = np.array([np.mean(proj_xx[:,i:i+window],1) for i in range(1,int(nhab-window))])\n",
    "proj_xx_smooth = np.concatenate((proj_xx[:,0:1].T,proj_xx_smooth),axis=0)\n",
    "blocks_xx = np.repeat(np.array([0,2])[np.newaxis,:],nsess,axis=0) # block label\n",
    "blocks_xx = np.reshape(blocks_xx,(nsess*2))\n",
    "\n",
    "proj_xy = np.reshape(Proj_cv_av[:,[1,3],:],(nsess*2,nhab,nbins))\n",
    "proj_xy = np.mean(proj_xy[:,:,b_start_glob:b_end_glob],-1)\n",
    "proj_xy_smooth = np.array([np.mean(proj_xy[:,i:i+window],1) for i in range(1,int(nhab-window))])\n",
    "proj_xy_smooth = np.concatenate((proj_xy[:,0:1].T,proj_xy_smooth),axis=0)\n",
    "blocks_xy = np.repeat(np.array([1,3])[np.newaxis,:],nsess,axis=0) # block label\n",
    "blocks_xy = np.reshape(blocks_xy,(nsess*2))\n",
    "\n",
    "v = np.linspace(window+1,nhab,proj_xx_smooth.shape[0]-1)\n",
    "v = np.concatenate((np.ones(1),v))\n",
    "\n",
    "## test difference\n",
    "# average the two same blocks from every session\n",
    "proj_xx = np.mean(Proj_cv_av[:,[0,2],:],1)\n",
    "proj_xx = np.mean(proj_xx[:,:,b_start_glob:b_end_glob],-1)\n",
    "proj_xy = np.mean(Proj_cv_av[:,[1,3],:],1)\n",
    "proj_xy = np.mean(proj_xy[:,:,b_start_glob:b_end_glob],-1)\n",
    "d = proj_xy-proj_xx\n",
    "res = linregress(np.arange(d.shape[1]),np.mean(d,0))\n",
    "    \n",
    "# PLOT\n",
    "# plot xx\n",
    "plt.figure(figsize=(5,4))\n",
    "av_xx = np.mean(proj_xx_smooth,1)\n",
    "sem_xx = np.std(proj_xx_smooth,1)/np.sqrt(nsess*2)\n",
    "plt.fill_between(v,av_xx-sem_xx,av_xx+sem_xx,color='#165aa9',alpha=.2,linewidth=0)\n",
    "plt.plot(v,av_xx,color='#165aa9')\n",
    "\n",
    "# plot xY\n",
    "av_xy = np.mean(proj_xy_smooth,1)\n",
    "sem_xy = np.std(proj_xy_smooth,1)/np.sqrt(nsess*2)\n",
    "plt.fill_between(v,av_xy-sem_xy,av_xy+sem_xy,color='#ee8c22',alpha=.2,linewidth=0)\n",
    "plt.plot(v,av_xy,color='#ee8c22')\n",
    "\n",
    "plt.plot(v,np.zeros(len(v)),'k--',linewidth=1,alpha=.4)\n",
    "plt.xticks([1,10,20,30,40,50])\n",
    "sb.despine(trim=True)\n",
    "plt.ylabel('Block type trajectory')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(os.path.join(out_path,'plots','LearningBlock_%s_%s.pdf'%(exp,soa)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats\n",
    "\n",
    "# stats\n",
    "print('Slope',res[0])\n",
    "CI = 1.96*res[4]\n",
    "print('CI slope',CI)\n",
    "z = res[0]/res[4]\n",
    "print('Z score',z)\n",
    "p = scipy.stats.norm.sf(abs(z))*2\n",
    "print('P',p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Blocktype effect for matched trials and post-global deviant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3373: RuntimeWarning: Mean of empty slice.\n",
      "  out=out, **kwargs)\n",
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:163: RuntimeWarning: invalid value encountered in true_divide\n",
      "  ret, rcount, out=ret, casting='unsafe', subok=False)\n",
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:234: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  keepdims=keepdims)\n",
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:195: RuntimeWarning: invalid value encountered in true_divide\n",
      "  arrmean, rcount, out=arrmean, casting='unsafe', subok=False)\n",
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:224: RuntimeWarning: invalid value encountered in true_divide\n",
      "  ret, rcount, out=ret, casting='unsafe', subok=False)\n",
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3373: RuntimeWarning: Mean of empty slice.\n",
      "  out=out, **kwargs)\n",
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:163: RuntimeWarning: invalid value encountered in true_divide\n",
      "  ret, rcount, out=ret, casting='unsafe', subok=False)\n",
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:234: RuntimeWarning: Degrees of freedom <= 0 for slice\n",
      "  keepdims=keepdims)\n",
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:195: RuntimeWarning: invalid value encountered in true_divide\n",
      "  arrmean, rcount, out=arrmean, casting='unsafe', subok=False)\n",
      "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:224: RuntimeWarning: invalid value encountered in true_divide\n",
      "  ret, rcount, out=ret, casting='unsafe', subok=False)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUMAAAEYCAYAAADGepQzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABuEklEQVR4nO2dd3xUVfr/3+dOTya9UwOEkkDogrCAioKArquiyFpA1wW7u+q6YhdX1y+6uqz956q4WFHBjgqCBbCiSO+9BdKTyfS55/fHTYZ0JiGAgfN+vfIic+fee84dkk/OOc/zfI6QUqJQKBQnO9rx7oBCoVD8FlBiqFAoFCgxVCgUCkCJoUKhUABKDBUKhQJQYqhQKBSAEkOFQqEAlBgqFAoFcBgxFELYhBA3CSG+EkLsE0LsFkIsFEJMFUJYj1UnFQqF4mjToBgKIU4DfgKygUeBM4Czgf8AvYHlQoiRx6KTCoVCcbQxN/Le74FhUsqyWsfXAR8LIeKA+4DFR6tztVB1g79dxPHugEJxpIiGapOFEFFSSvcx7k9jKDH87aLEUNHqaWzNcLkQoucx64lCoVAcRxoTw6eAr4QQVx+rzigUCsXxosFpMoAQohvwP2AbMOU4T5vVNPm3i5omK1o9jYohgBBCA17GiCTnVx2XUvY+ul2rQ6MdLS8vb/aNY2Jimn3tydZuAygxVLR6GosmI4SwAA8B44A7MEaICkVE/Pzzz6lms/lFoBcqwV9xfNGBNcFg8M8DBgw4WN8JDYqhEKIP8BpQCAyQUu4+On1UnKiYzeYX09PTs1NSUoo1TVPLHIrjhq7rIj8/PycvL+9F4Lz6zmnsr/W3wFvAGUoIFc2kV0pKSpkSQsXxRtM0mZKSUooxS6mXxqbJI6WUP9T3hhAiWkpZcaQdVJzwaEoIFb8VKn8WGxwANvhGlRAKIT4XQiRWHRdCDARWtGQnjxWhkt2ECrchQ4Hj3RWFQvEbI5JF7R+AX4QQQ4UQtwEfAvcc3W61PLqnGOkvQ4Y86KV7j3d3FK2Ytm3b5n7zzTdRtY9//PHHMV27dm12oYIQYsD+/fsbDWo2xpNPPpmUkZGRO2zYsK5ff/111KWXXtqhufc6HBMnTuy4ZMmSOp9Bc5k0aVKHW2+9tc2R3GPWrFkJgwYN6t7c6w8rhlLK+4C/AEuAWzGCKW83t8HjhazIP/R9sALdXXQce6NQtDxvvPFG0r333rt36dKlm1euXOnYv3//UXOWWrJkSeyJts3wYcVQCDEGoxrl38BK4E0hxFH7i3M0kH4XUq85NdYr8pF66Dj1SHE8eeqpp5LatWuXW1paqpWVlWmdO3fu+fTTTyc1dLy+ezz55JOpOTk52VlZWT1nzpxZ55zCwkLTH/7wh05du3bt2a1bt5xrr722XSBg/AwuXrw4unfv3j26du3aMycnJ/vDDz+skfi5a9cuc9euXXv+85//TKl931WrVtmGDh3atW/fvj3atGmTe+aZZ3Zxu93i6quvbr9q1arohx9+uO3tt9+e8cgjj7RZvny586KLLsoEeOONN+J69+7dIzs7O6d///49vvjii2iAW2+9tc2wYcO6duvWLecPf/hDp+ptBQIBLrvssg7dunXL6dmzZ/bYsWM7l5aWajfddFPbgwcPWiZNmtR58eLF0YsWLYoeOHBg9969e/fIyMjInTBhQkeAjRs3Wtu3b99r8uTJ7XNzc7M7dOjQ67///W8CQFFRkTZu3LjOmZmZvQYNGtR9y5Yt9qp233zzzbh+/fr16NWrV3ZGRkbuX/7ylzZgjLy7d++e069fvx7du3fP8Xg84q9//Wub9u3b98rNzc1+77334pvwY1CHSIbkLwOTpZQLAYQQtwLLgdQjabjyXgKYBayRUv7rSO/XELq3tvEOIIPopXswJXQ8Ws0qfqPcdNNNhQsXLoy94YYb2vn9fjF48GDXjTfeWAjQ0PHa2O12fd26det37Nhh6d+/f86wYcNqBBSnTJnSPjExMbhx48a1Pp9PjBo1Kuv+++9Pv//++w9ccsklXZ566qmdEydOLF2yZEnU1VdfnXnOOeesA9i+fbtl0qRJnW+77bb91113XZ3pyzPPPJNy+eWXF15//fVFPp9P5ObmZr/99ttxL7300u7Vq1dHXXfddQevuuqq4vbt2/vfe++9hHfffXfH6tWrbQ888EDbb775ZmN6enpo+fLl9rFjx3bbunXrGoC9e/daN23atNZisdRoa9GiRc7vvvsuZsuWLWs1TeO6665r++OPPzqeeuqpvfPmzUucPXv2thEjRrh///vfd3rggQf2nXvuueWlpaVaZmZm7pIlS/JTU1ODe/bssY0ZM6bsf//73+5XXnkl/q677mo/ZcqU4ttvv72t3W7Xt23btiYvL8/cv3//nFNPPdWl6zozZ85Mmz179vbc3Fzfjh07LFlZWb2nTZt2EGDLli2O9evXr+7WrZv/tddei//4448TVq1atc7pdOqjR4/OOpKfi0jEsI+UsnrlyRNCiCO27RJCZAPPAKcCa470fg0hpUT66hFDQAZchEr3IO0dQQ8hAx7QNIQjEaPwRtGSpF329oCj3caB1yf8HMl5//vf/3bm5ub2tNvt+qpVq9Yd7nhtbrjhhnyAzMzMwIgRI8o+++yz2L59+3qq3v/qq6/ivv766w2apuFwOOTUqVPzn3rqqbRx48aVaZrGxIkTSwGGDx/u3rRpU7id888/v2taWlrgmmuuqXcd59lnn93z/vvvx95zzz1pmzdvtufn51vKy8tNjT3rJ598Epufn285/fTTw+tpQgjWrl1rA+jfv39FbSEEGDhwoMdkMtGnT5/skSNHlk6YMKH4jDPOqFOS+8477+x455134qZNm5a+ceNGu9fr1crKyrTU1FTMZrOcMGFCKcCpp57qLi0tNQMsWbIk5rHHHtutaRpt2rQJjhkzphhA0zQWLFiwZc6cOXGvvPJK4saNGx1SSsrLyzWA9PR0f7du3fwACxcujDnnnHOKExISdIArr7yy4Jlnnklr7LNojEjE8D7gJiHEU1LKmwCklL82t8Fq3IAxKtzVAvdqmEAFyBCE/OApAkcimA4tpUhfKaGirTUuEZ4StMTOShBbmEiF6liwd+9ei9/vF36/37Rjxw5LTk6Ov7HjtTGbzeEFMyklFoulxgJa7fU0XdcJBoPCYrFIY0J0iJ9++snet29fL8DMmTN3zpgxI2P69Olp06dPP1C73fPOO69zMBgUF110UdH5559fumfPHuvh1u5CoZAYOnRo+SeffBKuINuyZYslMzMz8M477yRER0fr9V2XnJwcWr9+/dqFCxc6Fy5cGHv55Zd3mTp16oH777+/RgXHoEGDuvfs2dMzZsyY0ssvv7x4xYoVTimlALBYLNJkMrRaCFHjc6n+fdXnWVZWpvXu3Ttn7NixJSNGjCi/7rrrChYuXBhfdW5UVFS4r7XvV/3/pDlE8tv+u8p/hx1JQ7WRUt4opXy1Je9Zbzu+CggFYNk/YelDsOjv4Clu/JqgB73k6Gq04vjh8/nExIkTO0+bNm3f3//+930TJ07s7PP5REPH67vHf//732SAzZs3W5cuXRo7duzYGtOP4cOHl/373/9O1XUdj8cjXnzxxZTTTz+9rHfv3l4hhHzvvfdiAZYuXRp19tlndw+FQgLgjDPOqJg9e/b2f//73xk//fSTvXa733zzTewDDzywb8qUKcVCCFatWhVddW11zGazDAQCAuDss88uW7JkSeyKFSvsAHPmzInr379/T7fb3ejv/5tvvhk3dOjQ7meddZbriSee2DdhwoTCVatWRQGYTCbp9/tFfn6+ae3atdFPPvnknsmTJ5fs3LnTsnv3blso1Ph6/MiRI8teeuml5FAoRH5+vunzzz+PB1izZo3N5XKZZs6cuffSSy8tXbBgQYzf7xfBYLDOM44dO7bsww8/TCwoKDCFQiFee+21etd3I6XZYfxjgRBiKjAVICcnh7Vr1zb5HtJfDqtfBXsCDL8PNn8E3z0Kp02vMUKs7zrddQDN2exRt+I3ys0339w2NTU1cOuttxYAfPTRR/F/+ctf2kopqe/4888/v6f2Pbxer5aTk5MdCAS0GTNm7Ordu7dv165d4R+oF154YdeUKVM6dO/evWcgEBCnn3566SOPPLLfbrfLOXPmbL3lllva33nnne0sFot88803t9rt9vCopk+fPr5bb711/xVXXNH5119/XV/9vbvvvnvvxRdfnBUfHx+02+36Kaec4tqyZYutdv9OO+20ioceeqjtqFGjuixcuHDrzJkzd/7xj3/sLKXEbDbLt99+e0tsbGy9I8IqLr744tL58+fHdevWrWd0dLQeGxsbnDVr1k6AMWPGlFxxxRWdn3nmmZ033HDD/r59++YkJCQEExISgv369XNt2rTJ1qNHD19D937sscf2TZo0qUOXLl16JSUlBbp37+4BGDx4sOfMM88sycrK6hUbGxvMzMz0denSxbthwwabzWarMfK75JJLSletWuXo379/dmxsbKhnz57uwsLCuvP9CInEteYXKWV/IcQKKWW/5jbUyP1fIYIAysCBA+Xy5csbfL8+FxcZ8hPa/R18eTeM/jdYokBK+Ok/EJ0GPf8IQIwzusH7ajFt0KIa/oPTmHuMDLjRyw8gbLFo0XXvcaK71qxcuXJHnz59Clq6MYWiuaxcuTK5T58+mfW9d0IvikmfC7YtgA7DDSEEEAL6XAW7lkD+4eM2evk+9PL9ddaADtt2wEOoZBcy6EavyEP6VfWiQvFbJhIxrIqQ/Zb2Q4kI3VMEO7+GzmfXfMMWB6fcBD89A3t/IFS0Df+Gj/GteJVQPWuFursAvWgLMtjgqL8GUg+hl+02AjdV9yjPO6JnUSgUR5dIKlB+V/3flkZKeeVRyzHctQQSOkF0ndxVSM6GAdfB9kV4lvwLvWQnIPEsvIfArm/r9jPoJVS0NaIRnl62t079swx5jdQdhULxm+SwARQhhBMjveZsIAR8BPxTShnZMOk4IQNe2PUNtGtEw9N6Q1pvnNXWDM0dhuJZ/CCmpK5otUVUhgiV7MSU0AlhcdR7S9110Aja1PeepxhTA9cpFIrjSyTT5BeBtsAtwN+BHsCTR7NTLYFevg8KN0KbU5p0nSkpC0v27/F+93T964SVgihDddPPdG8puju/7jVVl3pLkbLRAJ5CoThORCKG/aSUl0kpF0spvwD+CIw4yv06YuS2hcbIz3woVUvY4hHW2MNea825EOkrJ7jty/pP0APoxdtrCKL0u9DL9h3mzjqyvtJAhUJx3IlEDA8IIZKrvY4GftPpElLqxhS57anhY5ozA1NcW0zx7dGcGY1eLzQT9lNvwPfL/9DL99ffRshPqHArekUhuusgoZLdGNssAHoISneCu+7HJL2lzX6u3wLKC1JxohJJ0nUe8LMQ4h0gCPwBQyCfBJBS3nwU+9csZOkeKNkOaX0AENYYtKiwP63xvR5sdEprSuqCtc9EPF8+hGPUQ2iOhHoaCqJXVIsSB72w9TPYthAsDvCXQ7uh0Osy0IyPWgYqkHoQof2m890bRPeVY6r2WSoUJwqRjAzXYjjXlGOk2bwFfImxUVS9jh7HG2OK3KeywkSgxaTXOUdzpiIszkbvY+02FnPn03HPv43gvsOYe+9fbpT6le2G390JZ/0LznoCyvfCuur2jxLpaZ2jQxn0IQ9TyngycLKau37zzTdRbdu2zW3u9ZdccknH999/v9kZ/0888UTy//3f/9WTGtIyHPaDl1JOP1qNHy3kjkXQxcgt1OwJiAbK7rT4DsY6nwbCbEdYHMiAG+l3hc+x9boYU1I3vN8/gyk1G/vAqxG2auuOIT+sfg0Orob+10JKzqH3rNEw4AZYPM1I/I5tb/TPWwj1VKT81tFdeaAHj3c3FA1QZe56/fXXFz355JNJR9PctTnMmTNn55Fcv2zZMmfPnj2PWn5aJOauq4UQq2p/Ha0OHSl68Q5jvS61DyAQjUzphBCY4tpiTumBKSETzZmGFp+JsMfXOM+c0Yfo3z+JsMVS8dHN+Dd8jF62j8Du7+GruyHghjMerimEVdjjoMcFsPbN8CEZCiB9hy+nk1JHdxcQKtpGKH+j8VW42XhdsgvdlW+kEB0jdFcdE5VWiTJ3jczcFWDGjBkpmZmZvXr16pX91FNP1fAwveOOO9JzcnKye/TokXPWWWd12bFjh2XVqlW2hISEPl6vVwAEg0HS0tJ6//zzz/ZBgwZ1nzVrVgLAtGnT0nNzc7O7deuW0759+16zZ8+Or+rP+PHjM4cNG9a1U6dOPQcMGNB9x44dltmzZ8d/8cUX8c8//3zaI488clRGh5EMyW+s9r0VOB84XNj0uKFv/gQyBoDJggx60cvz0GLSEeY6tez1IoRAi8kg5HOBPDQKEmY79oFXE+p8Bv417+Bf/6GRh9jjImgzyCjza4iOp8PG940pc0xbo5+uPGRSBrXtnMLP4Sk2xKeWQ3dVjEb6K4x8RvdBhMmBiEpE2JyGL6PfBbqOsMdG/NyHQw/5CZXuwZwa1+x73Pfar23+36ebwtGrD+49Yz3AH/7xZXbVsWvGdtv/4OV99+Vc+0HvwnKfBaBb21j3kkfHrL/m6e86vv/d7nAwb/nMc1b9tLkg6rpnfsiqfu3h+qHMXSMzd/32228djz76aJsVK1as7dChQ7D6tPvpp59OWrt2bdTKlSvXWywW/vWvfyVPnjy549dff70lKyvL++abb8ZfddVVxfPmzYtt166db8CAAeG/2ps2bbJ+9dVXsd99990Gp9MpX3jhhYSHHnqozaRJk0oAfvzxR+fKlSvXJSYm6iNHjsz6z3/+k/Lvf/973wcffFDSs2dPz5133tnwYv8REMk0+evqr4UQX2Dsqfzw0ehQU6jtKiNDftjxBeRMQK8oRNhikAEXeumeSn/CRgSrGkIzozlT6o0kmxI74xhxR/h1uSuCmmOTFTqdaQRX+l4d7qt05SFiDkW2q4xoZUU+Mhj5bECGPMjyvcaqbvXn8BYbz90CwRpZvJ3A+g/Qi7ZhTmnenjsPXt53X31iVZ/P4brn/1Bn9vH/bhyy8//dOKTGVKt9SnTphUM7NtknUZm7Ht7c9fPPP48dPnx4WYcOHYIA119/fcHXX38dBzB//vy4lStXRufm5uaA4dfo8Xg0gKuuuip/9uzZSVdddVXxK6+8kjx58uQaaRWVLtXbX3jhhaQtW7bYli9f7qyoqAg/w5AhQ8oTExN1gD59+riLiooafb6WojlGDUnAEe1i1VLoFQXo5XmGiEhJaPvX4CtHJmRByI8wG0smMuRFVjQtG0g4khpca2wWnc6Cvd9DtfXIqimwXr6fUOluQvkb0Et3HVYIpZRGtHzPd4ZXY0Pn6QH0svpTg5qCDPoIFu0kuOdHtPg6M6lWSZWJa1lZmWnHjh2Wwx2vTUubu1ZNoWfOnLlT0zQ5ffr0er3jzjvvvM4vvPBCSseOHf233377gZycHHek5q4bNmxYV/X17bffrj/llFM8AA2Zu9Y2T63+jKFQSNx88815VfdbuXLl+mXLlm0AuPLKK4t//fXX6F9++cX+ww8/xFx11VU1om5Lly6NGjZsWHZZWZlpzJgxZX/7299q/JBWtyyr7ENko5gjpKlrhquBrcCco9+1yNDd+YQKNhIq2AjbPkN2PB29oghhrzmd0z2FyCYs/gshENGRehlqaI5kTPGZmOIz0ZzpCK3W75EtDtL6GfXS1ZCBCnR3AdJbUmNa3hB62T48X9wHP/4HdiyGL26DioMNni/95REbTNR7vZSESnYS3PsDhAKYkro0+16/FZS5a2TmruPGjStbsmRJ7NatWy2VzxxeGx01alTp7Nmzk4uKijQw1vomTpzYCSAqKkqed955RVdeeWXm2LFji2NiYmqI7eLFi525ubkVDzzwwIFx48aVv//++wn1PUNjz3Q0aOqaoQTypZTrj1J/moceAG8x7P0eOWQaAhDWSu9CTyE4koAQsqIQEdO4wFVVlQiTFc0Rj/QW14gu14eIaWOszckQmO1o1mikIwG9eCcyWM3sp9OZ8MsLRqS7CVsKyICH4P4VBHd9T2j/Cqw9LyLUfqRxjy2fGsI44v4GzGol0l2EiG080bwh9LK9SHc+we1LENHJDdZktyaUuWtk5q6DBg3yPPDAA3vOPPPM7k6nM9S3b9/wmtAtt9xSsHfvXuspp5ySLYSgTZs2/tdff31H1fvXX399wamnnpr67LPP1rGBuuqqq4o++OCDhC5duvS0WCxy+PDhZaWlpabi4uJGfynGjh1bevvtt3cAeOSRR1rcBqpBc1chxKiqHfEavFiI0VLKBS3dqfqobu4q9SDen/6Lud1ghNkYgZUvfQopJXrHkWhRyQh7DLjyjGivI9FIcYnvjCm5a501tJiYGGTIj166Bxmo/P/WLGi2WLBEVaaU1J2OSj2Ey+UGkyVs1yU0C1pCJsJkRepB9OJth6o2pISv7oEeF0LGgEZNZcEQIt/qtwnu/gFTSnfMbU/B0vl0hNV5aK1SSlj+NFidhk9jvWiYkrshNFP4eSNBL89Dd+cTPLAGz6IHsHQ+E/vQv9S3ZqjMXRWtguaau44UQrwvhBgnhAirhxDCKoQ4WwjxMTC6hfsaMXrxNjyL7kMv349/+zfIfT+itznVGC2VbjcSnZ3pcM5/oevv4ceZ0MDaoSFaOw4JIRj1x55Cw5dQDyKD3homCzLgQboLQffDmteMERqV63QlO5FSGoGYmGrLq0JA9/ONyHIj6zxSD+Jb/Q7uz6ehxbbFeeGLRJ05HWuPcxHWWoniQhhBmQMrYffShj6tJpcB6hUF6O58ZChAYON8TB1+h7nLGU26h0LRmmhQDKWUdwKPADcBxUKITUKIrRhVJ38BHpZS/u3YdLMmQjMTdfYMtMQsKj79O/5fX0PveRkIgbbra8SK/wepvatOhg4jYPAtYLKie4rqrKHppXvQg15CxduNrUPrbC4vwWQHvwfpKTVSXnwV4EgwpqjufEjoDAEPlO8zIsWVoiusToS9WilfxgBDQPf9WO+zBfNW4/70dkIH1xE19nFsuRPqCmDV52Cyg2YxXLwH3wLr5hgJ4OV1s0ukpySizxZA+srRXcaadujgOoJ7f8bkTDVG1BFG5BWK1kaja4ZSyh+AsUKIOCALI8ttm5TyuNeTCaFh6TISc6fhhsNMyIrY9Q0EXHDaP4ypcXXiOhoBh5h26JYoTAmZAOjuQiTloAfxfv+sscdyKIC19yVYu55drT3AZkxrBUZwJrj7e+PNATeAZoK9P8Dmj+G06ejufIQjzlh7jE4h5C0FdEOc+06B7x8nlJKJKbGzsSa461v8mxcgvSXY+lyGOXN4o6lAwhKNFt8BpEQv3oqM62g896YPYOnDENcB+v4Jooz8VBnyIAOew675yaCPUOlu4/uAB+8Pz6LFtgn/AWlImBWK1k5ECWiV4veb2fMWIFSwkcDOpegFm4x9j/vfAJ3HIGLr1iGHMdlg1SvIhM6EiraD2YYs2oxn/ZvYh91G9NjHjHuX7IKAB6kHCe7+AXOHoXWESeohLO1PxRuXYwihMEObwbBlvuGY0/F0dNdBTHHtECYLWlTSIWOIxCzoPQnPovvB6kS6izCl98ba80LMbQeG1/YaQmhWtLg0Y19nAcKRjHTtB3s89J4MPS89tAvg8AeMskAM4TfFtWvwvlLq6KW7QIaQoQCerx5Bi+tA6MAarH0vR5gdLZtupFD8hmid1ikYKSZ62V5MKdlYcy/G5ReIqPj6T9YsRkCj3VDYsww2fYDsMd6Y3n73KKZuZyNMh1JhTPFGor3uLsK/5h1C+1ZgG3xdWKRkwIP78zuxD7sNzIkIayym+PbGOlvvScbUuf1wpK8UGUhGWOyIqCRjE/uqfVHaDSG662noZXvRYjIirxQRJrT49ghvCBlwG3Zh9njDLqwqyGOyGEEaXxmsegUG3mD021dWryktGCk0euluY21UD+Fd9m8wmTGlZiP9LjR7HMLmNKL0CsUJSKvdHc/S+QxsuZdgbj8IYYlCmKMQ5jppWQhLNKakLKP2OCYD+v4ZvCVGsvLX90FaH6zZ59XbhhaVSNToR9A9RXgWP1BZU1yGZ+njmJK7Yoo3jBc0p1GyKaISEUk9YPi9xmgR0CtzAIVmQqs1dRcmC6aEzCYIoRktti0y4CVUtL1yjXMXVBysc28Aek6EgnVGgjYAsl4DWmNEuNtYIgB8P7+E9JVj7XMZ5vTe2PpNAmFCOBIRFiWGihOTSJKu/1DPsSuOTneahhadYvySmqxosW2MYEI1hMmBFtfemE4CWlQSIraDEX01WWDoHYbXYNX5ZkelI82hj0VYHDhOvxtrj/MQJgven15Ec6ZiO2WK8b71UP2vEJoxAnQkG9HsoNdIeq40UxBRSSCaU1kkKl26o9BLd6GX7zXyFyuj27qnEExm6mS4mO3Q/QJY+1b4kAxUVG59Wnmttwy9cAvSZywD+7csJLR/Ffbf3YpesoPAzmXGNN8eV7+no0JxgtCgGAohfi+EuBB4QghxYbWvS4B/HrsuNowWlYi5TT8snU7DnJqNObkbpoTOaI4ktKgUtMROCM2Erku8fmN6qsW2Mdb3AKKrTDgEWmw7YwQZ3xFTSnc0RxJV4iI0E+Z2xl4q9sHXYT9lanjtTKvliiPs8YaRa8UB2PY5YJTdVd1HszdNUIQlGi2xC4T84ZGbDPlh5f/g15eMk6REuotrWotV0fF0cO2H4q3hQ7q7gFD+BoIH1qKX7kSGjOBIqHAL/hWvYj9tGiAJbvsKkIaA2+PqVPWcaEyaNKlDnz59egSDhyqBgsEg/fr163HzzTfXKEHduHGjNSoqql9997n11lvbTJo0qVlegkfqiQhw2WWXdWjbtm3uTTfd1PZoegAWFhaaTj311G4tec+uXbv2/Pjjj5vteQjG/+Ott97a5JLhxkaGfTHSalIr/636+jPwWDP62OKYkrtjimtfI8IpLA60mHTDvFUIQrrOlv1lbM8rx+sPGVPT2DaER1HChCm+Y41Rj9DMaLFtMCVm1QkYVI/GCrMDYa2ZOC00E5otDnqMN3IPA+4aa3UiOvLRoRadakSbK/LDlSx6+X7cn00DX6kRKDnwK/z8LDLorpzC1hodamboei5s/KDmcRkibIEDSJ8Lz5LHsA26BlNcO2MKXrQVc9tT0KJT0SyOE36K/Pzzz+92u92mu+66K1yuc9ddd2WYTCb5xBNP/Gadmmrz5ptvpnzzzTcbnnrqqb3Lli1zHq7srrkUFBSYVq1a1XjlQCuisTzDf0gpzwDukFKeUe1rFPD0sevikVFU5sMXCBHUdfKKDUERthhMiV3QnOmYEjsjbPX/IRIWO1pCJ4RWfwRVRNVv0CqiEiGmDaT3h6JNGCVxRq260MxoUY3/oRaaFVNCJzRnGrq3FOkrNVKApCSYtxpLl5Fwyk1GlDg5xxj15a9FesvqH3l2PB1KttUYHVZHhgJ4vpmBud0gLB1/Z1TjFGzE3GEoppg2CFt0zVzJVk5DvoVRUVHy9ddf3/bss8+mf/vtt45ly5Y5XnnllZR33nlnm9lcN9YopeSSSy7pmJ2dnZObm5u9aNGiOsKwfPly++DBg7t169Ytp3v37jnV/RFnzpyZlJWV1bNbt245p556arctW7bUWOf5/PPPnW3atMlduHBhnfu++eabcf369evRq1ev7IyMjNy//OUvbQAGDBjQXUrJ2LFju95+++0ZtT0A6/MgBBg0aFD30aNHd+nSpUvPhx9+uIZv4a5du8y/+93vuubk5GTn5ORkV7U1adKkTJ/Pp/Xo0SMnGAwyc+bMpCrPxDZt2uTOmDEjBQwH7jPPPLPLqFGjulR5OFbVXf/888/2Pn369OjWrVvOuHHjOlcX7sY8D6t7MBYVFWnjxo3rnJmZ2WvQoEHdt2zZUjd4EAGRRJP7CyGsUko/gBCiPfAmMKw5DR5LpJQUlB9KsHZ5ggSCOhazhjDbIgpcCJMVLb6Dkb6DrHbcZkyJXXXrloXZhrA4kf2mhJOUdW8JwplSua6YiPCW1L1OsyKikoz3hYbUg4RKduNfPQf/ho+JHvso1q5G0Y+vqhzPZDVGfls/R6b0RDgSEAEPMlTN+cZkheyLjITs4ffV/IyCXrxLn0BYo7D1v9I45nNhbjcIYYtDRBkiKBzxh/2sWguN+RkOGjTIc+edd+6dOnVqZigUEk8//fSOTp061WsN5PV6tbPOOqtszpw5O+fOnRt7+eWXd96+ffuaqvcDgQAXXHBB1kMPPbRn8uTJJTt27LAMHjw4u0ePHt7o6Gh9+vTp7X766ad1WVlZgQcffDD1vvvuy7j00kuLAT766KOYG264oeN77723ZfDgwTVsjHRdZ+bMmWmzZ8/enpub69uxY4clKyur97Rp0w7+/PPPG4UQA77++utNGRkZwW3bttmqPAAb8yAEiI+PD27dunVt7ed86qmnUjp27OhbtmzZ5rKyMu3SSy/NLCwsNM2ePXtHv379em7YsGFdaWmp9sorr6QsWLBgc3p6emjRokXR5513Xrc77rgjH+DHH3+M+fXXX9d26dIlMHny5PaPPPJI+rx583ZcfvnlnadOnXrwlltuKViwYEH0mDFjesDhPQ+rezBeffXV7e12u75t27Y1eXl55v79++eceuqpjRsK1EMkYmgFfhBCjAf6YYwKn2hqQ8eDYpefYMiYCoZ0iUmDEpePlPimmQ0YU+826OV7w8e0mIaNWcFYSwwFXLBtgZFo3ekspLfUECsh0BI6IvQSpN+Y3oqoRLDGoGmHBuuhwq14Ft8PMkT0uTPRao1EhT0BYbahtx9hGFVIiXTnIxIyoWxfzc3sO4wwNqra/CH0v9S4f9E2vN89iZbQCfvg68OpQ4GtiwGJLfdiox1bbIvkFxY/3mXAEd/kMCTctjWifNjGfAvvvvvug++//35Ct27dPBMmTGhwb9eYmJjQlClTigHGjx9fdsMNN4iVK1eGRyWrVq2y+3w+bfLkySVgeB+OGzeu+JNPPomLi4sLjRgxoiwrKysAcN999x0EY80wLy/PMmHChKwrrrgiv7YQAmiaxoIFC7bMmTMn7pVXXkncuHGjQ0pJeXm5lpHRsCFHYx6EAMOGDatXQM4999zS888/v+tpp51mPeOMM8off/zxPUlJSaGCgoLwek9cXJz+2WefbX733XfjN2/ebFu9enVU9VFer1693F26dAkA9O/f3/3BBx8k5OXlmTZt2uS4/vrrCwFGjx5dkZWV5YHDex5W92BcsmRJzGOPPbZb0zTatGkTHDNmTLM26onE3HWSEOJPwCqgBDhLSlnnr8dvkZIKPz9sKiLGYWb9Hhe6Lrn89A5NFkOoCpToSE8xIiq5wal1FcIWgzDZkbEd4NcXIXMk0l1klPBhTJdNcYY3oJQ6eskuZOkepD0WYYurLBv0YG4/GEvW6DqJ2Jo9Ea2aE43eYzz4y5FCQHkepvgO6OUH0D2VtdhCg1Nvg6UP4y4xRrl64RZs/Sdj7jwyLOx60Edg22Jsp1xzqK3oGrOmZhOpUB0LqnwL/X6/aceOHZacnJwaCZgdOnTwdenSpVHvM5PJVMfDsLrnn67XNYXRdZ1AICDMZnMNX0OXyyU2b95sA8Oqat68eZsnTJiQ9cc//rHojDPOcFe/R1lZmda7d++csWPHlowYMaL8uuuuK1i4cGF8JL6GN998c17VaM3j8Yj8/PzwD1Ztq60qTjvtNPeOHTtWf/jhh7GLFi2KGTJkSPacOXO2dOjQITxi3rp1q+V3v/td9hVXXJE/YsQI18SJE4uHDh0ajrjZ7fbwvat8Equev3q/q3wily5dGnXRRRdlXXvttQfGjBlTdvrpp5ffcsstHavOq+3BWN89mkokqTW5wM3AZ0A+cIcQ4oiiPccCfzDE7vwKnvl0OzaLiaE9Evl4eR4lFX78wdq1x5GhRSVjSuoacYqJiE6BpO7GNPXgasNkth47ML18f+UoTkd6S/CvfZeK96ZAyI+129g6QihssTWEUDgSjW1Jv/i7kc7jK0V35aPFpNVc63MkwmkPYOkyEkun04m+4L9YupxZY4QbylsFSMwZxiZowhpzQth2VacpvoWNUVJSYn7zzTfjwNhjxGaz6bm5uWF7+969e/ssFov83//+Fw+wY8cOy6effppw9tlnl40ePbp82bJlMTt37rQAPPHEEym33XZbO4Dk5OTgqFGjKh588ME9V111Vefy8vIav6dr1qyxuVwu08yZM/deeumlpQsWLIjx+/0iGAw26mvYmAdhY1x//fVtp02blnHFFVeUvPzyy7u7du3qWb9+vd1sNktd14Wu63z77bfRiYmJgRkzZuwfP3582fvvvx8HRjS+IdLS0kI5OTnumTNnJoMhgJs2bXJA0zwPR44cWfbSSy8lh0Ih8vPzTZ9//nn84Z6pPiKZJn8D/E1K+ZIQwgr8C1gNZDanwWNFSbmf937Yz+m9kumSbqw/9+wQw+JV+XRJjyEp9ug7iWv2WGSFHdntPKMaBMMNxlQt+q17y5Aew+E9VLAZ7/IXkb4yY9paR4SEYU/mTKuxVimEQEvsip7UzTCP7TwK3X0QzFa0mHT0QMWhyhNrDJaOv6u3v1JK8BRjzT6vMjfT2A/mRKMhP8P6fAsbIzExMTh37tyE6dOnt3U4HPqcOXO2VrfPt9ls8t13391y8803d3j44YfbhEIh8be//W3/73//+3KA6dOn7xk9enRXgNTU1MCrr766Y82aNeFp9k033VT43nvvJVx77bXtXn/99bAv4ODBgz1nnnlmSVZWVq/Y2NhgZmamr0uXLt4NGzbYevbsWWM0W90D8KGHHsprzIOwIaZNm3bg0ksv7dS1a9eeVqtVz8nJ8UyZMqXIYrHInJwcd5cuXXouXrx4U3p6eqBz5869HA6H3q9fv4qEhIRg1dYCDfHWW29tmzx5cubLL7+c0rFjR1/nzp290DTPw8cee2zfpEmTOnTp0qVXUlJSoHv37s3aQa9BP8PwCUL0kFJuqHXsHCnlJ81psLlU9zOsosztJzbKWMsqLz+0PqbrkvW7S7jvzXXcfG4XEp3GOdvyKthd6OHcgRlkph0a3Ebq71cf1dutD91TjF7lIiOl4awT2x7NHoszOopQ4WZDqPQA0l1EqHib4dMYHg1qCGuUsW5njwt7MdZuV4YChDZ9AsufgpH/ZzjZaBZMiV0g6CVUsuPQ8zbgoyiFFSEPzRa1qJR695yuB+VnqGgVNOZnGMnIcJMQ4m9ALobr9Y3Aoy3Xveazt9BDtN2MSav5x6Kw3EdI13lgYnaN453To+mUFkVphR9dl2ja0bejEvZ4REU+ctMHEPRBjwvQy/ZC0IvuC4EewL92HrJ8P/YhN4bFR9jjjRI7iyNcQdNoOyYLIr0fsvfkQ47XegBZkW/sDmhPQHobWVe2xeL/+WWENbpyam5FOFtmrVChaA1Ekoz5GNAbGIQxAhgD/PtodipSgiGd3fnuGounbl+QPfkV/O2VNXh8h9YG0+IctEmM4r8Ld7Jw5UHKPQ1vpNSSCCEQjiQjH3Dnl4axAnqlcaoPvWwvgfXvY+09MXyNFpNRmUweHZEQhttyxBt+iSU7YKexqaHuKUbqQaN+uva+LFXX2eNBCgIbPsKUVrlWGJPepLYVitZOJD/tZwJXAl4pZRmGu/Woo9mpplDq9rPzYAUuT4DCch87D7r4cm0+CU4rDpsx1Yx32khNcJAUa6dPZhw/by2h7BiJIYBwxEF8Z2Mvlv01p/q+X1/Hkv0HtGhjO2BhT0CLSq7vNodvxxZrlBraYmD924b7NbqxB4pmrll5U4kWnYrmSMS34hVMKT0MyzFzFNoJXnqnUNQmEjEMyGp+91JKHxD5NnPHgFK3n+0HytlXWIEvEOK97/dz/mBj4d8kBBkJhwIRYweks2ZXGcVlvjrbOR4thGY2hKrHBXU2grJ0G4u1x++NF5rliAIWQghj35boNBh4o7H5lLfYSNPRQwir0xhxmmyGm09CJ4QjEb10D8JkwdZvstGNyNYJFYoTikjWDNcIIW4ATEKI7sCtwMqj263mU+ENMqhrPD3bG0GRlDg7ZtMhAWqbFM05A9IpdvtxeYLERDW4NW6LojniCaX2NoIoZXvAmY5/43wsXc8OB0u0mIzDGrseDhGVCN4iSM6G7PHGPs32BKS7EOFMNUoRK/dlkXqIUNE2QgfXYO1xrnG9NaZOvbVCcTIQycjwL0B/IA1YBkRj5B0eEUKIcyr3Yt4ohHhHCFGP5UrTCOlGIueVIzsihMBs0kiKrVumeP04I8JcUO6t5y5HB2GNNgxk3fmw9CH46l6C+34JjxSF2d4iU1Nhth3KLcwcadiJVeQbo8Nqxq5SSvSSXfjXzcP74/PhfV80pxoVKk5OIhHD7lLKq6WUaVLKZCnl5Rji2GyEECnALGC8lLI7sA34v6bcY8XWIg6U1EwnWrQqn6c+OWRGkBrvqDdiXFIR4O7X1uHyBHD7jt2MX9gSDNuwIbdD7mU4TrsznPAc+Yb1h0dzph6yKdv3I/zyHOhBo8ol5DeMGIq3E8xfi3/VmziG3YbQTEZ5n6VZNe4KRaunMT/DfkKI/sDsqu8rvwYDLx5hu6OBn6SUmytfPwdcJhor9q3FW99s57lPNrI9z0V+iZc9hR7eXLKHi4a0xR+UBIOSaFv9U87ubWPZXeihoMzH/iL3sVs7dMQBAhK6QEqv8JRYmGxo9iMeGB9qRzMbWxcIE7QfZlSnHPgVGfIRKtxMqGAjMlCBf808bP2vRItta7hoq7VCxUlMYyPD64C5QGdgXuX3c4HXgHeOsN32wO5qr/cAsUCN7GchxFQhxHIhxPL8/PwaN7jnkt6s3VXKFyv3U+TyMfOjLZzaLR5Nk+zKd+EJhNi8r4wyd909P2wWE4O7JbBsQxFuX5CdByvQ9aMviMJkRVjqrsc1ZAV2RG1ZHMbmT5oZci6BDXNr7NUsQwHsQ2/C3HkkYJjeViV0n6woc9fIuO+++9LGjx+f2dzre/TokVPd5KGpTJw4seOSJUta3FyzMT/DqVLKTsBjUspO1b66SilvP0rt1igallK+IKUcKKUcmJJS8/8zJsrCzef1YOPuUrwBnUt+14aRuUZKikXTcNqNz3pPfkXYuaY65w9qSzBkiEOp28+GPaUUuxqty28RajtjI8yIo2SnL6xOIzqdPgD6Xxu2E/Nv+gzP4ulGlFsII50nwvXKYzWKPh4oc9djw4YNG9YlJyc3zyAAWLJkSezR+DmM5EO6VwhxmhDigur2/0fY7i6geg5JW6BYSlnRlJukJziYOrYbDquJBKc17Iaxbq+Ld7/dx+4CNyEpKS6vK3JDeiQxfkib8C93IKSzK7+CTXvrH022FIabzSE7LK3S4zASfIEQecUeCsq8EYuS5khAWGMgtp3hvP3p9fjXf4j91BsqT4g8nafCa+RytnaUuWvTzF19Pp+47LLLOnTo0KFXv379enz33Xfh4vrCwkLT+PHjM3v27JndrVu3nKuvvrp9IBDgX//6V/IZZ5yRVXXeihUr7Kmpqb2DwSBCiAH79+83l5WVaRdccEFmnz59emRmZvbq2bNn9sqVK21V/bnhhhvaDhw4sHvbtm1zzz///E6hUIibbrqp7cGDBy2TJk3qvHjx4hZNe4hkXvQKRuL1Fg65m0qMqXNzWQA8LoToWrlueC3wwWGuOSy+gM6rX+0mEJL0bB/LPW+sZ8aknljNpjq2XU6HmSXrCvl1eynTLukTPu7xB9l+wIXdYiIxxkZctAWruWVNHbRYw1ZfmO1GdUoEVHiDbM8rJ1QpgiG/hy7pMRGVFGqxGYQKt0KXsyFjINGJaWETCM2ZHlE6TzCks+NABSnxzQ+wuL96uI3v55fDyuu85K31AK45h+ombQP+tD/q9Lv3lTx7Sm/pKbIAaElZ7rgrP1/v+vgvHQMbPw5npMf++etVoX2/RFXMvyWr+rWH64cyd22aueujjz6asm3bNvumTZvW+v1+MWTIkO6xsbEegGuvvbZ9//793XPnzt0RDAa56KKLOk2fPj3ttttuy58+fXq7Xbt2mTt06BB84YUXkiZOnFhQ/Y/K3Llz4+Li4kIrV67cAHDppZd2eOKJJ1L/97//7QbYvn277YcffthYVlamde/evdf8+fNjnnrqqb3z5s1LnD179rYRI0a4a/f1SIhEDEcAPaSUTXaObQgp5UEhxFXAu5VOOFuBSUdyz3JPkFe+3E2HJAd/PS8LkyZIibPx+AdbmDGpJ+XuQI2cQpOmMbxnMq9+tYtftxXRt3PN6as3EGJfkZt9RUbittmsYTFp2K0mYhxmYhyWRs1dG0NYHJiSstDiEiO6h5SSPQUVYSEE8PqD7Mp31TCcaLA9k/XQJvbRKWEhFBYnWoQO1vuLPM22Pqsi6vS799UnVvX5HMZf/9Oq2sec5/5nJ+f+Z2f1Y6a4dqXW7POa7JOozF0jN3f98ssvYydMmFBkt9ul3W6XF198cdHq1asdAIsWLYr79ddfo1999dVkMP5AACQkJOjjxo0rfvHFF5Puu+++A3Pnzk365ptvahi+XHXVVcVdu3b1Pfzww6lbtmyxfffddzEDBgwIzw7POeecEpPJREJCgt6xY0dfYWHhUbWaimR+tqslhbAKKeV8KWUfKWW2lPJcKWVRU65/9pMNrN5RzI4DLr5ec4D/fLyNbm2iuX5cZ0yVo6VzB6ah65IV20oprWfqmxZrZ8roTGa8s5b9RQ27/oSkxBcI4fIGKCjzsv2AizU7S9hxwEV+qbdZ6TnCZI1YTAvLfXgDdYWo3BOgKMJpq4hONvIcDx2JeHrs8QcpLPfy95d/Zs2OZpkI/+aoMnctKyszVU0Vq3M8zF1XrFhhB8OD8MMPP9z0zjvvJH355Zd1AgVV5q6//PJL9IABA9wzZ87cbTabZaTmrhs2bFi3YcOGdStXrly/bNmysEA1ZO5a9WxVVDdP1XVdzJkzZ2vVPZcvX77+pZde2gUwderUgrfeeivp3XffjcvKyvL06NGjxi/hjBkzUqZMmZIZFRWlX3HFFUUXXnhhUfV2oqKiwi8ql8COqrNKJGK4TAjxlhDishZcMzwidF1S5g6weGUe7yzdybb95Vw6vC0XDGqDw3roj4cQgnMGpDP/57x6jRlioi2ckpXAXZf0IjXOxoFiT8RrcbqUlLr95BW72bq/jM37ynAdpXrnwrKGfycPFHsI1fNLVxshNLTYdlTVJmsx6RHnFB4o9rJudykVvhCpRzBN/q2gzF2bZu46evTo0jfeeCPJ7XYLt9st5s2bF55GDR8+vPSxxx5LqxxlijFjxmQ9+uijKQBnnnlmBcA///nPjD/96U91rNwWLlwYe+mllxbccsstBb179/Z+/vnn8Q0ZuFbHZDJJv9/f4sIYyTR5SOW/f6527EjXDI8ITRNMuziX7PZx+IOGEHjcLhJijMCESQjapUTj9gYZlpPE/77axe4CNx1Soom2HxoEOKxmLCaN3Ewjmvv0xxuJcVi45fxsLOamBeC8/iDbD5STHGMnLaH+ZO/mUOEN1jsqrCKo6xSVRbavi7BEYUrsginKdthtC6pw+4KUuv1kt4vjL+dlHxPbs6ONMndtmrnr3/72t/wtW7bYevTo0TM+Pj7YqVOncBsvvPDC7muuuaZ99+7dewYCATF8+PCyBx988EDV+5MmTSp4/PHHMy6//PKS2ve9/fbbD1x//fUdX3/99WSTySR79+7tXr9+/WF/kMeMGVNyxRVXdH7mmWd2XnjhhQ0uYzSVw5q7hk8Uwlx5/rGze6lGfeau3204yJ4CYw1VC3nJSDQ+x8zUmPD64K58F4+/v4nYKAvXjelCm6Sas479hW480jjX6w/xyNurAbjvj70xmQ4viBWuuuauZpNGSqydpFhbo1PhSExld+dXUFQt5ScY0lm2Lp/+mQ4EhMsOu7eNi1iommJmu+OAi89+2cvOAxX06RRP786J5LSPr32aMndVtAoaM3eNZA+UVCHEp0AF4BVCLBZCNHm3+qOB024hPcFOfLSVtMrpW4zDUiNQ0iYxiuHZSXy7obDedcPEWBui8nfZbjVxz8TeOGxm8kqaVrdcUhFgb6GHQEgnGNLZX+xmx0FXRFPYhpBSUlpxqM8rthbx5ye/45Of9lDuCfL85ztYtqGQYEinpKLl04G8/hAlFT5eW7wNs0nw7rJdfPRDkwZPCkWrIZK54NPA9xhGDanAEozyud8EsVFWUuPt4VFRWkLNUbbZpDGkezLlniA7DlZQ4a05sLVZTMQ7D+X9Wcwa0y7uRXqCg2/WHDjsGmKFN8hzn23nxhdW8q/3t+Dxh8LXuDwB9hQ2P/pf5g6EI8jbD7h49N21XDu2G/935QD0kGRQVjz/7/MdlHkC9Qr9kXKw1MvanaWUe4KkxtlYsbWIcwe3a/F2FIrfApGIYTcp5XQpZYmUslBKeT+QddirjgOJMTYc1rrLoCkJdoZ0T+LbDUX1jqDSE2rafIExKprzzQ5mvLu23sCIlJJQSPLDpmICQZ2np/Th31fnEmM388jcTWzYa0yfyyr8zR61VV1X5g7QKc3Ji38Zwildk9mV76LQ5SMu2kxWejRL1xZS4QnWW2nTXPzBECUuP9sPuDizTzqf/LSXsQPbklKPC5BCcSIQiRhahBDh3wAhRBSHkq9/M5g1jbQGgghWs4kze6fww6YiyioCdeqQrWYTndKcmKvtpRJtN/P4nwdis5iY/MS3HCzxsqeggpc+38zj89Zy9X++492lO+mQYmfcgFRKKvx4/TpCCM7ul8aMeZspdhlidrCk6Zt1hXSdMneAvYVupj75HWXuAA6rmb2F7nDQCGBojwQsZg1JzSn1kVJQ6kMiOS03jU7pTrbsL+fCoe2behtd1/XWH3FRnBBU/iw2OGKIRAzfAr4QQlwthLgao3rk3RbqX4tgEoIOqc46o7vqDM1OpqDcz74iT701yFE2M1ltYmpUm9itJm45P5vnbhgcLveLspvplB7D5ad3pk/nQzXFQV0nr9hDKCQZ0CWeUX1SePbT7cjKHMWmlviVVATQpWTWwi1cMLQDsVEW8ku94cjyviIvP20poUOKg44pdsrcgRbbyiCk6xSV+3jty218unwvi37N47TcNBJjGt31sT7W5OfnxylBVBxvdF0X+fn5ccCahs45bGqNlPIfQojdwFgM8ZwFvNxivTwCUuLsmDRBXLQFn6fxtblEp41BWQn8sLmI9ilR9Zq+2izGCHHLvrIa1R5VuXVtk6IYN7BdDR9FbyDEqh1ltE920CEpCokxSr34d21Zs7MsHE0uKPOFtzWNhKJyH7sOVrBqewm3XdgTlycYXhdc8Ms+Fq/cT4LTwrrd5SQ4LaTG2bnitA6EdL3OboFNpbDcjycQ4sPv93DlWZ35dVsR/54ysMkVN8Fg8M95eXkv5uXl9SKyP7wKxdFCB9YEg8E/N3TCYcVQCBEDZEgpLxFCdARuAaIwosvHlfRqwZLD1WFomuCM3im8+tUufn9KBoVl3noF0W410TY5il35dR+vxOXnYOmhKHOxy8cz83cQ4zDz2Yp8pozKpGNaNKlxdvYXe+jTKY7Fq/I5rVcyFd4A/mAoojpnjz+I2xckKdbG/Zf1xmbW2FFk9GfF1iKWrD3IXeO70SUjmjtfXUdqnJU1O8sIhoypdYKzySO4MFJKCkq9LFt7kLZJDjbvK6d/ViJtk42UpLgmbJMwYMCAg8B5ze6MQnEMieSv9StAlZtACcZ64X+PUn+OKiNz09hb6OVgiY8DxQ3X2iY4bUTbav6dqC2EgZDOfz/bQu/MWKaM6sBt53Vh1qKdlLkDFJT7SIqxIYBFq/P5Zq2RaldSHtlU+UCxlzJ3gG155fTsEE9BmY9gSOdgiZc53+zgmrHd6Fo5pT9vUAa78r3sOOim0OU/4i1Qi11+AiEdjz/E6blpLFl7kNH92mCzmIiLtmKzHNXyUIXiuBGJGHaVUv4NQEpZKqW8BTgi88njRWKMld/1SOKbdQWEpGR7ngtfA9Ud1ZOzi2sJIcC8ZbtIcFo5q3cy0XYLfTrFc2r3ROZ9vy8sshaTiQtPbcPHy/OQUlIcQYCjquJj4Yp9fPLTXvwBndKKAP6gzksLtnDOKW05pVsSZpMxZR3aI5H8Mh8DusRxoMiDy103QNQUDpZ6CYZ0RvZJp8IXIj3eTo/2htehiiQrTmQijSaHPemFEE4aqDj4rWPSNM7ul8ZXawuQUuIPhtiyr4z9Re46ZgtRNjNxURbyir3k1xLCHzcWsHFPGZef0QlNEyRXisT5gzNYtKoAjy9EYbmPaIeZfp3j8PhD7Djoxh8MNVq/bLjTGNsQfPLTXs45pS35pV50qfPW19vJSHTw+8HtiHEYU1W71YzDYmZ4TjLRdhN2m0a5N9jsQEpRuQ9fIMSMd9ay6Nf9LFl7kOG5acQ6zDisZqLtJ7cTtuLEJhIxnA38IIR4UAgxHfgOY+rcKhncLRGzJvh5awlgGC4cLPWyeV8Zq3cUs3lfGTsOuNiWV06xy18nSXvzvjLmfbuLP5+dhcNmJtFpI9pmom1SNJ3SnPTOjOWLVQeRSLz+EGZN47HJveiUZvhQFjXipn2gxIvHH+TXbcWYNEGnNCcVviDfrDHKDqeO6UpKnCG8NouJzulOurSJ4Xc9kli1s4xnPt1Bmbt5Cdi6Ltlf7GHHARcrtxdjt5rIL/UypHsyJpNGUmzz1yEVitbAYcVQSvkIcAcQBziBv0spnzjaHTtaxEVbufh3bZmzbG+d6hJdyvA0tdwTQKv0RKxie56Llxds4apRRo2z024mwWmlTVIUiTE22iRGcd4pGXzy8wFCupFSI4SRszj3u724vEHK3YF6k6NdnkA4St2zYzz3X9qHgjIfW/aV8dnPe/nrH7LpmBodjui2S47GpGlYzSZO6ZqIyxPkQImP/UUeSsp9TS4DzCv2EAzpvP7lds4b3I7vNuQzNDuFxBgbJiFIiI48Eq5QtEYiSneQUn4opfyLlPI2KeWnR7tTRxObxcTpvVIIBHW+WJV/2PNjo6wkOK2s2VnC//t0E5ef0Znu7eKIj7aSkejA6bCEU2bsVhODuiUS4zCHR56BgI6uS7YdcLN0XaExEq1V9+wPhthZGb3OL/Xy06YCoqwmCsq8zPpiK1ePzqJ3p4SwECY6bURVC/Akx9ro1zmeRKeFjftcFLv8FDVi+1UblydAfpnRp9xO8ZzSLYnlmwsZ0SuVaLuZeKf1hHCrUSga46TM/YqPsnL7+V15/evdfLk6n1AjAYcKb5C5y3bx9pId3HJ+NsN7pdIxNZrUeDtCCFLjala9JDptnDMgnY+X5wGgmQTegM7I3BQWVYpvUbkPr98IsoR0ne15rvBo8c2vt7NmZwkFZT5eWbiV32WncFa/NmExMmta2J2nCqfDwtDuSUiMDfDK3AEOlkbmzRjSjb1f/EGd5ZsLOfeU9ixbl0/3trF0SHEihFBTZMVJwUm5Ih4TZaVtkoO7L+7OS1/sZNbiXSQ4bZg0gdNhJinGhtNhodjl49dtxQzLSeW5G04lwVlzqui0W4iuFWB1Oiyc1iuZ17/Zzeb9LrpmGHvn9Gofi8sb5ECJl7R4O9v2l+GTFkqrTZv3F3n4evVB/nFFHz79eS8Al53RGWs1b8X0REe9lTan56bwn0+2ktsxhpCUFJT5KSr31ZtLWZ09BUZg55mPN1Lu9pOZ7uTLVXlMPC2TuGgL0TZzvfXeCsWJxkn5Ux5tN0xdu2Y4eeTyHALCRmlFgFClg3ZhuY8Kb5BeHeO5+bweDSYxp8XbkcG6Vl/JMXYuGJzBu9/u487x3XBYNYpdAZ78c++waWxIyjo7zc3/aQ9jB7ZhW56LJWsPcueEXiRVK4GLj7bWeF2jzVg7WenRvLFkL6dkxdMnM459hW7ioq0NlimWVppIvLN0Jxv3lPKvqwfy2c97cVhN9M6Mx2YxhQM2CsWJToNiKIT4iEYMGaSUrbqyICnGRl6JByEECU5bk6s2Yh0WnA4L5eV1xTAu2sqZvVN597t9bD9QEY4k6zq89MVOJvyuLTGOmh+9PxBiwvBM1u4q4fH31nHFyM50SnOGp8cxDgsdUhreGdHpMDMwK54vV+ezYnspuR1jySvxEh/joX1y3etCus6eQjeBoM66XaU8eHlfPL4Any433GkSnDasZiPRWqE4GWhszfBdYG4jX62aeOchU9emoglBRmKdfXrCRNvNOO1mxp/ahte/McxQnQ4LRRVGdcecpTUNUr/fkM9fXljO3kI3r3yxlWE5qfTqGB8WorhoK53SnI3WBlvNJoZnJ3Gw1MfmfRWUe4K4vAH2FbrrtRDbXeBm2/5yPP4gD1zWhyibmc9+2Y8QgoHdknA6LDUi6QrFiU6DYiil/F99Xxh5h98euy4eHSxmLbxnSlNpmxSF3dp4WVpclJXRfVPZle9m/Z5yrGaBrutcMLgNS9YXsm63sXXDV6vzePy9dVwyoiNzv91FlM3M2QPaEBtlMbY7jbWTmdq4EFbROSOGzNRourVxUhS2D/Oy+6CrRrJ3fqmXA8Ue7nvtV1btKMHtDbJ1fzkffL+bS0Z0JC3egdVsanBKrlCciERi+3+NEKJMCBESQoSAIIbbdasnIyGqybW26fGOiKysYqOtWMwaE4e149WvdiOlJDbKis8f5MaxXdhd4KHcE2DWgq387cKe7C/ysHZnCZPONLY6jY+2kh7vqLNnS2NE2y0M7ZGINxCifZIDjz9EUNc5WOpjxwEXu/Jd7DzoYl+Rmze/3k7XtrEM6prM3gI3s77YwsCuSZzSLRmH1URGoqPZ+0IrFK2RSFJrpgGjgE+AfsB9wHtHs1PHCk0TdMmIoW1SFBkJDjJTnWRlxNI5LYbUODv2akJps5jITHXW2VagIaLtZmwWE6f1SsblDbJ8awnRNhMhXdIm0cZZvVOJtpl5+Mp+gOSdJTu5enQWUTZjg/r2ydERt1WF025mSPdE9hR6+WptIbO/3G1sq+rxG1sDVPgpc/up8AZZuGI/U87OYl+Rmw9/3IM/oDN1TBbx0VbioqzEq7VCxUlGJNHkIinlD0KIX4E0KeXDQoifj3K/jhkmTSM+pm7ENCbKQkaiUaYmkc3yCExw2vAFQlxxentmf7mLfp3iiI22UOzys/2gC4fLqEd+/avtnN47nfYp0QgB3dvENlkIwRDspBgbQ3skkF/qQ5ewcGU+Z/dLJb/Ui9kUhcOqEW038+LNQyh2+fl+Qz4/bS7gqWsGkRhjbH/QLjny0ahCcaIQyW94QAiRAGwGBlUecx69Lv220DTRbLPUBKcVgWBgl3iSYqx8vuIgsQ4r1Wef363Pp8IbZFS/DADaJ0eTmR75Vp61iXdaGdU7hXV7XJw7IJUV20vZfsCNLiX7i9y8tWQPry3eRpHLz678Ct5espPpl/UhJd6OJgSZh3EMVyhOVCL5qX8B+BhjmnyNEGI5sP6o9uoEwWzSiKvcLuBPZ3bk7W/3UujyEeswpqD5pV4+/GEPV5xhrBM67WYGdEk6zF0bJ9pmJj3BwfDsRBauKuCmcZ3ITHWwfk85r329m0Wr8umcEYPLG+C1L7cxYXhHurczLLraJkUpZxrFSUskRg0vA6OllEXAEOAfwB+PdsdOFKo8ADukRDFhaFsee28zQkCFJ8jz8zdxzqB2tEmKwmLSGNIjFbP5yEZlToeFaLuZkbnJFJUH2LSvAiEEXr9O+2QHN4zNJDnWxpK1BwG4eHgmYCRtN2OPE4XihCES238NuFwIMQoIAJ9KKZu+3dtJit1qIjbKSpnbz7gBaZS4A9zy8mp0CWf1zWB4TyOQ0rNDXI39m5uLzWIiymYmJc7OxOFteHHhLtLibfTrHBc+J7/Uy6c/7eXfUwdi0gR2i4k2iU1fo1QoTiQimRM9DvQBXsMwdZ0ihOgqpbz3qPbsBCI13k6Z248QgstGtOecAelYTIKUpHh0aQhY+5SWW4aNi7Li9QfpnBrNHwal8/KiXUw+oz3tkx2Vm95v44+nd6JjqhOBoEM1azCF4mQlEjEcBfSVUgYBhBCvAz8DSgwjxGE1E++0UVJp7BofbThVCyEwCWNjq5a0yIqJsnCwVJCeGMUpWRKLWePFL3aRFmclr8THmX0zuHhYR8Cw/1JGDApFZGJYXnlelS++zm9gZ7zWRnq8g/IKf40tSMEYFSa2wPS4OlE2c2UyeYg2SdFomqBrRjTbD1TQNtHOoJ7GZvBmk0ZagjJiUCigcaOGWyu/zQO+EUK8CoSAiahocpOxmDXSEqPYV1jz70j75KijMkVNirGxr8iN1SxolxRFuTdIWrydWIclbAmWFm8/4j2WFYoThcZGhrmV/5ZUfvWvfL2VRtxsFA2TFGPD5w9RWOl0kxbvINoe+T7ETSHBaeNgiZegrqNpos5+x6r2WKGoSYNiKKW8qvrryg3kLVLKLS3RsDCGQ7OANVLKf7XEPVsDbZKiiHFYkEjSm1FlEinG/i129he7632/bZKqPVYoqhOJUUOWEGIt8CvwsxBiqxCix5E0KoTIBhYBE47kPq2VmKhD+6YcTZLj7Dhsdf/exTqOTfsKRWsikgWjp4FHpZQJUso44CHg2SNs9waMUeHbR3gfxWHomOKs4cxjs5jokNqwSaxCcbISiRimVfoYAiClnAWkHO4iIcQ4IUSwnq9JUsobpZSvRnCPqUKI5UKI5R6PyvNuDhazRpeMGDISokhPiKJzeowKmigU9RBJao1ZCJFYWY6HECKZCAIoUsr5Ed6/sXu8gFEbTSRtKurHpGkkq71MFIpGiUSsngK+F0LMqXx9CfDvo9clhUKhOPYcVgyllC8IITYDYzCm1ddLKb846j1TKBSKY0hjSdf9q70sBeZUf09K+cvR7JhCoVAcS4SU9S/FCSG2N3KdlFJ2PjpdarjNY9yeInJUwqKi1dOgGP4GaTUdPQlRYqho9TS6ZiiEaAdESyk3CiEeAmIxDBvuklLW3T1doVAoWikNJpwJIQYBvwADKw9dAhQCvTCSphUKheKEobE1w8+B/5NSfln5+hcpZX8hRAfgXSnloHovPHqoafJvFzVNVrR6Gpsmd64SwkoEgJRylxDieHjEN/oLJ4RYAxztqXsyUHACtGGXUvY6ym0oFK2KxsTQX+v18Ebe+y3glVIOPPxpzUcIsfxEaeNo3l+haI00VqRaXhlAAUBK6QIQQrRHOV0rFIoTjMbE8L/AG0KI1KoDlZvJzwKeO9odUygUimNJY+auLwkhugDbhRDrMAIYPYD/SCnfPFYdbAIvHP4U1cYxbEOhaFUcNum6cmQ4BCOA8YOUcv+x6JhCoVAcS1pTBYpCoVAcNZTLp0KhUKDEUKFQKAAlhgqFQgEoMVQoFApAiaFCoVAASgwVCoUCUGKoUCgUgBJDhUKhAJQYKhQKBXCEm7wfY1SpzG8XZe6qaPWokaFCoVCgxFChUCgAJYYKhUIBKDEM4/f7KSsrQ7n4KBQnJ0oMKyktLcXlclFeXn68u9JiuFwuDh48iK7rx7srCsVvHiWGgNvtJhAIAIaAVH3fmqka6QaDQUpKSo53dxSK3zwnvRiGQiFKS0trHDsRRofVBdDr9RIKhY5fZxSKVsBJL4b1rRN6vV6CweBx6tGR4/P56vTf7XYfp94oFK2Dk1oMQ6EQHo+n3vdcLleLtuXz+SgtLaW8vLzRNTyv13vEQZz6hE+JoULROCe1GFYXiF27dtV4z+Vy4ff7W6Sd4uJiCgsLqaiooLy8nIKCgnqnrcXFxRQVFVFUVNTstqSUeL3eOsdDoVC9xxUKhcFJK4ZSSioqKgBYsWIFQ4YM4YYbbsDj8eD3+ykoKGDXrl0UFRXh8/maveZWXl5eZ/QZDAYpKiqqMQL0eDzh83w+X7ODOB6Pp8GRZUuPdhWKE4mTVgw9Hk94uvrEE09w7733UlpayiuvvEJJSQm6ruP1evF6vRQWFnLgwAHKysoivn8oFKKgoKDBYEwgEKC0tBQpJcFgsE4Qp6HpeyTP1RB+v7/FRrsKxYnGSSuGVaPCVatWsW7dOq666iqmTZvGCy+8EJ5O+v3+GoEIl8uFz+eL6P7l5eWHFR63201BQQEFBQV11hGbM6XVdf2w/at6boVCUZOTUgyDwWB4GvrMM89w7bXXYrPZ6NKlC506deKzzz4Ln1t7pFU1mgManI42FpipTSAQqDegEgwGmxzRrt3m/v37efrpp2v00+v1qiRshaIeTkoxrBp1bd++nW+//ZZLL70UKSXl5eVMnDiRuXPnhs+tvQYXDAYpKCjA4/GQl5dHUVFRHXFxu90tUtbX1Clt9VHhzz//TExMDJ999hmPPPJI+LiUssHIsq7rFBQU1JmyV0flKypOVE5qMZw9ezaXXnop0dHRuFwudF1n6NChHDhwgM2bNwOGQNQecbndbrZs2UIgEAivKVYXP5/Ph67rRyyITRXDqvNXrFjB1KlT0TSN2bNn89prr5GXlxc+r7F0Ir/fT0VFRYPTbbXmqDhROenEUNd1/H4/Ukrmz5/PBRdcQDAYDI+WzGYzQ4YMYdasWeFrqoQSDCEtKioiEAhQXFxMKBQiEAiEr6+6f0lJCYWFhRGNpAKBAEVFRXVEqinCUzXdllLy4IMPcvvttxMVFUViYiJTpkxh9+7dNc6t3a/q0XWgXjEMhUJqZKg4YTnpxLDql3zNmjWYzWa6desWXgesyu8rKChgwYIFvPXWW4AhcEVFRZSUlFBSUhIe8YVCIYqLi8NTbCklPp8vPMKqL0pcG13XKS0txe/3h+9RRfW1zaY8V15eHhdffHH4HrfccgsDBw6sIWS1AzQVFRU12q5PiNWoUHEic9KJYdUv9Pz58xk3blzYpGHx4sVcc801hEIhnnrqKdq3b88rr7zCL7/8AhiiUl+ENxgMUlFRga7rHDx4kP3799fI5/P7/bjdbn788UcuvvhiBg4cyNtvvw0YYlpUVBQOlOi6XuPasrIy9uzZE1H1SNVz9ezZk7fffjucK1lQUEBZWRnTp0/n1VdfDZ9ftRbqdrspLy+nsLCwzv1qT/OVGCpOZE5KMZRS8vHHHzNmzBjcbjf79u3jn//8J88++ywOhwOAs88+mxEjRtCvXz9WrVrVqCC5XC48Hg/l5eX1OsS43W5++OEHLr30Up5//vnwvUpKSupEjCsqKsJTbLfbjdvtprS0tNEIcFVKjdvtZs6cOaSlpVFaWhq+t9vtZuTIkTzzzDPhEaTf72fPnj0UFxezZ8+eGqJc/bNq7LVCcSJxUomhrusEAgF+/fVXpJRkZ2ej6zrPPfccV199Nbm5uURFRQFw5plnsmzZMsAYRV5++eVs3bq1wXtX1R1XUVX+Nm3aNFauXMnUqVO54IILGDhwIFdeeSXvv/8+69atq/deXq83PD32+Xzh0WVDVI3yXn31VRYvXlzvud26daNjx468+uqr4ZK/goIC8vPzw9fXntJXXzesWhtVKE5UTioxrBrZvPvuu4wfPz487b377ru5/vrrAbDZbGiaRlZWFkIItm/fzrRp05g8eTI33ngjxcXFDd6/alpZXl7OyJEjOeecczCbzXTt2rXGmlwwGCQ/P5+77rorHDQpKSlh5cqV7Nu3r859XS4XLperweCF2+2moqKC5557jr/+9a8NRoKvueYaEhISwpUoUsoaI85AIFDj2oa+VyhORFrTVqFHTJUIfPjhh3z88ce4XC7++c9/MmPGDGw2GwBCCOx2O263myFDhvDtt9/SuXNn/vCHPxAMBg87VSwuLiYmJoYvvviCPXv20LFjRzRNw+/3U1hYSFRUFF6vl3PPPZfvvvuO+++/n+nTp/PSSy+xatUqdu/ezYUXXsiNN94YvqfP58Pr9VJSUkJSUhJAWMjKysoIBAJ8+umnDB06lM6dOzcYtMnNzSUQCFBYWBi+T20qKirCn0VVhFrTNCWGihOek2pkGAgE+Oabb+jcuTMZGRm89dZbVFRUEBcXV+M8u90OwJAhQ/juu+/Cx8ePH4/D4eDDDz+kvLycffv21QkyzJgxg48++giLxUKnTp3QtEMfcTAYpKysDL/fjxCCu+66i7Zt2yKlZPjw4cycOZN58+aRmppa575Vhg8ej4fS0lLy8vI4cOBAeGR50UUX8Z///OewlS8ffvghM2fObPB9v99fI1BUJYJKDBUnOiedGH7wwQecf/75BAIBPvvsM6688srw+1WjQpvNhslk4pRTTmHVqlU1BCYUCvHSSy8xfvx4LrvsMv70pz+FU3LcbjfffvstZ5xxRkT9cTqd/OlPf2LixIk8/vjjTJ48mY0bN5KUlMR1113Ho48+Gk6WDgQC5Ofns3PnzjppMHv27OG5557DZDLVGblarVaEOLTH+4gRI1i6dGmj638ulyt8f5/PF46WKxQnMieNGIZCISoqKvjiiy8499xzKSoqwmazceqppwKgaRopKSkkJiaSlpZGYmIiMTEx9OzZk6VLl4bvk5CQwKhRo/B4PIwfP55BgwYxZcoUPB4PX331FX379kUIwX/+8x/efvvtw47Unn32WQYPHsycOXP485//zPPPP8/rr7/O6tWr8Xg83HnnneEor5SyXouujz76iO3bt9dJ/bHb7SQmJpKUlBQWxJSUFDIzM/npp58a7FP1JHSPxxN263n88ceZPXt2JB+3QtHqOGnEMBAIsGDBAvr27UtKSgo2m41Zs2aFo8dxcXGYzcYSqqZppKWlYbPZOP/883nvvffC99m+fTvz5s3jhRdeoGvXrlx33XV07dqVWbNm0alTJ6644gruvfde9u3bxw8//MDkyZPZv39/vX1auXIlixcv5qabbgLgvPPOY9asWbz88ssMGDCAYcOGYbfbef3118PX1DZvDYVCvPnmm1xwwQX4fD6klDz55JOcf/75jBs3jn//+9+Yzebw1B9gypQpJCYmsnDhwhoR8qqKnC1btoQ3xqpK8QFYsGABWVlZR/pfoVD8JjlpAih+v5+5c+cyfvx4AoEAjzzyCFdddRXp6enYbLZwfmEVVquVhIQEzjrrLP71r3+xe/duMjIy+Mc//sE111xDdnY22dnZzJ8/n7Fjx3L//ffjdDpZs2YNHo+Hxx9/HLPZzAsvvMDDDz/MU089VWO66vF4eOCBB5g2bRqxsbF1+jtgwAB+/vln/v73vzN16lQuueSSsKD5/f5wf1evXk1ycjKDBw8mPz+f77//nm+++YZnnnkGp9PJ1KlTsdvtTJ06NbylwNChQ3n88cdZsmQJPp+PIUOGkJ2dzWuvvcYdd9xBp06d8Hg8NRKxXS4Xu3fv5pRTTjka/z0KxXHnpBkZ7tu3j59++omxY8fi8/n45JNPSExMBKhXjAAcDgeJiYlcdNFF3H333dx77704nc5wqRtAdHQ0DzzwAGazmc8++4z09HSefvppLBYLQgj+9Kc/sX//fr788svwNVJKHn30UXJzcxtcXxw5ciTDhg2jU6dO5OTkMH/+/PB7Vet9e/fupU+fPrz99tsEAoFwruGVV15J//796d69O2+88QZvvPEGjzzyCKFQiD179vCPf/yDDz/8kL59+zJ37lzKy8t57bXXeOKJJxg6dCgrV64MpxqBIb6LFi1iyJAhLeLGo1D8FjlpxHDevHmcddZZREdHs3nzZuLi4khJScFut2OxWOq9xmazERUVxU033cSIESPweDzMmDGjRoT4tNNO44033mDGjBm8+uqr3HrrreHUFACLxcI999zDI488wj//+U/27NnDiy++yNq1a7nzzjsb7G/79u3p378/Pp+PSZMmMWvWrHCpXjAY5J133mH06NFs2bIFs9mM3+9n8+bNbNu2jQsuuCA8Cm3bti0ffPABO3fu5LTTTuPKK68kNjaW119/nR9//JG//e1vXHrppWiaFt4HJjc3l127drFnz55wfxYvXkz//v1VFYrihEW0or/0ze5oKBRi1KhRTJw4kXPPPZe5c+fy0UcfMWvWLBITE+tMkatfd+DAAXRdp7i4+IgqMEpLS3n11VeZO3cuw4YN47bbbiM+Pr7Ra2699VZGjhzJueeey0MPPYTH4+Ghhx7i008/5dlnn+WNN96gW7duABQWFnLXXXfRsWNH7rjjDkwmU537uVyuGjXWHo+HL7/8ktzcXPLz87n77ruZO3cuUVFRPProo8THxzN16lQCgQBnn302r7/+OgMGDCA9Pb32rUWdxhSKVsZJIYYlJSW0b9+e77//HqfTSWFhITabjYSEBNLT02uM9Gpz4MABQqHQEQti1Rpk9XXDoqKiBkdaZrOZJUuW8MQTT/DWW2/h9/v54x//yJ133smSJUsYP348Q4YMAYyR4vr165kwYQKff/45HTp0aLAfVXuz1Pf/ft9995GcnMzNN9/Mli1b2LhxI+eccw5fffUVr732Gi+++CJt2rRRYqg4ITkppslLly6la9euJCQk4PV6mTdvHmVlZVgslkaFEAwRAyPCnJCQ0OCUujo2m61GdLo+IQSIiYmp8VrTNKxWK4mJiSQnJzNmzBgcDgdffPEFDoeD6667jieffJLbbruNTp06hcvzvF4v7733HmPHjm1UCAFMJlONaXx1brrpJj788ENWr15NVlYW48aNo6Kigk8//ZSxY8ce9rkVitbMSSGGixYtYvjw4YAxNXz22WexWCwNikJ1qsQQDglifVPQ6sTExISDL3FxcSQmJtYRQqBGH6KiosJ5jlVtOhwO7r33XnJycgBo06YNW7ZsCRtIVFRUhPdb+eSTT7j88ssj+DRocFkgJSWFu+++m2nTplFSUsKHH37IHXfcwffff89ZZ50V0b0VitZKi4uhEOIcIcQqIcRGIcQ7QojYWu9PEkL8Wu1ruxAiIIRIa+m+VLFkyRKGDRuGz+dj7969mEwmMjIyIhLD2iNBTdPqjOiqY7Vaa+QrNiQ8VTidTqxWK7GxsfUK5pAhQ0hKSuLhhx/mnnvuYfTo0eHIstvtJj8/n19//RWLxUK/fv0O+zxwyIyiPs444wxGjRrFvffeS2pqKt9//z1XX311nZJFheJEo0XFUAiRAswCxkspuwPbgP+rfo6UcraUsq+Usi9wCpAH3CilPNCSfanC4/Gwbt06+vTpg9vtZuXKlQwYMAAhRERT3qoUmerY7XacTicWi6XGyBEIJ3FHisVioW3btg2+bzabiY2NZdiwYfz1r3/lr3/9K8uWLath0/XJJ58wfvz4Gv0UQtQrrlU09ofghhtuwOPx8OCDDzJ27FjWrVuHlJKlS5eq1BrFCUtLJ12PBn6SUm6ufP0csFIIcYOs/7foDuCglPL/NbfBKqv96hUW1Vm5ciVt27bFarXicrkYPXo0Y8aMQdO0w64XAmHRrB3oiImJISMjg/LycioqKigvL0fTtBoiYzKZsFqt+P3+Bu23qtYTY2Nj8fl8lJeX1zk3Li6O0047Lfy6T58+fPXVV4wbN46ysjIWLFjAF198EX7fYrGQkpJCKBSipKSkXpMFu92O3+9H07Q6QSGLxcLTTz+NEAJd11myZAl5eXm89NJLjBo1ioyMjMN+bgpFa6Olp8ntgd3VXu8BYoE680ohRDJwG/DXhm4mhJgqhFguhFjes2fPOu/ruk5+fn7YrLQ+vf3pp5/o3bs3LpcLTdNYt25djeBGJNQntE6nk5iYGJxOJ9HR0URFReFwOMKjsejoaNLS0khISCAlJaXBUWjVlNtkMhEVFUVaWhrJyck1hLr2CHTixIk888wzVFRUMH/+fM4++2zatGkTfr8qidxkMpGYmFjvGqfD4SAlJYWkpKR6n6/KsMLhcDB69GgyMjJ4+eWXSUhIiOQjUyhaHS0thg3dr75h0VTgAynl9oZuJqV8QUo5UEo5sL61t+q2+V6vt47xqpSSn3/+mZycHPx+PxaLhWuvvRbgiMRQCEF0dDRgCE9VsKSqosXhcNRYY9M0jeTkZBITE2u021AQx2q11lmjq/78Q4cO5eyzz2bixIk8//zz4WcCY/pb/Z5CiDqRbCFEOLhT9QyRjpIPFzxSKForLT1N3gUMrva6LVAspayo59xLgJuPpLHa0z+v10tZWVl4ZOT1elm1ahVnnnlmeFTYvXt3HA5HROuFVZjNZkwmU3j6Gh0dXUM8qo+WpJT1rtVVtwcrLS3F5/M1OspyOBxhs4Sq116vN/zMt99+OxdddFGNvD8hRL2J3FarlZSUlPBeKnFxceFAT9VWBTExMYfdya/qGRSKE5GWHhkuAE4VQnStfH0t8EHtk4QQCUAW8G0Lt4/L5Qr7/RUVFbFhwwY6deqEw+Hg448/5uyzzwaaNjKEQyOz6qPC+mgsaFH1fnx8PKmpqYftQ23hiYuLC4/MnE4n/fv3r5EAHRMT0+DIzWw2k5SUREpKSp3cyarns1qt4RSf+uq1a3sjKhQnEi06MpRSHhRCXAW8K4SwAluBSUKIgcCLlRFkMIRwv5TyqOwwVFpaSkVFBdu3bycpKYmoqCjsdjuTJk0K2903VQxjYmLweDyNCk5TiERUbDZbjU2mqsSrKmG8OiaTqVGRbqyNxMREXC4XnTt3pqysDF3Xw9sTVA8cqVGh4kSmxS28pJTzgfm1DhcBfaud8xOGIB41gsEgW7ZsITMzEyEEu3fvJjY2ltTU1GatfQkhSE5OPqZrZlarFU3TarhMm81mUlNT62zr6XQ6mz1qs9vtYaHTNC3s3B0XFxcu3TOZTEoMFSc0J3QFysaNG+nYsSNms5l77rmHr776CqBObmCkHI/gQX0ClJiYWCNIYrPZmjUqbKi9qooZk8kUjnZHR0erKbLihOaENnfdsGEDvXr1YsWKFezYsSPsQ9iU4Mnxxul01kiwrgp8JCYm4vV6CYVCh61yaSpVglhUVITT6cRsNoeTz1vTZ6dQNIUTVgyllGzZsoXf//73PP/889x2223hEWFr+oU2m81ERUWFBbF6IKelRbA6NpuNtLQ0NE0LJ7ZXGUkoFCciJ6wYVgVQOnXqxMyZM2tUTbQmMQRj7a4qredoCmBtqtKHVEqN4mSg1a8ZVm2mXvvYjh07sFqtvPnmmzUiwEKIJkeSjzdCCJKSkkhKSoooOVqhUDSdVv2bVVJSwoQJE/j9738f3l8YjGTrdevWkZyczOLFi3E6neH31DRPoVDUR6sVw2AwyGWXXUZ2djZnnXUWF110UXizc5fLxZdffklBQQG33XZbjSiwmu4pFIr6aF3zxWrMmDGD2NhYpk+fHs4jvO+++7j77rsJBAJ8/fXXZGdnM2LEiBrXHcs1N4VC0XpolWLo8/mYN28eTz75JMFgEI/Hw2233cbFF1/MW2+9RVpaGklJSTz99NM11gcbMzVVKBQnN61SDG02G8uXL2flypXhjc5NJhMzZ87kz3/+M8XFxZx33nl1Kk2aaryqUChOHlrtMKm+aogOHTowc+ZMnE4nt99+O3AoPcRsNqspskKhaJBWK4YN8eWXX3LOOeeEAyVVI8PG9i1RKBSKVjlNboyuXbuSnZ0dfl1lva9GhQqFojFOGDFcv349Gzdu5Pzzzw8fq1ozrM/wVKFQKKpzQkyTV65cyY033lgnQGKxWEhKSmp1FScKheLY0+pVQkrJE088wW233cbo0aPDx6s2PFIVJwqFIhJa/cjQ7/dzxhlnMGbMmBrHo6OjlRAqFIqIafUjw4MHD3LllVfWOGaxWMIefAqFQhEJrXpkuGrVKm644QZsNluNvMOqtBo1MlQoFJHSasXQ5/Px3HPPccUVV5CQkBDe5a3KDFXTNLXHr0KhiJhWO02+5557iImJYcqUKQDhbS4dDoeyp1coFE1GSCmPdx8iYuDAgXL58uXh1y6Xi4MHD4anxFX7EVffKL2lNklSHBa1U5Si1dNqR4ZOp5Py8nKqxNzhcOBwODCZTLjdbiWECoWiSbRaMaxN1daZVqtVBU4UCkWTabUBlNpU30dYoVAomsoJIYZWq1WZtioUiiPihFAQNSpUKBRHygkhhmqNUKFQHClKDBUKhYKjIIZCiHOEEKuEEBuFEO8IIWLrOSdXCPGVEGKFEGK5EGJAc9uzWq31bgGgUCgUTaFFxVAIkQLMAsZLKbsD24D/q3VOFLAAeFRK2Q/4B/B6c9tUo0KFQtEStPTIcDTwk5Ryc+Xr54DLRM2h22hgq5RyfuXrD4EJzW1QBU8UCkVL0NJi2B7YXe31HiAWqL4bUzcgTwjxkhBiObCQZiZ/WywWJYYKhaJFaOkKlIbENVTtewswDjhDSvmDEOIPwHwhREcppa/6RUKIqcBUMLYBrY3T6WyRTisUCkVLjwx3ARnVXrcFiqWUFdWO7QM2SCl/AJBSfgCYgM61byalfEFKOVBKOTAlJaVOY1UmDQqFQnGktLQYLgBOFUJ0rXx9LfBBrXM+BTKrIshCiBGABLa3cF8UCoUiYlp0miylPCiEuAp4VwhhBbYCk4QQA4EXpZR9pZR5QojzgWeFENGAD7hQSultyb4oFApFU2i1foaK3xQq0VPR6jkhKlAUCoXiSFFiqFAoFCgxVCgUCkCJoUKhUABKDBUKhQJQYqhQKBSAEkOFQqEAlBgqFAoFoMRQoVAoACWGCoVCASgxVCgUCkCJoUKhUABKDBUKhQJQYqhQKBSAEkOFQqEAlBgqFAoFoMRQoVAoACWGCoVCASgxVCgUCkCJoUKhUABKDBUKhQJQYqhQKBSAEkOFQqEAlBgqFAoFoMRQoVAoACWGCoVCASgxVCgUCkCJoUKhUABKDBUKhQI4CmIohDhHCLFKCLFRCPGOECK2nnMeF0LsEkL8Wvk1p6X7oVAoFE2hRcVQCJECzALGSym7A9uA/6vn1KHARCll38qvS1qyHwqFQtFUWnpkOBr4SUq5ufL1c8BlQghRdYIQwgb0A/4mhFgphJgrhOjQwv1QKBSKJtHSYtge2F3t9R4gFoipdqwNsBi4E+gLfA98UF0wFQqF4lgjpJQtdzMh7gI6SCmvrXxtBgKAU0pZ0cA1AigF+kgpt9d6byowtfKlXUrZq5G21wDeI3+KRkkGCk6ANhr9LBWKkxFzC99vFzC42uu2QHF1IRRC9MYQvlernScwRLMGUsoXgBcibNsrpRzY9C5HjhBi+YnSxtG8v0LRGmnpafIC4FQhRNfK19cCH9Q6RweeFEJ0qnx9HbBKSrmnhfuiUCgUEdOiI0Mp5UEhxFXAu0IIK7AVmCSEGAi8WBk5XiOEuAn4SAhhwlhX/GNL9kOhUCiaSktPk5FSzgfm1zpchBEsqTrnNeC1Fm460um0auPYtKFQtCpaNICiUCgUrRVVjqdQKBS0MjGMsNTvsOe0QBtHXE4oDF4RQvytuf1ogTZUWaRCUUmrEcNISv2aUA7Y7DYqOaJyQiFENrAImHCE/Wh2G5WoskiFopJWI4ZEUOoX4TlH1EYLlRPegCF2bze3H0fahiqLVChq0prEMJJSv0jOOdI2jricUEp5Y62k8+b040jbUGWRCkU1WpMYNtTXUBPPOaI2pJTbpZTjpJQbpRGK/xfQBciMsI0W6ceRcoyeQ6FoNbQmMdwFZFR7XafUL8JzjqgNIURvIcQVta6rt5zwCDjS5zgsx+g5FIpWQ2sSw0hK/SI550jbOBblhEf6HJGgyiIVimq0eAXK0SLCUr96z2nhNo5KOWFLPkeEbaiySIWiGqoCRaFQKGhd02SFQqE4aigxVCgUCpQYKhQKBaDEUKFQKAAlhgqFQgEoMVQoFApAiaFCoVAASgwVCoUCgP8PWYR+pFPLaNEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 324x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "time_bins_seq = np.load(os.path.join(out_path,'time_bins_seq_gauss_%s.npy'%dates[0])) # time vector whole seq \n",
    "Proj_first3 = np.load(os.path.join(out_path,'Proj_first3_%s_%s.npy'%(exp,soa)))\n",
    "# plot block type effect for matched trials (precedet by xY or xx) in xY or xx blocks\n",
    "# for xY trials, get local deviance response and blocktype response in following trial\n",
    "nperm = 0\n",
    "nsess = len(np.unique(Session))\n",
    "\n",
    "last_stim_on = np.argmin(abs(time_bins_seq-3*.6))\n",
    "nbins = last_stim_on\n",
    "Proj_block = np.zeros((nsess,2,2,nbins))\n",
    "\n",
    "block_ax = 0 # which dimensions in projection matrix to use\n",
    "AUC = np.zeros((nsess,2,nbins))\n",
    "RandAUC = np.zeros((nsess,2,nbins,nperm))\n",
    "AUCdev = np.zeros((nsess,2,nbins))\n",
    "RandAUCdev = np.zeros((nsess,2,nbins,nperm))\n",
    "\n",
    "\n",
    "for sess in range(nsess):\n",
    "    # store trials to compute AUC later\n",
    "    AllProj_block = [[],[]]\n",
    "    AllBlock = [[],[]]\n",
    "    for c in range(2): # xx or xy block\n",
    "        for c2 in range(2): # previous trial xx or xY\n",
    "            # complicated because only should include trials preced by a completed trial within the same block\n",
    "            last_item = np.where((df.ItemID==3) & (df.date==dates[sess]) )[0]\n",
    "            previous_trial_completed = (df.ItemID.iloc[last_item-4]==3)\n",
    "            previous_trial_same_block = np.array(df.blockType.iloc[last_item-4])==np.array(df.blockType.iloc[last_item])\n",
    "            \n",
    "            # choose block of current trial and local deviance of previous trial\n",
    "            trial_block = Blocktype[Session==sess]==c\n",
    "            \n",
    "            # apply all criteria\n",
    "            inc = np.array((trial_block[1:]) & (Local[Session==sess][:-1]==c2) & previous_trial_completed[1:]  & previous_trial_same_block[1:]  & (Include[Session==sess][:-1]==1))\n",
    "            inc = inc[Include[Session==sess][1:]==1]\n",
    "            inc = np.where(np.concatenate((np.zeros(1).astype(int),inc)))[0] - 1            \n",
    "            #average\n",
    "            Proj_block[sess,c,c2,:] = np.mean(Proj_first3[sess,block_ax,(Session[Include==1]==sess),:][inc,:last_stim_on],0)\n",
    "            # all trials\n",
    "            AllProj_block[c2].append(Proj_first3[sess,block_ax,(Session[Include==1]==sess),:][inc,:last_stim_on])\n",
    "            AllBlock[c2].extend(c*np.ones(len(inc)))\n",
    "            \n",
    "\n",
    "    # compute AUC\n",
    "    for c2 in range(2):  # previous trial xx or xY\n",
    "        AUC[sess,c2,:] = getAUROC(np.array(AllBlock[c2]).astype(float),np.concatenate(AllProj_block[c2]))\n",
    "        # random permutations\n",
    "        for perm in range(nperm):\n",
    "            randind = np.random.permutation(len(AllBlock[c2]))\n",
    "            RandAUC[sess,c2,:,perm] = getAUROC(np.array(AllBlock[c2]).astype(float)[randind],np.concatenate(AllProj_block[c2]))\n",
    "    \n",
    "# 1. plot trajectories from above but for same global context, different global deviance condition in previous trial\n",
    "threshold = 3\n",
    "alpha = 0.05\n",
    "\n",
    "colors = ['#165aa9','#ee8c22'] # block type colors\n",
    "ls = ['-','--'] # standard or deviant before\n",
    "labels = [['xx block after standard','xx block after deviant'],['xY block after standard','xY block after deviant']]\n",
    "binsize = np.mean(np.diff(time_bins_seq))\n",
    "\n",
    "fig,ax = plt.subplots(2,1,figsize=(4.5,4),sharex=True)\n",
    "\n",
    "# PLOT PROBABILITY\n",
    "for c in range(2): # xx or xy block \n",
    "    for c2 in range(2): # previous trial standard or dev\n",
    "        av = np.mean(Proj_block[:,c,abs(c2-1*c),:],0)\n",
    "        sem = np.std(Proj_block[:,c,abs(c2-1*c),:],0)/np.sqrt(nsess)\n",
    "        \n",
    "        ax[0].fill_between(time_bins_seq[:last_stim_on]+binsize,av-sem,av+sem,alpha=.2,linewidth=0,color=colors[c])\n",
    "        ax[0].plot(time_bins_seq[:last_stim_on]+binsize,av,color=colors[c],linewidth=1,linestyle=ls[c2],label=labels[c][c2])\n",
    "ymin,ymax = ax[0].get_ylim()\n",
    "ax[0].legend(bbox_to_anchor=(1,1))\n",
    "for i in range(3):\n",
    "    ax[0].fill_between(np.linspace(i*.6,i*.6+.6/2,2),\n",
    "                               np.ones(2)*ymin,np.ones(2)*ymax,alpha=.05,color='k',linewidth=0)\n",
    "\n",
    "    \n",
    "# PLOT AUC BLOCK TYPE, MATCHED IN LOCAL PRE\n",
    "labels = ['preceded by xx','preceded by xY']\n",
    "for c2 in range(2): # previous trial xx or xY\n",
    "    av = np.mean(AUC[:,c2,:],0)\n",
    "    randav = np.mean(RandAUC[:,c2,:],0)\n",
    "    sem = np.std(AUC[:,c2,:],0)/np.sqrt(nsess)\n",
    "        \n",
    "    ax[1].fill_between(time_bins_seq[:last_stim_on]+binsize,av-sem,av+sem,alpha=.1,linewidth=0,color='k')\n",
    "    ax[1].plot(time_bins_seq[:last_stim_on]+binsize,av,color='k',linewidth=1,linestyle=ls[c2],label=labels[c2])\n",
    "    # stats:\n",
    "    p = clustertest(av,randav,threshold)\n",
    "    l = label(~np.isnan(p))\n",
    "    for j in np.unique(l):\n",
    "        if p[l==j][0]<=alpha:\n",
    "            print(p[l==j][0])\n",
    "            ax[1].plot(time_bins_seq[:last_stim_on][l==j]+binsize,np.ones(sum(l==j))*.6+.06*c2,color='k',\n",
    "                       linestyle = ls[c2],linewidth=1.5)\n",
    "            \n",
    "    ax[1].set_ylim(.5,None)\n",
    "\n",
    "sb.despine(trim=True)\n",
    "plt.tight_layout()\n",
    "ax[0].set_ylabel('Global context (up=xY)')\n",
    "#ax[2].set_ylabel('Effect of preceding global deviant (context AUROC)')\n",
    "#plt.savefig(os.path.join('plots','Context_BlockMatched_%s_%s_%s.pdf'%(exp,soa,animal_id)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Effect of distance from global deviant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ttest_relResult(statistic=2.3641092260419474, pvalue=0.01855779827187455)\n",
      "Ttest_relResult(statistic=-0.2352891859978196, pvalue=0.8141065301474748)\n",
      "Ttest_relResult(statistic=-1.5049951256765712, pvalue=0.1331265390401334)\n",
      "Ttest_relResult(statistic=-5.02123687109009, pvalue=7.575394870271354e-07)\n",
      "Ttest_relResult(statistic=-1.7728522760517513, pvalue=0.07697108177697776)\n",
      "Ttest_relResult(statistic=-0.06299914097394924, pvalue=0.9497968761999898)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOkAAADQCAYAAADiQkaGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa2ElEQVR4nO2de7xVZZnHvz+uB4lKEQeNEFTSMswUNXUybKZsUKdIG++IyWDmqFnKGDMUWqIC4mQ2IpkammleQdGwksrEy4ACokhqkh0VRLyCcbg988f7bs4+m332WWdf19nn+X4++7P3Wut9n/Wstfaz3tvzvo/MDMdx0kuXWivgOE5h3EgdJ+W4kTpOynEjdZyU40bqOCnHjdRxUk63WitQZXy8yUkrau2Al6SOk3LcSB0n5biROk7KcSN1nJTjRuo4Kaez9e52OMaNG8fKlSvp378/kydPrrU6Tg1wI005K1eu5JVXXqm1Gk4N8equ46QcL0mz8Kqlk0bcSLPwqqWTRtxIK0QaS+U06uS0Tc2MVNKRwKVAT2AJcLqZvZuT5grga8CbcddyMzsuHhsPjCJcw83ARZaitWDSWCqnUadyUq8voZoYqaR+wA3AoWb2vKTLgcuAb+YkPQQ43szm5+QfQTDe/YHNwFzgWeBXldbdSS/legmlzdhr1bv7ReD/zOz5uH0NcJKkrTMBJPUEPg2cL2mxpDslDYyHRwK3mNk6M1tPMPiTq6i/U8dkjH3lypW1VgWonZF+FPhb1nYj8EGgT9a+XYCHgO8C+wKPAbOiIefLPyDfiSSNlbRA0oK99967bBfgONWiVkba2nk3Z36Y2UtmNsLMlse25lRgd2BQK/k359mHmc0ws2FmNqxXr14lqu041adWHUcvAwdlbX8EeMvM1mV2SNoH+JSZ3ZSVTsDGmH/nnPyNlVO38vzlptPz7t/43qqt362l2e2Un1VML6c0ytG+rVVJ+iDwGUlD4vY3gFk5abYAV0kaHLfPBJaYWWNMe5Kk3rHtOhq4p+JaO047KUf7tiYlqZm9Luk04A5JPYAXgVGShgHXmdm+ZrZU0tnAvZK6EkrKE2L+eyUNBZ4AehCMdmYtrsVxKk3NxknN7H7g/pzdbxI6iTJpbiaMgebLPwmYVCn9HCctuMdRndFauxXabt/Wom2btjHJNNIpjbTYThrvoCk/7XVA6GgvoXLQKY20XHTGP0w90VF61H0+qeOkHDdSx0k5bqSOk3K8TZpydvxAtxbfTvqodN+EP/mUc94X+9daBafGeHXXcVKOG6njpByv7jpVwR1IiseNtBNR751Q9Xp99XU1TkHqvROqXq/PjTSLcr6J6/Wt3hlI27NLhxYpoZxv4np9q3cG0vbs2jRSSd/L2WXA+8BSM5tbEa2cVOPTy5JTjlI5Sc6hhPVv7yAs9jUSWAH8m6QDzewHRZ/d6ZDU+yLb5aQcpXKScdJ/APY3s3PN7NvAMEJp+lng30rWwHGcgiQx0r5mtnUVJTNbE/dtIKzc5zhOBUlS3f2LpEuBGYQlNU8HXpR0EK2sdes4SUlbT2oaSXJnTgN+DDwFbALuA8YAxwHnV041pzOQtp7UNNKmkZrZG8SlNHO4pvzqOI6TS5IhmM8BE4EdCNVdAMxsn8qp5ThOhiTV3Z8A1wNPEnp1y0LC+KQnAxfQPDZ7jpktiMcWAr2ADTH5L8xsSrn0c5y0kMRIN5jZtHKeNEl8Ukl7AlOA/czstRiT9C5goKTehOBN/czMe5iduiaJkS6VNNTMni7jefPFJ10s6aysaN1NwBgzey1uLwD6x7AUBwJrgTmSdgZ+C4w3s7+XUcdOT0dZ8rLeSTJOuhuwUNLzkpZkPiWet834pGa2wszmAMSYpNOA2XF8tg8wDzgWOAAYSKg6b0N2fNLVq1eXqLbjVJ8kJel/VeC8bcYnzRCrtjcSDPtLAGY2G5idlWYSoSr8rdz8ZjaDMMbLsGHDytamdpxq0WpJKmmv+PO9Vj6lkC++aIv4pFGHgcB8gvEebmZvx/1HSzosOynu/eTUKYVK0iuAI4E78xwzQjW4WB4ErpA0JLZLt4lPKmkH4A/AjWZ2UU7+AcCEODy0Afg2cFsJ+jhOamnVSM3syPg9uLU0xZIkPikhaPBAYKSkkVnZ/wm4lvCSeDJewzzg4nLr6ThpoFUjzTOPtAVmVpJRtBWf1MwuAS4pIOKC+HGcuqZQdbdf/N4L2BO4m+C7+2WC84HjOFWgUHX3bABJDxEcCt6I2z8kp/3oOE7lSDIEs3PGQCNvAztVRh2nI+DTy6pLkru8RNINwEya55M+XlGtnFTj08uqSxKPozGE0vNHBK+fRsKQieM4VSDJfNL3JE0gtEO7AY+7j6zjVI8k80kPIBjoKqArMEDSUWY2v9LKOY6TrLp7BXCSmX06TvQ+llDtdRynCiQx0g+a2bzMhpk9BGxXOZUcx8kmiZFukbRrZkPSIHyVQMepGkmGYC4GHpP027j9RbJWUHDSwcaNG2lsbGTdXqPJWoqqXSxbtqzF9qa9Titan3LJyshpaGhgwIABRevTkUnSu3uPpGXA5wkl7yVm9lzFNXPaRWNjI3369GHnPgMIc+TbT8++g1psN61ZUbQ+5ZLVs+8gzIw1a9bQ2NhYtD4dmTaru5L2AaaZ2TXAH4Fb4/pDTopYv349ffv2LdpA04wk+vbty/r162utSk1I0ia9BrgOIK5zNJEwVcxJGfVooBnq+draIomR9jazuzMbZnYPYT0ix3GqQJKOI8teLVDSx/He3Q7BiVPb52+ibgtbbNumpoLpbzn/kHbr5LSfJEY6AfijpMySnnsBJ1VOJcdxsknSu3tf7Cg6lDDp+3Eze73imjmOAyRrkwJ8PH53B/5R0lcrpI9TZ/z85z9nt912Y+3adaxb9z6fOvifuenWO9lr/8Na7Lv5tnzr3TmQzMF+JmGM9AWaY8EYYZ1bxynIqaeeyoMPPsj4iy6jacMGDv3MAZxy/DH8dt7DLfadfNwxtVY1tSRpk34W2MvM1lZaGac+mT59OvsM3ZteDQ08+rt7Abh66g85YPiIFvuc/CSp7r7sBuqUwqpVq2ha38Q777zLqytDHJnXV7+xzT4nP0lK0kck3QrcC2yd7G1mXt112mTjxo2ccMIJTLjwPLZs2cKpY8/hoTm3M2rsudvs6969e63VTSVJjPTg+D0ma1/JbdKE8UnzppHUlTCn9Yh4DVPNbHop+tQj7R3HLKfvbobx48fTv39/vn7K8QDc+8Bv+N4Pp/APO+3YYt/3L5nKpInfLfl89UiSIZjDy33ShPFJC6U5AxgCfJIQYe1RSU+a2RPl1tUpjSlTQlznjMHPuvWGbdLk2+c0k8TBvr+kOZL+LGknSXMllbpcXL74pCeppYNmoTQjgRvMbJOZvQXcCpxcok6Ok0qSdBz9L3APoT36FrAIKDVCbJvxSdtIk+9Y55xs6NQ9Sdqkg8zsp5K+aWYbgf/MchEsliTxSQulyXcsrz+xpLHAWICBAwcCMO7PIxIpmcsdOdvFyqmkrMb1Hypazu4522mQlSsnjfe8nDrlI+nyKVvTSeqTMF8hksQnLZQm37G8M4LNbIaZDTOzYf369cuXxHFSTRJjuwv4BfAhSWcADwG/KvG8DwKfkTQkbm8Tn7SNNLOAr0vqJunDwPGEKrnj1B1JencnSTqFYNBfIIS2v66UkyaJT9pamijiGkJNaDHQA7jWzP5Qik6Ok1YSRdwxs5uAm8p54rbikxZIg5ltAr5VTn3qkQ1zz2tX+r9079oy/8bC04Z7HHFlu3Vy2k+pbUvHaTcT/2scxxz9BTZvbn4JbN68ma99+QimTS4UN7pz4kbqVJ3vTriY999fx/Srm0vi6VdfSbdu3Tj3OxfWULN0ksSZ4SN59n2iMuo49UZmPum6dWt5//11HDH8IO6/7x7+5+rrmHHNVTz7zNM8s3QJv5h5PVde/VO6du3attBORqttUkk7xJ/3SxpO84rL3Qm9q0PyZHOcFmTmk15+yUQ2NDUx7MCDGXls8Nn9zrj/ZvwF57J582YunfIj+u+8S421TSeFOo5+SejNBViTtX8TcPe2yR0nP9OnT2fvTw6lZ0MDs+7fGlaIUV8fywNzZjPkY3vyuc9/oYCEzk2r1V0zO8LMugA3mlmXrE8PMzuuijo6HZxVq1bR1NTEe++8y6pVK1sc++jAXRm466DaKNZBSNJx9FD2hqSuki6rkD5OnZGZT3rudy7k7PMu4LyzxrBx48Zaq9WhSGKk4yVNl9RD0mDgEcIUMcdpk8x80uNOHMXxJ4/mw9vvwLTLf1hrtToUSZwZDgB+BCwEPgxcZGYleRw51aG9zga7fWT7FtsvvvJWyTpk5pNmZP3sppYepZOv/EnJ56h3kpSkBjQRAgd3AbZUVCPHcVqQxEifJszj/BTwOeAbknKd4R3HqRBJjPQiMzvFzNaa2QvAIcDyCuvlOE6kTSM1s5mSDpR0RpyNcoCZjauCbk47MbO2E3VQ6vna2iKJW+BowoJg4wgdR7Mk/Xtl1XLaS0NDA2vWrKnLP3Mm0ndDQ0OtVakJSXp3zyEs6/mHOMdzf+DXwE8rqpnTLgYMGEBjYyOvvLam6IC7G97drsX26rffL1qfcsnKyGloaGDAgAHAk0Xr1FFJYqSb41q3AJjZ3yRtqqxaTnvp3r07gwcP5oLri/8T3/GDlvFYJkwoPohSuWTlyumMJOk4elPSvsRgTZJOIkzOdhynCiQpSb8F3A7sLulVYD3w5Uoq5ThOM0mM9DnCGOnHgK6E4ZcPVFIpx3GaSVLdXWhmm81smZktjWvvPlJpxRzHCRSa9P07gt/udpKyAyl1BZ6qtGKO4wQKVXdHAjsA1wOnZe3fBLxWSaUcx2mmVSONYQjfBT4fPY1607yEyofxHl7HqQpJPI7OAd4B3sj6rC7lpJKOlLRE0nJJt0v6YCvpTpa0WNIiSfPj4tmZYwslPRuPLZJ0QSk6OU5aSdK7ey4hRmhZXD2SxCaN6fYEpgD7mdlrkkYQQl4MlNSbsIJ9v9iR5Th1S5Le3dfKZaCRJLFJIcxhHWNmmfbvAqB/rHofCKwF5kh6WtKVknqVUUfHSQ1JjPQ3ks6U9BFJO2Q+bWWSNELSptwPoQRsKzYpZrbCzOZEWQKmAbPNbENMOw84ltADPRC4NMG1OE6HI4mRXgj8hGBYidukZna/mXXL/RB6h/PRWnzR3oQobnsAY6Ls2XGO67tmth6YROiNzpd/rKQFkhasXl1SU9pxakKS+aS9cpb07GJmpSwzniQ2KQCSBgLzCQZ8uJm9HfcfLemw7KRA3rapxyd1OjpJene7SBonaZ6kP0n6nqRE0dhaIUls0swK+n8A7jKz483s71mHBwBTJfWS1BX4NnBbCTo5TmpJUt29FPg8YcXAaYTlU6YWe0Ize53gHHGHpGXAUOA7AJKGSVoUk55JaGuOzBpmWSSpL3AtwYCfJPgWrwUuLlYnx0kzSUrELwHDMkMdkuYQgvcWTYG4owuI8UnN7BKgUBy8C+LHceqaJCVpl+yxSDNropX2n+M45SeJkS6K45C7x880YEmlFXMcJ5DESM8Ctif0sj4G9APOrqRSjuM002abNDrajwaQ1BDHJR3HqRKF5pP2IKwIeLeZ3RN33yFpNfDvZuaLkXVS3ntmNlua3qNLzz702ftfa61O3VOounsxwV1vfta+MwhV34kV1MlJOVua3mPL+nfY0vRerVXpFBQy0qOAE+O4JgBm9gowilZc8BzHKT+FjHRDjpcPsLWN2lQ5lRzHyaaQkW6W1Cd3Z9zXvXIqOY6TTSEj/SVwXZyFAmydkXIdUPzS5o7jtItCRvo/hGVTVkp6TNITwErgLdxP1nGqRqGFyLYAYyVNAvYjRPh+wsxerZZyjuMkc2ZYAayouCaO4+SllHmhTp3TWkSzUaNm8cor77Bz3w8w06OeVZwkvruO49QQL0k7Ee7O1zFxI+1EZNz50kQaXxxp08mN1KkpaXxxpE0nN9IKkba3sdNxcSOtEGl7GzsdF+/ddZyU4yVpyvFqs+NGmkUaDcKrzclJ4/Mrh05upFm4QSSjf//+Lb7TQhqfXzl0qomRSjqSsDJ+T8LyoKfHyeS56a4AvkZzVPHlZnZcPDaesEpEN+Bm4CIzsyqon2pac+WD8rnzTZ48uei8aaKjuD1W3UiTBhGOHAIcb2bzc2SMIBjv/oRgTnOBZwnR15wU0lEMIo3UoiTNF0R4saSzsktCST2BTwPnS9odeAE4z8xeJqyxdEsmEpukG4CTqbKRVqPUcpyKDcGUGkQY2AV4CPguIT7MY8CsGFD4o3lkDGhFj63xSf/+922WbHKc1FOxkjQGZdpGfmxL5qNFEGEzewkYkZVvKjABGET+l0veIMRmNgOYkdlsS2/HSRu1cGZIFERY0j6STsnJmwkWnE9GYwV0dZyaU4s26YPAFZKGxHZp3iDChOVarpL0p1iqngksMbNGSbOA70uaAWwihMG4sSraVwjvWHFao+pGamavS8oEEe4BvEgYSkHSMOA6M9vXzJZKOhu4N0bzbgROiDLulTQUeALoQTDymUl1cIPouKSxs67SOtVknDRJEOG4fTNhDDSfjEnApAqp6DipwT2OOhFp9RRKG2m7T26knYg0egqlzSAgfffJjbRCpPHPl0bSZhBpxI20QvifzykXPunbcVKOG6njpByv7maRxnZkGnVKK2m8V+XQSZ1sCmanulinQ6HWDnh113FSjhup46QcN1LHSTlupI6TctxIHSfldLYhmFZ70LYmkJYC68t0vh2BN1Ikp5yy0qhTOWWVU6cGM/tksZk7m5EmYb2ZDSuHIEkLyiGrXHLqXadyyiq3TqXk9+qu46QcN1LHSTlupNsyo+0kVZflOlVfVmp06mxugY7T4fCS1HFSjhtpFpL+Q9IzkpZKmiVppxJknSxpsaRFkubHlRCLlSVJN0o6vwQZR0paImm5pNslfbCW+kQ55bxHZXt2Ud5XJG0TRKydMq6Q9HK8vkWSbitKkJn5J1T59wdWAB+K21OBa4uUtSfwGrBz3B4BvFykrI8Twm28D5xfpIx+wOvAkLh9OfC/tdKnAveobM8u5h9CiD20tsT/1KPAIaXIMDMvSTOY2ULCn/gdSQ2EVfHXFCmuCRhjZq/F7QVA/7jOcHs5ixCFrpRgVPmCZJ0U4+rUQh8o4z0q57OTtB1hGdlvF5M/S052wLHFku6UNLAYWW6kWZjZRklfISzEfRjhz1iMnBVmNgdC1RCYBsw2sw1FyPoPM7upGD2yyBfgKl+QrGrpU9Z7FOWV5dkB18bPkiLzZygUcKxduJHmYGb3mNmOwERgrqSi75Gk3oQSZw9gTHk0LIrWriFvkKtqUs57VOqzk/RNYJOZXV+KHlGXl8xshJktt1D3nUqIKDiovbI6tZFKujirUX+LpH/MOnw9sCuwfRGyLo5Vm/kEQzjczN4uRk47L6k1EgXJqjbF3qM8cvYo5dllMRo4QNIiQoSFXvE57FKEToUCjrWPUhu19fIBPkuoKu0Yt0cBi4uUtQPwEvD9Mup3I8V3HO0ErKK54+hS4IZa6VPue1TOZ5clcxAldBwBnwTeAgbH7W8CjxQjyx3sI2b2sKRLgN8rBDt+FfhKkeLOBAYCIyWNzNr/T2ZWbGdU0ViBIFk1pGz3qMzPrixYgYBj7cU9jhwn5XTqNqnjdATcSB0n5biROk7KcSN1nJTjRuo4KceNtAQkXZXleLAhzjDJbPfKSfsNSRe2IW+4wkJo+Y5dI+mlONRQEwrpl5POJO3YTtm/l3RsgnRLJQ1vj+ysvP8q6api8sb8gyXdWWz+YvFx0hIws3MyvyWtAE4ys7yLTpnZ9BJPdwYw0MwaS5TTaTGz2cDsEkTsSpi9U1W8JK0QkiZKmhvncN4ct6+Ox46K8ycXxPmGP2hD1sMEl7IHJH1W0gpJt0laJmmkpL0lzYvnWixpVMw3XNKjcQbGc5KelHS0pN/E817ZyvkOjLotkXR3zDc8J82H4nUtlfS0pMmSsl/6l8R8iyQdFfP0ljRT0mOS/ixpoaSCf3pJn4jpl0j6FdA769ghkh6O51mQdZ752aWypMskXS5ptKT74r7PSPqjpMfjvfhZ3D9I0ouSfizpCUkvSDouOiRcB+wuaW4hnctOqS5Z/tnqBrYCGJa1PRF4DuiWtX01wdjm0eyitwuwibDO63BgaSvyjWa3txXAhPi7G8GD6KtZ8hqBg6O8TcCn47EHCL6yPeL5NgC75JynG2HGzL/E7cOBLVHWVv2AnwM/itfTE5gLXJil639mucetIcxpPRa4Kutc04Efx9+/B47Nc91PAafH34dm6bI9sBwYlHXdfyN4MZ0G3Bf3Z7x9hhB8czP7fwkMj78/AKwmzEsdFPU/Kh47Bvhr/N3q86nkx6u7leUxM9uUvcPMTNLRwFGSTiRMohZZJURCHo7fHyMsvnxXlP9qbDd9ifAyeMnMnoppXwTesTAd7A2FlQd2ILjRZRga5TwQv+e10g79F+BQC//eJknTgW8Bl8Xj02P+pZKeBQ42szsk/SW6y+1B+NM/2toFSuoL7APMjLIeydLlYMKkgXvUPPvLYvpfAVMl9Qf2A14ws+clHZol/lRghKTxwF7AdgRjXUNwgr8/pnsy3qOa4UZaWdbm7lCYmvUUcDfB0K4n+Jm2d55hRna+JksXoHv83ZRzrK1ZGJvy6JJvSlvuebPPmZtHwEZJZwJjCTWKW4A3gcEFdMn4rGbrk3npdQWWmdlBW08SZqustjC39HbgRIIx/zSP7IeBxcCvCUZ9UNZ5NpjZliwdipkcXza8TVp9hhAmXP+3md0LfI5QXexapLzlwAZJX4Wtf9RjgN8UKW8ZoWT8UpR3IKF0zXXyngucpUBPgvFln3N0zL8f4ZofB44AbjSzn0W9j6bAdZvZm8BC4jzTKGtoPPwYMETSYfHYvsDzhGovBMMcDRwCtOiRlbQ9MIxQJb+LMHVvj0K6RDbR8kVUFbwkrT5LgPuA5yS9TVhL51nCnyS31GsTa16R4CpJEwnP9OJYTR1ehLxNko4Bpku6FPgzsJKwplF2lfwc4MfA04Q27q+B7OGh3SQ9RTDu483sTUlTgRkKM3I2EwxwKIU5AbghlsIvEF4imNnqqOcUhSVTugCnmNlf4/GFCjNi7jSzFrF9zOyteG1PSlpDiPnyCOEZvFhAl2eAzZKeAA6KVf2K47NgnG2QNAWYamarJH2UUC3czYqclO2UhpekTj7+CvxO0kZCe2yMG2jt8JLUcVKOdxw5TspxI3WclONG6jgpx43UcVKOG6njpBw3UsdJOf8PXC/6Q4Nlpl4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 216x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# get all global deviant trials and plot average global context response as a function of the distance from the global deviant trials\n",
    "trials_around = np.arange(-3,6)\n",
    "# selectrion criteria: global deviant trial and next 5 trials in same session & same block\n",
    "\n",
    "glob_criterion = Global[Include==1][abs(trials_around[0]):-trials_around[-1]]==1\n",
    "noglob_criterion = glob_criterion.copy()# check if no global deviants around\n",
    "for shift in trials_around:\n",
    "    if shift==0:\n",
    "        continue\n",
    "    elif shift<trials_around[-1]:\n",
    "        noglob_criterion = (noglob_criterion==1) & (Global[Include==1][abs(trials_around[0])+shift:-trials_around[-1]+shift]==0)\n",
    "    else:\n",
    "        noglob_criterion = (noglob_criterion==1) & (Global[Include==1][abs(trials_around[0])+shift:]==0)\n",
    "\n",
    "\n",
    "block_criterion = Blocktype[Include==1][:-trials_around[-1]+trials_around[0]] == Blocktype[Include==1][abs(trials_around[0])+trials_around[-1]:]\n",
    "sess_criterion = Session[Include==1][:-trials_around[-1]+trials_around[0]] == Session[Include==1][abs(trials_around[0])+trials_around[-1]:]\n",
    "block = Blocktype[Include==1][abs(trials_around[0]):-trials_around[-1]]\n",
    "session = Session[Include==1][abs(trials_around[0]):-trials_around[-1]]\n",
    "proj = Proj_first3[:,block_ax,:][:,:,int(abs(time_bins_seq[0])/binsize):last_stim_on]\n",
    "\n",
    "resp = [[],[]]\n",
    "for bl in range(2): # loop across contexts\n",
    "    for sess in range(nsess):\n",
    "        resp[bl].append([])\n",
    "        criterion = np.where(glob_criterion & block_criterion & sess_criterion & (block==bl) & (session==sess))[0]\n",
    "        criterion = criterion + abs(trials_around[0])\n",
    "        for i,shift in enumerate(trials_around):\n",
    "            av_resp = np.mean(proj[sess,criterion+shift,:],-1)\n",
    "            resp[bl][sess].append(av_resp)\n",
    "            \n",
    "# pool sessions\n",
    "context_resp = [[],[]]\n",
    "for bl in range(2):\n",
    "    for sess in range(nsess):\n",
    "        context_resp[bl].extend(np.array(resp[bl][sess]).T)\n",
    "\n",
    "        \n",
    "import scipy.stats as stats\n",
    "\n",
    "# STATS\n",
    "# ANOVA with variables trial number, block\n",
    "\n",
    "labels = ['xx','xY']\n",
    "fig,ax = plt.subplots(1,1,figsize=(3,3),sharey=True)\n",
    "colors = ['#165aa9','#ee8c22'] # block type colors\n",
    "for c in range(2):\n",
    "    \n",
    "    r = np.array(context_resp[c])\n",
    "    # ANOVA for each block type\n",
    "\n",
    "    #st = stats.f_oneway(r[:,3],\n",
    "    #           r[:,4],\n",
    "    #           r[:,5],\n",
    "    #           r[:,6],\n",
    "    #           r[:,7],\n",
    "    #           r[:,8],)\n",
    "    #print(st)\n",
    "    # pairwise t-test\n",
    "    for i in range(3):\n",
    "        T = stats.ttest_rel(r[:,np.where(trials_around==1+i)[0][0]],r[:,np.where(trials_around==0)[0][0]])\n",
    "        print(T)\n",
    "    # PLOT\n",
    "    ctx_df = pd.DataFrame(r,columns=trials_around)\n",
    "    #av = abs(np.mean(r,0))\n",
    "    #sem = np.std(r,0)/np.sqrt(r.shape[0])\n",
    "    sb.barplot(data=ctx_df,linewidth=1,label=labels[c],color=colors[c],ax=ax,alpha=.8)\n",
    "    #plt.fill_between(trials_around,av-sem,av+sem,alpha=.1,linewidth=0)\n",
    "ymin,ymax = plt.ylim()\n",
    "plt.xlabel('Trial from global deviant')\n",
    "plt.legend()\n",
    "ax.set_ylabel('Context encoding')\n",
    "sb.despine(trim=True)\n",
    "#plt.savefig(os.path.join(out_path,'plots','EffectOfDeviantOntoContext_%s_%s.pdf'%(exp,soa)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}