optimizeKaDensity_sameFit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis for generating traces with preset ROI having a certain amplitude AP\n",
    "\n",
    "This file runs optimizeKaDensity. It finds the KaDensity that brings the AP amplitude in the target ROI to a requested amplitude. \n",
    "It does this for the intact cell and then uses the same parameters for the cut cell.\n",
    "\n",
    "In this run, we're saving KaDensity, ApAmp, CaAmp, InputResistance (computed with a -10pA injection into the site), and EPSP size in dendrite and soma (using an alpha synapse with properties TBD). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload\n",
    "\n",
    "import numpy as np\n",
    "from neuron import h, gui\n",
    "\n",
    "from src.collection_uncageMapping import L23\n",
    "import src.morphologyFunctions as mfx\n",
    "import src.neuronFunctions as nfx\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "from matplotlib import cm\n",
    "\n",
    "from scipy.io import savemat, loadmat\n",
    "\n",
    "import pickle\n",
    "import time\n",
    "\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from contextlib import contextmanager\n",
    "import sys, os\n",
    "\n",
    "@contextmanager\n",
    "def suppress_stdout():\n",
    "    with open(os.devnull, \"w\") as devnull:\n",
    "        old_stdout = sys.stdout\n",
    "        sys.stdout = devnull\n",
    "        try:  \n",
    "            yield\n",
    "        finally:\n",
    "            sys.stdout = old_stdout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def determineKaDensity(kaDensity, targetAmplitude, cellID, cutExperiment, naDensity, idxROI):\n",
    "    # Create cell\n",
    "    for sec in h.allsec(): h.delete_section(sec=sec)\n",
    "    with suppress_stdout():\n",
    "        cell1 = L23(cellID=cellID,cutExperiment=cutExperiment,dendNa=[naDensity,None,None,False],dendK=[kaDensity,np.inf,True,False],dxSeg=1,fixDiam=None);\n",
    "\n",
    "    # Record response of AP at all desired sites\n",
    "    stim1 = nfx.attachCC(cell1.soma, delay=1, dur=1, amp=3.5, loc=0.5) # set stim up for somatic injection\n",
    "    \n",
    "    tv = h.Vector() # Time stamp vector\n",
    "    tv.record(h._ref_t)\n",
    "    \n",
    "    vsec = h.Vector()\n",
    "    vsec.record(getattr(cell1.sectionList[idxROI](cell1.segmentList[idxROI]),'_ref_v'))\n",
    "    \n",
    "    # Simulate Data\n",
    "    nfx.simulate(tstop=15,v_init=-75,celsius=35)\n",
    "\n",
    "    # Analyze Data\n",
    "    vData = np.array(vsec)\n",
    "    apAmp = np.amax(vData)\n",
    "\n",
    "    # Reset stim program\n",
    "    stim1 = None\n",
    "    \n",
    "    return np.abs(apAmp - targetAmplitude)\n",
    "\n",
    "def saveKaResults(fname, kaDensity, apAmp, caAmp, vTraces, cTraces, tv, cellID, idxROI, silentID, ires, cutExperiment, epspAmpDend, epspAmpSoma, epspDendTraces,epspSomaTraces,tvEpsp):\n",
    "    saveList = [kaDensity, apAmp, caAmp, vTraces, cTraces, tv, cellID, idxROI, silentID, ires, cutExperiment, epspAmpDend, epspAmpSoma, epspDendTraces,epspSomaTraces,tvEpsp]\n",
    "    fid = open(fname,'wb')\n",
    "    pickle.dump(saveList, fid)\n",
    "    fid.close()\n",
    "    return None\n",
    "\n",
    "def loadKaResults(fname):\n",
    "    fid = open(fname,'rb')\n",
    "    loadedData = pickle.load(fid)\n",
    "    fid.close()\n",
    "    kaDensity=loadedData[0]\n",
    "    apAmp=loadedData[1]\n",
    "    caAmp=loadedData[2]\n",
    "    vTraces=loadedData[3]\n",
    "    cTraces=loadedData[4]\n",
    "    tv=loadedData[5]\n",
    "    cellID=loadedData[6]\n",
    "    idxROI=loadedData[7]\n",
    "    silentID=loadedData[8]\n",
    "    ires=loadedData[9]\n",
    "    cutExperiment=loadedData[10]\n",
    "    epspAmpDend=loadedData[11]\n",
    "    epspAmpSoma=loadedData[12]\n",
    "    epspDendTraces=loadedData[13]\n",
    "    epspSomaTraces=loadedData[14]\n",
    "    tvEpsp=loadedData[15]\n",
    "    return kaDensity,apAmp,caAmp,vTraces,cTraces,tv,cellID,idxROI,silentID,ires,cutExperiment,epspAmpDend,epspAmpSoma,epspDendTraces,epspSomaTraces,tvEpsp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t1 \n",
      "\t1 \n",
      "Working on cell 1/8, ROI 0\n",
      "\t1 \n",
      "\t1 \n",
      "Finished. AP Amps: [-9.99952426 20.27466232 25.56092987], CaAmps: [0.81348358 0.8216668  0.82330701], Ires: [173.74753992 217.64762957 231.63432693], EpspDend: [ 34.91987026  95.66144398 100.96854096], EpspSoma: [131.11008593 131.10940659 131.10940659]\n",
      "\t1 \n",
      "Finished. AP Amps: [34.80829089 20.18377833 25.46381891], CaAmps: [0.82467975 0.82107923 0.8228251 ], Ires: [272.38098916 235.14477571 248.22940636], EpspDend: [110.20385678  95.57894627 100.87945051], EpspSoma: [131.13144755 131.1307522  131.1307522 ]\n",
      "Working on cell 1/8, ROI 1\n",
      "\t1 \n",
      "\t1 \n",
      "Finished. AP Amps: [-57.81763893 -10.00012036  -3.4261504 ], CaAmps: [0.00475691 0.60703157 0.74090155], Ires: [170.78257251 214.37695343 228.20682777], EpspDend: [16.50222887 27.13947229 24.37357165], EpspSoma: [130.25016232 130.24862515 130.24862515]\n",
      "\t1 \n",
      "Finished. AP Amps: [28.20195116 12.01110639 17.72430743], CaAmps: [0.80903662 0.80410579 0.80609091], Ires: [268.04747124 230.73292958 243.67006312], EpspDend: [103.80085121  87.73165802  93.476373  ], EpspSoma: [130.26228688 130.26077386 130.26077386]\n",
      "Working on cell 1/8, ROI 2\n",
      "\t1 \n",
      "\t1 \n",
      "Finished. AP Amps: [-57.81791575 -16.26853233  -9.99991988], CaAmps: [0.00466115 0.47268819 0.66024497], Ires: [170.78242954 214.37681065 228.20667971], EpspDend: [16.50200039 27.13864029 24.37274462], EpspSoma: [130.2501336  130.24859635 130.24859635]\n",
      "\t1 \n",
      "Finished. AP Amps: [28.2017761  12.01090473 17.72402957], CaAmps: [0.80903604 0.80410505 0.80609014], Ires: [268.04733205 230.73278791 243.66991663], EpspDend: [103.80066422  87.73149985  93.47615514], EpspSoma: [130.26225842 130.26074536 130.26074536]\n",
      "\t1 \n",
      "Working on cell 2/8, ROI 0\n",
      "\t1 \n",
      "Finished. AP Amps: [-9.99990941 21.63560101 25.25332634 25.60505812], CaAmps: [0.61894207 0.75939664 0.75749654 0.75854934], Ires: [178.89600574 235.65644327 241.84869241 236.92510742], EpspDend: [ 62.12754935  97.51366222 101.22777383 101.48024316], EpspSoma: [124.63540501 124.63247747 124.63247747 124.63247747]\n",
      "\t1 \n",
      "Finished. AP Amps: [26.00358258 21.75532295 25.51754074 25.83705217], CaAmps: [0.75856395 0.75697372 0.75675412 0.75822591], Ires: [273.73377146 258.56327113 263.61940852 259.11159778], EpspDend: [101.79572493  97.74858773 101.42819812 101.86005565], EpspSoma: [124.63533469 124.63242672 124.63242672 124.63242672]\n",
      "Working on cell 2/8, ROI 1\n",
      "\t1 \n",
      "Finished. AP Amps: [-46.98506744  -9.99987311  16.57912212  16.93488941], CaAmps: [0.00671071 0.50044846 0.68900285 0.6911854 ], Ires: [172.62736827 229.4150239  235.99242724 231.19102581], EpspDend: [28.14039121 45.24751953 92.82271283 93.07080275], EpspSoma: [123.23702291 123.2339407  123.2339407  123.2339407 ]\n",
      "\t1 \n",
      "Finished. AP Amps: [17.42788447 10.6468772  17.1925178  17.56100522], CaAmps: [0.69255549 0.67841679 0.69024405 0.69269899], Ires: [266.74423252 251.00936738 256.55715597 252.16665701], EpspDend: [93.69832912 86.64003957 93.3699524  93.77659103], EpspSoma: [123.22444039 123.22135467 123.22135467 123.22135467]\n",
      "Working on cell 2/8, ROI 2\n",
      "\t1 \n",
      "Finished. AP Amps: [-57.16033944 -47.47009399 -10.00044903  -3.57162251], CaAmps: [0.0023515  0.00617309 0.38841833 0.48510367], Ires: [165.5603445  222.41359967 229.52253309 224.86100772], EpspDend: [19.08662963 28.8186346  60.88600746 70.70692124], EpspSoma: [121.65638691 121.65207241 121.65207241 121.65207241]\n",
      "\t1 \n",
      "Finished. AP Amps: [ -1.17665818 -43.52458255  -0.42038707   1.53214544], CaAmps: [0.51544443 0.00922479 0.52761699 0.54831884], Ires: [259.10588194 242.74381644 248.82818999 244.57448595], EpspDend: [73.78283997 32.52446643 74.75879019 77.2350377 ], EpspSoma: [121.62704695 121.62272644 121.62272644 121.62272644]\n",
      "Working on cell 2/8, ROI 3\n",
      "\t1 \n",
      "Finished. AP Amps: [-57.53764098 -48.01856381 -22.3424012   -9.99982839], CaAmps: [0.00221369 0.00568894 0.12139053 0.3681013 ], Ires: [165.01500547 221.87009901 229.01818133 224.36810405], EpspDend: [18.76551974 28.32225406 50.60752995 62.3689425 ], EpspSoma: [121.52656566 121.5221555  121.5221555  121.5221555 ]\n",
      "\t1 \n",
      "Finished. AP Amps: [ -6.5501038  -44.62837012  -5.41074859  -1.70150338], CaAmps: [0.4346458  0.00844275 0.46171954 0.50930632], Ires: [258.4981854  242.08538931 248.21341989 243.97089337], EpspDend: [66.84227304 31.51166499 68.12659236 73.41637178], EpspSoma: [121.4961197  121.49169711 121.49169711 121.49169711]\n"
     ]
    }
   ],
   "source": [
    "numCells = 8\n",
    "naDensity = 5\n",
    "initKaDensity = 0.01\n",
    "kaBounds = [(0,None)]\n",
    "roiType = False # True means silent, False means active (will get all from each category)\n",
    "targetAmplitude = -10\n",
    "\n",
    "cellID = []\n",
    "cutExperiment = []\n",
    "idxROI = []\n",
    "silentID = []\n",
    "tv = []\n",
    "vTraces = []\n",
    "cTraces = []\n",
    "\n",
    "kaDensity = []\n",
    "apAmp = []\n",
    "caAmp = []\n",
    "\n",
    "ires = []\n",
    "\n",
    "epspAmpDend = []\n",
    "epspAmpSoma = []\n",
    "epspDendTraces = []\n",
    "epspSomaTraces = []\n",
    "tvEpsp = []\n",
    "\n",
    "for n in range(2): # Comment numCells back in to run them all! range(numCells):\n",
    "    # Create cell just to get silent IDs\n",
    "    for sec in h.allsec(): h.delete_section(sec=sec)\n",
    "    cell1 = L23(cellID=n,cutExperiment=0,dendNa=[naDensity,None,None,False],dendK=[initKaDensity,np.inf,True,False],dxSeg=1,fixDiam=None);\n",
    "    cSilentID = cell1.silentID\n",
    "    \n",
    "    # Trade out these two lines to use all or just active/silent type of interest\n",
    "    listTarget = [n for n in range(len(cSilentID))]\n",
    "    #listTarget = [n for n in range(len(cSilentID)) if cSilentID[n]==roiType]\n",
    "    \n",
    "    for r in listTarget:\n",
    "        print(f'Working on cell {n+1}/{numCells}, ROI {r}')\n",
    "        \n",
    "        # -- do it for the normal cell -- \n",
    "        \n",
    "        # Optimize kaDensity for this cell\n",
    "        results = minimize(determineKaDensity, initKaDensity, args=(targetAmplitude,n,0,naDensity,r),method='Nelder-Mead',bounds=kaBounds)\n",
    "        kaDensity.append(results.x) # Store optimal kaDensity\n",
    "\n",
    "        # Run cell with optimal kaDensity and record voltage/calcium traces for all ROIs in the cell\n",
    "        for sec in h.allsec(): h.delete_section(sec=sec)\n",
    "        cell1 = L23(cellID=n,cutExperiment=0,dendNa=[naDensity,None,None,False],dendK=[results.x,np.inf,True,False],dxSeg=1,fixDiam=None);\n",
    "        \n",
    "        # Record response of AP at all desired sites\n",
    "        stim1 = nfx.attachCC(cell1.soma, delay=1, dur=1, amp=3.5, loc=0.5) # set stim up for somatic injection\n",
    "        \n",
    "        # Record peak of AP in all the sites\n",
    "        vsec,ctv = mfx.recordSites(cell1.sectionList,cell1.segmentList)\n",
    "\n",
    "        # Record ica in all sites + soma\n",
    "        csec = mfx.recordSites(cell1.sectionList,cell1.segmentList,recordVariable='_ref_ica')[0]\n",
    "\n",
    "        # Simulate Data\n",
    "        nfx.simulate(tstop=15,v_init=-75,celsius=35)\n",
    "\n",
    "        # Convert calcium current to conductance\n",
    "        gca_sec = []\n",
    "        for ica,v in zip(csec,vsec):\n",
    "            gca_sec.append(nfx.conductanceFromCurrent(ica,v,cell1.Eca))\n",
    "        \n",
    "        # Store Data\n",
    "        cellID.append(n)\n",
    "        cutExperiment.append(0)\n",
    "        idxROI.append(r)\n",
    "        silentID.append(cSilentID[r])\n",
    "        vData = np.array(vsec)\n",
    "        gcaData = np.array(gca_sec)\n",
    "        apAmp.append(np.amax(vData,axis=1))\n",
    "        caAmp.append(np.amax(gcaData,axis=1))\n",
    "        vTraces.append(vData)\n",
    "        cTraces.append(gcaData)\n",
    "        tv.append(np.array(ctv))\n",
    "        \n",
    "        # --------- And also measure input resistance for sites!! ---------\n",
    "        stim = None\n",
    "        amplitude=-0.01\n",
    "        vsection,ctv,stim = mfx.injectSites(cell1.sectionList,cell1.segmentList,stim=stim,amplitude=amplitude)\n",
    "\n",
    "        # Reset stim program\n",
    "        stim = None\n",
    "\n",
    "        # Delay is 50ms, duration is 50ms\n",
    "        dvm = (np.array(vsection)[:,np.where(np.array(ctv)<=100)[0][-1]] - np.array(vsection)[:,np.where(np.array(ctv)<=50)[0][-1]])\n",
    "        ires.append(dvm/amplitude)        \n",
    "        \n",
    "        # --------- And also measure EPSP size and transfer function for sites! ---------\n",
    "        syn = None\n",
    "        stim = None\n",
    "        onset=50\n",
    "        tau=1\n",
    "        gmax=0.0005\n",
    "        tstop = 80\n",
    "        epspDendrite,epspSoma,epspTV,syn = mfx.injectAlphaSites(cell1.sectionList,cell1.segmentList,syn=syn,onset=onset,tau=tau,gmax=gmax,tstop=tstop)\n",
    "        vEpspDend = np.array(epspDendrite)\n",
    "        vEpspSoma = np.array(epspSoma)\n",
    "        cTvEpsp = np.array(epspTV)\n",
    "        cEpspAmpDend = np.amax(vEpspDend,axis=1) - vEpspDend[:,np.where(cTvEpsp>=onset-1)[0][0]]\n",
    "        cEpspAmpSoma = np.amax(vEpspSoma,axis=1) - vEpspSoma[:,np.where(cTvEpsp>=onset-1)[0][0]]\n",
    "        epspAmpDend.append(cEpspAmpDend)\n",
    "        epspAmpSoma.append(cEpspAmpSoma)\n",
    "        epspDendTraces.append(vEpspDend)\n",
    "        epspSomaTraces.append(vEpspSoma)\n",
    "        tvEpsp.append(cTvEpsp)\n",
    "        \n",
    "        print(f'Finished. AP Amps: {apAmp[-1]}, CaAmps: {caAmp[-1]}, Ires: {ires[-1]}, EpspDend: {epspAmpDend[-1]}, EpspSoma: {epspAmpSoma[-1]}')\n",
    "\n",
    "        # Reset stim programs\n",
    "        syn = None\n",
    "        stim = None\n",
    "        \n",
    "        \n",
    "        # -- now measure it again for the cut experiment -- (using the same fit parameters)\n",
    "        kaDensity.append(results.x) # Store optimal kaDensity\n",
    "\n",
    "        # Run cell with optimal kaDensity and record voltage/calcium traces for all ROIs in the cell\n",
    "        for sec in h.allsec(): h.delete_section(sec=sec)\n",
    "        cell1 = L23(cellID=n,cutExperiment=2,dendNa=[naDensity,None,None,False],dendK=[results.x,np.inf,True,False],dxSeg=1,fixDiam=None);\n",
    "        \n",
    "        # Record response of AP at all desired sites\n",
    "        stim1 = nfx.attachCC(cell1.soma, delay=1, dur=1, amp=3.5, loc=0.5) # set stim up for somatic injection\n",
    "        \n",
    "        # Record peak of AP in all the sites\n",
    "        vsec,ctv = mfx.recordSites(cell1.sectionList,cell1.segmentList)\n",
    "\n",
    "        # Record ica in all sites + soma\n",
    "        csec = mfx.recordSites(cell1.sectionList,cell1.segmentList,recordVariable='_ref_ica')[0]\n",
    "\n",
    "        # Simulate Data\n",
    "        nfx.simulate(tstop=15,v_init=-75,celsius=35)\n",
    "\n",
    "        # Convert calcium current to conductance\n",
    "        gca_sec = []\n",
    "        for ica,v in zip(csec,vsec):\n",
    "            gca_sec.append(nfx.conductanceFromCurrent(ica,v,cell1.Eca))\n",
    "        \n",
    "        # Store Data\n",
    "        cellID.append(n)\n",
    "        cutExperiment.append(1)\n",
    "        idxROI.append(r)\n",
    "        silentID.append(cSilentID[r])\n",
    "        vData = np.array(vsec)\n",
    "        gcaData = np.array(gca_sec)\n",
    "        apAmp.append(np.amax(vData,axis=1))\n",
    "        caAmp.append(np.amax(gcaData,axis=1))\n",
    "        vTraces.append(vData)\n",
    "        cTraces.append(gcaData)\n",
    "        tv.append(np.array(ctv))\n",
    "        \n",
    "        # And also measure input resistance for sites!!\n",
    "        stim = None\n",
    "        amplitude=-0.01\n",
    "        vsection,ctv,stim = mfx.injectSites(cell1.sectionList,cell1.segmentList,stim=stim,amplitude=amplitude)\n",
    "\n",
    "        # Delay is 50ms, duration is 50ms\n",
    "        dvm = (np.array(vsection)[:,np.where(np.array(ctv)<=100)[0][-1]] - np.array(vsection)[:,np.where(np.array(ctv)<=50)[0][-1]])\n",
    "        ires.append(dvm/amplitude)        \n",
    "        \n",
    "        # --------- And also measure EPSP size and transfer function for sites! ---------\n",
    "        syn = None\n",
    "        stim = None\n",
    "        onset=50\n",
    "        tau=1\n",
    "        gmax=0.0005\n",
    "        tstop = 80\n",
    "        epspDendrite,epspSoma,epspTV,syn = mfx.injectAlphaSites(cell1.sectionList,cell1.segmentList,syn=syn,onset=onset,tau=tau,gmax=gmax,tstop=tstop)\n",
    "        vEpspDend = np.array(epspDendrite)\n",
    "        vEpspSoma = np.array(epspSoma)\n",
    "        cTvEpsp = np.array(epspTV)\n",
    "        cEpspAmpDend = np.amax(vEpspDend,axis=1) - vEpspDend[:,np.where(cTvEpsp>=onset-1)[0][0]]\n",
    "        cEpspAmpSoma = np.amax(vEpspSoma,axis=1) - vEpspSoma[:,np.where(cTvEpsp>=onset-1)[0][0]]\n",
    "        epspAmpDend.append(cEpspAmpDend)\n",
    "        epspAmpSoma.append(cEpspAmpSoma)\n",
    "        epspDendTraces.append(vEpspDend)\n",
    "        epspSomaTraces.append(vEpspSoma)\n",
    "        tvEpsp.append(cTvEpsp)\n",
    "        \n",
    "        print(f'Finished. AP Amps: {apAmp[-1]}, CaAmps: {caAmp[-1]}, Ires: {ires[-1]}, EpspDend: {epspAmpDend[-1]}, EpspSoma: {epspAmpSoma[-1]}')\n",
    "\n",
    "        # Reset stim programs\n",
    "        syn = None\n",
    "        stim = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "timeStamp = time.strftime(\"optKaDensity_%Y%b%d_%H%M%S_sameFits\")\n",
    "fname = './'+timeStamp+'.pkl'\n",
    "saveKaResults(fname, kaDensity, apAmp, caAmp, vTraces, cTraces, tv, cellID, idxROI, silentID, ires, cutExperiment, epspAmpDend, epspAmpSoma, epspDendTraces, epspSomaTraces, tvEpsp)\n",
    "\n",
    "# Stack saved results in format matlab will like\n",
    "numFits = len(kaDensity)\n",
    "maxROI = 4\n",
    "NT = tv[0].shape[0]\n",
    "matVoltage = np.empty((NT,numFits,maxROI))\n",
    "matVoltage[:] = np.NAN\n",
    "matCalcium = np.empty_like(matVoltage)\n",
    "matCalcium[:] = np.NAN\n",
    "matTv = tv[0]\n",
    "matApAmp = np.empty((numFits,maxROI))\n",
    "matApAmp[:] = np.NAN\n",
    "matCaAmp = np.empty_like(matApAmp)\n",
    "matCaAmp[:] = np.NAN\n",
    "matIres = np.empty_like(matApAmp)\n",
    "matIres[:] = np.NAN\n",
    "matKaDensity = np.empty_like(matApAmp)\n",
    "matCellID = np.empty(numFits)\n",
    "matIdxROI = np.empty(numFits)\n",
    "matCutExp = np.empty(numFits)\n",
    "matEpspDend = np.empty_like(matApAmp)\n",
    "matEpspSoma = np.empty_like(matApAmp)\n",
    "matEpspDend[:] = np.NAN\n",
    "matEpspSoma[:] = np.NAN\n",
    "NT = tvEpsp[0].shape[0]\n",
    "matEpspDendTraces = np.empty((NT,numFits,maxROI))\n",
    "matEpspSomaTraces = np.empty((NT,numFits,maxROI))\n",
    "matEpspDendTraces[:] = np.NAN\n",
    "matEpspSomaTraces[:] = np.NAN\n",
    "for n in range(numFits):\n",
    "    cNumROI = vTraces[n].shape[0]\n",
    "    matVoltage[:,n,:cNumROI] = vTraces[n].T\n",
    "    matCalcium[:,n,:cNumROI] = cTraces[n].T\n",
    "    matApAmp[n,:cNumROI] = apAmp[n]\n",
    "    matCaAmp[n,:cNumROI] = caAmp[n]\n",
    "    matIres[n,:cNumROI] = ires[n]\n",
    "    matEpspDend[n,:cNumROI] = epspAmpDend[n]\n",
    "    matEpspSoma[n,:cNumROI] = epspAmpSoma[n]\n",
    "    matKaDensity[n] = kaDensity[n]\n",
    "    matCellID[n] = cellID[n]\n",
    "    matIdxROI[n] = idxROI[n]\n",
    "    matCutExp[n] = cutExperiment[n]\n",
    "    matEpspDendTraces[:,n,:cNumROI] = epspDendTraces[n].T\n",
    "    matEpspSomaTraces[:,n,:cNumROI] = epspSomaTraces[n].T\n",
    "\n",
    "matname = './'+timeStamp+'.mat'\n",
    "matdict = {\"matVoltage\":matVoltage, \"matCalcium\":matCalcium, \"matApAmp\":matApAmp,\"matCaAmp\":matCaAmp,\"matIres\":matIres,\"matKaDensity\":matKaDensity,\"matCellID\":matCellID,\"matIdxROI\":matIdxROI,\\\n",
    "           \"matTV\":matTv,\"matCutExp\":matCutExp,\"silentID\":np.array(silentID),\"epspAmpDend\":matEpspDend,\"epspAmpSoma\":matEpspSoma,\"epspDendTraces\":matEpspDendTraces,\"epspSomaTraces\":matEpspSomaTraces,\"tvEpsp\":tvEpsp}\n",
    "savemat(matname,matdict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "fname = './optKaDensity_2021Aug31_152955_sameFits.pkl'\n",
    "kaDensity,apAmp,caAmp,vTraces,cTraces,tv,cellID,idxROI,silentID,ires,cutExperiment,epspAmpDend,epspAmpSoma,epspDendTraces,epspSomaTraces,tvEpsp = loadKaResults(fname)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Everything below this is just for playing with the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[True,\n",
       " True,\n",
       " False,\n",
       " False,\n",
       " False,\n",
       " False,\n",
       " True,\n",
       " True,\n",
       " False,\n",
       " False,\n",
       " False,\n",
       " False,\n",
       " False,\n",
       " False,\n",
       " True,\n",
       " True,\n",
       " False,\n",
       " False,\n",
       " False,\n",
       " False,\n",
       " True,\n",
       " True,\n",
       " False,\n",
       " False,\n",
       " True,\n",
       " True,\n",
       " False,\n",
       " False,\n",
       " True,\n",
       " True,\n",
       " False,\n",
       " False,\n",
       " True,\n",
       " True,\n",
       " False,\n",
       " False,\n",
       " True,\n",
       " True,\n",
       " False,\n",
       " False]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "silentID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "'bool' object is not subscriptable",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-8-c87b2b7a3ca2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     10\u001b[0m     \u001b[0mnumROI\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvTraces\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mr\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnumROI\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 12\u001b[0;31m         \u001b[0;32mif\u001b[0m \u001b[0msilentID\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0;32mFalse\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0midxROI\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m!=\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;32mcontinue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     13\u001b[0m         \u001b[0mccol\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'b'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     14\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0msilentID\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mccol\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'k'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: 'bool' object is not subscriptable"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAFpCAYAAABj6bgoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAO70lEQVR4nO3cX4jld3nH8c/TXQP1T1XMKjZ/MC3RuBem6BilaBsrrUl6EQQvEsXQICyhRrw0FKoX3tSLgojRZZEg3piLGjSWaCgUa8GmzQY0GiWyjTTZRshGxUKEhtWnFzOWYZzd+c3kzGx8eL1gYH6/850zz3yZfeeXM+ec6u4AMMPvXOgBAFgdUQcYRNQBBhF1gEFEHWAQUQcYZMeoV9VdVfVUVX3vHLdXVX2qqk5V1cNV9cbVjwnAEkuu1D+f5Lrz3H59kis3Po4l+exzHwuAvdgx6t39zSQ/Pc+SG5N8odc9kORlVfXqVQ0IwHKreEz9kiRPbDo+vXEOgAN2eAX3Uduc2/a9B6rqWNYfosmLXvSiN1111VUr+PYAszz00ENPd/eRvXztKqJ+Osllm44vTfLkdgu7+0SSE0mytrbWJ0+eXMG3B5ilqv5rr1+7iodf7k1yy8azYN6a5Ofd/eMV3C8Au7TjlXpVfTHJtUkurqrTST6W5AVJ0t3Hk9yX5IYkp5L8Ismt+zUsAOe3Y9S7++Ydbu8kH1zZRADsmVeUAgwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMMiiqFfVdVX1aFWdqqo7trn9pVX11ar6TlU9UlW3rn5UAHayY9Sr6lCSO5Ncn+Rokpur6uiWZR9M8v3uvjrJtUn+vqouWvGsAOxgyZX6NUlOdfdj3f1skruT3LhlTSd5SVVVkhcn+WmSsyudFIAdLYn6JUme2HR8euPcZp9O8vokTyb5bpIPd/evVjIhAIstiXptc663HL8rybeT/H6SP0ry6ar6vd+4o6pjVXWyqk6eOXNm18MCcH5Lon46yWWbji/N+hX5ZrcmuafXnUryoyRXbb2j7j7R3WvdvXbkyJG9zgzAOSyJ+oNJrqyqKzb++HlTknu3rHk8yTuTpKpeleR1SR5b5aAA7OzwTgu6+2xV3Z7k/iSHktzV3Y9U1W0btx9P8vEkn6+q72b94ZqPdPfT+zg3ANvYMepJ0t33Jblvy7njmz5/MslfrHY0AHbLK0oBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBhF1gEFEHWAQUQcYRNQBBlkU9aq6rqoerapTVXXHOdZcW1XfrqpHqupfVjsmAEsc3mlBVR1KcmeSP09yOsmDVXVvd39/05qXJflMkuu6+/GqeuV+DQzAuS25Ur8myanufqy7n01yd5Ibt6x5b5J7uvvxJOnup1Y7JgBLLIn6JUme2HR8euPcZq9N8vKq+kZVPVRVt2x3R1V1rKpOVtXJM2fO7G1iAM5pSdRrm3O95fhwkjcl+csk70ryt1X12t/4ou4T3b3W3WtHjhzZ9bAAnN+Oj6ln/cr8sk3HlyZ5cps1T3f3M0meqapvJrk6yQ9XMiUAiyy5Un8wyZVVdUVVXZTkpiT3blnzlSRvr6rDVfXCJG9J8oPVjgrATna8Uu/us1V1e5L7kxxKcld3P1JVt23cfry7f1BVX0/ycJJfJflcd39vPwcH4DdV99aHxw/G2tpanzx58oJ8b4Dns6p6qLvX9vK1XlEKMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIIuiXlXXVdWjVXWqqu44z7o3V9Uvq+o9qxsRgKV2jHpVHUpyZ5LrkxxNcnNVHT3Huk8kuX/VQwKwzJIr9WuSnOrux7r72SR3J7lxm3UfSvKlJE+tcD4AdmFJ1C9J8sSm49Mb5/5fVV2S5N1Jjp/vjqrqWFWdrKqTZ86c2e2sAOxgSdRrm3O95fiTST7S3b883x1194nuXuvutSNHjiydEYCFDi9YczrJZZuOL03y5JY1a0nurqokuTjJDVV1tru/vJIpAVhkSdQfTHJlVV2R5L+T3JTkvZsXdPcVv/68qj6f5B8FHeDg7Rj17j5bVbdn/Vkth5Lc1d2PVNVtG7ef93F0AA7Okiv1dPd9Se7bcm7bmHf3Xz33sQDYC68oBRhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYBBRBxhE1AEGEXWAQUQdYJBFUa+q66rq0ao6VVV3bHP7+6rq4Y2Pb1XV1asfFYCd7Bj1qjqU5M4k1yc5muTmqjq6ZdmPkvxpd78hyceTnFj1oADsbMmV+jVJTnX3Y939bJK7k9y4eUF3f6u7f7Zx+ECSS1c7JgBLLIn6JUme2HR8euPcuXwgydeey1AA7M3hBWtqm3O97cKqd2Q96m87x+3HkhxLkssvv3zhiAAsteRK/XSSyzYdX5rkya2LquoNST6X5Mbu/sl2d9TdJ7p7rbvXjhw5spd5ATiPJVF/MMmVVXVFVV2U5KYk925eUFWXJ7knyfu7+4erHxOAJXZ8+KW7z1bV7UnuT3IoyV3d/UhV3bZx+/EkH03yiiSfqaokOdvda/s3NgDbqe5tHx7fd2tra33y5MkL8r0Bns+q6qG9Xhh7RSnAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDiDrAIKIOMIioAwwi6gCDLIp6VV1XVY9W1amqumOb26uqPrVx+8NV9cbVjwrATnaMelUdSnJnkuuTHE1yc1Ud3bLs+iRXbnwcS/LZFc8JwAJLrtSvSXKqux/r7meT3J3kxi1rbkzyhV73QJKXVdWrVzwrADtYEvVLkjyx6fj0xrndrgFgnx1esKa2Odd7WJOqOpb1h2eS5H+r6nsLvv90Fyd5+kIPcYHZg3X2YZ19SF631y9cEvXTSS7bdHxpkif3sCbdfSLJiSSpqpPdvbaraQeyD/bg1+zDOvuwvgd7/dolD788mOTKqrqiqi5KclOSe7esuTfJLRvPgnlrkp9394/3OhQAe7PjlXp3n62q25Pcn+RQkru6+5Gqum3j9uNJ7ktyQ5JTSX6R5Nb9GxmAc1ny8Eu6+76sh3vzueObPu8kH9zl9z6xy/VT2Qd78Gv2YZ19eA57UOs9BmACbxMAMMi+R91bDCzag/dt/OwPV9W3qurqCzHnfttpHzate3NV/bKq3nOQ8x2UJftQVddW1ber6pGq+peDnnG/Lfg38dKq+mpVfWdjD0b+na6q7qqqp8719O499bG79+0j639Y/c8kf5DkoiTfSXJ0y5obknwt6891f2uSf9/PmQ76Y+Ee/HGSl298fv20PVi6D5vW/XPW/4bzngs99wX6fXhZku8nuXzj+JUXeu4LsAd/k+QTG58fSfLTJBdd6Nn3YS/+JMkbk3zvHLfvuo/7faXuLQYW7EF3f6u7f7Zx+EDWn+c/zZLfhST5UJIvJXnqIIc7QEv24b1J7unux5Oku6ftxZI96CQvqapK8uKsR/3swY65/7r7m1n/2c5l133c76h7i4Hd/3wfyPp/mafZcR+q6pIk705yPHMt+X14bZKXV9U3quqhqrrlwKY7GEv24NNJXp/1FzF+N8mHu/tXBzPe88qu+7joKY3PwcreYuC32OKfr6rekfWov21fJ7owluzDJ5N8pLt/uX6BNtKSfTic5E1J3pnkd5P8W1U90N0/3O/hDsiSPXhXkm8n+bMkf5jkn6rqX7v7f/Z7uOeZXfdxv6O+srcY+C226Oerqjck+VyS67v7Jwc020Fasg9rSe7eCPrFSW6oqrPd/eWDGfFALP038XR3P5Pkmar6ZpKrk0yJ+pI9uDXJ3/X6A8unqupHSa5K8h8HM+Lzxq77uN8Pv3iLgQV7UFWXJ7knyfsHXY1tteM+dPcV3f2a7n5Nkn9I8tfDgp4s+zfxlSRvr6rDVfXCJG9J8oMDnnM/LdmDx7P+fyqpqldl/Q2uHjvQKZ8fdt3Hfb1Sb28xsHQPPprkFUk+s3GVeraHvaHRwn0Yb8k+dPcPqurrSR5O8qskn+vuMe9ouvB34eNJPl9V3836QxAf6e5x79xYVV9Mcm2Si6vqdJKPJXlBsvc+ekUpwCBeUQowiKgDDCLqAIOIOsAgog4wiKgDDCLqAIOIOsAg/wcmKmg9Y+iulwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1440x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot Somatic AP\n",
    "fig = plt.figure(figsize=(20,6))\n",
    "\n",
    "cmap = cm.get_cmap('jet')\n",
    "\n",
    "numFits = len(kaDensity)\n",
    "\n",
    "plt.subplot(1,3,1)\n",
    "for n in range(numFits):\n",
    "    numROI = vTraces[n].shape[0]\n",
    "    for r in range(numROI):\n",
    "        if silentID[n][r]==False and idxROI[n]!=r: continue\n",
    "        ccol = 'b'\n",
    "        if silentID[n][r]==False: ccol = 'k'\n",
    "        plt.plot(np.array(tv[n]),vTraces[n][r].T,c=ccol)\n",
    "\n",
    "plt.subplot(1,3,2)\n",
    "for n in range(numFits):\n",
    "    numROI = cTraces[n].shape[0]\n",
    "    for r in range(numROI):\n",
    "        if silentID[n][r]==False and idxROI[n]!=r: continue\n",
    "        ccol = 'b'\n",
    "        if silentID[n][r]==False: ccol = 'k'\n",
    "        plt.plot(np.array(tv[n]),cTraces[n][r].T,c=ccol)\n",
    "\n",
    "plt.subplot(1,3,3)\n",
    "for n in range(numFits):\n",
    "    numROI = cTraces[n].shape[0]\n",
    "    for r in range(numROI):\n",
    "        if silentID[n][r]==False and idxROI[n]!=r: continue\n",
    "        ccol = 'b'\n",
    "        if silentID[n][r]==False: ccol = 'k'\n",
    "        plt.scatter(apAmp[n][r],caAmp[n][r],c=ccol,s=80)\n",
    "    plt.plot([apAmp[n][0],apAmp[n][idxROI[n]]],[caAmp[n][0],caAmp[n][idxROI[n]]],c='k',linewidth=0.2,linestyle='dashed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cellID = 0\n",
    "cutExperiment = 0\n",
    "idxROI = 0\n",
    "\n",
    "naDensity = 6\n",
    "targetAmplitude = -20\n",
    "initKaDensity = 0.01\n",
    "\n",
    "results = minimize(determineKaDensity, initKaDensity, args=(targetAmplitude, cellID, cutExperiment, naDensity, idxROI),method='Nelder-Mead',bounds=(0,None))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this cell, I'll go through the cells and vary the potassium channel density until the requested ROI has an AP amplitude of a set voltage. \n",
    "\n",
    "numCells = 8 # number of cells to run through\n",
    "naDensity = 6 # Na channel density (in channels / um2)\n",
    "initKaDensity = 0.01 # K channel density (in units of S/cm2)\n",
    "kaMinimum = 0 # minimum k channel density \n",
    "kaMaximum = 0.2 # maximum k channel density\n",
    "stepSize = 0.00001 # ka S/cm2 per mV\n",
    "tolerance = 0.05 # tolerance in mV for AP\n",
    "\n",
    "# Index of ROI to set AP amplitude for\n",
    "targetROI = [0,0,0,0,0,0,0,0]\n",
    "targetAmplitude = 10 # target AP amplitude\n",
    "\n",
    "cAP = np.Inf\n",
    "cKa = initKaDensity\n",
    "while np.abs(cAP - targetAmplitude) > tolerance:\n",
    "    # Create cell\n",
    "    for sec in h.allsec(): h.delete_section(sec=sec)\n",
    "    cell1 = L23(cellID=1,cutExperiment=0,dendNa=[naDensity,None,None,False],dendK=[cKa,np.inf,True,False],dxSeg=1,fixDiam=None);\n",
    "\n",
    "    # Record response of AP at all desired sites\n",
    "    stim1 = nfx.attachCC(cell1.soma, delay=1, dur=1, amp=3.5, loc=0.5) # set stim up for somatic injection\n",
    "\n",
    "    # Record peak of AP in all the sites\n",
    "    vsec,tv = mfx.recordSites(cell1.sectionList,cell1.segmentList)\n",
    "\n",
    "    # Record ica in all sites + soma\n",
    "    csec = mfx.recordSites(cell1.sectionList,cell1.segmentList,recordVariable='_ref_ica')[0]\n",
    "\n",
    "    # Simulate Data\n",
    "    nfx.simulate(tstop=8,v_init=-75,celsius=35)\n",
    "\n",
    "    gca_sec = []\n",
    "    for ica,v in zip(csec,vsec):\n",
    "        gca_sec.append(nfx.conductanceFromCurrent(ica,v,cell1.Eca))\n",
    "\n",
    "    # Analyze Data\n",
    "    vData = np.array(vsec)\n",
    "    gcaData = np.array(gca_sec)\n",
    "    apAmp = np.amax(vData,axis=1)\n",
    "    gcAmp = np.amax(gcaData,axis=1)\n",
    "\n",
    "    # Reset stim program\n",
    "    stim1 = None\n",
    "\n",
    "    cAP = apAmp[targetROI[0]]\n",
    "    print(cAP)\n",
    "    \n",
    "    voltageError = (apAmp[targetROI[0]] - targetAmplitude)\n",
    "    if voltageError>1: \n",
    "        updateValue = voltageError**2 * np.sign(voltageError) * stepSize\n",
    "    else:\n",
    "        updateValue = voltageError * stepSize\n",
    "    newKa = cKa + updateValue\n",
    "    if newKa > kaMaximum or newKa < kaMinimum:\n",
    "        print(f'Out of range!! Ka:{newKa}')\n",
    "        break\n",
    "    cKa = newKa\n",
    "\n",
    "    print(cKa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t1 \n",
      "-19.862718519098134\n"
     ]
    }
   ],
   "source": [
    "testKa = 0.01871557\n",
    "for sec in h.allsec(): h.delete_section(sec=sec)\n",
    "cell1 = L23(cellID=0,cutExperiment=0,dendNa=[naDensity,None,None,False],dendK=[testKa,np.inf,True,False],dxSeg=1,fixDiam=None);\n",
    "\n",
    "# Record response of AP at all desired sites\n",
    "stim1 = nfx.attachCC(cell1.soma, delay=1, dur=1, amp=3.5, loc=0.5) # set stim up for somatic injection\n",
    "\n",
    "# Record peak of AP in all the sites\n",
    "vsec,tv = mfx.recordSites(cell1.sectionList,cell1.segmentList)\n",
    "\n",
    "# Record ica in all sites + soma\n",
    "csec = mfx.recordSites(cell1.sectionList,cell1.segmentList,recordVariable='_ref_ica')[0]\n",
    "\n",
    "# Simulate Data\n",
    "nfx.simulate(tstop=8,v_init=-75,celsius=35)\n",
    "\n",
    "gca_sec = []\n",
    "for ica,v in zip(csec,vsec):\n",
    "    gca_sec.append(nfx.conductanceFromCurrent(ica,v,cell1.Eca))\n",
    "\n",
    "# Analyze Data\n",
    "vData = np.array(vsec)\n",
    "gcaData = np.array(gca_sec)\n",
    "apAmp = np.amax(vData,axis=1)\n",
    "gcAmp = np.amax(gcaData,axis=1)\n",
    "\n",
    "# Reset stim program\n",
    "stim1 = None\n",
    "\n",
    "print(apAmp[targetROI[0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t1 \n"
     ]
    }
   ],
   "source": [
    "testKa = 0.01871557\n",
    "for sec in h.allsec(): h.delete_section(sec=sec)\n",
    "cell1 = L23(cellID=0,cutExperiment=0,dendNa=[naDensity,None,None,False],dendK=[testKa,np.inf,True,False],dxSeg=1,fixDiam=None);\n",
    "\n",
    "stim = None\n",
    "syn = None\n",
    "onset=30\n",
    "tau=2\n",
    "gmax=0.0002\n",
    "tstop = 100\n",
    "vsection,vsoma,tv,syn = mfx.injectAlphaSites(cell1.sectionList,cell1.segmentList,syn=syn,onset=onset,tau=tau,gmax=gmax,tstop=tstop)\n",
    "vsec = np.array(vsection)\n",
    "vsoma = np.array(vsoma)\n",
    "tv = np.array(tv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7AAAAFlCAYAAADS0QR3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABoN0lEQVR4nO3deXxU1f3/8deZyb4RAoGwJkDCJiibAi6gKLbu1n2tu1Vr1brUpd9aW/uz1lpba61WW+uuoKhg3TcUFUEQZF8iewIhkH1PZs7vj5kJARIIZJKbmXk/H4955M6du3wyjty855x7jrHWIiIiIiIiItLZuZwuQERERERERKQ1FGBFREREREQkJCjAioiIiIiISEhQgBUREREREZGQoAArIiIiIiIiIUEBVkREREREREJClNMFHIzu3bvbrKwsp8sQEZEwsXDhwh3W2nSn6whlujaLiEgwtXRtDskAm5WVxYIFC5wuQ0REwoQxZqPTNYQ6XZtFRCSYWro2qwuxiIiIiIiIhAQFWBEREREREQkJCrAiIiIiIiISEhRgRUREREREJCQowIqIiIiIiEhIaFOANcZMM8Ys9j82GGMW+9cf0WT998aYn7Sw/7PGmPVNth3VlnpEREREREQkfLVpGh1r7fmBZWPMX4BS/9NlwDhrbYMxphfwvTHmbWttQzOHucNa+3pb6hAREREREZHwF5R5YI0xBjgPmAJgra1q8nIcYINxHhEREdk3Y8w0YIj/aSpQYq0dZYzJAlYCq/2vfWOtva6Z/dOAaUAWsAE4z1pb3L5Vi4iItE6w7oE9Biiw1q4NrDDGjDfGLAeWAte10PoK8P+MMUuMMX81xsQGqR4REZGIZK0931o7ylo7CpgBvNHk5R8CrzUXXv3uAj6x1uYAn/ifi4iIdAr7DbDGmI+NMcuaeZzRZLMLgVea7metnWetPQQ4HLjbGBPXzOHvBob6t0kD7txHHdcaYxYYYxYUFha24lcTERGJXE16R72yv233cAbwnH/5OeDMIJYlIiLSJvvtQmytPWFfrxtjooCzgLEt7L/SGFMBjAAW7PHaVv9irTHmv8Dt+6jjKeApgHHjxqlLsoiIyL7t1TsKGGCMWQSUAf9nrZ3TzH49m1yftwE927lOERGRVgtGF+ITgFXW2i2BFcaYAf5gizEmE18r64Y9d/QP8BT4lvhMfIM/ibSLNWvWsHHjRqfLEBFps4PsHbUV6G+tHQ3cCrxsjEnZ13mstZZ9jGOh3lHSkXbs2MHs2bMpKChwuhQRcVAwAuwF7N096Wh8Iw8vBt4EbrDW7gAwxrxrjOnt3+4lY8xSfPfJdgf+EIR6RPbyxhtvMGzYMAYPHsyXX37pdDkiIm1irT3BWjuimcdM2K131LQm+9Raa3f6lxcCPwCDmzl8QZMvmHsB2/dRx1PW2nHW2nHp6enB+wVF9vD888/Tv39/jjvuOPr378+f/vQnfN+viEikafMoxNbay5tZ9wLwQgvbn9xkeUpbzy+yP16vl9tuu42cnByqqqq46aabWLhwIb6GfxGRsNRc76h0oMha6zHGDARygHXN7DsLuAx40P9zZgfUK9KiTz/9lMsvv5zJkydzxx138N///pe77vKNLXbnnS0OnyIiYSoo0+iIdGbz589nw4YNvPDCC9TU1HDNNdfw9ddfc9RRRzldmohIe2mud9Qk4PfGmHrAi2+GgCIAY8y/gSettQvwBdfpxpirgI34BoIScURtbS1XX3012dnZvPPOOyQkJPDjH/+Yiy++mHvuuYdJkyYxceJEp8sUkQ6kACthb+bMmURFRXHqqafidrv5+c9/zowZMxRgRSRstdA7aga+aXWa2/7qJss7gePbrTiRA/DSSy+xfv163nvvPRISEgBwuVw8/fTTfPnll9xwww18++23REXpT1qRSBGseWBFOq2vv/6acePGkZqaSnJyMlOnTmXmTPWIExER6cystfz5z39m9OjR/OhHP8Lj8fDxxx8zffp0ampqeOSRR1i8eDEvvvii06WKSAdSgJWw5vF4WLhwIYcffnjjuqlTp7Ju3To2b97sYGUiIiKyLwsWLGDVqlX8/Oc/p7i4mOOOO46pU6dy/vnnM2jQINxuN6NGjeKBBx7A4/E4Xa6IdBAFWAlrq1evprKyknHjxjWumzx5MgCff/65U2WJiIjIfrzyyitER0dz5plncvHFFzNv3jyefvpp5s6dy9ChQzn//PM59dRTWbt2La+99prT5YpIB1GAlbC2ePFiAMaOHdu4buTIkaSmpirAioiIdFLWWqZPn86Pf/xjPv74Y95//30eeeQRrr76aiZMmMBHH33EkCFDeOqpp8jMzOSf//yn0yWLSAdRgJWwtmbNGowxZGdnN65zu91MmDCB+fPnO1iZiIiItGTFihXk5eVx2mmn8X//938ceuihXHfddY2vp6Sk8Oqrr1JUVETPnj2ZM2cOy5Ytc7BiEekoCrAS1nJzc+nfvz+xsbG7rR87dizLly+npqbGocpERESkJR9//DEAMTEx5Obmcvfdd+N2u3fbZsSIEdx4440sWLCAqKgonnzySSdKFZEOpgArYW3t2rXk5OTstX7MmDF4PB6WLFniQFUiIiKyLx9//DE5OTm89dZbpKenc9ZZZzW73W9+8xuSkpLo3bs3L7/8MrW1tR1cqYh0NAVYCWu5ubm7dR8OGDNmDADfffddR5ckIiIi+1BfX8/s2bOZNGkS7777LpdccgkxMTHNbpuWlsb111/Pli1bKC4u5v333+/gakWkoynAStgqKiqiqKio2RbYzMxM0tLSFGBFREQ6mWXLllFRUUFiYiJ1dXUttr4G3HTTTURFRREXF8fLL7/cQVWKiFMUYCVsrV+/HoABAwbs9ZoxhlGjRvH99993dFkiIiKyD4FBFtevX096ejoTJ07c5/a9e/fmggsuwOPxMHPmTMrKyjqiTBFxiAKshK28vDwA+vbt2+zrw4cPZ+XKlVhrO7IsERER2Yd58+bRrVs3vvrqK0455ZS9Bm9qzjXXXEN9fT21tbW89dZb7V+kiDhGAVbCViDA9unTp9nXhw8fTnl5eeN2IiIi4rz58+czbNgwioqKOO6441q1z1FHHcWwYcOIiYlhxowZ7VyhiDhJAVbCVl5eHm63m549ezb7+rBhwwDfXHMiIiLivPLyclasWEFycjIAkyZNatV+xhiuvvpq6urqeP/996msrGzPMkXEQQqwErby8vLIyMhosevR8OHDAQVYERGRzuL777/HWktpaSn9+/cnKyur1ftedNFFuFyuxhArIuFJAVbCVl5e3m73v27btm23+eHS09Pp1q0bK1eudKI8ERER2cOyZcsAWLNmDZMnTz6gfTMyMjj22GNxuVy88cYb7VGeiHQCCrAStvLy8hrvf33sscfo1asXQ4cOJT8/H/B1Nxo+fLhaYEVERDqJZcuWkZSUxI4dO5gwYcIB73/RRRfh9XqZOXMm9fX17VChiDhNAVbCViDA5ufn86tf/YrDDz+c7du3c/vttzduM2zYMJYvX66RiEVERDqBpUuXNn75PHbs2APe/6yzziIqKorKykpmz54d5OpEpDNQgJWwVFNTQ2lpKRkZGbz66qvU1NTwwgsvcOONN/Lqq6+yefNmAIYOHUpxcTE7d+50uGIREZHIZq1l2bJlxMfH43a7OfTQQw/4GF27duXEE08E0HQ6ImFKAVbCUmFhIeC7z/X1119n1KhRDBkyhJ/97GdYa3nppZcAyMnJASA3N9exWkVERMQ3VkVRURE1NTUccsghxMfHH9RxLrzwQsAXYNXDSiT8KMBKWAoE2JSUFObPn8/JJ58MwMCBA5kwYULj4A7Z2dmAAqyIiIjTli5dCviC7MF0Hw445ZRTcLlc5Ofns2bNmmCVJyKdhAKshKVAgN25cycej4eJEyc2vnbSSSexYMECduzYwYABA3C5XKxdu9apUkVERARYvnw5ACUlJW0KsF27dm287r/77rtBqU1EOg8FWAlLO3bsAGDDhg0AjB8/vvG1E088EWstn3zyCbGxsfTv318tsCIiIg5bu3YtSUlJAAd1/2tTF1xwAQDTp09vc10i0rkowEpYCrTA5ubmkpmZSXp6euNrhx9+OCkpKXz22WeArxuxAqyIiIiz1q5dS9euXQEYPnx4m451xhlnADB//nzKy8vbXJuIdB4KsBKWCgsLcbvdrFu3bq+LoNvtZvz48XzzzTeAAqyIiEhnkJubS0xMDD169KBbt25tOla/fv0YPHgwXq+XTz75JEgVikhnoAArYamwsJDu3buzZs0ahgwZstfrEydOZOnSpZSXl5OTk0NRURFFRUUOVCoiIiK1tbVs2rSJ2traNre+Blx00UUAvP7660E5noh0DgqwEpYKCwtJTU2lurq62QA7YcIEvF4vCxYs0EjEIiIiDlu3bh1er5eioqKgBdizzjoLgHfeeUfT6YiEEQVYCUuFhYXExcUBMHTo0L1eDwzqNHfuXAVYERERhwWuwVVVVUELsCNGjCAtLY2SkpLGKXpEJPQpwEpY2rFjBy6X7+PdXAtsWloaAwcOZNGiRQwcOBBjjAKsiIiIQ5pOZ7e/AFtdDR99BO++C6WlLW9njOGUU04B4O233w5KnSLiPAVYCUs7duzAWktsbCwZGRnNbnPYYYfx/fffExcXR79+/TQXrIiIiENyc3OJj48HYNiwYS1u9/77MGAAnHginHIK9O4NjzwCXm/z25977rkAvPbaa0GvWUScoQArYcdaS0lJCXV1dfTt2xdjTLPbHXbYYeTm5lJZWcnAgQNZv359B1cqIiIi4GuBTU5OJikpiZ49eza7zUcfwamnQs+e8M478NlncPzxcNttcMMNzYfYKVOm4Ha7WbJkCWVlZe38W4hIR1CAlbBTUVGBx+OhqqqKfv36tbjdYYcdhrWWZcuWMWDAAAVYERERh/zwww9ERUUxaNCgZr943r4dLrgAhg2DL7+Ek0+GY4+FmTPhrrvgX/+C//f/9j5uYmIiY8aMwVrLp59+2v6/iIi0OwVYCTvFxcUAlJWV0bdv3xa3O+ywwwD4/vvvycrKIj8/n5qamg6pUURERHw8Hg9btmyhtraWQYMGNbvNr38NZWUwfTokJ+9abww88ABccgn89rfQ3JSvF154IQDTpk1rj/JFpIMpwErYKSkpafy5rxbYrKwsUlJS+P777xkwYAAAGzdu7IgSRURExG/btm3U19dTWlrabIBdvRr+8x+46SZfC+yejPG1wObkwLXXQlXV7q+fdtppAHzwwQeaTkckDCjAStgJBFiv17vPAGuM4dBDD90twKobsYiISMcKfHnc0NDQOLVdU48+CjExcOedLR8jIQGeegrWrYP779/9tezsbHr27ElxcTErV64MZuki4gAFWAk7gS7EwD4DLPjmiFuxYgVZWVmAAqyIiEhHa9r7ac8W2OJieO45uPhi6NFj38eZPBl++lP4619h8+bdXzv11FMBmDlzZlBqFhHnKMBK2Am0wML+A+ywYcMoLi7G7XYTGxvLhg0b2rc4ERER2c2+Auxrr/m6BP/85607VqD19d57d19//vnnA7oPViQcKMBK2GnaAturV699bjt06FAAVq9eTWZmplpgRUREOtimTZuIi4sjOjp6ry+eX30VhgyB0aNbd6z+/eEXv/C12q5atWv9pEmTiI6OZunSpVRUVASxehHpaAqwEnaatsB269Ztn9sGJktftWoVWVlZCrAiIiIdbOPGjcTFxZGVlYXb7W5cn58Ps2f7ps9pYUr3Zv3qVxAXBw89tGtdbGwshx9+OF6vl88++yx4xYtIh1OAlbBTXFxMTEwM6enpu10Im9O3b18SExNZuXKl5oIVERFxwMaNG7HW7tV9+M03wVpfgD0Q6elw9dXw4ou73wt70UUXAfDyyy+3tWQRcZACrISdkpISoqKi6LG/0R7wjUQ8dOjQxgC7c+dOysvLO6BKERERsdayceNGampqGDhw4G6vvf8+DBwI/rt9Dshtt/nC71//umvd6aefDsCHH36o6XREQpgCrISdkpISjDGtCrDg60a8atUqTaUjIiLSwUpKSqioqKC2tpbMzMzG9XV18Nln8OMfH9xxMzPhwgt9U+uUlvrW9evXj969e1NUVMTatWuDUL2IOEEBVsJOcXExXq/3gALs5s2b6dmzJ6AAKyIi0lGajkDcv3//xuWvvoLKSvjRjw7+2Dff7DvG88/vWnfaaacB8NZbbx38gUXEUQqwEnZKSkpoaGhodYANjERcX18PoKl0REREOkjTANt0BOIPPoCoKDjuuIM/9tixMH48PP64rzsx7LoP9pVXXjn4A4uIoxRgJewUFxdTX19/QC2wAFu3biUxMVEtsCIiIh1kc5NRlpq2wM6eDRMmQHJy245/442wejV88onv+cSJE4mJiWHp0qVUV1e37eAi4ggFWAk7gWl0WhtgBw0ahNvtZvXq1RqJWEREpAPl5eVhjMHtdjfO3V5dDd99B0cd1fbjn3uub1Tif/zD9zw6Oprx48fj8XiYPXt2208gIh1OAVbCitfrbZygvLUBNiYmhqysLNauXasAKyIi0oHy8/OJj4+nT58+REVFAbBgAdTXByfAxsbCNdfA22/Dpk2+dZdccgkAL774YttPICIdTgFWwkplZWXjcmsDLEBOTs5uAVbD64uIiLS/vLw8oqKidrv/9euvfT8nTgzOOa65BrxeeO453/PAdDoffPBBcE4gIh1KAVbCStM5XNPT01u9X3Z2Nrm5uWRlZVFRUcHOnTvbozwRERFpIj8/H4/Hs9cIxIMHQ/fuwTlHVhYcfzw884wvyGZkZNCvXz927tzJunXrgnMSEekwCrASVpoG2LS0tFbvl5OTQ3l5eeM+6kYsIiLS/vLy8qiurm4MsNb6WmCD0X24qauugg0bfHPLwq7pdGbMmBHcE4lIu2tzgDXGTDPGLPY/NhhjFu/xen9jTIUx5vYW9h9gjJlnjMn1HyumrTVJ5CorK2tcTk1NbfV+OTk5AI1dhzWVjoiISPuqrKykrKwMr9fb2IV43TrYudM3AnEwnXkmpKb6WmFh132wr776anBPJCLtrs0B1lp7vrV2lLV2FDADeGOPTR4B3tvHIf4E/NVamw0UA1e1tSaJXIEW2MTERNxud6v3CwTYwD20CrAiIiLtKz8/v3E50AK7aJHv+ZgxwT1XfDxcfDHMmAHFxXDEEUcQFxfHkiVLqK2tDe7JRKRdBa0LsTHGAOcBrzRZdyawHli+j32mAK/7Vz0HnBmsmiTyBAJsSkrKAe2XmZmJ2+0mPz+frl27qguxiISslnpGGWOyjDHVTV57soX9/2yMWWWMWWKMedMYk9qR9UvkyMvLa1xuGmDdbhgxIvjnu/JKqK2Fl18Gt9vNhAkTaGho4Isvvgj+yUSk3QTzHthjgAJr7VoAY0wScCfwu33s0w0osdY2+J9vAfoEsSaJMIEAeyDdh8E3L9yAAQNYu3YtWVlZaoEVkZC1n55RPwRes9Ze18IhPgJGWGsPBdYAd7dvxRKpWmqBHT4c4uKCf74xY2DUqF3diH/6058C8Oyzzwb/ZCLSbloVYI0xHxtjljXzOKPJZhfSpPUVuA9f1+CKYBRqjLnWGLPAGLOgsLAwGIeUMBS4B/ZABnAK2HMqHRGRUNZcz6jWsNZ+2OSL5W+AvsGuTQR2BdjExMTGL54XLYLRo9vvnFdeCd99B99/v2s6nY8++qj9TigiQdeqAGutPcFaO6KZx0wAY0wUcBYwrclu44GHjDEbgFuAe4wxN+5x6J1Aqn9/8F0k82iGtfYpa+04a+24A5keRSJLoAW2+0GMvZ+Tk0Nubi6ZmZls2LBBc8GKSKjbrWeU3wBjzCJjzOfGmGNacYwr2cc4FvpyWdoiLy8Pt9tN7969McawbRts29a+AfaiiyAmBp59Frp160ZmZiaFhYVs3ry5/U4qIkEVrC7EJwCrrLVbAiustcdYa7OstVnA34AHrLX/aLqT9SWEz4Bz/KsuA2YGqSaJQG0JsNnZ2VRUVNC9e3dqamooKCgIdnkiIkFxkD2jtgL9rbWjgVuBl40xLQ4YYIz5NdAAvNTSNvpyWdoiPz+fmJgYevXqBcDixb717Rlgu3WD00+HF1+EurpdrbDTp09vv5OKSFAFK8BewAF0UTLGvGuM6e1/eidwqzEmF989sf8JUk0SgdrahRhoHL1Y98GKSGd1MD2jrLW11tqd/uWFwA/A4OaOb4y5HDgVuNiqO4q0k8AgThkZGcCuEYhHjWrf815+OezYAe++C5dddhkAr7xyQD3tRcRBUfvfZP+stZfv5/X79nh+cpPldcARwahDpKSkBICuXbse8L6BAFtXVwfA+vXrmRDsiehERDrGXj2jjDHpQJG11mOMGQjkAOv23NEY82PgV8Bka21VRxUskSc/P5/6+vrGALtsGfTvD126tO95f/QjyMiA//4X3nxzNPHx8Xz//ffU1dURExPTvicXkTYL5ijEIo4rKioCDi7AZmZmEhUVRWlpKaAWWBEJac31jJoELPFPq/M6cJ21tgjAGPNvY8w4/3b/AJKBj/Y13Y5IW1hrycvLo6GhobEL8YoVcMgh7X/uqCi49FJ45x0oLHQxceJEGhoamDNnTvufXETaTAFWwkpxcTFwcAE2KiqKgQMHsnHjRtLT0xVgRSRkWWsvt9Y+uce6GdbaQ/xT6Iyx1r7d5LWrrbUL/MvZ1tp+rZhuR+SgFRUVNfZ4ysjIwOOBVatg2LCOOf/ll4PHAy+9tKsb8X//+9+OObmItElQuhCLdBaB1tODCbDgG8gpMBesptIRERFpH03ngO3VqxcbN0JNjW8O2I4wfDgccYSvG/GcOb6xzzSdjkhoUAushJXAKMQHG2ADU+kMGDBALbAiIiLtpOlI/xkZGaxY4VvuqAALvlbYZctg3bouDBgwgO3bt+8WrEWkc1KAlbBSWVkJtK0FtrKykvT0dDZu3IjX6w1meSIiIgJs3769cblpgO2oLsQAF1wAsbG+VtjAdDoajVik81OAlbBSXV0NtK0FFiAuLo66ujq2bt0atNpERETEJ9AC63a76d69OytWQK9ekJracTV07QpnngkvvwxXXum71full1qc9lhEOgkFWAkb1lpqamoA6HKQY/BnZ2cDNLa86j5YERGR4CsoKMAYQ48ePXC73axc2bHdhwOuuAKKimDNmqF06dKFpUuX0tDQ0PGFiEirKcBK2KipqcHr9RIdHU1U1MGNTxaYSqeqyjf1oe6DFRERCb6CggJiYmLIyMjAWt8UOk4E2BNOgD594Nln4bjjjqOhoUGDOYl0cgqwEjYC97/GxcUd9DGioqLIysqisLAQUAusiIhIeygoKMDlctGrVy+2bIGKCmcCrNsNP/0pvPceXHDBzwH417/+1fGFiEirKcBK2AgE2ISEhDYdJzs7mw0bNpCRkaEWWBERkXZQUFCA1+slIyOD1at964YMcaaWyy4Drxc2bJhCTEwMs2fPdqYQEWkVBVgJG8EKsE2n0lELrIiISPBt27aNuro6MjIyyM31rfOPo9jhhgyBiRPhuedcjB07ltLSUlauXOlMMSKyXwqwEjYCATYpKalNx8nOzqasrEwtsCIiIu3AWsv27dux1tKrVy9ycyEuDnr3dq6mK66AlSth6tSfAfDwww87V4yI7JMCrISNwMBLwQiwAMnJyWzatEmjEYqIiARRSUlJ47U1IyODtWth0CBwOfhX6XnnQXw8bN16IS6Xi//973/OFSMi+6QAK2Ej0AJ7sFPoBAQCrNvtxuPxkJeX1+baRERExCcwByzQ2ALrVPfhgC5d4KyzYPr0GIYNO4Tt27ezZcsWZ4sSkWYpwErYCFaAzcrKwuVyUVdXB2gkYhERkWBqGmB79Mjghx/A/92xoy6/HEpLYdSonwLw6KOPOluQiDRLAVbCRiDApqamtuk4MTExZGVlUVpaCmguWBERkWBqGmC93gxqaztHgJ0yBfr1g61br8EYw/Tp050uSUSaoQArYSMQONPS0tp8rOzsbPLz8zHGqAVWREQkiAIBNikpifz8RMD5LsTguwf3sstg9uwuZGYOYtOmTezcudPpskRkDwqwEjaKioqA4AXYH374gT59+qgFVkREJIgCAbZXr16sXetb1xlaYMHXjdjrhUGDLgXgb3/7m6P1iMjeFGAlbBQXFwPQrVu3Nh8rOzub0tJS+vbtqwArIiISRAUFBURHRzfOARsbC337Ol2Vz6BBcMwxsGHDzQC8/PLLDlckIntSgJWwEehC3NZ7YAFy/H2Zunbtqi7EIiIiQVRQUIAxpnEE4oEDnZ1CZ0+XXw4//NCFPn0GsW7dusYvyEWkc+hE/1yItE1ZWRngm7+1OdbCn/8MRx8Nf/qTr4tQSwJT6cTFxZGXl9c4IrGIiIi0zfbt2/F4PI0tsJ3h/temzj0XEhIgI+MyAB5++GGHKxKRphRgJWyUl5cDkJKS0uzrjzwCv/oVFBTAXXfB3Xe3fKwBAwZgjMHj8eD1etm8eXN7lCwiIhJxtm7disfjoWdPX4DtLPe/BiQnwznnwKpVt2CM4cUXX3S6JBFpQgFWwkYgwDbXAltUBPfeC6edBmvWwDXX+FpjFyxo/lixsbH079+/cWoe3QcrIiLSdtbaxkGc4uN7UV3tu++0s7niCqisTCYjYzCbNm1i27ZtTpckIn4KsBI2AmGzuRbY556Dqir4wx/AGHj4YejWDe67r+XjZWdnNw6fr/tgRURE2q6iooLa2lrANwcs+O6B7WwmTYKsLEhMvBqAP/3pT84WJCKNFGAlbFRVVQHNB9jnn4cJE+DQQ/FvAzfdBO+8A8uXN3+8nJwcNm3ahNvtVgusiIhIEARaXwFqa30BNivLoWL2weXyDeaUm3s9xrh49dVXnS5JRPwUYCVs1NTUAJCYmLjb+vx8WLwYzjxz9+2vvx6io+GZZ5o/XnZ2NkVFRfTp00ctsCIiIkHQNMBWVPQCIDPTqWr27ac/BUikZ89D2bZtG7m5uU6XJCIowEoYqa2tJSoqCtceY/G//77v50kn7b599+5wyinw0kvQ0LD38QIjEaenp6sFVkREJAgCAdblcrF9e3d69ID4eIeLasGAAXDcceDx3ALA73//e2cLEhFAAVbCSG1tLdHR0Xut//RT6NkTRo7ce5/LLvONSvzJJ3u/FgiwSUlJaoEVEREJgu3btwO+L4c3bXJ3yu7DTV1+ORQWXkpUVDQzZ850uhwRQQFWwkh9fT2xsbF7rZ8/33f/qzF77/PjH0NiIrz11t6vDfSPKuFyudi6dWtjF2URERE5OIEW2F69erFxY+e8/7Wps8+GpCQX6emTKSsrY/bs2U6XJBLxFGAlbDQ0NBAXF7fbupISWLsWDj+8+X3i4uBHP4JZs8Da3V+Lj4+nX79+1NXVAbBx48Z2qFpERCRyFBQU4Ha7ycjwBdjOev9rQGIinHceFBf/DoDf/e53DlckIgqwEhY8Hg9er3evFtjAPK8tBViA00/3DfS0cOHer2VnZ1NaWgpoLlgREZG2KigowBhDamovams7fwsswLXXQk3NkcTHd+XLL7+kvr7e6ZJEIpoCrISFwByw8XuMBBEIsOPGtbzvKaf4hsv/3//2fi07O5v8/HxAc8GKiIi01bZt2/B4PMTGdt4pdPY0frzvi/CYmHNpaGjgiSeecLokkYimACthIRBgExISdlu/fDn06wdpaS3v2707jB0LH3+892uBqXSio6PVAisiItJG+fn5WGsBX4Dt7F2IA37xCygtvQ+Av//9784WIxLhFGAlLFRVVQF7zwG7ahUMHbr//Y8/HubNg/Ly3dfn5OQAkJGRoRZYERGRNgqMQlxf37nngN3TeedBjx69SEwcwg8//MCmTZucLkkkYinASlgItMAmJSU1rrO29QH2hBN8c8F+8cXu6wNT6aSmpqoFVkREpA2qq6sbv3CurMyge3doctnu1GJjfffCVlbeCcBdd93lcEUikUsBVsJCcwE2Px8qKloXYI880ndx2nM+2MBUOrGxsQqwIiIibRBofQUoKsoImdbXgOuuA7f7ctzuON58801/V2gR6WgKsBIWysrKAEhOTm5ct2qV72drAmx8PBx99N73wSYmJtK7d2+8Xi/bt29vDMoiIiJyYAJzwPqWM0JiAKem+vSBs882uFxnUlNTwzPPPON0SSIRSQFWwsLOnTsB6NKlS+O6QIAdNqx1xzj+eFi6FAoLd1+fnZ3dGFzXrVvX5lpFREQiUSDAJiQksHlzUsi1wALcdBPU1/8ZgD/+8Y8OVyMSmRRgJSyUlJQAuwfYtWt999ZkZLTuGJMm+X5+/fXu63NyctixYwcAubm5bS1VREQkIgUCbPfuPamuDo0pdPZ05JEwYUJfoqOH8cMPP+jvAhEHKMBKWAh0IW4aYDds8I1uaEzrjjFunO8+2C+/3H19dnZ2YwuvLlQiIiIHJxBgU1N7A6EZYI2BO++E+vrfAvDLX/7S4YpEIo8CrISF5gLsxo0HNjx/bKxvovLmAixA165dFWBFREQOUkFBAS6Xi4QEX4ANxS7EAKefDjk55+FypfD+++9TW1vrdEkiEUUBVsJCIMB27dq1cd2BBljwDeS0cCH4R/kHdgXY7t27K8CKiIgcpEALrNsdWnPA7snlgl/9yuD1XkVDQwO//e1vnS5JJKIowEpYqKioACAtLQ2A8nIoLj64AFtfD99+u2vdoEGDAN+IxAqwIiIiByc/Px+v10tDQwZdu0KTTlMh59JLoUePBzDGzRNPPKEpdUQ6kAKshIXAKMGBLsQbN/rWH2iAPfJI38+m3YiTk5PJyMjAGMPmzZupqalpa7kiIiIRZ+vWrQBUV4feHLB7io2FW2+Nw9pTKCsr4/nnn3e6JJGIoQArYSEQYBMTE4GDD7Bdu8IhhzR/H2x1dTXWWk2lIyIichC2b98OQElJr5AcwGlP110HycmPAXDvvfc6XI1I5FCAlbAQrAALvm7EX38NHs+udRqJWERE5ODV19dTXl4OwPbtGWERYLt0gdtv7w+MZdOmTXzxxRdOlyQSERRgJSxU+UddahpgY2JaPwdsU0cdBWVlsGLFrnXZ2dkUFhYCCrAiIiIHKtD6ClBTE/pdiANuvhmSkv4JwI033uhwNSKRQQFWwkLgvtRAgN20Cfr29Y0UeKAmTPD9/OabXetycnIASElJUYAVERE5QIERiI1xAelh0QILvlbYO+88AhjM0qVLWbJkidMliYQ9BVgJCzU1NRhjcLvdAGzdCr17H9yxsrMhLW33ABuYSqdHjx4KsCIiIgcoEGATE7sC7rAJsAA33QRJSY8DcOWVVzpcjUj4U4CVsFBbW4urSXPr1q3Qq9fBHcsYXyvsvHm71g0ePBiAhIQEBVgREZEDtG3bNgASE3sCoTsHbHNSUuDuu08ABrBw4UKWL1/udEkiYU0BVsJCXV0dUVFRjc/bEmABxo/33QNbWup7npSURN++fbHWsnHjRurq6tpYsYiISOQItMDGxPQlJQVSU52tJ9h+8QtITfWNSHzFFWqFFWlPbQqwxphpxpjF/scGY8ziPV7vb4ypMMbc3sL+zxpj1jc5xqi21CORq2mArayE8vK2BdgJE8Ba+PbbXeuGDRtGRUUFXq+XDRs2tK1gERGRCFJQUIAxBuhDZqavt1M4SU6GP/7xFCCLb7+dr3thRdpRmwKstfZ8a+0oa+0oYAbwxh6bPAK8t5/D3BE4hrV2cVvqkchVX19PdHQ04Gt9hbYF2COO8P1seh/s0KFDG79BVjdiERGR1tu2bRvWWmprw2MO2OZcfTX06/c0AJdccqnD1YiEr6B0ITa+r9TOA15psu5MYD2gGwGk3TU0NBATEwMEJ8CmpsKwYbvfBzt06NDG6XoUYEWks2qpd5QxJssYU93ktSf3c5zbjDHWGNO9QwqXsLZlyxYAysrCYw7Y5kRFwT/+cQIwlKVLl/Dpp586XZJIWArWPbDHAAXW2rUAxpgk4E7gd63Y9/8ZY5YYY/5qjIkNUj0SYTweD7Gxvo+Pf5yIg5oDtqkJE3wtsNb6ng8bNgzwDeS0du3ath1cRKSd7Kd31A9Nej1d19IxjDH9gBOBTe1brUSKrf5vl8NpDtjmnHYajB79MgCXXXaFw9WIhKf9BlhjzMfGmGXNPM5ostmFNGl9Be4D/mqtrdjP4e8GhgKHA2n4Qm9LdVxrjFlgjFlQWFi4v7IlwjQNsMFogQXfQE47dsC6db7nQ4cOBaBbt24KsCLS6TXXO+oA/BX4FWCDWpRErF1/u/UK6wBrDDz99GjgKLZs2cSzzz7rdEkiYWe/AdZae4K1dkQzj5kAxpgo4CxgWpPdxgMPGWM2ALcA9xhjbmzm2FutTy3wX+CIfdTxlLV2nLV2XHp6+oH8jhIBrLXExcUBvgAbFQXdurXtmBMm+H4GuhFnZGSQkpJCXFwcq1atatvBRUTa3269o/wGGGMWGWM+N8Yc09xO/i+o86y13+/vBPpyWVqjoaGBsrIy/7NeDBjgaDntbuxYuOSSVwHDjTfehNfrdbokkbASjC7EJwCrrLVbAiustcdYa7OstVnA34AHrLX/2HNHY0wv/08DnAksC0I9EmGstXsF2IwMcLXx033IIZCYuGsgJ2MMw4YNw+PxsHHjxsb7YUVEOtpB9o7aCvS31o4GbgVeNsak7HHcBOAe4N7W1KEvl6U1dv9yI3wHcWrq0Uf7Eht7KZWV5fziFzc5XY5IWAlGgL2AA+ieZIx51xjT2//0JWPMUmAp0B34QxDqkQhTW1sL+O5NhbbPARsQFQWHH773SMTFxcUA6kYsIo45mN5R1tpaa+1O//JC4Adg8B6HHgQMAL7396LqC3xnjGnjqAISyQIj+EdFJZCYmEBamsMFdYC0NHjssX8DCTzxxBONg1iJSNu1OcBaay+31rY4kqG19j5r7cNNnp9src33L0+x1o70X3QvacU9syJ7qajwfWwCAbawEHr0CM6xx4+HxYuhpsb3vGmAVTdiEenE9uodZYxJN8a4/csDgRxgXdOdrLVLrbU9mvSi2gKMsdZu67jSJdxs84+uGB2dTlZW+M0B25Krropm0KDHsdbLKaec6XQ5ImEjWKMQizhm586dACQmJgK+gZe6B2nShwkToL4eFi3yPQ+MRGyMUYAVkc6sud5Rk4Al/ml1Xgeus9YWARhj/m2MGdexJUqkCLTAGhMZ3YcDXC6YOfNyYChLlixk2rRp+9tFRFpBAVZCXklJCQBJSUlY62uBDdatWOPH+34GuhEHRiLu3r07q1evDs5JRESCrLneUdbaGdbaQ/xT6Iyx1r7d5LWrrbULmjlOlrV2R0fULOErEGDr6/tHVIAF33gat932DmC47LIrqQl06RKRg6YAKyEvEGATExOpqvJ19w1WC2yvXpCZuSvADhw4kKioKJKTk9UCKyIi0gqBLsT19f3Cegqdljz44EB69LiZ2toqTjnlJ06XIxLyFGAl5JWWlgK+Ftgd/naCYA6GOWHCrgAbHR1NTk4OxhhWr16tofFFRET2Y9cARpHVhTggKgo++eQRoDeffvo+s2a9vd99RKRlCrAS8gIBNjk5mcBI/cFqgQVfgN20CfLzfc+HDh1KRUUFVVVVGlVQRERkP3ZdK3tHZIAFGDHCcNtt7wOGc8+9QFPxibSBAqyEvPLycsAXYAMtsMEMsBMn+n4GWmGHDx/eOKed7oMVERHZt0AXYugVkV2IAx56aCQDBtxGXV0VEyee4HQ5IiFLAVZCXiDApqSkNLbABrML8ahREBMDc+f6no8YMaKx67DugxUREdm3wGwBsbG9gnp9DjUuF3zzzZ9xuwezZMlc7r//T06XJBKSFGAl5DUNsO3RAhsbC2PH7mqBHTlyJOCbd1YBVkREpGUNDQ2UlZUBkJnZK2LmgG1Jjx4wffqXQAz33nsPS5YsdbokkZCjACshr6KiAoDU1FR27AC3G7p0Ce45JkyABQugrg5ycnKIjo4mNTVVAVZERGQfAt2HjYlh4MBkh6vpHM46K50LL3wZ8DJ+/DGaWkfkACnASsgLBNguXbpQWOhrfXUF+ZM9YYJvep4lSyAmJoYhQ4ZgjFGAFRER2Yf8wAiIdCMrK8KbX5t46aWzGTDgGmpqShkx4kinyxEJKQqwEvICI/mlpaWxY0dwuw8H7DmQ08iRI6msrCQ/P79xHloRERHZXSDAWpsRsSMQN8cYWLr0KeLiDuOHHxZx9tlXOV2SSMhQgJWQV11dDezqQtweAbZvX+jde/eBnALBdfny5cE/oYiISBjY1QLbTwF2D4mJsGjRNxjTlTfeeIY//OExp0sSCQkKsBLyAi2wgS7E7THCoTG+Vtg9B3ICWLZsWfBPKCIiEgby8vL8SwMjegqdlgwdGsebb34LxPCb39zMyy/PcrokkU5PAVZCXqAF1u12t1sLLPjug123DrZv97XAAsTFxSnAioiItGDDhg3+pUwGDXKyks7rjDMG8ec/fwgYLr74LD7//BunSxLp1BRgJeTV1tZijMHjgaKi9mmBhd3vg83MzCQpKYnU1FSWLtUQ+CIiIs1Zv349APHx/drtC+ZwcPvtk7npplcBD8cdN4n58/W3hUhLFGAl5NXU1GCMobQUvF5IS2uf84wZA1FRvvtgXS4XhxxyCMYYli1bhrW2fU4qIiISwgJdiPv16xfxc8Duz6OPnsullz6OtfVMmHC4QqxICxRgJeTV1tbidrspLvY979q1fc4THw+jR+9+H2xpaSk7d+6koKCgfU4qIiISwnbs2AHA4MH9HK4kNDz//A1cdNFfsbaWCRPG8fnnC50uSaTTUYCVkFdXV4fL5Wr3AAu++2Dnz4eGBt99sIEBpHQfrIiIyO5qamr810k3hxzS0+lyQsZLL93C5Zc/jrV1HHvsBF599WOnSxLpVBRgJeTV19cTFRXVIQF24kSoqoJly+DQQw9tXK8AKyIisrutW7f6l7qSk6M/OQ/Ef/97A7fe+izg4cILT+S++55yuiSRTkP/mkjICwRY/7Ss7d4CC777YA877DAAEhMTFWBFRET2sGsKnQyNQHwQ/vKXy3j88Q+AKH73u59x+uk3O12SSKegACshr6Ghgejo6A5pgc3Kgp49fffBpqWlkZWVpQArIiLSjE2bNvmXMsnOdrSUkHXDDVP55JMluFzJvP323+nT5wgqKqqdLkvEUQqwEvL2DLCpqe13LmN8rbBz5/qejx49mtraWpYvX47X622/E4uIiISYdevWAeB2D6F3b4eLCWFTpgxl69Y8UlOHk5//LV279uLjjxc5XZaIYxRgJeR5vV5iYmIoLoboaEhIaN/zTZwIa9fCjh2+AFtaWkpFRQUbN25s3xOLiIiEkFWrVgGQnj4Il/7ibJMePZIpKlrOMcf8jIaGUqZOHctZZ93pdFkijtA/JxLyPB4PsbGxFBf7ug+39zxzRx3l+/nVV74AG7B48eL2PbGIiEgIWbNmDQADBw50uJLwYAx88cWTPP742xgTx5tvPkRy8kAWLVrndGkiHUoBVkKetZbY2FhKStr3/teAww+H2FiYMwfGjBkDgMvl4rvvvmv/k4uIiISITZs2AzBqlEZwCqYbbjiVnTsL6d9/MhUV6xkzJocpU26goUG3MklkUICVkGetJS4urrEFtr3FxsIRR/gCbK9evejRowcpKSksWqT7UURERMB3bd6xoxAwjB+f6XQ5Yadr10Q2bpzNgw++issVz2efPUFcXFf++MdpTpcm0u4UYCWkNTQ0ABAfH99hARbgmGNg4UKorDSMHj0aY4wCrIiIiN/27dvxeDxAVw49NMbpcsLWnXeeT1VVMRMnXorHU84991xAUlI2L774hdOlibQbBVgJaeXl5cCuANueIxA3NWkSeDy+6XQCAznl5+dTUFDQMQWIiIh0Yhs2bPAv9WPIECcrCX+xsdF8/fXzLF78A716jaOy8gcuvXQyXbsewssvf+l0eSJBpwArIa3YP3dOQkJCh7bATpwILteu+2ADU+ioFVZERARyc3MBSEwcTHy8w8VEiMMOG0B+/rd88MG3pKUNp6RkBRdffAwJCVnce++LTpcnEjQKsBLSioqKAIiPT+iwQZwAUlJg1ChfgG06ErECrIiIyK7rYb9+Yx2uJPKceOI4du5czrvvziMjYyzV1Ru5//5LcbtTmDTpetav3+l0iSJtogArIa20tBSA2NhEvN6OC7Dguw927lzo23cgqampJCcnK8CKiIgA8+cvBGDMmFHOFhLBTjrpCLZuXcCqVZsYOfJ0vN5q5sx5koEDu5Oaegi33PIvams9TpcpcsAUYCWklZSUAOB2JwMdG2AnTYKaGli0yMURRxyBy+VSgBUREQFWrvTNATt16giHK5EhQ/qxZMlM6uqqufnmR0hOHkBp6QoeffQ64uJiSU8fwy23PEFVVb3TpYq0igKshLSysjIA3O4koGMD7NFH+37OmQPjx4+nrKyM3NzcxlZhERGRSOT1eikqKgBimTq1t9PliF90dBR/+9svKStbx9q1mznhhGuIje3Ojh2LePTRG0hMjCUxcRDHHns9H3641OlyRVqkACshLRBgwdcC21GjEAP06AFDhuwKsNZaABYvXtxxRYiIiHQy69atw+v14Hb3o08f43Q50ozs7L589NFT1NRsY/XqTZx00o0kJfWnqmo9n3/+JD/60aG4XPGkpR3KCSf8gmnTvsbrtU6XLQIowEqIC0yjY0wK0LEtsOC7D/bLL2Hs2CMa13377bcdW4SIiEgnMn/+fAAyMsY5XIm0xuDB/Xj33ccoL99AXV0dDz74IoMH/4ioqASKi5fyySf/4IILjsLtjiIurjcDBkzlnHPu5fXXv6Ghwet0+RKBFGAlpAUCrLXOBNjjjoOSEsjLS2fAgAEkJCQwb968ji1CRESkE5kx4wMApkw51eFK5EBFR0dx550Xs3r1+9TV7aSiooqHHnqRsWPPJSmpP7W1O9iw4WNmzLifc8+dSHR0FG53F1JThzNs2Bmcc85vePzxdygoKNv/yUQOUpTTBYi0RUVFBQAeTxeg4wPslCm+n5984utGPHPmTAVYERGJaHPmfAnAT396jMOVSFslJsZzxx0Xc8cdFzeuy83dzH/+8z8++eQLcnOXUFa2mdLSVZSWrmTVqlnMmAE33gjgxu1OIi6uGykpPcnI6M+gQdkceugwjjhiBBMnDiUlJdax301ClwKshLTKykoA6utTcbkgObljz5+RAYcc4guwJ500nldffZXNmzezdetWevXq1bHFiIiIOKy8vJzCwvW4XGkcf3x/p8uRdpCd3Y8//vF64PrGddZaFi9eyxtvfMHcuQtYu3YFRUVbqK7eSWXlRior17F161wWLYLXX296NAPE4HbHER2dRHx8FxIT00hJ6UZaWnfS03vSo0d3evfuQb9+GQwcmMHQoX3p0SMFY3R/daRSgJWQFgiwDQ2ppKSAE/+WHX88PP003HPP+MZ18+bN48wzz+z4YkRERBz05JPTAMvgwSc4ck0WZxhjGD16MKNHDwau3uv12to6vvpqBXPmfM/3369g8+aNFBbmU1a2g5qaUurqKqitLaSmJp/i4tYOFmXw3Q3pxpgoXK4oXK5oXK4Y3O4YoqJiiY6OJTo6jujoOGJifMuxsXHExMQRHR1NbGwssbFxxMXF+h++5cTEOOLj40hIiCc2Noa4uFiio6OJjw9sF0NMTDQxMW5iY6P9z6OIi4smNjaa2FjfclxcNFFRboXtIFOAlZBWVVUFQF1dV7p0caaG44+Hv/8damtHExUVhdfr5ZtvvlGAFRGRiOL1ernvvt8C8Le//dbhaqQziY2NYcqUUUyZMmq/29bW1rN69RZWrtzI+vV55OVtp6CgkJ07iygpKaKsrITKyjLq6qqoq6uioaEGj6cWr7cBj6eGhoZKrPUAHsD6H53dngF3X4G3uddas7/Z40ultp6TFoP52rU/MHBgxj6O1zYKsBLSAgG2piaFlBRnapg8GVwu+PLLOMaNG8fSpUt1H6yIiISFwsIyvvsul+XL15Obu5HNm/PYtm0rRUXbKSvbSVVVCfX15TQ0VGFtLeClT5+j+NGPhjtduoSo2NhoDj10AIceOiBox2xo8FBRUc3OnWUUFZVTXFxOeXkNFRVVVFTUUFlZTWVlDVVV1VRV1VBdXU1NTR21tTXU1dXh8TRQX19PfX09DQ31NDQ0YK0Xj8dDQ4MHr9f38Hg8eL1e/yPw3IO1vnWBn+DFWtvk4QV8P33TMtrG6RkDy76nvnWBbXYtN79u38+bBvvA+t2fN7XrGHu+vve2UVHtGzEVYCWkVVdXA1BREedYgO3SBQ4/3Hcf7KRJk5g/fz4LFizA4/HgdrudKUpERGQfKipqeOed+Xz99fesXr2WTZvWs3NnPuXlhdTVleDxVAMNrThSFC5XLNHRScTG9mbQoNF89dWL7V2+yAGJinKTmppEamoSgwY5XY20lQKshLSamhoAyssN6enO1XH88fDQQ3Dzzcfg9T5ERUUFK1asYOTIkc4VJSIiEc3j8fLWW3OZMeMjli79nq1bf6CsbCv19aVAfQt7ReN2J5CQ0IfExDSSk9NIS+tBz54Z9O/fl+zsTA45ZABjxmSTnu7QN8ciEtEUYCWk1dTU4HK5KCvD0W/Ujj8eHngArD2qcd3XX3+tACsiIh2isLCMv/3tNd577wPWr19GWdkWvN7yvbYzJp64uB6kpfUjM3MwQ4cOZfTo4UyePIpDDumH2+1yoHoRkdZTgJWQVldX1xhgnRrECeDIIyEuDr75pisjR45k9erVzJkzh5/97GfOFSUiImFr9uxl/PnP/2HBgi/YsWMNXm9Fk1cNbncyaWkjGDhwJJMnH8M550zh8MNzFFBFJOQpwEpICwTY0lIcuwcWfOH12GPhvffghBMmsWLFCj7//HOstRo6XURE2mz9+u3ce++TfPjhLAoLl2NtTeNrbncKPXuOY/z4Y7nqqnM45ZTDFVRFJGwpwEpIq6urw+2Oorra2QALcNJJcPPNcMMNk/B4HmfLli1s3LiRrKwsZwsTEZGQtHDhD9x664PMm/c/amu3Na6PiurKgAGTufDCC7n11nPp0iXBwSpFRDqWAqyEtPr6+sahup0OsCef7Auw5eXHNK774osvFGBFRKTVNm3awZVX3svnn0+joaHIv9ZN9+6Hcdpp53P//T+jT580R2sUEXGS+pdISGtoaMDt9gVYJ++BBcjO9j2+/roXgwYNIjo6mi+++MLZokREJCT89rfP0bXrMDIz0/nkkydoaCijd+8J3HvvszQ01FFYuJhnnrlb4VVEIp4CrIQ0X4CNBpxvgQVfK+ynn8Lkycfj9XqZM2eO0yWJiEgnVVFRw6mn3orbnczvf385JSWrSEoayPXX/4Xa2mry8ubyu99dpvtZRUSa0L+IEtI8Hg9udwzQeQJsTQ1kZEzF4/GwZs0atm3btv8dRUQkYhQWljFmzHkkJyfxzjt/xeut4bDDzmbFis2Ul//AP/95KzExustLRKQ5bQqwxphpxpjF/scGY8xi//osY0x1k9eebGH/NGPMR8aYtf6fXdtSj0Qer9fbqQLs5MkQHw/bt09pXPf55587WJGIiHQWpaVVHHnkZfTokcaiRa/hciVw1ll3UVVVyeLFrzNsWF+nSxQR6fTaFGCttedba0dZa0cBM4A3mrz8Q+A1a+11LRziLuATa20O8In/uUireb1eXK5YoHME2Lg4mDIFPv00jbFjx+J2u/n444+dLktEIkhbv1z2b/sLY8wqY8xyY8xDHVZ8GLv44t+RmprK3LnP43LFccUVD1BXV8KMGX8kPj7G6fJEREJGUPqnGN9El+cBU/a37R7OAI71Lz8HzAbuDEZNEhmstY0B1ulBnAJOOw3eeQeuuupEFi5cyAcffKD5YEWkw1hrzw8sG2P+ApQ2efkH/5fOLTLGHIfv+nyYtbbWGNOjXQqNEM8//wnXXHMRdXXbgRjOOeceXn31ft3XKiJykIL1r+cxQIG1dm2TdQOMMYuMMZ8bY45pYb+e1tqt/uVtQM8g1SMRoKGhAQBj4oDO0QILcMYZYAx4PFMB2Lx5M2vXrt3PXiIiwdXky+VXDnDX64EHrbW1ANba7cGuLRJUV9cxZMhJXHbZCdTVbWfo0FPYuXMnr732/xReRUTaYL//ghpjPjbGLGvmcUaTzS5k9wvkVqC/tXY0cCvwsjFmn/HCWmsBu486rjXGLDDGLCgsLNxf2RIBqqurAXC54nG5IKGTzOOekQFHHgkLFx5JXJwvXH/44YcOVyUiEehgv1weDBxjjJnn3+7w9i81vDz99HskJXVjzZr3iYvrzYcfLmLlyv+RlpbkdGkiIiFvvwHWWnuCtXZEM4+ZAMaYKOAsYFqTfWqttTv9ywuBH/BdEPdUYIzp5T9OL6DFb3mttU9Za8dZa8elp6cfyO8oYaqsrAwAa+NJSfG1enYWP/kJLF0ay/jxxxIdHc1HH33kdEkiEkba+cvlKCANmADcAUw3LdwDoS+Xd+fxeDnyyMu49tqT8XorOfHEG6mo2MzUqaOcLk1EJGwEow/LCcAqa+2WwApjTLoxxu1fHgjkAOua2XcWcJl/+TJgZhDqkQhRUlLiX0roNPe/BvzkJ76faWmnUF9fzyeffEJ9fb2zRYlI2GjnL5e3AG9Yn/mAF+jeQh36ctlv/frtpKbmMHfu80RFpfHhh9/xwQePqbuwiEiQBeNf1QvY+/6aScAS/8iHrwPXWWuLAIwx/zbGjPNv9yAw1RizFl8QfjAI9UiEKC31jUvi9SZ0mvtfAwYOhMMOg02bTgOgsrKSefPmOVyViESQtny5/BZwnH+7wUAMsKO9Cw5lL788m0GD+lNRsY7MzMmUlW1Vq6uISDtpc4C11l5urX1yj3UzrLWH+KfQGWOtfbvJa1dbaxf4l3daa4+31ub4v00uams9EjkCLbBeb1KnC7AAZ50F332XyZAhh2CM4d1333W6JBGJHG35cvkZYKAxZhnwKnCZf5wKacbdd/+biy+egrV1XHHFA2zYMFvT4oiItCP1a5GQFWiB9Xg6b4C1FgYMOBNrLW+++abTJYlIhGjjl8t11tpL/F2Sx1hrP+3o+kPFGWfcwYMPXgO4eeqpd3jmmbudLklEJOwpwErICgziVF/fOQPsiBG+R36+rxvxqlWr+OGHHxyuSkREguHoo69g1qyHcbmS+PLLJVxzzUlOlyQiEhEUYCVklZeXA1Bbm9TpBnEKuOgiWLLkcLp27QbA22+/vZ89RESksxs37kK++upZoqO7sX79eo46apjTJYmIRAwFWAlZFRUVANTWpnTKFliACy8EcJGdfRYul0vdiEVEQtxhh53NwoWvEhPTk02bcunfv9kBmkVEpJ0owErICgTY+vrOG2CzsuCoo2D79nPxer18+eWXFBVprDIRkVA0efI1LFnyBrGxvdi6NZeMjFSnSxIRiTgKsBKyAgEWunTaAAu+bsQbNx5HcnIqXq+Xd955x+mSRETkAJ177q/54ot/ExWVxubNq0hLS3K6JBGRiKQAKyGrsrLSv5TSae+BBTjvPIiKimLAgPMxxvDKK3vObCEiIp3Z7bc/weuvP4DLlciqVctJT+/E35qKiIQ5BVgJWVVVVf6lrp26BbZ7dzj5ZNi8+XystXz44YfqRiwiEiJmzPiKv/zl50A0X321kEGDMpwuSUQkoinASsjaFWBTO3WABbj6aigunkRSUjc8Ho8GcxIRCQHr12/n3HNPACxPPTWTCROGOF2SiEjEU4CVkFVdXe1f6pzzwDZ10knQu7ebtLSLMMbwwgsvOF2SiIjsg8fjZfjwMVhbw2WX/UHzvIqIdBIKsBKyampq/EtxnT7ARkXBlVfCpk1XYK3liy++oKCgwOmyRESkBaNGnU1NTR7Dhp3Ks8/+2ulyRETETwFWQtauAOvq1IM4BVx1FRgzmu7dh2KtZdq0aU6XJCIizbj//hdZtuwtYmN7sXTpTKfLERGRJhRgJWTV1tZijAHo9C2w4JsTdupUqK+/HoDHH38ca62zRYmIyG5WrtzCvfdeAbj56qs5uN36U0lEpDPRv8oSsmprawEXLhckJDhdTetcfz2Ull6My+VmzZo1LFy40OmSRESkiSOOmAw0cNttjzF27CCnyxERkT0owErIqqurwxg3KSngb4jt9E47DQYO7EZKyqkA/Otf/3K4IhERCbjggnupqFhH//6TePjh650uR0REmqEAKyErEGBD4f7XALcbbr4ZSkp+AcBLL71EZWWlw1WJiMiSJRuZNu3/YUwcixa943Q5IiLSAgVYCVn19fUYExUS9782dcUVkJw8hYSETKqrq5k+fbrTJYmIRLxjjpkKePn97/9NWlqS0+WIiEgLFGAlZDU0NAChF2CTk+GaawxVVXcB8OCDD2owJxERB1133V8oK1tLnz5H8n//d7HT5YiIyD4owErI8gXY6JALsAA33QRu96VERcWzZs0a5syZ43RJIiIRqbS0in/96x4gigUL1HVYRKSzU4CVkOXxeICYkAywmZlw+eWJeL3XAvDAAw84XJGISGQ66qgLgTrOO+8uMjJSnS5HRET2QwFWQpbX68XamJAaxKmpe+4BuA1w8cEHH7Bu3TqHKxIRiSwffbSY5ctnEROTzssv/87pckREpBUUYCVkeb1evN7YkGyBBRg4EH760364XBcAcP/99ztckYhIZDn77HMB+Ne/XsLt1p9EIiKhQP9aS8iy1mJtXMgGWPC1wlp7P2B44YUXyM/Pd7okEZGI8MgjMygvz6Vnz3FcfvlUp8sREZFWUoCVkOS7/xUgtANsTg5cccVAjDkDj8fDfffd53RJIiIR4e67fw7AO++86nAlIiJyIBRgJSRVV1f7l+JD9h7YgN//HmJiHgQMzzzzDFu3bnW6JBGRsHbrrf+krq6A7OwTGTt2kNPliIjIAVCAlZBUVlbmX4oP6RZYgD594I47hgDn4fF4uOuuu5wuSUQkbHk8Xv7+97sAFx9++ILT5YiIyAFSgJWQVFxc7F9KCPkAC/CrX0G3bn8F3Lz44ovk5uY6XZKISFi66qo/4vGUM3bseQwY0MPpckRE5AApwEpIKikp8S+FR4BNToYHH+wF3IzX6+Waa65xuiQRkbD04osPAW4+/PBpp0sREZGDoAArIam0tNS/lBQWARbgyith/Pj7gGRmz57NZ5995nRJIiJh5ZZb/oHHU8ahh55BWlqS0+WIiMhBUICVkLTrHtikkB/EKcDlgqefTsaYJwC49NJLaWhocLgqEZHw8c9/3gcY3n77X06XIiIiB0kBVkJSOLbAAowcCbfffhEwiry8PH73u985XZKISFi4//4Xqa/fSXb2VPr37+50OSIicpAUYCUklZeX+5eSSUhwtJSgu+8+Q1bWa4CLBx54gHXr1jldkohIyPvjH+8B4M03de+riEgoU4CVkBQIsPHxKRjjcDFBlpAAr72WjTG/wev1ctppp2GtdbosEZGQNWvWPKqrN9Oz51hGjOjvdDkiItIGCrASkioqKgBITAyj/sNNjBsHv/3tvcBgVqxYwX333ed0SSIiIev6628B4LHH/uZoHSIi0nYKsBKSAgE2KSlMRnBqxq9/7WLMmI+AKO6//34WLFjgdEkiIiEnL6+I/Px5xMb24txzj3a6HBERaSMFWAlJlZWVAKSkhG+AjYqCWbP6k5T0FNZapk49sfH3FhGR1jn33F8Clquv/pXTpYiISBAowEpIqqqqAiA1tavDlbSvPn3g7bevAM6jpKSYY489TvfDioi0ksfj5ZtvpmNMHI8+epPT5YiISBAowEpICgTYtLTwDrAAxx4Lf/7zS8AgFiz4lp/97GdOlyQiEhLuueffWFvD+PHn4nbrTx4RkXCgf80lJFVXVwPQrVuyw5V0jNtui+KKK+YCyTz99NM89thjTpckItLpPfnkwwA899yDDlciIiLBogArISkQYFNTox2upGMYA08/nc6UKfOAaG666SZee+01p8sSEem01qzJp6xsLSkpgxk8uLfT5YiISJAowEpIqq6uAaBLlzCbBHYf3G743/+GMXLk+4CL88+/gPfee8/pskREOqXLL78bgOuuu83hSkREJJgUYCUk1dTUAi5SwnMa2BbFx8NXX01h2LAZWGs55ZRTeeedd5wuS0Sk05k37w2MieWBB652uhQREQkiBVgJSTU1dYCLLuE7i06LkpNh3rwzGTJkOtZaTjvtNF5++WWnyxIR6TQee2wmXm8FhxzyYw3eJCISZvSvuoSk2tpawB1xLbABycmwYME5jBkzC2tdXHzxxfzxjxqkREQE4I9/fACAJ5/8o8OViIhIsCnASkiqra0DoiI2wAIkJcE335zKaad9DcRyzz13c/bZ52ueWBGJaHV1DWzdupCYmB4cddQwp8sREZEgU4CVkNTQUE+kB1iA6GiYOfMI7rgjF0jnjTemM2DAYHbs2OF0aSIijvjVr54EPEyadI7TpYiISDtQgJWQ1NDQAERHfIAF3xQ7Dz3Ul//9L4+oqGPYuDGXXr36MmvWLKdLExHpcM8//yQATz75G4crERGR9qAAKyHJ42kAYiNyEKeWnHJKNLm5X9C37x9oaKjjjDPO4PTTz8Hj8ThdmohIhygtraK4eAUJCf0ZNCjD6XJERKQdKMBKSPJ6PUAcyclOV9K5ZGbC+vW/5sYbvwfSefvtGXTt2oM5c+Y4XZqISLu74YY/A5ZTT/2p06WIiEg7UYCVkGStF4gnMdHpSjqfqCh47LGRzJ+/lZSUCykvL2LSpEkcf/yplJaWOl2eiEi7mTnzOcDwxBN3Ol2KiIi0EwVYCTm+UXYtUVHxGON0NZ3X4Ye7KSh4mauv/hrowaefvkPXrt255pqf++8hFhEJH0VFFVRWriclJZu0tCSnyxERkXaiACshp66uDoDo6ASHK+n84uLg6acnsmHDNkaOfABro/n3v/9JXFwSV111PVVVVU6XKCISFLff/hgAP/7x+Q5XIiIi7alNAdYYM80Ys9j/2GCMWexfn2WMqW7y2pMt7H+fMSavyXYnt6UeiQxFRUUAxMbqG/bWysw0LFlyN998U0pm5s/xeOCZZ54kKSmFqVPPIDc31+kSRSRIgnBtHmWM+ca/zQJjzBEd+gscpLfeehGARx+9zeFKRESkPbUpwFprz7fWjrLWjgJmAG80efmHwGvW2uv2cZi/Ntnu3bbUI5Fh586dAMTHK8AeqPHjo9mw4R/MnVvBsGH3Ym0iH388i5ycHHr2zOLhh/9KfX2902WKSBsE4dr8EPA7//73+p93anV1DRQXryIurg8ZGalOlyMiIu0oKF2IjTEGOA94JRjHE9mXQIBNSNAQxAdrwoQoVqz4HRs2lHDGGdNxu0eyfftG7rjjVmJj4xk6dDRPPvkvhVmRENaGa7MFArNsdwHyg1lXe/jNb54BvBx99BlOlyIiIu0sWPfAHgMUWGvXNlk3wBizyBjzuTHmmH3se6MxZokx5hljTNeWNjLGXOvvyrSgsLAwSGVLKAoE2KSklP1sKfuTmWl4661zqaxcwn/+U0hW1vVY253Vqxdz/fXXERsbR2bmUG6//Q42bNjgdLkicmAO9tp8C/BnY8xm4GHg7pZO0Fmuzc8//28A/vpXjT4sIhLu9htgjTEfG2OWNfNo+jXnhez+De9WoL+1djRwK/CyMaa5tPEEMAgY5d/nLy3VYa19ylo7zlo7Lj09ff+/mYStnTuLAUhN7eJwJeEjNhauvLI769f/k+3bt/GHP2ykb9+fYW1PNm1azV/+8jADBgwgLi6FsWOP5ve/v59169Y5XbZIxGrna/P1wC+ttf2AXwL/aamOznBt9ni8bNu2mOjobowY0d+RGkREpOPsN8Baa0+w1o5o5jETwBgTBZwFTGuyT621dqd/eSHwAzC4mWMXWGs91jep59NASAwUIc7autUXYNPSUp0tJEylp8Ovf92fzZufpLQ0n2efLefwwx8jJuYoams9fPfdV/z2t/cyaNAgoqJ8LbTnn38xTz/9NFu3bnW6fJGI0J7XZuAydt03+xqd/Nr8zDMfAPUceuhUp0sREZEOEBWEY5wArLLWbgmsMMakA0XWWo8xZiCQA+zVXGOM6WWtDfzF+xNgWRDqkTBXUFACQPfuLfY4lyBJSYHLLkvisstuxNobyc2FGTPyePXVV1i16j1qa5ewadNqNm1azfTpLwPgckXTtWsPBg0axOjRwzniiMOZPHkyAwcOxGjiXpGOctDXZnz3vE4GZgNTgLXNbNNp/P3vTwFw9903OlyJiIh0hGAE2AvYe4CIScDvjTH1gBe4zlpbBGCM+TfwpLV2AfCQMWYUvgEjNgA/C0I9Eua2bdsBQL9+vR2uJLIYAzk5cNddfbjrrtuB29m2DebOrWfmzG+YM+c9tmyZS13danbuLGDnzjzmz/+Cf/0rMFOHITY2kS5d0sjIyGDQoH4MHpzN8OHDGT16NIMHDyY2NtbJX1EknLTl2nwN8Ki/FbcGuLbjyj5wq1fPAWI4++yjnC5FREQ6gLHWOl3DARs3bpxdsGCB02WIQ4466hy+/noG06at4rzzhjhdjuyhqAhWrIB583Yye/bXLF06l+3bv6e6eh2wDSgHPM3ua4wLtzuGuLgEkpKS6dq1Kz16dKdPnwx6986gZ8+e9OzZk4yMDHr37k1GRgZdunQhKioY38VJJDPGLLTWjnO6jlDmxLV506YdZGam06PHWAoK9HeBiEg4aenarL/6JOQUFRUBMHJkP4crkeakpcHRR8PRR3fjtttOA04DoKEB8vJg/XpYvbqW775by8qVi9myZQVFRblUVm6moWEHDQ1lVFRUUlFRxrZtG1m5sjVnNbhcbtzuaKKifI/o6BhiYmKJi/M94uPjSUiIJykpgeTkRJKSEomPjyc+Pp64uLi9HoHX4uPjSUxMJC4ujoSEBBITffvFxMTgdrtxu924XK7G5cBzEWl/d931DwBOP/08hysREZGOogArnYrH46G6uprq6mrKysooKyujtLSU0tJSqqqqqKqqYt0637fsAwcmOFytHIioKMjM9D2OPTYWGOF/7FJfDwUFsHUrbN8ORUWWzZuLWLculy1bNlBYuI3S0kLKynZQVVVEfX0JDQ2leL3leL2VeL3V1NfX4ev1WIGvpdfrfzjJd+/vrnuAjX/ZNK7b8+eufXYt+7Znj2WzxzYt7bP7dk3vR969rj3P3fw+Lf2erbnNuaXj7Pl7tHSOphYunE///r32f1IJSx984Btr6oEHrnO4EhER6SgKsNImHo+HoqIi1q1bx4YNG9i0aROFhYUUFxdTXFxCcXEZpaXllJeXU1lZRU1NNbW1tdTX1+HxNGCtB6/Xi28g6gPpzh6HbpcMP9HR0Lev7+FjgG7+x/gW9/N6oaoKKip2PcrLd/2sqYGqKg8VFdWUlpZTUlJMWVkxFRXlVFdXUlVVTW1tjf/zWUNdXQ0NDTU0NNTi8dTg8dT6H3V4vXX+n/VY2+D//HqwNvDT4/88e9k9QHubrN992few/kfTsB1Yt6/n7Gfdnsfbc/uW9m/p+Z6CeRvKgR+rpKSO/po5JSJ5PF6KilYSE9OT9HTNCy4iEikUYGUv1loKCgpYsmQJK1euZO3ataxbt4FNm/IpLNxBZWU5tbU1eDx1/j/CW8vtf8T4H4lNlmOAWFyuWNzuOFyueKKifI/o6ESioxOIiUkgJiaOmJh4Tjvt+KD/3hK6XC5ISvI9WuYGkvyPjmmxs9b38Hp3Pfb1/GBe83h2P9+ey82t29/rB7NPex1zf4YNa/22El6effYjoIGRI49zuhQREelACrARqr6+nm+//Za5c+eyePH3LF++hi1b8igt3UldXfU+9owC4oAuQBLGpBIT0424uB4kJ2fQtWsfunfvRffuPUhP70mPHj1IT08hNTWG5GRITvaFjMDPhASIi4PYWF8IEQknxvge+myLBN/jj/8HgFtv1QQGIiKRRAE2AqxZs4ZZs2YxZ87XfP/9CrZt20JtbWUzW0YDKcBAYmP7kpLSnx49csjMHMawYSMYOrQfffu6yciA7t2hWzeIj+/gX0ZERARYufIrIIqLLjrW6VJERKQDKcCGmfLycl5//XVmzXqHb75ZyPbtW/B6G5psYfCF1JGkpIwgM3MMI0aMZ8KEsQwblsCAAdC/P8TEOPQLiIiI7EddXQM1NVtJScl2uhQREelgCrAhrq6ujjfeeINnn32Jr7/+hvLyHU1edQE9iY09jEGDjuWYY37McceNZORIFzk5vgFzREREQs1f/vIaYBk3TuMhiIhEGgXYELRz507+8pdHeOGFaWzZso5dI3fGYMxI+vQ5gSlTzuPMM49g/HgXvXs7Wa2IiEhwvfTSqwDcdttVDlciIiIdTQE2RJSWlvL739/Ps8++QlFRfpNXMunR40eccsrVXHLJOI480hAX51iZIiIi7W7t2m+AGE4+eZzTpYiISAdTgO3ErLW88MIL3H//n8nNXY6vpdWFMYcwfPgl/OIXP+ecc5Lp1s3pSkVERDpGaWkVdXXb6dr1EKdLERERByjAdkIVFRXceuuveO6556irq/Kv7c/Qoddw7723c+aZcRr9V0REItKDD74MwJFHnuhwJSIi4gQF2E4kLy+Piy++ki+++BhrvUA0Xbqcxy23PMiNNw6ge3enKxQREXHW669PB+Duu69xuBIREXGCAmwnsHXrVi688Ao+//xDfN2EuzFs2M089thdTJkSjTFOVygiItI5bNz4HcbEctRRw5wuRUREHKAA66CamhrOP/8yZs3yTQcA3Zk8+SGeffYKsrIcLk5ERKSTqa6uo75+J127Dne6FBERcYjL6QIikbWWBx98mOTkVGbNmg50YdKkp9m4sZDZsxVeRUREmvPYY28BMHr0JGcLERERx6gFtoOtXr2aSZN+xPbtG4EoBg26g//970GGDtV3CSIiIvsybdobAFx77YUOVyIiIk5Rauog1lquv/5Whg4dxvbtG0lImMyMGQXk5j6k8CoiItIKq1Z9A7g455yjnS5FREQcohbYDpCXl8fo0UdRWLgRSODcc5/jxRfPISbG6cpERERCR1XVFuLje+N264tfEZFIpStAO3vxxen07z+AwsKNJCVNYd68QqZPV3gVERE5ELNmzQM8ZGcf7nQpIiLiIAXYdmKt5eyzr+DSS8/H6/Vw7LGPsnPnJxxxRILTpYmIiIScf/7zJQDOPvtMZwsRERFHqQtxO6ivr2fEiCNZs2YBxqTxj398wQ03HOJ0WSIiIiFrwYLZANxyyznOFiIiIo5SgA2ygoICBg8eRVnZNmJjR7B48bcMHRrndFkiIiIhragol6iornTpop5MIiKRTF2Igyg3dx39+g2irGwbGRlnsX37EoVXERGRNlq/fjvWVtOrl3oziYhEOgXYIFmyZDlDhx5CfX0lY8feyZYtM0hJMU6XJSIiEvIee+x1ACZOPNbZQkRExHEKsEEwd+5CRo0ag8dTw/HHP8C33z6I2+10VSIiIuHho48+AuDqq892uBIREXGaAmwbLV26iqOOOhJr6zj77L/z8cd3Y9TwKiIiEjTr1i0CXEydOsrpUkRExGEKsG2wYcNmRo8eg7V1nHvuo7z++i+cLklERCTsVFXlExvbw+kyRESkE1CAPUg7dhQxePAIPJ5qpk69j+nTb3K6JBERkbCzbNkmoJ4+fTSAk4iIKMAeFI/Hw8CBh1FfX8YRR/yCDz/8rdMliYiIhKV//tM3gNPRRx/rbCEiItIpKMAehGHDjqO8fAv9+p3MN9/83elyREREwtZnn30CwPXXn+twJSIi0hkowB6gM874OWvXziEhIZu1a9/WgE0iIiLtaOPGpUAUEyYMcboUERHpBBRgD8Df/vYKs2b9E5erC6tXf0dsrN4+ERGR9lRdvZW4uJ5OlyEiIp2EElgrrV69kV/+8qeAi/fe+5q+fZOdLklERCSszZ+/Fmigf/+RTpciIiKdRJTTBYQCr9fLmDFHAg3ccMM/OfHE4U6XJCIiEvaefNI3gNPkyVMcrkRERDoLtcC2wrhx51BVlU929qk8/vj1TpcjIiISEebM+QyAG288z+FKRESks1CA3Y9//OMNFi16k5iYnixf/pbT5YiIiESMzZuXAdEcemim06WIiEgnoQC7DyUlFdx00yWA4YMPPiUmxu10SSIiIhGjtnY7CQm9nS5DREQ6EQXYfTj00B9hbTU/+cndHHus7nsVERHpKJ9+ugTwMGDAYU6XIiIinYgCbAv+7//+w+bNX5OcnM0bb/w/p8sRERGJKP/+9xsAHH/8CQ5XIiIinYkCbDNKSip54IGfA26+/vpTp8sRERGJOHPnfg7ATTed63AlIiLSmSjANuOII36CtbWcd949jBjRz+lyREREIk5+/gqMiWXQoAynSxERkU5EAXYPL774GWvXfkRcXG+mTfu90+WIiIhEHI/HS13dDhIT+zhdioiIdDIKsE14PF6uvNLXVWnGjJkOVyMiIhKZ3ntvIeBl0KDRTpciIiKdjAJsE6effhv19Ts57LCzOfnkcU6XIyIiEpH++1/fAE4nnfRjhysREZHORgHWb9u2Et599zGMiWPu3JedLkdERCRizZ8/B4AbbzzL4UpERKSziXK6gM7imGPOBzxce+2fiI+PcbocERGRiFVQsApj4ujTJ83pUkREpJNRCyzw6adLyc39kNjYDJ588janyxEREYlYHo+X+voikpP7O12KiIh0QgqwwE9+4hu46cknn3O4EhERkcg2bdoXgGXw4LFOlyIiIp1QxAfYxx9/m7Ky1aSnj+byy090uhwREZGI9vzzvgGczjrrNIcrERGRzijiA+wdd9wAwKxZrzhciYiIiCxa9BUAN9xwhsOViIhIZ9SmAGuMmWaMWex/bDDGLG7y2qHGmLnGmOXGmKXGmLhm9k8zxnxkjFnr/9m1LfUcqD/9aRrV1Vvo1+9oJkwY0pGnFhERaRehfm3esWMtbncyXbokdORpRUQkRLQpwFprz7fWjrLWjgJmAG8AGGOigBeB66y1hwDHAvXNHOIu4BNrbQ7wif95h/ntb28BDP/734sdeVoREZF2E8rX5sLCMrzectLSsjvqlCIiEmKC0oXYGGOA84BAP9wTgSXW2u8BrLU7rbWeZnY9AwiMnPQccGYw6mmN3/zmv9TWbmPQoOM59NDMjjqtiIhIhwjFa/M99zwJwJQpp3bUKUVEJMQE6x7YY4ACa+1a//PBgDXGfGCM+c4Y86sW9utprd3qX94G9GzpBMaYa40xC4wxCwoLC9tccGFhES5XEu++q5GHRUQkLIXctbmurg6XK5mHHrqpzccSEZHwZKy1+97AmI+BjGZe+rW1dqZ/myeAXGvtX/zPbwd+DhwOVOHrgvR/1tpP9jh2ibU2tcnzYmvtfu+1GTdunF2wYMH+NhMREWkVY8xCa+04p+toLV2bRUQk3LV0bY7a347W2hP2c+Ao4Cyg6YRtW4AvrLU7/Nu8C4zBd7FsqsAY08tau9UY0wvYvr96REREIp2uzSIiEqmC0YX4BGCVtXZLk3UfACONMQn+i+hkYEUz+84CLvMvXwbMDEI9IiIikU7XZhERCUvBCLAXsGuACACstcXAI8C3wGLgO2vtOwDGmH8bYwJNwQ8CU40xa/FdbB8MQj0iIiKRTtdmEREJS/u9B7Yz0n02IiISTKF2D2xnpGuziIgEU0vX5mCNQiwiIiIiIiLSrhRgRUREREREJCQowIqIiIiIiEhIUIAVERERERGRkKAAKyIiIiIiIiFBAVZERERERERCggKsiIiIiIiIhAQFWBEREREREQkJCrAiIiIiIiISEoy11ukaDpgxphDYGIRDdQd2BOE44U7vU+vofWodvU+to/dp/4L5HmVaa9ODdKyIpGtzh9P71Dp6n1pH71Pr6H3av3a/NodkgA0WY8wCa+04p+vo7PQ+tY7ep9bR+9Q6ep/2T+9ReNJ/19bR+9Q6ep9aR+9T6+h92r+OeI/UhVhERERERERCggKsiIiIiIiIhIRID7BPOV1AiND71Dp6n1pH71Pr6H3aP71H4Un/XVtH71Pr6H1qHb1PraP3af/a/T2K6HtgRUREREREJHREegusiIiIiIiIhIiIDbDGmB8bY1YbY3KNMXc5XU9nYYzpZ4z5zBizwhiz3Bhzs399mjHmI2PMWv/Prk7X6jRjjNsYs8gY8z//8wHGmHn+z9Q0Y0yM0zU6zRiTaox53Rizyhiz0hgzUZ+lvRljfun//22ZMeYVY0ycPk9gjHnGGLPdGLOsybpmPz/G5+/+92uJMWaMc5XLwdK1uXm6Nreers37p2tz6+ja3LzOcG2OyABrjHEDjwMnAcOBC40xw52tqtNoAG6z1g4HJgA/9783dwGfWGtzgE/8zyPdzcDKJs//BPzVWpsNFANXOVJV5/Io8L61dihwGL73S5+lJowxfYCbgHHW2hGAG7gAfZ4AngV+vMe6lj4/JwE5/se1wBMdVKMEia7N+6Rrc+vp2rx/ujbvh67N+/QsDl+bIzLAAkcAudbaddbaOuBV4AyHa+oUrLVbrbXf+ZfL8f2j1gff+/Ocf7PngDMdKbCTMMb0BU4B/u1/boApwOv+TfQeGdMFmAT8B8BaW2etLUGfpeZEAfHGmCggAdiKPk9Ya78AivZY3dLn5wzgeevzDZBqjOnVIYVKsOja3AJdm1tH1+b907X5gOja3IzOcG2O1ADbB9jc5PkW/zppwhiTBYwG5gE9rbVb/S9tA3o6VVcn8TfgV4DX/7wbUGKtbfA/12cKBgCFwH/93bn+bYxJRJ+l3Vhr84CHgU34Lo6lwEL0eWpJS58f/bse+vTfsBV0bd6nv6Fr8/7o2twKujYfsA69NkdqgJX9MMYkATOAW6y1ZU1fs76hqyN2+GpjzKnAdmvtQqdr6eSigDHAE9ba0UAle3RJivTPEoD/PpEz8P1R0RtIZO+uOdIMfX4k0uja3DJdm1tN1+ZW0LX54HXE5ydSA2we0K/J877+dQIYY6LxXSBfsta+4V9dEGjy9//c7lR9ncBRwOnGmA34urhNwXc/Saq/mwnoMwW+b9m2WGvn+Z+/ju+iqc/S7k4A1ltrC6219cAb+D5j+jw1r6XPj/5dD336b7gPujbvl67NraNrc+vo2nxgOvTaHKkB9lsgxz+SWAy+m7JnOVxTp+C/X+Q/wEpr7SNNXpoFXOZfvgyY2dG1dRbW2ruttX2ttVn4PjufWmsvBj4DzvFvFtHvEYC1dhuw2RgzxL/qeGAF+iztaRMwwRiT4P//L/A+6fPUvJY+P7OAn/pHPJwAlDbpziShQdfmFujavH+6NreOrs2tpmvzgenQa7PxtfJGHmPMyfjulXADz1hr/5+zFXUOxpijgTnAUnbdQ3IPvnttpgP9gY3AedbaPW/gjjjGmGOB2621pxpjBuL71jcNWARcYq2tdbA8xxljRuEbTCMGWAdcge+LM32WmjDG/A44H99Io4uAq/HdIxLRnydjzCvAsUB3oAD4LfAWzXx+/H9g/ANfF68q4Apr7QIHypY20LW5ebo2Hxhdm/dN1+bW0bW5eZ3h2hyxAVZERERERERCS6R2IRYREREREZEQowArIiIiIiIiIUEBVkREREREREKCAqyIiIiIiIiEBAVYERERERERCQkKsCIiIiIiIhISFGBFREREREQkJCjAioiIiIiISEj4/2PxsmOlK5CPAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(16,6))\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot(tv,vsec[0].T,c='b')\n",
    "for n in range(len(vsec)-1):\n",
    "    plt.plot(tv,vsec[n+1].T,c='k')\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(tv,vsoma[0].T,c='b')\n",
    "for n in range(len(vsoma)-1):\n",
    "    plt.plot(tv,vsoma[n+1].T,c='k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t1 \n"
     ]
    }
   ],
   "source": [
    "testKa = 0.07\n",
    "cellID = 7\n",
    "cutExp = 0\n",
    "\n",
    "for sec in h.allsec(): h.delete_section(sec=sec)\n",
    "cell1 = L23(cellID=cellID,cutExperiment=cutExp,dendNa=[naDensity,None,None,False],dendK=[testKa,np.inf,True,False],dxSeg=1,fixDiam=None);\n",
    "\n",
    "syn = None\n",
    "syn = None\n",
    "onset=50\n",
    "tau=1\n",
    "gmax=0.0005\n",
    "tstop = 120\n",
    "epspDendrite,epspSoma,epspTV,syn = mfx.injectAlphaSites(cell1.sectionList,cell1.segmentList,syn=syn,onset=onset,tau=tau,gmax=gmax,tstop=tstop)\n",
    "vEpspDend = np.array(epspDendrite)\n",
    "vEpspSoma = np.array(epspSoma)\n",
    "tvEpsp = np.array(epspTV)\n",
    "cEpspAmpDend = np.amax(vEpspDend,axis=1) - vEpspDend[:,np.where(tvEpsp>=onset-1)[0][0]]\n",
    "cEpspAmpSoma = np.amax(vEpspSoma,axis=1) - vEpspSoma[:,np.where(tvEpsp>=onset-1)[0][0]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fdb91c3ab38>]"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIYAAAFlCAYAAACeB+OLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAACDAUlEQVR4nO3dd3xUVf7/8ddJQhI6JIQO0kFNKCEoqIAKoisKivsDsYCrrqjo6u6qi7vqqvtVV11RcS0L9oaKigU7ihVEAgSIKFUQpIUSCKQn5/fHnYEQJqRN5k55Px+P+8jkTrmfS0hu8p7zOcdYaxERERERERERkcgT5XYBIiIiIiIiIiLiDgVDIiIiIiIiIiIRSsGQiIiIiIiIiEiEUjAkIiIiIiIiIhKhFAyJiIiIiIiIiEQoBUMiIiIiIiIiIhEqxu0CymrRooXt1KmT22WIiASlxYsX77TWJrldh5t0nRAR8U3XCIeuEyIivh3tOhFUwVCnTp1IT093uwwRkaBkjNnodg1u03VCRMQ3XSMcuk6IiPh2tOuEWslERERERERERCKUgiERERERERERkQilYEhERERERIKSMaanMSajzLbPGHOj23WJiISToJpjSERERERExMtauwroC2CMiQZ+A2a7WZOISLjRiCEREREREQkFw4B11lpNtC0i4kcKhkREREREJBRcCMwsv9MYc5UxJt0Yk56VleVCWSIioU3BkIiIiIiIBDVjTCwwCphV/j5r7XRrbZq1Ni0pKSnwxYmIhDgFQyIiIiIiEux+Byyx1m53uxARkXCjYEhERERERILdeHy0kYmISO0pGBIRERERkaBljGkInAG87XYtIiLhSMvVi4iIiIhI0LLWHgAS3a5DRCRcacRQkFi7di27du1yuwwREQlx27dvZ8mSJVhr3S5FRERERPyksLCQ119/vU5eW8FQEFi5ciXHHXccQ4YM0S/yIiJSY6tWraJHjx7079+ff/zjH26XIyIiIiJ+sGLFCk488UQuvPBCvv/+e7+/voKhIPD8889TVFTEypUr6+SLLCIikeHmm28mOjqaUaNG8e9//5uMjAy3SxIRERGRGiopKeH+++8nLS2NLVu28M477zBw4EC/H0fBUBBYuHAhPXr0AOCrr75yuRoREQlFW7Zs4YMPPuCaa67hhRdeoGHDhjz88MNulyUiIiIiNbB27VqGDBnClClTOOecc8jMzGT06NF1ciwFQy6z1pKRkcHw4cPp2rUr6enpbpckIiIh6JNPPqG0tJQLL7yQZs2aMX78eN5++23y8vLcLk1EREREqshayxNPPEGfPn1YuXIlL7/8Mm+++SZJSUl1dkwFQy7bs2cP+/bto2vXrvTv358lS5a4XZKIiISguXPn0qpVK5KTkwEYO3Ys+/fv5+OPP3a5MhERERGpis2bN3PmmWcyefJkTjnlFFasWMHFF1+MMaZOj6tgyGUbN24EoFOnThx77LFs2LCBgoICl6sSEZFQ8/XXX3Paaacd/MXh1FNPpUmTJnzyyScuVyYiIiIiR2Ot5aWXXiI5OZnvvvuOJ598ko8//pj27dsH5PgKhly2YcMGAI455hi6deuGtZZffvnF3aJERCSk7Nq1i82bN9O/f/+D+2JiYhg6dChz5851sTIREREROZodO3ZwwQUXMGHCBJKTk1m+fDlXX311nY8SKkvBkMvKBkNdu3YFYN26dS5WJCIioWb58uUA9OnT57D9w4YNY926dQevNSIiIiISPGbPnk1ycjIffPABDzzwAF999dXBXCCQYgJ+RDnM9u3biYmJITExEWst4Mw+LiIiUlXLli0DoHfv3oftHzp0KADz58+nU6dOgS5LRERERHzIzs7mhhtu4MUXX6Rfv3588cUXB+eJdINGDLls165dtGjRAmMMLVq0oHHjxhoxJCIi1bJ8+XJatmxJq1atDtufnJxMgwYNWLhwoUuViYiIiEhZn332GSkpKbzyyivcfvvtfP/9966GQqARQ67buXMnLVq0AMAYQ4cOHdi8ebPLVYmISChZtWoVxx577BH7Y2Ji6N+/v4IhEREREZcdOHCAv/3tbzz++OP06tWLBQsWMGDAALfLAjRiyHVlgyGAdu3a8dtvv7lYkYiIhJr169dX2I9+4oknsnTpUgoLCwNclYiIiIiA09bft29fHn/8cW688UaWLFkSNKEQKBhynYIhERGpjdzcXLZt20aXLl183n/CCSdQWFh4cB4iEREREQmMgoICpkyZwuDBgykuLmbevHk8/PDD1K9f3+3SDqNgyGXeOYa82rdvz7Zt2ygpKXGxKhERCRW//PILQIXBkPfdqCVLlgSsJhEREZFIl5GRwYABA7j//vu5/PLLWb58OaeeeqrbZflUq2DIGPO6MSbDs20wxmSUue9WY8xaY8wqY8yZta40DJWWlh4RDLVr146SkhK2b9/uYmUiIhIq1q9fD1QcDB1zzDE0adLk4JL2IiIiIlJ3iouLueeeezjhhBPIyspizpw5zJgxg8aNG7tdWoVqNfm0tXac97Yx5iFgr+f2ccCFwPFAW2CuMaaHtVbDYMrIzs6mtLT0iGAIYPPmzbRt29at0kREJERUFgwZY+jdu7eCIREREZE6tmrVKiZOnMjChQsZN24cjz/+OImJiW6XVSm/tJIZYwwwFpjp2TUaeM1aW2Ct/QVYC5zgj2OFk127dgEc9h/FGwxpniEREamKX3/9lfr16x/2JkN53mDIWhvAykREREQiQ2lpKdOmTaNfv36sWbOG1157jddeey0kQiHw3xxDg4Ht1to1ns/bAZvK3L/Zs0/K2Lt3LwBNmzY9uM8bDG3ZssWVmkREJLRs2bKFdu3a4bxH41vv3r3Zt28fv/76awArExEREQl/GzduZPjw4dxwww2cdtppZGZmMm7cuMqfGEQqDYaMMXONMZk+ttFlHjaeQ6OFqsUYc5UxJt0Yk56VlVWTlwhZOTk5AIf1GrZo0QJjDDt27HCrLBERCSFbtmyhTZs2R31M7969AdROJiIiIuIn1lqee+45UlJSWLRoETNmzGDOnDmV/l4WjCoNhqy1w621yT62dwGMMTHAGOD1Mk/7DehQ5vP2nn2+Xn+6tTbNWpuWlJRU8zMJQfv27QOgSZMmB/dFR0fTokULBUMiIlIlW7ZsqXROuuTkZEDBkIiIiIg/bNu2jdGjR3P55ZeTmprK8uXLufLKK486gjuY+aOVbDjws7V2c5l97wEXGmPijDGdge7AD344VljxNWIIoGXLlgqGRESkUtZatm7dWmkw1LhxY7p06aJgSERERKSWZs2aRXJyMp9++ilTp07liy++oHPnzm6XVSu1WpXM40LKtZFZa380xrwBrASKgclakexIvkYMgYIhEQkNxpjXgZ6eT5sB2dbavsaYTsBPwCrPfd9ba6/28fwHgXOBQmAd8AdrbXYdlx1WcnJyOHDgQJVWsTzuuOP4+eefA1CViIiISPjZvXs31113HTNnziQtLY0XX3yRY4891u2y/KLWI4astZdZa5/ysf8ea21Xa21Pa+1HtT1OODraiKHt27e7UZKISJVZa8dZa/taa/sCbwFvl7l7nfc+X6GQx2dAsrW2N7AauLVuKw4/3oUKqtLL3qtXL1atWkVJid6nEREREamOjz76iJSUFGbNmsXdd9/N/PnzwyYUAv+tSiY1sG/fPqKjo6lfv/5h+1u1aqURQyISMozTTD2Wai5CYK391Fpb7Pn0e5z56KQatm7dClClEUO9evWioKCAjRs31nVZIiIiImFh//79TJo0ibPPPpvmzZuzcOFCbr/9durVq+d2aX6lYMhFOTk5NG7c+IgJqlq2bMnevXspKChwqTIRkWoZDGy31q4ps6+zMWapMeYrY8zgKrzG5YDP0aWRvHplZbwjhqoSDHnf1VI7mYiIiEjlvvnmG3r37s2MGTO4+eabSU9PJzU11e2y6oSCIRfl5OQcMb8QOMEQgP4AEhG3GWPmGmMyfWyjyzxsPIePFtoKdLTW9gP+ArxqjDnyh92hY/wDZz66V3zdH8mrV1bGO2KoKq1kPXs600EpGBIRERGpWH5+PjfddBNDhw7FGMPXX3/NAw88QHx8vNul1Rl/TD4tNbRv374j5heCQ8HQjh07aN9enRUi4h5r7fCj3W+MiQHGAP3LPKcAKPDcXmyMWQf0ANJ9PP8y4BxgmLXW+q/yyJCVlUV8fLzPa0l5iYmJJCUlKRgSERERqcDixYuZMGECK1eu5Oqrr+bBBx+kUaNGbpdV5zRiyEUVjRhq1aoVgCagFpFQMBz42Vq72bvDGJNkjIn23O4CdAfWl3+iMeYs4BZglLU2N0D1hpWdO3eSmJh4REtyRXr16qVgSERERKScoqIi7rrrLgYOHEh2djYff/wxTz75ZESEQqBgyFVVGTEkIhLkLuTISaeHAMuNMRnAm8DV1trdAMaYp40xaZ7H/RdoDHxmjMkwxhyxwqUc3c6dO2nRokWVH69gSERERORwK1euZNCgQdx5552MGzeOzMxMzjzzTLfLCii1krkoJyfHZ6uYNxjSiCERCXbW2st87HsLZ/l6X4+/ssztbnVXWWSoSTCUlZXFrl27SExMrMPKRERERIJbSUkJjzzyCP/4xz9o3Lgxb775JhdccIHbZblCI4ZclJubS8OGDY/Y37BhQ+Lj49m5c6cLVYmISKioSTAEmoBaREREItv69es5/fTTuemmmzjrrLPIzMyM2FAIFAy5Kjc3lwYNGhyx3xhDYmIiu3btcqEqEREJFTUNhn766ae6KklEREQkaFlrmT59Or179yYjI4Pnn3+e2bNnH5znN1KplcxFFQVDgIIhERE5qpKSEvbs2VOtlrBjjjmG2NhY1qxZU4eViYiIiASfLVu2cOWVV/LRRx8xbNgwnn32WTp27Oh2WUFBI4ZcYq0lLy+P+vXr+7xfwZCIiBzNnj17sNZWa8RQdHQ0Xbp0UTAkIiIiEcNay8yZM0lOTubLL7/kscce49NPP1UoVIaCIZcUFRVRUlJy1BFDmmNIREQq4r1GVCcYAujevTtr166ti5JEREREgsrOnTsZN24cF110ET179iQjI4PrrruOqChFIWXpX8Mlubm5AGolExGRGqlpMNStWzfWrl1LaWlpXZQlIiIiEhTmzJlDcnIy77zzDvfeey/ffPMNPXr0cLusoKRgyCV5eXkAFbaStWjRgt27d2OtDWRZIiISIrzBUHWXne/evTt5eXls3bq1LsoSERERcdW+ffu44oorOPfcc2nVqhWLFi3i1ltvJSZGUyxXRMGQS6oyYqikpIS9e/cGsiwREQkR3lGlNRkxBGieIREREQk78+bNo3fv3jz//PPceuut/PDDD/Tp08ftsoKegiGXVCUYAtROJiIiPtVmjiFA8wyJSEgwxjQzxrxpjPnZGPOTMWaQ2zWJSPDJy8vjxhtv5PTTTyc2NpZvv/2We++9l7i4OLdLCwkKhlziDYaOtioZoAmoRUTEp507dxIfH1/hGwwV6dChA7GxsQqGRCRUPAp8bK3tBfQBfnK5HhEJMgsXLqRfv348+uijXHfddSxdupRBg5QhV4eCIZd45xjSiCEREamJPXv2kJCQUO3nacl6EQkVxpimwBDgGQBrbaG1NtvVokQkaBQWFnLbbbdx0kknkZuby9y5c3nsscdo2LCh26WFHAVDLqmslczbGqBgSEREfMnOzqZZs2Y1eq53ZTIRkSDXGcgCnjPGLDXGPG2MOeIvPmPMVcaYdGNMelZWVuCrFJGAW7FiBSeeeCL33HMPEyZMYMWKFQwbNsztskKWgiGXVLWVTMGQiIj4UptgqHv37qxdu1YrX4pIsIsBUoEnrbX9gAPAlPIPstZOt9amWWvTkpKSAl2jiARQSUkJ999/P2lpaWzZsoV33nmH5557jqZNm7pdWkhTMOSSylrJmjZtSlRUlOYYEhERn7Kzs2v8S1C3bt3Izc3VkvUiEuw2A5uttQs9n7+JExSJSARau3YtQ4YMYcqUKZx77rlkZmYyevRot8sKCwqGXFLZiKGoqCgSEhI0YkhERHzau3dvrUYMgZasF5HgZq3dBmwyxvT07BoGrHSxJBFxgbWWJ554gj59+rBy5UpefvllZs2ahUYI+o+CIZdUNscQOO1kCoZERMSX2s4xBFqyXkRCwvXAK8aY5UBf4F53yxGRQNq0aRNnnnkmkydPZvDgwaxYsYKLL74YY4zbpYWVGLcLiFSVtZKBgiEREfHNWlurVrIOHTpQr149jRgSkaBnrc0A0tyuQ0QCy1rLyy+/zPXXX09RURFPPvkkkyZNUiBURzRiyCW5ublERUVRr169Ch/TokULBUMiInKEvLw8iouLazxiKCYmhi5dumjEkIiIiASdHTt2cMEFFzBhwgSSk5NZvnw5V199tUKhOqRgyCW5ubk0aNDgqP+5ExMTNfm0iIgcITs7G6DGwRA47WQaMSQiIiLBZPbs2SQnJ/PBBx/wwAMP8NVXX9G1a1e3ywp7CoZc4g2GjkatZCIi4os3GKrN0qxasl5ERESCRXZ2NhMnTmTMmDG0b9+exYsXc/PNNxMdHe12aRFBwZBL8vLyKlyRzCsxMZH8/PyDE1WLiIiAsyIZ1H7EkJasFxEREbd99tlnpKSk8Morr3D77bfz/fffk5yc7HZZEUXBkEtyc3OrFAwBGjUkIiKH8UcrmXfJ+nXr1vmhIhEREZHqOXDgAJMnT2bEiBE0atSIBQsWcPfddxMbG+t2aRFHwZBLCgoKiI+PP+pjFAyJiIgv/mgl05L1IiIi4pb58+fTt29fnnzySf785z+zZMkSBgwY4HZZEUvBkEsKCgqIi4s76mMUDImIiC/+aCXr2LEjMTExCoZEREQkYAoKCpgyZQqDBw+muLiYL774gqlTp1baTSN1K8btAiJVVYKhhIQEAHbv3h2IkkREJET4o5UsJiaGzp07KxgSERGRgMjIyGDChAmsWLGCK6+8kqlTp9K4cWO3yxI0Ysg1GjEkIiI1lZ2dTb169SptSa5Mt27dFAyJiIhInSouLuaee+7hhBNOICsrizlz5jBjxgyFQkFEwZBLNGJIRERqKjs7m2bNmmGMqdXreIMhLVkvIiIidWHVqlWcfPLJ3HbbbYwZM4bMzExGjhzpdllSTq2CIWPM68aYDM+2wRiT4dmfaIyZZ4zZb4z5r18qDTNVCYbi4+Np0KCBRgyJiMhh9u7dW6uJp726devGvn372Llzpx+qEhEREXGUlpYybdo0+vbty9q1a3nttdd47bXXDnbFSHCp1RxD1tpx3tvGmIeAvZ5P84HbgWTPJuVUJRgCp51MI4ZERKSsffv20aRJk1q/TtmVyZKSkmr9eiIiIiIbN27kD3/4A/PmzePss8/m6aefpk2bNm6XJUfhl1Yy44xlHwvMBLDWHrDWfosTEIkPVQ2GEhISNGJIREQOk5OT45e+fC1ZLyIiIv5ireW5554jJSWFRYsWMWPGDObMmaNQKAT4a46hwcB2a+0aP71e2NOIIRERqSl/BUOdOnUiKipKwZCIiIjUyrZt2xg9ejSXX345qampLF++nCuvvLLW8yFKYFQaDBlj5hpjMn1so8s8bDye0ULVZYy5yhiTboxJz8rKqslLhKT8/HyNGBIRkRrZv3+/X4Kh2NhYjjnmGAVDIiIiUmOzZs0iOTmZTz/9lKlTp/LFF1/QuXNnt8uSaqh0jiFr7fCj3W+MiQHGAP1rUoC1djowHSAtLS1ilkWpzoghBUMiIlJWTk4OjRo18stracl6ERERqYndu3dz3XXXMXPmTAYMGMALL7zAscce63ZZUgP+aCUbDvxsrd3sh9eKCNZaCgsLq9VKpqWERUTEy1+tZKBgSERERKrvo48+IiUlhVmzZnH33Xczf/58hUIhzB/B0IX4aCMzxmwApgKXGWM2G2OO88OxwkJhYSFAlVvJSkpK2LdvX12XJSIiIaCkpITc3Fy/BkO7d+/WfHYiIiJSqZycHCZNmsTZZ59N8+bNWbhwIbfffjsxMbVa8FxcVuuvnrX2sgr2d6rta4ergoICAOLj4yt9bGJiIuAM02vatGmd1iUiIsFv//79AH4NhgDWrVtHQkKCX15TREREws8333zDxIkT2bBhAzfffDN33313lf6mleDnr1XJpBq8wVBVRwwBmmdIREQA55068H8wpHYyERER8SU/P5+bbrqJoUOHYozh66+/5oEHHlAoFEYUDLmgOsFQ2RFDIiLBxBjzujEmw7NtMMZkePZ3MsbklbnvqUpe56/GGGuMaRGQwkOcv0cMdenSBWMMa9as8cvriYiISPhIT08nNTWVhx56iEmTJrFs2TJOOeUUt8sSP1MjoAs0YkhEwoG1dpz3tjHmIWBvmbvXWWv7VvYaxpgOwAjgV78XGKa8I4b8tSpZfHw87du314ghEREROaioqIh77rmH//u//6N169Z8/PHHnHnmmW6XJXVEwZALNGJIRMKJMcYAY4HTa/D0h4FbgHf9WlQY83crGWhlMhERETlk5cqVTJgwgcWLF3PJJZcwbdo0mjdv7nZZUofUSuaC6gRD3m9AjRgSkSA2GNhurS3bi9TZGLPUGPOVMWawrycZY0YDv1lrlwWkyjChYEhERETqQklJCQ899BCpqals3LiRN998k5deekmhUATQiCEXVCcYqlevHk2aNFEwJCKuMMbMBVr7uOsf1lrvKJ/xwMwy920FOlprdxlj+gPvGGOOt9buK/O6DYC/47SRVVbDVcBVAB07dqzZiYSRugqGsrKy2Lt3r1bAFBERiUDr16/nsssu45tvvmH06NH873//o1WrVm6XJQGiYMgF1QmGwGknUyuZiLjBWjv8aPcbY2KAMUD/Ms8pAAo8txcbY9YBPYD0Mk/tCnQGljmdaLQHlhhjTrDWbitXw3RgOkBaWpqt7TmFuroKhsBZsj41NdVvrysiIiLBzVrLjBkz+Mtf/kJ0dDTPP/88EyZMwPP7mUQItZK5oLrBUEJCgkYMiUiwGg78bK3d7N1hjEkyxkR7bncBugPryz7JWrvCWtvSWtvJWtsJ2Ayklg+F5Ej+XpUMtGS9iIhIJNqyZQsjR45k0qRJDBw4kBUrVjBx4kSFQhFIwZALNGJIRMLIhRzeRgYwBFjuWb7+TeBqa+1uAGPM08aYtMCWGF68I4YaNGjgt9fs2rUroGBIREQkElhrmTlzJsnJyXz55Zc89thjfPrpp2rZj2BqJXNBTUYMrV+/vvIHiogEmLX2Mh/73gLequDxV1awv5NfCwtjOTk5NGrUiKgo/72307BhQ9q0aaNgSEREJMzt3LmTa6+9llmzZjFw4EBeeOEFevTo4XZZ4jKNGHKBRgyJiEhN5eTk+LWNzEsrk4mIiIS3OXPmkJyczDvvvMO9997LN998o1BIAAVDrsjPzweqN2Joz549lJSU1GVZIiISAhQMiYiISHXs27ePK664gnPPPZdWrVqxaNEibr31VmJi1EAkDgVDLqjJiCFrLdnZ2XVYlYiIhIK6DIa2bt3KgQMH/P7aIiIi4o558+aRkpLC888/z6233soPP/xAnz593C5LgoyCIRfUZI4hQO1kIiLC/v376ywYAmfJehEREQltubm53HDDDZx++unExcXx3Xffce+991b5b1CJLAqGXOANhuLj46v0+MTERAAtWS8iIgcnn/Y3LVkvIiISHhYuXEhqairTpk3j+uuvJyMjg4EDB7pdlgQxBUMuqEkrGWjEkIiI1F0rmZasFxERCW2FhYXcdtttnHTSSeTm5jJ37lymTZtGgwYN3C5Ngpxmm3JBQUEBMTExVV5q2NtKphFDIiKyf//+Ohkx1LRpU5KSkhQMiYiIhKAVK1YwYcIEMjIyuOyyy3jkkUdo2rSp22VJiNCIIRcUFBRUq7dTI4ZERMTrwIEDdRIMgVYmExERCTUlJSXcf//9pKWlsWXLFt59912ee+45hUJSLQqGXFDdYKhp06YYYzRiSEQkwllrOXDgQJ0NCVcwJCIiEjrWrl3LkCFDmDJlCueeey6ZmZmMGjXK7bIkBCkYckF1g6Ho6GiaN2+uEUMiIhEuPz8fay0NGzask9fv1q0bmzZtIi8vr05eX0RERGrPWssTTzxBnz59WLlyJS+//DKzZs0iKSnJ7dIkRCkYckFhYSGxsbHVek5CQoJGDImIRLjc3FyAOg2GAH755Zc6eX0RkZowxmwwxqwwxmQYY9LdrkfETZs2beLMM89k8uTJDB48mMzMTC6++GKMMW6XJiFMwZALioqKqFevXrWek5iYqGBIRCTCHThwAKj7YEjtZCIShE6z1va11qa5XYiIG6y1vPTSS6SkpDB//nyefPJJPvroI9q1a+d2aRIGFAy5oLi4mJiY6i0Il5CQoFYyEZEI5w2G6nKOIVAwJCIiEkx27NjBBRdcwIQJE0hJSWHZsmVcffXVGiUkfqNgyAUaMSQiIjVR1yOGEhISaN68uYIhEQk2FvjUGLPYGHNV+TuNMVcZY9KNMelZWVkulCdSd2bPnk1ycjIffPABDz74IF9++SVdu3Z1uywJM9UbtiJ+UdNgSCOGREQiW13PMQRamUxEgtIp1trfjDEtgc+MMT9ba7/23mmtnQ5MB0hLS7NuFSniT9nZ2fzpT3/ipZdeIjU1lXnz5nH88ce7XZaEKY0YckFxcXG1g6GEhAT27dtHUVFRHVUlIiLBrq5HDIGCIREJPtba3zwfdwCzgRPcrUikbn322WekpKTw6quvcscdd/D9998rFJI6pWDIBUVFRdWeYygxMRGAPXv21EVJIiISAgIVDG3cuJHCwsI6O4aISFUZYxoaYxp7bwMjgEx3qxKpGwcOHGDy5MmMGDGCRo0asWDBAu66665qDyoQqS4FQy6oSStZQkICgOYZEhGJYHU9+TQ4wVBpaSkbNmyos2OIiFRDK+BbY8wy4AfgA2vtxy7XJOJ38+fPp2/fvjz55JP8+c9/ZsmSJQwYMMDtsiRCaI4hFxQVFREfH1+t53hHDGmeIRGRyBWoOYbAWZmsR48edXYcEZGqsNauB/q4XYdIXSkoKOCf//wnDz74IB07dmTevHkMHTrU7bIkwigYckFNl6sHjRgSEYlkgWolAy1ZLyIiUtcyMjKYMGECK1as4Morr2Tq1Kk0btzY7bIkAqmVzAU1XZUMFAyJiESyQLSSJSUl0bhxYwVDIiIidaS4uJh77rmHE044gaysLObMmcOMGTMUColrNGLIBbWZY0itZCIikSs3N5e4uDiio6Pr7BjGGK1MJiIiUkdWrVrFhAkT+OGHHxg3bhyPP/74wUEAIm7RiCEX1KSVrEmTJsTExGjEkIhIBDtw4ECdtpF5KRgSERHxr9LSUqZNm0bfvn1Zu3Ytr732Gq+99ppCIQkKCoZcUJMRQ8YYEhISNGJIRCSCBTIY+uWXXyguLq7zY4mIiIS7jRs3Mnz4cG644QaGDRtGZmYm48aNc7sskYNqFQwZY143xmR4tg3GmAzP/jOMMYuNMSs8H0/3S7VhoibBEDjtZBoxJCISuQIVDPXo0YPi4mJ++eWXOj+WiIhIuLLW8uyzz5KSksKiRYt4+umnef/992nTpo3bpYkcplZzDFlrD8acxpiHgL2eT3cC51prtxhjkoFPgHa1OVY4KS4urlEwlJiYqBFDIiIRLDc3NyDBUK9evQD4+eef6d69e50fT0REJNxs27aNP/7xj8yZM4ehQ4fy/PPP06lTJ7fLEvHJL61kxhgDjAVmAlhrl1prt3ju/hGob4yJ88exwkFRUVG15xgCjRgSEYl0Bw4cqNMVybx69uwJOBNkioiISPXMmjWL5ORk5s6dy8MPP8wXX3yhUEiCmr/mGBoMbLfWrvFx3wXAEmttga8nGmOuMsakG2PSs7Ky/FROcKtpK1liYqKCIRGRCBaoVrLmzZvTsmVLfv755zo/loiISLjYvXs3F110EWPHjqVLly4sXbqUG2+8kagoTe0rwa3S/6HGmLnGmEwf2+gyDxuPZ7RQueceD9wPTKro9a210621adbatKSkpJqcQ8ipaSuZJp8WEYlsgQqGwGknUzAkIiJSNR999BHJycnMmjWLu+++m/nz5x9szRYJdpX2M1lrhx/tfmNMDDAG6F9uf3tgNjDBWruuNkWGm5q2kiUmJpKbm0t+fj7x8fF1UJmIiASzQM0xBE472ezZswNyLBERkVCVk5PDTTfdxPTp0zn++OOZM2cOqampbpclUi3+GNM2HPjZWrvZu8MY0wz4AJhirf3OD8cIG9baWrWSARo1JCISoQI1xxA4I4Z27tzJzp07A3I8ERGRUPP111/Tp08fZsyYwc0330x6erpCIQlJ/giGLuTINrLrgG7AHWWWs2/ph2OFvJKSEoAat5IBmmdIRCRCBbqVDDQBtYiISHn5+fn89a9/5dRTT8UYw9dff80DDzygrg4JWbVarh7AWnuZj33/B/xfbV87HBUXFwPUuJUMNGJIRCQSlZaWBryVDJxg6OSTTw7IMUVERIJdeno6EyZM4KeffuKaa67hgQceoFGjRm6XJVIrmh49wIqKigCNGBIRkerJz88HCFgw1KlTJ2JjYzUBtYiICM7fcXfeeScDBw5k3759fPzxxzzxxBMKhSQs1HrEkFRPbYIhjRgSEYlcBw4cAAIXDEVHR9OjRw8FQyIiEvFWrlzJhAkTWLx4MZdccgnTpk2jefPmbpcl4jcaMRRgtWkl04ghEZHI5Q2GAjX5NDjtZJpjSEREIlVJSQkPPfQQqampbNy4kbfeeouXXnpJoZCEHQVDAVabEUMNGjQgLi5OwZCISAQK9IghcCagXrduHYWFhQE7poiISDBYv349p512GjfddBNnnXUWmZmZjBkzxu2yROqEgqEAq00wZIwhISFBrWQiIhEoNzcXCHwwVFJSwrp16wJ2TBERETdZa5k+fTq9e/dm2bJlPP/888yePZtWrVq5XZpInVEwFGDeVrKaBEPgzDOkEUMiEgyMMa8bYzI82wZjTIZnfydjTF6Z+546ymtcb4z52RjzozHmgYAVH4LcaiUDLVkvIiKRYcuWLYwcOZJJkyYxcOBAVqxYwcSJEzHGuF2aSJ3S5NMB5h0xVJM5hsAJhjRiSESCgbV2nPe2MeYhYG+Zu9dZa/se7fnGmNOA0UAfa22BMaZlnRQaJtwYMeQNhjQBtYiIhDNrLa+99hqTJ08mPz+fxx57jGuvvZaoKI2jkMig/+kBVptWMnAmoNaIIREJJsZ5G20sMLOaT70G+Le1tgDAWrvD37WFk7y8PADq168fsGM2adKEtm3bKhgSEZGwtXPnTsaNG8dFF11Ez549ycjI4LrrrlMoJBFF/9sDrLbBkEYMiUgQGgxst9auKbOvszFmqTHmK2PM4Aqe1wMYbIxZ6HncAF8PMsZcZYxJN8akZ2Vl+bv2kOFGMARamUxERMLX+++/T3JyMu+88w733Xcf3377LT169HC7LJGAUzAUYLVZrh4OjRiy1vqzLBERn4wxc40xmT620WUeNp7DRwttBTpaa/sBfwFeNcY08fHyMUACMBC4GXjD+Gjit9ZOt9amWWvTkpKS/HZuocbbShboYKhXr1789NNPuu6IiEjY2LdvH5dffjmjRo2iVatWpKenM2XKFKKjo90uTcQVmmMowPwxYqiwsJDc3NyAzjMhIpHJWjv8aPcbY2KAMUD/Ms8pALztYYuNMetwRgell3v6ZuBt6yQOPxhjSoEWQOQOCzoKt0YMHX/88ezdu5ctW7bQrl27gB5bRETE3+bNm8dll13G5s2b+fvf/84///lPYmNj3S5LxFUaMRRg/phjCNA8QyISLIYDP1trN3t3GGOSjDHRnttdgO7Aeh/PfQc4zfO4HkAssLOuCw5V3mAokKuSgRMMAfz4448BPa6IiIg/5ebmcsMNN3D66acTFxfHd999xz333KNQSAQFQwFX21ayxMREQMGQiASNCzly0ukhwHLP8vVvAldba3cDGGOeNsakeR73LNDFGJMJvAZMtOpXqlBeXh7GGOLi4gJ6XAVDIiIS6hYuXEhqairTpk3j+uuvJyMjg4EDB7pdlkjQUCtZgPlrxJAmoBaRYGCtvczHvreAtyp4/JVlbhcCl9RZcWEmLy+P+Ph4fEzDVKeSkpJISkoiMzMzoMcVERGprcLCQu6++27uu+8+2rVrx9y5cxk2bJjbZYkEHQVDAeaPOYZAI4ZERCJNbm5uwOcX8kpOTtaIIRERCSkrVqxgwoQJZGRkcNlll/HII4/QtGlTt8sSCUpqJQswbytZbYMhjRgSEYkseXl5AZ9fyOv4449n5cqVWplMRESCXklJCffffz9paWls2bKFd999l+eee06hkMhRKBgKMO+IodosVw8aMSQiEmny8vJcGzF0/PHHk5OTw6ZNm1w5voiISFWsXbuWIUOGMGXKFM4991wyMzMZNWqU22WJBD0FQwFW21ayuLg4GjZsqBFDIiIRxu1gCDQBtYiIBCdrLU888QR9+vRh5cqVvPzyy8yaNYukpCS3SxMJCQqGAqy2wRA4o4Y0YkhEJLK4OceQgiEREQlWmzZt4swzz2Ty5MkMHjyYzMxMLr744oAv1iASyhQMBVhtl6sHZ54hBUMiIpHFzTmGEhISaNOmjVYmExGRoGGt5cUXXyQlJYX58+fz1FNP8dFHH9GuXTu3SxMJOQqGAsxfI4bUSiYiElncbCUDZ9SQRgyJiEgw2LFjB2PGjGHixImkpKSwbNkyJk2apFFCIjWkYCjA/BEMacSQiEjkCYZgaOXKlZSWlrpWg4iIyOzZs0lOTubDDz/kwQcf5Msvv6Rr165ulyUS0hQMBZg/WslatGjBzp07/VWSiIiEgGAIhnJzc9m4caNrNYiISOTKzs5mwoQJjBkzhg4dOrBkyRJuuukmoqOj3S5NJOQpGAowf4wYatmyJbt37z4YMomISPhzc/Jp0ATUIiLins8++4yUlBReffVV7rjjDr7//vuD1yURqT0FQwHmj2AoKSkJa63ayUREIoibk08DHHfccQCagFpEAs4YE22MWWqMmeN2LRJYBw4cYPLkyYwYMYJGjRqxYMEC7rrrrlr9LSUiR1IwFGDFxcUYY2o15LFly5aAM+maiIhEBrdbyZo1a0b79u0VDImIG24AfnK7CAms+fPn07dvX5588kn+/Oc/s2TJEgYMGOB2WSJhScFQgBUVFdVqfiFQMCQiEmlKSkooLCx0NRgC6N27N8uXL3e1BhGJLMaY9sBI4Gm3a5HAKCgoYMqUKQwePJji4mLmzZvH1KlTXb8GioQzBUMBVlRUVOuhjwqGREQiS15eHoDrvxT36dOHn376iYKCAlfrEJGI8ghwC6AlESNARkYGaWlp3H///Vx++eUsX76coUOHul2WSNhTMBRgxcXFfguGsrKy/FGSiIgEOW8w5OYcQwB9+/aluLiYlStXulqHiEQGY8w5wA5r7eJKHneVMSbdGJOu349DU3FxMffccw8DBgxg586dzJkzhxkzZtC4cWO3SxOJCAqGAswfrWTNmzcnOjpaI4ZERCJEMI0YAli2bJmrdYhIxDgZGGWM2QC8BpxujHm5/IOstdOttWnW2rSkpKRA1yi1tGrVKk4++WRuu+02fv/735OZmcnIkSPdLkskoigYCjB/tJJFRUXRokULBUMiIhEiWIKhbt26Ub9+fQVDIhIQ1tpbrbXtrbWdgAuBL6y1l7hclvhJaWkpjz76KH379mXt2rW8/vrrzJw5k8TERLdLE4k4tRu6ItXmj2AInHYyBUMiIpEhNzcXcD8Yio6Opnfv3mRkZLhah4iIhLaNGzfyhz/8gXnz5jFy5EhmzJhBmzZt3C5LJGJpxFCAFRcX17qVDBQMiYhEkmAZMQROO9myZcuw1rpdiohEEGvtl9bac9yuQ2rHWsuzzz5LSkoKixYt4umnn+b9999XKCTiMgVDAebPEUOaXE9EJDIEy+TT4ARDe/bsYdOmTW6XIiIiIWTbtm2MGjWKK664gtTUVFasWMEVV1yBMcbt0kQiXq2CIWPM68aYDM+2wRiT4dl/Qpn9y4wx5/ul2jCgVjIREakuf44YWroUxo+HtDS45BJIT6/e8/v27QtoAmoREam6WbNmkZyczNy5c3n44Yf54osv6NSpk9tliYhHrYIha+04a21fa21f4C3gbc9dmUCaZ/9ZwP+MMZrPCP+2ku3bt4/8/Hw/VCUiIsHMX8HQ66/DCSfAp59Cixbw4YfO5w89BFXtDEtJSQEUDImISOV2797NRRddxNixY+nSpQtLly7lxhtvJCpKjSsiwcQv35HGGf83FpgJYK3NtdYWe+6OBzQRgYe/Rgx5l+JUO5mISPjzx+TTS5Y4I4QGDoQ1a+Djj+GXX+CCC+Cmm+DRR6v2Oo0bN6Zbt26agFpERI7qo48+Ijk5mVmzZvGvf/2L+fPn06tXL7fLEhEf/BXVDga2W2vXeHcYY040xvwIrACuLhMUHcYYc5UxJt0Ykx4JIYc/W8kAtZOJiESA2s4xVFwMl10GSUnw3nuQkODsb9rUGUU0Zgz85S8wd27VXs87AbWIiEh5OTk5TJo0ibPPPpuEhAQWLlzIbbfd5peuCRGpG5UGQ8aYucaYTB/b6DIPG49ntJCXtXahtfZ4YABwqzEm3tfrW2unW2vTrLVp3lEw4ay4uNivwVAkhGkiIpGutq1kr7wCK1bAY49B8+aH3xcVBS++CD17wh/+ANnZlb9enz59WLduHTk5OTWqR0REwtPXX39Nnz59mDFjBrfccguLFy8mNTXV7bJEpBKVBkPW2uHW2mQf27sAnrmDxgCvV/D8n4D9QLI/Cw9VRUVFfptjCDRiSEQkEtQmGCopgXvugb59nZFBvjRsCC+8AFu3wt/+Vvlr9u3bF2sty5cvr3Y9IiISfvLz8/nrX//KqaeeSlRUFN988w33338/cXFxbpcmIlXgj1ay4cDP1trN3h3GmM7eyaaNMccAvYANfjhWyFMrmYiIVFdubi5RUVE1un58/rkzp9CUKXC0FYFPOAEmT4ann3ZGFx2N993fJUuWVLseEREJL+np6aSmpjJ16lSuvvpqMjIyOPnkk90uS0SqwR/B0IWUayMDTgGWeZavnw1ca63d6YdjhTx/BUONGjUiPj5ewZCISATIy8ujfv36mKMlOxV45hlITITzzqv8sf/8pzPv0F/+cvRVytq2bUvr1q1Jr+5a9yIiEjaKioq48847GThwIPv27eOTTz7hiSeeoFGjRm6XJiLVVOueJmvtZT72vQS8VNvXDkf+Wq7eGEPLli0VDImIRIC8vLwaTTy9Zw+88w5ccw1UZTR/QgLccQf8+c/w1Vdw6qm+H2eMIS0tTcGQiEiEWrlyJRMmTGDx4sVccsklTJs2jeblJ7ETkZDhr1XJpIr8NWIInCXrNfm0iEj4844Yqq4PPoDCQrjooqo/Z9IkaN0a/u//jv64tLQ0fvrpJ/bv31/tukREJDSVlJTw0EMPkZqaysaNG3nrrbd46aWXFAqJhDgFQwHmz2CoZcuWbN++3S+vJSIiwaumwdB770GbNpCWVvXn1K8PN93kzE20YEHFj0tLS8Nay9KlS6tdl4iIhJ7169dz2mmncdNNN3HWWWeRmZnJmIpWNRCRkKJgKMD81UoG0Lp1a7Zt2+aX1xIRkeCVm5tb7WCooAA+/hjOPddZkr46rr7amZfoaKOG+vfvD6B2MhGRMGetZfr06fTu3Ztly5bx/PPPM3v2bFq1auV2aSLiJwqGAsyfI4batGnD9u3bKS0t9cvriYhIcKrJHEPffgs5OXDOOdU/XsOG8Kc/wYcfws8/+35M69atad++vYIhEZEwtmXLFkaOHMmkSZMYOHAgK1asYOLEiTVaDEFEgpeCoQDzdzBUXFzMzp1a8E1EJJzVpJXsyy8hOhqGDq3ZMSdNgthY+O9/K36MJqAWEQlP1lpeffVVkpOT+fLLL/nvf//Lp59+SseOHd0uTUTqgIKhACsuLvZbMNS6dWsAtm7d6pfXExGR4FTTYKh/f2jSpGbHbNUKxo+H55+HvXt9PyYtLY3Vq1ezt6IHiIhIyNm5cydjx47l4osvpmfPnmRkZDB58mSiqtuXLCIhQ9/dAVZUVOS3OYbatGkDoHmGRETCXHXnGMrNhR9+qPloIa/rr4cDB+C553zfn+aZ1XrJkiW1O5CIiASF999/n+TkZN59913uu+8+vv32W3r06OF2WSJSxxQMBZi/W8lAI4ZERMJddecY+v57Z5n6U0+t3XH794eTT3bayXxNZ6cJqEVEwsO+ffu4/PLLGTVqFK1atSI9PZ0pU6YQHR3tdmkiEgAKhgJMwZCIhAtjzOvGmAzPtsEYk+HZ38kYk1fmvqcqeH5fY8z3nsekG2NOCOgJhJDqtpJ98w0Y44Q6tXXttbBuHcybd+R9LVq0oFOnTgqGRERC2Lx580hJSeGFF17g73//O4sWLaJ3795ulyUiAeSfniapktLSUkpLS/3WSla/fn2aNm2qYEhEXGGtHee9bYx5CCg70cw6a23fSl7iAeAua+1HxpizPZ+f6u86w0F1g6FFi+DYY6Fp09ofe8wYaN4cnn4ahg078v60tDR++OGH2h9IREQCKjc3l1tvvZVp06bRo0cPvvvuOwYOHOh2WSLiAo0YCqDi4mIAv40YAmcCagVDIuIm46xZOxaYWc2nWsA7NXJTYIs/6won1QmGrIX0dBgwwD/Hjo+HSy+Ft98GX4tgDhw4kA0bNmi+OxGRELJw4UL69evHtGnTuP7661m6dKlCIZEIpmAogIqKigD/BkNt2rTRL+Mi4rbBwHZr7Zoy+zobY5YaY74yxgyu4Hk3Ag8aYzYB/wFu9fUgY8xVnlaz9KysLL8WHgqKioooLi6ucjD022+wfTt45oX2iyuvdOYsevnlI+8bNGgQAAsWLPDfAUVEpE4UFhZy2223cdJJJ5GXl8fcuXOZNm1ateaxE5Hwo2AogLwjhvzVSgZOMKQRQyJSV4wxc40xmT620WUeNp7DRwttBTpaa/sBfwFeNcb4WjT9GuDP1toOwJ+BZ3zVYK2dbq1Ns9amJSUl+efEQkheXh5AlX9pX7TI+eivEUMAKSlw4okwY4YzIqms1NRU6tWrp2BIRCTIrVixghNPPJF77rmHiRMnsmLFCob56hEWkYijYCiA6mrE0NatW7Hlf1MXEfEDa+1wa22yj+1dAGNMDDAGeL3Mcwqstbs8txcD6wBfa91OBN723J4FaPJpH7zBUFVHDKWnQ0wM9Onj3zquvBJWrnRWPCsrPj6e1NRUBUMiIkGqpKSE+++/n/79+7Nlyxbeffddnn32WZr6YyI6EQkLCoYCqC6CodatW5Obm0tOTo7fXlNEpBqGAz9bazd7dxhjkowx0Z7bXYDuwHofz90CDPXcPh1Y4+MxEa+6wdCiRc4In/h4/9Zx4YXQqJEzCXV5J510Eunp6RQWFvr3oCIiUitr165lyJAhTJkyhdGjR/Pjjz8yatQot8sSkSCjYCiA6mLyae+S9ZpnSERcciFHTjo9BFjuWb7+TeBqa+1uAGPM08YY7+w3fwQeMsYsA+4FrgpMyaElNzcXqHowtGwZ9Ovn/zoaNYL/9/9g1izwlHTQoEGDyM/PZ9myZf4/sIiIVFtpaSmPP/44ffr0YeXKlbzyyiu88cYbtGjRwu3SRCQIKRgKIO+IIX/PMQRoniERcYW19jJr7VPl9r1lrT3eWtvXWptqrX2/zH1XWmvTPbe/tdb2t9b2sdae6Gk7k3KqM8fQjh3OlpJSN7VMmAA5OfDOO4fv1wTUIiLBY9OmTZx55plcd911DB48mMzMTC666CKcRURFRI6kYCiA6mqOIVAwJCISrqrTSvbjj87H5OS6qWXIEDjmGHjxxcP3t2/fnvbt2ysYEhFxkbWWF198kZSUFBYsWMBTTz3FRx99RLt27dwuTUSCnIKhAKqLVrLWrVsDCoZERMJVdYKhzEznY10FQ1FRcOml8NlnsGXL4fcNGjRIwZCIiEt27NjBmDFjmDhxIikpKSxbtoxJkyZplJCIVImCoQCqi1ay5s2bExcXp2BIRCRMVWeOocxMSEyEVq3qrp5LL4XSUnj11cP3Dxo0iI0bN+p6JCISYLNnzyY5OZkPP/yQBx98kC+//JKuXbu6XZaIhBAFQwFUF61kxhjatm3LlvJv3YqISFio7oih5GSoyzeIe/SAgQPhhRfA2kP7vfMMfffdd3V3cBEROSg7O5sJEyYwZswYOnTowJIlS7jpppuIjo52uzQRCTEKhgKoLoIhcOZ22LRpk19fU0REgkNVgyFrDwVDdW3CBOdYGRmH9qWmptKgQQO++eabui9ARCTCffbZZ6SkpPDqq69yxx138P3333P88ce7XZaIhCgFQwHknWPIn61k4ARDmzdv9utriohIcMjPzwcqD4Y2b4Z9+yAQfxeMGwexsYdPQh0bG8ugQYP46quv6r4AEZEIdeDAASZPnsyIESNo1KgRCxYs4K677vL7G88iElkUDAVQXY0Y6tChA5s3b8aWHdMvIiJhoaojhn7+2fl47LF1XREkJMC55zrzDHkubQAMHTqU5cuXs2fPnrovQkQkwsyfP5++ffvy5JNP8pe//IUlS5YwYMAAt8sSkTCgYCiA6rKVrLCwkKysLL++roiIuK+qwdDq1c7Hnj3ruiLHhAmwYwd8+umhfUOGDMFay7fffhuYIkREIkBBQQFTpkxh8ODBFBcXM2/ePB566KEqzT0nIlIVCoYCqC5byQC1k4mIhKG8vDyio6MrfVNh9Wpo1Ahatw5MXWedBS1aOJNQe5144onExsby9ddfB6YIEZEwl5GRQVpaGvfffz9XXHEFy5cvZ+jQoW6XJSJhRsFQANXliCFQMCQiEo7y8vKq9K7wmjXQvXvdrkhWVmwsXHQRvPceZGc7++Lj4znxxBM1z5CISC0VFxfzf//3fwwYMIBdu3bxwQcfMH36dBo3bux2aSIShhQMBVBdzjEEaGUyEZEwVNVgaPVqZyn5QJowAQoK4I03Du0bOnQoS5YsIScnJ7DFiIiEiVWrVnHyySdz++238/vf/57MzEzOPvtst8sSkTCmYCiAvK1k/g6GWrZsSUxMjEYMiYiEoby8POLj44/6mMJC+OUXZ8RQIKWmwnHHHb462ZAhQygpKWH+/PmBLUZEJMSVlpby6KOP0rdvX9auXcvrr7/OzJkzSUhIcLs0EQlzCoYCyDtiyN9zDEVFRdGuXTsFQyIiYagqI4Z++QVKSwM/YsgYZ9TQd9/BunXOvpNOOomYmBi1k4mIXxhj4o0xPxhjlhljfjTG3OV2TXVh48aNDB8+nBtvvJFhw4aRmZnJ2LFj3S5LRCKEgqEAqqtWMnDaydRKJiISfqoSDHlXJAt0MARw8cVOQPTSS87nDRs2pH///nz55ZeBL0ZEwlEBcLq1tg/QFzjLGDPQ3ZL8x1rLs88+S0pKCosWLeLpp5/m/fffp02bNm6XJiIRRMFQANVlMNS+fXuNGBIRCUNVCYbWrHE+BrqVDKB9exg2zGkns9bZN2zYMH744Qf27dsX+IJEJKxYx37Pp/U8m3WxJL/Ztm0bo0aN4oorrqB///6sWLGCK664AhOoVQRERDwUDAVQXS1XD4eCIWvD4jopIiIeVR0xlJgIbk1DMWGC08723XfO52eccQYlJSUaNSQifmGMiTbGZAA7gM+stQtdLqnWZs2aRXJyMnPnzuWRRx7h888/p1OnTm6XJSIRSsFQAFU2Ysha+M9/4IQTYMoUZ6WXqmrfvj0FBQXs3LnTH6WKiEiQqEowtHYtdOsWoIJ8OP98aNjw0CTUgwYNokGDBsydO9e9okQkbFhrS6y1fYH2wAnGmOSy9xtjrjLGpBtj0rOyslypsap2797NRRddxNixY+nSpQtLly7lhhtuICpKf5aJiHtq9RPIGPO6MSbDs23wJPll7+9ojNlvjLmpVlWGicqCoSeegJtvhgMH4P774Y9/PDQsvzJasl5EJDxVdfLpLl0CVJAPjRrB738Pr78OeXkQFxfHkCFD+Oyzz9wrSkTCjrU2G5gHnFVu/3RrbZq1Ni0pKcmV2qrio48+Ijk5mVmzZvGvf/2L+fPn06tXL7fLEhGpXTBkrR1nre3rSfDfAt4u95CpwEe1OUY4OVorWV4e3HEHDB8OmZnwz386E3nOmVO11z7mmGMAZ0UDEREJH/n5+UcNhkpK4Ndfwe0OhAkTYN8+eO895/Phw4fz888/a/47EakVY0ySMaaZ53Z94AzgZ1eLqqacnByuuuoqzj77bBITE/nhhx+47bbb6mR6CRGRmvDLmEXjzJA2FphZZt95wC/Aj/44Rjg42oihN96A3bvh7393Vnf5xz+gVy+46Sbnl/7KeHuSN2zY4MeKRUTEbZWNGPrtNygudj8YOvVU6NDhUDvZGWecAaB2MhGprTbAPGPMcmARzhxDVXzr1H1ff/01ffr04emnn+aWW24hPT2dfv36uV2WiMhh/NXMOhjYbq1dA2CMaQT8DbirsieGUk9wbRUVFREdHe1zpYE334TOnZ1frAHq1YO77nImFH3//cpfOyEhgcaNG/PLL7/4t2gREXFVZcGQ9/2Azp0DU09FoqLgkkvgk09g2zZISUmhZcuWCoZEpFastcuttf2stb2ttcnW2rvdrqkq8vPz+etf/8qpp55KVFQU33zzDffffz9xcXFulyYicoRKgyFjzFxjTKaPbXSZh42nzGgh4E7g4TJLS1YoVHqC/aG4uNjnaKGCApg3D373O2e0kNeYMc47wFOnVv7axhg6deqkEUMiImGmsmDI+36A2yOGAC691BnlOnOmc10aPnw4c+fO1YqZIhJR0tPTSU1NZerUqVx99dVkZGRw8sknu12WiEiFKg2GrLXDPel8+e1dAGNMDDAGeL3M004EHjDGbABuBP5ujLnO/+WHlqKiIp+9xN9950w4fdZZh++PiYFrroFvvoE1ayp//c6dO2vEkIhIGCktLaWgoID4+PgKH7Nhg/OmQseOgaurIsceCwMGHN5Otn37dlasWOFuYSIiAVBUVMSdd97JwIED2bdvH5988glPPPEEjRo1crs0EZGj8kcr2XDgZ2vtwdklrbWDrbWdrLWdgEeAe621//XDsUJaUVGRzxFD8+c7H4cOPfI5l1ziDM9/4YXKX987YkjvzIqIhIf8/HyASkcMtW0LwdKdMGECZGTA8uUwYsQIAD788EN3ixIRqWMrV65k0KBB3HXXXVx00UWsWLHi4M9AEZFg549g6EIObyOTClTUSpaeDj16QJMmRz6nbVsYMcIJhkpLj/76nTt3Zv/+/ezatctPFYuIiJvy8vKAowdDGza4P79QWRde6MyT9+KL0LZtW/r168cHH3zgdlkiInWipKSEhx56iNTUVDZu3Mhbb73Fiy++SPPmzd0uTUSkymodDFlrL7PWPnWU+++01v6ntscJBxW1ki1eDP37V/y8Sy+FzZsPjSyqiFYmExEJL1UJhn75JTjmF/Jq0QJGjoRXXnFWSxs5ciTz589n9+7dbpcmIuJX69ev57TTTuOmm27irLPOIjMzkzFjxrhdlohItflrVTKpAl+tZDt2OKHP0YKhc86B2Fh4++2jv35nz1vGCoZERMJDZcFQUZFzDQmmEUPgtJNt2wZz5zrBUGlpKZ988onbZYmI+IW1lunTp9O7d2+WLVvGCy+8wOzZs2nVqpXbpYmI1IiCoQDyFQx55+Ps27fi5zVpAmec4QRDR5s+yDtiSBNQi4iEh8qCoU2bnDbjYBoxBHD22ZCQ4LSTDRgwgBYtWqidTETCwpYtWxg5ciSTJk1i0KBBZGZmMmHCBEzZpYVFREKMgqEAKi4uPqKVzLvaWM+eR3/uBRfAxo1O21lFmjZtSvPmzTViSEQkTFQWDHnfBwi2EUNxcc5cQ7Nnw4ED0fzud7/j448/pqSkxO3SRERqxFrLq6++SnJyMl9++SX//e9/+eSTT+jQoYPbpYmI1JqCoQDyNWJo9WqoX9+ZZPpoRo2C6Gh4662jP05L1ouIhI/KgiHv+wDBNmIInHay/Hx44w2nnWzXrl0sXLjQ7bJERKpt586djB07losvvphevXqxbNkyJk+eTFSU/pQSkfCgn2YBVFEw1L27syT90SQmwqmnwnvvHf1xnTp1UjAkIhImKguGNm4EY6B9+0BWVTUnnADHHgtPPw1nnnkm0dHRaicTkZDz/vvvk5yczLvvvst9993HN998Q/fu3d0uS0TErxQMBZCv5epXr3aWqq+KkSNh5cpDrQO+dOvWjfXr12u4vohIGMjPzwcqDoZ++w1at3aWhw82xsAf/wgLF8KmTc04+eSTmTNnjttliYhUyb59+7j88ssZNWoUrVu3Jj09nSlTphAdHe12aSIifqdgKIDKL1dfVATr11cvGAI42huuPXr0oLCwkF9//bUWlYqISDCobMTQb79Bu3aBrKh6Lr3UWVXz6afhvPPOY/ny5axbt87tskREjmrevHmkpKTwwgsv8Pe//50ffviB3r17u12WiEidUTAUQOVbyTZuhJISp5WsKnr0cB5bWTAEsHr16tqUKiIiQcAbDMXHx/u8f/Pm4A6GWrSA88+Hl16Cs88+H4DZs2e7XJWIiG+5ubnccMMNnH766cTHx/Pdd99xzz33EBsb63ZpIiJ1SsFQAJVvJdu40fl4zDFVf42RI2HePDhwwPf93p7nNd7lzkREJGSF+oghgCuvhD17YPHiTqSmpvL222+7XZKIyBEWLlxIv379mDZtGn/6059YunQpAwcOdLssEZGAUDAUQOVbyTZtcj5WZ5XLkSOhoAC++ML3/a1ataJx48YaMSQiEgaOFgzl5kJ2dvAHQ6efDp07O+1kY8aMYcGCBWzZssXtskREDrLWcuONN5Kfn8/nn3/Oo48+SoMGDdwuS0QkYBQMBVD5VrLNm52P1VlNZsgQaNSo4nYyYww9evRQMCQiEgaOFgz99pvzMRhXJCsrKgquuMIZ7dq//xgA3nnnHXeLEhEpwxjDa6+9xvLlyzn99NPdLkdEJOAUDAVQ+WBo0yZISoIKpo7wKTYWzjjDCYas9f2Y7t27q5VMROqcMeZ1Y0yGZ9tgjMkoc19vY8wCY8yPxpgVxpgjftIZYxKMMZ8ZY9Z4PjYP6AmEgLy8POrVq+dzFRxvMBTsI4YA/vAHJyD66qtj6dWrl9rJRCToHHPMMTRt2tTtMkREXKFgKICKi4uPaCWrThuZ18iRzmij5ct939+jRw82bNhAQUFBDSsVEamctXactbavtbYv8BbwNoAxJgZ4GbjaWns8cCpQ5OMlpgCfW2u7A597Ppcy8vLyjjq/EIRGMNS2rXPteu45GD16DF9++SW7du1yuywRERERQcFQQPkaMVSTYOjss52PFbWT9ejRg9LSUtavX1+DKkVEqscYY4CxwEzPrhHAcmvtMgBr7S5rbYmPp44GXvDcfgE4r45LDTnhEgwB/PGPsH07JCaOoaSkhHfffdftkkREREQEBUMB5SsYqsncEG3aQP/+FQdDWplMRAJsMLDdWuv9odMDsMaYT4wxS4wxt1TwvFbW2q2e29uAVnVdaKipLBhq3NjZQsHvfue8GfLxx6l07dqVmTNnVv4kEREREalzCoYCqGwrWU4O7N1bsxFD4AzJ//572LnzyPt69OgBwKpVq2paqogIAMaYucaYTB/b6DIPG8+h0UIAMcApwMWej+cbY4Yd7TjWWgv4nDnNGHOVMSbdGJOelZVVuxMKMUcLhjZvDp3RQgAxMXD11fDFF4YRIy7iiy++YOvWrZU/UURERETqlIKhACo7Ysj7u3BNf6kfORJKS+Hjj4+8r1mzZrRu3ZqffvqphpWKiDistcOttck+tnfh4HxCY4DXyzxtM/C1tXantTYX+BBI9fHy240xbTyv0wbYUUEN0621adbatKSkJH+eXtCrbMRQKAVDAFde6SyisHfveEpLS3njjTfcLklEREQk4ikYCqCywdD27c6+VjVsnEhLg5YtK24nO/7448nMzKzZi4uIVN1w4Gdr7eYy+z4BUowxDTzB0VBgpY/nvgdM9NyeCGjSmXLy8vKIr2Dpyt9+C/6l6str2RLGjoX33z+WPn368eqrr7pdkoiIiEjEUzAUQMXFxQeDoR2e98VbtqzZa0VFOZNQf/wxFBcfeX9ycjI//vgjpaWlNaxWRKRKLuTwNjKstXuAqcAiIANYYq39AMAY87QxJs3z0H8DZxhj1uAETP8OVNGhIj8/3+eIoZISZ+RpqI0YApg82Wmn7tLlIn744QfNhyciIiLiMgVDAVRUVHRwjiFvMFTTEUPgtJNlZ8OCBUfel5ycTG5uLhs3bqz5AUREKmGtvcxa+5SP/S9ba4/3tJ3dUmb/ldbadM/tXdbaYdba7p6Wtd2BrD0UVNRKtmOHEw6FYjB04onOAgqZmRdijNEk1CIiIiIuUzAUINZan61kLVrU/DXPOMOZzHPOnCPvS05OBlA7mYhICKsoGAq1perLMsYZNbRmTXv69BnKK6+8gjP3uIiIiIi4QcFQgJSUlAAc1kqWmOgEOzXVtCkMGeJ7nqHjjjsOUDAkIhLKwjEYArjwQkhIgKioS1i9ejULFy50uyQRERGRiKVgKECKPRMBeVvJtm+vXRuZ18iR8OOPUL5jrEmTJnTs2JEff/yx9gcRERFXhGswVL++s0LZ0qVjadCgIc8884zbJYmIiIhELAVDAVJUVAQcPmKophNPlzVypPPR16ih5ORkjRgSEQlhFQVDmzc7I079cR1xy/XXQ3R0Yzp3Hsdrr73G/v373S5JREREJCIpGAqQ8sGQv0YM9egBXbtWPM/QTz/9dHC0koiIhJajjRhq0waio10oyk/at4fx42HduivYv38/b7zxhtsliYiIiEQkBUMBUr6VzF8jhoyBc86BefMgN/fw+5KTkyksLGT16tW1P5CIiARUcXExxcXFFQZDodpGVtZf/wr5+YNISuqldjIRERERlygYCpCyI4YKCmDvXv+1AIwcCfn58MUXh+/v168fAEuWLPHPgUREJGDy8vIAwjoY6tMHhg83FBRcyfz58/npp5/cLklEREQk4igYCpCywdCOHc4+f7SSgbMyWcOGR84z1KtXL+rXr8/ixYv9cyAREQmYSAiGAG66Cfbtu5To6Biefvppt8sRERERiTgKhgLE20pWNhhKSvLPa8fFwYgRzjxD1h7aHxMTQ9++fRUMiYiEIG8wFB8ff9j+nBxnC5dgaMQISE5uSaNGY3j22Wc5cOCA2yWJiIiIRBQFQwHiHTEUExPD7t3OvhYt/Pf655zjrFKTkXH4/v79+7N06VJKS0v9dzAREalzFY0YCvWl6sszxhk1tHfv9WRnZ/PKK6+4XZKIiIhIRFEwFCBlW8m8wVBCgv9e/9xzISoK3nrr8P39+/dn//79moBaRCTEVBYMtW8f6Irqzvjx0KHDyTRo0JfHHnsMW3b4q4iIiIjUKQVDAVK2lcwbDDVv7r/XT0qCoUPh7bcP39+/f38AtZOJiISYioKhzZudj+EyYgggNhZuvdWQm3s9mZmZfPXVV26XJCIiIhIxFAwFSNlWsj17nH3+DIYAxoyBn35yNq9jjz1WE1CLiISg/Px8IPxbybz+8Ado02Y8MTGJPPbYY26XIyIiIhIxFAwFSPlWsgYNoNx8orV2/vnOx7KjhmJiYujTpw/p6en+PZiIiNSpo7WSNW8OPhYrC2nx8fC3v9WnuPhK3nnnHTZu3Oh2SSIiIiIRoVbBkDHmdWNMhmfbYIzJ8OzvZIzJK3PfU36pNoSVD4b8Ob+QV7t2MGjQkfMMDRo0iEWLFlFYWOj/g4qISJ04WjAUbqOFvP74R0hMnIy1UTz88MNulyMiLjPGdDDGzDPGrDTG/GiMucHtmkREwlGtgiFr7ThrbV9rbV/gLaDsDDfrvPdZa6+uzXHCgXeOIe+qZP5uI/MaMwaWLoVffjm07+STTyY/P5+lS5fWzUFFRMTvIjEYatAAbrmlA9ZezP/+N4OdO3e6XZKIuKsY+Ku19jhgIDDZGHOcyzWJiIQdv7SSGWMMMBaY6Y/XC0dlRwzt2VM3I4bACYbg8Hayk08+GYBvv/22bg4qIiJ+l5ubC0CDBg0O2x/OwRDAtddC8+a3kJ+fy3//+1+3yxERF1lrt1prl3hu5wA/AWH8E1BExB3+mmNoMLDdWrumzL7OxpilxpivjDGDK3qiMeYqY0y6MSY9KyvLT+UEn0C0kgF06QL9+sHrrx/a17p1a7p166ZgSEQkhPgKhoqLYfv28FqqvrxGjeCf/zwOGMXUqY+xf/9+t0sSkSBgjOkE9AMWulyKiEjYqTQYMsbMNcZk+thGl3nYeA4fLbQV6Git7Qf8BXjVGNPE1+tba6dba9OstWlJSUm1OZegVn65+roKhgAuuggWLYI1ZWK6k08+me+++w5rbd0dWERE/MZXMLRtG5SWhveIIYCrr4ZWraaQk7ObGTOedrscEXGZMaYRzrQVN1pr9/m4PyLeaBYRqSuVBkPW2uHW2mQf27sAxpgYYAzwepnnFFhrd3luLwbWAT3q5hRCQ9nl6us6GBo/HoyBV145tO+UU04hKyuLNWXTIhERCVq5ubnExMQQGxt7cN/mzc7HcA+G4uLgvvsGAUO4556HKCgocLskEXGJMaYeTij0irX2bV+PiZQ3mkVE6oo/WsmGAz9bazd7dxhjkowx0Z7bXYDuwHo/HCtkeYOhkpJ65OfX3eTT4PzBcNppTjDkHSB0yimnAPDNN9/U3YFFRMRvDhw44HN+IQj/YAjg0kuhQ4fb2LVrM089NcPtckTEBZ55TJ8BfrLWTnW7HhGRcOWPYOhCjpx0egiw3LN8/ZvA1dba3X44VsjytpLl5tYD6nbEEMDFF8PatfDDD87nPXv2pHXr1sydO7duDywiIn6Rm5sb0cFQTAw8/PBwYAj//Oc9B1vrRCSinAxcCpxujMnwbGe7XZSISLipdTBkrb3MWvtUuX1vWWuP9yxVn2qtfb+2xwl13hFDOTkxQN0HQxdc4AzF97aTGWM444wzmDt3LqWlpXV7cBERqbWKgqHYWGjRwqWiAmzMGMNxx/2LvXu38fDDT7hdjogEmLX2W2utsdb29vxd0dda+6HbdYmIhBt/rUomlfAGQ/v3B2bEUNOmcO658Npr4Dk0I0aMYOfOnWRkZNTtwUVEpNZyc3Np2LDhYfu8S9Ub41JRAWYMPPfcEGAE99xzPzk5OW6XJCIiIhJ2FAwFSPlgqC7nGPK69FLIyoKPPnI+Hz58OACfffZZ3R9cRERqpaI5hiKhjaysE06As8/+F3l5O/nnPx91uxwRERGRsKNgKEC8cwzt3RuYVjKAs8+GNm1g+nTn89atW9O7d28+/fTTuj+4iIjUSkWtZJEWDAHMmHEC0dHn8dhjD7B9+3a3yxEREREJKwqGAuTQHEOBGzEUEwNXXOGMGPr1V2ffiBEj+Pbbb9m/f3/dFyAiIjVWPhiyNnKDobZt4U9/up/i4jyuuOJ2t8sRERERCSsKhgKkfCtZ48aBOe4VVzh/TDz7rPP5yJEjKSws5JNPPglMASIiUiPlg6HsbMjNjcxgCODee3vQpMl1fPDBM6SnL3e7HBEREZGwoWAoQLytZDk5UTRuDFEB+pfv1AlGjIBnnoHiYjjllFNITExk9uzZgSlARERqpPzk05G0VL0v8fEwffodQDPGjfsL1lq3SxIREREJCwqGAqSoqIh69eqxb5+hadPAHvuqq2DzZqelLCYmhlGjRjFnzhwKCwsDW4iIiFRZ+cmnIz0YAhg3rjn9+t3J+vWf89RT77tdjoiIiEhYUDAUIIeCIWjSJLDHPvdcZ36Gxx5zPj///PPZu3cv8+bNC2whIiJSZeVbybzBUPv2LhUUJN5552qioo7jz3++nv37D7hdjoiIiEjIUzAUIMXFxa4FQ/XqwXXXwWefwYoVcMYZZ9CwYUPefvvtwBYiIiJVUlpaSl5ens9gqG1bl4oKEh071uPPf/4fBQW/MnbsXW6XIyIiIhLyFAwFSFFRETExMezdS8BbyQAmTYIGDeCRRyA+Pp5Ro0Yxa9YsCgoKAl+MiIgcVX5+PsARwVBSEsTGulVV8HjggVNo2fIKPvpoKnPnaiJqERERkdpQMBQgbraSASQkwMSJ8PLLsH07TJgwgT179jBnzpzAFyMiIkd14IDTIlV+8ulInl+orKgoeO+9B4AE/t//m0RJSanbJYmIiIiELAVDAeINhtwaMQRw441QWAiPP+60k7Vp04YXX3zRnWJERKRCubm5wOEjhjZvVjBU1oknJnDppQ+Rnf09F1/8uNvliIiIiIQsBUMBUlxcTExMjGsjhgB69IDzznMmod6/P5pLLrmEDz/8kKysLHcKEhERn3wFQxoxdKTnnruExMTf8frrf2Pu3NVulyMiIiISkhQMBYh3xNCBA+4FQwC33w7Z2TBtGkycOJHi4mKef/559woSEZEjlA+GCgogK0srkpUXHW34+ONnMKY+558/gfz8YrdLEhEREQk5CoYCpKioiKioeoB7rWQAqanO8vUPPwwdOhzPkCFDeOKJJygpKXGvKBEROUz5YGjrVme/RgwdKS2tDddd9wT79y/kzDPvd7scERERkZCjYChAiouLMSYGcHfEEMA//wl79jgtZX/605/YsGGDJqEWkWozxrxujMnwbBuMMRll7uttjFlgjPnRGLPCGBPv4/kPGmN+NsYsN8bMNsY0C2T9waz85NPepeoVDPk2bdo4unQZx9df38mjj6a7XY6IiIhISFEwFCDBMmIIoH9/OOcc+M9/4JRTRtOhQwceffRRd4sSkZBjrR1nre1rre0LvAW8DWCcFPxl4Gpr7fHAqUCRj5f4DEi21vYGVgO3BqLuUFB+xNDmzc5+tZJV7Ouvn6Bevbb85S9jyczMdrscERERkZChYChAioqKMMYJhtweMQRw332wbx/8+98xXHfddcybN48ffvjB7bJEJAQZYwwwFpjp2TUCWG6tXQZgrd1lrT2iX9Va+6m11jspzPeAYg8PBUPV165dAi+++BqlpZsYOvRyCgqs2yWJiIiIhAQFQwFSXFwMBE8wlJwMV1zhLF0/YsQ1JCQkcPfdd7tdloiEpsHAdmvtGs/nPQBrjPnEGLPEGHNLFV7jcuAjX3cYY64yxqQbY9IjZRVFX8FQo0bBcf0IZhdeOIiLL/43u3fP5vTTp2GVDYmIiIhUSsFQgBQVFWGtM8eQ261kXnffDbGx8K9/Neavf/0rH3zwAYsWLXK7LBEJIsaYucaYTB/b6DIPG8+h0UIAMcApwMWej+cbY4Yd5Rj/AIqBV3zdb62dbq1Ns9amJSUl1fqcQoE3GPLOMbR5szNayBg3qwoNL730F7p1O5f582/iL3/52u1yRERERIKegqEAcYKhWCB43vFt3Rr+9jd4+2047rjraN68OXfccYfbZYlIELHWDrfWJvvY3oWD8wmNAV4v87TNwNfW2p3W2lzgQyDV1+sbYy4DzgEutlbjO7z2798PHD5iSG1kVWOMYeHCF2nYsCuPPHIBL7/8i9sliYiIiAQ1BUMBUlhYiLXBMfl0WTffDD16wF//2oSbb/47H3/8MR995LObQ0TEl+HAz9bazWX2fQKkGGMaeIKjocDK8k80xpwF3AKM8gRI4pGTk0NsbCyxsc4bCgqGqichoRlff/0eUVHFTJw4isWLc9wuSURERCRoKRgKkKKiIkpL6xEdDZ43gINCfDw89RSsXw979/6J7t278+c//5miIl8LCImIHOFCDm8jw1q7B5gKLAIygCXW2g8AjDFPG2PSPA/9L9AY+Myz5P1TAas6yO3fv5/GjRsDUFwMW7cqGKqu1NQePPfcG5SW/sTgwRfz669HzH8uIiIiIigYChhvMNSkSfDNEXHaaTBxIjz0UCx/+tNUVq1apeXrRaRKrLWXWWuPCHSstS9ba4/3tJ3dUmb/ldbadM/tbtbaDt4l7621Vwey9mCWk5NDo0aNANi+HUpKFAzVxIQJZ/C3vz1KXt779O07md271a0oIiIiUp6CoQApLCykpCQ2aOYXKu8//4GEBJg+fSQjR57L7bffzurVq90uS0QkIuXk5BwcMfTbb84+BUM18+9/T2b8+FvZs+d/9OnzT3LVtCgiIiJyGAVDAVJUVERxcb2gDYZatIDnnoMVKwzt2j1F/fr1+cMf/kBJiYbei4gEWtlWss2e2ZsUDNXcK6/cw+mnX8Hmzf+iX7/HyMtzuyIRERGR4KFgKEC8wVAwTTxd3tlnw7XXwvTpbbnqqmnMnz+f+++/3+2yREQiTtlWMgVDtWeM4ZNPniI19TxWr/4TqanTFQ6JiIiIeCgYCpCioiKKioK3lczrwQfh2GPhmWcuZtSoC7n99tv5/PPP3S5LRCSilG0l27zZWSggIcHlokJcTEwM3303k969z+bnnyeRmvqUwiERERERFAwFTGFhIUVFwT1iCJwV095+GwoKDJs3z6BHj56MHz+eX3/91e3SREQiRvlWsvbtg2/hglAUHx/PDz+8TZ8+5/Dzz9fQp88T7N3rdlUiIiIivhUUwMaN8P33zt/pjz8OOTn+P06M/19SfCkqKsLa4J1jqKxeveDFF+H88xtx3nlv8/nnJ3D22WfzzTff0Lx5c7fLExEJe+VbydRG5j9xcXEsXPgmJ500liVLJnPssXtZvHgKbdooeRMREZHAyM+HrVudbcuWI297P+7adeRzhwyBlBT/1qNgKACstRQVFVFSEhrBEMB558Edd8Ddd/di4sR3ePXVsxg9ejSffvop8fHxbpcnIhLWyo8YOvlklwsKM3FxcSxYMIuzzvoD8+b9nZ49N/PDD9Po1Sva7dJEREQkhOXmHgp5fAU93o979hz53JgYaN0a2raFrl1h8GBo08bZ2rY99DEpyf91KxgKgJKSEqy1lJTEBn0rWVl33ukMW3vhhdO58soXefrp8Zx33nnMnj2b+vXru12eiEhYKigooKioiEaNGlFa6ixX366d21WFn9jYWObOfYmJE9vy8sv/oU+frbz//iuMGKHrm4iIiBzuwAHfAU/5fb5a1OvVOxTq9OwJp53mO/BJTIQolyb7UTAUAEVFRZ5b9fC8ARwSjIEZMyArC5599kKuvvoA//vfHznnnHN47733aNiwodslioiEnRxP43jjxo3Ztg0KC+GYY1wuKkxFRUXx0ksP0rFjO+699y+ceeZQ7rnnbf7+d/XuiYiIRIKcnMrbubZs8T2vT1zcoVDnuONg2LDDgx5v+JOYGPxzRdYqGDLGvA709HzaDMi21vb13Ncb+B/QBCgFBlhr82tzvFBVNhgKpRFD4KSbb7wBZ54JM2ZcwTXXxPLUU5dx+umn8+6779K6dWu3SxQRCSv79+8HnGBowwZnX6dOrpUTEe6550aOP74zEyZcwj/+0Z+FC99i1qxTiI11uzIRERGpLmth377KA5+tW8Hza9dh4uMPBTu9ezt/C5cNery3mzcP/sCnqmoVDFlrx3lvG2MeAvZ6bscALwOXWmuXGWMSgSLfrxL+DgVDwb9cvS8NG8JHH8E558BTT13K9dc3ZsaMiznhhBOYM2cOvXv3drtEEZGw4R0x1KhRo4PBUOfO7tUTKS66aDQpKQsZOvQ83nvvNLp3f5i5cyfTvXuY/MYnIiIS4qyF7OyjBz3e27m5Rz6/QYNDoU6/fjBy5JHtXG3aQNOm4RP4VJVfWsmMMQYYC5zu2TUCWG6tXQZgrfUxl3bkKCws9NwKncmny2vcGD78EEaNgmnTzuNPf/qGWbPO5aSTTuJ///sfF198sdslioiEhbKtZIsXO/vUShYYKSnHsX79D5xxxqWkp1/Pccd9xn//+wyTJrVwuzQREZGwZa0zGXNl8/ds3eqs5lVeo0aHgp0BA3y3c7Vt6/xNG2mBT1X5a46hwcB2a+0az+c9AGuM+QRIAl6z1j7gp2OFnLKtZKEaDIEzcmjOHJg4ER59NJVLLlnEL7+M45JLLmHevHlMmzaNBg0auF2miEhIKz9iqFUr0Hz/gdOsWTMWLnyXO++cxv/939+4+uo+vP32S7z++uk0a+Z2dSIiIqHDWme59crm79m2DQoKjnx+kyaHQp1Bg3y3c7VpQ0jN4xusKg2GjDFzAV8TyfzDWvuu5/Z4YGa51z0FGADkAp8bYxZbaz/38fpXAVcBdOzYsXrVh4hQnXzal/r14bXXnOXz/v3vtgwbNo8bbriDRx+9j4ULFzJr1ix69erldpkiIiErOzsbcAKKDRs0v5AboqKiuPvuGxk16lTOOutCPv10GO3aTeaZZ+7jwgtD/EIuIiJSS6WlsHPn0Zdj925FPiaUadbsULjjXZLd1ygfrXUUOJUGQ9ba4Ue73zOf0Bigf5ndm4GvrbU7PY/5EEgFjgiGrLXTgekAaWlptsqVh5BDrWShOcdQeVFRcN990L07XHttDKtW3cvUqUO4995L6d+/P/feey/XXnst9erVc7tUEZGQ4w2GmjdvzoYN0L//UR8udSgtrS8bNy7myiv/wWuvTWP8+Pd46qmneOONs2nZ0u3qRERE/KukxFmRurL5e7Ztg+LiI5/fvPmhYKdnT9/tXG3aaCR0MPJHK9lw4Gdr7eYy+z4BbjHGNAAKgaHAw344VkgKl1ay8i6/HPr2hd//Hm655Sz+8Y8MFiy4nBtvvJGnnnqKqVOnctZZZ2HUyCkiUmV79uwBoGnT5mzc6PyMFfc0bNiQmTMf4ZprxvH731/JV1+NpH37cfz97w/wj390RO+BiNQtY8yzwDnADmttstv1iISi4mLYsaPy+Xu2b3fCofISEw+FOscff2TQ07YttG7trOYlockfwdCFHN5GhrV2jzFmKrAIsMCH1toP/HCskFQ2GAq34XCpqbBkCVx2Gdx1VztOO+1j/ve/OfznP3/l7LPP5uSTT+aOO+7gjDPOUEAkIlIFe/bsITY2luzs+hQVqZUsWAwZMohNm5bw17/+myef/Dd33fUejz9+M889dwvnnBNmF3eR4PI88F/gRZfrEAk6RUVOmFPZKl07djjtX+UlJR0Kdnr39t3O1bo1xMUF/twksGodDFlrL6tg/8s4S9ZHPG8w1KBBbFjOgt6sGcyeDTNmwM03G77//lzuuONMGjV6hgceuI8zzzyTE088kVtuuYVRo0YRE+OvOc9FRMLPnj17aN68Ob/84lwwFAwFj7i4OP77339y002Xcemlf+Pbb+/m3HOf4bjj7uLppycwaJCGD4n4m7X2a2NMJ7frEAmkoiKnXauyVbqyspwJnssyBlq2PBTw9Ovnu52rVSuIjXXn/CT46C/0APDOMdSwYfj+wmgMXHUVjBwJ114Lt94aS8+e1/Doo5ezY8cL/Pvf93HBBRfQsWNHrr32Wq688koSExPdLltEJOhkZ2fTvHlz1njW+eze3d165EidOh3DN9+8xhdfXM9ll/2FlSuv5KST7iM19Q6efvoi+vXTr1ciInKkgoKqBT47dx753KgoJ8xp0wbat3eWZfc1aXOrVqD34aW69F8mALwjhsI5GPJq1w7eeQfefx9uvhnGjIlj2LCrePXVy9m27X0ee+wxpkyZwp133sn555/PhAkTGD58uEYRiYh4eEcMrVrlvJN3zDFuVyQVOf30k9m48XtmzfqAG2+8gyVLJpKaeg99+vyN//znIoYNiw/LkcIiwSYSVjmW4Jaff/Tl2L337dp15HOjo512rTZtnGv+oEG+J21u2dJ5rEhd0F/jARBJwRA4o4dGjYLf/Q6eegruugtOOimGM888n3/963yaNFnBU089xcyZM5k5cyatW7fm4osv5sILL6R///6ai0hEItqePXto1aoVq1dDt276JTDYGWMYO/Yc/t//G8lLL73DLbfcxbJlV3DGGbfSvv113H33NVx6aQu9eytShyJhlWNxR25u5fP3bN0KnnUjDhMT4wQ+bdtC166HlmUvP8qnRQtd68V9+jUlALytZI0aRUYw5FWvHlx/PUycCE8+CQ89BKecAkOGpHDDDY/zwANT+fTTD3nxxReZNm0aDz30EB06dOC8887jvPPOY8iQIRpJJCIRZ8+ePfTq1YulS6FHD7erkaoyxjBhwvlceul5fPjhF9xyy0OsXHkHl19+L5MnX8i4cVdy550nccwxevNDRMRtBw5U3s61dSvs3Xvkc+vVOxTq9OwJp57qe9LmFi2c9i+RUKC/ugPAGww1aRKZ6/c1aQJ/+5sTEs2YAVOnwgUXQPv2cUyadD5PPXU+MTG7mDNnDrNnz2bGjBk89thjNG/enOHDhzN8+HDOOOMMOnfu7PapiIjUuezsbJo2bcbatXDOOW5XI9VljGHkyGGMHDmMFStWcvPN05g79xWef/55nn++F8ceeyU33XQp48e3pH59t6sVCX7GmJnAqUALY8xm4J/W2mfcrUqCVU7O0YMe7+2cnCOfGxd3KNQ5/ngYPtz3pM2JiahNWMKOgqEAyMvLA6Bp08gMhrwaNIAbboDrroM5c+Dxx+H22+HOO+HMMxO55JKJvPrqRKw9wKeffsr777/PZ599xqxZswDo0qULw4cPZ/DgwQwaNIguXbqo7UxEwkpxcTF79uwhOjqRwkKNGAp1KSnH8fHHT7F//3944ok3eOyxZ/jpp5u44oq/cdVVp3PCCWO58cbzOf/8ROpF1qBikSqz1o53uwZxl7Wwb1/l8/ds2eKMBCovPv5QqNO7N5x5pu9Jm5s3V+AjkUvBUADk5+cDCoa8oqNh9GhnW7UKnn8eXnkFLroIGjWC885ryHnnnc+jj55Po0aW1atX89lnn/HZZ5/x2muvMX36dABatmzJoEGDOOmkkzjhhBPo06cPzZs3d/fkRERqISsrC2stRUVtADjuOJcLEr9o1KgRt9xyObfccjnLl6/k3/9+hTlzXmfBgj+yYME11Ks3nLS085kw4XdcdFEHmjRxu2IRkbpnLWRnVz5/z5Yt4Hmf/TANGhwKdfr1c1ZH9hX4NG2qwEekMgqGAiAvzwmGmjfXmPHyevaE++6De+6Bb76Bl1+Gt95yPsbGwrBhhtGje3LuuT257rrrKCkpYeXKlSxYsID58+czf/583n333YOv17FjR/r27XtwS05OpnPnzpqrSERCwrZt2wDYs6c1xjjvbEp46d37OF599R6s/T8WLlzKAw+8zqefvsGCBZNYsACuuSaFTp1Gcu65Z/OHPwykb996+oNGREKKtbB7d+Xz92zd6qzmVV6jRoeCnQEDfLdztW0LjRsr8BHxF/21HADZ2U7E3by5RgxVJCoKhg51tiefhG+/hffeg3ffhauvdh7TpQucfno0w4alMHp0CldddRXgvMO+dOlSMjIyDm5z5syhtLQUgHr16tG1a1d69ux52Na5c2dat25NlGaFE5EgsXXrVs/H1vToAQ0bulyQ1BljDAMHpvL226lY+29WrPiJJ5/8gA8++JANG/7DY4/9m8cea0i9eifRs+dQRowYysSJA0hJidMfQiLiitJSZ7n1qqzS5Zli9TBNmhwKdcovyV42/GncOPDnJhLpFAwFwL593hFDcS5XEhpiYpzZ/U891VnJbOVK+PxzZ3vjDXj6aedx3bvDiSfCwIFJnHjiCG68cQSxsc59ubm5ZGZmsnLlSlatWnVw++ijjw5OBg4QGxtLhw4dOOaYYw7bOnbsSOvWrWnTpg3NmjXTXEYiEhDeEUPr1rXmlFNcLkYCxhhD797H8eSTxwE3s3fvXl59dS5vvfUlixd/RWbmbWRmwtSp8cTEnECHDgPo338Av/vdAEaN6kyLFrpGiUjNlZbCzp2Vt3Nt2wZFRUc+v1mzQ6GOd0l2X6N89GaHSPBSMBQATjAUQ7Nm+ueuLmOcVQGOPx7+9CcoLoalS+GLL+D772HuXKftDJyJ5fr2dVovUlIa0Lv3CYwefQJlpx0qLi5m48aNrFq1ig0bNrBx48aD28cff3zw3fqyYmNjadWqFa1btz5sa9GiBQkJCUdszZo1U+uaiNSINxj67bdW9OvncjHimqZNm3LNNRdwzTUXAJCVtZNZs77h7be/YtmyBWzY8Bi//FLIm28CJFK/fhodOvTl+OOTOemkZM46qxfHHhtPdLSrpyEiLispgayso7dzbdkC27c7v2OX17z5oVDntNN8t3O1aYNWWBQJA/rrNQCcVrL6JCS4XUnoi4lxeo0HDHA+txY2bYKFC52gaPFiePNN8MxPDUD79pCc7Iww6t49hu7du9KzZ1dGjHBer6yCggI2bdrEpk2b2L59O9u2bTts+/XXX/nhhx/YsWMH1toK62zatCkJCQk0b96cJk2a0LhxYxo1akTjxo2PujVo0ID69esTHx9P/fr1D96Oj4/XqCWRCLBhwwYaN04kJ6c+gwe7XY0Ei6SkFlx77flce+35ABQWFvL99yuYPXsR3367iNWrF7F69eesXl3M7Nlw881RGNONpk2Tad/+WLp160ZKSlcGDerGiSe2JiFB1xORUFZcDDt2VD5/z/btTjhUXmLioVDnuOOODHratoXWrZ03XUUkMigYCoDs7P1AAxIT3a4k/BgDHTs62//7f84+a52L4fLlzrZihdOO9t13kJNz6LkxMc68RZ06Oc/v0AE6dIijQ4dudOjQjRNPdFY78KWkpITs7Gx2795d6ZaTk8OmTZvYv38/OTk55OTkkJubW+1zjYuLOywsKvsxNjaW2NhY6tWrd3Ar/3lV9kVHRxMdHU1UVFSlt2vyOGNMjbaoqKgaP7eiTSQYrVu3joYNu1NUBGlpblcjwSo2NpYhQ/ozZEh/wJmIr7CwkOXL1/DZZz+yYEEmP/2UyW+/rSAz810yM0t45x3vsxsQHd2VJk26kpTUmbZt29OpU3t69GhPSkoH+vVrQ5s2MWj6PZHAKypywpzK5u/ZscNp/yovKelQsNO7t+92rtatIU6zW4hIOUEVDGVkZNCsWTO/vJa//vDzx+vs27cfOFbBUIAY41z82raFs846tN9a50K6Zo2zrV7tfNy4EZYtcy7E5SUmQqtWztaypbM5t6Np1SqRli0TSUpyRiM1bUqVh+2XlJQcFhR5b+fm5pKXl0deXh75+flVvl1YWMiBAwcoKio6uBUWFh72edl9pb5+m4hAvoKi8t/zVb2vOo+tzTEkvK1Zs4YDB4YyeDAH50wTqYrY2FjS0o4nLe14YOzB/UVFRaxatZHvvlvLkiXr+Omndfz661qyslaxZs2nrF5d/o2KKKA1cXHtadCgNY0bJ5GQkERSUktat06iffuWdOqURLduSfTokUSrVnHUqxfIMxUJPYWFzvw8la3SlZXl/L5aljHO75/ecCc11Xc7V6tWum6ISM0FVTCUmJjI+PHja/06R2vxceN1liyB774boWDIZcYcCnl8TepaUAC//Qa//uq0p3m3HTuc0GjJEuf23r0VH6NpU2cCvubND3303m7WzFlloWFDaNQomoYNm9KoUVMaNXL2tWzJwdtxcXW7/GZpaanPAKmkpITS0lJKSkqOuH20+6r6OGttjbfS0tJaPf9om1f57/mq3ledx9bmvkcffdT3F1TCwvbt29m0aRPQm9//3u1qJFzUq1eP5ORuJCd3O+I+ay179mSzYsVmli/fzM8/b2b9+s389ttmsrI2s3//RrZsWcSvv2YBPiYgAaAhxjQlOroZsbHNiItrRoMGzWjYsClNmjSjadNmJCQ0o2nTJjRt2pBmzRrSvLmzJSQ0IDGxIUlJDUlKakCDBlFacU1CSkHBoZCnovl7tm51JnYuLyrK+Z20TRtn2oMTTvA9aXPLlih8FZE6F1TBUIcOHXjkkUfcLsPvbr3VmQNHSy8Gt7g4p7WsS5ejPy4/33lHZ/t2JyjKyoI9eyA72/lY9vaaNYc+P3Cg6rVERzsBUXz8oS0u7vCPR7sdG+v8EhETU9HHKOrViyMmJq6Sxx3aoqKcuqKijtyqs1+/9NecgqHwNm/elwA0aHAyY8ce/bEi/mCMISGhOUOHNmfo0JQKH2etZdeubNauzWLNmix++WUHmzZlsXXrDvbsyWbfvmxycrLJzc0mP38HWVmr2bo1m9LSbCoOlHypjzENiYrybnFERztbTEwcMTHxB69dsbFx1KvnfIyLO3JzWqRjDm6xsc5W9nZsbD3i4g59Hh/vfIyLO7R5Hx8dHUV0tCE6OspzDY06eDsmxrnPe7vsPmN03QtF+flHD3q8t3fvPvK50dFO4NO2rTNdwUknVRz4aIJ4EQkWQRUM/fSTs/x4oATqQr1wofNOgH4xCA/x8d75iKr3vKIiJxzav//Qx/K3y39eUOD8cuL96L194IDzy0jZfWVv+1pZIlgYU3GI5P0eKfuLtPd2ZZ+7+djy51fZvpo+T2DJkh+JizvWxz2+/7F8t+JV9A975P6KW/mq/rpVraGgYD1wDLfddiJ+6qoW8QtjDC1aNKdFi+YMHNijys+z1rJ/fx6bN2ezY8c+du48wO7dB9izZz979hxg374D7N17gJycA+zff4ADBw6Qm3uAvLwDFBQcoLi4gKKiAoqLCyguzqagoICSkgJKSw9t1hYABYCPNbSDgsFpz/O1GYw5fJ/zuTlsc36GmDKvV3Zf2c+lNiZPhi+/dAKf7Owj74+JcebnadsWunZ1Rp/7mrS5RQsFPiISeoIqGIqJIWArd/mpS6xKzjwTzj03cMeT4FSv3qGWsrpWUuKEQ96tqKh2H611XrO09MitOvuP9ljvqhnWHvr+9N6u7HM3H1tWVfb583mRKCamPq1a9T5sX0Vtv773l2/X873/aK/rfezhd1f82CNfxvdjo6O7M3Hi35gyRX9RSHgwxtC4cQOOPbYBxx7btk6PVVpaSm5uIfv2FZCbW0B+fjH5+cUUFDgfCwu9t4soKDj0ufdjUZFzu7Dw0O2iomKKi4spLi6itNRpaS4ttZSUlHpul3rapUsPtjuX/fzQZrG2/L7Dn+O931rv5m11PtTyfHgLdPmWaMvGjWvq9N843DVoAD17wqmnHhn4tGnjBD6alF1EwlVQBUPdu8NHH7ldhUjoi452Nq06EV70hjD07t2F9PTX3S5DRIJMVFQUjRrF06hR5K6vbcxrbpcQ0h580O0KRETco9xbRERqxBjzujEmw7NtMMZklLmvtzFmgTHmR2PMCmNMhX+tGWP+aoyxxpgWASlcREREREQOCqoRQyIiEjqsteO8t40xDwF7PbdjgJeBS621y4wxiVQwAYgxpgMwAvi17isWEREREZHyNGJIRERqxTizno4FZnp2jQCWW2uXAVhrd1lrSyp4+sPALVQ8SY+IiIiIiNQhBUMiIlJbg4Ht1lrvzKc9AGuM+cQYs8QYc4uvJxljRgO/eQMkEREREREJPLWSiYhIhYwxc4HWPu76h7X2Xc/t8RwaLQTOteUUYACQC3xujFlsrf28zOs2AP6OM7qoshquAq4C6NixY01OQ0REREREKqBgSEREKmStHX60+z3zCY0B+pfZvRn42lq70/OYD4FU4PMyj+kKdAaWOZ1otAeWGGNOsNZuK1fDdGA6QFpamlrORERERET8SK1kIiJSG8OBn621m8vs+wRIMcY08ARHQ4GVZZ9krV1hrW1pre1kre2EEyallg+FRERERESkbikYEhGR2riQw9vIsNbuAaYCi4AMYIm19gMAY8zTxpi0QBcpIiIiIiK+qZVMRERqzFp7WQX7X8ZZsr78/isreHwnvxYmIiIiIiJVohFDIiIiIiIiIiIRSsGQiIiIiIiIiEiEUjAkIiIiIiIiIhKhFAyJiIiIiIiIiEQoY611u4aDjDE5wCq366gjLYCdbhdRR8L13ML1vEDnFqp6Wmsbu12Em3SdCFnhem7hel6gcwtFEX+NADDGZAEba/j0cP2/UZFIOt9IOlfQ+Yaz2pzrMdbaJF93BNuqZKustWG5jLExJl3nFlrC9bxA5xaqjDHpbtcQBHSdCEHhem7hel6gcwtFukY4KvqjpyrC9f9GRSLpfCPpXEHnG87q6lzVSiYiIiIiIiIiEqEUDImIiIiIiIiIRKhgC4amu11AHdK5hZ5wPS/QuYWqcD63qgrnfwOdW+gJ1/MCnVsoCtfzCqRI+zeMpPONpHMFnW84q5NzDarJp0VEREREREREJHCCbcSQiIiIiIiIiIgEiOvBkDEm2hiz1Bgzx/N5Z2PMQmPMWmPM68aYWLdrrAljzAZjzApjTIZ3lQhjTIIx5jNjzBrPx+Zu11ldxphmxpg3jTE/G2N+MsYMCpPz6un5Wnm3fcaYG8Ph3ACMMX82xvxojMk0xsw0xsSHw/eaMeYGzzn9aIy50bMvJL9mxphnjTE7jDGZZfb5PBfjmOb52i03xqS6V3nd0jUi9Og6EZLnFpbXCNB1IhKuEzVhjDnLGLPK8+8zxcf9cZ7/92s93wedXCjTL6pwrn8xxqz0/D/53BhzjBt1+ktl51vmcRcYY6wxJqRXsqrK+Rpjxnq+xj8aY14NdI3+UoX/yx2NMfM8vzcuN8ac7Uad/uLrZ365+/36c971YAi4AfipzOf3Aw9ba7sBe4ArXKnKP06z1vYts5zcFOBza2134HPP56HmUeBja20voA/O1y7kz8tau8rzteoL9AdygdmEwbkZY9oBfwLSrLXJQDRwISH+vWaMSQb+CJyA83/xHGNMN0L3a/Y8cFa5fRWdy++A7p7tKuDJANXoBl0jQo+uEyEkXK8RoOsEkXOdqBZjTDTwOM6/0XHAeGPMceUedgWwx/P//2Gc74eQU8VzXYrz/d8beBN4ILBV+k8VzxdjTGOc3y8WBrZC/6rK+RpjugO3Aidba48Hbgx0nf5Qxa/tbcAb1tp+ONexJwJbpd89z5E/88vy6895V4MhY0x7YCTwtOdzA5yO80MJ4AXgPFeKqxujcc4JQvDcjDFNgSHAMwDW2kJrbTYhfl4+DAPWWWs3Ej7nFgPUN8bEAA2ArYT+99qxwEJrba61thj4ChhDiH7NrLVfA7vL7a7oXEYDL1rH90AzY0ybgBQaQLpGhN656ToRsucWjtcI0HUi7K8TNXQCsNZau95aWwi8hvPvVVbZf9c3gWGea1CoqfRcrbXzrLW5nk+/B9oHuEZ/qsrXFuBfOGFffiCLqwNVOd8/Ao9ba/cAWGt3BLhGf6nKuVqgied2U2BLAOvzuwp+5pfl15/zbo8YegS4BSj1fJ4IZHsu3gCbgXYu1OUPFvjUGLPYGHOVZ18ra+1Wz+1tQCt3SquxzkAW8JxniN7TxpiGhP55lXchMNNzO+TPzVr7G/Af4FecX/b3AosJ/e+1TGCwMSbRGNMAOBvoQBh8zcqo6FzaAZvKPC4Uv35V8Qi6RoQaXSdC7NzC+BoBuk5EwnWiJqryb3PwMZ7vg70416BQU93/B1cAH9VpRXWr0vP1tNt0sNZ+EMjC6khVvr49gB7GmO+MMd8bY442AiWYVeVc7wQuMcZsBj4Erg9Maa7x689514IhY8w5wA5r7WK3aqhjp1hrU3GGeE02xgwpe6d1loMLtSXhYoBU4EnPEL0DlBt+HaLndZBx5lAYBcwqf1+onptnvoHROH+wtQUacvRhiSHBWvsTzrs9nwIfAxlASbnHhOTXzJdwOpeq0DUiZL/euk6E2LmF6zUCdJ0QqQ5jzCVAGvCg27XUFWNMFDAV+KvbtQRQDE6r0anAeGCGMaaZmwXVofHA89ba9jhvBLzk+ZpLFbj5D3UyMMoYswFnKNjpOPMSNPMMZQZnKONv7pRXO5534LzD9WbjDH/b7h3e5fkYakP5NgObrbXeftw3cf4ACPXzKut3wBJr7XbP5+FwbsOBX6y1WdbaIuBtnO+/kP9es9Y+Y63tb60dgjMHxmrC42vmVdG5/IbzrrdXSH79KqFrRGj+39V1IvTOLWyvEaDrRJnHhezXsA5U5d/m4GM83wdNgV0Bqc6/qvT/wBgzHPgHMMpaWxCg2upCZefbGEgGvvT8fjEQeC+EJ6Cuytd3M/CetbbIWvsLzs/A7gGqz5+qcq5XAG8AWGsXAPFAi4BU5w6//px3LRiy1t5qrW1vre2EMyT7C2vtxcA84Peeh00E3nWpxBozxjT0TGqGZwj9CJzhzO/hnBOE4LlZa7cBm4wxPT27hgErCfHzKmc8h9oDIDzO7VdgoDGmgac/3vt1C4fvtZaejx1x5o14lfD4mnlVdC7vARM8qxEMBPaWaSUIC7pGhOa56ToRkucWttcI0HUinK8TtbAI6G6clfdica4x75V7TNl/19/jXINCcTRWpedqjOkH/A8nFArlkBQqOV9r7V5rbQtrbSfP7xff45x3ujvl1lpV/i+/gzNaCGNMC5zWsvUBrNFfqnKuv+JcwzDGHIsTDGUFtMrA8u/PeWut6xvOf9Y5nttdgB+AtTjDtOPcrq8G59MFWObZfgT+4dmfiLNixBpgLpDgdq01OLe+QDqwHOcHTfNwOC/PuTXEeTeoaZl94XJudwE/4/zx+RIQFybfa9/g/AGzDBgWyl8znD80twJFOO/uXFHRuQAGZ2WGdcAKnNVEXD+HOvy30TUihDZdJ0Lv3ML1GuE5N10nIuA6UYN/y7NxRk6sK/Mz+G6ckACcPyhnef7//wB0cbvmOjzXucB2nFbLDJzRJa7XXVfnW+6xX4b690YVvr4Gp31upednwYVu11yH53oc8J3n530GMMLtmmt5vr5+5l8NXF3ma+u3n/PG86IiIiIiIiIiIhJhNBmTiIiIiIiIiEiEUjAkIiIiIiIiIhKhFAyJiIiIiIiIiEQoBUMiIiIiIiIiIhFKwZCIiIiIiIiISIRSMCQiIiIiIiIiEqEUDImIiIiIiIiIRCgFQyIiIiIiIiIiEer/A9eDUdYu4I4oAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1440x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize=(20,6))\n",
    "plt.subplot(1,3,1)\n",
    "plt.plot(tvEpsp,vEpspDend[0].T,c='b')\n",
    "for n in range(len(vEpspDend)-1):\n",
    "    plt.plot(tvEpsp,vEpspDend[n+1].T,c='k')\n",
    "plt.xlim(40,100)\n",
    "plt.subplot(1,3,2)\n",
    "plt.plot(tvEpsp,vEpspSoma[0].T,c='b')\n",
    "for n in range(len(vEpspSoma)-1):\n",
    "    plt.plot(tvEpsp,vEpspSoma[n+1].T,c='k')\n",
    "plt.xlim(40,100)\n",
    "plt.subplot(1,3,3)\n",
    "plt.plot(range(len(vEpspDend)),cEpspAmpDend,c='k')\n",
    "plt.plot(range(len(vEpspDend)),cEpspAmpSoma,c='b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "neuron",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.21"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}