YiqunSibilant.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.integrate import cumtrapz\n",
    "from scipy.special import lambertw\n",
    "from scipy.special import gammaln\n",
    "import matplotlib.pylab as plt\n",
    "import os\n",
    "import pymc as pm\n",
    "#run if missing: conda install pymc\n",
    "\n",
    "class AttributeDict(dict):\n",
    "    #extend dictionary object to set atributes as for Matlab objects.\n",
    "    def __getattr__(self, attr):\n",
    "        return self[attr]\n",
    "    def __setattr__(self, attr, value):\n",
    "        self[attr] = value\n",
    "\n",
    "class topologies_and_genes:\n",
    "    def __init__(self):\n",
    "        \"\"\"\n",
    "        Documentation: This is the master class for computing transition, marker genes and topologies\n",
    "        \"\"\"\n",
    "        pass\n",
    "    def set_data(self,tfdata,idxs,labs,iunique):\n",
    "        \"\"\"\n",
    "        inputs:\n",
    "        #tfdata --> log-transformed data. rows = genes. columns = replicates (several per cell type). \n",
    "        #idxs --> which cell types are we interested in out of all of the cell types in tfdata?\n",
    "        #labs --> cell type names\n",
    "        #iunique --> what is the cell type identity of each column of tfdata?\n",
    "        #iunique has same length as the # of columns in tfdata, and the max value\n",
    "        #of iunique is the number of cell types.\n",
    "        \"\"\"\n",
    "        self.tfdata = tfdata\n",
    "        self.idxs = idxs\n",
    "        self.labs = labs\n",
    "        self.iunique = iunique\n",
    "    def start_Sibilant(self):\n",
    "        \"\"\"\n",
    "        #compute integrals and calculate topologies\n",
    "        #inputs grabed from self (after running self.set_data)\n",
    "        #tfdata --> log-transformed data. rows = genes. columns = replicates (several per cell type). \n",
    "        #idxs --> which cell types are we interested in out of all of the cell types in tfdata?\n",
    "        #labs --> cell type names\n",
    "        #iunique --> what is the cell type identity of each column of tfdata?\n",
    "        #iunique has same length as the # of columns in tfdata, and the max value\n",
    "        #of iunique is the number of cell types.\n",
    "        \"\"\"\n",
    "        \n",
    "        print 'Sibilant started.'\n",
    "        #get data from self\n",
    "        tfdata = self.tfdata\n",
    "        idxs = self.idxs\n",
    "        labs = self.labs\n",
    "        iunique = self.iunique\n",
    "        \n",
    "        if not hasattr(self,'Params'):\n",
    "            Params = self.get_parameters()\n",
    "            self.Params = Params\n",
    "        else:\n",
    "            Params = self.Params\n",
    "\n",
    "        loggenemeans, loggenestds = self.calc_log_mean_std(tfdata, iunique,Params.dist)\n",
    "        # internalize\n",
    "        self.loggenemeans,self.loggenestds = loggenemeans, loggenestds\n",
    "\n",
    "        #the function calc_log_mean_std has been modified to be able to calculate\n",
    "        #sigma and mu in ZIP distribution. However, these values will be used only\n",
    "        #if kde is on or filter gene by expression is on. \n",
    "        #######verified up to here\n",
    "        print 'Beginning integration.'\n",
    "        Integrals, best_mu, best_sigma,combinations = self.calculate_integrals(tfdata, idxs, iunique, Params)\n",
    "        #internalize\n",
    "        self.Integrals,self.best_mu,self.best_sigma,self.combinations = Integrals, best_mu, best_sigma,combinations\n",
    "\n",
    "        #combiantions contains all possible combinations of 3\n",
    "\n",
    "        ncells = len(idxs)\n",
    "        ngenes = tfdata.shape[0]\n",
    "\n",
    "        #ncandid = # of likely topologies\n",
    "        #tlikely = most likely root of the triplet\n",
    "        #plikely = probability of T being the root\n",
    "\n",
    "        # Initialize holders\n",
    "        ncombinations = len(combinations)\n",
    "        ncandid = np.zeros(ncombinations,dtype=int);\n",
    "        tlikely = np.zeros(ncombinations,dtype=int);\n",
    "        plikely = np.zeros(ncombinations);\n",
    "        pT_g_sum = np.zeros([ncombinations,5]);\n",
    "\n",
    "        makeplot = Params.get('makeplot',False)#make plots? specify in Params and default to false\n",
    "        plot_dir = Params.get('plot_dir',None) #save plots somewhere if in Paramas, default to not save\n",
    "        #find most likely topologies for each triplet\n",
    "\n",
    "        # Add some progress reporting\n",
    "        print 'Starting trajectories.'    \n",
    "        for j in range(ncombinations):\n",
    "            icomb = j\n",
    "            iii = combinations[icomb,:]\n",
    "            trip_data = loggenemeans[:,idxs[iii]]\n",
    "            tlikely_ind, plikely_ind, ncandidates, pT_g_sum_j = self.process_one_triplet(icomb, Integrals, Params, trip_data, makeplot=makeplot, plot_dir=plot_dir, labs=labs[iii]);\n",
    "            tlikely[j] = tlikely_ind\n",
    "            plikely[j] = plikely_ind\n",
    "            ncandid[j] = ncandidates\n",
    "            pT_g_sum[j,:] = pT_g_sum_j\n",
    "        #internalize\n",
    "        self.ncandid,self.tlikely,self.plikely,self.pT_g_sum = ncandid,tlikely,plikely,pT_g_sum\n",
    "\n",
    "        print 'Starting gene selection.'\n",
    "        #find union of marker and transition genes over all triplets\n",
    "        good_genes = np.zeros(ngenes)\n",
    "        transition_probs = np.zeros([ngenes, ncombinations])\n",
    "        marker_probs = np.zeros([ngenes, 3, ncombinations])\n",
    "        for i in range(ncombinations):\n",
    "            if plikely[i]>Params.plikely_thresh and ncandid[i]>=1:\n",
    "                icomb = i\n",
    "                iii = combinations[icomb,:]\n",
    "                trip_data = loggenemeans[:,idxs[iii]]\n",
    "                prob_asym_topol, prob_marker_topol = self.calculate_gene_probabilities(icomb,Integrals, Params,trip_data,tlikely[icomb]);\n",
    "                transition_probs[:,i] = prob_asym_topol\n",
    "                marker_probs[:,:,i] = prob_marker_topol\n",
    "                genestatus = self.find_classes(prob_asym_topol, prob_marker_topol, tlikely[icomb], trip_data, Params)\n",
    "                good_genes = (good_genes + abs(genestatus))>0\n",
    "        self.good_genes,self.transition_probs,self.marker_probs = good_genes,transition_probs,marker_probs\n",
    "    def find_classes(self,pasym_gi_T, pmar_gi_T, tlik, data, Params):\n",
    "        \"\"\"\n",
    "        Documentation...\n",
    "        \"\"\"\n",
    "        N_thresh_asym = Params.N_thresh_asym; \n",
    "        N_thresh_sym = Params.N_thresh_sym; \n",
    "        p_thresh_asym = Params.p_thresh_asym;\n",
    "        p_thresh_sym = Params.p_thresh_sym;\n",
    "        ngenes = pasym_gi_T.shape[0]\n",
    "        status = np.zeros(ngenes)\n",
    "        # asym\n",
    "        asymgenes = np.where(pasym_gi_T > p_thresh_asym)[0] #find genes with prob(beta_i | g_i, T) > p_thresh_asym\n",
    "        asym_order = np.argsort(pasym_gi_T[asymgenes])[::-1]\n",
    "        asymgenes = asymgenes[asym_order[:min(N_thresh_asym,len(asymgenes))]]#keep top N_thresh_asym of these genes\n",
    "        ind_notlik = np.setdiff1d(range(3),tlik)\n",
    "        self.asymgenes,self.data,self.ind_notlik = asymgenes,data,ind_notlik\n",
    "        group = np.argmax(data[asymgenes,:][:,ind_notlik],axis=1)\n",
    "        status[asymgenes[group==0]] = 1 #mark status in output vector: +/- 1\n",
    "        status[asymgenes[group==1]] = -1 #mark status in output vector: +/- 1\n",
    "\n",
    "        #markers for root cell type\n",
    "        psym_gi_T = pmar_gi_T[:,tlik]\n",
    "        symgenes = np.where(psym_gi_T > p_thresh_sym)[0] #find genes with prob(alpha_i = 1| g_i, T) > p_thresh_sym\n",
    "        sym_order = np.argsort(psym_gi_T[symgenes])[::-1]\n",
    "        symgenes = symgenes[sym_order[:min(N_thresh_sym,len(symgenes))]]#keep top N_thresh_sym of these genes\n",
    "        status[symgenes] = 2; #mark status in output vector: +2\n",
    "\n",
    "        # markers for leaf cell types\n",
    "        daught = np.setdiff1d(range(3),tlik)\n",
    "        daught1 = daught[0] # code: +3\n",
    "        daught2 = daught[1] # code: -3\n",
    "        markergenes1 = np.where(pmar_gi_T[:,daught1] > p_thresh_sym)[0] #find genes with prob(alpha_i = 1| g_i, T) > p_thresh_sym\n",
    "        markergenes2 = np.where(pmar_gi_T[:,daught2] > p_thresh_sym)[0] #find genes with prob(alpha_i = 1| g_i, T) > p_thresh_sym\n",
    "        mar1_order = np.argsort(pmar_gi_T[markergenes1, daught1])[::-1]\n",
    "        mar2_order = np.argsort(pmar_gi_T[markergenes2, daught2])[::-1]\n",
    "        markergenes1 = markergenes1[mar1_order[:min(N_thresh_sym,len(markergenes1))]]#keep top N_thresh_sym of these genes\n",
    "        markergenes2 = markergenes2[mar2_order[:min(N_thresh_sym,len(markergenes2))]]#keep top N_thresh_sym of these genes\n",
    "        status[markergenes1] = 3;\n",
    "        status[markergenes2] = -3;\n",
    "\n",
    "        return status\n",
    "    def calculate_gene_probabilities(self,icomb,Integrals, Params,trip_data,tlik):\n",
    "        \"\"\"\n",
    "        Documentation:\n",
    "        #See Equations (27) and (28)\n",
    "        \"\"\"\n",
    "        odds0 = Params.odds0\n",
    "        psla = Params.psla\n",
    "        pgi_mas_T, pgi_sla, pgi_sym = self.calc_pgi_topo(icomb,Integrals, Params, trip_data,tag = 'calculate_gene_probabilities') #get probability of data given different hypotheses\n",
    "        # transition genes\n",
    "        odds_asym_topol = odds0*pgi_mas_T[:,tlik]/pgi_sla; # p(beta_i = 1 | g_i, T)/p(beta_i = 0 | g_i, T)\n",
    "        prob_asym_topol = odds_asym_topol/(1.+odds_asym_topol); # p(beta_i = 1 | g_i, T)\n",
    "        prob_asym_topol[np.isnan(prob_asym_topol)] = 0  #get rid of NaNs  \n",
    "\n",
    "        p_mas = odds0/(1.+odds0) # prior p(beta_i=1)\n",
    "        p_sla = 1./(1.+odds0) # prior p(beta_i=0)\n",
    "        pg_T = p_mas*pgi_mas_T[:,tlik] + p_sla*pgi_sla; # p(g_i | T)\n",
    "\n",
    "        #markers\n",
    "        #pg_T as above\n",
    "        prob_marker_topol = np.zeros(pgi_sym.shape)\n",
    "        for t in range(3):#1:3\n",
    "            psym = pgi_sym[:,t] # p(g_i | T, gene i marker for cell type t) i.e. p(g_i | T, beta_i = 0, alpha_i = 1)\n",
    "            pnotsym = pg_T-p_sla*psla[t+1]*psym; # p(g_i | T, gene i NOT marker for cell type t). psla(tlik+1) is p(mu_A > mu_BC | b_i = 0) (or equiv) \n",
    "            prior_sym_odds = (p_sla*psla[t+1])/(1-p_sla*psla[t+1]); # p(alpha_i = 1 for cell type t)/p(alpha_i = 0 for cell type t);\n",
    "            odds_marker_topol = prior_sym_odds*psym/pnotsym; # p(alpha_i = 1 | g_i, T)/p(alpha_i = 0 | g_i, T)\n",
    "            prob_marker_topol[:,t] = odds_marker_topol/(odds_marker_topol+1.); # p(alpha_i = 1 for cell type t | g_i, T)\n",
    "            #deal with nonsense\n",
    "            prob_marker_topol[odds_marker_topol == np.inf,t] = 1 ## changing because think it's providing spurious markers\n",
    "            prob_marker_topol[np.isnan(prob_marker_topol[:,t]),t] = 0 ## changing because think it's providing spurious markers\n",
    "        return prob_asym_topol, prob_marker_topol\n",
    "    def process_one_triplet(self,icomb, Integrals, Params, tripdata, makeplot=False, plot_dir=None, labs=['x','y','z']):\n",
    "        \"\"\"\n",
    "        Documentation:\n",
    "        Calls: calc_pgi_topo with tag process_one_triplet\n",
    "        \"\"\"\n",
    "        allodds = Params.oddsrange\n",
    "        topology_prior = Params.topology_prior\n",
    "        plikely_thresh = Params.plikely_thresh\n",
    "        #find p(g_i|T,b_i) for each gene starting from integrals\n",
    "        pgi_mas_T, pgi_sla, _ = self.calc_pgi_topo(icomb,Integrals, Params, tripdata,tag='process_one_triplet')\n",
    "\n",
    "        pT_g = np.zeros([len(allodds),4]) #pT_g as a function of prior odds\n",
    "        odds_mas_sla_T = pgi_mas_T/np.expand_dims(pgi_sla,-1) # See Equation (18)\n",
    "        #deal with nonsense\n",
    "        odds_mas_sla_T[np.isnan(odds_mas_sla_T)] = 1 #set nans to 1\n",
    "        odds_mas_sla_T[odds_mas_sla_T==np.inf] = 1 #set infs to 1\n",
    "        for i in range(len(allodds)):\n",
    "            prior_odds = allodds[i] #p(beta_i = 1)/p(beta_i = 0) Equation (21)\n",
    "            odds = 1+ prior_odds*odds_mas_sla_T\n",
    "            sum_log_odds = np.sum(np.log(odds),axis=0)\n",
    "            prododds = np.exp(sum_log_odds - np.max(sum_log_odds)) #divide by max(prod(odds)) to avoid Inf\n",
    "            pT_g[i,:] = prododds*topology_prior/np.sum(topology_prior*prododds); # this is Equation (20)\n",
    "        \n",
    "        self.odds_mas_sla_T = odds_mas_sla_T\n",
    "        self.pT_g = pT_g\n",
    "        \n",
    "        greaterthan0 = pT_g[:,:3] > plikely_thresh\n",
    "        ncandidates = np.sum(np.sum(greaterthan0,axis=0)>0); #how many topologies have probability greater than plikely_thresh\n",
    "        #ncandicates is a number in [0,1,2,3]\n",
    "\n",
    "        #find likely topology\n",
    "        tlikely_ind, plikely_ind, report_ind = self.find_Tlikely(pT_g, icomb)\n",
    "        self.report_ind = report_ind\n",
    "        self.tlikely_ind = tlikely_ind\n",
    "        # Report the best odds for a triplet\n",
    "        pT_g_max = np.concatenate([[allodds[report_ind]],pT_g[report_ind,:]])\n",
    "\n",
    "        if makeplot:\n",
    "            ###Not tested \n",
    "            #make a dot plot\n",
    "\n",
    "            pgi_mas_Tmin = np.column_stack([3*Integrals.IABC_Amin[:,icomb],3*Integrals.IABC_Bmin[:,icomb],3*Integrals.IABC_Cmin[:,icomb]])#find p(g_i | b_i=1, cell T is min)\n",
    "            imax = np.argmax(pgi_mas_Tmin,axis = -1); #find cell type that gene i votes against\n",
    "            f=plt.figure(figsize=(15,15))\n",
    "            plt.semilogx(odds_mas_sla_T,4-imax,'o')\n",
    "            #plot boundaries\n",
    "            plt.xlim([1e-10,10*np.max(allodds)])\n",
    "            plt.ylim([0,4])\n",
    "            #weird labels\n",
    "            plt.yticks(range(1,4), [labs[i_lab] for i_lab in [2,1,0]])\n",
    "            if plot_dir is not None:\n",
    "                f.savefig(plot_dir+os.sep+'topvotes-'+str(icomb)+'.png')\n",
    "                plt.close(f)\n",
    "            else:\n",
    "                plt.show()\n",
    "\n",
    "            #plot p(T|{g}) as a function of p(b_i=1)/p(b_i=0)\n",
    "            f = plt.figure(figsize=(15,15));\n",
    "            plt.semilogx(allodds, pT_g, '-'); \n",
    "            #legend(labs{1},labs{2},labs{3},'null', 'Location', 'northwest'); \n",
    "            plt.xlim([np.min(allodds),np.max(allodds)])\n",
    "            plt.ylim([-0.05,1.05])\n",
    "            if plot_dir is not None:\n",
    "                f.savefig(plot_dir+os.sep+'pT_g-'+str(icomb)+'.png')\n",
    "                plt.close(f)\n",
    "            else:\n",
    "                plt.show()\n",
    "        return tlikely_ind, plikely_ind, ncandidates, pT_g_max\n",
    "    def find_Tlikely(self,pT_g, triplet):\n",
    "        \"\"\"\n",
    "        Finds the most likely topology of triplet given pT_g(odds_prior,Topology)\n",
    "        tlikely_ind - most likely Topology 4-null \n",
    "        plikely_ind - maximum probability for that topology\n",
    "        report_ind - index of mostlikely odd\n",
    "        \"\"\"\n",
    "        m = np.max(pT_g,axis=-1)\n",
    "        itmax = np.argmax(pT_g,axis=-1) #find most likely topology for each of the odds\n",
    "\n",
    "        #this is not a fair treatment of topologies\n",
    "        ## if any of the non-nul topologies have a maximum for some odd then exclude nulltopology\n",
    "        itmaxnot4 = itmax!=3#bool not null topology across odds\n",
    "        if np.sum(itmaxnot4)!=0:\n",
    "             m[itmaxnot4==False]=-np.inf # set all null topology values to - infinity\n",
    "\n",
    "        report_ind = np.argmax(m)\n",
    "        tlikely_ind  = itmax[report_ind] #ensures the convention null topology is 4\n",
    "        #tlikely_ind  = (itmax[report_ind]+1)%4 #ensures the convention null topology is 0\n",
    "        plikely_ind = np.max(m)\n",
    "\n",
    "        return tlikely_ind, plikely_ind, report_ind\n",
    "    def calc_pgi_topo(self,icomb,Integrals, Params, data,tag=''):\n",
    "        \"\"\"\n",
    "        Documentation...\n",
    "        This appears in 2 files: , identical forms but the filter_min condition is different.\n",
    "        process_one_triplet.m,calculate_gene_probabilities.m\n",
    "        Currently the process_one_triplet.m is implemented\n",
    "        \"\"\"\n",
    "        log_asym_min = Params.log_asym_min\n",
    "        log_sym_min = Params.log_sym_min\n",
    "        filter_min = Params.filter_min\n",
    "        psla = Params.psla\n",
    "\n",
    "        pgi_sla_1 = Integrals.IABCsame[:,icomb] #p(gi | alpha_i = 0, beta_i = 0) (Equation 9)\n",
    "        pgi_sla_2Amax = 2*Integrals.IA_BC_Amax[:,icomb] \n",
    "        pgi_sla_2Amax[pgi_sla_2Amax<0.] = 0. #(this is unlikely: should all be > 0 but I did find one floating point error)\n",
    "        pgi_sla_2Bmax = 2*Integrals.IA_BC_Bmax[:,icomb]\n",
    "        pgi_sla_2Bmax[pgi_sla_2Bmax<0.] = 0.\n",
    "        pgi_sla_2Cmax = 2*Integrals.IA_BC_Cmax[:,icomb] \n",
    "        pgi_sla_2Cmax[pgi_sla_2Cmax<0.] = 0. #p(gi | alpha_i = 1, beta_i = 0) Equation (8)\n",
    "        \n",
    "        #internalize for troubleshooting:\n",
    "        self.pgi_sla_1 = pgi_sla_1\n",
    "        \n",
    "        \n",
    "        pgi_sym = np.column_stack([pgi_sla_2Amax,pgi_sla_2Bmax,pgi_sla_2Cmax])\n",
    "\n",
    "        pgi_mas_Tmin = 3.*np.column_stack([Integrals.IABC_Amin[:,icomb],Integrals.IABC_Bmin[:,icomb],Integrals.IABC_Cmin[:,icomb]]) #p(gi | beta_i = 1, cell type T is a leaf)\n",
    "        pgi_mas_T = (1./2.)*np.column_stack([pgi_mas_Tmin[:,1]+pgi_mas_Tmin[:,2],\n",
    "                                             pgi_mas_Tmin[:,0]+pgi_mas_Tmin[:,2],\n",
    "                                             pgi_mas_Tmin[:,1]+pgi_mas_Tmin[:,0]]); #p(gi | beta_i = 1, topology T) Equation (4)\n",
    "        pgi_sla = psla[0]*pgi_sla_1 + psla[1]*pgi_sla_2Amax + psla[2]*pgi_sla_2Bmax + psla[3]*pgi_sla_2Cmax; # p(g_i | beta_i = 0) See Equation (14)\n",
    "\n",
    "        if filter_min: #some additional coarse filtering if option chosen\n",
    "            expressed = data>log_asym_min\n",
    "            sym_expressed = data>log_sym_min\n",
    "            pgi_mas_T = pgi_mas_T*expressed*((np.expand_dims(np.sum(expressed,axis=1)>1,-1))*np.ones([1,3])) #transition genes above expression threshold in at least 2 cell types\n",
    "            pgi_sym = pgi_sym*sym_expressed #marker genes above expression threshold \"marked\" cell type\n",
    "            if tag=='calculate_gene_probabilities':\n",
    "                imax = np.argmax(data,axis=-1)\n",
    "                for t in range(3):\n",
    "                    pgi_sym[:,t] = pgi_sym[:,t]*(imax==t); #marker gene most highly expressed in \"marked\" cell type\n",
    "            \n",
    "        pgi_mas_0 = (1./3.)*np.sum(pgi_mas_T,axis=1) #p(g_i | b_i = 1, T = null) Equation (6), Equation (7)\n",
    "        pgi_mas_T = np.column_stack([pgi_mas_T,pgi_mas_0]) #p(g_i | b_i = 1, T)\n",
    "        self.pgi_mas_T = pgi_mas_T\n",
    "        self.Integrals = Integrals\n",
    "        return pgi_mas_T, pgi_sla, pgi_sym\n",
    "    def calculate_integrals_one_gene(self,g_dist, i, ncells, combinations, log2tfdata, \n",
    "                                 iunique, idxs, P, M, S, logS, lM, lS, firstcut, dm, ds,normed=True):\n",
    "        ncombinations = len(combinations)\n",
    "        tempIA = np.zeros([1,ncells])\n",
    "        tempIAB = np.zeros([ncells,ncells,1])\n",
    "        tempIABCsame = np.zeros([1, ncombinations])\n",
    "        tempIA_BC_Amax = np.zeros([1, ncombinations])\n",
    "        tempIA_BC_Bmax = np.zeros([1, ncombinations])\n",
    "        tempIA_BC_Cmax = np.zeros([1, ncombinations])\n",
    "        tempIABC_Amin = np.zeros([1, ncombinations])\n",
    "        tempIABC_Bmin = np.zeros([1, ncombinations])\n",
    "        tempIABC_Cmin = np.zeros([1, ncombinations])\n",
    "        Isigma = np.zeros(ncells,dtype=object); #pre-allocate vector for integration over sigma for one cell type\n",
    "        Isigma2 = np.zeros([ncells,ncells],dtype=object); #pre-allocate matrix for integration over sigma for two cell types\n",
    "        diags = np.zeros([ncells,ncells],dtype=object); #pre-allocate cell for integration over mu for two cell types\n",
    "        data = log2tfdata[i,:]; #vector of gene(i) expression levels in all cells\n",
    "\n",
    "        best_mu_gene = np.zeros(ncells); #pre-allocate vector for the parameters that give the highest likely hood for each gene in each cell type\n",
    "        best_sigma_gene = np.zeros(ncells);\n",
    "        for j in range(ncells): #for each cell cluster j\n",
    "            #N1 = sum(iunique==idxs(j));#integrals over one cell type (to use in integrals for Equation 9)\n",
    "            xk1 = data[iunique==idxs[j]] #gene(i) expression levels in all cells in cluster j\n",
    "\n",
    "            F1 = self.logDa(g_dist,M,S,xk1,logS); #returns a sbin by mbin matrix of logDa(gi|mu,sigma) for gene(i) in cluster j. See equation (11) (or the equivalent for ZIP).\n",
    "            # Save most likely distribution for each gene, per cluster.\n",
    "            i1,i2 = np.unravel_index(np.argmax(F1),F1.shape)\n",
    "            best_mu_gene[j],best_sigma_gene[j]= M[0,i2],S[i1,0]\n",
    "            #if normed: F1-=np.max(F1) ##normalize\n",
    "            F1 = np.exp(F1)\n",
    "            F1 = P*F1  # multiply by prior, but flat.\n",
    "            F1[np.isnan(F1)] = 0 # correct for NaN?\n",
    "\n",
    "            self.F1=F1#internalize for troubleshooting\n",
    "\n",
    "            tempIA[0,j] = dm*ds*np.trapz(np.trapz(F1,axis=0),axis=0); #integrate Da over all prior mu and sigma for gene(i) in cluster j. Can be used in Eq.7 for T=0.\n",
    "            #if normed: F1 = F1/tempIA[0,j]\n",
    "            \n",
    "            Isigma[j] = ds*np.trapz(F1,axis=0);#store integrals over sigma for use later. Isigma{j} is a vector of length mbins.\n",
    "            for k in range(j+1,ncells):#integrals over 2 cell types (to use in integrals for Equation 9)\n",
    "                #N2 = sum(iunique==idxs(k)); N12 = N1+N2;\n",
    "                xk2 = data[iunique==idxs[k]]#(ismember(iunique,idxs(k)));\n",
    "                xk12 = np.concatenate([xk1,xk2])\n",
    "                F2 = self.logDa(g_dist,M,S,xk12,logS)\n",
    "                #if normed: F2-=np.max(F2) ##normalize\n",
    "                F2 = np.exp(F2)\n",
    "                F2 = P*F2  # multiply by prior, but flat.\n",
    "                F2[np.isnan(F2)] = 0 # correct for NaN?\n",
    "                tempIAB[j,k,0] = dm*ds*np.trapz(np.trapz(F2,axis=0),axis=0); #integrate Da over the hyperprior space for gene(i) in cluster j+k.\n",
    "                tempIAB[k,j,0] = tempIAB[j,k,0]; #tempIAB is a symmetric 3by3 matrix.\n",
    "                #if normed: F2 = F2/tempIAB[j,k,0]\n",
    "                Isigma2[j,k] = ds*np.trapz(F2,axis=0); #Isigma2 is a symmetric 3by3 cell.\n",
    "                Isigma2[k,j] = Isigma2[j,k];\n",
    "        for j in range(ncells):\n",
    "            for k in range(j+1,ncells):\n",
    "                Isigmaj,Isigmak = np.meshgrid(Isigma[j],Isigma[k])\n",
    "                AB = Isigmaj*Isigmak\n",
    "                #AB = repmat(reshape(Isigma{j},[lM 1]),[1 lM]).*repmat(reshape(Isigma{k},[1 lM]),[lM 1]);\n",
    "                AB_diag = np.diag(cumtrapz(cumtrapz(AB[::-1,::-1],axis=0),axis=1))[::-1]#AB_diag is a vector of length mbins.\n",
    "                #Each element is the integral of [Integral(Dj over sigma)*Integral(Dk over sigma)]\n",
    "                #over mu values greater than a certain mu (in the prior space).\n",
    "                diagjk = np.array(AB_diag[:])#delete the first element of integration of Dj and Dk over the complete hyperprior space.\n",
    "                                             #done internally by cumtrapz in python\n",
    "                prev_firstcut = 0 if firstcut==0 else firstcut-1\n",
    "                diagjk[:prev_firstcut] = diagjk[firstcut] #replace integrations over mu>mu_replace for mu_replace<mu_cutoff with the integrations over mu>mu_cutoff\n",
    "                diags[j,k] = diagjk #diags is a 3by3 symmetric cell with each element being a (mbins-1)by1 vector\n",
    "                diags[k,j] = diags[j,k]\n",
    "\n",
    "        for j in range(ncombinations):\n",
    "            ind123 = np.sum([iunique==idx_ for idx_ in idxs[combinations[j,:]]],axis=0)>0\n",
    "            xk123 = data[ind123] #computation of integrals over 3 cell types assuming they are from a single distribution\n",
    "            #N123 = length(xk123);\n",
    "            F3 = self.logDa(g_dist, M,S,xk123,logS)\n",
    "            #if normed: F3-=np.max(F3) ##normalize\n",
    "            F3 = np.exp(F3)\n",
    "            F3 = P*F3  # multiply by prior, but flat.\n",
    "            #if normed: F3 = F3/np.sum(F3) ##normalize\n",
    "            F3[np.isnan(F3)] = 0 # correct for NaN?\n",
    "            tempIABCsame[0,j] = dm*ds*np.trapz(np.trapz(F3,axis=0),axis=0); # Equation (10). Integration of Da over the entire hyperprior space\n",
    "            \n",
    "            #troubleshoot\n",
    "            self.Isigma=Isigma\n",
    "            self.diags = diags\n",
    "            self.Isigma_t = Isigma[combinations[j,0]][:-1]\n",
    "            self.diags_t = diags[combinations[j,1],combinations[j,2]]\n",
    "            self.combinations = combinations\n",
    "            self.j = j\n",
    "            \n",
    "            \n",
    "            tempIABC_Amin[0,j] = dm*dm*dm*np.trapz(Isigma[combinations[j,0]][:-1]*diags[combinations[j,1],combinations[j,2]]); # Equation (6)\n",
    "            tempIABC_Bmin[0,j] = dm*dm*dm*np.trapz(Isigma[combinations[j,1]][:-1]*diags[combinations[j,0],combinations[j,2]])# if cutoff=0, this is exactly Equation(6)\n",
    "            tempIABC_Cmin[0,j] = dm*dm*dm*np.trapz(Isigma[combinations[j,2]][:-1]*diags[combinations[j,0],combinations[j,1]])\n",
    "\n",
    "            IMesh1,IMesh2 = np.meshgrid(Isigma[combinations[j,0]],Isigma2[combinations[j,1],combinations[j,2]])\n",
    "            A_BC = IMesh1*IMesh2\n",
    "            IMesh1,IMesh2 = np.meshgrid(Isigma[combinations[j,1]],Isigma2[combinations[j,0],combinations[j,2]])\n",
    "            B_AC = IMesh1*IMesh2\n",
    "            IMesh1,IMesh2 = np.meshgrid(Isigma[combinations[j,2]],Isigma2[combinations[j,0],combinations[j,1]])\n",
    "            C_AB = IMesh1*IMesh2\n",
    "            A_BC = np.tril(A_BC,0);\n",
    "            B_AC = np.tril(B_AC,0);\n",
    "            C_AB = np.tril(C_AB,0);\n",
    "            tempIA_BC_Amax[0,j] = dm*dm*np.trapz(np.trapz(A_BC,axis=0),axis=0); # Equation (9) but why tril(X,1) instead of tril(X,0)??\n",
    "            tempIA_BC_Bmax[0,j] = dm*dm*np.trapz(np.trapz(B_AC,axis=0),axis=0);\n",
    "            tempIA_BC_Cmax[0,j] = dm*dm*np.trapz(np.trapz(C_AB,axis=0),axis=0);\n",
    "        return [tempIA, tempIAB, tempIABCsame, tempIABC_Amin, tempIABC_Bmin, tempIABC_Cmin, tempIA_BC_Amax, tempIA_BC_Bmax, tempIA_BC_Cmax, best_mu_gene, best_sigma_gene]\n",
    "    def calc_log_mean_std(self,data,iunique, g_dist):\n",
    "        \"\"\"Given:\n",
    "        # data # an expression array m x n\n",
    "        # iunique # a vector of cell identies of length n\n",
    "        # g_dist # a string id for the distribution N - (log)normal distribution, P - ZeroInflatedPoisson (ZIP)\n",
    "        Returns: \n",
    "        means and stds IF g_dist = N\n",
    "        max_likelihood of means of ZIP and zero probs IF g_dist = P \"\"\"\n",
    "        #cast as numpy arrays for easy computations\n",
    "        data = np.array(data)\n",
    "        iunique = np.array(iunique)\n",
    "\n",
    "        cell_types = np.unique(iunique)\n",
    "        sx,sy = data.shape\n",
    "        loggenemeans = np.zeros([sx, len(cell_types)])\n",
    "        loggenestds = np.zeros([sx, len(cell_types)])\n",
    "\n",
    "        # iterate through cell types\n",
    "        for i,cell_type in enumerate(cell_types):\n",
    "            datai= data[:,iunique==cell_type]#data of cell type\n",
    "            splmean = np.mean(datai,axis=-1)\n",
    "            splstd = np.std(datai,axis=-1)\n",
    "            if g_dist == 'N':\n",
    "                loggenemeans[:,i] = splmean\n",
    "                loggenestds[:,i] = splstd\n",
    "            elif g_dist == 'P':\n",
    "                #maximum likelihood \n",
    "                # for lambda we want to solve: xbar (1-exp(-lambda)) = lambda (1-n0/n)\n",
    "                # https://en.wikipedia.org/wiki/Zero-inflated_model\n",
    "                # the analytical solution to 1-exp(-x)= b x is: (1+b*lambertw(-np.exp(-1./b)/b))/b\n",
    "                n0 =np.sum(datai==0,axis=-1)\n",
    "                n = float(datai.shape[-1])\n",
    "                b = (1.-n0/float(n))/splmean\n",
    "                max_likelihood_lambda = ((1.+b*lambertw(-np.exp(-1./b)/b))/b).real\n",
    "                max_likelihood_pi = 1.-splmean/max_likelihood_lambda\n",
    "                non_sense_1 = np.isnan(max_likelihood_lambda)&(max_likelihood_pi>1)\n",
    "                non_sense_2 = max_likelihood_pi<0\n",
    "                #deal with nonsense\n",
    "                max_likelihood_lambda[non_sense_1]=splmean[non_sense_1]\n",
    "                max_likelihood_pi[non_sense_1]=1.\n",
    "                max_likelihood_lambda[non_sense_2]=splmean[non_sense_2]\n",
    "                max_likelihood_pi[non_sense_2]=0.\n",
    "                loggenemeans[:,i] = max_likelihood_lambda\n",
    "                loggenestds[:,i] = max_likelihood_pi\n",
    "        return loggenemeans, loggenestds  \n",
    "    def simulate_ZIP(self,mu_true=12,psi_true=0.5,nx=7,ny=2):\n",
    "        \"\"\"Simulates zero inflated poisson\n",
    "        with mean mu_true, and sigma = 1-psi_true (prob for 0) nx - number of sets, ny is number of observations/set\n",
    "        \"\"\"\n",
    "        # True parameter values\n",
    "        #mu_true,psi_true,nx,ny = 12,0.5,7,2\n",
    "        # Simulate some data\n",
    "        datai = np.array([[pm.rpoisson(mu_true)*(np.random.random()<psi_true) for i in range(ny)] for i_ in range(nx)])\n",
    "        return datai\n",
    "    def test_likelihood(self):\n",
    "        \"\"\"\n",
    "        Test logDa function using simulate_ZIP function \n",
    "        \"\"\"\n",
    "        datai = self.simulate_ZIP(mu_true=35, psi_true=0.5, nx=1, ny=30)\n",
    "        xk1 = datai[0]\n",
    "        m = np.linspace(30,40,100)\n",
    "        s = np.linspace(0.1,0.9,98)\n",
    "        #s = np.linspace(3,800,100)\n",
    "        M,S = np.meshgrid(m,s)\n",
    "        log_vals = self.logDa('P',M,S,xk1,logS=None)\n",
    "        i1,i2 = np.unravel_index(np.argmax(log_vals),log_vals.shape)\n",
    "        \n",
    "        import matplotlib.cm as cm\n",
    "        import matplotlib.pylab as plt\n",
    "        plt.figure()\n",
    "        plt.title('Max likelihood vals:'+str([m[i2],s[i1]]))\n",
    "        plt.imshow(np.exp(log_vals),interpolation='nearest')\n",
    "        plt.show()\n",
    "        plt.figure()\n",
    "        plt.hist(xk1)\n",
    "        plt.show()\n",
    "    def get_parameters(self):\n",
    "        \"\"\"\n",
    "        Sets default paramaters\n",
    "        \"\"\"\n",
    "        Params = AttributeDict()\n",
    "        #numeric integration parameters\n",
    "        Params.use_kde = 0; # do NOT turn this on until function is fixed to deal with n=[sbins,mbins]. 1 -> use kde-estimated empirical prior;  0 -> use uniform prior\n",
    "        Params.use_parallel_integrals = 1; #compute integrals in parallel?\n",
    "        Params.use_parallel_topology = 0;\n",
    "        Params.use_parallel_genes = 0;\n",
    "        Params.use_parallel = any([Params.use_parallel_integrals, Params.use_parallel_topology, Params.use_parallel_genes]);\n",
    "        Params.parallel_pool_size = 16; #number of processors for parallelization\n",
    "\n",
    "        Params.dist = 'N'; #distribution for gene expression. Default is 'N' for normal. Can also use 'P' for ZIP.\n",
    "        Params.mumin = 0; # min and max for mean prior in normal distribution, or for lambda prior in ZIP. \n",
    "        Params.mumax = 7;\n",
    "        Params.sigmin = 0; # min and max for std prior in normal distribution, or for sigma prior in ZIP.\n",
    "        Params.sigmax = 3;\n",
    "        Params.mbins = 2**8; #number of bins for numeric integration over distribution parameters\n",
    "        Params.sbins = 2**8;\n",
    "\n",
    "        Params.cutoff = np.log2(10); # expression threshold for numeric integration\n",
    "\n",
    "        Params.verbose = 2; # 0 = default (no) reporting\n",
    "                            # 1 = extra results\n",
    "                            # 2 = debugging\n",
    "        Params.save_integrals = 1; # save integrals for later dot-plotting\n",
    "        Params.save_gene_probs = 1; # save gene probabilities for later\n",
    "        Params.save_best_dists = 1; # save info about 'best' mu and sigma per gene/cluster\n",
    "\n",
    "        Params.psla = [1./2, 1./6, 1./6, 1./6]; #vector of prior probabilities for non-transition genes\n",
    "        Params.topology_prior = [1./4, 1./4, 1./4, 1./4]; #topology prior [p(A) p(B) p(C) p(0)]\n",
    "        Params.plikely_thresh = 0.6; # threshold for probability of topology given data\n",
    "        Params.filter_min = 0; #ignore genes with values under threshold?\n",
    "        Params.log_asym_min = np.log2(2); #threshold for transition gene expression\n",
    "        Params.log_sym_min = np.log2(2); #threshold for marker gene expression\n",
    "        Params.oddsrange = np.logspace(-4,2,50); #range to vary prior odds p(beta_i = 1)/p(beta_i = 0)\n",
    "\n",
    "        #parameters for gene probabilities\n",
    "        Params.N_thresh_sym = 1000; # max number of genes in gene classes\n",
    "        Params.N_thresh_asym = 1000; \n",
    "        Params.p_thresh_asym = 0.8; # probability thresholds for defining gene classes\n",
    "        Params.p_thresh_sym = 0.8; \n",
    "        Params.odds0 = 0.05; # prior transition \n",
    "        #do checks\n",
    "        if Params.dist == 'P':\n",
    "            if Params.sigmax > 1 or Params.sigmax < 0:\n",
    "                raise ValueError('Params.sigmax for Zero Inflated Poisson distribution must be in the range [0,1]')\n",
    "        self.Params=Params\n",
    "        return Params\n",
    "    \n",
    "    def calculate_integrals(self,log2tfdata, idxs, iunique, Params):\n",
    "        \"\"\"Given:\n",
    "        #log2tfdata# expresion matrix of shape m x n\n",
    "        #idxs# a vector with a subset of values from iunique ...\n",
    "        #iunique# a \n",
    "        This returns:\n",
    "        [Integrals, best_mu, best_sigma,combinations]\n",
    "        The integrals used in this paper, the best_mu and best_sigma given data and the combinations matrix used.\n",
    "        \"\"\"\n",
    "        #cast as numpy arrays for easy computation\n",
    "        log2tfdata = np.array(log2tfdata)\n",
    "        idxs = np.array(idxs)\n",
    "        iunique = np.array(iunique)\n",
    "        #internalize paramaters\n",
    "        nbins = [Params.sbins, Params.mbins]\n",
    "        g_dist = Params.dist\n",
    "        mumin = Params.mumin\n",
    "        mumax = Params.mumax\n",
    "        sigmin = Params.sigmin \n",
    "        sigmax = Params.sigmax\n",
    "        cutoff = Params.cutoff\n",
    "\n",
    "        use_parallel_integrals = Params.use_parallel_integrals\n",
    "        #compute statistics\n",
    "        [loggenemeans, loggenestds] = self.calc_log_mean_std(log2tfdata, iunique, g_dist)\n",
    "\n",
    "        #uniform prior over means and standard deviations #kde not implemented\n",
    "        dm = (mumax-mumin)/nbins[1]\n",
    "        ds = (sigmax-sigmin)/nbins[0] # will be zero is sigmax==sigmin\n",
    "        if sigmax==sigmin: ds=1\n",
    "        M,S=np.meshgrid(np.arange(mumin,mumax,dm),np.arange(sigmin,sigmax,ds))\n",
    "        density = np.ones(M.shape)\n",
    "        P = density/(dm*ds*np.sum(density))\n",
    "\n",
    "\n",
    "        ################ Compute integrals\n",
    "\n",
    "        #initialize paramaters\n",
    "        ngenes = log2tfdata.shape[0]\n",
    "        ncells = len(idxs)\n",
    "        from itertools import combinations as combnk\n",
    "        combinations = np.array(list(combnk(range(ncells),3)))\n",
    "        ncombinations = len(combinations)\n",
    "        lM,lS = M.shape\n",
    "        Mu = M[0,1:]###Mu = M(1,2:end); what? stupid?\n",
    "        firstcut = np.sum(Mu<=cutoff)###\n",
    "        ## Compute paramters over integration space required for likelihood of data given hyperprior\n",
    "        if g_dist=='N':\n",
    "            logS = np.log(2*np.pi*S**2);\n",
    "        elif g_dist=='P':\n",
    "            logS1 = np.log(1-S);\n",
    "            ExpM = np.exp(-M); #exp(-mu). sbins by mbins matrix.\n",
    "            logSM = np.log(S+(1-S)*ExpM); #log[sigma+(1-sigma)exp(-mu)]. sbins by mbins matrix.\n",
    "            logS = [logS1,logSM];\n",
    "        else:\n",
    "            raise ValueError('Invalid distribution parameter. Acceptable input: \"N\" or \"P\"')\n",
    "\n",
    "        #initialize to 0 target variables\n",
    "\n",
    "        IA = np.zeros([ngenes, ncells]);#integrals over one cell type (to use in integrals for Equation 8)\n",
    "        IAB = np.zeros([ncells,ncells,ngenes]);#integrals over 2 cell types (to use in integrals for Equation 8)\n",
    "        IABCsame = np.zeros([ngenes, ncombinations]);#integrals over 3 cell types (Equation 9) muA=muB=muC\n",
    "        IA_BC_Amax = np.zeros([ngenes, ncombinations]);#integrals with muA < muB=muC (Equation 8)\n",
    "        IA_BC_Bmax = np.zeros([ngenes, ncombinations]);#integrals with muB < muA=muC (Equation 8)\n",
    "        IA_BC_Cmax = np.zeros([ngenes, ncombinations]);#integrals with muC < muB=muA (Equation 8)\n",
    "        IABC_Amin = np.zeros([ngenes, ncombinations]);#integrals with muA < muB and muA < muC (Equation 5)\n",
    "        IABC_Bmin = np.zeros([ngenes, ncombinations]);#integrals with muB minimum (Equation 5)\n",
    "        IABC_Cmin = np.zeros([ngenes, ncombinations]);#integrals with muC minimum (Equation 5)\n",
    "\n",
    "        best_mu = np.zeros([ngenes, ncells]);\n",
    "        best_sigma = np.zeros([ngenes, ncells]);\n",
    "\n",
    "        #itterate through each cell type in idx\n",
    "        for i in range(ngenes):\n",
    "            [tempIA, tempIAB, tempIABCsame, \n",
    "             tempIABC_Amin, tempIABC_Bmin, tempIABC_Cmin, \n",
    "             tempIA_BC_Amax, tempIA_BC_Bmax, tempIA_BC_Cmax, \n",
    "             best_mu_gene, best_sigma_gene] = self.calculate_integrals_one_gene(g_dist, i, ncells, \n",
    "                                                                           combinations, log2tfdata, iunique, idxs, \n",
    "                                                                           P, M, S, logS, lM, lS, firstcut, dm, ds)\n",
    "            IA[i,:] = tempIA\n",
    "            IAB[:,:,i] = np.squeeze(tempIAB)\n",
    "            IABCsame[i,:] = tempIABCsame\n",
    "            IABC_Amin[i,:] = tempIABC_Amin\n",
    "            IABC_Bmin[i,:] = tempIABC_Bmin\n",
    "            IABC_Cmin[i,:] = tempIABC_Cmin\n",
    "            IA_BC_Amax[i,:] = tempIA_BC_Amax\n",
    "            IA_BC_Bmax[i,:] = tempIA_BC_Bmax\n",
    "            IA_BC_Cmax[i,:] = tempIA_BC_Cmax\n",
    "            best_mu[i,:] = best_mu_gene\n",
    "            best_sigma[i,:] = best_sigma_gene\n",
    "            if i%100==0:\n",
    "                print str(i)+' genes integrated.'\n",
    "        #Save to Integrals object\n",
    "        Integrals = AttributeDict()\n",
    "        Integrals.IA = IA;#not used\n",
    "        Integrals.IAB = IAB;#not used\n",
    "        Integrals.IABCsame = IABCsame;\n",
    "        Integrals.IABC_Amin = IABC_Amin;\n",
    "        Integrals.IABC_Bmin = IABC_Bmin;\n",
    "        Integrals.IABC_Cmin = IABC_Cmin;\n",
    "        Integrals.IA_BC_Amax = IA_BC_Amax;\n",
    "        Integrals.IA_BC_Bmax = IA_BC_Bmax;\n",
    "        Integrals.IA_BC_Cmax = IA_BC_Cmax;\n",
    "        return [Integrals, best_mu, best_sigma,combinations]\n",
    "\n",
    "    def logDa(self,g_dist,M,S,xk1,logS=None):\n",
    "        \"\"\"\n",
    "        ##This function returns a sbins by mbins matrix of logDa(gi|mu,sigma) for gene(i) in cluster j.\n",
    "        #input arguments:\n",
    "        # g_dist: Assumption for gene expression.'N' for normal, or 'P' for ZIP (zero inflated poisson).\n",
    "        # M,S: meshgrid components (each is a sbin by mbin matrix) created with the prior mu and sigma values.\n",
    "        # xk1: gene(i) expression levels in all cells in cluster j\n",
    "        # logS: Terms independent of gi in logDa: Equals log[2pi(sigma)^2] for\n",
    "        # 'N', and {log(1-sigma), log[sigma+(1-sigma)exp(-mu)]} for 'P'.\n",
    "        # lM, lS: number of prior mu and sigma values. Should be roughly mbins, sbins\n",
    "        Don't understand why are you passing these. They should be exactly mbins and sbins\n",
    "        \"\"\"\n",
    "        xk1 = np.array(xk1)\n",
    "        M = np.array(M)\n",
    "        S = np.array(S)\n",
    "        if logS is None:\n",
    "            if g_dist=='N':\n",
    "                logS = np.log(2*np.pi*S**2);\n",
    "            elif g_dist=='P':\n",
    "                logS1 = np.log(1-S);\n",
    "                ExpM = np.exp(-M); #exp(-mu). sbins by mbins matrix.\n",
    "                logSM = np.log(S+(1-S)*ExpM); #log[sigma+(1-sigma)exp(-mu)]. sbins by mbins matrix.\n",
    "                logS = [logS1,logSM];\n",
    "\n",
    "        N1 = len(xk1); #number of cells in cluster j\n",
    "        if g_dist == 'N':\n",
    "            F1 = -N1*logS/2; #-(N/2)log[2pi(sigma)^2], sbins by mbins matrix, same value for each row.\n",
    "            MM,XK = np.meshgrid(M[0,:],xk1); #MM, XK are Ncells by mbins matrices. MM contains prior mu values column-wise. XK contains gene(i) levels row-wise.\n",
    "            F1a = np.tile(np.sum(-(MM-XK)**2,axis=0),[S.shape[0],1]); #sum[-(gi-mu)^2]. sbins by mbins matrix. Identical rows.\n",
    "            F1 = F1 + F1a/(2*S*S); #logDa = -(N/2)log[2pi(sigma)^2]+sum[-(gi-mu)^2]/(2sigma^2). sbins by mbins.\n",
    "            #F1 = F1a/(2*S*S); #logDa = -(N/2)log[2pi(sigma)^2]+sum[-(gi-mu)^2]/(2sigma^2). sbins by mbins.\n",
    "        elif g_dist == 'P':\n",
    "            logS1 = logS[0]; #log(1-sigma). sbins by mbins matrix. Identical rows.\n",
    "            logSM = logS[1]; #log[sigma+(1-sigma)exp(-mu)]. sbins by mbins matrix.\n",
    "            xk1p = xk1[xk1>0]; #gene(i) levels in all cells expression gene(i) in cluster j. Vector.\n",
    "            Np = len(xk1p); #number of cells with positive gene(i) expression.\n",
    "            N0 = N1-Np; #number of cells with zero gene(i) level.\n",
    "            sum_g = np.sum(xk1p); #single value\n",
    "            sum_gamma = np.sum(gammaln(xk1p+1)); #for gi>0, sum{log[gamma(gi+1)]}. Single value.\n",
    "            sum_Mg = np.log(M[0,:])*sum_g; #vector of length mbins. log(mu)*Sum(gi) for gi>0.\n",
    "            F1 = Np*logS1+N0*logSM-Np*M; #Np*log(1-sigma)+N0*log[sigma+(1-sigma)exp(-mu)]-Np*mu. sbins by mbins matrix.\n",
    "            F1 = F1 + np.expand_dims(sum_Mg, 0); #F1 as above + log[mu*Sum(gi) for gi>0]. sbins by mbins matrix.\n",
    "            F1 = F1 - sum_gamma; #logDa = Np*log(1-sigma)+N0*log[sigma+(1-sigma)exp(-mu)]-Np*mu+log[mu*Sum(gi)]-sum{log[gamma(gi+1)]}, for gi>0.\n",
    "        return F1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sibilant started.\n",
      "Beginning integration.\n",
      "0 genes integrated.\n",
      "Starting trajectories.\n",
      "Starting gene selection.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\bogda\\Anaconda2\\lib\\site-packages\\ipykernel\\__main__.py:703: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAz4AAAEOCAYAAABILsJcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XmYJFWVsPH3dDcCDdIouxuKCC4g0OCK+y7ihoiCoyDu\n4oyD+onrAI7juA0wwuCIC+0GfrhvCIKCGyKftAuIgsim0qBNY7M0Db2c74/IgqjsrKpcbldmRb2/\n54mnIm/eOHEyMioqT0XkjchMJEmSJKnJ5gw7AUmSJEla3yx8JEmSJDWehY8kSZKkxrPwkSRJktR4\nFj6SJEmSGs/CR5IkSVLjWfhIkiRJajwLH0mSJEmNZ+EjSZIkqfEsfCRJIyEizo2Ic4adx6iJiKMi\nYu2w85Ckmc7CR5JmgIg4OCLWTjCtiYhHDjvHAhJo9Af8iHhnRDy/x8Uav10kaTpEZg47B0nSFCLi\nYOAzwHuBqzp0OSMzl01rUoVFxDyAzFw97FzWl4i4GfhyZh7awzJzgHmZecf6y0ySmm/esBOQJPXk\njMxcPOwkImLjzLytZMwmFzz9iIj5mbkiM9cCFj2SNCAvdZOkBml9H2RNRDy5rf2kiLg9InZtPX5i\n6zK5AyLiAxGxJCJuiYhvRsR92pY9NyJ+GxELI+LHEXEr8B+155/dar8lIm6KiO9ExEPbYmwTESdH\nxJ8jYmVEXBsR34iI+7Wt54dty20VEZ+OiOsi4raI+HVEvKKtz/at1/KWiHhNRFzeWscFEbFXF9ts\n7DLCvSPiYxHxt4i4MSL+NyLmRcSCiPhcRCxrTR/qEONtEfGziFgaESsi4pcR8aK2PmuB+cAhtcsU\nP1N739ZGxEMi4pSIWAb8pP5cLc7Y8oe0xX9Xq/1ZU71mSZqNPOMjSTPLgojYoq0ta5e5vR94LvDp\niNg1M2+NiGcCrwbenZkXtS37bqrvj3wQ2Bo4HDgrInbPzNvH4gNbAqcDXwI+B1wPEBEvBxYBZwBv\np/pg/wbgJxGxR2Ze04rxNeAhwMeAq1vrejpwP+Ca2nruFBEbAT8CdgCOp7rE78XAoohYkJnHt72W\nlwGbAv/binUE8NWI2CEz13TYlu2OB5YA/wY8GngN8A/gsa2c3wnsA7wtIi7KzC/Ulv0X4JvAF4C7\nAS8FTouIfTPze60+/wR8GvgFcFKr7U9tr/3LwGWtdUXtuTu3TWYuioj9gGMi4qzM/GuroP034JOZ\neUYXr1WSZp/MdHJycnIa8Qk4mKpA6TStaOv7MGAl8AlgAfAX4HxgTq3PE1vLXgPMr7Xv32p/U63t\nHGAN8Oq29WwCLAM+3ta+FXAj8L+txwtaMd8yxWs8B/hh7fGbW+t9aa1tLvAzYDmwSatt+1b8vwGb\n1fo+t7X8Pl1u2++2tf+stfwJtbY5rW32w7a+G7Y9ngv8Fjirrf1m4DMdcjiylcPnJ3huTVvbNsBS\nqoJzA2AxcAWw6bD3VScnJ6dRnbzUTZJmjqQ6m/K0tunZ4zpl/o7qw/JrgDOBewIHZ/VdkXafzcwV\ntWW/QnXWY5+2frdTndmpezpVUfOliNhibGrl+Qtg7HK726i+o/KkiNi8h9f7bOC6zPxSLb81VGeN\nNqUq3uq+lJk31R7/hOqsyQ5drCupBo+o+0Xr553trW34y/aYedfZMVqv8R6t9S/sYt31HD7RVcfM\n64HDgGe01vNw4NDMvKWH9UnSrOKlbpI0s/y/7G5wg49QXW71COBdmXnpBP0un6Dt/m1tf811Bx94\nEFVh0eneOwncBJCZd0TEEcBHgesj4nzgO8DnWh/gJ7I98McO7b9vrXf7tvY/j0sg8x8RAVUR0o1r\n2h4v7xS31T4uZkTsS3XZ4O7AhrWneh2G+spuO2bm/42IfwKeA3wiM8/tcV2SNKtY+EhSMz2QqjAB\n2LVAvE4juM2hKnD+idZ3ftrcWShl5n9HxLeAFwDPBN4HvDMinpyZvymQH1SXpXUSE7R3u3yn9jtj\nRsTjqb7fcy7VGbklwCrgUODALtc9puuR8iLinsBeVO/BQ6foLkmznoWPJDVMVKc5FlGdmTgWeHdE\nfCUzv9Gh+4M6tO0IdFOM/ImqAPh7Zv5wqs6ZeWUrn2Mj4oGtdbwVeMUEi1xN56LtIbXnR8F+VAXL\nM+tnxSLiVR36lrx53olUl/y9E/hgRPxrZh5XML4kNYrf8ZGk5nkrd41K9m/AecDHW2cI2r0iIjYd\nexARLwa2oxrBbSpnUl3O9q6xm4/WRcSWrZ8bR8SGbU9fSfVF//b2utOBbSPiJbWYc4F/bi37oy5y\nnA5rqAqaO7dBRNwfeH6HvrcCvXzPqaOI2B84ADgiMz9MNdre+yNix0FjS1JTecZHkmaOAPaJiId0\neO68zLyy9dz7gJMz83So7vsC/Br4OPCStuWWAT+NiJOBbalGUrsM+NRUyWTmzRHxBqrhrRdHxJeA\nv1MNUf0c4KdUwzzvBPwgIk4DLqG6BG4/qiGtT51kFScBr6Mavnov7hrO+jHAmzPz1qly7EG3l8N1\n8l3gLcCZEXEK1Yhrb6T6ftLD2/peCDwtIg4HrgWuzMwLeko0Ymuq9/IHmXliq/lNVINJfBbYu98X\nIklNZuEjSTNHAkdP8NwrI+Jqqkvc/kZ1P55qoczLI+KdwHERsX9r5LaxeB+g+nD+DuDuwFnAYZm5\nssO6100o89SI+Gtr+bdRncH5K9VIYye3uv0ZOAV4KtX3gVYDfwBe3OHyu/r9alZGxBOp7jH0CmAz\n4FLgkMz8fIflOuU4UXunfr2o53lORBxKtQ2OpTqb9XbgAaxb+LyFauS2fwc2pipUuil86vmdSPX3\n+5W1HJZFxGuBb0TE2zLzoz2+HklqvMgsebmxJGkmaBUU5wD7Z+bXhp2PJEnrm9/xkSRJktR4Fj6S\nJEmSGs/CR5JmL691liTNGn7HR5IkSVLjzbhR3SJiC6q7fl8FtI86JEmSJGn22Ai4P3BmZt4wWccZ\nV/hQFT1fHHYSkiRJkkbGy6hunTChmVj4XAXw6UWw886D3G8Ozr7b0wqkAyuYXyTOh/fcpkicRRee\nUyTOhw6/gYOPfdjAceaypkA2cBX3LxJnLqsLxVlbJM7pB72WU07ZgoMOuoFLLx00t1K/0stHK852\ne5SJs6TMew+fLhTn6YXinFUozrMLxfkI8J8F4txeIAbAV6bu0pXXFYpzS6E4pX5Pv1MozmsLxSl1\nX9pNB46w885wyilw0EErufTSQb8asMHA+VRKXuwyt1CczxWKU+p37Popnj8aOLKLOO23GevXvoXi\nLCgUZ/Dfjco1heIM+rf5CuAIaNUIk5mJhc9KqIqe3fcYrPC5fMN7FEnoZu5eJA5sVyTKzgvLFGLz\nF9zEDgsH/yUrVfisYusicUoVPvMKva758x/OwoXbMH/+9cCqAaOV+sM66ZniHiwrE2bDhWXiDLx9\nx5T5XYWHFopzcaE4g/+jo7IA2L1AnNsKxAD4WaE4pfbDmwrFKfV7emGhOKW2z82F4mw2cIT582Hh\nQpg//zYGHwuk1PG51O8FlPsYWOqfL6X2ob9M8fzdgV27iPOLArkA7FIozhaF4gz+u1Ep9fm31N/m\nqf8r4KhukiRJkhpvZAqfiDgsIq6MiNsi4vyIeMSwc5IkSZLUDCNR+ETES4D/orrgcg/gN8CZEbHl\nUBOTJEmS1AgjUfgAhwOfyMzPZeYfgNcDK4BDh5vW7Lb3gfcadgqSZqT9h52AJE3i+cNOQEMy9MIn\nIjYA9gR+MNaW1V1VzwYeM6y8BI878N7DTkHSjGThI2mUvWDYCWhIhl74AFtSjafYPvbg9cC205+O\nJEmSpKaZicNZA3DE/0kWLBg/tOT+BwQHvGSwIa4lSZIkjaLvAqe3tXU/BP4oFD5LgTVA+907twGu\nm2ihD31k8Pv4SJIkSZopntOa6i4BXtzV0kO/1C0zV1HdPe2pY20REa3H5w0rL0mSJEnNMQpnfACO\nARZFxIXABVSjvM0HFg0zKUmSJEnNMBKFT2ae1rpnz/uoLnH7NfDMzPz7cDOTJEmS1AQjUfgAZOaJ\nwInDzkOSJElS8wz9Oz6SJEmStL5Z+EiSJElqvJG51K1XG1+TbLJpTt1xEgdfc1qZZH5ZJsz98ylF\n4uy5/DdF4jz6b2XixE5HFolzWl5YJM6LT/1OkTi/O3CHInG+wG7A64BvAEuKxBzUr/OUInF2u/iP\nReJwa5kwPKBMmCx05JyzRZk4S/K4InG2Pf3oInH+vs+mReJstd8tReLE18scg5bk5kXibBf/WiRO\nKXl+mfd95S5l4mz0rSJhiINKvO/bUR2fzwKWDRjrqoGzqexdKA6se+/4fj2tUJxLCsUp5YBCcS4v\nFOdHheJcVSjOkwrFOXfA5bv/7OQZH0mSJEmNZ+EjSZIkqfEsfCRJkiQ1noWPJEmSpMYbicInIh4f\nEd+KiL9GxNqIeN6wc5IkSZLUHCNR+ACbAL8G3ggMNlSbJEmSJLUZieGsM/MM4AyAiIghpyNJkiSp\nYUbljI8kSZIkrTcWPpIkSZIar69L3SLiYGBpZn639fjDwGupbrl7YGZeXS7Fzg7/MCy4+/i2A58N\nB+6zvtcsSZIkafpdBFzc1ray66X7/Y7Pu4A3AETEY4DDgMOBfYFjgf36jNu1Y98OCx+6vtciSZIk\naTTs2prqlgAndbV0v4XPfYHLW/MvAL6amSdFxM+Ac/uMKUmSJEnrRb/f8bkF2KI1/wzgrNb8SmDj\nXoNFxCYRsVtE7N5q2qH1+L595idJkiRJd+r3jM9ZwKci4lfATsDprfaHAVf1EW8v4Byqe/gk8F+t\n9s8Ch/aZoyRJkiQB/Rc+hwHvp7rk7UWZeUOrfU/g1F6DZeaPcIQ5SZIkSetJX4VPZv4DeFOH9iMH\nzkiSJEmSCuv7LEtEPD4ivhAR50XEvVttL4+Ix5VLT5IkSZIG1+99fF4EfB74IrAQ2LD11AKqoa7X\n+9109nzJa4HtBoqRhxxdJJdYVOZE1x3HlMnnmOVFwnDEU8rEyReXeV3sVCYMDy4T5ssHXVEm0B6b\ntWY2A24fMNjCAZev7B47Folz1+CPTbN3oTjLikTZLk4oEgdWF4rTPtRov35fKE4Z28W/For0pEJx\nzi0SJR59TpE4sKpQnOsLxbm5QIx7AvCYxZuxFTlQpO/HVQXygTnXPbZIHIC1t8wvE2jHMn/nt8+X\nFolz9afL/KG/z6v+WCTOX37z9CJx2P2oMnHOKBRnl+7vnzOp+5xbJk4X+j3j8x7g9Zn5GsYf6X5G\nqU9ekiRJklRIv4XPzsCPO7QvBzbvPx1JkiRJKq/fwuc6oNO1MI8DCl3/I0mSJEll9Fv4fBL474h4\nFNV9d+4VES8DPgp8vFRykiRJklRCv/fx+SBV0fQDYD7VZW+3Ax/NzON7CRQR7wReSPWV89uA84Aj\nMvOyPnOTJEmSpHH6OuOTlf+gGu5kF+DRwFaZ+d4+wj0eOB54FPA0YAPg+xGxcT+5SZIkSVK7fs/4\nAJCZdwCXDBhj3NDXEXEI8DdgT+Cng8SWJEmSJOj/Pj6bAO8AngpsTduZo8zcYYCcNqf63lCZG1xI\nkiRJmvX6PePzKeCJVDcxXQID3tWrJSICOA74aWYOdCZJkiRJksb0W/g8G3hOZv6sZDLAicBD6eq2\n6GcAG7W17UK5O4ZLkiRJGh0XARe3ta3seul+C58bKXwpWkScAOwDPD4zl0y9xLOA7UqmIEmSJGlk\n7cq6JzmWACd1tXS/9/F5L/C+iJjf5/LjtIqe5wNPzsxrSsSUJEmSpDH9nvF5K/BA4PqIuApYVX8y\nMxd2GygiTgQOBJ4H3BoR27SeWp6Z3Z+7kiRJkqQJ9Fv4fKNgDq+nGhzh3Lb2VwKfK7geSZIkSbNU\nX4VPZh5dKoHM7PdyO0mSJEnqSt9FR0RsHhGvjoj/jIh7ttoWRsS9y6UnSZIkSYPr9wamDwfOBpYD\n9wc+STXK237A/YBXFMpPkiRJkgbW7xmfY4BFmfkgxg+efTrwhIGzkiRJkqSC+h3c4BHA6zq0/xXY\ntv90ptdLTl5UJM5JJ59dJM7z+WqROPfi2iJxIpYWibNbPqtInF25qEicm9m0SJxS+bDwQa2ZBzH4\nvak2HnD5MdtM3aUrqwvFuWehOLcVirNZoTgbFIpzU6E4iwvF2bFQnPsUirN1oTil3q9Snl0ozh8L\nxXlIoTil7o3+sgIxNgTg599+Er8a9O6FP33+4OkAa88oEqbyh0Jxzj2qSJirv1skDDxt1dR9uvCX\nbz5o6k7d2CXLxDnuqCJhdnjm74rEuYO7FYnzlyJRutPvGZ/b6fyXfyfg7/2nI0mSJEnl9Vv4fAv4\nt4gY+/dXRsT9gA9BodMWkiRJklRIv4XPW4FNgb9RXVvzI+By4Bbg3WVSkyRJkqQy+r2Pz3Lg6RGx\nN7AbVRG0ODN7/rJLRLweeAPV6HAAvwPel5klr2KVJEmSNIv1dMYnIjaOiH1rTftSfYt1W2CfiPhw\nRGzUYw5/Bo4AFgJ7Aj8EvhkRpb4lKUmSJGmW6/WMz8HAc4DvtB6/ieoMzdhQSQ8GlgDHdhswM9vH\n8HhPRLwBeDTw+x7zkyRJkqR19Podn5cBJ7W1HZSZT87MJwP/Bzig32QiYk5EvBSYD/y83ziSJEmS\nVNdr4bMjjLt5yUpgbe3xBcBDe00iInaJiJuphsk+EXhhZpYaXV6SJEnSLNfrpW6bM3Y3LyAzt2p7\nfk79+R78gWqQhAXA/sDnIuIJkxc/ZwDtXyfaBdi1j9VLkiRJGm0XARe3ta3seuleC5+/UFUXl07w\n/MPp4wasmbkauKL18FcR8UjgzVSjvU3gWQx+l3tJkiRJM8OurHuSYwnrfhOns14vdTsdeF+nkdsi\nYmPgSKB9sIJ+9HvmSJIkSZLW0esZnw9QDV5waUScAFzWat+ZaoS3ea0+XYuIDwDfA64B7k41gMIT\ngWf0mJskSZIkddRT4ZOZ10fEY4GPAx8EYuwp4CzgjZl5fY85bA18luq6teXAb4FnZOYPe4wjSZIk\nSR31esaHzLwSeFZE3JNqlDeAyzNzWT8JZOar+1lOkiRJkrrVc+EzplXoXFAwF0mSJElaL3od3ECS\nJEmSZhwLH0mSJEmN1/elbkN3wOtg64UDhTjtxWVSOe2qg8sE+uVRZeKMmN+8/NFl4nzhjCJx+OBR\nRcJ887gDi8TZY7tPAHsBF1KNRT8C9j2qTJx1Br7v0/0LxTm3UJzu75U2uU0LxTmkUJxPPa9MnN3L\nhCnmU0eVifOdQnFK/X4dUijOou+VicPlZcKce1SZOCX+ZNyr9fPIT8CvBjw+v+OoQbOprC4TBoCP\nHlUmzrxCcUr56QZl4vyjTBhecHSZOO85qkiYK977sCJxOLtMmOnkGR9JkiRJjWfhI0mSJKnxLHwk\nSZIkNd7IFT4R8Y6IWBsRxww7F0mSJEnNMFKFT0Q8Angt8Jth5yJJkiSpOUam8ImITYEvAK+m3Dga\nkiRJkjQ6hQ/wP8C3M/OHw05EkiRJUrOMxH18IuKlVHd+2GvYuUiSJElqnqEXPhFxH+A44GmZuarr\nBX96ONxtwfi2nQ6sJkmSJEkNcxFwcVtb93cVH3rhA+wJbAUsjohotc0FnhARbwI2zMxcZ6nHHQtb\nL5y+LCVJkiQN0a6tqW4JcFJXS49C4XM2676CRcDvgQ92LHokSZIkqQdDH9wgM2/NzEvqE3ArcENm\n/n7Y+c1uFw07AUkz0Z9OHXYGkjSxSzxGzVZDL3wm4FmekdB+DaUkdcHCR9Ios/CZtUbhUrd1ZOZT\nhp2DJEmSpOYY1TM+kiRJklSMhY8kSZKkxhvJS92msBEANxYY9+DGwUMA1VAMRSwpFaiQlRTJ6YbF\ng8cAim2fvxbKp/u7Tk1qxYprWLx4MStWXAMsLRN0UMsLbaPbyoRh40JxSv2u3lEoTkzdpSvXFIpT\nbPssh6Wlfu9LKHTsuHzEjmWjdmwt5bJCr+u6wUOsWAOLF1Pm+Hxdode1pkyYSqH3fsko/b4z9Ta6\nfXl378fNRbKhsdt5ZD7/3vm7udFUPWOmjRYdEQcBXxx2HpIkSZJGxssy85TJOszEwmcL4JnAVfRy\nq1ZJkiRJTbMRcH/gzMy8YbKOM67wkSRJkqReObiBJEmSpMaz8JEkSZLUeBY+kiRJkhrPwkeSJElS\n41n4SJIkSWo8Cx91FBGPjojVEfHtYeciabRFxMkRsbY2LY2I70XErsPOTZLGRMQ2EXF8RPwpIlZG\nxNUR8a2IeMqwc9P0sPDRRF4FfAx4QkRsO+xkJI287wHbANsCTwFWA/7jRNJIiIjtgcXAk4C3ArsA\nzwLOAU4YXmaaTt7HR+uIiE2AJcCewNHAbzPzg8PNStKoioiTgQWZuV+tbW/gx8DWU91QTpLWt4g4\nnarY2SkzV7Y9t1lm3jSczDSdPOOjTl4C/D4z/wh8kersjyR1JSI2BV4O/NGiR9KwRcQ9gGcCJ7QX\nPQAWPbPHvGEnoJF0KPD51vwZwGYR8YTM/PEQc5I02p4bETe35jcBrgX2HWI+kjRmRyCAS4ediIbL\nMz4aJyJ2Bh4JfAkgM9cAp+FZH0mT+yHwcGA34BHAmcAZEXHfoWYlSVXRI1n4aB2vAuYCSyJiVUSs\nAl4PvCgi7j7c1CSNsFsz88rMvCIzLwReQ3Xm5zVDzkuS/ggk8OBhJ6LhsvDRnSJiLtV1+W+h+q9t\nfboWOHB42UmagRLYeNhJSJrdMvNGqrPQh0XEOsekiFgw/VlpGCx8VPdcYHPgM5l5SX0Cvga8erjp\nSRphG7bukbFNRDwYOB6YD3xryHlJEsBhVFe0XBAR+0XEjhHx4Ij4F+C8IeemaWLho7pDgbMy8+YO\nz30V2DMidpnmnCTNDM+iOjN8LXA+1XD4+2fmT4aalSQBmXklsJDqvj0fBS4Cvg88g+pKF80C3sdH\nkiRJUuN5xkeSJElS41n4SJIkSWo8Cx9JkiRJjWfhI0mSJKnxLHwkSZIkNZ6FjyRJkqTGs/CRJEmS\n1HgWPpIkSZIaz8JHkiRJUuNZ+EiSZryIODIiflV7fHJEfG2YOUmSRouFjyRp6CJim4g4PiL+FBEr\nI+LqiPhWRDylhzC53hKUJM1484adgCRpdouI7YHzgGXAW4GLgQ2AZwEnAA8dXnaSpKbwjI8kadg+\nDqwBHpGZ38jMyzPz95l5LPBogIhYEBGfioi/RcTyiDg7Ih7e7QoiYv+I+G1ErIiIpRHx/YjYeD29\nHknSCLLwkSQNTUTcA3gmcEJmrmx/PjNvas1+Bdii1XchsBg4OyI272Id2wKnAJ8CHgw8EfgaECVe\ngyRpZvBSN0nSMO1IVYBcOlGHiNgb2AvYOjNXtZrfHhEvBPanKmgmsx0wF/h6Zv651fa7gbKWJM04\nFj6SpGHq5qzLbsDdgWUR47pvBDywi+V/A/wAuDgizgS+D3wlM//RY66SpBnMwkeSNEx/pBqN7cHA\nNyfosylwLdUlau2F0pTFS2auBZ4REY8BngH8M/D+iHhUZl7db+KSpJnF7/hIkoYmM28EzgQO6zTY\nQEQsoPo+z7bAmsy8om1a1sO6fp6ZRwN7AKuAF5Z5FZKkmcDCR5I0bIdRfQfngojYLyJ2jIgHR8S/\nAOdl5tnA+cA3IuLpEbF9RDw2It4fEQunCh4Rj4yId0bEnhFxX+BFwJbAJevzRUmSRouXukmShioz\nr2wVMO8GPko1GMHfgd8Cb2l1ezbwH8BngK2A64AfA9d3sYqbgCcAbwY2A64G3pKZ3y/4MiRJIy4y\nvdG1JEmSpGbzUjdJkiRJjWfhI0mSJKnxLHwkSZIkNZ6FjyRJkqTGs/CRJEmS1HgWPpIkSZIaz8JH\nkiRJUuNZ+EiSJElqPAsfSZIkSY1n4SNJkiSp8Sx8JEmSJDWehY8kSZKkxrPwkSRJktR4Fj6SJEmS\nGs/CR5IkSVLjWfhIkiRJajwLH0mSJEmNZ+EjSZIkqfEsfCRJkiQ1noWPJEmSpMaz8JEkSZLUeBY+\nkiRJkhrPwkeSJElS41n4SJIkSWo8Cx9JkiRJjWfhI0mSJKnxLHwkSZIkNZ6FjyRJkqTGs/CRJEmS\n1HgWPpIkSZIaz8JHkiRJUuNZ+EiSJElqPAsfSZIkSY1n4SNJkiSp8Sx8JEmSJDWehY8kSZKkxrPw\nkSRJktR4Fj6SJEmSGs/CR5IkSVLjWfhIkiRJajwLH0mSJEmNN2/YCfQjIu4HbDnsPCRJkiQN3dLM\nvGaqTpGZ05FMMVXRs8HVsGrYqUiSJEkavhXAQ6YqfmbiGZ8tq6Jnv2p23Euoz88t1N7N/PpeVzdx\n+ll+ovkedZPeRPPd9K9fkNnNSyzVf6I8u+nf6+vtZ5lS23d9x+912V7fv3HztX/kzF1TW3btnbMx\n5672ufNq83Pv6jNn7tR9AObUY1Fbftz86tpLmKhPPU7nmPX2Xped01V7qT6dX+Pc2nYol1v/eU72\n3LyutvV0bsfetsX4faXX7djr/tpjnmtq7ePms9aHWvtd8zH+149aqox7a+vz9WUm6r92gvZu+g8S\nc6I+g/Qf5fn1Fbeb92OQ92897Dera+2ra/1X1/qsWdO5/0Srmmi1E/Xppn/72zfIOnpdXy/rWgp8\nDeZTFQaNK3xatgTuxfiXsEFtfn23T+e6umkvmV8XYoL5+ofWbj7kTpTObOs/6DpKLTvj49cKn1rB\nUp+PuXcdNudsUPuAWC9wxs3XPgjOG/9nYO64wmei+foHyfXRfzrXNcq59bbsdK9vlLfF+Pm5XfSf\nM0H7XSb8P0WtwJlb+5Q3b1w7Hedjsk9hg3wCLPUp0fVO3md9xR217ThngvZa//o1S6tzgva1E7TT\nub2bPoO0j0oek+XXDQc3kCRJktR4Fj6SJEmSGs/CR5IkSVLjWfhIkiRJajwLH0mSJEmNZ+EjSZIk\nqfEsfCRJkiQ1noWPJEmSpMaz8JEkSZLUeBY+kiRJkhrPwkeSJElS41n4SJIkSWo8Cx9JkiRJjWfh\nI0mSJKnxLHwkSZIkNZ6FjyRJkqTGs/CRJEmS1HgWPpIkSZIaz8JHkiRJUuNZ+EiSJElqPAsfSZIk\nSY1n4SPinqYvAAAJKklEQVRJkiSp8Sx8JEmSJDWehY8kSZKkxrPwkSRJktR4Fj6SJEmSGs/CR5Ik\nSVLjWfhIkiRJajwLH0mSJEmNZ+GjKZw67ARmh8vcztMlv37asFOYNf506i+HncKscPapNww7hVnj\n1LOGncHscOo1w85g9vj+sBOYZhY+moIfyKeFhc+0ya9/ZdgpzBpXWPhMix9Y+EybU88edgazg4XP\n9LHwkSRJkqSGsfCRJEmS1HgWPpIkSZIab96wE+jf0tbP+kuoz88t1N7N/PpeVzdx+ll+ovm65cDi\ndZtzgvm1tfk1E4TMCdonWrZenq+uzU/0EudM0KfX/hNttm76TzQ/Uf/bl8P1izs/1+s6el12fcfv\nddle379x87Wda25tJ5pX27luWk7+9tcArJ13V5+YW+szbtnafL0PsHbOXc+tre3Aa2rzc2o77fj2\n+vxdceaOm1/bsb3XZed01V6qz13tdyxfydLFf261r+5p2e7W23+ekz03r6ttPX3bcaplb12+hksX\nr+jYp/fteNf7NK+rbd3ja1xT26fXZG2eWh86tsf4X7/xfw/WTDBfX2ai/hP97enQf/ktsPiyQjEn\n6jNI/1Ge72GZ5atg8Y1d9u/mPR5knyiw37T3X71mgvlan9qvx4QhV08w302fsflbgEu7iDPIOrpd\nttf2sfmldC8yJ/oUOpoi4n7AZcCGw85FkiRJ0tCtAB6SmZMOjTHjCh+4s/jZcth5SJIkSRq6pVMV\nPTBDCx9JkiRJ6oWDG0iSJElqPAsfSZIkSY1n4SNJkiSp8Sx8JEmSJDWehY8kSZKkxptxhU9EHBYR\nV0bEbRFxfkQ8Ytg5zWQR8c6IuCAiboqI6yPi6xGxU1ufkyNibdt0+rBynqki4sgO2/GStj7vi4hr\nI2JFRJwVETsOK9+ZqnV8aN/OayPi+Nbzi9yf+xMRj4+Ib0XEX1vb7Xkd+ky6D0fEhhHxPxGxNCJu\njoivRMTW0/cqRt9k2zki5kXEhyLitxFxS6vPZyNiu7YY57bt42si4sTpfzWjbap9upu/f+7TU+ti\nO4/to+3b+q21Pu7TU+jmM12r36w9Ts+owiciXgL8F3AksAfwG+DMiPCePv17PHA88CjgacAGwPcj\nYuO2ft8DtgG2bU0HTmeSDXIx47fj48aeiIgjgDcBrwUeCdxKtX/fbQh5zmR7cdf23RZ4OpDAaa3n\nE/fnfm0C/Bp4I9V2HKfLffg44DnAi4AnAPcCvrp+055xJtvO84HdgaOp/g6+ENgZ+GZbvwRO4q79\nfDvg7esv5Rlr0n26Zarjhfv01KbazmP76Ng2PhRYC3yl1sd9empTfqab9cfpzJwxE3A+8N+1xwH8\nBXj7sHNrykR1Y9i1wONqbScDXxt2bjN9oirYF0/y/LXA4bXHmwG3AQcMO/eZPFEdwC+rPXZ/LrNd\n1wLPa2ubdB9uPb4deGGtz86tWI8c9msaxanTdu7QZy9gDXCfWts5wDHDzn8mTRPs05MeL9yny2zn\nDn2+AZzV1uY+3fu27vSZblYfp2fMGZ+I2ADYE/jBWFtW78bZwGOGlVcDbU71X5Vlbe1Pap02/UNE\nnBgR9xxCbk3woNap/j9FxBci4r4AEfEAqv9g1ffvm4Bf4P7dt9Zx42XAp9uecn8urMt9eC9gXluf\nS4FrcD8fxNhx+x9t7S+LiL9HxEUR8YEOZ/LVncmOF3viPl1U65KqfYBPdXjafbo34z7TeZyuXthM\nsSUwF7i+rf16qkpUA4qIoPrv+E8zs/7dk+9RneK8Engg8J/A6RHxmFbxqe6cDxwCXEp1iv4o4McR\nsQvVgSjpvH9vO30pNs4LgQXAZ2tt7s/rRzf78DbAHa0/tBP1UQ8iYkPgg8ApmXlL7akvAldT/Xf3\n4cCHgZ2A/ac9yZltquPFtrhPl3YIcBPw9bZ29+keTPCZbtYfp2dS4aP170TgocDe9cbMPK328HcR\ncRHwJ+BJVKee1YXMPLP28OKIuIDqIH4A8IfhZNV4hwLfy8zrxhrcn9UUETEP+DLVB5k31p/LzPp/\ny38XEUuAH0TEAzLzymlMc0bzeDEUrwS+kJl31Bvdp3vW8TPdbDdjLnUDllJdw7xNW/s2wHXrdlcv\nIuIEqlPLT8rMJZP1bR1glgKOODaAzFwOXEa1Ha+j+s6a+3chEXE/qi93fnKyfu7PxXSzD18H3C0i\nNpukj7pQK3ruCzyj7WxPJxdQvT/u5wPocLxwny4oIh5PdRan02Vu7dynJzDJZ7pZf5yeMYVPZq4C\nLgSeOtbWOo33VOC8YeXVBK1fkOcDT87Ma7rofx9gC2DSAkmTi4hNqQ7Y17b+mF7H+P17M6qRWdy/\n+3Mo1an5SYeqdn8uo8t9+EJgdVufnYH7AT+ftmRnuFrRswPw1My8sYvF9qA6M+R+PoAOxwv36bJe\nBVyYmRd30dd9uoPJPtN5nJ55l7odAyyKiAupKv3DqYb2XDTMpGay1hj4BwLPA26NiLH/AizPzJUR\nsQnVaGRfpfpl2RH4ENWZijM7hNQEIuIjwLepLm+7N9VwtKuAL7W6HAe8JyIuB64C/p1q1ML2YWo1\nhdY/RQ4BFmXm2lq7+/MAWttvR6r/GALsEBG7Acsy889MsQ9n5k0R8WngmIi4EbgZ+Bjws8y8YFpf\nzAibbDtTfcj7KtWQ1vsCG9SO28syc1VE7AAcRFX03wDsRvX380ddfqCcNabY1suY4njhPt2dLo4d\nYx/A96f6bNe+vPt0F6b6TNean93H6WEPK9frRHUd81VUQ+/9HNhr2DnN5IlqeMI1HaZXtJ7fCDiD\n6qC/ErgC+Diw1bBzn2kTcCrVweU2qtFRTgEe0NbnKKovbq6g+sO647DznokT1b171rRvP/fngbfr\nEyc4Znym1mfSfRjYkOo+E0up/qB+Gdh62K9tlKbJtjOwfYfnxh4/obX8fYBzgb+33odLqb6Uv+mw\nX9uoTVNs666OF+7Tg23nWp/XALcAd++wvPt0d9t50s90tX6z9jgdrRcoSZIkSY01Y77jI0mSJEn9\nsvCRJEmS1HgWPpIkSZIaz8JHkiRJUuNZ+EiSJElqPAsfSZIkSY1n4SNJkiSp8Sx8JEmSJDWehY8k\nSZKkxrPwkSRJktR4Fj6SJEmSGu//A2tt33eXQmYFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xd9f4fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAADyCAYAAAB6QJUSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xm4JFV9//H3Z2ZwAWVER1kSgRANiwKCigtxT1QQNEaj\nJv4iARFcosYl4ApEjCgaJYAo4gYuE4lxQRBQBAUBBdk3ZVVQ1mEZ9mVmzu+PqsvUNN23+/bt232X\n9+t56pnuU6eqvnX63J7+dp06nVIKkiRJkqTO5o06AEmSJEma7kycJEmSJKkLEydJkiRJ6sLESZIk\nSZK6MHGSJEmSpC5MnCRJkiSpCxMnSZIkSerCxEmSJEmSujBxkiRJkqQuTJwkaY5Isk+SFaOOYxCS\nrEiy16jj6EeSDer43zTqWPoxk9tekibDxEmSBiTJTvWHyq2HdLxNk+ydZP0eNynAjEmckmyXZO8O\nq0u9TOXx35Zkp6k8RpfjT/T1lSRNIRMnSRqsKf0w32IzYG9gwx7r7wusPmXRDN72QKcrG48E/nOK\nj/92YGSJExN/fSVJU2jBqAOQJPUtTCBRK6WsAO6funAGLp1WlFJm0nn0a0KvryRpannFSZKGKMlq\nST6W5DdJbktyZ5KTk7ywTd031PVuT7I0yflJ3lmv2wk4sq7683qI4PIkzx/n2A+5x6ne7sAkr0py\nQZJ7k1yY5GWNOq+p6z2vzT53r9dt1ijbOMl3k9yc5J4kZybZsWW7BfUwtEvrOkuSnJLkJfX6r1Fd\n8RmLcUWS5S1x79WyzxfW7XVPksuS7NbhnHdO8rMkN9Tne1GSt7bUuQp4CvDCxvFPbKxfmOSAJFfX\n+7gsyR5J0rKfhUm+Xr/Wt9bn9ZhOr1Fju66vb5K316/VvUn+lOTgJAtb9vPzut9sneTUJHcnuTLJ\n7m2O+fgkX0lyfd2G5/Z6H1aSrZIcW/fTO5KckORZbeptkeQXdRzXJPlw/XqsGBuSWLfXTUnmt9n+\nJ0ku6SUmSRo0rzhJ0nCtCewCLAa+BDwaeDNwXJJtSinnAyT5W+DbwE+BPeptNwW2BQ4CTgYOBN4J\nfBz4bV1nvA+Vne4Leh7w98AhwB3Au4DvJlm/lHIrcAxwJ/A64JSWbV8HXFhKubiO+ynAL4E/AvsB\nd9V1fpDk70spP6y3+w/gA3UbnFm3yzOArYGfAV8E1gP+Bngj41x9qo+7FXAscC3wUar/3z4KLGlz\nzm8FLgR+CCwDdgQOSZJSyhfqOu8GDq7b4+P18W+oj/VIqvZft47zGuC59fmuA7y3cayj6nVfoHqN\nXg0c3iamVr9gnNc3yT5Uwxh/QvW6bUyVaD4jyballLEkswCPpXoNj6TqU68DvpDkvlLK1+v9PaI+\n5kZU/ev3wD8AX0+ysJRyUKdA66T5ZGAp8EmqNt2dKuF7finlzLreesBJwHKqYZZ3A7tSXQVttsc3\ngH8GXgb8uHGctYEXUQ1flKThK6W4uLi4uAxgobofZjmw9Th1AixoKVsTuA44rFH2OeDWLsd7TX28\n5/cY397A8payFcA9wIaNss3r8rc3yr5Vx5hG2dpUH5I/1Cg7ATinzTn+Evht4/k5wFFd4j2oNd6W\nuPdqPD+KKslZu1G2EdWH8tZzfnib/R0LXNZSdgFwYpu6HwFuBzZqKf9Efbw/q5+/qo7zvS2v/y/q\n1+1N/by+wCLgXuDHLeVvr+vv1CgbS1Te3ShbDTi7fj3n12Xvruu9oVFvPnAqVUK0xjht//26D23Q\nKFun3u6kRtmBdX/ZvFH2GKrkdjmwfqONrga+3XJ+76m336C0aS8XFxeXqV4cqidJQ1QqywBSWQt4\nGPAbqqstY24D1mgOmZtCPy2l/L4R4wXUiUGjzneAJwAvbJT9A9WH3CMB6nN5EfC/wMIkjxtbqK6M\nPDnJuvW2twFPSfKkyQafZB7wEuAHpZQbGudxJVVCtIpSyn2Nbdes4zsZ2CjJo3s45GuprrwtbTnH\nn1Fd6RobTrc98ADVVamxYxeqhHDcK2hd/A1V8nNAS/lhVMnjK1rKl1Fd2RuL4QHgUKrX8+l18XbA\n9aWU/2nUW06V7DwKeEG7QOq2/1vg+6WUPzS2vZ7q6tZfJ3lUXfwy4PS6f43Vu40qKadRVuqyVyZZ\no7Hqn4DTmseRpGEycZKkIUs1bfl5VFcNbgZupPqw27w/5RDgUuDH9b0gX5nCJOqaNmW3Ams1nh9H\nlUy9vlH2OuDcUsrl9fMnUSUE+wI3tSz71HWeUP+7F9XVhkvre3D2T7J5n/E/gWqWvcvbrHtIWZJt\n63tw7qRK4G5i5Qx9C1vrt/Fk4OU89Bx/SjXkbOwc1weuK6Xc3bL973o4xng2qP+9tFlYJ0RXNtaP\nubaUck9L2aVUr9WGjX1e1uZYl9T1Wvc55vFUMzVe2mbdJVSfM57YOEZPrxFwRL3fV0N13xxVkndE\nhzgkacp5j5MkDVGS/wd8DfgesD9V0rQc+BCNKzyllJuSPI3qW/rt6mXnJIeXUnYecFjLO5Q/eFWk\nlHJ/kh8Ar07ydqr7e7aluk9pzNiXcZ8Bju+wz8vr/Z2S5C+phrO9lOo+r/ck2b2U8tW+z6SLJBtR\nDSe8hGro1zVUw+teAfwbvX2hOI8qSfoU7a8ctUsiNAGllEuSnAX8P+Cb9b/3UV3NlKSRMHGSpOF6\nDXBFKeW1zcIkH2utWA/pO6ZeSPIFYLck+9bD0IY9VfV3gDdRDYt7Sl12ZGP9lfW/D5RSTqSLepjW\n4cDhSVanGv62DzCWOPV6fjdSXb1rN+zvyS3Pd6QaGrljKeVPY4WpZ/NrDbHD8a4AHlVKOalLXH8A\nXpxk9ZarTpt02a7b8ceGqm1MNYkDUM3YCPwFVVLXtF6SR7Zcddq43v9VjX22u+K3acsxW91ENcnD\nxh22XcHKK5p/oLfXaMwRwH8lWQf4R+CYUsrSDnUlaco5VE+ShushV3fqaZuf01L22Dbbjt0b8vD6\n37uornh0nd56QE6gGsL3Bqpheme03NdyE/BzYPf6w+4qkixqPF7l/OrE4nJWnhtU50eSNccLqlS/\nT3UC8HfN49b3T728pfpY+89r1FsI/EubXd9F+7Y9EnhOkpe2rkg1/fjYvn9MdS/S2xrr51HNlNdL\nUtjp9T2B6t6pd7WU70o10cjRLeULqGYSHIthNapZ726imiRiLNZ1kry+UW9+HesdVBNaPETd9j8B\nXjU2nXi97dpUyc4ppZQ76+Ljqdpti0a9x1Ldu9TO4vrf/6ZKCL/RoZ4kDYVXnCRpsAK8Ocl2bdYd\nQPWh9u/rYW/HUA3P2x24iOom/DFfrj9Unkg1tfeGwL8C55RSxqYcP5cqEdgzyWOohjL9rJSyZOBn\nRXUFLMn3qBKn1YH3tan2DqorRxckOYzqKtTaVInhnwFb1fUuTvJz4CzgFuCZVJMuHNjY11lU7XlQ\nkuOpZsf7Tofw9qEa8ndafWVuQR3LhcCWjXo/oUo6jk5yKNV08LtSTTXemuydBbw1yYepkrob66tM\nnwZeWe/j63W9NYAtqKZ137A+px9RzUr3ySR/AVxcr+9lAgoY5/VNsh+wV5LjqGYU3IQqQTuDlskW\nqKZo3yPJhlTDCN9Qx/qWsnLa8i9R9cOvJ3kGK6cjfw7VjHx3jRPnR6gmrDg1ySF1zLtRXdnbo1Fv\nf6ohdyckOYgqMdyV6krUWrQkk/V5HlfHcSuNqcklaSRGPa2fi4uLy2xZWDkdeadlvbrenlQJxd1U\ns+ltR3Xf0xWNfb2aaka466imer4K+DzwhJZj7kJ1U//9dJmanGo68mUtZcuB/25T90rgK23KX1Jv\n88DY+bSps2F9Pn+iGkJ3NdVvJr26UeeDwOlUk2PcSZU47kk9PXZdZx5Vsnk91cxwy1vi/mjLcV9Y\nt+c9dZvsSpXk3NVS7xVU06HfRTXs7n1UV5wenBK7rvcEqqTktnrdiY11q1P9vtLv6uPdQJUw/lvL\nOTwG+DrVB/9b6nbZgh6mI+/2+lIlShfVbXwt1Wx9a7ZsfxJwPlXCemp9zlcCb21zrEXAl+tzuYcq\ncfvnNvXatf2WVInNUqorVD8Ftmmz7RZUVyXvrvvFB6muai0HHt+m/muphvsdMuq/bxcXF5eUMuwh\n8pIkDUeS7wOblVLa3YMz6yU5CXhcKWWLrpVHJMkBwFuo7hsrLeteSfU7Uc8rpZw2ivgkaYz3OEmS\nZoUkj2h5/mSq31LqNomDhqTNa/Q4quF7p7QmTbXdgCtNmiRNB97jJEmaLa6s7zm6kmq44FuphrF9\neoQxaVWn1/e2XUJ1T9kuVPd87duslGTsPqzteOgkGJI0EiZOkqTZ4liqiQ/WoZpI4TTgQ6WUK0Ya\n1ehNpzH5x1Ddt/QWqrjOAnYupZzaUu/bVPdKfRn4wlAjlKQOvMdJkiRJkrqYk1ec6jHVL6OabvXe\n0UYjSZIkaYQeQTXE+/hSys2dKs3JxIkqaWr9nQtJkiRJc9cbqYYKtzXhxCnJc0opp08qpNH7PVS/\nQrioQ4XjeOjPzUtTwb6mYbGvaVjsaxoW+5oGYQnwverh78er188Vp1OSXAp8FfhGKeWGPvYxavdC\nlTSt26HCI8ZZJw2SfU3DYl/TsNjXNCz2NQ3YuLfw9PM7TusD36D6RfZrkvwwySuTzO8nukFJ8uwk\ny5L8aJRxSJIkSZpeBjEd3oQTp1LKtaWU/UopmwAvBG4EjgD+mGT/JJsMIK5+vBk4EHh+knVGFIMk\nSZKkaeA+4PSFCzl6o404+elP5+iNNuL0hQu5r8/99XPF6UH1L3l/kOrHBdei+pG6C5P8JMnGk9n3\nRCRZA3g91W89HAP8y7COLUmSJGl6uQ84dt11efe3vsWZl1/OSb/5DWdefjnv+uY3OXbddftKnvpK\nnJLMS7JDku8Bf6L6Mbv3U/3o4GbAUuC7/ey7T68HLimlXEY1W96bJ7vDp046JKk39jUNi31Nw2Jf\n07DY19TJ2QsX8onDDmP7V7yCJAAk4RU77MDHDz2UcxYunPA+J5w4Jdkf+CNVgnIjsG0pZatSysGl\nlNtKKZdSXXnabMLR9G8XqvuuoJpgZc0kz5/MDjefdEhSb+xrGhb7mobFvqZhsa+pk5sf9zi22377\ntutescMOLFnUaW7tzvqZVW9b4MPAd0opd3eocxPQPtIBq4cEbgP8HUApZXmSI6muOp083rbHUc3G\n0vRU/COUJEmSZqoCPHqttR680tS0ePFiFi9ezO+XLGFxXTbuVHoNKWUQc0yMTn0F7H3AipZV9wHr\nllLuaLPN1sBZu+EUlpIkSdJsc/RGG3Hm5Ze3TZ5KKTzzSU9ihyuvBOA64EvVqqeXUs7utM9+fgD3\npR1WFaqE7YpSyrUT3W8/6inQ/xl4L/DTltU/AP6RB9tBkiRJ0lzwuJtv5sfHHMMrdtjhIeuOOfpo\nFt1884T32c9QveNYORX6WArXfF6SnAC8rpSytI/9T8SOwGOAr7ZeWaonrtgVEydJkiRpTtl66VI+\nvNtulEMP5RU77EASSikcc/TRfGT33dlu6cTTlAkP1UvyEuBTwN7AGXXxNsA+wL7AXcDngV+UUnaf\ncEQTi+UogFLKK9useybwK2DLUsqFLescqidJkiTNYvcB5yxcyJJFi3j0woXcsXQpi5YsYaulS3l4\no96UDdUD/gt4ZynllEbZMUnuAA4upWyR5N0M4UpPu4Spse5MYP5UxyBJkiRp+nk48OylS2HpUgor\nh8r1q5/fcdoYaDco8GbgyfXj3wKP7zcoSZIkSRqUySZN0N8Vp3OB/ZLsVEq5DSDJWsB+wDl1nb+k\n+mHcae1L/BrYagRHfmAExxyzbITHnux5z+TYR3XsudpmM/m8Z3LskzGZ8x5l3DB3XzPfmybO8+7P\nTI59Mjzv4WkM1htHP4nTrsBRwLVJrqzLNqqPODZ07rHAJ/vYtyRJkiRNOxNOnEopF9U/OrsD8Fd1\n8e+AY0opy+o63x1ciJIkSZI0Wv1ccaJOkH4w4Fj6kuRrwE6NoluAM4E9SikXjCYqSZIkSbNJP5ND\nkORZSf43yYX1cmSSbQYd3AQcC6wNrAO8mGpg5o9GGI8kSZKkWWTCiVOS1wGnAg8DjqiXhwOnJvmH\nwYbXs/tKKTeVUm4spZxPdX/VE5M8bkTxSJIkSZpF+hmqtzfw4VLKp5qFSfak+hHc/x1AXH1L8ijg\nn4HLSintpk2XJEmSpAnpZ6jek4D/a1P+f1TTkI/CjknuqH+E93aqiSveMKJYJEmSJM0y/SROfwKe\n36b8BYzut5tOBLYAtgSeCRwPHJfkiSOKR5IkSdIs0s9QvQOAzyfZHDitLtsW2A3Yc1CBTdBdpZSr\nxp4keQuwFHgLsFfnzd4PrNlS9nq8WCVJkiTNRhcAF7aU3dvTlv38jtOBSW4C3keVmAD8Fti5lPKd\nie5vChXgkeNX+Qyw1TBikSRJkjRym9dL03XAl7pu2e/vOC0GFvez7RR5eJK168drAe8EVgeOGl1I\nkiRJkmaLvhKnaejlwLX14zuoroC9tpRyyuhCkiRJkjRb9JQ4JbmOauhbV6WU9SYV0QSVUnYGdh7m\nMSVJkiTNLb1ecdpnKoOQJEmSpOmsp8SplHLoVAciSZIkSdNVX/c4JQmwPbBpXXQRcFwppafhfNPF\n23kWQx1XKEmSpAmZLTfkD9Nqow5gRPrtK9dQzbU98P0n2RA4GngycEVd/JfA75LsWEr5w0T3KUmS\nJEnT2bw+tjkIuB5Yv5SyWSllM2AD4CbgwEEGJ0mSJEnTQT+J04uA95ZSbhgrKKVcT/WDuC8eVGAT\nkWTtJAcluSLJvUn+kOSoJCOJR5IkSdLs0s9QwGXAI9uUP6JeN1RJNgBOA26hSt4upBra+XLgYGCz\nYcckSZIkaXbpJ3H6MfDFJP9SSjkPIMnTgC8AxwwyuB59AVgOPLOUcm+j/JIkXxlBPJIkSZJmmX6G\n6r0TuBE4J8mdSe4EzqK67+ldgwyumyRrAS8DDm5JmgAopdw+zHgkSZIkzU4TvuJUSrkZeFmSzYFN\n6uJLSikXDjSy3jwJCPC7ERxbkiRJ0hzR99T4pZQLgAsGGEs/MpmNf0x1Y1bT5sCWk9mpJEmSpGnp\nLODslrJ7ety23x/AfSPwHlZecfotcEAp5Zv97G8SLgNKHccPJ7rx9uAP4EqSJElzxNPrpanXH8Cd\n8D1OST4CHAacAuxSL6cAh9brhqaUcitwPPCOJA+Z6S/JwmHGI0mSJGl26ueK07uA3Usp32iUHZnk\nHKpk7eMDiax37wB+CZyRZG/gfKrzeimwO/CUIccjSZIkaZbpJ3F6OPCrNuWn89BbhqZcKeWqJFsD\nH6ZK3NYFbqJKoN477HgkSZIkzT79JE6LgTcDH2gp3xn4n0lH1IdSyg1UV8KGOh26JEmSpLmhn8Tp\nLuBtSf6W6ioTwLOBJwOHJfnEWMVSyocmH6IkSZIkjVY/idNzgYvrx2OTUjxQl23bqFcmEZckSZLE\nslEHMAPZZhNzd4/1+vkB3OdMdBtJkiRJmskmPB15U5JFSRYNKpg+Y3hBkhVJ1qyf75Tk1lHGJEmS\nJGl26ed3nJJkjyQ3AjcANyS5Mcm/J8ngQ+xJ67BAhwlKkiRJGph+7nH6D6rfTvo4cGpd9tdU04Gv\nAewzkMgkSZIkaZroJ3F6M7BrKeX7jbIzkvwB+G8mmDglOYnqN5fuBXYF7ge+WEr5jyQbAFcBTyul\nnF/XXwjcCrywlHJyH/FLkiRJ0oT0c4/T44CL2pRfUK/rx5uAO4FtgD2AvZK8pF7nsDtJkiRJI9VP\n4nQhsFub8t3rdf04v5SybynlilLKN4DfAGOJ06jum5IkSZIkoL+heh8AflRfETqtLnsusDGwQ59x\nnN/y/DrgCX3uS5IkSZIGqp/fcTohySbAO4FN6+KfAX9XSvlDn3E80HoYqqthK+rnzatOq/V5jIf4\nMfCIlrLNgS0HdQBJkiRJ08Z5VPcXNd3b47Y9J05J9gI+U0q5u06Q3t/rtpNwU/3vulTnCbAVA7rv\naXtgvQ7rzsMESsNhX9Ow2Nc0LPY1DYt9TRO1JQ/tM9cCh/Sw7UTucdobeNQE6k9aKeVe4FfAB5Js\nkuQFwL5tqg78PqjWTFSaKvY1DYt9TcNiX9Ow2Nc0TBNJnKZqkoZuV492oboy9hvgs1S/FzXRfUiS\nJElS3yZ6j9PAE5RSyovblL268fi3VD+w2zS/sf4XLc8PBw4fdJySJEmS5q6JJk6XJhk3eSqlPHYS\n8UiSJEnStDPRxGlvYOlUBDJkj4CVM0+0cy/VjWLSVLOvaVjsaxoW+5qGxb6mQWjkBK0Tbq8ipfQ2\n+i7JCmCdUsqNkwlsOkjyT8C3Rh2HJEmSpGnjjaWUb3daOZErTrNpAobjgTcCv6f3qdslSZIkzT6P\nADakyhE6mpNXnCRJkiRpInpOnCRJkiRprprI7zhJkiRJ0pxk4iRJkiRJXZg4SZIkSVIXJk6SJEmS\n1IWJU4sk70hyVZJ7kvwqyTNHHZNmtiTPS3JUkj8lWZHklW3qfCzJtUnuTvLTJE8aRaya2ZJ8MMkZ\nSW5PckOS7yf5qzb17G+alCRvTXJekqX1clqSl7fUsZ9poJJ8oP5/9LMt5fY1DYWJU0OS1wP/BewN\nbAWcBxyfZNFIA9NMtwZwLvB22vweWpI9gX8FdgO2Ae6i6ncPG2aQmhWeBxwEPAv4G2A14CdJHjlW\nwf6mAbkG2BPYGng6cCLwwySbgv1Mg1d/kb0b1WezZrl9TUPjdOQNSX4F/LqU8u76eaj+cziwlLL/\nSIPTrFD/HtrflVKOapRdC3y6lPK5+vmawA3ATqWUI0cTqWaD+kufG4Hnl1J+WZfZ3zQlktwMvL+U\n8jX7mQYpyaOAs4C3AR8FzimlvLdeZ1/T0HjFqZZkNapvzX42VlaqrPIE4DmjikuzW5K/ANZh1X53\nO/Br7HeavMdQXeW8BexvmhpJ5iV5A7A6cJr9TFPg88CPSiknNgvtaxq2BaMOYBpZBMyn+pai6QZg\n4+GHozliHaoPtu363TrDD0ezRX3F/ADgl6WUi+ti+5sGJslTgdOBRwB3AK8upfwuyXOwn2lA6qT8\nacAz2qz2PU1DZeIkSbPTIcBmwLajDkSz1m+BLYGFwGuBI5I8f7QhaTZJ8udUXwD9TSnlgVHHIzlU\nb6UlwHJg7ZbytYHrhx+O5ojrgWC/0wAlORjYHnhhKeW6xir7mwamlLKslHJlKeWcUsqHqW7afzf2\nMw3O04HHA2cneSDJA8ALgHcnuZ/qypJ9TUNj4lSrv8k4C3jJWFk91OUlwGmjikuzWynlKqo392a/\nW5NqVjT7nSasTppeBbyolHJ1c539TVNsHvBw+5kG6ARgc6qhelvWy2+AbwJbllKuxL6mIXKo3qo+\nC3w9yVnAGcB7qG52/foog9LMlmQN4ElU34oBbJRkS+CWUso1VMMQPpLkcuD3wL7AH4EfjiBczWBJ\nDgH+EXglcFeSsW9hl5ZS7q0f2980aUk+ARwLXA08Gngj1ZWAl9ZV7GeatFLKXcDFzbIkdwE3l1Iu\nqYvsaxoaE6eGUsqR9fS9H6O6zHsu8LJSyk2jjUwz3DOAk6huYC1UvxUGcDiwSyll/ySrA4dSzYJ2\nCrBdKeX+UQSrGe2tVH3s5y3lOwNHANjfNCBPoHoPWxdYCpwPvHRs1jP7mabQKr+jY1/TMPk7TpIk\nSZLUhfc4SZIkSVIXJk6SJEmS1IWJkyRJkiR1YeIkSZIkSV2YOEmSJElSFyZOkiRJktSFiZMkSZIk\ndWHiJEmSJEldmDhJkmaMJBskWZFki1HHAr3Fk+QFdZ01hxmbJGmwFow6AEmSJuBqYB1gyagDaSgD\nqiNJmsa84iRJmhGSrFYqN5ZSVkxmP4OMC8iA9ydJmoZMnCRJQ5fkpCQH1cttSW5K8rGWOlcl+UiS\nw5MsBQ5tNzSuHgr36yT3Jrk2yX5J5rU51ueS3AQc1yGmJNkryTX1vs5J8rKWOtskOTvJPUnOALai\n5WpSku2T/C7J3Ul+BmzYsn79JEcluSXJnUkuSPLy/lpSkjQsJk6SpFF5E/AA8EzgXcB7k7y5pc77\ngHOBpwH71mUPJipJ1gOOAX4NbAG8FXgz8JE2x7oPeG5dp51/A94DvBfYHDgeOCrJX9bHWgP4EXAh\nsDWwD/CZ5g6SPBH4P+CHwJbAl4FPthznEOBhwF8DTwX2BO7sEJMkaZpIKQ67liQNV5KTgMeXUp7a\nKNsP2HGsLMlVwFmllNc26mwAXAU8rZRyfpL/BF5dStmsUedtwCdLKQsbx3p0KeUZXWL6I3BQKeVT\njbJfA2eUUt6ZZDfg48Cfl1Lur9fvTpUIbVXH84n6HDZvOa89gLVKKbcnOQ/4billXyRJM4ZXnCRJ\no/KrluenA09O0rxn6Kwu+9ik3q7pVOBRSf681/0keTSwHnBam31t2jjW+WNJUyPmZrybUF39oqVO\n04HAR5P8Msk+STZHkjTtmThJkqazu6bZfiatlPIV4C+AI6iG6p2Z5B2jjUqS1I2JkyRpVJ7V8vw5\nwGVlYmPIL6m3a/pr4I5Syh973Ukp5Q7gWmDbllXbAhc1jrVFkoe1xNyM9xJgm5Z9tMZHKeVPpZQv\n1cMQPwu8pddYJUmjYeIkSRqV9ZN8JslfJflH4F+BAya4j0OAJ9az5m2c5FVUkzb8Vx/xfBrYM8nr\n6pg+STXBw4H1+m9TJUlfTrJpku2pJq9o+iLVcMP96338E7BTs0I9u99Lk2yYZGvgRcDFfcQrSRoi\nfwBXkjQqRwCPBM4AlgGfK6V8ubG+05WnB8tLKdfWCcynqWbfuwU4DPjPHvbT6kBgTaqZ8p5Alczs\nWEq5oj7WXUl2pEqOzq7X70E1i95YPNckeQ3wOapE8Azgg8BXG8eZDxwM/DlwO3As1Ux+kqRpzFn1\nJElDV890d04pxYRBkjQjOFRPkiRJkrowcZIkjYLDHSRJM4pD9SRJkiSpC684SZIkSVIXJk6SJEmS\n1IWJkyRJkiR1YeIkSZIkSV2YOEmSJElSFyZOkiRJktSFiZMkSZIkdWHiJEmSJEldmDhJkiRJUhcm\nTpIkSZInqe1QAAAMjUlEQVTUhYmTJEmSJHVh4iRJkiRJXZg4SZIkSVIXJk6SJEmS1IWJkyRJkiR1\nYeIkSZIkSV2YOEmSJElSFyZOkiRJktSFiZMkSZIkdWHiJEmSJEldmDhJkiRJUhcmTpIkSZLUhYmT\nJEmSJHVh4iRJkiRJXZg4SZIkSVIXJk6SJEmS1IWJkyRJkiR1YeIkSZIkSV2YOEmSJElSFyZOkiRJ\nktSFiZMkSZIkdWHiJEmSJEldmDhJkiRJUhcmTpIkSZLUhYmTJEmSJHVh4iRJkiRJXZg4SZIkSVIX\nJk6SJEmS1IWJkyRJkiR1sWDUAYxKkvWBRaOOQ5IkSdLILSmlXD1ehZRShhXMtFElTav9AR4YdSiS\nJEmSRu9uYNPxkqe5esVpUZU0/T2wTqO42RzzB1Te2sSd1k3F8Qa1n4kea4J6Ca3T417qNwek9nKK\nrQNYJ7pNL03XS/3JnPNk6k9026ne/0S37ef1W+Vx48uk+csb26548GHmrSyfv6DxeP7KOvPm91Cn\nuR8a5as8XtYIv1Od5n7a77NZPtFt5/VUPpg6D13XPIdlHepM5tiDb4sFPbX11LbjZNph1b4ysW2r\n5xPtsxOMdXmjfJXHpVGHRvnKx2mGuqzxeHmHx73UX9GhvJf6k9lnpzoTrT/e9nPh8VT3iSnoN8ta\nXr9ljW2WNeo1/jxW2abT4TodulOdXupPtHxUx1oCfA9WpxqNZuLU3iJgvcbz1RqPF0xBeeu6qT7e\nZMons20P0uFx80NvLx+SO4UzqPrDOMZ0qD+obafb/vs6RiNxaiQ8zceZv/Ltdt5qjQ+ZzQRplceN\nD5LN8lUSp06Pmx9Cp6L+MI/Vvf5o45utxxpebNXz+T1sM69D+Uodv+doJEjzG58SF6xSTtvHGdSn\nwan4hDmq447y2NNhn9PtuPM6lDfqt46XWlbar2smUauU0768lzpTXT6qGHrl5BCSJEmS1IWJkyRJ\nkiR1YeIkSZIkSV2YOEmSJElSFyZOkiRJktSFiZMkSZIkdWHiJEmSJEldmDhJkiRJUhcmTpIkSZLU\nhYmTJEmSJHVh4iRJkiRJXZg4SZIkSVIXJk6SJEmS1IWJkyRJkiR1YeIkSZIkSV2YOEmSJElSFyZO\nkiRJktSFiZMkSZIkdWHiJEmSJEldmDhJkiRJUhcmTpIkSZLUhYmTJEmSJHVh4iRJkiRJXZg4SZIk\nSVIXJk6SJEmS1IWJkyRJkiR1YeIkSZIkSV2YOEmSJElSFyZOkiRJktSFiZN68N1RBzAzXbt41BHM\nPCfbZv1Yuvi4UYcwI120+MJRhzAjHbd46ahDmJEW/3TUEcxMi68YdQQz03eXjzqC2cnEST0wcerL\ndSYBE2bi1JfbTZz6crGJU1+OW3z7qEOYkRafMOoIZiYTp/6YOE0NEydJkiRJ6sLESZIkSZK6MHGS\nJEmSpC4WjDqA0VrCqk3QfDx/QOWtTTzM4w1qP0uBcydwrB6UDo9XNB53Gp9bOpR32rb59cCyxuNO\nzdD6dUIvTTqvTZ0HlsJtZ7fftl391n12ejzV9Se67SD3f/dSuOLsycXWz+u3yuNGB5vf6EgLVnaw\nMm9l+YoFKx9nfqMTrrJt43GjzormfhodeHnj8bxGp121vFF/6Z3cc/Yl9ams3Of8Rp1m+bxV9rO8\nbZ35HetMZtvudR66rnkOyzrU6e/Y9y69j+vOvn4S59O+fEFPbT217TjxbTv1lYdue+fS5Vxy9r1t\nt62er3ydFnQ4xqRiXd6IdXlpPG7E0PyTazxOM9Tm/wfLOzzupX6n/3ta6i+9E87+3QD32anOROuP\nt/00eLz0fjh7yRQea6r7xCT7Tbv6y1pev+bzZXW9pQXOa/x31mm3yzo87qVOL/UnWj6qYy2hNyml\n06fQ2SvJ+sAlwOqjjkWSJEnSyN0NbFpKubpThTmZOMGDydOiUcchSZIkaeSWjJc0wRxOnCRJkiSp\nV04OIUmSJEldmDhJkiRJUhcmTpIkSZLUhYmTJEmSJHVh4iRJkiRJXZg4iSQbJPlykiuT3J3ksiT7\nJFmtUWeLJN9OcnVd56Ik72qznxUty/Ik2wz/rKZeL+1W13tikmOS3JXk+iT7J5nXUmeLJCcnuSfJ\nH5L8+3DPZriSfCjJqXWb3NJm/U6N/tPanxbVdeZUf4Pu7VbXadcmr2upM2f6Ww99zfe2Nnrsa763\ndZHkBR3ey1YkeXqjXte/27kkye/btMceLXW69r+5ZAKfSexrk7Bg1AFoWtgECPAW4ArgqcCXqX4g\neOyN6unADcAbgWuA5wKHJVlWSjmksa8CvAS4uFF285RGPzpd261+E/8xcC3wbGA94BvA/cBH6jqP\nBo4HfgLsDmwOfC3JraWULw/xfIZpNeBI4HRglzbr/wc4tqXscOBhpZTmD3zPpf4G3dttzE7AcVT9\nE+C2sRVzsL91azPf29obt918b+vZqcA6LWUfB15cSjmrpbzj3+0cVKj60WGsbI87xlb20v/moF4+\ny42xr/WrlOLi8pAFeD9weZc6BwMnNJ5vAKwAthh1/NOl3YDtgAeARY2y3YFbgQX187cBS8ae12X7\nAReP+nyG0F47Abf0UG8RcB/wT42yOdvfxmu3uk1eOc62c7K/9drX6rq+t3VpN9/b+m7PBVSJ+oda\nysf9u51rC3AV8K5x1nftfy7tP8vZ1ya3zNlLmurqMUDb4RkNCzvUOSrJDUlOSbLj4EOb1lrb7dnA\nBWXVqyTHU7XdUxp1Ti6lLGups3GShVMZ7AyyE3AX8H9t1s3l/tbJ55PclOTXSXZuWWd/6873tu58\nb+vPq4DHAl9vs268v9u56ANJliQ5O8n7k8xvrOul/6nzZzn7Wp8cqqeHSPIk4F+B945T57nA64Dt\nG8V31tucSvWNxmuBHyR5VSnl6KmLeHro0G7rUH272HRDY9159b9XjlNn6WAjnZF2Ab5VSrmvUTan\n+9s4PgqcCNwNvBQ4JMkapZSD6/X2t3H43tYz39v6swtwfCnl2pbybn+3c81/A2dTfeh/LvBJqj7z\n/np9L/1vThvns5x9bRJMnGaxJPsBe45TpQCbllIubWzzZ1T3lnynlPLVDvt9KvADYJ9Sys8e3Fkp\nNwMHNKqelWRd4N+BGfPhYqrabbbrp9163O9zqMZuv3GVnc3h/jaeUsp/Np6el2QNqjaZNf8pTmFf\n872tj3abiybx/8TLqBLvVSv7dwuNNiulNP/eLkxyP3Bokg+WUh6Y0kCnmUF/JpkLfW0qmTjNbp8B\nvtalzoPfBiZZj+pbiF+WUnZvVznJZsAJwBdLKfv1EMMZwN/2Fu60Mch2ux54ZkvZ2o11Y/+u3aXO\nTDChdpuAXYFzSynn9lB31ve3PpwBfDTJavUHjtnQ3wbeZr63PajXdptL723t9NOWu1Dd8/WjHvbf\n+nc7G0ym/51B9Zl1Q+Ayeut/s8XAP8u1mI19bcqYOM1i9bekPc36VH87cSJwJh1m60ryFOBnwNdK\nKXv1GMZWwHU91p0WBtxupwMfSrKoMRb7pVRDVC5u1Pl4kvmllOWNOr8rpcyYoSwTabde1d+E/QPj\nf9vWNKv7W5+2Am5t/Ic44/vboNvM97a+zJn3tnb6bMt/AQ5vtMV4Wv9uZ7xJ9r+tqIbJ3lg/76X/\nzQqD/izXxqzra1PJxElj3078nGoWmz2AJyTVDJWllBvqOk+l+mM8Fjggydg3O8vH3rSSvIlqKtBz\n6nWvofqP4s3DOI9h66XdqKbhvRj4RpI9gXWBfYGDG29S3wb2Ar6a5FNUU/a+C3j3cM5k+JI8keoG\n6Q2A+Um2rFddXkq5q1H1DcB84Ftt9jGn+ht0b7ckO1B96/or4F6qDxIfBPZv7GZO9bce2sz3tjZ6\n+Bv1vW0CkryE6mrJV9qs6+Xvds5I8mzgWcBJVFOQPxf4LPCNRsLdS/+bU3r8LGdfm6xRT+vnMvqF\nasay5S3LCqoPDmN19m5TZzlwZaPOm4CLqN7obqX6RujVoz6/UbZbXe+JVPdB3El18+qngHktdZ4K\n/ILqZs2rgfeP+vymuO2+1qE/Pb+l3qnAER32Maf6Wy/tRnX/xNlU37reXj/etc1+5kx/66HNfG/r\no93qOr639d6e36KaYbDdup7+bufKQnUF5HSqiSHuAi6kSgRWa6nXtf/NpYXePsvZ1ya5pG5ISZIk\nSVIH/o6TJEmSJHVh4iRJkiRJXZg4SZIkSVIXJk6SJEmS1IWJkyRJkiR1YeIkSZIkSV2YOEmSJElS\nFyZOkiRJktSFiZMkSZIkdWHiJEmSJEldmDhJkiRJUhf/H46/gljzB1oTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xcc14588>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\bogda\\Anaconda2\\lib\\site-packages\\ipykernel\\__main__.py:108: RuntimeWarning: divide by zero encountered in log\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOkAAAFzCAYAAADbvS5yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJztnXm4XEWZ/z/f3LAjoCCbQBAJiGwzREFcWERZhAH3MYPD\nNi4IqD+HkTgMGgyOgcgoKqAgDohiBAd+oBgIu7IJDkEJIWFNAhgSSIAACUnIve/8UdXkpNN9eknf\nrtP3vJ/nOc+9p05Vnfd09/dU1VubzAzHcYrLsNQGOI6Tj4vUcQqOi9RxCo6L1HEKjovUcQqOi9Rx\nCo6L1HEKjovUcQqOi9RxCk5pRCpprqTzm4z7J0mTBtum1MTP5IoO5neQpAFJH24i7q8lTc+crxXT\nnpIJOz6GbdrCvfds/wmKSdMijR9Ao6Nf0j6DafBqMAC8PgZS0q6SxkraskZcy8YdwgzGMzabpxG+\nk0ZxVspP0pckHbma9+4phrcQ9zNV50cDH4zhyoRPp5iMAPoz57sBY4HrgDlVcd/PEP3Cu4AaRwFW\n/d3U4kLgYjNblgn7MvAocFk2oplNlrROVdwhQdMiNbNfZc8l7Q180MwmNpNe0tpmtqRF+zqGmb1W\nFSTqCNHMlg++RZ0n9WfcCmbW30QcA5oW3VAUKAxSmzTTPviopLMk/Q14RdKakjaR9H1JD0p6RdKL\nkn4n6R118jhc0umS/iZpsaTJkkZUxX27pKtjG+tVSU9K+qWkdTJxXm+TSvoCcGm89KdMVX3PeH2V\nNqmkzSRdIunZeI/7JY2uirNjzOuEeDwe494lafcmPrdKG2wvST+T9Hz8fH4m6Q1VcedKukLSoZLu\nk7QEOCpeGy5pnKQnJC2Nf0+XVPOlLOnDkv4abZ0q6bCq6019ZxED1pD03WjjK5KulLRFVZ4rtUkb\nfB6bxvNngO2AgzNNrEnxWiVsz6o83ivpRkkLoy0314izoaRzJc2StETSPEnXS9o5z75u0Up1tx3O\nABYBZwHrEaqbOwIHA/8DzAa2AI4HbpP0DjObX5XHWGApcCawMXAKcAmwP4TSA7iR0L75PvAssDVw\nOLA+8GrMJ1tq3gz8ON73m8ATMfyxGnGRtB5wB/AW4IfA08A/ApdJWt/Mflpl878AawPnAX3AGOB/\nJO1g+XMDK9cuBJ4DTgN2jna+hfC5ZePuBvwcOB/4CTAtXvsl8ElgInAn8N74nCOB6vbcLsAvoq3P\nA58FrpL0ATO7I8Zp5TsTMI5QAv5ntPv/AdtLGpWppTTT7q+Oc0J81rnAhHivOXXiIulg4Brg7vj8\nxOe7TdK7zeyBGPbf8fl+CDwCbALsE597Gqkxs7YO4EdAf51rBxFEMw0YXnVtzRrxtycI8eQaeUwB\n+jLhXyOIfbt4vleMd0gDe58Bzs+cHxnz2bNG3LuBSZnzMTHuRzJhw4H/BRYAa8ewHaMtfwPWy8T9\nZEz/gQY2fiGmvx0Ylgk/Lab/YNXz9APvq8pjz5jHOVXhP4zx96qRx0GZsI0IL7o7VuM7e6zymcTw\nz8Twz2bCJgIPZc7XinFOqfo8+oFNM2GPZr+bqnu//n0Saokzgauq4q0LPAlcnQlbBExoVwuDfQx2\nF8x/W1X7zjLtBkl9kt4EvBg/0D1q5HGRrdx+uZ3wBt0unr8Y/x4iaa2OWb4yhwCzzezqSkB8rh8R\nftTvqYp/mZktypxX25yHAT8xs6zn89yYvrprY7qtKO0qfDjm8f2q8P+KeRxaFT7TzCa/fnOzFwlO\nmb0lbRjDWv3OLraV28a/IrzMGnbNdJA9Cc7CiZI2rhwEkd5KrIlFXiI872ZdtK9pBluks6oDJA2T\ndIqkxwlv4vmEN/dIYMMaeTxVdf5C/PtGADN7mFBVOxFYIGlSbMus35lHAMKX/UiN8OmEH/6IqvBc\nm5vgsexJFM5zwLZV8WbWSDsCWGZms6vymE2o+lfb+miNPCrPOgLa+s6q7R8gNCmq7R9MRsa/lxM+\nu8rxLKFkXy/zUv834J3A05LulvSNar9HSga7TfpqjbBxwKmENtSthB/wAKGNWOulUc8L+Lr73sy+\nJOmnhHbogQTRnhLbHc+2b37bNLS5Q9T6fAeDVr+zIjCMUKP4MvW7BZcBmNllkm4FPgp8iNC8GSPp\nH8zs1m4Ym8dgi7QWHye0KU7IBsYq1OPtZmrBCfAA8G1J+wG3EJwE36mXpIXsZwM71AjfKeYzu8a1\n1WEkcE/lRNJGwJupUTOpwWxgTUkjsqWppG2AdVjV1pGsyo6ZvKD172xkVbxhhKr+7U3Y34hmv7fH\nCS/FhWZ2S8NMzeYQXu7nxWrvX4F/J7yUkjKYb8F6H2Y/VSWKpH8meG6bzSObdoP4I8gyNf7Na6Mu\ninZs1OgewCRghKQjMvcdDpxEaJvd2UQezSLg+Kpn+hLhs7iuSVtF8KhmOTnm8fuq8LdKOuT1m0tv\nBP4JuNvMFsbgVr4zgGMkrZs5PzLG7cRQy0U09539idDsOEWZrrgKkjaJf4dXN43MbB4wj/zfT9cY\nzJK0XtXuWuBrki4E/gzsTujOmNVCHlkOASZI+g2hfbUWYTTUEuCqnHT3E360p0nanNDWusHMXqgR\n9zxCqfwrSecSvvxPE5wmx5vZ0ibsbIX1gRslXUXoIvk8cJOZ3dgooZndK+ly4MuS3syKLph/Aiaa\n2T1VSWYAv1ToQ14Q77URoRSp0Mp3BvAK8EdJlwJbEaqc0wjdRavLfcBRkr5OaJM/Y2Z/jNeyTaDl\nkj5H6IKZGm2ZE+35IMED/4+El8cj8fczFVhM6I7ZhdDlk5523cIEz+byOtcq7vAP17i2NsHz+Dfg\nZUK19O+Bu4DfN8qDUBXrBz6V6Qr4GcFZsYjgGJjMql0Tc4DzqsKOJ1SLlrGy+/7urC0xbDPg4pj/\nq4SuoU/Xse2EqvC1YvjXGnymlS6HvYCLCKJ5MT7fG2o8z+V18hkOnE5w1iyJf8eS6crK5kF40T0Q\nn+sB4LDV/M6OIPSNz43xrwS2qMpzIjAt7zOidhfMloTawEvx2qSqe+9ZdZ+/J7ysnyMI8HFCP/L7\nMs/2XeAv8bNeSOhaO7ZdbXT6UDTUKQAKI6HOB3Y1s4dS2+MUg6J65hynI0g6UdLMOOTxT5Le1SD+\nforDLCU9IunobtlaDxepM2SR9I+EQRxjCdXevwKTK06jGvG3JbS/bya0u38AXCTpQ92wtx4uUmco\n81XgAjO71MxmEHwQi4Hj6sT/IvCEmZ1iZg+b2XmE8cpf7Y65tXGRFggzu8DM+rw9uvpIWgMYRSgV\ngdenvt0E7F0n2bvj9SyTc+J3hRSDGXqOOObzIEKXQ0/M1+wQaxOG8k02swWDcYM4yKJm9bMJ5pvZ\nk3WubUKYgTSvKnweKwZrVLN5nfgbSFrLOt/V1hQu0uY4iKqVAErGkYRB8h1F0jZrwOzq2fgtsDRO\n/6sn1CGBi7Q5ZgF8jPZf+b3IfF4fDTJrkG6xyWu097lG29aKSWuJdD6h37R6ZstmhP7bWsytE/+l\nVKUouEibZQmEX8MWDSIOUQa1ir85YYRCKzT64ZrZa5LuAw4AfgsgSfH8h3WS3U0Y2JHlwBieDBep\nk5zhwBptpGmC7wGXRLHeS/DSrktY2QNJ44EtzazSF/oT4ERJZxFWazgA+ATdnQe7Ci5SJzl9tP5D\n7GsijpldEftExxGqrX8hrELxXIyyOWGpnUr8WZIOJQyB/DJhmZx/MbNqj29XcZE6yRnEkhQzO58w\n1LLWtWNrhP2R0HVTGLyf1HEKjpekTnIGq7o7VHCROskZzOruUKBMz+oUFC9J83GROsnxkjSfMj2r\nU1CG0/oPsUw/XPfuOk7BKdMLySkoXt3Np0zP6hQUF2k+ZXpWp6C4dzcfF6mTHC9J83HHkeMUnDK9\nkJyC4tXdfFykTnK8uptPmZ7VKShekubjInWS4yVpPmV6Vqeg+LDAfNy76zgFp0wvJKegeHU3nzI9\nq1NQ3HGUT6mru61ui+cMDpWStJWjTKVLaUXa6rZ4zuDhIs2ntCKl9W3xnEGiUt1t5fDq7hCnzW3x\nHCcJZao1ZGlnWzxnkBjeB2uoxTRG2I6pBJRVpG1xPWHDziy7ALsmsKXTTAUerArr1kasfX0wvMU6\nXd8ALtIhTjvb4nEwQ3dXtV1Z9WXzDHBhF+49fBis0WIjs0w/3FK2Sc3sNaCyLR6w0rZ4d6Wyq6wM\nH97eURZK9KirkLstntM9hvfBGi3+Esv0wy3Ts65EE9viOU4hKK1IIX9bPKeLDKP1js+BwTCkmJSy\nTeoUjMSjGSSdKulOSYskPd9kmoslDVQdkzpn1QpKXZI6BaGdCaWdLUnXAK4A7qa1EWfXAccAlV7e\npR21KuIiddLTzjSYDvaRmtm3ACQd3WLSpd3wYXh113HaZz9J8yTNkHS+pDcNxk28JHXS047jKH3x\nch1wJTATeBswHpgkae84DrxjuEid9DSo7k5cCBNfWjlsYYPqrqTxwJicKAbsZGaPNGdkVWKzKzKn\n0yRNBR4H9gNubSfPerhInfQ0cByN3jgcWaa8CqMey831bODiBnd+ohnzmsHMZkqaD2yPi9QZcgxC\nddfMFgAL2rSoZSRtBWxMGPLcUdLX7B0nfT/p1pJ2B0YAfZJ2j8d6mTgzJB0R/19P0gRJe0kaIekA\n4GrgEWBy5ywLeEnqOGFo6FGZ8ynx7/7AH+P/I4EN4//9wG4xzUbAHII4vxknb3QUF6mTnsSrY5vZ\nscCxDeL0Zf5fQpi52BVcpE56erMLpmu4SJ30+MK7ubhInfS4SHNxkTrp8R2bcilRzd5xepMSvY+c\nwuKOo1xcpE56vE2ai4vUSY+LNBcXqZOePloXnYvUcbqIl6S5lKj57Ti9iZekTnq8JM3FReqkx9uk\nubhIW+CQ+zZl5z3WTG1G15g2ZRkXjnp28G/kJWkuLlInPS7SXNxx5DgFx0tSJz1ekubiInXS446j\nXFykTnq8JM3FReqkx0Wai4vUSY9Xd3Nx767jFBwvSZ30eHU3Fxepkx4XaS4uUic9LtJcXKROetxx\nlIuL1EmPl6S5uHfXcQqOl6ROerwkzcVF6qTH26S5uEid9HhJmou3SZ30pN/p+xpJsyW9KmmOpEsl\nbdFEunEx/mJJN0ravnNWrcBF6qSnss1EK0dnf7m3AJ8EdgA+BrwN+E1eAkljgJOAzwN7AouAyZI6\nvr6OV3ed0mNmP8icPiXpTOD/S+ozs/46yb4CnGFm1wJIOgqYB3wEuKKT9nlJ6qSn1apuO1slNomk\nNwFHAnfWE6iktwKbAzdXwszsJeAeYO9O2+QiddKTuE0KIOlMSa8A84GtCSViPTYHjFByZpkXr3UU\nr+466Wng3Z14RziyLFycn6Wk8cCYnCgG7GRmj8TzCcBFwAhgLPAL4LD8u3QHF6mTngb7k47eNxxZ\npjwOo76Wm+vZwMUN7vxE5R8zex54HnhM0gxC23QvM7unRrq5gIDNWLk03Qy4v8E9W8ZF6gxJzGwB\nsKDN5JVXxlp18p4paS5wAPAAgKQNgL2A89q8Z11cpE562nEEdeiXK2lP4F3AHcALwPbAOOBR4O5M\nvBnAGDO7JgadA5wm6TFgFnAG8DRwDR2mlI4jSe+X9FtJf5M0IOnw1DaVmrSOo8WEvtGbgBnAT4G/\nAPuZ2WuZeCOBDSsnZjYB+BFwAcGruw5wiJkt65hlkbKWpOsRvoifAVcltsVp0Catm6YDmNmDhGpr\no3irWGhmpwOnd8aS+pRSpGZ2PXA9gCQlNsfxsbu5lFKkTsFI2CbtBUrZJnWcXqJE76PV5ztffZE3\nbLjye+2w0ety2Oh1E1nUOa6duJhrJ648QuDlhQPduXnCNmkv4CJtgVO/v9GQ3US41stm2pRlfNQ3\nEU6Oi9RJj7dJcynRo65A0nqETuuKZ3c7SbsDz5vZU+ksKylekuZSSpEC7wRuJQyyNuC/YvjPgeNS\nGVVavE2aSylFamZ/oFRfs9PLlFKkTsHw6m4uLlInPe44yqVEj+oUFm+T5uIiddLj1d1cXKROelyk\nuZSo0uA4vYmXpE563HGUS4ke1SkqNgxWnVLdOE1ZcJE6yenvg/4Wf4n9JWqTukid5Ay0IdKBEom0\nRJUGx+lNvCR1ktPfJ5b3tbbUVH9fZW7E0MdF6iSnv6+P/uGtVer6+waA5YNjUMFwkTrJGejro7+v\nNZEO9AkXqeN0iX6G0d/iEKJ6m4YORVykTnL66WO5i7Qu7t11nILjJamTnAH66G/xp9ilxUYLgYvU\nSU57bdLyyNRF6iQnlKStiXTAReo43WOgjZJ0oESuIxepk5zlDGvZu7u8RD7P8jyp4zRA0pqS/hI3\nlt6tQdyLY7zsMWkw7PKS1EnOAMPb8O4OSnV3AvA0sGuT8a8DjmHFTghLB8EmF6mTnvbapJ2tBEo6\nBPgQ8HHgw00mW2pmz3XUkBq4SFvgulHP8pfURnSRZ7p0n/a6YDonUkmbARcChwOvtpB0P0nzgBeA\nW4DTzOz5jhkWcZE6yWlvWGBHZ31fDJxvZvdLGtFkmuuAK4GZwNuA8cAkSXubWUfn0LlIneQ0GnF0\nw8QXuHHiCyuFvbIwv00qaTwwJieKATsBBwPrA2dVkja2GMzsiszpNElTgceB/QibgXUMF6lTeA4c\n/UYOHP3GlcIenrKYY0Y9kpfsbEIJmcdMYH9gb2CptJI+/1fSZWZ2bDM2mtlMSfMJW2q6SJ2hRX8b\nI44axTezBcCCRvlI+hLwH5mgLYHJwKeAe5u1R9JWwMYMQlPeReokJ6V318yezp5LWkSo8j5hZnMy\n4TOAMWZ2TdyEeiyhTTqXUHqeBTxCEHhHcZE6yUnt3a1BLcfPSGDD128PuwFHARsBcwji/KaZvdZp\nY1ykTnIK4N19HTObTY2dZsxWLN9tZksIDqeu4MMCHafgeEnqJKe9Sd/lWR3bReokp4Bt0kLhInWS\n096kby9JHadr9Lcxn7RMJWnPPKmkoyUdmjmfIOlFSXe1MN7SKSD9sU3a2lGekrRnRAqcSpyhIGlv\n4ETgFGA+8P2EdjnOoNJL1d2tgcfi/x8BrjSzCyXdCdyWzCpntfE2aT69VJK+QhgbCXAgcGP8fwmw\nThKLnI5Q8e62dvTST3f16KWS9EbgIkn3AzsAlfVkdgZmpTLKWX2KNOKoiPSSSE8Evk2o9n48znIA\nGAVMTGaVs9r4YIZ8ekakZvYicFKN8LEJzHE6iA9myKennlTS+yX9Mna7vCWG/bOk96W2zXEGi54R\nqaSPE6YDvQrsAawVL21I6J5xepSBlp1GfaWq7vaMSIHTgOPN7HNAds7enQTROj3KQBve3U4v6Vlk\neqZNCuwI/LFG+ELCxFunR1nehne31fi9TC+JtLJMxayq8PcBT3TdGqdjuHc3n14S6U+BH0g6jrC8\nxZZxeODZwBlJLXNWC/fu5tNLIj2T0Ia+GViXUPVdCpxtZj9qJSNJ/w58FHg7wRF1F2GRqdw1Ih0n\nBT0j0rgq+H9K+i6h2rs+8JCZvdJGdu8HfgT8L+EzGA/cIGknM2tlmwGnA/jY3Xx6RqQVzGwZ8NBq\n5rHShjySjgGeJYxeumN18nZax+eT5tMzIo1rnX4dOADYlKruIzPbbjWy34jQzu34ZjtOY/rbcBz5\n2N1ichGwL/ALwirhHdkUR2FvgXOAO8xstUpopz28uptPL4n0EOBQM7uzw/meD7wDeG+jiNcDa1eF\n7ULzO84WmanAg1VhS1IY4qxCL4n0BTpcHZV0LmHD2PebWcM9PA4GtuikAQViV1Z92TxD2LRzsCnC\nJsJFppee9BvAOEnrdiKzKNAjgP3N7MlO5Om0h0/6zqeXStKTCZu1zpM0i5XH72JmTY/flXQ+MJqw\ns/OiuNMzwMK4hYDTRXzSdz69JNKrO5jX8QTH021V4ccCl3bwPk4T+LDAfHpGpGb2rQ7mVZ66Ug/g\nwwLz6aknlbSRpM9KGi/pTTFsj8oEcMcZivRMSSppN+AmwtS0bQkD7p8HPgZsQ9gr0ulBvJ80n14q\nSb8HXGJmI1m5C28SsE8ak5xOUBkW2MrR6equpEMl/UnSYknPS7qqiTTjJM2JaW6UtH1HjYr0kkjf\nBVxQI/xvwOZdtsXpIK1vMTG8ZUdTHnFpnkuBnxG6i98D/KpBmjGEhfE+D+wJLAImS1qzY4ZFeqa6\nS5iWtkGN8B2A57psi9NBUg5mkNRHGBZ6spldkrk0o0HSrwBnmNm1MZ+jgHmE3RWu6IhxkV4qSX8L\nfFPSGvHcJG0DnAVcmc4sZ3VJPJhhD2BLAElTYvV1kqSd6yWQ9FZC7e3mSpiZvQTcA+zdKcMq9JJI\nTybMIX2WsK3EHwh7w7wC/EdCu5zeZjtAwFhgHHAoYQjqbZLqrZ21OaGffV5V+DwGoenVM9VdM1sI\nfEjSe4HdCYKdYmY3pbXMWV0aeXdnTbyH2RPvWSls2cL8ufmSxgNjcqIYsBMrCqpvm9nVMe2xwNPA\nJwm9CEkpvEglrQMcUKn7A4exYs3dD0s6EPimD+frXRpN+t5q9HvYavR7Vgp7Ycosbhx1el62ZwMX\nN7j1E8SqLjC9EmhmyyQ9Qejaq8VcQum7GSuXppsB9ze4Z8sUXqTA0YQqSEWkJwHTiHuVEtYpegbf\no7RnGYxJ33GvoAW5kQBJ9xGckjsS1roi+j22BWbXyXumpLmEBQgeiGk2APYCzmv2GZqlF9qkR7Lq\njKl/MrP9zWx/4GvAp7pvltMpUq5gb2YvAz8BviXpQ5J2AH5MqA7/phJP0gxJR2SSngOcJukfJO1K\n6MJ5GrimI4Zl6IWSdHvCnOQKS4CBzPm9DMLbyykV/0aYVXUpwSl5D/CB6AepMJKwpQkAZjYhTpu8\ngLD8zu3AIXENro7SCyLdiBVtUMzszVXXh2WvO71H6gH2ZtYPnBKPenFWMdDMTgdO75ghdegFkT5N\nWKXk4TrXd4txnB7F55Pm0wtt0kmEFRmqlxeqeH7HAr/vulVOx6jMJ23lKNMA+14oSb9DcAw9HJc8\nqawyvyPB0zs8xnF6lNTV3aJTeJGa2TxJ7yF43M4k9E9B8L7dCJxgZtUjP5wewqeq5VN4kULolwIO\njhO9K9OBHjMzX8zaGfL0hEgrRFHem9oOp7P4NhP59JRInaGJbzORj4vUSY63SfNxkTrJ8RXs83GR\nOsnpb6MkLVN1tzyvI8fpUbwkdZLj3t18XKROcty7m4+L1EmOe3fzcZE6yXHvbj4uUic5yxlGX4si\nXV4ikZbnSR2nR/GS1EnOQBvbRgyU6Kdbnid1Cou3SfNxkTrJ6WcYw7yftC4uUic5AwN99A+0WJK2\nGL+XKc/ryHF6FC9JneT09w+D5S1Wd/vLU764SJ3k9C/vg+UtDgtsUdS9jIvUSc5Af1/LJelAv4vU\ncbpGf/8wrGWRenXXcbpG//I+Bl5rTaStirqXKc/ryHF6FC9JneTYQB/W3+JPsUT9pC5SJz3LW++C\nYXl5KoHleVKnuFS8u60cHfLuStpX0oCk/vg3e4zKSXdxjfiTOmJUFV6SOunpFyxX43jVaTrDncDm\nVWHfJmwifF+DtNcBx7Bif6KlnTIqi4vUSU8/sLyNNB3AzJYDz1bOJQ0HjgB+0ETypWb2XGcsqY9X\ndx1nZY4A3gRc0kTc/STNkzRD0vlxQ7GO4yWpk56EJWkNjgMmm9mcBvGuA64EZgJvA8YDkyTtbWbW\nSYNcpE56ltO6SBvElzQeGJMTxYCdzKyyKTWS3gIcBHyi0e3N7IrM6TRJU4HHgf2AWxulbwUXqZOe\n5cBrOddvmBiOLK8sbJTr2cDFDeI8UXV+HDAf+F2jzKsxs5mS5hP2z3WROkOMAfKrrweMDkeWh6fA\ncXV7SDCzBcCCFi05Bvi5mbVcmZa0FbAx8EyraRtRSseRpOMl/VXSwnjcJeng1HaVlkqbtJWjw21S\nSQcA2wI/q3N9hqQj4v/rSZogaS9JI2Laq4FHgMmdtaykIgWeIrRX9gBGAbcA10jaKalVTkqOA+7M\ntlGrGAlsGP/vB3YDrgEeBn4K/BnYx8zyKu5tUcrqrpn9viroNElfBN4NTE9gUrkZBMdRq5jZkQ2u\n92X+XwJ0reZVSpFmkTQM+BSwLnB3YnPKSbG6YApHaUUqaReCKNcGXgY+amYz0lpVUlykuZRWpMAM\nYHdCO+MTwKWS9skT6vUERWfZBdh10EzsHlOBB6vClnTr5i7SXEor0jhms9JPdr+kPYGvAF+sl+Zg\nYIsu2JaCXVn1ZfMMcGECW5yVKa1IazAMWCu1EaXES9JcSilSSd8hjL18EngDcCSwL3BgSrtKS6MR\nR/XSlIRSihTYFPg5ofa6EHgAONDMbklqVVnpp/WS0UvSoY2ZfTa1DU4Gr+7mUkqROgXDRZpLWYcF\nOk7P4CWpkx4vSXNxkTrpKcDY3SLjInXS4yVpLi5SJz0u0lxcpE56fDBDLu7ddZyC4yWpkx4fcZSL\ni9RJj7dJc3GROulxkebiInXS4yLNxUXqpMe9u7m4d9dxCo6XpE563Lubi4vUSY+3SXNxkTrpcZHm\n4iJ10uOOo1zcceQ4BcdLUic97jjKxUXqpMfbpLm4SJ30uEhzcZE66XHHUS7uOHLSM8CKdmmzx0Dn\nbi9pD0k3SHpB0nOSLpC0XhPpxkmaI2mxpBslbd85q1bgInVKjaQtgBuBR4A9Cfty7Qxc0iDdGOAk\n4PMx3SJgsqQ1O22jV3db4EI+z9DdV60WXdpXLe1qgYcBy8zspEqApOOBByRtZ2ZP1En3FeAMM7s2\npjkKmAd8BLiiY9bhJalTBCqOo1aOzjmO1gKWVYVVtmZ9X60Ekt4KbA7cXAkzs5eAe4C9O2ZZxEXq\npKfiOGrl6FxJeguwuaR/k7SGpDcC4wGjfrVp83h9XlX4vHito3h110lPxXFUj3kTw5Fl+cLcLCWN\nB8bkRDFgJzN7SNLRwPcI4lwO/BB4lo66p9rHReqkp1E/6cajw5HllSlw/6i8XM8GLm5w5ycAzOzX\nwK8lvZngAAI4mRU7wVczFxCwGSuXppsB9ze4Z8u4SJ0hiZktABa0mOY5AEnHAa8SvL614s2UNBc4\ngLC3LZIDKuy3AAAI+0lEQVQ2APYCzlsNs2viInXSk3gvGEknAncBrxB2e58AnBKdQZU4M4AxZnZN\nDDoHOE3SY8As4AzgaeAaOoyL1ElP+hFHewKnA+sDM4DPmdmvquKMBDasnJjZBEnrAhcAGwG3A4eY\nWbWneLVxkTrpaeQ4qpemQ5jZ0U3E6asRdjpB3IOKi9RJjw+wz8VF6qTH9yfNxQczOE7B8ZLUSU96\nx1GhcZE66UnsOCo6LlInPe44ysVF6qTHRZqLO44cp+B4Seqkpx0nkDuOHKeL9BPmlLSapiS4SJ30\ntCM4F6njdJF+whTsVvAuGMfpIstpvbrbqqh7GPfuOk7B8ZLUSU87jqMSlaQuUqcYlEh0reLVXccp\nOC5SQNLXJQ1I+l5qWxynmtKLVNK7CPt5/DW1LY5Ti1KLVNL6wC+BzwIvJjbHcWpSapES1kj9nZnd\nktqQcpN2n4miU1rvrqRPA38HvDO1LY4vcpRHKUUqaSvC4sYfNLMWFu64Hli7KmwXYNeO2ZaOqcCD\nVWFLakUcBHz9lDxKKVJgFPBmYIqkSjd6H7CPpJOAtcysRs/dwQzd/Ul3ZdWXTZf2J/VZ37mUVaQ3\nseov8hJgOnBmbYE6g4eXpHmUUqRmtgh4KBsmaRGwwMymp7HKcWpTSpHWwUvPZHhJmoeLNGJmH0ht\nQ3nxNmkeLlKnAHhJmoeL1CkAXpLmUfYRR45TeFykTgFIOyxQ0qmS7pS0SNLzdeJsLen3Mc5cSRMk\n5epH0m1xdlXl6Jd0fqv2eXXXKQDJhwWuAVwB3A0cV30xinESMAd4N7Al8AtgGXBaTr5GGA3yDVas\nPbG4VeNcpE4BSOs4MrNvAUiqt+P3QcDbgf3NbD4wVdI3gDMlnW5mecYsNrPnVsc+r+46BaDiOGrl\n6Krj6N3A1CjQCpOBDYGdG6Q9UtJzkqZK+o6kdVq9uZekTgEofBfM5sC8qrB5mWv1Fgy4DJhNqCbv\nBkwAdgA+0crNXaROD3AzUD3l95XcFJLGA2Nyohiwk5k9snq25dzA7KLM6TRJzwA3S3qrmc1sNh8X\nqVMAGvWT7huPLI8CJ+RlejZwcYMbP9HQtMBc4F1VYZtlrjXLvQQH0vaAi9TpJTpf3TWzBcCCdi2q\n4m7gVEmbZNqlBwILqZqo0YC/J5Tgz7Ryc3ccOQUgreMo9oHuDowA+iTtHo/1YpQbCGL8haTdJB0E\nnAGcW1k0QNKWkqZLemc8307SaZL2kDRC0uHAz4E/mFn17PpcvCR1CkByx9E44KjM+ZT4d3/gj2Y2\nIOkw4MfAXcAiwvzjsZk0axCcQuvG82XAB4GvAOsBTwG/Af6zVeNcpE4BSN5PeixwbIM4TwGH5Vyf\nTVjdo3L+NLBfJ+zz6q7jFBwvSZ0CkHxYYKFxkToFIHmbtNC4SJ0C4PNJ83CROgXAS9I8XKROAfCS\nNA/37jpOwfGS1CkAXt3Nw0XqFACv7ubhInUKgJekebhImyNupTY/P9aQ4/Xnrd5KrsPMpXXRlee7\ncJE2x7bhz1VJjUjItoSB5Z1mPrAYrlq3YczaLKYEapVvINYYSRsTFqOaRfc27SwCaxMEOjnOz+w4\nkrYBNmkz+Xwze7KT9hQRF6njFBzvJ3WcguMidZyC4yJ1nILjInWcguMidZyC4yLtEnFXrcNzro+I\ncXbrgi1H19s9zCkepROppM0k/UjS45KWSJot6beSPpDYtCcJWxa0tNxjIyTNlPTlquBfE1a2c3qA\nUo04kjSCMHLmeeBkgiDWAA4GzgXekco2Cx3Wz3bpXkuBpd24l9MBzKw0B2GPySeBtWtc2yD+3Rq4\nBniZsEL55cCmmXhjgfsJS0DOjvHOJdRKTiGsTj4POLUq/wHg+GjDYuBx4OOZ6yNinN3i+b7x/APA\nnwlrvd4J7JBJsx1wNWHw68uEbQwOyFy/NebRX/kbw48BXqiy74vAYwTxTgc+U8P+fyGMjVwEPAL8\nQ+b6RoQNip6Nz/cwcHTq73woHMkN6NqDwhvjj/WUnDiKAvwD8HeE/T/+DNySiTMWeCmK9+3AoYSh\ngtcB5wAjowgGgHdl0g3EH/CxhL1AxhGmfuwYr4+I9lWL9C7gffFefwBuz+S5G/A5YCfgbcC3ooC2\nyjzzk8CpwKaVlw1wNPB8Jp+PRnF+Idr21WjbvlX2zwY+FV8O58TPYaN4/VzgPsJWCtvEl8uhqb/3\noXAkN6BrDxoENwAckRPnQ4SVx7fMhO0U042K52MJpda6mTjXAY9X5TU9+0KIeZxbFefuShi1S9J+\nYL9M/ENi2Jo5zzAVOCFzPhP4clWcapHeAfy4Ks7lwO+q7D89c75uDDswnl8DXJT6ex6KR5kcR2oc\nhbcDT5nZnEqAmU0HXiSItcIsM8tuqz6PVTfumUcovbL8qer87qp8azE1839lo59NASStJ+lsSQ9J\nekHSy/EZtmmQZzU7seoslztr2Pa6LfH5X2LFM/4YGC3pfklnSdq7RRucOpRJpI8SdrR6ewfyqp6h\nbHXCOvH5ZvOtzIao5PtfwBHA1wlV4t0JzrA1O3DfRrZU7BkGYGbXE14O3wO2AG6SNGGQ7CgVpRGp\nmb1A2EL9xFpbokvakFBF3VrSWzLh7yA4RaZ1wIx31zifvhr5vQe4xMx+a2bTCG3ebaviLCOzR0kd\npgPvrQp7L61t64eZLTCzX5jZUYR27edbSe/UplRdMMCJhPbXvZLGAg8QPoMDgS+Y2c6SHgQuk/RV\nQvfMecCtZnZ/B+7/SUn3RRs+Q2gn520UVKuKng17FPiYpGvj+bgaaWYB+0i6HFhqteeFfhe4XNJf\ngJuAwwnOpAPyHydjlPQtguNoGmEe6mG0KHKnNqUpSQEsbIG+B6Fr4mxCG+sGgkj/NUY7HHiB4Em9\ngdAt8el2blfjfGzM668EkX7azB5ukCYv33+Ntt5JcNxcz4pt+yp8k1C6Pk6dflgzu4awRV+l7/hz\nwDFmdnsTtlTClwHfITzbbYT1UEbXup/TGj7p23EKTqlKUsfpRVykjlNwXKSOU3BcpI5TcFykjlNw\nXKSOU3BcpI5TcFykjlNwXKSOU3BcpI5TcFykjlNw/g/QGLp/fgS7hgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xde58da0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAGHCAYAAAAOf30gAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XmcHFW5//HPN0MEBAwgEMDLchHZZBEiCqgg+iMI4gK4\n3AiKiF5AFC/IIqgQEGQRRWRzAYJyb6K4AVeBsBuBKBcDyBKQJQkIEshCWBIgmXl+f5xqqFR6Zmp6\nuqd6Zr7v16teM111qurp6up+6pw6VaWIwMzMrJ2NqDoAMzOz3jhZmZlZ23OyMjOztudkZWZmbc/J\nyszM2p6TlZmZtT0nKzMza3tOVmZm1vacrMzMrO0NmmQlaQNJXZKOrDqWgSBpfPZ+V2/iMmdKuqpE\nuV2yde+cG3eppBmFcl2STujDui/pe9SN60t8JZe3zDaw1+W+o5+rOpbuSLpF0s1Vx9Ff9b6jg4mk\nz2fxr192nj4lK0kHZCvokrRTN2WeyKb3+qNoPYpsaPYyGy1bL56lxknaUdKJkt5UZ3ldfVx/MzR7\nG7biMwFA0jrZttu6ZPmetrV1L0j74qAg6VBJB3QzuaXfJ0krZvtYKxJin79LyzW4okXAZ4Db8yMl\n7QK8BXi5weVaG4iIP0laMSJe7aXoisCS3OudgBOACcDzhbKbMvA/EsX42tm6wInADODvJcr3tK2t\ne7tVHUAffRl4Fvh5fmQfvqP98UbSPhnAlBaup5RGmwGvBj4pqTj/Z4A7gaf7FVWOkuWbtbxe1tUh\naWSLlv3GViy3Vcp8CSLi1YjIJyD1UHZxRHQ2JbiS6sTXzrrddk0qb0BELImIwXIA06MWJypos32s\nkWQVwCTgzeSOUrIf+U8AE6nzJiUdJek2SXMkLZR0p6R965TrkvQjSZ+RdB+plrZ7d8FI+qmklyV9\nPDdulKQfSno8m/awpGMkKVfmtXNgkr4m6ZFsXZv3sK58bA9KWpS9j/cVytXON20uaaKkecCfc9M/\nIOnPkl6UNF/SFZI262a1a0q6XNKCbNv9sJi8JR0o6UZJs7P3e7+kQ3p4H7tJuiuL/35Jexeml2oP\nz58TknQicGY2aWY2rbPWJl3vnFWZzykr9x/Zdn4+2w5/l3R4T7EV48te1z6Xtyqdf5ov6TlJl0ha\nobfldbOOsvv1btlnPl/SC9n+c2o2bRfgDtJ369Lctqt77qfEtu6Q9G1Jj2TbdYakUyW9obCcmZKu\n6m1/yMr+u6RfS5or6SVJUyXtWXIbldrfJb0/236Lsn3hP2ufWa7MLZLu7mY9D0m6ppdYbpF0U+51\nbV//pKRvKp3GWCTpBklvzZU7N/vcltlPJE2S9FR+v5W0h6Qp2Xt+XtIfJG1RmG+0pAnZOl/OlnFF\n7nOcAbwdeL9eP/1yUyHunQvLPEzSo9m++BdJ763znkdKOjnb1s9lMU6R9P5cmQ2AZ0j75Pjc+vPf\np00l/SbbJxZJ+j9JH6mzfbaQdFMW0xOSvkkjuSciSg/AAUAnsB1wK3BpbtrHgMXA2qSmjKsK8z4O\nnAscCnwNmJota49CuS7gflLt7FvAIcDWwAbZtCOzciNIVeOFwIdy868I3JNt6JOBL5GaSjqBH+TK\n1ZZ3H/AwcDRwOPBvPbz/LlITzWzgm8BRwGPAi8AWuXIn5pb9O+Bg4JBs2v8DXgWmA1/P3uMzwBxg\n/TrLuAe4IttuP8/GXVqI66/AxVn8XwauycodWig3A3gQmAucmn0Od5Oayj6YK7dLtr12zo2bADxW\nZ3uckP2/FfA/2XxfJdWyPwOsmFv3JQ18Trtl65mc7QuHAOcAvyyxv74WX2Gb/g34dfa5/CRb52kl\nlldvG/S6XwNbkA6E/gJ8JXuvZwA3Z9PXyvaDLuDC3LbbsJs4tuxlW1+aLeuX2faakL3+bYP7w1qk\n7+NzwElZuWlZuY/V+U59Ljeu7P6+Len0wqOk79U3gCeAu4DOXLmDsve9ReG9bJ+t+zO9fIY3AzcV\n9vXaPnEH6Tv0bdJ3emqu3Huz9e5bWN6KwAvAOblxn83K/pH0fTwqe19zC+/5NmAeMB44EDgWuAF4\nbzb9o9n+dT8wLvuMP9jDd/TQ7L3cDBwGfC/bzg8X3vObgX9m0/8z+1weIO2jW2dl3phN6wJ+w+v7\n2JbZ9LcD84F7s/d3aLbezsI+MTr3eX8LOJK0z92dlV2/p89rqW1dtmC24nyy+jJp510+m/Yr4Ibc\nl6CYrJYvvO4g/fBfX+cHZjGwaWH8a8kqm/eX2Q71wUK5b5Ha8DcqjP8u6UvzlsLy5gOrl3z/Xdn7\nf0du3HqkhPmbOj+Kl9VZxl3Av4BRuXFbkb74E+os43eF+c/LYtiyu22bjbsGeLjOj1NxZ1oFeBK4\ns/AF7lOyyl5/vbsdkGWTVdnP6Wxgfl/20x7iq23TnxbK/RZ4psTy6m2DXvdr0o97J7BaD8seQ+GH\nvpdY6m5r0oFdF/Djwvgzs/K7NLA/nJ2V2zE3biXSD/Cjdb6j+WRVdn+/ivSjPzo3bqNsX8gnqzeR\nvm/fLby/c7L9acVetlt3yeo+oCM3/qsUkiIpeV5eWN4ns3LvyW2XecCFhXJrkn5rfpy9HkXu4LuH\neO/Nx1uI+7XvKDCSdG5rKjAiV+6z2Xry71nAcoXlvSn7nH6WG/dmCt+h3LQbss+2uJxbgQfr7Dtj\nCsudX2//7WnoT9f1y0nZdy9JKwN7kY726oqIV2r/S1oVWI3UNLZdneK3RMRD3SzqDaRMvyfp6PXG\nwvRPZMtdIOnNtQG4kdShpNi09ZuImNdd3HXcHhGvNUNExBPAlcDu+WYAUvX5J/kZJa0NbEP6ki7I\nLeNe4PrsPeUFcH5h3Lmkne21soVt+6bs/U4BNpK0SmH+pyLiyty8LwC/ALaVtFZPb7zJyn5OzwEr\nSeq2KbiPlvlcsjjenO3HfVtYuf36uezv3oV9pBX2JL3Hswvjv0/abz5cGF9mf9gDuCMipubKvQT8\nFNiw2LxVU3Z/Vzr3/UHgioiYnSv3GOmgi9y450nft3G59YwAPgX8PiIW1YulhEti6XOqfyZtr41y\n434N7Kmlzz9/GngyIm7LXu9GSkS/LOzXQWoB2TUrt4iUiN+f7Tf99U5SEvhZLH2ediIpMbwmkiXw\nWp+A1Ui/q3dS//d4KVn5XUnbY1ThfV4HvE3SOlnxPYC/RMTfcuufSw+5ojsNJ6uImEPKrp8B9smW\n9ZvuykvaK2vnXkQ68niGVHUcVaf4zB5WfTypevyJiPhznelvAz5EOsrID9eTdpjiD3JP66rnkTrj\n/kFK3GsWxs8ovN4gV75oOrCGpBV7Wd+jpKOdDWsjJL0na2N/kfTD+CypWQeW3b7dxU9+mQOg7Od0\nQRbf1Vl798VNSFyPF17Xvsyr9XVBJffrX5GafH4GzM7OcXyyRYmrVrtZ6nPOksBzvL4P1pTZHzYA\n6h08Ts9N7y6W/PKK89b297VIzWn1Yqk37hfA+pLem73eLVvGZd3EUcYThdf19olfkb7nHwWQtBLp\nx/jyXJm3kZLczSy9Xz+Ti5NInSOOzeafLelPko6WNLrB+DcgfW8ezY/MEvDMYmGly5DuITX9zc3i\n+zD1f4+LNia9x++w7Pd3fFam9v3dgNQMWdRdZaRbjXZdr5lI+gKuA1yTHZUtQ6kDwpXALaQv8r9I\nTX1fIHeElNPT0dG1pB+5YyTdEsv2iBlB+sE7g/q9WYpfnEaPxMpoxbIj/0LSRqSDhunAEaQv3auk\nHe+/aN8Lv0t9ThHxrKR3kDrZ7JENB0r6eUQc2OC6u+uV2KfkUXa/joiXgZ0l7Ur6XD5EOiK/UdLY\nyNpGmqwVy2wXk0k/rvuTmp32J51TK7ay9EWv+0RE/FXSTFIt7pekpLUCSyerEaRtvz/p3HbRaz0R\nI+IcpetRP07av08GjpO0a0Tc0/hb6Zmk/UlN2r8jNQ8/Q3r/x7N0TbI7td+Us0ifRT31DjL6pb/J\n6vekJpV3k7583dmH9MO9e+S6jUo6qIF1/gX4Menk5a8l7V2o9j4KrBwRNzew7DLeVmfcpqR29Gd7\nmXdWrnzRZsCcOs0Yb8vNB+moZgSv19o+SqrCfyQinqwVkvTBbmLYuM64Wjwzu428nL78QJb+nLJ9\n5o/ZgKQLgf+U9J2sqagqfdqvs/d6M3CUpOOAU0jNKTfR9+TSXflZpP3jbeSOXrMmvVVZel+CcvvD\nLOrvs5vnpncXS355ea/t75JeIR3h14tlme9bRHRJmggcIOkbpM5dP2lR0i+6HDg8azL+NDAzIu7I\nTX+UlOCejYib6i0gLyJmkJpsz1bqfXgP6XxkrSdo2fc0K1vvxsCfaiMldZBqyPnkty/pXOMn8guQ\ndHIxvG7WVfvOLS7xHmdR/zezu97P3erXUXfWbn0Iqer3vz0U7SS98deSo6QNSTtZI+u9CfgP0lF2\nsep/ObCjpLHF+ZS6Snc0ss6cHSVtm1vmeqSEMbm3L0tEPE3qBXOAcncekLQlMJbsxzgfMqlXT97h\npG15bfa69iP52mcpaRTw+W7CWFe5rslZHJ8F7oqIZ3qKv4SXsr9l2uBLfU6qf7upe7O/A3L9XQ9K\n7ddZG3/RPaTPt/Ye+rLteip/dbbc/yqM/3oWa3EfK7M/XA28S9K7c+VWIvUWmxERD9QLsOz+nh1s\n3gB8PDvPVSu3MakWWs9lwOqkg+WVaOAcSIN+RfrMPk+qDf2qMH0yqaPH8ZKWqQxIWiP7u6KWvX50\nBqmTSX78S5TbJ+4kNed9SUtf/7o/yzZvL1OLzD7bHQujF2Z/l1p/RDxLak04OP955Za1Ru7l1cAO\nkt6Zm74m6fRRnzRSs1qqqSQiyrQT/5HUi29ydkQ0mtSb8GFS76U+i4grJR0I/ELSCxFRu67oe6Tk\n8QdJl5K6pK6UrWcf0lFGXzpUFN0HXCvpXFJz26Fk1yKUnP9o0gf4F0kXk9rAv0JqIz+pTvl/l3Ql\nKTntBOwH/Hd2khrSCc3FpPf7E1Jvri+SmiCW2ZFIzWsXSdo+K3MQqX35gEK5Rs6n/C2b77uSfpnF\ndVU3J73Lfk4XZQnrJlJ32w1J2+uuiJi+zFIHVtn9+gSl62H+SDrSHE3abx4nNWNBOiJ/DjgkO/f4\nEvDXiJjZzbq729Z/l/RzUs1zNdJR9rtJR+q/i4g/FZZTZn84ndSsea2kH5E+l8+Tzkfs08s2Kru/\njyclsNuzmvNypAO1+0idNJYSEXcrXYf5SeCBfKenVoqIuyQ9Sjon/AaWbgIkIl6QdCjpvNq07LN5\nFlif1AR8K+mAcxNSM/DlpG7jS0jbci3Sdaw1fyPtE98kNa09k2uNyDdRLpY0HvgRcHO23A1JXeIf\nYela0h+AfSRdQdonNyJdxnE/sHJumS9LegD4tKSHSZ/7fRFxP+mz+TNwr6SfkWpbo0kJ7y2kSxEg\nNTN+lvQdOYeUAL9EqrX37be/bLfBrNJwAFnX9V7KPQZcWRj3eVL/+oXZRvkcqStxZ6FcJ7lrFnLj\nN8imHVEYf0g2/ozcuDeSmlgeIjXTzM427H+RdU/tbnm9vK8u0s4wLlv2QuD/gPcVyp2YLbtul3hS\n088UUtf7+aTm1GJX/RNJO/CmpC/Ec6RrFX4IvKFQ9sOkbqQvkX70vp5t76W6htY+F9K1L3fnPou9\nC8vrruv6o4VyncC3C+OOJ/0IL86vP1v3xYWyZT6nvUk9wv6VlZlB6iG5VonPa6n4uvtceH2/7rEb\nbTfb4PP0sl8D7yedH3giew9PkGoGby0say9SrfGVLJ4eu7H3sK1HkC4NqF3oPpN0MnxkYf4ZpC7j\nPe4PWdkNSbWIudl+NpXc9Y2F79TnCuN73d9z2+nObBs9TDro+h7wUjfv/yjSd/KYPnyHbwZurLOv\n71PmvWTTvpNNe7CH9exMStLzsu31D9K1kNtm01cn/ZbcT6qJzSPdvq4Yx1rZZ/Rcts6buvuOZuMP\nI33XFpJ6H+5E+o36Y6Hcsblyd5Jaqert3+8mXX+2KFtf/lKQDbN5nsz2s8dJvy8fLyzj7aSDzZey\nMseRkmifuq4rW5iVoHQl/XkR0evdE8zandIdEu6NiI9WHUt3JP2edK3TMue9JH2N1CV/w4j454AH\nNwhkPU6fJV0QfnDV8fRHu/YUM7NhRoVbGUl6G+larO464XyBdE2mExVQ5xwYpFaD1el+Gw4a/e0N\naGbWLI9l5y8fIzUxHUJqXvperYDSBbkfIzUtbkl2zZMBqSPD2aSLdeeS7oryBdIdVbq9BnawcLLq\nm2BoX79iw0u77c/XkHr5rk06b3c7cHxE5C90XZPU828+cGpEFHs3DmczSeeEvkqqTc0j3SfyuBgC\nd5r3OSszM2t7PmdlZmZtz82AQ5zSzSV3JzUR+AnONlStQDrPNTnSjVKbSukZU2v0WrB7cyKieE9K\n6wMnq6Fvdwbu6n6zqu1Humdp00hafyTMWty/xSyUtLkTVuOcrIa+menPPvTvwLBRtfsODzxdV/dB\nuwMiTjgWnXxGdes/o3jz/gFyzxGwTfHpJAPg+enwf/tD/+9vWc8ai0nPtGnkGTrPAL9JF8CvwbJ3\n/LeSnKyGvqzpbw3SzfEH2goVrRe09TsqWS9AvGlUtetfbaVqVjxyFKzW6yORWqllTd3rku4j1Ff+\nkW0Od7AwM7O256RvZlZCB439YPb3MQ+WOFmZmZWwHDCywfms/7wdrcW2rDqASmjvT/ReaChar96D\nv4cG16yq5WRlLbZV1QFUQnt/quoQqrH+0E1WrllVy9vRzKwE16yq5d6AZmbW9lyzMjMrwc2A1fJ2\nNDMrYTka+8H0j2xzeDuamZXgmlW1vB3NzEpwsqqWO1iYmVnbc9I3MyvBXder5WRlZlaCmwGr5e1o\nZlaCa1bVcrIyMyvBNatqeTuamZXgmlW13BvQzMzanmtWZmYluBmwWt6OZmYl+HZL1fJ2NDMrwTWr\nank7mpmV4A4W1XIHi0FI0mGSZkhaJOkvkravOiazoa5Ws+rr4BpBczhZDTKSPg18HzgR2Ba4B5gs\naY1KAzOzASHpw9lB6kJJ8yT9rjB9PUl/lPSSpKclnSlp0P/WO+kPPkcAP4mIXwBIOgT4MPAF4Mwq\nAzMbytrhnJWkfYGfAt8AbspC2jI3fQRwNfAUsAOwLnAZ8CrwrSaGMuCcrAYRSSOBMcB3a+MiIiTd\nAOxYWWBmw0DV56wkdQA/BL4eEZfmJj2Y+393YDNg14iYA9wr6dvA6ZLGR8SSJoUz4AZ91XCYWYO0\n788ujJ8NrD3w4ZgNH8t1wMjl+j4s17weFtuRakpImibpKUlXS3p7rswOwL1ZoqqZDIwC8uUGHScr\nM7MSOjpgueX6PnQ0L1ltBIh0vvpkUvP/fOAWSatmZdam/sFsbdqg5WbAwWUO0AmMLowfDTzd86zX\nAisUxm0JbNWk0MwGyOOT4IlJS49bvKDlq11uBIzsJfFcviQNeQui53kknQYc20ORADbn9crFKRFx\nRTbvgcA/gU8CP+t5TYObk9UgEhGLJf0N+CBwFYAkZa9/1PPcHwLWaXGEZgNg/XFpyJs/DW4cU008\nOZ9aLg15d3XCe17ucbazgAm9LPoxsiZAYHptZES8KukxYP1s1NNA8VKW0blpg5aT1eDzA+DSLGnd\nQeod+Ebg0iqDMhvqlmvw/NNy6nl6RMwF5va2nOw7/wqwKXB7Nm4ksCEwKys2FThe0hq581ZjgQXA\nA32Pvn04WQ0yEXF5dk3VyaQjpruB3SPi2WojMxvaah0s+jxfk9YfES9I+jFwkqR/khLUMaRmwl9n\nxa4jJaXLJB1Lak75DnBeRCxuUiiVcLIahCLiAuCCquMwG1ZG0Fg/9K6mRnEUsBj4BbAi8FfgAxGx\nACAiuiTtBVxIqn29RGp1ObGpUVTAycrMrIxGL7RqYrKKiE5SbeqYHso8AezVvLW2B3ddNzOztuea\nlZlZGY0+0Kq5zYDDlpOVmVkZjTYDdjY7kOHJycrMrIxGO1j4ZEtTOFmZmZVR9Z1shzknKzOzMho9\nZ+Vf2aZwBdXMzNqec76ZWRk+Z1UpJyszszJ8zqpSTlZmZmX4nFWlvBnNzMpwM2ClnKzMzMpwM2Cl\nnPPNzKztuWZlZlaGa1aVcrIyMyvDHSwq5c1oZlaGO1hUysnKzKwMNwNWysnKzKwMJ6tKuYJqZmZt\nzzUrM7MyOmisluSaVVO4ZmVmVkatGbCvQxOTlaTjJd0m6SVJ8+pM31rSREmPS1oo6X5Jh3dTboqk\nRZJmSTq6eVG2hmtWZmZltMc5q5HA5cBU4At1po8BZgP7AU8AOwE/k7QkIi4AkLQKMBm4DjgY2AqY\nIGl+RFzU1GibyMnKzKyMNmgGjIiTACQd0M30CYVRMyXtBOwDXJCN25+U9A6KiCXAdEnbAkcCTlZW\nrb/97WC22267qsOwATL+N6o6hAH1L+CnVQfRvkYB+SbDHYApWaKqmQwcI2lURCwY0OhKcrIyMyuj\nPZoB+ySrVX0K2DM3em3gsULR2blpbZms3MHCzKyMFnWwkHSapK4ehk5Jm/Q1XElbAlcA4yPixr7O\n325cszIzK6NEzWrSP9KQt+DVXpd8FlA811RUrAn1SNIWwA3AjyPitMLkp4HRhXGjc9PakpOVmVkZ\nJTpYjNs8DXnTnoExk7qfJyLmAnP7G16NpLcDNwITIuKEOkWmAqdI6oiIzmzcWOChdj1fBW4GNDMr\npz2us1pP0jbABkCHpG2yYaVs+pbAzaQOEz+UNDob1sgtZiLwKnCJpC0kfRo4HPh+8yJtPteszMwG\nj5OBz+VeT8v+7gpMAfYF3kzqnr5/rtwsYCOAiHhe0ljgfOBOYA7pvNbFrQ29f5yszMzKaIPegBFx\nIHBgD9NPAk4qsZz7gF2aF1nrOVmZmZXRBhcFD2dOVmZmZbRBzWo4c7IyMyvDyapSTlZmZmU4WVXK\nXdfNzKztuWZlZlaGO1hUysnKzKwMNwNWysnKzKwMJ6tKOVmZmZXhZsBKOVmZmZXhmlWl3BvQzMza\nnmtWZmZluGZVKScrM7MyRtBY4nH7VVM4WZmZlVF7PlUj81m/eTOamZXhZsBKuYJqZmZtzzUrM7My\nXLOqlJOVmVkZ7mBRKScrM7My3MGiUs75g4yk90m6StKTkrokfbTqmMyGhVozYF8HNwM2hZPV4LMS\ncDfwZSAqjsVs+Kg1A/Z18K9sU3gzDjIRcW1EnBARVwKqOh4zGziSjpd0m6SXJM3rpsz2km6QNF/S\nPEnXStq6UGZrSVMkLZI0S9LRA/MOGudkZWZWRns0A44ELgcurDdR0krANcBM4F3Ae4AXgGsldWRl\nVgEmAzOA7YCjgfGSvtjUSJvMp/7MzMpogw4WEXESgKQDuimyGbAacGJEPJmVPQm4B9gAeAzYn5T0\nDoqIJcB0SdsCRwIXNS/a5nLNysysjMFxzuohYC5wkKSRklYEvgg8QKptAewATMkSVc1kYFNJowYy\n2L5wzWqYOOKIIxg1aun9cNy4cYwbN66iiMwacy9wX2HcywOx4kFwUXBEvChpV+AK4IRs9D+A3SOi\nK3u9NqmGlTc7N21BywNtgJPVMHH22Wez3XbbVR2GWb9tlQ15/wJ+2uoVt6gZUNJpwLE9FAlg84j4\nR2+rkrQCcDFwK/DpbO1HAVdLemdEvFIy6rbjZDXIZCdQN+b1noAbSdoGmBcRT1QXmZlNuhYmTV56\n3IIXe53tLGBCL2WKNaHu7AdsEBE71EZI2g+YD3yM1DnjaWB0Yb7a66dLrmfAOVkNPu8EbiYdbQXw\n/Wz8z4EvVBWU2ZBXohlw3F5pyJs2Hcb00NoeEXNJ55maYUWgqzCu9ltRO3s2FThFUkdEdGbjxgIP\nRURbNgGCO1gMOhHxp4gYEREdhcGJyqyV2qCDhaT1spaUDYAOSdtkw0pZkeuB1SSdL2kzSW8n1doW\nkw5yASYCrwKXSNpC0qeBw3n9wLctuWZlZlZGe3SwOBn4XO71tOzvrqQefg9J+ghwInA7qZZ1F6mD\nxWyAiHhe0ljgfOBOYA4wPiIubmqkTeZkZWZWRntcZ3UgcGAvZW4EbuylzH3ALs2LrPXcDGhmZm3P\nNSszszL8PKtKOVmZmZXRHueshi0nKzOzMpysKuVkZWZWRht0sBjOvBnNzEqIERAN1JLC56yawpvR\nzMzanmtWZmYldHZAZwO/mJ0+Z9UUTlZmZiV0NZisupysmsLJysyshM4OsaRDvRdcZr7afWStP5ys\nzMxK6OzooHO5vp/m7+zoApb0Ws565mRlZlZCV0cHnR19T1ZdHcLJqv/cG9DMzNqea1ZmZiV0MoLO\nBm5H0dl7ESvBycrMrIROOljiZFUZJyszsxK66KCzgZ/M4jPmrTFOVmZmJTTeDOh01QxOVmZmJaSa\nVd+TVZeTVVO4N6CZmbU916zMzEroarAZsMtdLJrCycrMrIQljGioN+ASN2A1hbeimVkJXSxHZwND\nV5PqBJI2kHSRpMckLZT0sKTxkkZ2U351Sf+U1CnpTYVpW0uaImmRpFmSjm5KkC3kmpWZWQmNNwM2\nrU6wGSDgS8CjwJbARcAbgWPqlL8YuBtYJz9S0irAZOA64GBgK2CCpPkRcVGzgm02Jysbsl56Zfg+\nm2F8DK+7fE+bNo2fjhlTdRgtFRGTSUmmZqaks4BDKCQrSYcCo4DvAHsUFrU/MBI4KCKWANMlbQsc\nSUp+bcnNgGZmJdSus+r70NKf2VWBefkRkrYAvgV8lvrXJO8ATMkSVc1kYFNJo1oVaH+5ZmVmVkLj\nt1tqTQ1f0sbAV0g1otq4NwATgaMi4smsTNHawGOFcbNz0xa0INx+c83KzKyE2u2W+t7BoudkJek0\nSV09DJ2SNinM8xbgGuBXEXFJbtLpwAMRMalWtPB30HLNysyshM4Sd7C4ftI8bpi0VKscLy7o9Tqr\ns4AJvZR5rSYkaV3gJuDWiDi4UG5XYEtJn6wVz4ZnJZ0aEScBTwOjC/PVXj/dW7BVcbIyMyuhTG/A\nD4xbkw+MW3Opcf+Y9hJfGnN/t/NExFxgbpkYshrVTcD/AV+oU2QfYMXc63eRegW+l9cT3lTgFEkd\nEVHLpGOzfQM2AAAdNUlEQVSBhyKiLZsAwcnKzGxQyGpUtwAzSL3/1pJS615EzM7+zijMsyapZvVg\nRDyfjZ4InABcIukMUtf1w4Gvtf5dNM7JysyshMbvut60rgG7ARtlwxPZOAEBPQa21HUMEfG8pLHA\n+cCdwBxgfERc3KxAW8HJysyshKp7A0bEz4Gf93GeP1EnkUXEfcAuTQlsgDhZmZmV0PjDF4fvxenN\n5GRlZlZCGzQDDmtOVmZmJTT+8EXXrJrBKd/MzNqek1WLSTpA0odzr8+U9Jyk2yVtUGVsZlZeZ/Y8\nq74ObgZsDm/F1jseWAQgaUfgMNI1EnOAsyuMy8z6oLPB2y216t6Aw43PWbXeesAj2f8fB34bET+V\ndBvpAj8zGwR8zqparlm13ovAm7P/xwLXZ/+/zNK3RTGzNtamjwgZNlyzar3rgYsk3QVsAlydjX87\nMLOqoMysb6q+KHi4c8pvvcNIN45cE9g3u2klwBhgUrdzmZnZa1yzarGIeI70gLTi+BMrCMfMGuQ7\nWFTLNasBIOl9kv47667+lmzcZyW9t+rYzKwcn7Oqlrdii0naF5hM6r6+HbB8NmkUqVu7mQ0CXQ0l\nqg7XrJrEyar1vgUcEhFfAhbnxt9GSl5mNgh0NViz6vLPbFN4K7bepsCUOuMXAKsOcCxmZoOSO1i0\n3tPAxizbTT3/mGkza3O12yc1Mp/1n5NV6/0MOEfSF0hP7Fw3u+3SWcB3Ko3MzEpzb8BqOVm13umk\n5tYbgTeSmgRfAc6KiHP7siBJxwF7A5uROmzcDhwbEf9oasRmtgw/z6paTlYtFhEBnCrpe6TmwJWB\nByLixQYW9z7gXOBO0md3GnCdpM0jYlGzYjazZfnegNVyshogEfEq8EA/l7Fn/rWkzwPPkO6GcWt/\nlm1mPas9IqSR+az/nKxaTNJKwDeADwJrUeiBGREb9WPxq5LOg83rxzLMzNqek1XrXQTsAlwG/IuU\nXPpNkoAfArdGRL9qbGbWu84GO1g080a2kq4E3kE68J0P3EA6b/2vXJn1gB8D7wdeAH4BfCMiunJl\ntgbOA7Yntc6cFxHfa1qgLeBk1Xp7AB+OiNuavNwLgC2A9zR5uWZWR5ucs7oJOJV04PsW4PvAr0mX\nwiBpBOnJDk8BOwDrkg6UXyXdoABJq5DuqnMdcDCwFTBB0vyIuKiZwTaTk1XrzafJzXSSzgP2BN6X\nP6LqyRFHHMGoUaOWGjdu3DjGjRvXzNDMWm7SpElMmrT0AwsWLFjQ8vV2NdgbsJl3sIiIc3Ivn5B0\nOvB7SR0R0QnsTuotvGtEzAHulfRt4HRJ4yNiCbA/MBI4KHs9XdK2wJGklqC25GTVet8GTpZ0QEQs\n7O/CskT1MWCXiHi87Hxnn302223nuzvZ4FfvIGvatGmMGTOmpettt67rklYH9gNuyxIVpNrUvVmi\nqpkMXEh6ht49WZkpWaLKlzlG0qiIaH3mb4CTVet9HXgrMFvSTJa+PyARUTqDSLoAGAd8FHhJ0uhs\n0oKIeLk54ZpZPe3y8MWsNvUV0nWbU4G9cpPXBmYXZpmdm3ZP9rd495x8GSerYeqKJi7rEFIHjVsK\n4w8knUQ1swpNn/R3Hpz096XGvbKg5+NISacBx/ZQJIDNcxf/n0lqrtsAOJF0TmqvbuYdMpysWiwi\nTmrisnzBhllFytxuaZNx27HJuKUbS56Z9iQTx5zX02xnARN6Wf1rNaGImEc6D/6IpAdJ567eHRF/\nJd2LdPvCvLUWmKdzf0f3UqbtOFkNAEmrAp8gNQd+LyLmSdoOmB0RT1YbnZmV0apzVhExF5jbYFi1\ngGrPyZsKHC9pjdx5q7Gkpr0HcmVOyXXKqJV5qF3PV4GTVctl1zPcQNpZNiTd2HYesA+wPvC5yoIz\ns9Kq7rou6V2kWtOtpF7GGwMnAw+TEhCk7ugPAJdJOhZYh3TD7PMiona+fCJwAnCJpDNIXdcPB77W\nlEBbxM1KrfcD4NKIeBuQb7y+Gti5mpDMrK9qt1vq69DE3oALSQe5NwAPkg587wbeX0tE2YW/ewGd\npBtd/wK4lHRui6zM86Sa1Iak+4x+DxgfERc3K9BWcM2q9bYnXXhX9CSp542ZDQKdLNfgHSya8zMb\nEfeRbtvWW7kn6KXDRbasXZoS2ABxzar1XgHeVGf8JsCzAxyLmdmg5GTVelcBJ0gamb0OSesDZwC/\nrS4sM+uL2h0s+jo08w4Ww5m3Yut9nfQMq2eAFYE/AY8ALwLfrDAuM+uDzgaTlR8R0hw+Z9ViWVfQ\n3SS9B9iGlLimRcQN1UZmZn1RdW/A4c4pv0UkrSgpf5JzL1JX07WBPSWdKWmFaqIzs75qg96Aw5pr\nVq1zAPBh4A/Z668A9wO1x89vRrrN/9kDH5qZ2eDiZNU6+5Hu4ZX3mYh4DEDS/sBhOFmZDQrt8PDF\n4czJqnU2Bu7NvX4Z6Mq9vgM4f0AjMrOG+ZxVtZysWmdVXr9fFxGxZmH6iPx0M2tv7fY8q+HGyap1\n/glsCTzUzfStszJmNgi0y/Oshisnq9a5mvSE4D8WH4woaUXSvbr+WElkZtZnZR4R0t181n9OVq3z\nXeBTwEPZo+hrD07blNQzcLmsjJmZ9cLJqkUiYraknYALgdMB1SYB1wNfjoji46fNrE35nFW1nKxa\nKCJmAB+StDqpdyDAI9mTPs1sEHFvwGo5WQ2ALDndUXUcZta42h0sGpnP+s/JysysBF8UXC0nKzOz\nEtwMWC3XT83MrO25ZmVmVkJXg70B/fDF5nCyMjMrobPBZkCfs2oOp3wzsxLa4XlWkq6UNEvSIklP\nSfqFpHXqlPu8pHuyck9LOrcwfWtJU7LpsyQd3bQgW8Q1q2FizJg/MdxuRfibpW5yP7x8gquqDmGA\nPdryNbRJb8CbgFNJz8J7C/B94NfAe2sFJB0JHAEcRbpkZiVgw9z0VYDJwHXAwcBWwARJ8yPiomYG\n20xOVmZmJbRDb8CIOCf38glJpwO/l9QREZ2SVgW+A3w4Im7Jlb0v9//+wEjgoIhYAkyXtC1wJNC2\nycrNgGZmg1B2Z5z9gNsiojMbvRvp1m7rSXpA0hOSfiXp33Kz7gBMyRJVzWRgU0mjBiT4BjhZmZmV\nUOsN2Neh2b0BJZ0u6UVgDrAe8PHc5I2ADuA44HBgX2B14HpJtZa0tYHifUln56a1JScrM7MSljTY\nwWJJLz+zkk6T1NXD0Clpk9wsZwLvINWiOoHLctNGkE7vfDUiboiIO4BxwNuAXZu5PQaaz1mZmZXQ\nxXK9drBYOOl/WTjpD0vPt+CF3hZ9FjChlzKP1f7J7jU6D3hE0oOkc1fvjoi/kjpeAEzPlZ8jaQ6w\nfjbqaWB0Yfmjc9PakpOVmVkJZS4KXn7cx1l+3MeXGrd42n3MGfPRbueJiLnA3AbDqgW0fPb3tuzv\npsBT8Nq5rTWAmdm0qcAptU4Z2bixwEMRsaDBOFrOzYBmZoOApHdJOkzSNpLWl/QBYCLwMCkBEREP\nA1cB50jaUdKWwM+BB4BbskVNBF4FLpG0haRPk85vfX9g31HfuGZlZlZCJyMYUe0jQhYC+wDjSddO\n/Qu4Bjg1Ihbnyn0WOBv4A9BFSlJ71GpREfG8pLHA+cCdpI4a4yPi4mYF2gpOVmZmJXR1ddDZ1cB1\nVg3MU09E3Ad8sES5F4EvZUNPy9qlKYENECcrM7MSOjtHwJIGaladPtvSDE5WZmYldC7pgCUN3G6p\ngQRny3KyMjMroauzo6GaVVenk1UzuH5qZmZtzzUrM7MSOjtHEA3VrFwnaAYnKzOzEjqXdNC1uO/J\nqpEEZ8tysjIzKyG6OojOBn4ym9R1fbhzsjIzK2NJY13XWeJmwGZwsjIzK6PB3oC4N2BTOOWbmVnb\nc83KzKyMTsESNTaf9ZuTlZlZGZ3Akl5L1Z/P+s3JysysDCerSjlZmZmVsYTGklUj89gynKzMzMpY\nAizutVT9+azf3BvQzMzanpPVICLpEEn3SFqQDbdL+lDVcZkNC12k8099HbqqCHbocTPg4PIEcCzw\nMCDg88CVkt4REdOrDMxsyHMHi0o5WQ0iEfHHwqhvSToU2AFwsjJrJXewqJST1SAlaQTwKeCNwNSK\nwzEb+lyzqpST1SAjaUtScloBeAHYOyIerDYqM7PWcgeLwedBYBvgXcCFwC8kbVZtSGbDQK1m1deh\nBTUrSW+QdLekLklb58ZvLWmipMclLZR0v6TD68y/taQpkhZJmiXp6OZH2VyuWQ0yEbEEeCx7eZek\ndwFfAw7tec6LSS2GeTtng9lgMiUb8ha2frXt1Qx4JvBPYKvC+DHAbGA/UoesnYCfSVoSERcASFoF\nmAxcBxycLWOCpPkRcVFLom0CJ6vBbwSwfO/FDgLe2upYzAZAvYOsR4EjW7vaNklWkvYAdgP2BfbM\nT4uICYXiMyXtBOwDXJCN2x8YCRyUHfxOl7QtaQM6WVn/SfoucA3wOLAK6ehpF2BslXGZDQttcAcL\nSaOBnwIfBRaVnG0UMC/3egdgSpaoaiYDx0gaFRELmhJskzlZDS5rAT8H1gEWAH8HxkbETZVGZTYc\n1C7ybWS+5pkAXBARd0naoLfCWa3qUyxdA1ub108l1MzOTXOysv6JiC9WHYOZNZek00gX+3cngM2B\nDwErA2fUZu1luVsCVwDjI+LGJoRaKScrM7MyypyzmjoJ/jJp6XGLeq2onEWqMfVkBrArsCPwirRU\nnrpT0v9ExIG1EZK2AG4AfhwRpxWW9TQwujBudG5aW3KyMjMro0yy2n5cGvJmTYOTx3Q7S0TMBeb2\ntnpJXwW+mRu1Lulc06eAO3Ll3g7cCEyIiBPqLGoqcIqkjoioNVKOBR5q1/NV4GRlZlZOxb0BI+Kf\n+deSXiI1BT4WEU9l47YEbiJ1xPph1iEDoDMi5mT/TwROAC6RdAap6/rhpEtg2paTlZlZGe15b8Ao\nvN4XeDOpe/r+ufGzgI0AIuJ5SWOB84E7gTmk81oXtzTSfnKyMjMro02us6qJiFlAR2HcScBJJea9\nj3TZy6Dh2y2ZmVnbc83KzKyMNqtZDTdOVmZmZbTBHSyGMycrM7My2uMOFsOWk5WZWRluBqyUk5WZ\nWRlOVpVyb0AzM2t7rlmZmZXhmlWlnKzMzMpwb8BKOVmZmZXh3oCVcrIyMyvDzYCVcgcLMzNre65Z\nmZmV4ZpVpZyszMzKcAeLSjlZmZmV4Q4WlXKyMjMrw82AlXKyMjMrw8mqUu4NaGZmbc81KzOzMtzB\nolJOVmZmZXTRWJNeV7MDGZ6crIaJL3Mk61YdxAB7pOoAKnQ6H6s6hAH1JHBuq1eyhMZqSa5ZNYXP\nWZmZlVHrYNHXoQUdLCS9QdLdkrokbV2Ytr2kGyTNlzRP0rV1ymwtaYqkRZJmSTq6+VE2l5OVmVkZ\ntXNWfR1aU7M6E/gnEPmRklYCrgFmAu8C3gO8AFwrqSMrswowGZgBbAccDYyX9MWWRNokTlZmZoOI\npD2A3YCjABUmbwasBpwYEQ9HxHTgJGA0sEFWZn9gJHBQREyPiMuBHwFHDkT8jXKyMjMro9bBoq9D\nEztYSBoN/JSUcBbVKfIQMBc4SNJISSsCXwQeINW2AHYApkREvs43GdhU0qjmRdtcTlZmZmW0xzmr\nCcAFEXFXvYkR8SKwK/BZUjJ7ARgL7BkRtbS5NjC7MOvs3LS25GRlZlZGI4mqRA9CSadlHSW6Gzol\nbSLpcGBl4IzarHWWtQJwMXAr6ZzVTsB9wNWSlu/X+6+Yu66bmZVR5qLguZNg3qSlx3Uu6G3JZ5Fq\nTD2ZQaox7Qi8Ii2Vp+6U9D8RcSCwH7BBROxQmyhpP2A+8DHgcuBp0jmsvNrrp3sLtipOVmZmzfLm\ncWnIe2kaTB/T7SwRMZd0nqlHkr4KfDM3al3SuaZPAXdk41Zk2bNkkQ21lrSpwCmSOiKi1kg5Fngo\nInrNrFVxM6CZWRkVd7CIiH9GxAO1AXiY1BT4WEQ8lRW7HlhN0vmSNpP0dlKtbTFwc1ZmIvAqcImk\nLSR9Gjgc+H5zIm0N16zMzMpoz7uuL3WdVUQ8JOkjwInA7aRUeRewe0TMzso8L2kscD5wJzAHGB8R\nF7c00n5ysjIzK6PNbrcUEbOAjjrjbwRu7GXe+4BdWhNZazhZmZmV4buuV8rJysysDN91vVLuYGFm\nZm3PNSszszLas4PFsOFkZWZWhpNVpZyszMzKaLSjhDtYNIWTlZlZGZ3UuRtfyfms35yszMzKaDTp\nOFk1hXsDmplZ23PNysysjE4KNzcqyddZNYWTlZlZGUto7JxVIwnOluFkZWZWRqMdLJysmsLJysys\nLCeeyriDhZmZtT0nq0FM0jckdUn6QdWxmJm1kpPVICVpe+A/gXuqjsXMrNWcrAYhSSsD/w18EXiu\n4nDMzFrOyWpwOh/434i4qepAzIaP2tMX+zr45oDN4N6Ag4yk/wDeAbyz6ljMhpc2e679MONkNYhI\n+jfgh8D/i4hGHrBtZjYoOVkNLmOANYFpkmqXJ3YAO0v6CrB8RNS9EuRqYIXCuK2AbVoVqVmL3M2y\nvYpeHpA115oBG5nP+svJanC5gZRj8i4FpgOnd5eoAPYE1m1dXGYD5h3ZkPckcG7L11z90xclzQTW\nz40K4LiIOLNO2dWBvwPrAKtFxPO5aVsD5wHbA88A50XE95oWaAs4WQ0iEfES8EB+nKSXgLkRMb2a\nqMyGi7aoWQXwLeBnvH7zpxe6KXsxqSK6Tn6kpFWAycB1wMGkA+AJkuZHxEXNDLaZnKwGP98AxmxA\ntEWyAngxIp7tqYCkQ4FRwHeAPQqT9wdGAgdFxBJguqRtgSOBtk1W7ro+yEXEByLiyKrjMBv6as2A\nfR2a/vTFb0iaI2mapKMkdeQnStqCVPv6LPUfULIDMCVLVDWTgU0ljWp2sM3impWZ2eBxDjANmAfs\nBJwOrA0cBSDpDcBE4KiIeFLSxnWWsTbwWGHc7Ny0BS2Iu9+crMzMSinTDPiHbMjr7pRSIuk04Nge\nigSweUT8IyJ+mBt/n6RXgZ9IOi67nOV04IGImFRbfOHvoOVkZWZWSpnegB/KhrwHgE/0NNNZwIRe\nFlysCdXcQfod3xB4GNgV2FLSJ7PpyoZnJZ0aEScBTwOjC8upvX66lzgq42RlZlZKazpYRMRcYG4j\nEQHbks5LPZO93gdYMTf9XaRege/l9YQ3FThFUkdE1E6ojQUeioi2bAIEJyszs5Kqvd2SpB2AdwM3\nk9oWdwJ+AFxWSzIRMaMwz5qkmtWDueusJgInAJdIOoPUdf1w4GtNCbRFnKzMzEqpvOv6K8B/ACcC\nywMzgO8DZ/cy31KXt0TE85LGkm6IfScwBxgfERc3K9BWcLIyMxsEIuIuYMc+zvMn0i3ZiuPvA3Zp\nUmgDwsnKzKyU6m+3NJw5WZmZlVJ5M+Cw5mRlZlaKa1ZVcrIyMyvFNasqOVmZmZXimlWVfCNbMzNr\ne65ZmZmV4mbAKjlZmZmV4mRVJScrM7NSqr3d0nDnZGVmVoprVlVyBwszM2t7rlmZmZXirutVcrIy\nMyvFzYBVcjOgtdQ9VQdQkburDqAiQ/t912pWfR1cs2oGJytrqXurDqAiwzVJD+33XatZ9XVwzaoZ\n3AxoZlaKz1lVyTUrMzNre65ZmZmV4g4WVXKyGvpWANjpv/+bzTfffMBXft0RR7DH2WcP+HqrdvUR\nR7Cb3/eAmT59Oufuvz9k+3trPE1jiWdOswMZlhQRVcdgLSTpM8D/VB2H2QDZLyImNnOBktYHpgNv\n7MdiFgKbR8TjzYlq+HGyGuIkvRnYHZgJvFxtNGYtswKwITA5IuY2e+FZwlqjH4uY40TVP05WZmbW\n9twb0MzM2p6TlZmZtT0nKzMza3tOVmZm1vacrKwlJO0gaYmk/606loEgaYKkrtwwR9I1kraqOraB\nIGm0pHMlPSrpZUmzJF0l6QNVx2ZDg5OVtcpBwI+AnSWtXXUwA+QaYDSwNvAB0hWkQz5ZS9oAmAa8\nH/g6sCXwIeBm4LzqIrOhxF3XrekkrQT8CxgDnAT8PSJOrzaq1pI0ARgVEfvkxr0HmAKs1Yprf9qF\npKtJCWqTiHi5MO1NEfF8NZHZUOKalbXCp4HpEfEw6e4ZB1Ucz4CTtDLwWeDhIZ6oViNddH5eMVEB\nOFFZs/jegNYKXwAuy/6/FniTpJ0jYkqFMQ2Ej0h6Ift/JeApYK8K4xkIGwMCHqo6EBvaXLOyppK0\nKfAu4JcAEdEJXM7wqF3dBGwNbANsD0wGrpW0XqVRtZaqDsCGBycra7aDgA7gX5IWS1oMHALsK2mV\nakNruZciYkZEPBYRfwO+RKphfaniuFrpYSCAzaoOxIY2JytrGkkdpPM0R5JqF/nhKWBcddFVJoAV\nqw6iVSJiPqkGeZikZd6npFEDH5UNRU5W1kwfAVYFLomIB/ID8Dvgi9WG13LLZ9cbjZa0GXAu6bES\nV1UcV6sdRqpN3yFpH0kbS9pM0uHA7RXHZkOEk5U10xeA6yPihTrTfguMkbTlAMc0kD5EqkE+BfyF\n1HX/ExHx50qjarGImAFsR7qu6izgXuA6YCyplm3Wb77OyszM2p5rVmZm1vacrMzMrO05WZmZWdtz\nsjIzs7bnZGVmZm3PycrMzNqek5WZmbU9JyszM2t7TlZm3ZB0oqS7cq8nSPpdlTGZDVdOVjYkZffn\nO1fSo5JeljRL0lWSPtDHRZW+xYukAyTN7+PyzawEP3zRhhxJG5BuoDoP+DpwHzCSdO++84AtWrVq\n+pDczKw816xsKLoQ6AS2j4grIuKRiJgeEWcDO9QKSRol6SJJz0haIOkGSVs3skJJuwCXAKMkdUnq\nlHSCpG9LurdO+bslnZT9P0HS77PytVgulLRcrrwkHSfpMUkLJd0lad/c9FUl/U82/0JJD0k6oJH3\nYtaOXLOyIUXSasDuwHER8XJxekQ8n3v5G+DFrPzzwMHADZI2iYjn+rjq24D/Ak4CNiHVsl4EVgNO\nkDQmeyAjkrYFtgQ+lpv/g8AiYBdgQ+BSYA7w7Wz68cBngP8EHgF2Bi6T9Ex2V/dTSA9A3B2YS3rc\n/JB9jpYNP05WNtRsTEoUD/VUSNJ7gHcCa0XE4mz0MZL2Bj4BXNSXlUbEEkkL0r/xbG7SQknXAQcC\nf8vGHQj8KSJm5cq9AhwYEa8A0yWdAJwJfFvSG4DjgA9GxF+z8jMlvY+UYP8MrAfcFRG1DiGP9yV+\ns3bnZGVDjUqW2wZYBZgnLTXLCsBbmxzTz4CLJR1JOqc1Dvhaocw9WaKqmQqsLGm9LM43Atdr6WBH\nAtOy/y8EfitpDOlZUldExNQmvw+zyjhZ2VDzMCkhbAZc2UO5lUkPSdyFZRNcX5sAe/O/pJrT3sBi\n0vfut32Yf+Xs756kmPNeAYiIayWtn5XZjdSceX5EHNOfwM3ahZOVDSkRMV/SZOAwST+KiEX56ZJG\nRcQCUo1kbaAzIprVZPYq6fHuxZg6Jf2C9CTlV4FfFmpRANtIWj43fkfgxYh4IusO/wqwQUTc2t3K\nI2IucBnpXNatpGZEJysbEpysbCg6DLgVuEPSicDfSfv6WNI5nrdHxA2SpgJXSDoW+AfwFlLN5HcR\nMa3+ons0k9R09wHgHmBhLlleBEwn1freU2feN5CaCk8F/h0YD5wLEBEvSjoLOFtSR/beRmXLWRAR\nl2U9C/8G3E9qytwLeKCB92DWlpysbMiJiBmStgO+CZwFrAM8S0paR+aK7gmcSupyvibwNDAFmN3g\neqdK+jHwK2B1Us/Ak7Npj0i6HVgtIv6vzuw3kpowp5AS18Rs/tqyvy3pGeAbwEakpsppwHezIq9m\n/29I6lX4Z9K5MbMhQRG+htFsIEh6GDgvIs4pjJ8AjIqIfaqJzKz9uWZl1mKS1iDVckaTrp8ysz5y\nsjJrvWdIzZBfyjp3mFkfuRnQzMzanu8NaGZmbc/JyszM2p6TlZmZtT0nKzMza3tOVmZm1vacrMzM\nrO05WZmZWdtzsjIzs7bnZGVmZm3v/wP/lScksHPStQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xcd162e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#inputs:\n",
    "#tfdata --> log-transformed data. rows = genes. columns = replicates (several per cell type). \n",
    "#idxs --> which cell types are we interested in out of all of the cell types in tfdata?\n",
    "#labs --> cell type names\n",
    "#iunique --> what is the cell type identity of each column of tfdata?\n",
    "#iunique has same length as the # of columns in tfdata, and the max value\n",
    "#of iunique is the number of cell types.\n",
    "\n",
    "\n",
    "obj = topologies_and_genes()\n",
    "\n",
    "\n",
    "#simulate some data:\n",
    "#3 cell types A,B,C with 5 genes g1,g2,...,g4\n",
    "# g1,g2,g3 have clear minimum in B and similar values in A,C and g4 has a clear maximum in A and similar in B,C \n",
    "#for T=𝒜 and βi=1, our model is that the distribution of the expression levels of gA,B,Ci\n",
    "# in the three cell types A, B and C has the smallest mean value in either B or C but not in A\n",
    "\n",
    "nA = 10 #no. cells of type A\n",
    "nB = 20 #no. cells of type B\n",
    "nC = 13 #no. cells of type C\n",
    "idxs = np.arange(3,dtype=int)\n",
    "labs = np.array(['A','B','C'])\n",
    "iunique = np.concatenate([np.ones(nA)*idxs[0],np.ones(nB)*idxs[1],np.ones(nC)*idxs[2]])\n",
    "iunique = np.array(iunique,dtype=int)\n",
    "\n",
    "g1_A = obj.simulate_ZIP(mu_true=120, psi_true=0.8, nx=1, ny=nA)\n",
    "g1_B = obj.simulate_ZIP(mu_true=20, psi_true=0.7, nx=1, ny=nB)\n",
    "g1_C = obj.simulate_ZIP(mu_true=20, psi_true=0.75, nx=1, ny=nC)\n",
    "#min in C\n",
    "\n",
    "g2_A = obj.simulate_ZIP(mu_true=140, psi_true=0.8, nx=1, ny=nA)\n",
    "g2_B = obj.simulate_ZIP(mu_true=140, psi_true=0.7, nx=1, ny=nB)\n",
    "g2_C = obj.simulate_ZIP(mu_true=10, psi_true=0.75, nx=1, ny=nC)\n",
    "#min in C\n",
    "\n",
    "g3_A = obj.simulate_ZIP(mu_true=190, psi_true=0.8, nx=1, ny=nA)\n",
    "g3_B = obj.simulate_ZIP(mu_true=10, psi_true=0.7, nx=1, ny=nB)\n",
    "g3_C = obj.simulate_ZIP(mu_true=70, psi_true=0.75, nx=1, ny=nC)\n",
    "#min in B\n",
    "\n",
    "g4_A = obj.simulate_ZIP(mu_true=100, psi_true=0.8, nx=1, ny=nA)\n",
    "g4_B = obj.simulate_ZIP(mu_true=10, psi_true=0.7, nx=1, ny=nB)\n",
    "g4_C = obj.simulate_ZIP(mu_true=70, psi_true=0.75, nx=1, ny=nC)\n",
    "#min in B\n",
    "\n",
    "g5_A = obj.simulate_ZIP(mu_true=45, psi_true=0.8, nx=1, ny=nA)\n",
    "g5_B = obj.simulate_ZIP(mu_true=50, psi_true=0.7, nx=1, ny=nB)\n",
    "g5_C = obj.simulate_ZIP(mu_true=50, psi_true=0.75, nx=1, ny=nC)\n",
    "#min in A\n",
    "\n",
    "#plt.hist(np.ravel(g1_A),alpha=0.75)\n",
    "#plt.hist(np.ravel(g1_B),alpha=0.75)\n",
    "#plt.hist(np.ravel(g1_C),alpha=0.75)\n",
    "#plt.show()\n",
    "\n",
    "tfdata= np.vstack([np.column_stack([g1_A,g1_B,g1_C]),\n",
    "                   np.column_stack([g2_A,g2_B,g2_C]),\n",
    "                   np.column_stack([g3_A,g3_B,g3_C]),\n",
    "                   np.column_stack([g4_A,g4_B,g4_C]),\n",
    "                   np.column_stack([g5_A,g5_B,g5_C])])\n",
    "\n",
    "\n",
    "obj.set_data(tfdata, idxs, labs, iunique)\n",
    "\n",
    "#set the paramters\n",
    "obj.get_parameters();\n",
    "#modify the Params to have the right range and distribution type\n",
    "obj.Params.cutoff = 0.#no cutoff\n",
    "obj.Params.mumax = 250.\n",
    "obj.Params.mumin = 0.\n",
    "obj.Params.sigmin = 0.\n",
    "obj.Params.sigmax = 1.\n",
    "obj.Params.dist = 'P'\n",
    "\n",
    "obj.start_Sibilant()\n",
    "\n",
    "#plot last topology prob as a function of prior odds\n",
    "plt.figure(figsize=(10,7))\n",
    "for x_line in np.where(np.diff(obj.iunique))[0]:\n",
    "    plt.plot([x_line+0.5,x_line+0.5],[-1,obj.tfdata.shape[0]+1],'w')\n",
    "plt.imshow(obj.tfdata,interpolation='nearest')\n",
    "plt.xticks([np.mean(np.where(obj.iunique==idx_)[0]) for idx_ in obj.idxs],obj.labs)\n",
    "#plt.yticks([0,1,2,3],['A','B','C','null'])\n",
    "plt.xlabel('Cells')\n",
    "plt.ylabel('Genes')\n",
    "plt.colorbar(orientation='horizontal')\n",
    "plt.title('Expresion matrix')\n",
    "plt.show()\n",
    "\n",
    "#plot last topology prob as a function of prior odds\n",
    "plt.figure(figsize=(10,7))\n",
    "plt.title('Last investigated topology')\n",
    "plt.plot([obj.report_ind],[obj.tlikely_ind],'wo')\n",
    "plt.imshow(np.log(obj.pT_g).T,interpolation='nearest')\n",
    "plt.yticks([0,1,2,3],['A','B','C','null'])\n",
    "\n",
    "plt.xlabel('prior odds')\n",
    "plt.ylabel('Topology')\n",
    "plt.colorbar(orientation='horizontal')\n",
    "plt.show()\n",
    "\n",
    "#plot transition probabilities\n",
    "plt.figure()\n",
    "plt.imshow(np.log(obj.transition_probs),interpolation='nearest')\n",
    "plt.xticks([],[])\n",
    "plt.title('Transition probabilities')\n",
    "plt.ylabel('Genes')\n",
    "plt.xlabel('Combinations')\n",
    "plt.colorbar()\n",
    "\n",
    "plt.figure()\n",
    "plt.imshow(np.log(obj.marker_probs[:,:,-1]),interpolation='nearest')\n",
    "plt.xticks([],[])\n",
    "plt.title('Marker probabilities in last topology investigated')\n",
    "plt.xlabel('Cell types')\n",
    "plt.xticks([0,1,2],['A','B','C'])\n",
    "plt.ylabel('Genes')\n",
    "plt.colorbar();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}