05_MCMC_fitting_autonomous.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'emcee'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-1-f1cd78d02e76>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0;32mimport\u001b[0m \u001b[0memcee\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      7\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mnumpy\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mpandas\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'emcee'"
     ]
    }
   ],
   "source": [
    "from __future__ import print_function, division\n",
    "import os\n",
    "import sys\n",
    "import time \n",
    "\n",
    "import emcee\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import math\n",
    "from scipy.integrate import odeint\n",
    "import matplotlib.pyplot as plt\n",
    "from multiprocessing import Pool\n",
    "import pickle\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "from toggle_equations import *\n",
    "from attractor_stats import *\n",
    "from steady_states import *\n",
    "from sim_funcs import *\n",
    "from plotting_funcs import *\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Importing data to initialise run"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`mcmc_initialise.csv`: Excel file that contains the normalised values of Sox2, Tbxta and TCF_mRNA of an 18ss TCF-GFP tailbud. \n",
    "\n",
    "Filename of sample: TCF-GFP_18ss_B1tbxta594_B4sox2647_B2egfp546_4_Series009.ims\n",
    "\n",
    "NMp number = 182"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### df = pd.read_csv(\"normalised_exp.csv\")\n",
    "df = pd.read_csv(\"files/csv/mcmc_initialise.csv\")\n",
    "data = np.load(\"files/npy/loess_fit_vals_050421.npy\")\n",
    "#filt = (df[\"stage\"] == \"18ss\") & (df[\"ID\"] == 7)\n",
    "nmp_numbers = len(df)\n",
    "\"\"\"\n",
    "groups = df.loc[filt, ].groupby(\"ID\")\n",
    "for _ , group in groups:\n",
    "    nmp_list = group[[\"sox2_norm\", \"tbxta_norm\", \"tcf_norm\"]].values.tolist()\n",
    "\"\"\"\n",
    "\n",
    "nmp_list = df[['sox2_norm', 'tbxta_norm']].values.tolist()\n",
    "nmp_list_tcf = df[['sox2_norm', 'tbxta_norm','tcf_norm', 'TCF_group']].values.tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "60"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mylist = [l for l in nmp_list_tcf if l[3] == 0]\n",
    "len(mylist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0.360546782606608, 0.33999913174345897, 0.49255643564356394, 1.0],\n",
       " [0.136614642376552, 0.142272516786293, 0.44234851485148496, 1.0],\n",
       " [0.0022842615479472398, 0.289887705487381, 0.361341584158416, 1.0],\n",
       " [0.0167444172357735, 0.13319923593424401, 0.305214356435644, 1.0],\n",
       " [0.0892862232725647, 0.18254514934012497, 0.29348861386138597, 1.0],\n",
       " [0.0484215442668186, 0.0673694721000232, 0.420704455445545, 1.0],\n",
       " [0.20611309994994398, 0.21186183144246398, 0.534039603960396, 1.0],\n",
       " [0.0201958124205222, 0.127264702477425, 0.585712871287129, 1.0],\n",
       " [0.0282417336785062, 0.20639326232924302, 0.551816831683168, 1.0],\n",
       " [0.0287867960881521, 0.21359689742996105, 0.42941386138613896, 1.0],\n",
       " [0.0152667480426509, 0.28986744616809396, 0.40976435643564396, 1.0],\n",
       " [0.0965945600771289, 0.134669483676777, 0.286610891089109, 0.0],\n",
       " [0.0953894220888292, 0.0411481245658717, 0.237076237623762, 0.0],\n",
       " [0.121361895937085, 0.094022053716138, 0.358061881188119, 0.0],\n",
       " [0.13895341016546398, 0.347962491317435, 0.24055792079207897, 0.0],\n",
       " [0.0569000150517235, 0.207844697846724, 0.7270841584158421, 0.0],\n",
       " [0.168928842352449, 0.00085812688122251, 0.43578514851485206, 0.0],\n",
       " [0.153186039801557, 0.0890686501504978, 0.285781683168317, 1.0],\n",
       " [0.0810042749894863, 0.0242243574901598, 0.35708861386138596, 1.0],\n",
       " [0.0628281938281934, 0.197752662653392, 0.454590099009901, 1.0],\n",
       " [0.0911044314574019, 0.13492417226209802, 0.6809306930693071, 0.0],\n",
       " [0.0462447950290308, 0.0210783746237555, 0.32227871287128695, 0.0],\n",
       " [0.0209844027141108, 0.0427037508682565, 0.47135693069306894, 0.0],\n",
       " [0.08258995655002499, 0.0441161148413985, 0.21316534653465302, 0.0],\n",
       " [0.529096081501332, 0.0205588677934707, 0.19806930693069297, 0.0],\n",
       " [0.187726494683641, 0.0, 0.24159158415841603, 0.0],\n",
       " [0.23349723543345702, 8.24843713822778e-05, 0.30019702970296996, 0.0],\n",
       " [0.0230676412449226, 0.11278363047001598, 0.25424504950495, 1.0],\n",
       " [0.0289118104022911, 0.0882973489233619, 0.5147970297029699, 1.0],\n",
       " [0.522690348044851, 0.20315755962028198, 0.236316831683168, 1.0],\n",
       " [0.129595838723534, 0.24266323222968303, 0.42026485148514897, 1.0],\n",
       " [0.0773343547836228, 0.0677182218106043, 0.6207524752475251, 1.0],\n",
       " [0.17990609924836398, 1.0, 0.4455445544554461, 1.0],\n",
       " [0.00948008546978631, 0.34608561009492894, 0.40515198019802007, 1.0],\n",
       " [0.0190221780393855, 0.18828287798101398, 0.5740346534653471, 1.0],\n",
       " [0.191699949644234, 0.202646735355406, 0.27535495049505, 1.0],\n",
       " [0.0132550176995266, 0.19628820328779803, 0.34396138613861404, 1.0],\n",
       " [0.127533602597497, 0.287135332252836, 0.180358910891089, 1.0],\n",
       " [0.20835635680285397, 0.22255441074322801, 0.28778663366336604, 1.0],\n",
       " [0.0712211548222272, 0.19912450798795997, 0.53220297029703, 1.0],\n",
       " [0.0366441957604146, 0.21260997916184302, 0.465258910891089, 1.0],\n",
       " [0.0381303659268986, 0.20834683954619101, 0.31549504950495, 1.0],\n",
       " [0.133806320823734, 0.29976846492243603, 0.48537524752475203, 1.0],\n",
       " [0.0011861358125505401, 0.304522169483677, 0.38694356435643606, 1.0],\n",
       " [0.0139981027827686, 0.15507785366983098, 0.31501782178217796, 1.0],\n",
       " [0.14642976620823098, 0.29600601991201697, 0.525316831683168, 1.0],\n",
       " [0.0932751800081109, 0.141984545033573, 0.41048514851485207, 1.0],\n",
       " [0.00811042864407974, 0.211805394767307, 0.429431188118812, 0.0],\n",
       " [0.0, 0.0709539245195647, 0.507955445544555, 0.0],\n",
       " [0.0660035574073231, 0.123804700162075, 0.53110396039604, 1.0],\n",
       " [0.0422503376636625, 0.31314685112294505, 0.55640099009901, 0.0],\n",
       " [0.0146241744679766, 0.24687427645288298, 0.284507425742574, 1.0],\n",
       " [0.0198202694208487, 0.103837693910627, 0.310052475247525, 1.0],\n",
       " [0.3445349492516889, 0.148483445241954, 0.21202821782178197, 0.0],\n",
       " [0.20357280908664002, 0.174218569113221, 0.214376237623762, 0.0],\n",
       " [0.0731978811573925, 0.0986585436443621, 0.410580198019802, 0.0],\n",
       " [0.0318926517086206, 0.09508422088446401, 0.383867821782178, 0.0],\n",
       " [0.0631417297280539, 0.41568650150497793, 0.242291584158416, 0.0],\n",
       " [0.0468823680311396, 0.208821486455198, 0.197910891089109, 0.0],\n",
       " [0.171115592735368, 0.211149861078953, 0.272534653465347, 0.0],\n",
       " [0.321847351521749, 0.05249334336651999, 0.23384207920792102, 0.0],\n",
       " [0.0654209907034356, 0.163586767770317, 0.298889108910891, 0.0],\n",
       " [0.116909886181968, 0.22131714517249398, 0.28299603960396, 1.0],\n",
       " [0.0720842536470426, 0.228245832368604, 0.272206435643564, 1.0],\n",
       " [0.0279857043631496, 0.1979943273906, 0.5160841584158421, 1.0],\n",
       " [0.18894863461866399, 0.154594524195416, 0.20582574257425698, 1.0],\n",
       " [0.168264266258487, 0.26749536929844897, 0.312978712871287, 1.0],\n",
       " [0.0386859295389322, 0.25264962954387604, 0.464595049504951, 0.0],\n",
       " [0.0657700306685115, 0.160466832600139, 0.338978712871287, 1.0],\n",
       " [0.0177655341536606, 0.24931263023848102, 0.31520693069306893, 1.0],\n",
       " [0.0427528952065012, 0.37890426024542695, 0.514613861386139, 1.0],\n",
       " [0.0464298162139565, 0.116703808752026, 0.265221782178218, 1.0],\n",
       " [0.0170639538227127, 0.0300170757119704, 0.58360396039604, 0.0],\n",
       " [0.0182940946738402, 0.0616592382495948, 0.268670792079208, 0.0],\n",
       " [0.197272587711293, 0.335864783514703, 0.231080693069307, 0.0],\n",
       " [0.24262828093816696, 0.0801733618893262, 0.34787722772277196, 0.0],\n",
       " [0.073832453815962, 0.0671639847186848, 0.305692574257426, 0.0],\n",
       " [0.0494741647918687, 0.315006367214633, 0.44578019801980195, 0.0],\n",
       " [0.0788565290725788, 0.013374044917804999, 0.320264851485149, 0.0],\n",
       " [0.081261304419356, 0.172681755035888, 0.222585643564356, 0.0],\n",
       " [0.0156377905270154, 0.1813498495022, 0.401507920792079, 0.0],\n",
       " [0.0308205289505648, 0.0929454156054642, 0.289335643564356, 0.0],\n",
       " [0.15336706052843, 0.12757293354943303, 0.273534653465347, 0.0],\n",
       " [0.344264918333149, 0.0753979509145636, 0.37492673267326704, 0.0],\n",
       " [0.0658100352490361, 0.113019506830285, 0.139182178217822, 0.0],\n",
       " [0.0644753824312884, 0.195583468395462, 0.437741584158416, 0.0],\n",
       " [0.0447081190796347, 0.12603901366057002, 0.238646534653465, 0.0],\n",
       " [0.0321381798215896, 0.0704127112757583, 0.22448168316831696, 0.0],\n",
       " [0.0415767605390818, 0.0452983908312109, 0.157473762376238, 0.0],\n",
       " [0.182807431450901, 0.0440596781662422, 0.181358910891089, 0.0],\n",
       " [0.272081653349309, 0.0463171451724936, 0.282897524752475, 0.0],\n",
       " [0.314806545349443, 0.00175821949525349, 0.16219950495049498, 0.0],\n",
       " [0.0523219908679544, 0.0620325885621672, 0.25003762376237604, 1.0],\n",
       " [0.0846721949663237, 0.0380426603380412, 0.273808415841584, 1.0],\n",
       " [0.0570290298239149, 0.278080863625839, 0.29272524752475204, 1.0],\n",
       " [0.0936002172248722, 0.33466224820560303, 0.213514851485149, 1.0],\n",
       " [0.16820225915867398, 0.20774050706181998, 0.35803861386138597, 1.0],\n",
       " [0.190642828603875, 0.273338735818476, 0.287323267326733, 1.0],\n",
       " [0.22683547266162002, 0.33946370687659194, 0.360224257425743, 1.0],\n",
       " [0.0914179673572624, 0.251732171799028, 0.337907425742574, 1.0],\n",
       " [0.0362341488100388, 0.418637126649687, 0.43169752475247497, 1.0],\n",
       " [0.0471313965449044, 0.19634029868025002, 0.434568811881188, 1.0],\n",
       " [0.00427148908550029, 0.10673477656864999, 0.31471287128712894, 0.0],\n",
       " [0.0372687672738529, 0.186103554063441, 0.22536831683168304, 0.0],\n",
       " [0.0543252202377172, 0.00586652002778421, 0.151782178217822, 0.0],\n",
       " [0.10327432491020198, 0.0436863278536698, 0.25810099009901, 0.0],\n",
       " [0.103615363959173, 0.13785453808752, 0.20799752475247502, 1.0],\n",
       " [0.10113708019568199, 0.127386258393147, 0.203863861386139, 1.0],\n",
       " [0.0316466235383951, 0.232805626302385, 0.17176386138613903, 1.0],\n",
       " [0.122808561580301, 0.167569171104422, 0.30507821782178196, 1.0],\n",
       " [0.127592609353771, 0.267284093540171, 0.259610891089109, 1.0],\n",
       " [0.13345078011432301, 0.302229972215791, 0.158433168316832, 1.0],\n",
       " [0.10518954420281099, 0.34762676545496596, 0.252309405940594, 1.0],\n",
       " [0.15296651466592903, 0.264971636952998, 0.38410396039604, 1.0],\n",
       " [0.0415657592794375, 0.12668441768928002, 0.37310544554455394, 1.0],\n",
       " [0.0633852576119966, 0.18742040981708702, 0.552707920792079, 1.0],\n",
       " [0.117988509684359, 0.312123755498958, 0.40402376237623794, 1.0],\n",
       " [0.174240950588842, 0.37209423477656894, 0.46784158415841604, 1.0],\n",
       " [0.103098304755895, 0.0585045728177819, 0.29395544554455394, 1.0],\n",
       " [0.13500045755239, 0.37720826580226896, 0.198109900990099, 1.0],\n",
       " [0.0388819519835021, 0.221353322528363, 0.67330198019802, 1.0],\n",
       " [0.00843546586084106, 0.233566797869877, 0.376307425742574, 1.0],\n",
       " [0.13598006971798302, 0.181138573743922, 0.260884158415842, 1.0],\n",
       " [0.00146016718914317, 0.2093916415837, 0.4099366336633661, 0.0],\n",
       " [0.00124714279785035, 0.322370629775411, 0.273260396039604, 1.0],\n",
       " [0.0207193723681362, 0.25123437138226395, 0.3820272277227721, 0.0],\n",
       " [0.0316886283479458, 0.0582918499652697, 0.149335643564356, 1.0],\n",
       " [0.0121808947124446, 0.124026105579995, 0.202216831683168, 1.0],\n",
       " [0.0932431763436914, 0.0970435864783515, 0.14474158415841598, 1.0],\n",
       " [0.23936290705285798, 0.156843308636258, 0.21065742574257396, 1.0],\n",
       " [0.207616272063151, 0.26693968511229393, 0.17805445544554502, 1.0],\n",
       " [0.10652869753586801, 0.139159817087289, 0.15239603960396, 1.0],\n",
       " [0.21778743666149802, 0.338394304237092, 0.20623069306930703, 1.0],\n",
       " [0.263930720067448, 0.30255412132438103, 0.18666336633663397, 1.0],\n",
       " [0.14593170918070097, 0.144022053716138, 0.191820297029703, 1.0],\n",
       " [0.0467368513694818, 0.25780417920814996, 0.4157851485148521, 1.0],\n",
       " [0.276837197859155, 0.11320184070386699, 0.21806138613861398, 1.0],\n",
       " [0.0321181775313273, 0.32970884463996303, 0.44969207920792104, 1.0],\n",
       " [0.11489115503725199, 0.12375405186385699, 0.5112227722772279, 1.0],\n",
       " [0.11849156728445401, 0.33613104885390105, 0.44650742574257396, 1.0],\n",
       " [0.0781354465086252, 0.22536900902986803, 0.290509405940594, 1.0],\n",
       " [0.0951358930597553, 0.23601094003241502, 0.30204455445544603, 1.0],\n",
       " [0.0250358666067265, 0.218923651308173, 0.208515841584158, 1.0],\n",
       " [0.206442637682015, 0.168645809215096, 0.29371732673267303, 1.0],\n",
       " [0.0389874640646354, 0.0998581847649919, 0.378040594059406, 1.0],\n",
       " [0.0338278732914919, 0.0174143320213012, 0.15147772277227697, 0.0],\n",
       " [0.218965571557943, 0.19828664042602498, 0.33343366336633695, 1.0],\n",
       " [0.0522724851995553, 0.15922522574670098, 0.193209405940594, 0.0],\n",
       " [0.0716051987952621, 0.11854451261866199, 0.249275247524752, 0.0],\n",
       " [0.0932076722784759, 0.0726412363973142, 0.23280693069306896, 0.0],\n",
       " [0.13573504166227, 0.124086883537856, 0.196831188118812, 1.0],\n",
       " [0.16278313866937802, 0.23477222736744602, 0.131705445544554, 1.0],\n",
       " [0.41565309227906605, 0.20867243574901603, 0.12108168316831698, 1.0],\n",
       " [0.112057830621606, 0.148073917573512, 0.350274752475248, 1.0],\n",
       " [0.18310396540403898, 0.243142220421394, 0.17122673267326702, 1.0],\n",
       " [0.0367122035473062, 0.122908948830748, 0.17647277227722802, 1.0],\n",
       " [0.00683528263986225, 0.2986990622829361, 0.26623415841584197, 1.0],\n",
       " [0.0451156657437277, 0.119586420467701, 0.257719801980198, 1.0],\n",
       " [0.10382488794967, 0.250924693216022, 0.260478712871287, 1.0],\n",
       " [0.098005221597873, 0.3365203172030561, 0.315158415841584, 1.0],\n",
       " [0.17862745284335102, 0.0915098981245658, 0.215881188118812, 1.0],\n",
       " [0.10240872579910401, 0.17352396388052801, 0.258798514851485, 1.0],\n",
       " [0.0232381607694081, 0.139591051169252, 0.197383663366337, 1.0],\n",
       " [0.381239151882891, 0.119066913637416, 0.25529851485148497, 1.0],\n",
       " [0.33840924785888, 0.204212491317435, 0.18935495049505, 1.0],\n",
       " [0.0889666866856255, 0.0279593077101181, 0.215853465346535, 0.0],\n",
       " [0.21442755195469898, 0.0213938411669368, 0.165903465346535, 0.0],\n",
       " [0.0393890100416498, 0.131549548506599, 0.158992079207921, 0.0],\n",
       " [0.13663814506761, 0.0405302153276221, 0.205662871287129, 0.0],\n",
       " [0.10082604458210499, 0.183879370224589, 0.127921287128713, 1.0],\n",
       " [0.108502423527494, 0.218149455892568, 0.169359405940594, 1.0],\n",
       " [0.219698655496054, 0.139336362583931, 0.133959900990099, 1.0],\n",
       " [0.00777589033944387, 0.0323946515397082, 0.173274752475248, 1.0],\n",
       " [0.0853932775302772, 0.07505788376939099, 0.17333217821782199, 1.0],\n",
       " [0.102404225283795, 0.448995716601065, 0.19809306930693105, 1.0],\n",
       " [0.0649044315574134, 0.254993922204214, 0.31421287128712894, 1.0],\n",
       " [0.17537858084750696, 0.0153826117156749, 0.182453465346535, 1.0],\n",
       " [0.15332655589064997, 0.0548231650845103, 0.19545742574257396, 1.0],\n",
       " [0.10160063327251, 0.0628617735586941, 0.17707772277227696, 1.0],\n",
       " [0.0376423100445001, 0.17813295901829101, 0.123503465346535, 1.0],\n",
       " [0.0530400730883687, 0.0558375781430887, 0.187141089108911, 1.0],\n",
       " [0.0239652440204404, 0.0875202593192868, 0.23245, 1.0]]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nmp_list_tcf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([26., 40., 37., 24., 18., 16., 11.,  6.,  1.,  3.]),\n",
       " array([0.12108168, 0.18168193, 0.24228218, 0.30288243, 0.36348267,\n",
       "        0.42408292, 0.48468317, 0.54528342, 0.60588366, 0.66648391,\n",
       "        0.72708416]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAPbUlEQVR4nO3dfYxldX3H8ffHBapVDCij2QB2qFkfqNGlHSkJiVHUBqFVTGkDaQ0mtKuNthpN66pNqn1I17ZKm9SYLELZJj4GNVBQW0SIoSnYQZeHdasobluEsOMDUdqUFvz2jznbjrMze8/cx/mt71dyM+f87u/O+f44dz57OPf87klVIUlqz+NmXYAkaTgGuCQ1ygCXpEYZ4JLUKANckhp1zDQ3dtJJJ9X8/Pw0NylJzbv99tu/XVVzq9unGuDz8/MsLi5Oc5OS1Lwk/7pWu6dQJKlRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqN6B3iSLUm+nOS6bv20JLcluSfJx5IcN7kyJUmrbeQI/E3A/hXr7wEuq6ptwPeAS8dZmCTpyHoFeJJTgPOBD3brAc4Bru667AEumESBkqS19Z2J+ZfA7wHHd+tPBR6qqke79fuAk9d6YZIdwA6AZzzjGcNXOkPzO6+fyXYP7Dp/JtuV1IaBR+BJfhE4WFW3r2xeo+uat/apqt1VtVBVC3Nzh03llyQNqc8R+NnAK5OcBzweeDLLR+QnJDmmOwo/Bbh/cmVKklYbeAReVW+vqlOqah64CPh8Vf0acBNwYdftEuCaiVUpSTrMKNeBvw14S5Kvs3xO/IrxlCRJ6mNDXydbVTcDN3fL9wJnjr8kSVIfzsSUpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGbWgmpqbLr7GVdCQegUtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RG9bmp8eOTfDHJHUn2JXl3135Vkm8m2ds9tk++XEnSIX2uA38EOKeqHk5yLHBLks90z/1uVV09ufIkSesZGOBVVcDD3eqx3aMmWZQkabBe58CTbEmyFzgI3FBVt3VP/UmSO5NcluQnJlalJOkwvQK8qh6rqu3AKcCZSZ4HvB14DvBC4Cks36X+MEl2JFlMsri0tDSmsiVJG7oKpaoeYvmu9OdW1QO17BHgb1jnDvVVtbuqFqpqYW5ubuSCJUnL+lyFMpfkhG75CcDLgH9JsrVrC3ABcPckC5Uk/ag+V6FsBfYk2cJy4H+8qq5L8vkkc0CAvcDrJ1inJGmVPleh3AmcsUb7OROpSJLUizMxJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVF97on5+CRfTHJHkn1J3t21n5bktiT3JPlYkuMmX64k6ZA+R+CPAOdU1QuA7cC5Sc4C3gNcVlXbgO8Bl06uTEnSagMDvJY93K0e2z0KOAe4umvfw/Kd6SVJU9LrHHiSLUn2AgeBG4BvAA9V1aNdl/uAk9d57Y4ki0kWl5aWxlGzJImeAV5Vj1XVduAU4EzguWt1W+e1u6tqoaoW5ubmhq9UkvQjNnQVSlU9BNwMnAWckOSY7qlTgPvHW5ok6Uj6XIUyl+SEbvkJwMuA/cBNwIVdt0uAayZVpCTpcMcM7sJWYE+SLSwH/ser6rokXwE+muSPgS8DV0ywTknSKgMDvKruBM5Yo/1els+HS5JmwJmYktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNarPDR02hfmd18+6BEnaVDwCl6RG9bkn5qlJbkqyP8m+JG/q2t+V5FtJ9naP8yZfriTpkD6nUB4F3lpVX0pyPHB7khu65y6rqr+YXHmSpPX0uSfmA8AD3fIPkuwHTp50YZKkI9vQOfAk8yzf4Pi2rumNSe5McmWSE9d5zY4ki0kWl5aWRipWkvT/egd4kicBnwDeXFXfBz4APBPYzvIR+nvXel1V7a6qhapamJubG0PJkiToGeBJjmU5vD9UVZ8EqKoHq+qxqvohcDlw5uTKlCSt1ucqlABXAPur6n0r2reu6PZq4O7xlydJWk+fq1DOBl4D3JVkb9f2DuDiJNuBAg4Ar5tIhZKkNfW5CuUWIGs89enxlyNJ6quZqfSanll+bcGBXefPbNtSa5xKL0mNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqVJ97Yp6a5KYk+5PsS/Kmrv0pSW5Ick/388TJlytJOqTPEfijwFur6rnAWcAbkpwO7ARurKptwI3duiRpSgYGeFU9UFVf6pZ/AOwHTgZeBezpuu0BLphUkZKkw23oHHiSeeAM4Dbg6VX1ACyHPPC0dV6zI8liksWlpaXRqpUk/Z/eAZ7kScAngDdX1ff7vq6qdlfVQlUtzM3NDVOjJGkNvQI8ybEsh/eHquqTXfODSbZ2z28FDk6mREnSWvpchRLgCmB/Vb1vxVPXApd0y5cA14y/PEnSeo7p0eds4DXAXUn2dm3vAHYBH09yKfBvwK9MpkRJ0loGBnhV3QJknadfOt5yJEl9ORNTkhplgEtSowxwSWqUAS5JjTLAJalRfS4jlKZmfuf1M9nugV3nz2S70ig8ApekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEY5E1NidjNAwVmgGp5H4JLUqD73xLwyycEkd69oe1eSbyXZ2z3Om2yZkqTV+hyBXwWcu0b7ZVW1vXt8erxlSZIGGRjgVfUF4LtTqEWStAGjnAN/Y5I7u1MsJ67XKcmOJItJFpeWlkbYnCRppWED/APAM4HtwAPAe9frWFW7q2qhqhbm5uaG3JwkabWhAryqHqyqx6rqh8DlwJnjLUuSNMhQAZ5k64rVVwN3r9dXkjQZAyfyJPkI8GLgpCT3AX8AvDjJdqCAA8DrJlijJGkNAwO8qi5eo/mKCdQiSdoAZ2JKUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUqIHfBy5psuZ3Xj+T7R7Ydf5Mtqvx8Qhckho1MMCTXJnkYJK7V7Q9JckNSe7pfp442TIlSav1OQK/Cjh3VdtO4Maq2gbc2K1LkqZoYIBX1ReA765qfhWwp1veA1ww5rokSQMMew786VX1AED382nrdUyyI8liksWlpaUhNydJWm3iH2JW1e6qWqiqhbm5uUlvTpJ+bAwb4A8m2QrQ/Tw4vpIkSX0MG+DXApd0y5cA14ynHElSX30uI/wI8E/As5Pcl+RSYBfw8iT3AC/v1iVJUzRwJmZVXbzOUy8dcy2SpA1wJqYkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNGvh94JKOTvM7r5/Ztg/sOn9m2z6aeAQuSY0a6Qg8yQHgB8BjwKNVtTCOoiRJg43jFMpLqurbY/g9kqQN8BSKJDVq1AAv4B+S3J5kx1odkuxIsphkcWlpacTNSZIOGTXAz66qnwVeAbwhyYtWd6iq3VW1UFULc3NzI25OknTISAFeVfd3Pw8CnwLOHEdRkqTBhg7wJE9McvyhZeAXgLvHVZgk6chGuQrl6cCnkhz6PR+uqs+OpSpJ0kBDB3hV3Qu8YIy1SJI2wMsIJalRBrgkNcoAl6RGGeCS1Ci/TlbSj42j7St0PQKXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1Cin0kuaullOaT+aeAQuSY0aKcCTnJvkq0m+nmTnuIqSJA02yk2NtwDvB14BnA5cnOT0cRUmSTqyUY7AzwS+XlX3VtV/Ax8FXjWesiRJg4zyIebJwL+vWL8P+PnVnZLsAHZ0qw8n+eqqLicB3x6hjs3CcWwujmPzOVrGMtQ48p6RtvlTazWOEuBZo60Oa6jaDexe95cki1W1MEIdm4Lj2Fwcx+ZztIxlM41jlFMo9wGnrlg/Bbh/tHIkSX2NEuD/DGxLclqS44CLgGvHU5YkaZChT6FU1aNJ3gj8PbAFuLKq9g3xq9Y9vdIYx7G5OI7N52gZy6YZR6oOO20tSWqAMzElqVEGuCQ1aioBPmjKfZIXJflSkkeTXDiNmobRYxxvSfKVJHcmuTHJmtdubgY9xvL6JHcl2Zvkls06y7bv1zkkuTBJJdkUl3+t1mN/vDbJUrc/9ib5jVnUOUif/ZHkV7u/k31JPjztGvvqsU8uW7E/vpbkoakXWVUTfbD8Aec3gJ8GjgPuAE5f1WceeD7wt8CFk65pguN4CfCT3fJvAR+bdd0jjOXJK5ZfCXx21nUPM46u3/HAF4BbgYVZ1z3k/ngt8NezrnUM49gGfBk4sVt/2qzrHuW9taL/b7N8IcdU65zGEfjAKfdVdaCq7gR+OIV6htVnHDdV1X92q7eyfG38ZtRnLN9fsfpE1piktQn0/TqHPwL+DPivaRa3AUfL11L0GcdvAu+vqu8BVNXBKdfY10b3ycXAR6ZS2QrTCPC1ptyfPIXtjttGx3Ep8JmJVjS8XmNJ8oYk32A5/H5nSrVtxMBxJDkDOLWqrptmYRvU9731y93puauTnLrG87PWZxzPAp6V5B+T3Jrk3KlVtzG9/967U6WnAZ+fQl0/YhoB3mvKfQN6jyPJrwMLwJ9PtKLh9f0ahPdX1TOBtwG/P/GqNu6I40jyOOAy4K1Tq2g4ffbH3wHzVfV84HPAnolXtXF9xnEMy6dRXszyUesHk5ww4bqGsZHcugi4uqoem2A9a5pGgB8tU+57jSPJy4B3Aq+sqkemVNtGbXSffBS4YKIVDWfQOI4HngfcnOQAcBZw7Sb8IHPg/qiq76x4P10O/NyUatuIPu+r+4Brqup/quqbwFdZDvTNZiN/Ixcxg9MnwFQ+xDwGuJfl/8U49GHAz6zT9yo274eYA8cBnMHyBx/bZl3vGMaybcXyLwGLs657lPdW1/9mNueHmH32x9YVy68Gbp113UOO41xgT7d8EsunKZ4669qHfW8BzwYO0E2KnHqdU/qPcR7wtS7c3tm1/SHLR6kAL2T5X7z/AL4D7Jv1DhxyHJ8DHgT2do9rZ13zCGP5K2BfN46bjhSMm3kcq/puygDvuT/+tNsfd3T74zmzrnnIcQR4H/AV4C7golnXPMp7C3gXsGtWNTqVXpIa5UxMSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIa9b8cDqNDsVri1QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(df.tcf_norm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Defining parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_ss = 100\n",
    "dt = .1  # Time step.\n",
    "T = 10  # Total time.\n",
    "tolerance = 1e-8\n",
    "max1 = 1\n",
    "max2 = 1\n",
    "n = int(T / dt)\n",
    "t = np.linspace(0., T, n)  # Vector of times.\n",
    "sim_stages = np.linspace(0,n-1,num_ss, dtype = int)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-23-a240a0655fd5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mx_vals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0my_vals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mx\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[0m\u001b[1;32m      3\u001b[0m \u001b[0mcoefs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpolyfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_vals\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0my_vals\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mxf\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mnew_line\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpolyval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcoefs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mxf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'x' is not defined"
     ]
    }
   ],
   "source": [
    "x_vals = [0, n]\n",
    "y_vals = [x, 0]\n",
    "coefs = np.polyfit(x_vals,y_vals,5)\n",
    "xf=np.linspace(0,n,n)\n",
    "new_line = np.polyval(coefs, xf)\n",
    "plt.scatter(xf,new_line,c='g', marker='^', s=5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def em_cle_tcf_alpha1 (n_traj, x_init, y_init, tolerance, dt , T, alpha2, a, c, b, d, max1, max2, lambda1, lambda2, tcf_init, tcf_group, seed = 20):\n",
    "    np.random.seed(seed)\n",
    "    n = int(T / dt)\n",
    "    t = np.linspace(0., T, n)\n",
    "    sqrtdt =  np.sqrt(dt/omega) \n",
    "    x = np.zeros((n_traj,2,n))\n",
    "    x[:,:,0]= [x_init,y_init]\n",
    "    diverging_time = n/3\n",
    "    \n",
    "    # Euler-Marayama method\n",
    "    for k in range(n_traj):\n",
    "        if tcf_group == 0:\n",
    "            x_vals = [0, n]\n",
    "            y_vals = [tcf_init, 0]\n",
    "            coefs = np.polyfit(x_vals,y_vals,5)\n",
    "            xf=np.linspace(0,n,n)\n",
    "            alpha_vals = np.polyval(coefs, xf)\n",
    "            \n",
    "            for i in range(n - 1):\n",
    "                pp = protein_production([x[k][0][i], x[k][1][i]],alpha_vals[i],alpha2,a,c,b,d, max1, max2)\n",
    "                de = decay([x[k][0][i], x[k][1][i]],lambda1,lambda2)\n",
    "                detT = deterministic_term(pp,de)\n",
    "                stochastic_term = np.array([np.sqrt(pp[0] + de[0]), np.sqrt(pp[1] + de[1])])\n",
    "                x[k][0][i + 1] = x[k][0][i] + dt * detT[0] + stochastic_term[0] * sqrtdt * np.random.randn() #Chemical langevin equation\n",
    "                x[k][1][i + 1] = x[k][1][i] + dt * detT[1] + stochastic_term[1] * sqrtdt * np.random.randn()\n",
    "\n",
    "                # Tolerance - setting small values and negative values to 0 to account for negative values from the standard normal\n",
    "                if(x[k][0][i + 1] < tolerance):\n",
    "                    x[k][0][i + 1] = 0\n",
    "                elif(x[k][1][i + 1] < tolerance):\n",
    "                    x[k][1][i + 1] = 0\n",
    "                    \n",
    "            elif tcf_group == 0:\n",
    "                x_vals = [0, n]\n",
    "                y_vals = [tcf_init, 1]\n",
    "                coefs = np.polyfit(x_vals,y_vals,5)\n",
    "                xf=np.linspace(0,n,n)\n",
    "                alpha_vals = np.polyval(coefs, xf)\n",
    "                \n",
    "                           \n",
    "                for i in range(n - 1):\n",
    "                    pp = protein_production([x[k][0][i], x[k][1][i]],alpha_vals[i],alpha2,a,c,b,d, max1, max2)\n",
    "                    de = decay([x[k][0][i], x[k][1][i]],lambda1,lambda2)\n",
    "                    detT = deterministic_term(pp,de)\n",
    "                    stochastic_term = np.array([np.sqrt(pp[0] + de[0]), np.sqrt(pp[1] + de[1])])\n",
    "                    x[k][0][i + 1] = x[k][0][i] + dt * detT[0] + stochastic_term[0] * sqrtdt * np.random.randn() #Chemical langevin equation\n",
    "                    x[k][1][i + 1] = x[k][1][i] + dt * detT[1] + stochastic_term[1] * sqrtdt * np.random.randn()\n",
    "\n",
    "                    # Tolerance - setting small values and negative values to 0 to account for negative values from the standard normal\n",
    "                    if(x[k][0][i + 1] < tolerance):\n",
    "                        x[k][0][i + 1] = 0\n",
    "                    elif(x[k][1][i + 1] < tolerance):\n",
    "                        x[k][1][i + 1] = 0\n",
    "                        \n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stoch model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`batch`: Total number of rounds to iterate the entire run. Equivalent to an entire HCR time-course experiment from 18ss to 28ss.\n",
    "\n",
    "`num_nmps_time`: Number of NMps / nmp_numbers (total number of cells / embryo) / batch calculated at `num_ss` intervals.\n",
    "\n",
    "`num_nmps_embryos_time`: `num_nmps_time` of all embryos\n",
    "\n",
    "`batch_nmp_num`: `num_nmps_embryos_time` for all batches.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Autonomous "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stoch_model(alpha1, alpha2, a, c, b, d, lambda1, lambda2, omega, exit_nmp_val):\n",
    "    batch_num = 0\n",
    "    batch = 1\n",
    "    batch_nmp_num = [[] for i in range(batch)]\n",
    "\n",
    "    while batch_num < batch:\n",
    "        trajs_list = []\n",
    "        i = 0\n",
    "        while i < nmp_numbers:\n",
    "            xinit = nmp_list[i][0]\n",
    "            yinit = nmp_list[i][1]\n",
    "            single_cell_sim = em_cle(1, xinit, yinit, tolerance, dt , T, alpha1, alpha2, a, c, b, d, max1, max2, lambda1, lambda2, omega, seed = batch_num)\n",
    "            trajs_list.append(single_cell_sim)\n",
    "            i+=1\n",
    "        #print(\"Embryo simulation complete: {} cells\".format(nmp_numbers))\n",
    "\n",
    "        num_nmps_time = []\n",
    "        for i in range(num_ss) :\n",
    "            nmp_counter = 0\n",
    "            traj = np.squeeze(np.array(trajs_list))[:,:,sim_stages[i]]\n",
    "            for cell in traj:\n",
    "                if not ((cell[0] <= exit_nmp_val and cell[1] >= exit_nmp_val) or (cell[0] >= exit_nmp_val and cell[1] <= exit_nmp_val)):\n",
    "                    nmp_counter += 1\n",
    "            num_nmps_time.append(nmp_counter)\n",
    "        #print(\"num_nmps_time\", num_nmps_time)\n",
    "        batch_nmp_num[batch_num].append(num_nmps_time)\n",
    "        batch_num +=1\n",
    "    return batch_nmp_num\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TCF model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stoch_model_tcf(alpha1, alpha2, a, c, b, d, lambda1, lambda2, omega, exit_nmp_val):\n",
    "    batch_num = 0\n",
    "    batch = 1\n",
    "    batch_nmp_num = [[] for i in range(batch)]\n",
    "\n",
    "    while batch_num < batch:\n",
    "        trajs_list = []\n",
    "        i = 0\n",
    "        while i < nmp_numbers:\n",
    "            xinit = nmp_list_tcf[i][0]\n",
    "            yinit = nmp_list_tcf[i][1]\n",
    "            tcf_init = nmp_list_tcf[i][2]\n",
    "            tcf_group = nmp_list_tcf[i][3]\n",
    "            single_cell_sim = em_cle_tcf_alpha1(1, xinit, yinit, tolerance, dt , T, alpha1, alpha2, a, c, b, d, max1, max2, lambda1, lambda2, omega, tcf_init = tcf_init, tcf_group = tcf_group, seed = batch_num)\n",
    "            trajs_list.append(single_cell_sim)\n",
    "            i+=1\n",
    "        #print(\"Embryo simulation complete: {} cells\".format(nmp_numbers))\n",
    "\n",
    "        num_nmps_time = []\n",
    "        for i in range(num_ss) :\n",
    "            nmp_counter = 0\n",
    "            traj = np.squeeze(np.array(trajs_list))[:,:,sim_stages[i]]\n",
    "            for cell in traj:\n",
    "                if not ((cell[0] <= exit_nmp_val and cell[1] >= exit_nmp_val) or (cell[0] >= exit_nmp_val and cell[1] <= exit_nmp_val)):\n",
    "                    nmp_counter += 1\n",
    "            num_nmps_time.append(nmp_counter)\n",
    "        #print(\"num_nmps_time\", num_nmps_time)\n",
    "        batch_nmp_num[batch_num].append(num_nmps_time)\n",
    "        batch_num +=1\n",
    "    return batch_nmp_num\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    " md = np.mean(batch_nmp_num, axis = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(95,)"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(((md - data)/0.25)**2)[0][5:].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3848738.8647658397"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum((((md - data)/0.25)**2)[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3769907.693793953"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum((((md - data)/0.25)**2)[0][5:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MCMC functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import data\n",
    "\n",
    "def logposterior(theta,sigma, x):\n",
    "    lp = logprior(theta)\n",
    "    if not np.isfinite(lp):\n",
    "        return -np.inf\n",
    "    return lp + loglikelihood(theta, sigma,x)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loglikelihood(theta, sigma,t):\n",
    "    alpha1, alpha2, a, c, b, d, lambda1, lambda2, omega, exit_nmp_val = theta\n",
    "    \n",
    "    batch_nmp_num = stoch_model(alpha1, alpha2, a, c, b, d, lambda1, lambda2, omega, exit_nmp_val)\n",
    "    md = np.mean(batch_nmp_num, axis = 0)\n",
    "    return -0.5*np.sum((((md - data)/0.25)**2)[0])\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add evidence of activation or repression\n",
    "def logprior(theta):\n",
    "    alpha1, alpha2, a, c, b, d, lambda1, lambda2, omega, exit_nmp_val = theta\n",
    "    if 0 < alpha1 < 1 and \\\n",
    "       0 < alpha2 < 1  and \\\n",
    "       0  < a < 1  and \\\n",
    "       0 < c <1  and \\\n",
    "       0  < b < 1  and \\\n",
    "       0  < d < 1  and \\\n",
    "       0 < lambda1 < 5 and \\\n",
    "       0  < lambda2 < 5  and \\\n",
    "       10  < omega < 1000  and \\\n",
    "       0 < exit_nmp_val < 0.3:\n",
    "        lp = 0.\n",
    "    else:\n",
    "        lp = -np.inf\n",
    "    return lp\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run MCMC model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nens = 100\n",
    "# Initial (Nens) points\n",
    "alpha1_ini = np.random.uniform(0,1, Nens) \n",
    "alpha2_ini = np.random.uniform(0,1, Nens) \n",
    "a_ini = np.random.uniform(0,1, Nens) \n",
    "c_ini = np.random.uniform(0,1, Nens) \n",
    "b_ini = np.random.uniform(0,1, Nens) \n",
    "d_ini = np.random.uniform(0,1, Nens) \n",
    "lambda1_ini = np.random.uniform(0,5, Nens) \n",
    "lambda2_ini = np.random.uniform(0,5, Nens) \n",
    "omega_ini = np.random.uniform(10,1000, Nens) \n",
    "exit_nmp_val_ini = np.random.uniform(0,0.3, Nens) \n",
    "\n",
    "inisamples = np.array([alpha1_ini, alpha2_ini, a_ini, \\\n",
    "                       c_ini, b_ini, d_ini,\\\n",
    "                       lambda1_ini, lambda2_ini, omega_ini,\\\n",
    "                       exit_nmp_val_ini]).T\n",
    "\n",
    "nwalkers, ndim = inisamples.shape\n",
    "ndims = inisamples.shape[1] \n",
    "Nburnin = 0   \n",
    "Nsamples = 100\n",
    "dsigma = 0.25\n",
    "argslist = (dsigma, sim_stages)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Backend saving"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "filename = \"test_saving.h5\"\n",
    "backend = emcee.backends.HDFBackend(filename)\n",
    "backend.reset(nwalkers, ndim) #In case file name exists already- this clears it.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multiprocessing "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████████████████████████████████████████████████████████████| 3/3 [00:26<00:00,  8.87s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "One run took 36.3 seconds\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "with Pool() as pool:\n",
    "    start = time.time()\n",
    "    sampler = emcee.EnsembleSampler(Nens, ndims, logposterior, args = argslist, pool = pool, moves = emcee.moves.StretchMove(a = 1.3))\n",
    "    sampler.run_mcmc(inisamples, Nsamples+Nburnin, progress=True)\n",
    "    end = time.time()\n",
    "    time_taken = end - start\n",
    "    print(\"One run took {0:.1f} seconds\".format(time_taken))\n",
    "postsamples = sampler.chain[:, Nburnin:, :].reshape((-1, ndims))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MPI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from schwimmbad import MPIPool\n",
    "\n",
    "with MPIPool() as pool:\n",
    "    if not pool.is_master():\n",
    "        pool.wait()\n",
    "        sys.exit(0)\n",
    "    np.random.seed(42)\n",
    "\n",
    "    inisamples = np.array([alpha1_ini, alpha2_ini, a_ini, \\\n",
    "                           c_ini, b_ini, d_ini,\\\n",
    "                           lambda1_ini, lambda2_ini, omega_ini,\\\n",
    "                           exit_nmp_val_ini]).T\n",
    "\n",
    "    nwalkers, ndim = inisamples.shape\n",
    "    ndims = inisamples.shape[1] \n",
    "    Nburnin = 0   \n",
    "    Nsamples = 3\n",
    "    dsigma = 0.25\n",
    "    argslist = (dsigma, sim_stages)\n",
    "    np.random.seed(42)\n",
    "    initial = np.random.randn(32, 5)\n",
    "    nwalkers, ndim = initial.shape\n",
    "    nsteps = 100\n",
    "\n",
    "    start = time.time()\n",
    "    sampler = emcee.EnsembleSampler(Nens, ndims, logposterior, args = argslist, pool = pooL)\n",
    "    sampler.run_mcmc(inisamples, Nsamples+Nburnin, progress=True)\n",
    "    end = time.time()\n",
    "    time_taken = end - start\n",
    "    print(\"One run took {0:.1f} seconds\".format(time_taken))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ndims"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Results "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MCMC properties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Viewing parameter distributions by accessing chain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'sampler' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-2de1bb96f0a7>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0msamples\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msampler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_chain\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mflat\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mhist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msamples\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"k\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhisttype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"step\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mxlabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mr\"$\\theta_1$\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mylabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mr\"$p(\\theta_1)$\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgca\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_yticks\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\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;31mNameError\u001b[0m: name 'sampler' is not defined"
     ]
    }
   ],
   "source": [
    "samples = sampler.get_chain(flat=True)\n",
    "plt.hist(samples[:, 3], 100, color=\"k\", histtype=\"step\")\n",
    "plt.xlabel(r\"$\\theta_1$\")\n",
    "plt.ylabel(r\"$p(\\theta_1)$\")\n",
    "plt.gca().set_yticks([]);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(72, 10)"
      ]
     },
     "execution_count": 149,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sampler.get_chain(flat=True).shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Viewing Mean acceptance fraction of chain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean acceptance fraction: 0.222\n"
     ]
    }
   ],
   "source": [
    "print(\"Mean acceptance fraction: {0:.3f}\".format(np.mean(sampler.acceptance_fraction)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### View integrated autocorrelated time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean autocorrelation time: nan steps\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/porcupine/miniconda3/lib/python3.7/site-packages/emcee/autocorr.py:36: RuntimeWarning: invalid value encountered in true_divide\n",
      "  acf /= acf[0]\n",
      "/home/porcupine/miniconda3/lib/python3.7/site-packages/emcee/autocorr.py:41: RuntimeWarning: invalid value encountered in less\n",
      "  m = np.arange(len(taus)) < c * taus\n",
      "/home/porcupine/miniconda3/lib/python3.7/site-packages/emcee/autocorr.py:99: RuntimeWarning: invalid value encountered in greater\n",
      "  flag = tol * tau_est > n_t\n"
     ]
    }
   ],
   "source": [
    "print(\n",
    "    \"Mean autocorrelation time: {0:.3f} steps\".format(\n",
    "        np.mean(sampler.get_autocorr_time())\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Corner plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n",
      "WARNING:root:Too few points to create valid contours\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Embryo simulation complete: 90 cells\n",
      "num_nmps_time [35, 44, 65, 74, 76, 78, 79, 79, 78, 78, 78, 78, 77, 77, 75, 75, 74, 74, 74, 74, 74, 73, 73, 73, 73, 73, 73, 73, 70, 70, 69, 69, 68, 68, 67, 67, 67, 67, 66, 66, 65, 64, 63, 63, 61, 60, 59, 60, 58, 57, 57, 57, 57, 56, 57, 54, 51, 50, 49, 46, 45, 44, 43, 43, 43, 43, 41, 36, 36, 35, 35, 35, 35, 35, 34, 33, 33, 33, 30, 29, 28, 28, 28, 27, 25, 25, 24, 22, 22, 21, 20, 20, 18, 16, 16, 15, 15, 14, 13, 12]\n",
      "Embryo simulation complete: 90 cells\n",
      "num_nmps_time [35, 44, 71, 80, 81, 82, 81, 82, 82, 82, 82, 81, 81, 81, 81, 81, 80, 79, 78, 78, 78, 78, 78, 78, 77, 76, 76, 76, 75, 74, 74, 73, 73, 73, 72, 72, 71, 71, 71, 70, 69, 67, 67, 67, 66, 65, 65, 65, 65, 65, 65, 65, 64, 63, 63, 63, 62, 62, 62, 62, 61, 59, 59, 59, 56, 55, 53, 51, 50, 50, 48, 46, 45, 44, 45, 44, 44, 44, 42, 39, 38, 36, 34, 34, 34, 34, 34, 34, 31, 31, 26, 25, 25, 25, 24, 24, 23, 23, 21, 20]\n",
      "Something else went wrong. Cannot plot\n"
     ]
    }
   ],
   "source": [
    "# Plot corner plot\n",
    "try:\n",
    "    import matplotlib as mpl\n",
    "    mpl.use(\"Agg\") # force Matplotlib backend to Agg\n",
    "    import corner \n",
    "except ImportError:\n",
    "    sys.exit(1)\n",
    "\n",
    "fig = corner.corner(postsamples,\\\n",
    "                    labels=[r\"$alpha1$\",r\"$alpha2$\",r\"$a$\",\\\n",
    "                            r\"$c$\",r\"$b$\",r\"$d$\",\\\n",
    "                            r\"$lambda1$\",r\"$lambda2$\",r\"$omega$\",\\\n",
    "                            r\"$exit_nmp_val$\"],\\\n",
    "                    truths=None)\n",
    "fig.savefig('emcee.png' )\n",
    "\n",
    "# Export the best set of parameters\n",
    "samples = sampler.flatchain\n",
    "params = samples[np.argmax(sampler.flatlnprobability)]\n",
    "np.save('140121_fitted_params.npy', params)    \n",
    "###    \n",
    "# If this emcee is not fed into the loop,\n",
    "# Need to output the results by uncommenting the following \n",
    "\n",
    "try:\n",
    "    plt.close(fig)\n",
    "    y_values = stoch_model(alpha1, alpha2, a, c, b, d, lambda1, lambda2, omega, exit_nmp_val)\n",
    "    f2 = plt.figure(1)\n",
    "    plt.plot(sim_stages, data, label = 'data')\n",
    "    plt.plot(sim_stages, y_values, 'b-', label = 'model')\n",
    "    plt.ylabel('NMp numbers')\n",
    "    plt.xlabel('Simulation_stages')\n",
    "    plt.legend(loc=\"best\")\n",
    "    plt.close(f2)\n",
    "    f2.savefig('plot.png')\n",
    "except:\n",
    "    print(\"Something else went wrong. Cannot plot\")\n",
    "    \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load fitted params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### From jupyter notebook example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = sampler.flatchain\n",
    "params = samples[np.argmax(sampler.flatlnprobability)]\n",
    "np.save('files/npy/060421_fitted_params.npy', params)  \n",
    "fitted_params = np.load(\"files/npy/060421_fitted_params.npy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Saved from nohup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10,)"
      ]
     },
     "execution_count": 134,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fitted_params = np.load(\"files/results/map_params_emcee.txt\" ,  allow_pickle=True)\n",
    "alpha1, alpha2, a, c, b, d, lambda1, lambda2, omega, exit_nmp_val = fitted_params\n",
    "fitted_params.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAD4CAYAAADSIzzWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAKLElEQVR4nO3dX4ite13H8c/XNhaB/aGtEGnzJCh08CY5hN5YYYScC8+NhIFUcAg06KarwJuou6CCQFAvpD9QWV3UJuqmUgzpWPugqQnGyWbZIUmlOjdRKf66WGuPy2lm72dmr/V8Z816vWBgzczas36/tZ55z+95nrXWrjFGAFjeS7oHAHCsBBigiQADNBFggCYCDNDkzlWufPfu3TFN056GAnA7Pffcc18eY7z8/NevFOBpmnL//v3djQrgCFTV6qKvOwQB0ESAAZoIMEATAQZoIsAATQQYoIkAAzQRYIAmAgzQRIB5bNM0papSVfFSdZhPgHlsq9UqY4yMMbJaXfiKS5GGC1zpvSDguh5EOkmqqnk0cDNYAQM0EWCAJgIM0ESAAZoIMEATAQZoIsAATQQYoIkAAzQRYIAmAgzQRIABmggwQBMBBmgiwABNBBigiQADNBFggCYCDNBEgAGaCDBAEwEGaCLAAE0EGKCJAAM0EWCAJgIM0ESAAZoIMEATAQZoIsAATQQYoIkAAzQRYIAmAgzQRIABmggwQBMBBmgiwABNBBigiQADNBFggCYCDNBEgAGaCDBAEwEGaCLAAE0EGKCJAAM0EWCAJgIM0ESAAZoIMEATAQZoIsAATQQYoIkAAzQRYIAmAgzQRIABmggwQBMBhkbTNKWqUlWZpql7OCzsTvcA4JitVquMMZIkVdU8GpZmBQzQRIABmggwQBMBBmgiwABNBBigiQADNBFggCYCDNBEgIEzXhq9LC9FBs54afSyrIABmggwQBMBZqdOTk4cQ7wm993xcQyYnTo9PT277Bji1bjvjo8VMEATAQZoIsAATQQYoIkAAzQRYIAmAgzQRIABmggwQBMBBmgiwABNBBigiQADNBFggCYCDNBEgAGaCDBAEwEGaCLAAE0EGKCJAAM0EWCAJgIM0ESAAZoIMEATAQZoIsAATQQY9mSaplRVqirTNHUPhxvoTvcA4LZarVYZYyRJqqp5NNxEVsAATQQYoIkAAzRxDJi9OTk5OTv2eXJy0jwauHkEmL05PT3tHgLcaA5BADQRYIAmAgzQRIABmggwQBMBBmgiwABNBBigiQADNBFggCYCDNBEgAGaCDBAEwEGaCLAAE0EGKCJAAM0EWCAJgIM0ESAAZoI8A01TVOqKlWVaZq6hwPsgf8V+YZarVYZYyTJ2X/tDtwuVsDcSPYAOAZWwNxI9gA4BlbA+cbVlhUXsBQr4Hzjaiux4gKWYQUM0ESAj5STXNDPIYgj5SQX9LMCBmgiwABNBBigiQADNBFggCYCDNBEgAGaCDBAEwEGaCLAAE0EGKCJAANtjv1NobwZD9Dm2N8UygoYOEi3YfVsBQwcpNuwerYCBmgiwABNBBigiQBf0204AQD0chLumm7DCQCglxUwQBMBBmgiwABNBBigiQADNBFggCYCDNBEgAGaCDBAEwEGaCLAAE0EGKCJAAM0EWCAJgJ8BdvvAXxyctI9HODAeT/gK9h+D2CAx2UFDNBEgAGaCDBAEwEGaCLAAE0EGKCJAAM0EWCAJgIM0ESAAZoIMLAo76nydd4LAliU91T5OitggCYCDDtk93q/btv96xAE7JDd6/26bfevFTBHaXslNU1T93A4UlbAHKXtlVRVNY+GY2UFzGxWjbBbVsBHZJqmrFarJLnWCQyrRtgtK+Aj8iCgY4ycnp52D4c9sadyOAQYboHt6CY5+0P7YI+ne0z+EFxMgGFhV30u65yQ3cS9m+0xdf4huMkcA+ZaTk5OzlZbt+EJ8Uu66nNZHXu/vQSYa1lylbUd+/NfvymrvUc59j9Yxz7/yzgEcYEHG8uxH7va3vXtfOnn6enp2a7s9sch7dZuz2HOH43tbXAf9/vSx2evOv9jsViAD+mA/PbGcki/5Ltw2cmc2/CLM+fY603ZTi8L1lUXB5ddf87xWQuR/VsswId6QH7OSmTfG+plP39OLOZEZ/vnJ9l7cLdv77KPOau+OT9n++OyuV02/8u20/N7BkvuJWyHOckjb3fOYuKybfyyf3vZdve4e0xztvPbdviirnIyoKq+lOSq9byb5MtX/De3wTHO25yPxzHO+3HmfDLGePn5L14pwNdRVffHGE/u9UZuoGOctzkfj2Oc9z7m7CQcQBMBBmiyRIDfv8Bt3ETHOG9zPh7HOO+dz3nvx4ABuJhDEABNBBigyc4CXFVvqarPVtXzVfULF3z/m6vqg5vvf6yqpl3ddpcZc/75qvpMVX2yqv6yqm7Fs8gfNe+t672tqkZVHfzTlebMuap+fPN4/0NV/e7SY9y1Gdv391bVh6rq45tt/KmOce5SVX2gqr5YVZ++5PtVVb+xuU8+WVWvf6wbvOg19lf9SPJNSf4pyauTvDTJ3yd54tx1fjbJezeX357kg7u47a6PmXP+kSTfurn8rkOf89x5b673siQfSfJskie7x73AY/2aJB9P8p2bz1/RPe4F5vz+JO/aXH4iyWn3uHcw7zcleX2ST1/y/aeS/HmSSvKGJB97nNvb1Qr4B5M8P8b43Bjjf5P8fpKnz13n6SS/tbn8R0neXA9e+3mYHjnnMcaHxhj/tfn02SSvXHiM+zDnsU6SX07yK0n+e8nB7cmcOf9MkveMMf4jScYYX1x4jLs2Z84jybdtLn97kn9dcHx7Mcb4SJJ/f8hVnk7y22Pt2STfUVXffd3b21WAvyfJv2x9/sLmaxdeZ4zx1SQvJvmuHd1+hzlz3vZM1n85D90j511VP5DkVWOMP11yYHs057F+bZLXVtVHq+rZqnrLYqPbjzlz/sUk76iqF5L8WZKfW2Zora76e/9Qu3o/4ItWsuef3zbnOodk9nyq6h1JnkzyQ3sd0TIeOu+qekmSX0/y00sNaAFzHus7WR+G+OGs93T+uqpeN8b4zz2PbV/mzPknkvzmGONXq+qNSX5nM+ev7X94bXbasV2tgF9I8qqtz1+Z/787cnadqrqT9S7Lw5b6N92cOaeqfjTJu5O8dYzxPwuNbZ8eNe+XJXldkg9X1WnWx8nuHfiJuLnb95+MMb4yxvjnJJ/NOsiHas6cn0nyB0kyxvibJN+S9RvW3Gazfu/n2lWA/y7Ja6rq+6rqpVmfZLt37jr3kvzU5vLbkvzV2BzVPlCPnPNmV/x9Wcf30I8JPvDQeY8xXhxj3B1jTGOMKetj328dY9zvGe5OzNm+/zjrk66pqrtZH5L43KKj3K05c/58kjcnSVV9f9YB/tKio1zevSQ/uXk2xBuSvDjG+MK1f9oOzx4+leQfsz5z+u7N134p61++ZP3g/GGS55P8bZJXd5/xXGDOf5Hk35J8YvNxr3vMS8z73HU/nAN/FsTMx7qS/FqSzyT5VJK3d495gTk/keSjWT9D4hNJfqx7zDuY8+8l+UKSr2S92n0myTuTvHPrcX7P5j751ONu216KDNDEK+EAmggwQBMBBmgiwABNBBigiQADNBFggCb/B7Qdv67bZcVwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "samples = np.load(\"files/results/samples_emcee.txt\" ,  allow_pickle=True)\n",
    "plt.hist(samples[:, 3], 100, color=\"k\", histtype=\"step\")\n",
    "#plt.xlabel(r\"$\\theta_1$\")\n",
    "#plt.ylabel(r\"$p(\\theta_1)$\")\n",
    "plt.gca().set_yticks([]);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(240000, 10)"
      ]
     },
     "execution_count": 168,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "samples.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([423.17156578, 583.93395984, 600.45008394, 446.47465559,\n",
       "       600.43979306, 449.43914778, 605.65964543, 400.22024519,\n",
       "       529.04925286, 411.48472675])"
      ]
     },
     "execution_count": 136,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "autocorr_time = np.load(\"files/results/autocorr_time_emcee.txt\" ,  allow_pickle=True)\n",
    "autocorr_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acceptance_fraction = np.load(\"files/results/acceptance_fraction_emcee.txt\" ,  allow_pickle=True)\n",
    "acceptance_fraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-5597042.8429107 , -6058490.15452306, -8619627.28374162, ...,\n",
       "        -109148.62934078,   -74867.23749369,  -301388.08785582])"
      ]
     },
     "execution_count": 144,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_probs = np.load(\"files/results/log_probs.txt\" ,  allow_pickle=True)\n",
    "log_probs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha1, alpha2, a, c, b, d, lambda1, lambda2, omega, exit_nmp_val = fitted_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f1e674324a8>]"
      ]
     },
     "execution_count": 171,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD4CAYAAADsKpHdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAARqUlEQVR4nO3dfbBcdX3H8feXxMQHkAdzoZk8mKDRGh1b4jXC0PowoiY4JWNL26TTgi1jahXbjlomFJs6dBgR2zq1jUpAxupUMD7VtF6MrdCxrYK5DBAJaco1RnMNY26AUpUKBL79Y38Jy2b33k3uPtxz7/s1c2fP+Z3fnv3+cpb9cB52T2QmkiSd0O8CJElTg4EgSQIMBElSYSBIkgADQZJUzO7XC8+bNy+XLFnSr5eXpEq64447DmbmQDfW3bdAWLJkCcPDw/16eUmqpIj4frfW7SEjSRJgIEiSCgNBkgQYCJKkYsJAiIgbIuJARNzTYnlExEciYiQidkTEis6XKUnqtnb2ED4JrBpn+WpgWflbD3xs8mVJknptwkDIzG8AD47TZQ3wqay5DTglIuZ3qkBJUm904nsIC4B9dfOjpe3+xo4RsZ7aXgSLFy8+7hd85LFDLN+47Wlt9121mmfM8pSIJB2vTnyCRpO2pjdZyMzNmTmYmYMDA8f/RbvGMAD4lb/9j+NenySpM4EwCiyqm18I7O/Aeo/JyIGftFz26KEnuOP7D/WwGkmqnk4EwlbgonK10dnAw5l51OGiXtr34CMs2fAVfvb4EwC8+H1f5dc+9k2WbPhKP8uSpCmtnctObwS+Bbw4IkYj4pKIeHtEvL10GQL2ACPAdcA7ulZtm375mlsB+Pk/+2qfK5Gk6pjwpHJmrptgeQLv7FhFkqS+qNxlOT999FDT9qZnsYHrvrGne8VI0jRSuUB47+fubtre7FIngKuGdnWvGEmaRioXCPeNczWRJOn4VS4QJEndYSBIkoBpFAjR6iSCJKkt0yYQstVlRpKktkybQJAkTc60CQQPGUnS5EybQJAkTU7lAiE9WSBJXVG5QGjFnJCkyenEHdN66slJfvDX/wS2d1mTpKdU7tPwewd/2rF1Lbvi5o6tS5KqrnKB0IpXGUnS5EyfQGj5e6eSpHZMm0CQJE3OtAmEbHmLHElSO6ZNIDz+RPodBUmahMpddjqepZcPsffqNx/Tc+ovQ+2Exac9m29c9rqOrlOSemHa7CEc9q4b7+zr6//gwUfYsn1fX2uQpOMxrfYQAP7p7v39LoHLvrCDy76wo+uvc6x7Q5I0nmm3hyBJOj4GgiQJmIaHjATfHDnIR265jxPK17fPeO4zueotL+PZc9zcklrzE2KaeeLJ5Leuv/2o9i/d+UPPOUgal4eMppknJvtzsJJmLANBkgR4yKjSjvVLdeP193CSJPcQBMAjjx3qdwmS+sxAEDD5O9FJqr62AiEiVkXE7ogYiYgNTZYvjohbI+LOiNgREed3vlRJUjdNGAgRMQvYBKwGlgPrImJ5Q7f3AVsy8yxgLfDRTheqmefLd/2QJRu+wp0/eKjfpUgzQjsnlVcCI5m5ByAibgLWAPfW9UnguWX6ZKD/PyikY/KyP9/W7xJaestHv9nvEtRDEZD51GN9W2OfVuadOJcIjtxHMalNH35cdsaJXH/RK3nWnFldGUNVtRMIC4D6n+8cBV7V0Of9wNci4l3Ac4Dzmq0oItYD6wEWL158rLVKmgEOf9DXf+A3fvhPdOuTgz95dNzlB378KC/Z+NXjqK4zzjnzedy4/uy+vX4r7ZxDaHaz4sbNsQ74ZGYuBM4HPh0RR607Mzdn5mBmDg4MDBx7tZI0DXxrzwP9LqGpdgJhFFhUN7+Qow8JXQJsAcjMbwHPBOZ1okBJUm+0EwjbgWURsTQi5lA7aby1oc8PgNcDRMRLqAXCWCcLlSR114SBkJmHgEuBbcAualcT7YyIKyPigtLtPcDbIuJu4EbgrekNjiWpUtr66YrMHAKGGto21k3fC5zb2dIkSb3kN5UlSYCBIEkqDARJEmAgSJIKA0GSBBgIkqTCQJAkAQaCJKkwECRJgIEgSSoMBEkSYCBIkgoDQZIEGAiSpMJAkCQBBoIkqTAQJEmAgSBJKgwESRJgIEiSCgNBkgQYCJKkwkCQJAEGgiSpMBAkSYCBIEkqDARJEmAgSJIKA0GSBBgIkqTCQJAkAW0GQkSsiojdETESERta9PmNiLg3InZGxGc6W6YkqdtmT9QhImYBm4A3AKPA9ojYmpn31vVZBlwOnJuZD0XE6d0qWJLUHe3sIawERjJzT2Y+BtwErGno8zZgU2Y+BJCZBzpbpiSp29oJhAXAvrr50dJW70XAiyLiPyPitohY1WxFEbE+IoYjYnhsbOz4KpYkdUU7gRBN2rJhfjawDHgtsA64PiJOOepJmZszczAzBwcGBo61VklSF7UTCKPAorr5hcD+Jn2+nJmPZ+b3gN3UAkKSVBHtBMJ2YFlELI2IOcBaYGtDn38EXgcQEfOoHULa08lCJUndNWEgZOYh4FJgG7AL2JKZOyPiyoi4oHTbBjwQEfcCtwJ/kpkPdKtoSVLnTXjZKUBmDgFDDW0b66YTeHf5kyRVkN9UliQBBoIkqTAQJEmAgSBJKgwESRJgIEiSCgNBkgQYCJKkwkCQJAEGgiSpMBAkSYCBIEkqDARJEmAgSJIKA0GSBBgIkqTCQJAkAQaCJKkwECRJgIEgSSoMBEkSYCBIkgoDQZIEGAiSpMJAkCQBBoIkqTAQJEmAgSBJKgwESRJgIEiSirYCISJWRcTuiBiJiA3j9LswIjIiBjtXoiSpFyYMhIiYBWwCVgPLgXURsbxJv5OAPwRu73SRkqTua2cPYSUwkpl7MvMx4CZgTZN+fwFcA/ysg/VJknqknUBYAOyrmx8tbUdExFnAosz85w7WJknqoXYCIZq05ZGFEScAHwbeM+GKItZHxHBEDI+NjbVfpSSp69oJhFFgUd38QmB/3fxJwMuAf4uIvcDZwNZmJ5Yzc3NmDmbm4MDAwPFXLUnquHYCYTuwLCKWRsQcYC2w9fDCzHw4M+dl5pLMXALcBlyQmcNdqViS1BUTBkJmHgIuBbYBu4AtmbkzIq6MiAu6XaAkqTdmt9MpM4eAoYa2jS36vnbyZUmSes1vKkuSAANBklQYCJIkwECQJBUGgiQJMBAkSYWBIEkCDARJUmEgSJIAA0GSVBgIkiTAQJAkFQaCJAkwECRJhYEgSQIMBElSYSBIkgADQZJUGAiSJMBAkCQVBoIkCTAQJEmFgSBJAgwESVJhIEiSAANBklQYCJIkwECQJBUGgiQJMBAkSYWBIEkC2gyEiFgVEbsjYiQiNjRZ/u6IuDcidkTE1yPi+Z0vVZLUTRMGQkTMAjYBq4HlwLqIWN7Q7U5gMDNfDnweuKbThUqSuqudPYSVwEhm7snMx4CbgDX1HTLz1sx8pMzeBizsbJmSpG5rJxAWAPvq5kdLWyuXADc3WxAR6yNiOCKGx8bG2q9SktR17QRCNGnLph0jfhsYBD7UbHlmbs7MwcwcHBgYaL9KSVLXzW6jzyiwqG5+IbC/sVNEnAdcAbwmMx/tTHmSpF5pZw9hO7AsIpZGxBxgLbC1vkNEnAVcC1yQmQc6X6YkqdsmDITMPARcCmwDdgFbMnNnRFwZEReUbh8CTgQ+FxF3RcTWFquTJE1R7RwyIjOHgKGGto110+d1uC5JUo/5TWVJEmAgSJIKA0GSBBgIkqTCQJAkAQaCJKkwECRJgIEgSSoMBEkSYCBIkorKBcLpJ83tdwmSNC1VLhBmndDs9gySpMmqXCBIkrrDQJAkARUMhBPCQ0aS1A2VC4RnzDIQJKkbKhcIkqTuqFwgZL8LkKRpqnKBIEnqDgNBkgQYCJKkonKBkJ5EkKSuqFwgSJK6w0CQJAEVDIT0wlNJ6orKBYIkqTsMBEkSUMFA+MsLf6HfJUjStDS73wUcq1ed+bym7W9YfgbXXTQIwP0P/x/nfOCWXpYlSZVXuUBotGLxKbzppT/H+lefeaRt/snPYu/VbwbghX86xKEnPREtSROp3CEjgM+87VVHpr/4jnP5/de8gGhxn4Qtbz+nV2VJUqW1tYcQEauAvwFmAddn5tUNy+cCnwJeATwA/GZm7u1sqU+Zd+LctvuuWHzqkb2Fw8Z+/CivvOpfO12WJFXahHsIETEL2ASsBpYD6yJieUO3S4CHMvOFwIeBD3a60HrLTj9xUs8fOKn9QJGkmaKdQ0YrgZHM3JOZjwE3AWsa+qwB/r5Mfx54fbQ6htMBXVy1JM1Y7QTCAmBf3fxoaWvaJzMPAQ8DR10OFBHrI2I4IobHxsaOr+Ji79VvPupQ0LH498tex5teegbrVi460nbeS86YVE2S1I4b3jrY7xKaauccQrP/HW+8bKedPmTmZmAzwODgYF8v/Vl02rO59ndqG+UDv/ryfpYiSVNCO3sIo8CiuvmFwP5WfSJiNnAy8GAnCpQk9UY7gbAdWBYRSyNiDrAW2NrQZytwcZm+ELgl0zsXSFKVTHjIKDMPRcSlwDZql53ekJk7I+JKYDgztwKfAD4dESPU9gzWdrNoSVLntfU9hMwcAoYa2jbWTf8M+PXOliZJ6qVKflNZktR5BoIkCTAQJEmFgSBJAiD6dXVoRIwB3z/Op88DDnawnKqZyeOfyWOHmT1+x17z/Mwc6MaL9C0QJiMihjNzan73uwdm8vhn8thhZo/fsXd/7B4ykiQBBoIkqahqIGzudwF9NpPHP5PHDjN7/I69yyp5DkGS1HlV3UOQJHWYgSBJAioYCBGxKiJ2R8RIRGzodz2TERF7I+I7EXFXRAyXttMi4l8i4r7yeGppj4j4SBn3johYUbeei0v/+yLi4rr2V5T1j5Tn9u3eoxFxQ0QciIh76tq6PtZWr9FrLcb//oj4Ydn+d0XE+XXLLi9j2R0Rb6prb/r+Lz9Pf3sZ52fLT9UTEXPL/EhZvqQ3I35KRCyKiFsjYldE7IyIPyrt0377jzP2qbntM7Myf9R+fvu7wJnAHOBuYHm/65rEePYC8xrargE2lOkNwAfL9PnAzdTuTnc2cHtpPw3YUx5PLdOnlmXfBs4pz7kZWN3Hsb4aWAHc08uxtnqNKTL+9wPvbdJ3eXlvzwWWlvf8rPHe/8AWYG2Z/jjwB2X6HcDHy/Ra4LN9GPt8YEWZPgn47zLGab/9xxn7lNz2fflwmMQ/7jnAtrr5y4HL+13XJMazl6MDYTcwv+7NtLtMXwusa+wHrAOurWu/trTNB/6rrv1p/fo03iU8/QOx62Nt9RpTZPytPhSe9r6mdi+Sc1q9/8uH4EFgdmk/0u/wc8v07NIv+vw++DLwhpm2/RvGPiW3fdUOGS0A9tXNj5a2qkrgaxFxR0SsL21nZOb9AOXx9NLeauzjtY82aZ9KejHWVq8xVVxaDovcUHc441jH/zzgfzLzUEP709ZVlj9c+vdFOWxxFnA7M2z7N4wdpuC2r1ogNDsGXuXrZs/NzBXAauCdEfHqcfq2GvuxtlfBTBnrx4AXAL8I3A/8VWnv5PinzL9NRJwIfAH448z83/G6Nmmr9PZvMvYpue2rFgijwKK6+YXA/j7VMmmZub88HgC+BKwEfhQR8wHK44HSvdXYx2tf2KR9KunFWFu9Rt9l5o8y84nMfBK4jtr2h2Mf/0HglIiY3dD+tHWV5SdTu81tT0XEM6h9IP5DZn6xNM+I7d9s7FN121ctELYDy8pZ9TnUTpRs7XNNxyUinhMRJx2eBt4I3ENtPIevnriY2jFHSvtF5QqMs4GHyy7wNuCNEXFq2e18I7VjiPcDP46Is8sVFxfVrWuq6MVYW71G3x3+oCreQm37Q63mteUqkaXAMmonTZu+/7N2kPhW4MLy/MZ/y8PjvxC4pfTvmbJNPgHsysy/rls07bd/q7FP2W3fzxMsx3lS5nxqZ+q/C1zR73omMY4zqV0pcDew8/BYqB3j+zpwX3k8rbQHsKmM+zvAYN26fg8YKX+/W9c+WN5o3wX+jj6eTARupLZr/Di1/3O5pBdjbfUaU2T8ny7j21H+451f1/+KMpbd1F0d1ur9X95P3y7/Lp8D5pb2Z5b5kbL8zD6M/ZeoHarYAdxV/s6fCdt/nLFPyW3vT1dIkoDqHTKSJHWJgSBJAgwESVJhIEiSAANBklQYCJIkwECQJBX/D0rscjAjV5sWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, samples.shape[0], num = samples.shape[0])\n",
    "plt.plot(x, samples[:,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing output (running stoch model with best parameters) against data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Embryo simulation complete: 182 cells\n",
      "num_nmps_time [91, 96, 99, 109, 121, 125, 129, 133, 137, 137, 134, 136, 136, 132, 130, 141, 136, 120, 118, 126, 123, 116, 119, 114, 120, 119, 121, 114, 115, 113, 111, 111, 107, 104, 105, 104, 105, 104, 102, 99, 99, 100, 97, 98, 98, 98, 98, 96, 95, 95, 95, 95, 92, 90, 90, 90, 90, 90, 90, 90, 89, 89, 88, 87, 85, 84, 82, 79, 71, 70, 70, 65, 63, 62, 60, 57, 54, 51, 48, 48, 46, 45, 42, 33, 27, 26, 25, 21, 14, 11, 10, 6, 3, 2, 2, 1, 1, 1, 1, 1]\n"
     ]
    }
   ],
   "source": [
    "num_ss = 100\n",
    "sim_stages = np.linspace(0,n-1,num_ss, dtype = int)\n",
    "dt = .1  # Time step.\n",
    "T = 10  # Total time.\n",
    "tolerance = 1e-8\n",
    "\n",
    "batch_num = 0\n",
    "batch = 1\n",
    "batch_nmp_num = [[] for i in range(batch)]\n",
    "\n",
    "while batch_num < batch:\n",
    "    trajs_list = []\n",
    "    i = 0\n",
    "    while i < nmp_numbers:\n",
    "        xinit = nmp_list[i][0]\n",
    "        yinit = nmp_list[i][1]\n",
    "        single_cell_sim = em_cle(1, xinit, yinit, tolerance, dt , T, alpha1, alpha2, a, c, b, d, max1, max2, lambda1, lambda2, omega, seed = batch_num)\n",
    "        trajs_list.append(single_cell_sim)\n",
    "        i+=1\n",
    "    print(\"Embryo simulation complete: {} cells\".format(nmp_numbers))\n",
    "\n",
    "    num_nmps_time = []\n",
    "    for i in range(num_ss) :\n",
    "        nmp_counter = 0\n",
    "        traj = np.squeeze(np.array(trajs_list))[:,:,sim_stages[i]]\n",
    "        for cell in traj:\n",
    "            if not ((cell[0] <= exit_nmp_val and cell[1] >= exit_nmp_val) or (cell[0] >= exit_nmp_val and cell[1] <= exit_nmp_val)):\n",
    "                nmp_counter += 1\n",
    "        num_nmps_time.append(nmp_counter)\n",
    "    print(\"num_nmps_time\", num_nmps_time)\n",
    "    batch_nmp_num[batch_num].append(num_nmps_time)\n",
    "    batch_num +=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[91,\n",
       " 96,\n",
       " 99,\n",
       " 109,\n",
       " 121,\n",
       " 125,\n",
       " 129,\n",
       " 133,\n",
       " 137,\n",
       " 137,\n",
       " 134,\n",
       " 136,\n",
       " 136,\n",
       " 132,\n",
       " 130,\n",
       " 141,\n",
       " 136,\n",
       " 120,\n",
       " 118,\n",
       " 126,\n",
       " 123,\n",
       " 116,\n",
       " 119,\n",
       " 114,\n",
       " 120,\n",
       " 119,\n",
       " 121,\n",
       " 114,\n",
       " 115,\n",
       " 113,\n",
       " 111,\n",
       " 111,\n",
       " 107,\n",
       " 104,\n",
       " 105,\n",
       " 104,\n",
       " 105,\n",
       " 104,\n",
       " 102,\n",
       " 99,\n",
       " 99,\n",
       " 100,\n",
       " 97,\n",
       " 98,\n",
       " 98,\n",
       " 98,\n",
       " 98,\n",
       " 96,\n",
       " 95,\n",
       " 95,\n",
       " 95,\n",
       " 95,\n",
       " 92,\n",
       " 90,\n",
       " 90,\n",
       " 90,\n",
       " 90,\n",
       " 90,\n",
       " 90,\n",
       " 90,\n",
       " 89,\n",
       " 89,\n",
       " 88,\n",
       " 87,\n",
       " 85,\n",
       " 84,\n",
       " 82,\n",
       " 79,\n",
       " 71,\n",
       " 70,\n",
       " 70,\n",
       " 65,\n",
       " 63,\n",
       " 62,\n",
       " 60,\n",
       " 57,\n",
       " 54,\n",
       " 51,\n",
       " 48,\n",
       " 48,\n",
       " 46,\n",
       " 45,\n",
       " 42,\n",
       " 33,\n",
       " 27,\n",
       " 26,\n",
       " 25,\n",
       " 21,\n",
       " 14,\n",
       " 11,\n",
       " 10,\n",
       " 6,\n",
       " 3,\n",
       " 2,\n",
       " 2,\n",
       " 1,\n",
       " 1,\n",
       " 1,\n",
       " 1,\n",
       " 1]"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(num_nmps_time)[0:100]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'x' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-33-5a9b1d933e03>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnum_nmps_time\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlabel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"Simulated\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'r'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlabel\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"Data\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'x' is not defined"
     ]
    }
   ],
   "source": [
    "plt.plot(x[5:100],(num_nmps_time)[5:100], label = \"Simulated\")\n",
    "plt.plot(x[5:100],data[5:100], color = 'r', label = \"Data\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f1e675ed390>"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3hUxdfA8e9JTyAkQOgt9CYQIHRUELG9ICAgUpQqKCoq1p+9ixUFBEVEUBQpCiJiQQRFmvQuRYqEXkNJQtq8f8wmJJCEtM1ukvN5nvske/fuvWcj7tm5M3NGjDEopZRSAB6uDkAppZT70KSglFIqmSYFpZRSyTQpKKWUSqZJQSmlVDIvVweQEyEhISY0NNTVYSilVL6ydu3aE8aYUmk9l6+TQmhoKGvWrHF1GEopla+IyP70ntPbR0oppZJpUlBKKZVMk4JSSqlk+bpPQSlVcMTFxREREUFMTIyrQykw/Pz8qFixIt7e3pl+jSYFpZRbiIiIIDAwkNDQUETE1eHke8YYTp48SUREBFWrVs306/T2kVLKLcTExFCyZElNCLlERChZsmSWW15OSwoiMllEjonIljSee1xEjIiEOB6LiIwRkd0isklEmjgrLqWU+9KEkLuy8/d0ZkthCnDL5TtFpBLQEfgvxe5bgZqObSgwwYlxuYc5c+Dff10dhVJKpeK0pGCM+RM4lcZTo4EngZQLOXQBvjDWSiBYRMo5KzaXu3ABevaEe+91dSRKqRQ8PT0JCwujfv36NGrUiPfff5/ExMQMX7Nv3z6+/vrrPIrQ+fK0T0FEbgcOGmM2XvZUBeBAiscRjn0F0/r1kJAAixeDzshWym34+/uzYcMGtm7dysKFC1mwYAEvv/xyhq/RpJBNIhIAPAu8kNbTaexLc0k4ERkqImtEZM3x48dzM8S8s3q1/RkQAO+849pYlFJpKl26NBMnTmTcuHEYY9i3bx/XXnstTZo0oUmTJixfvhyAp59+mqVLlxIWFsbo0aPTPS6/yMshqdWBqsBGR+dHRWCdiDTHtgwqpTi2InAorZMYYyYCEwHCw8Pz51qif/8NlSpB797w7ruwZw9Uq+bqqJRyGy//sJVth87m6jnrlS/Gi53rZ+k11apVIzExkWPHjlG6dGkWLlyIn58fu3btonfv3qxZs4ZRo0bx7rvvMn/+fACioqLSPC6/yLOWgjFmszGmtDEm1BgTik0ETYwxR4B5wD2OUUgtgUhjzOG8ii3PrV4NzZrBww+DpyeMHu3qiJRS6Uhaxz4uLo57772XBg0a0LNnT7Zt25bm8Zk9zl05raUgItOBdkCIiEQALxpjPkvn8AXAbcBuIAoY6Ky4smTiRHt7x/GPgnLl4JtvoEIOujtOnrSjjoYMgfLloW9f+OwzePFFCAnJnbiVyuey+o3eWfbs2YOnpyelS5fm5ZdfpkyZMmzcuJHExET8/PzSfM3o0aMzdZy7clpSMMb0vsrzoSl+N8ADzoolW/bvt9/k69SB+o5/oHPnwl13we+/Q8pp4z/9BNu3X3p87bW2JZCWpGZk8+b25+OPw5QpMH48vJBWd4tSyhWOHz/Offfdx4MPPoiIEBkZScWKFfHw8GDq1KkkJCQAEBgYyLlz55Jfl95x+YWWuUjPk0+CCMybZ+//A0yfDn36wDPP2BZEYiI89xy8+Wbq1wYEwM6dabcokjqZmza1P+vXhxtugNmzNSko5WLR0dGEhYURFxeHl5cXd999NyNHjgRg+PDhdO/enVmzZtG+fXuKFCkCQMOGDfHy8qJRo0YMGDAg3ePyDWNMvt2aNm1qnGLpUmPAmBdfvPK54cPtc998Y0zv3vb3oUONOX3amMhIY7ZuNcbHx5i770773J07G1OnTup9jz5qjL+/MQkJuf5WlMovtm3b5uoQCqS0/q7AGpPO56rWPrpcYiI88ghUrGhbC5d7/30ID7e3kaZPt62Ejz+G4GAoVgzq1YORI+HLL2HVqtSvNeZSJ3NKNWtCdDQcLrh960qp/EGTwuWmToW1a+Gtt+xtoMv5+sKsWdC2LXz9NTz9tL3NlNIzz0DZsja5mBSjZg8ehCNH0k4KALt25e57UUqpLNKkkNK5c/YDvWVLO4cgPaGhsHRp+scEBsIbb8DKlTZxJPn7b/szqZM5SY0a9qcmBaWUi2lSSOmNN+w3+Q8/vPLbf1b17287k598Ev5z1P5bvRq8vKBRo9THVqoEPj6aFJRSLqdJIcnevba/4O67r/wmnx0eHvDJJ7b4XYsWsG6dbSk0bAiXj1v29ITq1WH37pxfVymlckCTQpInnrDf4i8fXpoTTZvCsmW2FXDddbBiRfoJp0YNbSkopVxOkwLAH3/At9/aTuOczFZOS/36tm+hdm07wii9SW01a9qWwlXK9CqlnOf111+nfv36NGzYkLCwMFatWsWQIUNyrVRFaGgoJ06cyPCYN954I8vnnTJlCg8++GB2w0pFJ68lJNhRQpUr29nFzlCunE08X39tJ7+lpWZNiImBQ4fscFilVJ5asWIF8+fPZ926dfj6+nLixAliY2OZNGlSnsbxxhtv8Mwzz+TpNVPSlsLkybBhgx2C6u/vvOsULQpDh17Zn5BEh6Uq5VKHDx8mJCQEX19fAEJCQihfvjzt2rVLrnJatGhRnnrqKZo2bcqNN97I33//Tbt27ahWrRrz5s0DrvzW3qlTJ5YsWXLF9bp27UrTpk2pX78+EydOBGwZ7qRZ1X379gVg2rRpNG/enLCwMIYNG5ZcNuPzzz+nVq1aXH/99SxbtizX/g6Fu6UQGQnPPgtt2kCvXq6NJeWw1PbtXRuLUq72yCP2y1puCguDDz5I9+mbbrqJV155hVq1anHjjTfSq1cvrr/++lTHXLhwgXbt2vHWW2/RrVs3nnvuORYuXMi2bdvo378/t99+e6bDmTx5MiVKlCA6OppmzZrRvXt3Ro0axbhx49jgeO/bt29nxowZLFu2DG9vb4YPH85XX31Fx44defHFF1m7di1BQUG0b9+exo0bZ+/vcpnCnRRefx2OH4cFC3I+BDWnKlWyE+N0BJJSLlG0aFHWrl3L0qVLWbx4Mb169WLUqFGpjvHx8eGWW+zS8w0aNMDX1xdvb28aNGjAvn37snS9MWPGMGfOHAAOHDjArl27KFmyZKpjFi1axNq1a2nm6IuMjo6mdOnSrFq1inbt2lGqVCkAevXqxc6dO7Pztq9QeJPC7t32W8OAAbZshat5eNhhqRndPlq/3rYoAgPzLi6lXCGDb/TO5OnpSbt27WjXrh0NGjRg6tSpqZ739vbGsUgYHh4eybeaPDw8iI+PB8DLyyvVus4xMTFXXGfJkiX89ttvrFixgoCAANq1a5fmccYY+vfvz5uXjYqcO3duchy5rfD2KTzxhP1mno2efqfJaFjqqVN2vkPXrrZzXCmVq3bs2MGuFP//bdiwgSpVqmT5PKGhoWzYsIHExEQOHDjA30mVDFKIjIykePHiBAQE8M8//7By5crk57y9vYmLiwOgQ4cOzJ49m2PHjgFw6tQp9u/fT4sWLViyZAknT54kLi6OWbNmZTnO9BTOpPD773ZthGeesSOD3EXNmnYBnrSGpf76K8TF2dhfeinPQ1OqoDt//jz9+/enXr16NGzYkG3btvFSNv5fa9OmDVWrVqVBgwY8/vjjNGnS5IpjbrnlFuLj42nYsCHPP/88LVu2TH5u6NChNGzYkL59+1KvXj1ee+01brrpJho2bEjHjh05fPgw5cqV46WXXqJVq1bceOONaV4ju8SY/LnMMdg1mrO19umWLTBqFEyalP5oIFf45BO47z5bFqNSpdTP3XMP/PgjdOkCn39u+0FuvdU1cSrlBNu3b6du3bquDqPASevvKiJrjTFp3jcvnC2Fa66BadPcKyFA+oXxEhPh55/hllvgo49s7aR+/ezqcEoplYsKZ1JwV+nNVVi3zo6SuvVWO5di9myIj4cH3GsFU6VU/qdJwZ1UrJj2sNSkIbM332wf16gB998Pv/xiO6CVKiDy8+1sd5Sdv6fTkoKITBaRYyKyJcW+d0TkHxHZJCJzRCQ4xXP/E5HdIrJDRG52VlxuLb1hqT/9ZGsmOcYkA9Cjh20tOGZR5khCgu3gVsqF/Pz8OHnypCaGXGKM4eTJk/hl8Ta5M+cpTAHGAV+k2LcQ+J8xJl5E3gL+BzwlIvWAu4D6QHngNxGpZYxxytjLP3Ye55UftlLUz5uivp4U8fGiqK8XRRxbUV/P5N8DU+wP9HM87+NFEV9PvDydkFNr1kydFE6csMt6vvBC6uOaNoUqVeytpAEDsn+98+ftbO6ffrLrPTRtmv1zZUN8QiKR0XGciY7jTFQsZ6LiOB0VR3RsPHEJhoREQ4CvJzfUKU25ICeWIVEuV7FiRSIiIjh+/LirQykw/Pz8qJjFWmpOSwrGmD9FJPSyfb+meLgS6OH4vQvwjTHmIrBXRHYDzYEVzoitqK8ndcoW49zFeC5cjOfk+SjOO34/f9F+GGWGr5dHqmRiE4inI7GkTjJFfb0p4uuZvL9oimOK+nrh5+1hJ6PUrGk/oPfuhapV7VBUY+C221JfXMS2FsaMseU6goKuHvBjj9mRTQ8+aEt5HzkCnTrZcgIeHjBzZraTQkKi4azjw/10VCyRUfbn6ag4IqNiHfvtB3+k45gzUXGci4nP9DUaVw7m/xqUo1/LKvh5e2YrTuW+vL29qVq1qqvDKPRcOaN5EDDD8XsFbJJIEuHYdwURGQoMBahcuXK2Lty0SgmaVimR7vMX4xO4cDEhOUlcuBifnEDsvoQUv8enSignzsey/+SlJHMhNnONHQ+BIr5eNLpYiwmePsQ3Cue9Ee/SafEsrgkM5p2DfgT8/E+qZFKhcTtaxb3Hvs+/Ia53H5uEEuIouvR3PDt3tov3JDlwAEaPtslk9mxbB+bUKTh5En74AcaOhdmzSXzjTc5ejHd8Y49N9Q3+jONDPekDPvnD/kIsZzP4cBeBIH9vgv29CQrwoUQRH6qFFCE4wIfgAG+KB/jY5x2/Fw/wIcDXE28PDzw9hSORMfyy9QgLNh/mtR+38+26g3zUpzHVShXN9H9zpVTmuCQpiMizQDzwVdKuNA5L8+u6MWYiMBHsPAVnxOfr5Ymvlyclivjk+FyJiYYLsfFcuJhwKVFclmQu/Z7AhYtleb/C5wx/dwTPvTOcBPFgad3WzN10hAsX44lPvPSWxRiWBYawc/wUhh6pCMbw3o/v033rYv7X+VEWNr81uSXSb9E0ehvDq29+Q41/t9Dup6/xj4rhg8fHs2FvEG0Cr+GJPT/Tecg4tpapnu77KebnRfEiPskf8KEhRZJ/Lx5gP9iD/e2HfbBjX6CfN54e2Z+SX6N0UWqUrsED7Wuw+J9jPDpzA53H/sWo7g3p3Kh8ts+rlLpSnicFEekPdAI6mEs9ShFAytlaFYFDeR2bM3h4CIF+9oMx8xpB/7bQuTOsWcMtTwzklj43YYzhYnxicgI5fzEej+Pd6fDVFMbfXpOQ+XNovnUx8T6+jNw4D4/+/TkXm0hUdCztl89nU43G/OFdhkU1SzGu6vV4SSJBRf0J8vPiRIdbSPh2NC/EbGNLp87JH/BB/kkf9vbbfE4+3HND+zqlWTDiWh78eh0PTV/PzqPnGNmxltPqwChV6BhjnLYBocCWFI9vAbYBpS47rj6wEfAFqgJ7AM+rnb9p06amQDt/3pjp042Ji0v/mD//NAaMefJJY3x9jbn5ZmOmTrX75s+3xyxZYh9/8UXG1+vQwZiaNY1JTMy99+AksfEJ5slZG02Vp+abd3/5xyTmg5iVchfAGpPO56ozh6ROx3YU1xaRCBEZjB2NFAgsFJENIvKxIzFtBWY6EsbPwAPGSSOP8pUiReCuu+za0elp3RrKloW337ZDVqdNg969bZmMd96xx0yebCurdu+e8fV69LAjn7Zsyfg4N+Dt6cGbdzSgV3glxv6+mw9+08WJlMoNTksKxpjexphyxhhvY0xFY8xnxpgaxphKxpgwx3ZfiuNfN8ZUN8bUNsb85Ky4ChxPT+jZ0yaOmTMhJAS8veHRR+0SoL/9BrNm2UQREJDxubp1u9QRnQ94eAhv3tGAO8Mr8uGiXYxdpIlBqZwqnAXxCprz5+HgQahd+9K+c+dsa8HT044yWrnSlt6+mnbt7NyIzZttRdYJE6BxY1tR1k3v2ycmGp6YvYlv10Xw8u316d861NUhKeXWtCBeQVe0aOqEAPZ20fDhNiHUqwfNm2fuXD16wNatULcu3HijLcT33HO2SmtsbO7Hngs8PIS3ujegY70yvDhvK99vOOjqkJTKtzQpFGQjRkCxYrZwXma/5XfvbvsyfH1tX8SJE/Dqq7av4uab4fRp58acTV6eHozt3ZiW1Urw2MyNLN5xzNUhKZUv6e2jgi4qylZWzcqtn7ReM20aDBpki/H99JMtseGGzsXE0fvTlew+dp4vB7egWWj6kxSVKqz09lFhFhCQ9b6AtF7Tr5+tynroELRsCWvX5l6MuSjQz5upA5tTPtifQZ+vZsvBSFeHpFS+oklBZV779rB8Ofj42NpJM2fCzp12i4hwdXTJShb1ZdrgFhTz96b/5L/59/h5V4ekVL6hSUFlTb16diRTnTq2umrt2narVAk++8zV0SUrH+zPtCEtEIG7J60i4nSUq0NSKl/QPgWVPRcu2L6FpBFJ48bZFsPu3RAcnPFr89C2Q2e5a+IKggK8mTmslZbfVgrtU1DOUKSIHb7ap4/dxo+3w19fecXVkaVSr3wxvhzcgjMX4ujz6SqOno1xdUhKuTVNCip3hIXBkCG2BPeOHa6OJpVGlYKZMqgZx87G0OfTlRw/d9HVISnltjQpqNzz2mt2KOtjj7k6kis0rVKCyQOacehMDL0/XckxbTEolSZNCir3lC4Nzz8PP/5oZ0K7mRbVSvL5wGYcOhPNXRNXciRSE4NSl9OkoHLXiBFQqxYMHgzH3G9WcctqJZk6qDlHz8bQa+IKDp2JdnVISrkVTQoqd/n6wowZttO5Tx9IcL8K6M1CS/DlkBacOh9Lz49XsO/EBVeHpJTb0KSgcl9YmB2iumgRvPyy3RcdDZMm2TpMy5ZByqHQGzbAI4/ABx/A2bN5EmKTysX5+t6WRMXG0/OTFew4ci5PrquUu9N5Cso5jLG1kqZOtbeS5syBkyftbOjYWAgPt2s8/PADLFli14CIi7MVXwcNsutBhIY6PcxdR8/R77NVXIxPZOrA5jSq5D5zLJRyFp2noPKeCHz0ETRoYGc6X3edXfTn1Cm7RsP583aU0t69doW4o0dh9Wro2tU+HxYGixc7PcyaZQKZfV9rivl50+fTlfy164TTr6mUO9OWgnKuyEibACpUSL0/MdHOgK5R48rlRvfuhU6d7NKgkyfbYnxOdvRsTHKdpNG9wujUsLzTr6mUq2hLQblOUNCVCQHAw8PWT0pr/emqVW2/Q9u2cPfd8MYbTg+zTDE/ZgxrRVilYB6avp4vVuxz+jWVckdOSwoiMllEjonIlhT7SojIQhHZ5fhZ3LFfRGSMiOwWkU0i0sRZcal8IjjYznXo2xeefdau5+BkQf7efDm4BR3qlOGF77fyzi//kJ9b0kplhzNbClOAWy7b9zSwyBhTE1jkeAxwK1DTsQ0FJjgxLpVf+PjAlCm2P2LYMLtMqJP5eXvycb8m9G5eiY8W/8tjszYSG5/o9Osq5S6clhSMMX8Cpy7b3QWY6vh9KtA1xf4vjLUSCBaRcs6KTeUjXl7wzTd2zekePWz/xOVOnIC334bZs3Pnkp4evNGtASM71uK7dQcZPHU152LicuXcSrm7NG7oOlUZY8xhAGPMYREp7dhfATiQ4rgIx77Dl59ARIZiWxNUrlzZudEq91CuHEyfDjfeaIe3PvGE3R8VBV9+aW8txcTYBFK+PLRuneNLiggjOtSkbJAf//tuMz0/XsHnA5tp6W1V4LlLR3Na60WmeTPXGDPRGBNujAkvVaqUk8NSbqN9e1twb+ZMaNbMbtdfD199Bf37w4oVdt3oO++E48dz7bJ3hldi8oBmRJyOpttHy9l2KG8m1ynlKnndUjgqIuUcrYRyQFJxnAigUorjKgKH8jg25e6efhquvRbOnLGPRex60SVL2sezZ9vHffvaBYDOnbNDWr/99tJiQBmpXdvOrQgKSrX7+lqlmDmsFYOmrKbnx8sZ26cxN9Qpk8tvTin3kNdJYR7QHxjl+Pl9iv0Pisg3QAsgMuk2k1LJROww1fQklde4917o0AHWrLErxDVrZiu4ZiQx0dZs2rgRFiywy4umUK98MeY+0IZBU1YzaMoaBretypO31MbXyzMX3phS7sNpk9dEZDrQDggBjgIvAnOBmUBl4D+gpzHmlIgIMA47WikKGGiMueqsNJ28pq5gjF3s58svbRmNhx+GJpkc4fzbb9C9uy218eOPNslcJiYugTcXbGfqiv3UL1+MMb0bU71U0Vx+E0o5V0aT13RGsyp4EhNtJ3TRbHxYb94Mt91my3Hcey889BBUr37FYb9uPcKT324iKjaBYddV4/521QnwyeuGt1LZo0lBqaw4dAiefNJ2asfH2yRRs+YVh12IjWf13lP8e/wCRXw9adCsHvXeeh7x8XFB0EplniYFpbLj8GFbnO/zzzMs6R1vDLGxCQRcjGLujX0InfwRYVptVbkxTQpKOVlComFPr/7UnP0l997xHAE97+Dxm2pTqUSAq0NT6goZJQW9CapULvD0EGpO+5SEvdsY+8uH/F/ZqnTYfIS7W1XhgfY1KFFEbymp/MFdJq8plf/5+uI5ezZ+Pl789NcYetQrwefL9nL924v58LddnNVSGSof0KSgVG4KDYUvvsBn00bemPgkC/tfQ6vqJRn9206ufWsx437fpXWUlFvTPgWlnGH6dBgwAKpVgwUL2Oxdgg9+28mif45RzM+Le1qFMqBNKCFFfV0dqSqEctTRLCJFgGhjTKKI1ALqAD8ZY1z+dUeTgnJrf/5plxf19ob77gMPD45ExrB63yl2Hz+Ppwj1yxcjrFIwJbOSHEqWhD59oEQJ58WuCrScJoW1wLVAcWAlsAaIMsb0ze1As0qTgnJ7//xjE8OOHbl7Xn9/uOceuP9+W0UWbBmQkBD7U6kM5HT0kRhjokRkMDDWGPO2iKzP3RCVKqDq1IHt29N9+uT5i8xcE8G0lfs5eCaaEkV86BJWnh5NK1K/fFDaL9q8GcaMsQsQffJJ6ueuuQYeecS2JPy1zLfKusy0FNYDw4HRwGBjzFYR2WyMaZAXAWZEWwqqoEhINPy58ziz10awcNtRYhMSqVM2kG6NK9AlrAJlg/yufNHx4zBvHly8aB9HRdm1JTZutC2GYcNg+HC7xoRSKeT09tF1wOPAMmPMWyJSDXjEGDMi90PNGk0KqiA6ExXLDxsP8d36g6z/7wwitnz34LZVaVsjBMno9pAx8Mcf8MEHNmF4eto1Jh57LPOFAVWBl9Ok0NMYM+tq+1xBk4Iq6PaeuMCcdRFMX32A4+cuUrtMIMOur0a3xhUyTg4A//5rS4l/9pltRXz8sa0gqwq9nCaFdcaYJlfb5wqaFFRhcTE+gR82HmbS0j38c+Qc7WqX4u0eDSkdmMZtpctFRsJdd8HPP8Ozz8Krr2pndCGXraQgIrcCtwF3AjNSPFUMqGeMaZ7bgWaVJgVV2Bhj+HLlfl7/cTtFfL14844G3Fy/7NVfGBdn+xcmTbKd0KNHX33hIVVgZZQUMprRfAg7/DQGWJtimwfcnNtBKqWuTkS4p1UoP45oS/lgP4Z9uZZpK/df/YXe3jBxIrz+Onz9NVSuDIMHw6ZNzg9a5SuZuX3k7Q4T1dKiLQVVmMXGJ3L/tLUs3nGMj/s15abMtBjAzp0YMwamTrV9DSNGwPvv205pVShkt6WQpLmILBSRnSKyR0T2isieXI5RKZVFPl4ejO3TmAYVgxnxzXrW/Xc6cy+sUwfGj4cDB+DBB22C6N7dJghV6GUmKXwGvA+0BZoB4Y6fSikXC/Dx4rP+4ZQp5sfgKavZe+JC5l9cogSMHQsffmiHr7ZvD8eOOS9YlS9kJilEGmN+MsYcM8acTNpyclEReVREtorIFhGZLiJ+IlJVRFaJyC4RmSEiWoBeqUwIKerL1IF23MfD36wnLiExaycYMQK++87OlL75ZoiOdkKUKr/ITFJYLCLviEgrEWmStGX3giJSARgBhBtjrgE8gbuAt4DRxpiawGlgcHavoVRhExpShNe7NWBTRCTjF/+b9RN07WrXpN6wAR5+OPcDVPlGZpJCC+wtozeA9xzbuzm8rhfgLyJeQABwGLgBmO14firQNYfXUKpQua1BObqElWfs77vYcjAy6yfo1An+9z/49FP44ovcD1DlCy5ZT0FEHgZeB6KBX4GHgZXGmBqO5ythy3Nfk8ZrhwJDASpXrtx0//5MDMdTqpCIjIrjpg/+IMjfm3kPtsXPO4sjiuLjoWNHWLXKbg1cXuJMOUGORh+JyAtpbTkIpjjQBagKlAeKALemcWia2coYM9EYE26MCS9VqlR2w1CqQAoK8GZU94bsPHqe0Qt3Zv0EXl52gaCgIOjb19ZSUoVKZm4fXUixJWA/wENzcM0bgb3GmOOO+Q/fAa2BYMftJICK2MlzSqksal+7NL2bV2bi0j2s3X8q6ycoWxaef952PO/MRmJR+dpVk4Ix5r0U2+tAO6BCDq75H9BSRALEVvTqAGwDFgM9HMf0B77PwTWUKtSe/b+6VAj257GZG4mKjc/6CW67zf5csCB3A1NuLzMthcsFANWye0FjzCpsh/I6YLMjhonAU8BIEdkNlMTOj1BKZUNRXy/e6dGIfSejePvnbKz6FhoKdevCTz/lemzKvV115TUR2cyl+/ueQCnglZxc1BjzIvDiZbv3AC4vsqdUQdGqekkGtA5lyvJ93FS/DK2rh2TtBLfdZie3nT8PRYs6J0jldjLTUugEdHZsNwHljTHjnBqVUipXPHVLHaqGFOGJWZs4ExWbtRffeivExsLixc4JTrmlzPQp7AeCsUmhG1DP2UEppXKHv48no3uFcexcDCNnbiQxMQujidq2tS0EvVzh11IAAB2ZSURBVIVUqGRmSOrDwFdAacf2lYg85OzAlFK5I6xSMM93qsfv/xxjwh9ZmO3s6wsdOtjOZh2aWmhk5vbRYKCFMeYFY8wLQEvgXueGpZTKTXe3rELnRuV579cdLP/3ROZfeOutsH+/LbetCoXMJAXBzk9IkuDYp5TKJ0SEUXc0oGpIEUZMX8/BM5ksenerY16p3kIqNDKTFD4HVonISyLyErASHS6qVL5TxNeLj/s15WJ8Ind/toqT5y9e/UWVK0P9+jpfoRDJTEfz+8BA4BS2eulAY8wHzg5MKZX7apYJ5LP+zTh4OpqBU1Zz/mImJrbddhv8+Sf89pv2LRQCmVmOsyWw1RhzzvE4EKjnmITmUrocp1LZ89u2owybtpaW1UoweUAzfL0yKJy3fbtdgOfoUbjmGrtaW82al55v0gSCg50ftMo1OV2OcwJwPsXjC459Sql86sZ6ZXire0OW7T7J0C/WZlwKo25d2LcPPv/cruN83312VFLS1q9fnsWtnC9THc0mRXPCGJNIJmZCK6XcW4+mFRl1RwOW7jrO3Z/9TWRUXPoH+/nBgAGwfr3d/vjDbiNGwI8/wpYteRa3cq7MJIU9IjJCRLwd28PYkhRKqXzuruaV+ahPEzZHRHLnJys4ejYm4xeIQFgYXHed3V58EQIC4N2crrul3EVmksJ92NLWB4EI7EpsQ50ZlFIq79zaoByTBzTjwOkoun20jK2HsrBqW4kSMGQIfPUVREQ4L0iVZzIz+uiYMeYuY0xpY0wZY0wfY8yxvAhOKZU32tYMYeawVhigx4QV/LzlcOZf/OijdlTShx86LT6Vd7JTOlspVQBdUyGI7x9oQ+2ygdw3bR1jF+0iU8v1hobCnXfCJ59AZDbWhlZuRZOCUipZ6WJ+fDO0JV3DyvPewp0M+3It52Iy6IBO8sQTcO6cTQwqX9OkoJRKxc/bVlZ9vlM9Fv1zjC4fLWPX0XMZv6hxY7jxRnjtNVi0KG8CVU6RmSqpJUVkrIisE5G1IvKhiJTMi+CUUq4hIgxuW5WvhrTgbHQcXT9axo+brtLP8Pnn9lbSLbfA1Kl5EqfKfZlpKXwDHAO6Y9dQPg7McGZQSin30LJaSeY/dC21ygbywNfreP3HbcQnJKZ9cMWKsHQptGtn5zQ8/zzEXGWIq3I7mUkKJYwxrxpj9jq217CL7iilCoGyQX7MGNqKe1pV4dOle+k7aRXHz6VTTC8oyE5mGzDA3kqqXBleeAEOZ2E0k3KpzCSFxSJyl4h4OLY7gR9zclERCRaR2SLyj4hsF5FWIlJCRBaKyC7Hz+I5uYZSKvf4eHnwSpdrGN2rERsjztBp7FLW/Xc6nYN9YPJk27fQqpVNDlWqwDhdxTc/yExSGAZ8DVwEYrG3k0aKyDkROZvN634I/GyMqQM0ArYDTwOLjDE1gUWOx0opN9KtcUW+u78Nvl6e9PpkBdNW7k972KoI3HADfP897NwJN98MDz0Ejz8OiencflJu4apVUnP9giLFgI1AtZQ1lURkB9DOGHNYRMoBS4wxtTM6l1ZJVco1IqPieHjGepbsOE7PphV5tes1+HlnUGk1IQEeecS2Fnr0gC++AH//vAtYpZJRldR0C9uJSJOMTmqMWZfNeKphO6s/F5FGwFrgYaCMMeaw49yHRaR0Ns+vlHKyoABvJvdvxge/7WTM77vZefQcH9/dlHJB6XzQe3rCmDFQtaptLcTFwZw5tkWh3Eq6LQURSQS2Yj/AIfUSnMYYc0O2LigSjl29rY0xZpWIfAicBR4yxgSnOO60MeaKfgURGYqj9lLlypWb7t+/PzthKKVyyS9bjzByxgb8fTwZ37cpzauWyPgFH3xgS2O8/bad9KbyXHbXU3gMiASisUtydjbGtHds2UoIDhFARIpFemYDTYCjjttGOH6mWV/JGDPRGBNujAkvVapUDsJQSuWGm+uX5fsH21DMz5s+n65k2sqrfFF7+GF7C+l//7NDWJVbSTcpGGNGG2PaAg8ClYBFIjJTRMJyckFjzBHggIgk9Rd0ALYB84D+jn39ge9zch2lVN6pUTqQOQ+0oW3NEJ6bu4Xn5m4mLr35DCLw2WdQrRr06mVXdFNuIzNVUvdiP6B/BZoDtXLhug8BX4nIJiAMeAMYBXQUkV1AR8djpVQ+EeTvzWf9mzHs+mpMW/kf/Sat4tSF2LQPLlYMZs+G06ehd2/bEa3cQkZ9CtWAu4AuwAHsUNT5xhi3maKoo4+Uck/fbzjIE7M3UaaYL5PuaUbtsoFpHzhlCgwcCM8+a+czqDyRUZ/C1TqaN2FbCWeBVAcaY97P5TizTJOCUu5rw4Ez3PvFGqIuxjOmd2M61C2T9oGDB9vJbgsWwK235m2QhVR2O5pfAeYAiUBRIPCyTSml0hVWKZh5D7ahWqmiDPliDZ/+uSftiW7jxkHDhtCvH/z3X94HqlLJ88lruUlbCkq5v+jYBB6btYEFm49wV7NKvNLlGny8Lvs+umsXhIdDnTp2RJKPj2uCLSSyO3nthQzOaYwxr+Y4MqVUgefv48m43k0YXWonY3/fzf6TUUzo14TggBQf/DVr2tLb3bvbyW1jxrgu4EIuo9tHF9LYAAYDTzk5LqVUAeLhITx2U23ev7MRa/efptv45ew7cSH1QXfcASNHwtixMEOr87tKpm4fiUggthTFYGAm8J4xJs3JZXlJbx8plf+s3neKoV/Y/28n3hNOs9AUM6Dj4ux6DJs2wZo1UDvD8mcqm7Lb0YyjnPVr2FFIXkATY8xT7pAQlFL5U7PQEswZ3obiAT70/XQVc9cfvPSkt7dtJfj52VtJFy6kfyLlFOkmBRF5B1gNnAMaGGNeMsakU0BdKaUyLzSkCN8Nb03jysE8MmMDYxftujQyqWJF+Ppr2LYN7r8f8vFgmPzoarWPygPPAYdE5Kxjy8k6CkopBUBwgA9fDG5Ot8YVeG/hTp6cvelSaYyOHeHFF+HLL2HSJNcGWsikO/rIGJOZBXiUUirbfL08ef/ORlQqEcCYRbs4FBnNhH5NKebnDc89B8uX28V5wsOhcWNXh1so6Ae/UsqlRISRHWvxbs9GrNpzih4TlnPwTLRdg2HaNChVylZVXbUK1q2zW2Skq8MusDQpKKXcQo+mFZk6qDmHz8TQ7aNlbDkYaRPCzJl2pnPLltC0qd1q14a1a10dcoGkSUEp5Tba1Ahh9v2t8fIQ7vxkBYt3HINWrewQ1e+/t1vS6KTrroP5810dcoGjSUEp5VZql7VrM1QNKcKQqWv45u//oG5duP12u915J6xcafd16QIffgjx8a4Ou8DQpKCUcjtlivkxY1gr2tYI4envNvPerztSF9MrWxaWLIHbboNHHrFrP48aBSdPuizmgkKTglLKLRX19WJS/3B6hVdi7O+7eWzWRmLjU6zmVrSovZ00bx7UqmWX96xSRZf4zCFNCkopt+Xt6cGo7g0Y2bEW3607yKApqzkbE3fpAA8P6NwZFi2CjRuhfHld4jOHNCkopdyaiDCiQ03e6dGQlXtOcufHKzgSmcYCkA0bXlris08fXeIzmzQpKKXyhZ7hlZg8oBkRp6PpNn4ZO46cu/Kghg1hwgT4/Xd46aU8j7Eg0KSglMo3rqtVipnDWpFoDD0+Xs7yf09cedCAATBokF3z+fff8zzG/M5lSUFEPEVkvYjMdzyuKiKrRGSXiMwQEV16SSl1hXrli/Hd8DaUC/Kj/+S/U1dZTTJuHISEwJQpeR5ffufKlsLDwPYUj98CRhtjagKnsWs3KKXUFSoE+zPrvtY0rVKcR2ZsYPyS3amHrPr7Q/v2sHixVlnNIpckBRGpCPwfMMnxWIAbgNmOQ6YCXV0Rm1Iqfwjy92bqoOZ0blSet3/ewfPfbyEhMUUCaN8eIiLg339dF2Q+5KqWwgfAk0DSoOOSwBljTNK0xAigQlovFJGhIrJGRNYcP37c+ZEqpdyWr5cnH/YKY9j11Zi28j/um7aW6FjHqKP27e1P7VfIkjxPCiLSCThmjElZzUrSODTNNp8xZqIxJtwYE16qVCmnxKiUyj88PIT/3VqXl2+vz2/bj9L705WcPH/RFs0rV87eQlKZ5oqWQhvgdhHZB3yDvW30ARAsIknrO1QEDrkgNqVUPtW/dSgT+jZl++GzdJ+wnH0no7RfIRvyPCkYY/5njKlojAkF7gJ+N8b0BRYDPRyH9Qe+z+vYlFL52y3XlOXre1sSGR3HHROWs79hczu7+Z9/XB1avuFO8xSeAkaKyG5sH8NnLo5HKZUPNa1SnO+GtyHQz4sh/wXanXoLKdNcmhSMMUuMMZ0cv+8xxjQ3xtQwxvQ0xlx0ZWxKqfyrakgRvr2/NQF1a3GwWCn2z/rB1SHlG+7UUlBKqVwTUtSXb4a2Yl+D5gSuWs6b87eSmKh9C1ejSUEpVWD5+3jSanAPSkSf5Y/vFjPim/XExGmhvIxoUlBKFWgeHW4A4PkiR5m/6TB3f7aKM1GxLo7KfWlSUEoVbJUrQ/XqtNmzjrG9G7PxQCR3TFjOfyejXB2ZW9KkoJQq+O64A375hc5lPZk2pAUnz8dyx4RlbDxwxtWRuR1NCkqpgm/gQLvozrRpNK9agm/vb42/jye9Jq7g161HXB2dW9GkoJQq+OrWhVatYPJkMIYapYvy3f1tqF0mkGHT1jJl2V5XR+g2NCkopQqHQYNg2zb4+28ASgXaIas31i3DSz9s45UftqWuslpIaVJQShUOd94JAQG2teDg7+PJx/2aMrBNKJOX7WX4VymqrBZSmhSUUoVDsWLQsydMnw5Rl0YeeXoIL3auzwud6vHrtqPc9elKjp8rvAUVNCkopQqPgQPh3Dn49tsrnhrUtiqf9GvKjiNn6TZ+GbuPnXNBgK6nSUEpVXhcdx1Ur27XcD527Iqnb6pflpnDWhETl0i38ctZvvuEC4J0LU0KSqnCQwSefNJ2NleubDufN25MdUjDisHMfaA15YL8uGfy38xcc8BFwbqGJgWlVOEydChs324TwowZEBZmF+P5/ns7lwGoWDyA2fe3plX1kjw5exPv/PJPoSmmp0lBKVX41KkD48dDRAS8/Tbs2QNdu0KtWrBmDQDF/LyZPKAZvZtX4qPF//JQISmmp0lBKVV4FS8OTzwB//4Ls2ZBTAwMGZLcYvD29OCNbg14+tY6/LjpML0/XcmJ8wV7ZJImBaWU8vKCHj1g9Gjbx5BiLoOIcN/11fm4XxO2Hz5L14+WsetowR2ZpElBKaWS9OwJbdvCs89CZGSqp265phwzhrbiYnwid4xfzp87j7soSOfSpKCUUklE4IMP4MQJeO21K55uVCmYuQ+0oUJxfwZOWc2XK/e7IEjnyvOkICKVRGSxiGwXka0i8rBjfwkRWSgiuxw/i+d1bEopRdOmMGAAfPgh7Np1xdMVgv2ZfX9r2tUqxfNzt/DSvK3EJyTmfZxO4oqWQjzwmDGmLtASeEBE6gFPA4uMMTWBRY7HSimV9954A3x94eWX03y6qK8XE+8JZ3DbqkxZvo8hX6zhbExcHgfpHHmeFIwxh40x6xy/nwO2AxWALsBUx2FTga55HZtSSgFQtiz07Qtz50J0dJqHeHoIz3eqxxvdGvDXrhN0H18wVnNzaZ+CiIQCjYFVQBljzGGwiQMonc5rhorIGhFZc/x4wezoUUq5gR494MIF+OWXDA/r06IyXwxqzrFzF+k6fhl/7z2VRwE6h8uSgogUBb4FHjHGnM3s64wxE40x4caY8FKlSjkvQKVU4Xb99VCyJMyefdVDW9cIYe4DbQj296bvpJX5ujSGS5KCiHhjE8JXxpjvHLuPikg5x/PlgCurVSmlVF7x9oYuXWDePLh49QlrVUOKMGd4G1pUtaUxXpufPxftccXoIwE+A7YbY95P8dQ8oL/j9/7A93kdm1JKpdKjhy21vXBhpg4PCvBmysBm9G9VhUl/7WXI1NWcy2cd0K5oKbQB7gZuEJENju02YBTQUUR2AR0dj5VSynU6dICgoEzdQkri5enBy12u4bWu17B01wm6jV/OvhMXnBhk7hJj8l/zJkl4eLhZ4yhepZRSTtG/v72FdPQo+Phk6aXL/z3B8K/WYQxM6NuE1jVCnBRk1ojIWmNMeFrP6YxmpZTKSI8ecOYMLF6c5Ze2rh7C9w+0oXSgL3dP/pspy/bi7l/ENSkopVRGOnaEwMAs3UJKqUrJInw3vDXta5fipR+28fS3m7kY774luDUpKKVURvz84Pbb7YI8e/dm6xSBft5MvDucB9vXYMaaA/T5dBXHz7lnCW5NCkopdTWvvgqenvZWUkxMtk7h4SE8fnNtxvVpzNZDkXQe+xebIs7kcqA5p0lBKaWupmpVmDoV1q2DRx/N0ak6NSzPt/e3xtND6PHxCr5bF5FLQeYOTQpKKZUZt98OTz4JH38MX32Vo1PVLx/EvAfb0KRyMCNnbuTV+dvcptKqJgWllMqs11+H666De++Fxx+HffuyfaqSRX35cnALBrQO5bO/9nLP5L85dSE292LNJk0KSimVWV5eMHOmbTV88AFUrw7du8Phw9k6nbenBy/dXp93ejRkzf7TdB77F1sORl79hU6kSUEppbKiTBn45hs7EumJJ2DBAnjooRydsmd4JWYNa0WiMXSfsJw5613Xz6BJQSmlsqNSJRg1Cp55Br79Fv74I0ena1QpmB8eaktYpWAenbGRl+ZtJc4F/QyaFJRSKicefxwqV4ZHHoGEnE1KCynqy7QhLRjUxq7o1vfTVRw7l70hsNmlSUEppXLC3x/efhs2bIDPP8/x6bw9PXihcz0+6BXGpoNn6Dz2L9buz7uFezQpKKVUTt15J7RpA88+C2czvWZYhro2rsB397fB18uTuyau5IsV+/KkbpImBaWUyikROxrp2DF47rlcO2298sX44cG2XFuzFC98v5WRMzcSHevcukmaFJRSKjeEh9tRSGPHwpw5uXbaoABvJt0TzsiOtZi74SDdxi9jrxPXZ9CkoJRSueWdd6BZMxgwAP79N9dO6+EhjOhQkykDm3PkbAy3j/2LX7YeybXzp7qWU86qlFKFka8vzJp1qXhedHSunv76WqWY/1BbqpYqwjEnVVnVpKCUUrmpShX48ks7GqlvX9izJ1dPX7F4AN/e35p+LSrn6nmTeDnlrEopVZj93//ZiW3PPQdz59qyGIMG2fWewZbLaNYsy8t7JvH2dN73ebdrKYjILSKyQ0R2i8jTro5HKaWy5amnYP9+O+N52TLo0gXatbNb27YQGgqvvQbHj7s40NTEndYLFRFPYCfQEYgAVgO9jTHb0jo+PDzcrFmzJg8jVEqpbIiOhtWrL814PnUKJk2Cn3+2/RDVq2f9nIMHw8iR2QpHRNYaY8LTes7dbh81B3YbY/YAiMg3QBcgzaSglFL5gr+/LbmdUvfusH27XZ/h0KGsn7NMmdyJ7TLulhQqAAdSPI4AWqQ8QESGAkMBKld2TkeLUkrlibp14cMPXR1FKu7WpyBp7Et1f8sYM9EYE26MCS9VqlQehaWUUoWDuyWFCKBSiscVgWy0q5RSSmWHuyWF1UBNEakqIj7AXcA8F8eklFKFhlv1KRhj4kXkQeAXwBOYbIzZ6uKwlFKq0HCrpABgjFkALHB1HEopVRi52+0jpZRSLqRJQSmlVDJNCkoppZK5VZmLrBKR48B+V8fhIiHACVcH4UKF/f2D/g30/Wf//VcxxqQ50StfJ4XCTETWpFe7pDAo7O8f9G+g7985719vHymllEqmSUEppVQyTQr510RXB+Bihf39g/4N9P07gfYpKKWUSqYtBaWUUsk0KSillEqmScHNXW3NahEZKSLbRGSTiCwSkSquiNNZMrtmt4j0EBEjIgVqiGJm3r+I3On4N7BVRL7O6xidKRP//iuLyGIRWe/4f+A2V8TpLCIyWUSOiciWdJ4XERnj+PtsEpEmOb6oMUY3N92wlWL/BaoBPsBGoN5lx7QHAhy/3w/McHXcefn+HccFAn8CK4FwV8edx//9awLrgeKOx6VdHXcev/+JwP2O3+sB+1wddy7/Da4DmgBb0nn+NuAn7AJlLYFVOb2mthTcW/Ka1caYWCBpzepkxpjFxpgox8OV2IWJCoqrvn+HV4G3gZi8DC4PZOb93wt8ZIw5DWCMOZbHMTpTZt6/AYo5fg+igC3KZYz5EziVwSFdgC+MtRIIFpFyObmmJgX3ltaa1RUyOH4w9ltDQXHV9y8ijYFKxpj5eRlYHsnMf/9aQC0RWSYiK0XkljyLzvky8/5fAvqJSAS25P5DeROa28jqZ8RVud16CiqVq65ZnXygSD8gHLjeqRHlrQzfv4h4AKOBAXkVUB7LzH9/L+wtpHbYVuJSEbnGGHPGybHlhcy8/97AFGPMeyLSCvjS8f4TnR+eW8j0Z0RmaUvBvWVqzWoRuRF4FrjdGHMxj2LLC1d7/4HANcASEdmHvac6rwB1Nmfmv38E8L0xJs4YsxfYgU0SBUFm3v9gYCaAMWYF4IctFFdY5Pq69poU3NtV16x23D75BJsQCtL9ZLjK+zfGRBpjQowxocaYUGyfyu3GmDWuCTfXZWbN8rnYwQaISAj2dtKePI3SeTLz/v8DOgCISF1sUjiep1G61jzgHscopJZApDHmcE5OqLeP3JhJZ81qEXkFWGOMmQe8AxQFZokIwH/GmNtdFnQuyuT7L7Ay+f5/AW4SkW1AAvCEMeak66LOPZl8/48Bn4rIo9jbJgOMY1hOQSAi07G3BkMc/SYvAt4AxpiPsf0otwG7gShgYI6vWYD+fkoppXJIbx8ppZRKpklBKaVUMk0KSimlkmlSUEoplUyTglJKqWSaFJTbE5GSIrLBsR0RkYMpHi930jUbi8gkZ5w7q0QkNKlKpoiEi8iYXDjnSyLyuOP3d0XkhpyeUxUMOk9BuT3HuPswsB9mwHljzLtOvuwzwGs5PYnYySOSW2UXHBPzcnty3ljgU+D3XD6vyoe0paDyNRE57/jZTkT+EJGZIrJTREaJSF8R+VtENotIdcdxpUTkWxFZ7djapHHOQKChMWaj4/FLjrr2S0Rkj4iMSHHsSBHZ4tgecewLFZHtIjIeWAdUEpHzIvKWiKwVkd9EpHmK892e4nVLRWSdY2udRmztRGS+4/frU7SY1jviRkSecLy3TSLycorXPit2bYLfgNpJ+40x+4GSIlI2x/9BVL6nLQVVkDQC6mJLDe8BJhljmovIw9jqmY8AHwKjjTF/iUhl7GzZupedJxy4fFGTOthyEoHADhGZADTEziBtgS1MtkpE/gBOYz90BxpjhgOISBFgiTHmKRGZg22FdMSuATAVW67gGNDRGBMjIjWB6Y5Y0vM48IAxZpmIFAViROQmbO2j5o6Y5onIdcAFbJmIxtj/79cBa1Ocax3QBvg2g+upQkCTgipIVifVfRGRf4FfHfs346gPBNwI1HOUBAEoJiKBxphzKc5Tjivr5/zoKDZ4UUSOAWWAtsAcY8wFxzW/A67FfsDvd9S3TxIL/JwinovGmDgR2QyEOvZ7A+NEJAxbsqLWVd7vMuB9EfkK+M4YE+FICjdhF94BWwKlJjaZzUlae0NELi8Rcgwof5XrqUJAk4IqSFJWiE1M8TiRS//WPYBWxpjoDM4TjS2slt65ExznS6tscZILlz2OS1GTJzk2Y0yiiCTF9ihwFNvi8eAqiwYZY0aJyI/Y2jcrxVbLFeBNY8wnKY913NrKqKaNH/Z9q0JO+xRUYfMr8GDSA8e38sttB2pk4lx/Al1FJMBxe6gbsDQHsQUBhx2d0ndji8ClS0SqG2M2G2PewnY+18HeDhvkuJ2EiFQQkdKOWLuJiL+j76HzZaerxZW3zFQhpC0FVdiMAD4SkU3Yf/9/AvelPMAY84+IBKVxW4nLjlsnIlOAvx27Jhlj1otIaDZjGw98KyI9gcVc2dq43CMi0h7bctkG/GSMuSi2hPQKxy2y80A/R6wzgA3AflIkLxHxxibBglJyXOWAVklVKg2OUsznjDFuMVfBmUSkG9DEGPO8q2NRrqe3j5RK2wRS9yMUZF7Ae64OQrkHbSkopZRKpi0FpZRSyTQpKKWUSqZJQSmlVDJNCkoppZJpUlBKKZXs/wGyj+G37EPhowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0,1,100)\n",
    "plt.plot(x[5:100],data[5:100],label=\"Data\")\n",
    "plt.plot(x[5:100],(num_nmps_time)[5:100],label=\"Simulated\", color = 'r')\n",
    "plt.xlabel('Time (normalised)')\n",
    "plt.ylabel('NMp counts')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = sampler.get_chain(flat=True)\n",
    "print('samples')\n",
    "print(samples.shape)\n",
    "print(type(samples))\n",
    "print(' ')\n",
    "with open(\"samples_emcee.txt\", \"wb\") as fp:   #Pickling\n",
    "    pickle.dump(samples, fp)\n",
    "\n",
    "log_probs = sampler.flatlnprobability\n",
    "print('log_probs')\n",
    "print(log_probs.shape)\n",
    "print(type(log_probs))\n",
    "print(' ')\n",
    "with open(\"log_probs.txt\", \"wb\") as fp:   #Pickling\n",
    "    pickle.dump(log_probs, fp)\n",
    "    \n",
    "map_params = samples[np.argmax(sampler.flatlnprobability)]\n",
    "print('map_params')\n",
    "print(map_params.shape)\n",
    "print(type(map_params))\n",
    "print(' ')\n",
    "with open(\"map_params_emcee.txt\", \"wb\") as fp:   #Pickling\n",
    "    pickle.dump(map_params, fp)\n",
    "    \n",
    "acceptance_fraction = sampler.acceptance_fraction\n",
    "print('acceptance_fraction')\n",
    "print(acceptance_fraction.shape)\n",
    "print(type(acceptance_fraction))\n",
    "print(' ')\n",
    "with open(\"acceptance_fraction_emcee.txt\", \"wb\") as fp:   #Pickling\n",
    "    pickle.dump(acceptance_fraction, fp)\n",
    "    \n",
    "autocorr_time = sampler.get_autocorr_time(quiet=True)\n",
    "print('autocorr_time')\n",
    "print(autocorr_time.shape)\n",
    "print(type(autocorr_time))\n",
    "print(' ')\n",
    "with open(\"autocorr_time_emcee.txt\", \"wb\") as fp:   #Pickling\n",
    "    pickle.dump(autocorr_time, fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.05844943e-02, 9.65632033e-01, 4.52272889e-02, 3.30898025e-01,\n",
       "       4.10382923e-01, 4.17411003e-01, 4.92825227e+00, 4.68364994e+00,\n",
       "       1.70012427e+02, 2.91903166e-01])"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myparams"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "341.325px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}