ParticleFiltering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm as norm\n",
    "from scipy.stats import poisson as poisson\n",
    "from scipy.stats import nbinom\n",
    "import numpy.random as random\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.patches as patches\n",
    "import matplotlib.ticker as plticker\n",
    "from matplotlib import rc #for latex use on the plots\n",
    "from matplotlib import ticker  # for labels formattings\n",
    "from scipy.stats import poisson\n",
    "import datetime\n",
    "from os import getcwd,path\n",
    "import pandas as pd\n",
    "from numpy import pi as pi\n",
    "from time import time as timer\n",
    "%matplotlib inline\n",
    "from joblib import Parallel, delayed\n",
    "import multiprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Plot Configs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "hidden": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       ".output_png {\n",
       "    display: table-cell;\n",
       "    text-align: center;\n",
       "    vertical-align: middle;\n",
       "}\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.style.use('ggplot')\n",
    "perso_linewidth = 0.6\n",
    "# This will change your default rcParams\n",
    "def init_plotting():\n",
    "    plt.rcParams['figure.figsize'] = (40,20)\n",
    "    plt.rcParams['font.size'] = 30\n",
    "    plt.rcParams['font.family'] = 'serif'\n",
    "    plt.rcParams['axes.labelsize'] = plt.rcParams['font.size']\n",
    "    plt.rcParams['axes.titlesize'] = 1.5*plt.rcParams['font.size']\n",
    "    plt.rcParams['legend.fontsize'] = plt.rcParams['font.size']\n",
    "    plt.rcParams['xtick.labelsize'] = plt.rcParams['font.size']\n",
    "    plt.rcParams['ytick.labelsize'] = plt.rcParams['font.size']\n",
    "    #plt.rcParams['savefig.dpi'] = 2*plt.rcParams['savefig.dpi']\n",
    "    plt.rcParams['axes.linewidth'] = perso_linewidth\n",
    "    plt.rcParams['savefig.dpi'] = '300'\n",
    "    plt.rcParams['savefig.format'] = 'pdf'\n",
    "    plt.rcParams['axes.facecolor'] = '#D3D3D3'\n",
    "    plt.rcParams['axes.edgecolor'] = '0'\n",
    "    plt.rcParams['axes.grid'] = True\n",
    "    plt.rcParams['grid.color']='white'\n",
    "    plt.rcParams['grid.linestyle'] = '-'\n",
    "    plt.rcParams['grid.linewidth'] = '0.4'\n",
    "    plt.rcParams['axes.axisbelow'] = True\n",
    "    plt.rcParams['legend.edgecolor'] = 'black'\n",
    "    plt.rcParams['legend.facecolor'] = 'white'\n",
    "    plt.rcParams['lines.markersize']= 18 \n",
    "    plt.rcParams['lines.markeredgewidth']= '0.1'\n",
    "    plt.rcParams['lines.color']= 'r' \n",
    "    plt.rcParams['lines.marker']= 'o' \n",
    "    plt.rcParams['lines.linestyle']= '-' \n",
    "    plt.rcParams['xtick.color']= '0'\n",
    "    plt.rcParams['ytick.color']= '0'\n",
    "    #plt.rcParams['axes.color_cycle']= ['#3778bf', '#feb308', '#a8a495', '#7bb274', '#825f87']\n",
    "    plt.gca().spines['right'].set_color('none')\n",
    "    plt.gca().spines['right'].set_visible('False')\n",
    "    plt.gca().spines['top'].set_visible('False')\n",
    "    plt.gca().spines['top'].set_color('none')\n",
    "    plt.gca().xaxis.set_ticks_position('bottom')\n",
    "    plt.gca().yaxis.set_ticks_position('left')\n",
    "    plt.rcParams['ytick.minor.size']= 1.5\n",
    "    plt.rcParams['ytick.major.width']= perso_linewidth\n",
    "    plt.rcParams['ytick.minor.width']= perso_linewidth\n",
    "    plt.rcParams['xtick.major.width']= perso_linewidth\n",
    "    plt.rcParams['xtick.minor.width']= perso_linewidth\n",
    "\n",
    "init_plotting()\n",
    "\n",
    "plt.close() # this line to avoid an empty plot showing up\n",
    "\n",
    "### Centering the plots\n",
    "from IPython.core.display import HTML\n",
    "HTML(\"\"\"\n",
    "<style>\n",
    ".output_png {\n",
    "    display: table-cell;\n",
    "    text-align: center;\n",
    "    vertical-align: middle;\n",
    "}\n",
    "</style>\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stochastic Simulation Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initial state and constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simple SIS Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false
   },
   "source": [
    "##### Simple SIS model\n",
    "\n",
    "\\begin{itemize}\n",
    "\\item $E_1:X \\rightarrow X+1,Y \\rightarrow Y-1$\n",
    "\\item $E_2:X \\rightarrow X-1,Y \\rightarrow Y+1$\n",
    "\\end{itemize}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 403,
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "editable": false,
    "hide_input": true,
    "run_control": {
     "frozen": true
    }
   },
   "outputs": [],
   "source": [
    "initial_state=[30,70]\n",
    "\n",
    "### Define events increments for Gillepsie's Direct and FR: +1/-1 increment added to each compartment when Ei fires\n",
    "E1=[-1,1]\n",
    "E2=[1,-1]\n",
    "E=np.vstack((E1,E2))\n",
    "\n",
    "\n",
    "### Define state-change vector for tau-Leap\n",
    "v1=[-1,1]\n",
    "v2=[1,-1]\n",
    "v=np.vstack((v1,v2))\n",
    "\n",
    "\n",
    "simulation_period=10  ### Years\n",
    "\n",
    "dt=1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "source": [
    "#### Boarding School SIR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### Boarding School\n",
    "\n",
    "##### Actual case count\n",
    "csv_file_path=path.join(getcwd(), 'data','Boarding.csv')\n",
    "df=pd.read_csv(csv_file_path,index_col=0,header=0)\n",
    "\n",
    "t_reported=df.day\n",
    "Y_reported=df.flu\n",
    "\n",
    "\n",
    "\n",
    "# ### Define events increments for Gillepsie's Direct and FR: +1/-1 increment added to each compartment when Ei fires\n",
    "#E1=[1,0,0]  ## New birth\n",
    "E1=[-1,1,0] ## Infection\n",
    "#E3=[-1,0,0] ## Death in S\n",
    "E2=[0,-1,1] ## Recovery\n",
    "#E5=[0,-1,0] ## Death in I\n",
    "#E6=[0,0,-1] ## Death in R\n",
    "\n",
    "\n",
    "E=np.vstack((E1,E2))\n",
    "\n",
    "\n",
    "### Define state-change vector for tau-Leap: Signs before each rate\n",
    "v1=[-1,0,0]\n",
    "v2=[1,-1,0]\n",
    "v3=[0,1,-1]\n",
    "v4=[0,0,1]\n",
    "v=np.vstack((v1,v2,v3,v4))\n",
    "\n",
    "simulation_period=1/16  ### Years\n",
    "\n",
    "dt=0.01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Ratefunc(state,Beta,Gamma,Delta):\n",
    "    \n",
    "    \n",
    "    S=state[0]\n",
    "    I=state[1]\n",
    "    C=state[2]\n",
    "    R=state[3]\n",
    "    \n",
    "    N=S+I+C+R\n",
    "    #######\n",
    "    \n",
    "    Rates=np.zeros(0)\n",
    "    Rates=np.append(Rates,[(Beta*I*S)/N])\n",
    "    Rates=np.append(Rates,[Gamma*I])\n",
    "    Rates=np.append(Rates,[Delta*C])\n",
    "    \n",
    "        \n",
    "    return Rates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $\\tau$-Leap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TauLeap(state,Beta,Gamma,v,dt): \n",
    "    ### Returns the new state (Xn)\n",
    "    \n",
    "    Rates=Ratefunc(state,Beta,Gamma,dt)\n",
    "    \n",
    "    Lambda=dt*Rates\n",
    "    Lambda_nbin=dt*Rates\n",
    "    \n",
    "    if state[0]>0:\n",
    "        Lambda_nbin[0]=Lambda_nbin[0]/state[0]\n",
    "    if state[1]>0:\n",
    "        Lambda_nbin[1]=Lambda_nbin[1]/state[1]\n",
    "    if state[2]>0:\n",
    "        Lambda_nbin[2]=Lambda_nbin[2]/state[2]\n",
    "    \n",
    "      \n",
    "    probs=1-np.exp(-Lambda_nbin)\n",
    "    dMs_nbin=np.zeros(3)\n",
    "    \n",
    "    \n",
    "    \n",
    "    temp=np.ones(len(state))*-1\n",
    "    N=sum(state)\n",
    "    counter=1\n",
    "    while any(value<0 for value in temp):\n",
    "            \n",
    "            temp=np.copy(state)\n",
    "            \n",
    "#             dMs=np.random.poisson(Lambda)\n",
    "#             dstate=v@dMs.T\n",
    "                      \n",
    "            \n",
    "\n",
    "            for i in range(3):\n",
    "                 dMs_nbin[i]=np.random.binomial(state[i], probs[i], size=1)\n",
    "             \n",
    "            dstate=v@dMs_nbin.T   \n",
    "        \n",
    " \n",
    "            \n",
    "            \n",
    "            temp[0]=temp[0]+dstate[0]\n",
    "            temp[1]=temp[1]+dstate[1]\n",
    "            temp[2]=temp[2]+dstate[2]\n",
    "            temp[3]=N-temp[0]-temp[1]-temp[2]\n",
    "            \n",
    "        \n",
    "            counter=counter+1\n",
    "            \n",
    "        \n",
    "    \n",
    "    state=np.asarray(temp)\n",
    "     \n",
    "    return state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "code_folding": [
     0
    ],
    "collapsed": true,
    "deletable": false,
    "editable": false,
    "run_control": {
     "frozen": true
    }
   },
   "outputs": [],
   "source": [
    "def TauLeapSim(initial_state,v,dt,simulation_period):   ### dt:days, simulation_period:years\n",
    "    state=np.array(initial_state)\n",
    "    n_vars=len(initial_state)\n",
    "    state.shape=(1,n_vars)\n",
    "    \n",
    "    simulation_period_days=simulation_period*365\n",
    "    \n",
    "    n=int(simulation_period_days/dt)\n",
    "    \n",
    "    \n",
    "    t=np.array([0])\n",
    "    \n",
    "    \n",
    "    for i in range(n):\n",
    "        \n",
    "        R=Ratefunc(state[-1])\n",
    "        new_state=TauLeap(state[-1],R,v,dt)\n",
    "        state=np.vstack((state,new_state))\n",
    "        t=np.append(t,t[-1]+dt)\n",
    "#         if any(state[-1]<=0):\n",
    "#             break\n",
    "        \n",
    "        \n",
    "    return state,t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Weight Function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def measure(Counts,Y,arg1,arg2):\n",
    "    \n",
    "    tau=arg1\n",
    "    nbr=1/tau/tau\n",
    "    Counts[Counts==0]=1\n",
    "    w=nbinom.pmf(k=Y, n=nbr, p=1/(1+Counts/nbr))\n",
    "    #var=X\n",
    "    #var[var==0]=1\n",
    "    #mean=X\n",
    "    #w=norm(loc=mean,scale=var).pdf(Y)\n",
    "    #w=poisson.pmf(mu=Counts,k=Y)\n",
    "    #w=nbinom.pmf(k=Y, n=C, p=Rho)\n",
    "    \n",
    "    return w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bootstraping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def BootStrap(S,I,C,R,Y_obs,arg1,arg2):\n",
    "    \n",
    "    \n",
    "    pobs=measure(C,Y_obs,arg1,arg2)\n",
    "    \n",
    "    w=pobs/sum(pobs)\n",
    "    \n",
    "    \n",
    "    prediction=np.vstack((S,I,C,R))\n",
    "    idx=np.arange(prediction.shape[1])\n",
    "    filtered_idx=np.random.choice(idx,n_particles,p=w)\n",
    "    filtered=prediction[:,filtered_idx]\n",
    "    S=filtered[0,:]\n",
    "    I=filtered[1,:]\n",
    "    C=filtered[2,:]\n",
    "    R=filtered[3,:]\n",
    "    \n",
    "\n",
    "    return S,I,C,R,w,pobs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Direct Method Simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=764\n",
    "S0=763\n",
    "I0=1\n",
    "C0=0\n",
    "R0=0\n",
    "\n",
    "\n",
    "n_particles=50000\n",
    "particles=np.zeros((n_particles,3))\n",
    "Beta=3.2\n",
    "Gamma=1.5\n",
    "Delta=0.5\n",
    "tau=0.1\n",
    "\n",
    "\n",
    "\n",
    "obs_time=t_reported.values\n",
    "time=t_reported.values\n",
    "time=np.insert(time,0,0)\n",
    "dt=time[1]-time[0]\n",
    "dt=1/5\n",
    "t0=1/5\n",
    "time=np.arange(0,14+dt,dt)\n",
    "\n",
    "\n",
    "Y_obs=Y_reported.values\n",
    "#Y_obs=np.insert(Y_obs,0,0)\n",
    "\n",
    "\n",
    "\n",
    "counter=0\n",
    "\n",
    "\n",
    "S=np.zeros((n_particles,len(time)+1))\n",
    "I=np.zeros((n_particles,len(time)+1))\n",
    "C=np.zeros((n_particles,len(time)+1))\n",
    "R=np.zeros((n_particles,len(time)+1))\n",
    "\n",
    "Likelihood_t=np.zeros(14)\n",
    "logll=np.zeros(14)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "S[:,0]=S0\n",
    "I[:,0]=I0\n",
    "C[:,0]=C0\n",
    "R[:,0]=R0\n",
    "\n",
    "\n",
    "w=np.zeros(n_particles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Serial processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tau-leaping for 50000.000000 particles time: 62.061847\n",
      "Tau-leaping for 50000.000000 particles time: 66.084810\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-87-d682ede86c6b>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     15\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     16\u001b[0m         \u001b[0mS\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mstep\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnew_state\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 17\u001b[1;33m         \u001b[0mI\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mstep\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnew_state\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     18\u001b[0m         \u001b[0mC\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mstep\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnew_state\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     19\u001b[0m         \u001b[0mR\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mstep\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnew_state\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "#fig=plt.figure(figsize=(20,10))\n",
    "#ax=plt.subplot(111)\n",
    "#observed=plt.plot(obs_time,Y_obs,'-ok',linewidth=3,label='Actual reported cases',zorder=10)\n",
    "\n",
    "for step,t in enumerate(time):#time:\n",
    "    \n",
    "    if any(value<0 for value in S[:,step]):\n",
    "        break\n",
    "    time_before = timer()\n",
    "    for Particle in range(n_particles):\n",
    "        \n",
    "        current_state=[S[Particle,step],I[Particle,step],C[Particle,step],R[Particle,step]]\n",
    "        \n",
    "        new_state=TauLeap(current_state,Beta,Gamma,v,dt)\n",
    "        \n",
    "        S[:,step+1]=new_state[0]\n",
    "        I[:,step+1]=new_state[1]\n",
    "        C[:,step+1]=new_state[2]\n",
    "        R[:,step+1]=new_state[3]\n",
    "        \n",
    "                \n",
    "\n",
    "    time_after = timer()\n",
    "    func_time=time_after-time_before\n",
    "    print(\"Tau-leaping for %f particles time: %f\"%(n_particles,func_time))\n",
    "    C_prediction=C[:,step+1]\n",
    "    \n",
    "    \n",
    "    \n",
    "    if t in obs_time:\n",
    "        time_before = timer()\n",
    "        S[:,step+1],I[:,step+1],C[:,step+1],R[:,step+1],w,pobs=BootStrap(S[:,step+1],I[:,step+1],C[:,step+1],R[:,step+1],Y_obs[int(t)-1],arg1=tau,arg2=0)\n",
    "        time_after = timer()\n",
    "        func_time=time_after-time_before\n",
    "        print(\"Bootstraping time for %f particles: %f\"%(n_particles,func_time))\n",
    "        logll[int(t)-1]=np.log(np.mean(pobs))\n",
    "        #P=measure(C[:,step+1],Y_obs[int(t)-1],tau,0)\n",
    "       \n",
    "        #Likelihood_t[int(t)-1]=np.mean(P)\n",
    "        scale=1000/max(w*w*w)\n",
    "        print(\"time,obs,mean process:\",t,Y_obs[int(t)-1],np.mean(C[:,step+1]))\n",
    "        #PredcitionParticles=ax.scatter(np.ones(n_particles)*(t),C_prediction,marker='o',c='grey',label='Prediction Particles (Particle sizes scaled to their relative weights)',s=scale*(w*w*w),zorder=9)   ### Prediction particles\n",
    "\n",
    "        \n",
    "        time=np.arange(0,14+2*dt,dt)\n",
    "        \n",
    "LL=sum(logll)  \n",
    "print(sum(logll))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parallel processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tau-leaping for 50000.000000 particles time: 3.676531\n",
      "Tau-leaping for 50000.000000 particles time: 2.664212\n",
      "Tau-leaping for 50000.000000 particles time: 2.702446\n",
      "Tau-leaping for 50000.000000 particles time: 2.636783\n",
      "Tau-leaping for 50000.000000 particles time: 2.710520\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-88-cbd2ca9902e3>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     16\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     17\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 18\u001b[1;33m     \u001b[0mnew_state\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mParallel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn_jobs\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnum_cores\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdelayed\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mTauLeap\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcurrent_state\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mParticle\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mBeta\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mGamma\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdt\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mParticle\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mn_particles\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     19\u001b[0m     \u001b[0mnew_state\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnew_state\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     20\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\joblib\\parallel.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, iterable)\u001b[0m\n\u001b[0;32m    960\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    961\u001b[0m             \u001b[1;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_backend\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mretrieval_context\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 962\u001b[1;33m                 \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mretrieve\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    963\u001b[0m             \u001b[1;31m# Make sure that we get a last message telling us we are done\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    964\u001b[0m             \u001b[0melapsed_time\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_start_time\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\joblib\\parallel.py\u001b[0m in \u001b[0;36mretrieve\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    863\u001b[0m             \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    864\u001b[0m                 \u001b[1;32mif\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_backend\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'supports_timeout'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 865\u001b[1;33m                     \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_output\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mjob\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    866\u001b[0m                 \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    867\u001b[0m                     \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_output\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mjob\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\joblib\\_parallel_backends.py\u001b[0m in \u001b[0;36mwrap_future_result\u001b[1;34m(future, timeout)\u001b[0m\n\u001b[0;32m    513\u001b[0m         AsyncResults.get from multiprocessing.\"\"\"\n\u001b[0;32m    514\u001b[0m         \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 515\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mfuture\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mresult\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    516\u001b[0m         \u001b[1;32mexcept\u001b[0m \u001b[0mLokyTimeoutError\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    517\u001b[0m             \u001b[1;32mraise\u001b[0m \u001b[0mTimeoutError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\joblib\\externals\\loky\\_base.py\u001b[0m in \u001b[0;36mresult\u001b[1;34m(self, timeout)\u001b[0m\n\u001b[0;32m    424\u001b[0m                 \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__get_result\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    425\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 426\u001b[1;33m             \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_condition\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    427\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    428\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_state\u001b[0m \u001b[1;32min\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mCANCELLED\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mCANCELLED_AND_NOTIFIED\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Continuum\\anaconda3\\lib\\threading.py\u001b[0m in \u001b[0;36mwait\u001b[1;34m(self, timeout)\u001b[0m\n\u001b[0;32m    293\u001b[0m         \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m    \u001b[1;31m# restore state no matter what (e.g., KeyboardInterrupt)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    294\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mtimeout\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 295\u001b[1;33m                 \u001b[0mwaiter\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0macquire\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    296\u001b[0m                 \u001b[0mgotit\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    297\u001b[0m             \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "#fig=plt.figure(figsize=(20,10))\n",
    "#ax=plt.subplot(111)\n",
    "#observed=plt.plot(obs_time,Y_obs,'-ok',linewidth=3,label='Actual reported cases',zorder=10)\n",
    "\n",
    "\n",
    "num_cores = multiprocessing.cpu_count()\n",
    "\n",
    "for step,t in enumerate(time):#time:\n",
    "    \n",
    "    if any(value<0 for value in S[:,step]):\n",
    "        break\n",
    "    time_before = timer()\n",
    "    \n",
    "    current_state=np.vstack((S[:,step],I[:,step],C[:,step],R[:,step])).T\n",
    "    \n",
    "    \n",
    "\n",
    "    new_state=Parallel(n_jobs=num_cores)(delayed(TauLeap)(current_state[Particle,:],Beta,Gamma,v,dt) for Particle in range(n_particles))\n",
    "    new_state=np.asarray(new_state)\n",
    "    \n",
    "    S[:,step+1]=new_state[:,0]\n",
    "    I[:,step+1]=new_state[:,1]\n",
    "    C[:,step+1]=new_state[:,2]\n",
    "    R[:,step+1]=new_state[:,3]\n",
    "        \n",
    "        \n",
    "\n",
    "    time_after = timer()\n",
    "    func_time=time_after-time_before\n",
    "    print(\"Tau-leaping for %f particles time: %f\"%(n_particles,func_time))\n",
    "    C_prediction=C[:,step+1]\n",
    "    \n",
    "    \n",
    "    \n",
    "    if t in obs_time:\n",
    "        time_before = timer()\n",
    "        S[:,step+1],I[:,step+1],C[:,step+1],R[:,step+1],w,pobs=BootStrap(S[:,step+1],I[:,step+1],C[:,step+1],R[:,step+1],Y_obs[int(t)-1],arg1=tau,arg2=0)\n",
    "        time_after = timer()\n",
    "        func_time=time_after-time_before\n",
    "        print(\"Bootstraping time for %f particles: %f\"%(n_particles,func_time))\n",
    "        logll[int(t)-1]=np.log(np.mean(pobs))\n",
    "        #P=measure(C[:,step+1],Y_obs[int(t)-1],tau,0)\n",
    "       \n",
    "        #Likelihood_t[int(t)-1]=np.mean(P)\n",
    "        scale=1000/max(w*w*w)\n",
    "        print(\"time,obs,mean process:\",t,Y_obs[int(t)-1],np.mean(C[:,step+1]))\n",
    "        #PredcitionParticles=ax.scatter(np.ones(n_particles)*(t),C_prediction,marker='o',c='grey',label='Prediction Particles (Particle sizes scaled to their relative weights)',s=scale*(w*w*w),zorder=9)   ### Prediction particles\n",
    "\n",
    "        \n",
    "        time=np.arange(0,14+2*dt,dt)\n",
    "        \n",
    "LL=sum(logll)  \n",
    "print(sum(logll))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "# for Particle in range(n_particles):\n",
    "#     filtered_particles=ax.plot(time,C[Particle,:],'--or',alpha=0.5,zorder=1,label='Filtered Particles',markersize=3)\n",
    "\n",
    "# plt.xlim((0,15))\n",
    "# plt.ylim((-10,500))\n",
    "# plt.legend(handles=[filtered_particles[0],PredcitionParticles,observed[0]],loc=0,fontsize=20)\n",
    "# plt.xlabel('Time (days)')\n",
    "# plt.ylabel('Convalescents')\n",
    "# plt.title(r'Particle filtering for boarding school ($\\beta$=%.2f, $\\gamma$=%.2f, $\\delta$=%.2f) - LogLikelihood=%.2f' %(Beta,Gamma,Delta,LL),fontsize=20)\n",
    "# plt.show()\n",
    "#LogLikelihood=sum(np.log(Likelihood_t))\n",
    "#print(LogLikelihood)\n",
    "# print(sum(logll))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "logll"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python [default]",
   "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.6.3"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "oldHeight": 760,
   "position": {
    "height": "782px",
    "left": "1245px",
    "right": "20px",
    "top": "78px",
    "width": "600px"
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "varInspector_section_display": "block",
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}