From eb039d3022936411cafe094bf42167da2fb86177 Mon Sep 17 00:00:00 2001 From: Phil McCown Date: Fri, 3 Nov 2017 10:51:41 -0400 Subject: [PATCH 01/15] Initial start to the Python scripts --- Ex10Q1.py | 3 +++ Ex10Q2.py | 6 ++++++ 2 files changed, 9 insertions(+) create mode 100755 Ex10Q1.py create mode 100755 Ex10Q2.py diff --git a/Ex10Q1.py b/Ex10Q1.py new file mode 100755 index 0000000..5cee4e1 --- /dev/null +++ b/Ex10Q1.py @@ -0,0 +1,3 @@ +import numpy +import pandas +import scipy diff --git a/Ex10Q2.py b/Ex10Q2.py new file mode 100755 index 0000000..b5a1851 --- /dev/null +++ b/Ex10Q2.py @@ -0,0 +1,6 @@ +import numpy +import pandas +import scipy +#SIR model in epidemiology +#N=S+I+R (susceptible-->B-->infected-->g-->resistant); unidirectional +#R0=((B*(S+I+R))/g) ; R0 is how bad the disease is in the population; see YY's notes for where max daily incidence, max daily prevalence, percent affected, and basic reproduction number are found in a graph From 058ab91cb321accfe213dbb49b6c6e4e766df61f Mon Sep 17 00:00:00 2001 From: Phil McCown Date: Fri, 3 Nov 2017 11:24:20 -0400 Subject: [PATCH 02/15] Q1 almost done sans ggplot --- Ex10Q1.py | 33 +++++++++++++++++++++++++++++++++ Ex10Q2.py | 1 + 2 files changed, 34 insertions(+) diff --git a/Ex10Q1.py b/Ex10Q1.py index 5cee4e1..c955425 100755 --- a/Ex10Q1.py +++ b/Ex10Q1.py @@ -1,3 +1,36 @@ import numpy import pandas import scipy +import scipy.integrate as spint +from plotnine import * +def ddSim(y,t0,r,K): + N=y[0] + dNdt=r*(1-N/K)*N + return [dNdt] +times=(range(0,1000)) +rs=[-0.1,0.1,0.4,0.8,1] +store_rs=pandas.DataFrame({"time":times,"rl":0,"r2":0,"r3":0,"r4":0,"r5":0}) +for i in range(0,len(rs)): + K=100 + N=10 + pars1=(rs[i],K) + sim=spint.odeint(func=ddSim,y0=N,t=times,args=pars1) + store_rs.iloc[:,i]=sim[:,0] + +ks=[10,50,100] +store_ks=pandas.DataFrame({"time":times,"kl":0,"k2":0,"k3":0}) +for i in range(0,len(ks)): + N=1 + pars2=(0.2,ks[i]) + sim2=spint.odeint(func=ddSim,y0=N,t=times,args=pars2) + store_ks.iloc[:,i]=sim2[:,0] + +ns=[1,50,100] +store_ns=pandas.DataFrame({"time":times,"nl":0,"n2":0,"n3":0}) +#problem child is below +for i in range(0,len()): + pars3=(0.1,50) + sim3=spint.odeint(func=ddSim,y0=N[i],t=times,args=pars3) + store_ns.iloc[:,i]=sim3[:,0] + + diff --git a/Ex10Q2.py b/Ex10Q2.py index b5a1851..07fb63b 100755 --- a/Ex10Q2.py +++ b/Ex10Q2.py @@ -1,6 +1,7 @@ import numpy import pandas import scipy +from plotnine import * #SIR model in epidemiology #N=S+I+R (susceptible-->B-->infected-->g-->resistant); unidirectional #R0=((B*(S+I+R))/g) ; R0 is how bad the disease is in the population; see YY's notes for where max daily incidence, max daily prevalence, percent affected, and basic reproduction number are found in a graph From c837d3cfd64ef1c39c3a1d0053df142d857dcef9 Mon Sep 17 00:00:00 2001 From: kkilgoreND <31992787+kkilgoreND@users.noreply.github.com> Date: Mon, 6 Nov 2017 20:30:58 -0500 Subject: [PATCH 03/15] updated Q2 with codes, dataframe and model --- Ex10Q2.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/Ex10Q2.py b/Ex10Q2.py index 07fb63b..cac2be6 100755 --- a/Ex10Q2.py +++ b/Ex10Q2.py @@ -5,3 +5,33 @@ #SIR model in epidemiology #N=S+I+R (susceptible-->B-->infected-->g-->resistant); unidirectional #R0=((B*(S+I+R))/g) ; R0 is how bad the disease is in the population; see YY's notes for where max daily incidence, max daily prevalence, percent affected, and basic reproduction number are found in a graph + +#state vairables=I/S=state parameters=r/B +#basic reproduction only initial values + +#equations are dS/dt=-BIS, dI/dt=BIS-RI and dR/dt=RI +def epiSim(y,t,B,r): + I=y[0] + S=y[1] + + dSdt= (-B*I*S) + dIdt= (B*I*S)-(r*I) + dRdt= (r*I) + + return[dSdt,dIdt,dRdt] + +#generate a dataframe containing a list of the B and r parameters +parB=[0.0005,0.005,0.0001,0.00005,0.0001,0.0002,0.0001] +parr=[0.05,0.5,0.1,0.1,0.05,0.05,0.06] +par_array=numpy.column_stack([parB,parr]) +par_df=pandas.DataFrame(par_array,columns=['B','r']) + +#extract params from dataframe +for i in range(0,6): + I0=[1] + S0=[999] + params=(par_df.B[i],par_df.r[i]) + times=range(0,500) + #creat a model using odeint + modelSim=spint.odeint(func=epiSim,y0=(I0,S0),t=times,args=params) + From 16bd80d6b274284edfa36f0783fd271549fd88ed Mon Sep 17 00:00:00 2001 From: Phil McCown Date: Thu, 9 Nov 2017 15:25:23 -0500 Subject: [PATCH 04/15] Ex10Q1 is done; Ex10Q2 is half done, I'm going over the max prevalence and max infectivity rates --- Ex10Q1.py | 15 ++++++++++++--- Ex10Q2.py | 12 ++++++++---- 2 files changed, 20 insertions(+), 7 deletions(-) diff --git a/Ex10Q1.py b/Ex10Q1.py index c955425..bd77b16 100755 --- a/Ex10Q1.py +++ b/Ex10Q1.py @@ -8,6 +8,7 @@ def ddSim(y,t0,r,K): dNdt=r*(1-N/K)*N return [dNdt] times=(range(0,1000)) + rs=[-0.1,0.1,0.4,0.8,1] store_rs=pandas.DataFrame({"time":times,"rl":0,"r2":0,"r3":0,"r4":0,"r5":0}) for i in range(0,len(rs)): @@ -27,10 +28,18 @@ def ddSim(y,t0,r,K): ns=[1,50,100] store_ns=pandas.DataFrame({"time":times,"nl":0,"n2":0,"n3":0}) -#problem child is below -for i in range(0,len()): +for i in range(0,len(ns)): + N=1,50,100 pars3=(0.1,50) sim3=spint.odeint(func=ddSim,y0=N[i],t=times,args=pars3) store_ns.iloc[:,i]=sim3[:,0] - +#Placing models into dataframe +modelOutputRS=pandas.DataFrame({"t":times,"N":sim[:,0]}) +modelOutputKS=pandas.DataFrame({"t":times,"N":sim2[:,0]}) +modelOutputNS=pandas.DataFrame({"t":times,"N":sim3[:,0]}) + +#Plot simulation outputs into ggplot +ggplot(modelOutputKS,aes(x="t",y="N"))+geom_line()+theme_classic() +ggplot(modelOutputRS,aes(x="t",y="N"))+geom_line()+theme_classic() +ggplot(modelOutputNS,aes(x="t",y="N"))+geom_line()+theme_classic() \ No newline at end of file diff --git a/Ex10Q2.py b/Ex10Q2.py index cac2be6..bdf0dee 100755 --- a/Ex10Q2.py +++ b/Ex10Q2.py @@ -1,6 +1,7 @@ import numpy import pandas import scipy +import scipy.integrate as spint from plotnine import * #SIR model in epidemiology #N=S+I+R (susceptible-->B-->infected-->g-->resistant); unidirectional @@ -29,9 +30,12 @@ def epiSim(y,t,B,r): #extract params from dataframe for i in range(0,6): I0=[1] - S0=[999] + S0=[999,1,0] params=(par_df.B[i],par_df.r[i]) times=range(0,500) - #creat a model using odeint - modelSim=spint.odeint(func=epiSim,y0=(I0,S0),t=times,args=params) - + #create a model using odeint + modelSim=spint.odeint(func=epiSim,y0=S0,t=times,args=params) +#maxDailyIncidence is the maxInfectivityRate +#maxDailyPrevalence is the half-max of the equation +maxDailyIncidence= +maxDailyPrevalence= \ No newline at end of file From 26ebd1577bb24f743e3e2b127254273ac5dfdd4e Mon Sep 17 00:00:00 2001 From: Michelle Corley Date: Thu, 9 Nov 2017 16:01:21 -0500 Subject: [PATCH 05/15] some of first plot done --- Ex10Q1.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Ex10Q1.py b/Ex10Q1.py index c955425..afaba97 100755 --- a/Ex10Q1.py +++ b/Ex10Q1.py @@ -16,6 +16,8 @@ def ddSim(y,t0,r,K): pars1=(rs[i],K) sim=spint.odeint(func=ddSim,y0=N,t=times,args=pars1) store_rs.iloc[:,i]=sim[:,0] +a=ggplot(store_rs,aes(x="time",y="r2")) +a+geom_point()+coord_cartesian() ks=[10,50,100] store_ks=pandas.DataFrame({"time":times,"kl":0,"k2":0,"k3":0}) From e35bfd07401a14b6780a4371c342a89dd6a2a7bd Mon Sep 17 00:00:00 2001 From: Phil McCown Date: Thu, 9 Nov 2017 16:06:29 -0500 Subject: [PATCH 06/15] All done with one odd error somewhere in Q2 --- Ex10Q2.py | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/Ex10Q2.py b/Ex10Q2.py index bdf0dee..ef06c79 100755 --- a/Ex10Q2.py +++ b/Ex10Q2.py @@ -14,28 +14,32 @@ def epiSim(y,t,B,r): I=y[0] S=y[1] - dSdt= (-B*I*S) dIdt= (B*I*S)-(r*I) dRdt= (r*I) - return[dSdt,dIdt,dRdt] #generate a dataframe containing a list of the B and r parameters parB=[0.0005,0.005,0.0001,0.00005,0.0001,0.0002,0.0001] parr=[0.05,0.5,0.1,0.1,0.05,0.05,0.06] -par_array=numpy.column_stack([parB,parr]) -par_df=pandas.DataFrame(par_array,columns=['B','r']) - +DataOut=pandas.DataFrame(numpy.zeros([7,6]),columns=['MaxIncidence','MaxPrevalence','PercentAffected','R0','Beta','Gamma']) +times=range(0,500) #extract params from dataframe for i in range(0,6): I0=[1] S0=[999,1,0] - params=(par_df.B[i],par_df.r[i]) + parms=(parB[i],parr[i]) times=range(0,500) - #create a model using odeint - modelSim=spint.odeint(func=epiSim,y0=S0,t=times,args=params) -#maxDailyIncidence is the maxInfectivityRate -#maxDailyPrevalence is the half-max of the equation -maxDailyIncidence= -maxDailyPrevalence= \ No newline at end of file + sim=spint.odeint(func=epiSim,y0=S0,t=times,args=parms) + inc=sim[1:500,1]-sim[0:499,1] + maxinc=numpy.max(inc)#create a model using odeint + prevalence=((sim[:,1])/(sim[:,0]+sim[:,1]+sim[:,2])) + maxprev=numpy.max(prevalence) + affected=((sim[499,1]+sim[499,2])/(sim[499,0]+sim[499,1]+sim[499,2])) + R0=(parB[i])*1000/parr[i]/parr[i] + DataOut.iloc[i,0]=maxinc + DataOut.iloc[i,1]=maxprev + DataOut.iloc[i,2]=affected + DataOut.iloc[i,3]=R0 + DataOut.iloc[i,4]=parB[i] + DataOut.iloc[i,5]=parr[i] From 7f76b28d00009fd1b402e426f7bbd81377f0a3f3 Mon Sep 17 00:00:00 2001 From: Phil McCown Date: Thu, 9 Nov 2017 16:08:07 -0500 Subject: [PATCH 07/15] Try parts of 2-1 --- Ex10Q2-1.py | 45 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100755 Ex10Q2-1.py diff --git a/Ex10Q2-1.py b/Ex10Q2-1.py new file mode 100755 index 0000000..ef06c79 --- /dev/null +++ b/Ex10Q2-1.py @@ -0,0 +1,45 @@ +import numpy +import pandas +import scipy +import scipy.integrate as spint +from plotnine import * +#SIR model in epidemiology +#N=S+I+R (susceptible-->B-->infected-->g-->resistant); unidirectional +#R0=((B*(S+I+R))/g) ; R0 is how bad the disease is in the population; see YY's notes for where max daily incidence, max daily prevalence, percent affected, and basic reproduction number are found in a graph + +#state vairables=I/S=state parameters=r/B +#basic reproduction only initial values + +#equations are dS/dt=-BIS, dI/dt=BIS-RI and dR/dt=RI +def epiSim(y,t,B,r): + I=y[0] + S=y[1] + dSdt= (-B*I*S) + dIdt= (B*I*S)-(r*I) + dRdt= (r*I) + return[dSdt,dIdt,dRdt] + +#generate a dataframe containing a list of the B and r parameters +parB=[0.0005,0.005,0.0001,0.00005,0.0001,0.0002,0.0001] +parr=[0.05,0.5,0.1,0.1,0.05,0.05,0.06] +DataOut=pandas.DataFrame(numpy.zeros([7,6]),columns=['MaxIncidence','MaxPrevalence','PercentAffected','R0','Beta','Gamma']) +times=range(0,500) +#extract params from dataframe +for i in range(0,6): + I0=[1] + S0=[999,1,0] + parms=(parB[i],parr[i]) + times=range(0,500) + sim=spint.odeint(func=epiSim,y0=S0,t=times,args=parms) + inc=sim[1:500,1]-sim[0:499,1] + maxinc=numpy.max(inc)#create a model using odeint + prevalence=((sim[:,1])/(sim[:,0]+sim[:,1]+sim[:,2])) + maxprev=numpy.max(prevalence) + affected=((sim[499,1]+sim[499,2])/(sim[499,0]+sim[499,1]+sim[499,2])) + R0=(parB[i])*1000/parr[i]/parr[i] + DataOut.iloc[i,0]=maxinc + DataOut.iloc[i,1]=maxprev + DataOut.iloc[i,2]=affected + DataOut.iloc[i,3]=R0 + DataOut.iloc[i,4]=parB[i] + DataOut.iloc[i,5]=parr[i] From 85f5d397877c5556fd848d545af77e64abcda3e8 Mon Sep 17 00:00:00 2001 From: kkilgoreND <31992787+kkilgoreND@users.noreply.github.com> Date: Thu, 9 Nov 2017 16:23:15 -0500 Subject: [PATCH 08/15] updated q2 --- Ex10Q2-1.py | 50 ++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 38 insertions(+), 12 deletions(-) diff --git a/Ex10Q2-1.py b/Ex10Q2-1.py index ef06c79..7f74274 100755 --- a/Ex10Q2-1.py +++ b/Ex10Q2-1.py @@ -30,16 +30,42 @@ def epiSim(y,t,B,r): S0=[999,1,0] parms=(parB[i],parr[i]) times=range(0,500) + #simuate model sim=spint.odeint(func=epiSim,y0=S0,t=times,args=parms) - inc=sim[1:500,1]-sim[0:499,1] - maxinc=numpy.max(inc)#create a model using odeint - prevalence=((sim[:,1])/(sim[:,0]+sim[:,1]+sim[:,2])) - maxprev=numpy.max(prevalence) - affected=((sim[499,1]+sim[499,2])/(sim[499,0]+sim[499,1]+sim[499,2])) - R0=(parB[i])*1000/parr[i]/parr[i] - DataOut.iloc[i,0]=maxinc - DataOut.iloc[i,1]=maxprev - DataOut.iloc[i,2]=affected - DataOut.iloc[i,3]=R0 - DataOut.iloc[i,4]=parB[i] - DataOut.iloc[i,5]=parr[i] + #put output into dataframe + sim_df=pandas.DataFrame({"t": times, "dSdt": sim[:,0], "dIdt": sim[:,1], "dRdt": sim[:,2]}) + #plot S, I, and R against time + plots[i]=ggplot(sim_df, aes(x="t"))+geom_line(aes(y="dSdt"))+theme_classic()+xlab("Time")+ylab("N")+geom_line(aes(y="dIdt"), color = "red")+geom_line(aes(y="dRdt"), color = "blue")+ggtitle(plotnames[i]) + #calculate maximum incidence + incidence=range(500) + #creates a list to store incidence values + for j in range(0, 500): + #calculate incidence at every time step + incidence[j]=sim_df.iloc[j,0]-sim_df.iloc[j-1,0] + DataOut.iloc[i, 1]=max(incidence) #get the max + #calculate max prevalence + prevalence=sim_df.iloc[:,0]/1000 + #calculate prevalence at every time step + DataOut.iloc[i, 2]=max(prevalence) + #calculate percent affected + DataOut.iloc[i, 3]=(sim_df.iloc[499,0]+sim_df.iloc[499,1])/100 + #calculate R0 + DataOut.iloc[i,4]=(betas[i]*1000)/gammas[i] + +### Look at the results! +print(plots) +print(DataOut) + + + #inc=sim[1:500,1]-sim[0:499,1] + #maxinc=numpy.max(inc)#create a model using odeint + #prevalence=((sim[:,1])/(sim[:,0]+sim[:,1]+sim[:,2])) + #maxprev=numpy.max(prevalence) + #affected=((sim[499,1]+sim[499,2])/(sim[499,0]+sim[499,1]+sim[499,2])) + #R0=(parB[i])*1000/parr[i]/parr[i] + #DataOut.iloc[i,0]=maxinc + #DataOut.iloc[i,1]=maxprev + #DataOut.iloc[i,2]=affected + #DataOut.iloc[i,3]=R0 + #DataOut.iloc[i,4]=parB[i] + #DataOut.iloc[i,5]=parr[i] From cf94a9b108ed79a73982c700e61a2a099869b32e Mon Sep 17 00:00:00 2001 From: kkilgoreND <31992787+kkilgoreND@users.noreply.github.com> Date: Thu, 9 Nov 2017 16:23:49 -0500 Subject: [PATCH 09/15] updated q2 --- Ex10Q2-1.py | 1 - 1 file changed, 1 deletion(-) diff --git a/Ex10Q2-1.py b/Ex10Q2-1.py index 7f74274..21667ad 100755 --- a/Ex10Q2-1.py +++ b/Ex10Q2-1.py @@ -52,7 +52,6 @@ def epiSim(y,t,B,r): #calculate R0 DataOut.iloc[i,4]=(betas[i]*1000)/gammas[i] -### Look at the results! print(plots) print(DataOut) From 90282f8c7a462faff3a84a9a78c3a36b95529dbf Mon Sep 17 00:00:00 2001 From: Phil McCown Date: Thu, 9 Nov 2017 16:35:20 -0500 Subject: [PATCH 10/15] GGplot not working, otherwise everything is fine, I think --- Ex10Q2-1.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Ex10Q2-1.py b/Ex10Q2-1.py index 21667ad..573da72 100755 --- a/Ex10Q2-1.py +++ b/Ex10Q2-1.py @@ -35,7 +35,7 @@ def epiSim(y,t,B,r): #put output into dataframe sim_df=pandas.DataFrame({"t": times, "dSdt": sim[:,0], "dIdt": sim[:,1], "dRdt": sim[:,2]}) #plot S, I, and R against time - plots[i]=ggplot(sim_df, aes(x="t"))+geom_line(aes(y="dSdt"))+theme_classic()+xlab("Time")+ylab("N")+geom_line(aes(y="dIdt"), color = "red")+geom_line(aes(y="dRdt"), color = "blue")+ggtitle(plotnames[i]) + ggplot(sim_df, aes(x="t"))+geom_line(aes(y="dSdt"))+theme_classic()+xlab("Time")+ylab("N")+geom_line(aes(y="dIdt"), color = "red")+geom_line(aes(y="dRdt"), color = "blue") #calculate maximum incidence incidence=range(500) #creates a list to store incidence values @@ -50,9 +50,9 @@ def epiSim(y,t,B,r): #calculate percent affected DataOut.iloc[i, 3]=(sim_df.iloc[499,0]+sim_df.iloc[499,1])/100 #calculate R0 - DataOut.iloc[i,4]=(betas[i]*1000)/gammas[i] + DataOut.iloc[i,4]=(parB[i]*1000)/parr[i] -print(plots) +print(ggplot) print(DataOut) From 62311251c19e86c6c5689f91d25ae2655dfb2914 Mon Sep 17 00:00:00 2001 From: kkilgoreND <31992787+kkilgoreND@users.noreply.github.com> Date: Thu, 9 Nov 2017 19:15:00 -0500 Subject: [PATCH 11/15] fixed q2 with forgotten list designations --- Ex10Q2-1.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/Ex10Q2-1.py b/Ex10Q2-1.py index 573da72..c26526b 100755 --- a/Ex10Q2-1.py +++ b/Ex10Q2-1.py @@ -22,6 +22,10 @@ def epiSim(y,t,B,r): #generate a dataframe containing a list of the B and r parameters parB=[0.0005,0.005,0.0001,0.00005,0.0001,0.0002,0.0001] parr=[0.05,0.5,0.1,0.1,0.05,0.05,0.06] +#create list for parameters +plotnames=["Parameter 1", "Parameter 2","Parameter 3","Parameter 4","Parameter 5","Parameter 6", "Parameter 7"] +#create list to hold plots +plots=["Plot1", "Plot2","Plot3","Plot4","Plot5","Plot6","Plot7"] DataOut=pandas.DataFrame(numpy.zeros([7,6]),columns=['MaxIncidence','MaxPrevalence','PercentAffected','R0','Beta','Gamma']) times=range(0,500) #extract params from dataframe @@ -35,7 +39,7 @@ def epiSim(y,t,B,r): #put output into dataframe sim_df=pandas.DataFrame({"t": times, "dSdt": sim[:,0], "dIdt": sim[:,1], "dRdt": sim[:,2]}) #plot S, I, and R against time - ggplot(sim_df, aes(x="t"))+geom_line(aes(y="dSdt"))+theme_classic()+xlab("Time")+ylab("N")+geom_line(aes(y="dIdt"), color = "red")+geom_line(aes(y="dRdt"), color = "blue") + ggplot(sim_df, aes(x="t"))+geom_line(aes(y="dSdt"))+theme_classic()+xlab("Time")+ylab("N")+geom_line(aes(y="dIdt"), color = "red")+geom_line(aes(y="dRdt"), color = "blue")+ggtitle(plotnames[i]) #calculate maximum incidence incidence=range(500) #creates a list to store incidence values @@ -52,7 +56,7 @@ def epiSim(y,t,B,r): #calculate R0 DataOut.iloc[i,4]=(parB[i]*1000)/parr[i] -print(ggplot) +print(plots) print(DataOut) From 1756c46039d1320a677c9aae63c6786ebc164246 Mon Sep 17 00:00:00 2001 From: Phil McCown Date: Thu, 9 Nov 2017 21:12:50 -0500 Subject: [PATCH 12/15] Added a Readme file for additional thoughts as required for Q2 --- Readme-for-thoughts-on-Q2.txt | 6 ++++++ 1 file changed, 6 insertions(+) create mode 100644 Readme-for-thoughts-on-Q2.txt diff --git a/Readme-for-thoughts-on-Q2.txt b/Readme-for-thoughts-on-Q2.txt new file mode 100644 index 0000000..cb26295 --- /dev/null +++ b/Readme-for-thoughts-on-Q2.txt @@ -0,0 +1,6 @@ +For Q2, we came to the conclusion that, with a high beta value within +the models provided, there is a large amount of incidence of infection, +as well as a higher rate of infection and higher rates of, ultimately, +resistant hosts. This makes conceptual sense since beta is used to +determine the rates at which individuals are converted from the +susceptible hosts pool into the infected hosts pool. From ee0aeb50670234951be23ad4c353a59aae81dd12 Mon Sep 17 00:00:00 2001 From: Phil McCown Date: Thu, 9 Nov 2017 21:13:50 -0500 Subject: [PATCH 13/15] Tidying up an extra file --- Ex10Q2-1.py | 74 ----------------------------------------------------- Ex10Q2.py | 53 +++++++++++++++++++++++++++++--------- 2 files changed, 41 insertions(+), 86 deletions(-) delete mode 100755 Ex10Q2-1.py diff --git a/Ex10Q2-1.py b/Ex10Q2-1.py deleted file mode 100755 index c26526b..0000000 --- a/Ex10Q2-1.py +++ /dev/null @@ -1,74 +0,0 @@ -import numpy -import pandas -import scipy -import scipy.integrate as spint -from plotnine import * -#SIR model in epidemiology -#N=S+I+R (susceptible-->B-->infected-->g-->resistant); unidirectional -#R0=((B*(S+I+R))/g) ; R0 is how bad the disease is in the population; see YY's notes for where max daily incidence, max daily prevalence, percent affected, and basic reproduction number are found in a graph - -#state vairables=I/S=state parameters=r/B -#basic reproduction only initial values - -#equations are dS/dt=-BIS, dI/dt=BIS-RI and dR/dt=RI -def epiSim(y,t,B,r): - I=y[0] - S=y[1] - dSdt= (-B*I*S) - dIdt= (B*I*S)-(r*I) - dRdt= (r*I) - return[dSdt,dIdt,dRdt] - -#generate a dataframe containing a list of the B and r parameters -parB=[0.0005,0.005,0.0001,0.00005,0.0001,0.0002,0.0001] -parr=[0.05,0.5,0.1,0.1,0.05,0.05,0.06] -#create list for parameters -plotnames=["Parameter 1", "Parameter 2","Parameter 3","Parameter 4","Parameter 5","Parameter 6", "Parameter 7"] -#create list to hold plots -plots=["Plot1", "Plot2","Plot3","Plot4","Plot5","Plot6","Plot7"] -DataOut=pandas.DataFrame(numpy.zeros([7,6]),columns=['MaxIncidence','MaxPrevalence','PercentAffected','R0','Beta','Gamma']) -times=range(0,500) -#extract params from dataframe -for i in range(0,6): - I0=[1] - S0=[999,1,0] - parms=(parB[i],parr[i]) - times=range(0,500) - #simuate model - sim=spint.odeint(func=epiSim,y0=S0,t=times,args=parms) - #put output into dataframe - sim_df=pandas.DataFrame({"t": times, "dSdt": sim[:,0], "dIdt": sim[:,1], "dRdt": sim[:,2]}) - #plot S, I, and R against time - ggplot(sim_df, aes(x="t"))+geom_line(aes(y="dSdt"))+theme_classic()+xlab("Time")+ylab("N")+geom_line(aes(y="dIdt"), color = "red")+geom_line(aes(y="dRdt"), color = "blue")+ggtitle(plotnames[i]) - #calculate maximum incidence - incidence=range(500) - #creates a list to store incidence values - for j in range(0, 500): - #calculate incidence at every time step - incidence[j]=sim_df.iloc[j,0]-sim_df.iloc[j-1,0] - DataOut.iloc[i, 1]=max(incidence) #get the max - #calculate max prevalence - prevalence=sim_df.iloc[:,0]/1000 - #calculate prevalence at every time step - DataOut.iloc[i, 2]=max(prevalence) - #calculate percent affected - DataOut.iloc[i, 3]=(sim_df.iloc[499,0]+sim_df.iloc[499,1])/100 - #calculate R0 - DataOut.iloc[i,4]=(parB[i]*1000)/parr[i] - -print(plots) -print(DataOut) - - - #inc=sim[1:500,1]-sim[0:499,1] - #maxinc=numpy.max(inc)#create a model using odeint - #prevalence=((sim[:,1])/(sim[:,0]+sim[:,1]+sim[:,2])) - #maxprev=numpy.max(prevalence) - #affected=((sim[499,1]+sim[499,2])/(sim[499,0]+sim[499,1]+sim[499,2])) - #R0=(parB[i])*1000/parr[i]/parr[i] - #DataOut.iloc[i,0]=maxinc - #DataOut.iloc[i,1]=maxprev - #DataOut.iloc[i,2]=affected - #DataOut.iloc[i,3]=R0 - #DataOut.iloc[i,4]=parB[i] - #DataOut.iloc[i,5]=parr[i] diff --git a/Ex10Q2.py b/Ex10Q2.py index ef06c79..c26526b 100755 --- a/Ex10Q2.py +++ b/Ex10Q2.py @@ -22,6 +22,10 @@ def epiSim(y,t,B,r): #generate a dataframe containing a list of the B and r parameters parB=[0.0005,0.005,0.0001,0.00005,0.0001,0.0002,0.0001] parr=[0.05,0.5,0.1,0.1,0.05,0.05,0.06] +#create list for parameters +plotnames=["Parameter 1", "Parameter 2","Parameter 3","Parameter 4","Parameter 5","Parameter 6", "Parameter 7"] +#create list to hold plots +plots=["Plot1", "Plot2","Plot3","Plot4","Plot5","Plot6","Plot7"] DataOut=pandas.DataFrame(numpy.zeros([7,6]),columns=['MaxIncidence','MaxPrevalence','PercentAffected','R0','Beta','Gamma']) times=range(0,500) #extract params from dataframe @@ -30,16 +34,41 @@ def epiSim(y,t,B,r): S0=[999,1,0] parms=(parB[i],parr[i]) times=range(0,500) + #simuate model sim=spint.odeint(func=epiSim,y0=S0,t=times,args=parms) - inc=sim[1:500,1]-sim[0:499,1] - maxinc=numpy.max(inc)#create a model using odeint - prevalence=((sim[:,1])/(sim[:,0]+sim[:,1]+sim[:,2])) - maxprev=numpy.max(prevalence) - affected=((sim[499,1]+sim[499,2])/(sim[499,0]+sim[499,1]+sim[499,2])) - R0=(parB[i])*1000/parr[i]/parr[i] - DataOut.iloc[i,0]=maxinc - DataOut.iloc[i,1]=maxprev - DataOut.iloc[i,2]=affected - DataOut.iloc[i,3]=R0 - DataOut.iloc[i,4]=parB[i] - DataOut.iloc[i,5]=parr[i] + #put output into dataframe + sim_df=pandas.DataFrame({"t": times, "dSdt": sim[:,0], "dIdt": sim[:,1], "dRdt": sim[:,2]}) + #plot S, I, and R against time + ggplot(sim_df, aes(x="t"))+geom_line(aes(y="dSdt"))+theme_classic()+xlab("Time")+ylab("N")+geom_line(aes(y="dIdt"), color = "red")+geom_line(aes(y="dRdt"), color = "blue")+ggtitle(plotnames[i]) + #calculate maximum incidence + incidence=range(500) + #creates a list to store incidence values + for j in range(0, 500): + #calculate incidence at every time step + incidence[j]=sim_df.iloc[j,0]-sim_df.iloc[j-1,0] + DataOut.iloc[i, 1]=max(incidence) #get the max + #calculate max prevalence + prevalence=sim_df.iloc[:,0]/1000 + #calculate prevalence at every time step + DataOut.iloc[i, 2]=max(prevalence) + #calculate percent affected + DataOut.iloc[i, 3]=(sim_df.iloc[499,0]+sim_df.iloc[499,1])/100 + #calculate R0 + DataOut.iloc[i,4]=(parB[i]*1000)/parr[i] + +print(plots) +print(DataOut) + + + #inc=sim[1:500,1]-sim[0:499,1] + #maxinc=numpy.max(inc)#create a model using odeint + #prevalence=((sim[:,1])/(sim[:,0]+sim[:,1]+sim[:,2])) + #maxprev=numpy.max(prevalence) + #affected=((sim[499,1]+sim[499,2])/(sim[499,0]+sim[499,1]+sim[499,2])) + #R0=(parB[i])*1000/parr[i]/parr[i] + #DataOut.iloc[i,0]=maxinc + #DataOut.iloc[i,1]=maxprev + #DataOut.iloc[i,2]=affected + #DataOut.iloc[i,3]=R0 + #DataOut.iloc[i,4]=parB[i] + #DataOut.iloc[i,5]=parr[i] From 45c22f6ae76bdf25e0fa33dbe444f649bc517279 Mon Sep 17 00:00:00 2001 From: kkilgoreND <31992787+kkilgoreND@users.noreply.github.com> Date: Thu, 9 Nov 2017 21:19:19 -0500 Subject: [PATCH 14/15] final updates for q2 --- Ex10Q2.py | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/Ex10Q2.py b/Ex10Q2.py index c26526b..92ccb41 100755 --- a/Ex10Q2.py +++ b/Ex10Q2.py @@ -59,16 +59,7 @@ def epiSim(y,t,B,r): print(plots) print(DataOut) - - #inc=sim[1:500,1]-sim[0:499,1] - #maxinc=numpy.max(inc)#create a model using odeint - #prevalence=((sim[:,1])/(sim[:,0]+sim[:,1]+sim[:,2])) - #maxprev=numpy.max(prevalence) - #affected=((sim[499,1]+sim[499,2])/(sim[499,0]+sim[499,1]+sim[499,2])) - #R0=(parB[i])*1000/parr[i]/parr[i] - #DataOut.iloc[i,0]=maxinc - #DataOut.iloc[i,1]=maxprev - #DataOut.iloc[i,2]=affected - #DataOut.iloc[i,3]=R0 - #DataOut.iloc[i,4]=parB[i] - #DataOut.iloc[i,5]=parr[i] +#Gleaned information: +#increases in the gamma value appear to have a inverse-correlation with the calculated values +#increases in the beta values appear to have a direct-correlation with the calculated values +#Thus...increase beta to increase the spread of the disease From 718b058c940c98e62590a92f0fc438604c9d8e7e Mon Sep 17 00:00:00 2001 From: kkilgoreND <31992787+kkilgoreND@users.noreply.github.com> Date: Thu, 9 Nov 2017 21:22:35 -0500 Subject: [PATCH 15/15] Full script for exercise 10 --- Ex10FullScript.py | 116 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 116 insertions(+) create mode 100644 Ex10FullScript.py diff --git a/Ex10FullScript.py b/Ex10FullScript.py new file mode 100644 index 0000000..c720a7e --- /dev/null +++ b/Ex10FullScript.py @@ -0,0 +1,116 @@ +#---Part 1--- +import numpy +import pandas +import scipy +import scipy.integrate as spint +from plotnine import * +def ddSim(y,t0,r,K): + N=y[0] + dNdt=r*(1-N/K)*N + return [dNdt] +times=(range(0,1000)) + +rs=[-0.1,0.1,0.4,0.8,1] +store_rs=pandas.DataFrame({"time":times,"rl":0,"r2":0,"r3":0,"r4":0,"r5":0}) +for i in range(0,len(rs)): + K=100 + N=10 + pars1=(rs[i],K) + sim=spint.odeint(func=ddSim,y0=N,t=times,args=pars1) + store_rs.iloc[:,i]=sim[:,0] +a=ggplot(store_rs,aes(x="time",y="r2")) +a+geom_point()+coord_cartesian() + +ks=[10,50,100] +store_ks=pandas.DataFrame({"time":times,"kl":0,"k2":0,"k3":0}) +for i in range(0,len(ks)): + N=1 + pars2=(0.2,ks[i]) + sim2=spint.odeint(func=ddSim,y0=N,t=times,args=pars2) + store_ks.iloc[:,i]=sim2[:,0] + +ns=[1,50,100] +store_ns=pandas.DataFrame({"time":times,"nl":0,"n2":0,"n3":0}) +for i in range(0,len(ns)): + N=1,50,100 + pars3=(0.1,50) + sim3=spint.odeint(func=ddSim,y0=N[i],t=times,args=pars3) + store_ns.iloc[:,i]=sim3[:,0] + +#Placing models into dataframe +modelOutputRS=pandas.DataFrame({"t":times,"N":sim[:,0]}) +modelOutputKS=pandas.DataFrame({"t":times,"N":sim2[:,0]}) +modelOutputNS=pandas.DataFrame({"t":times,"N":sim3[:,0]}) + +#Plot simulation outputs into ggplot +ggplot(modelOutputKS,aes(x="t",y="N"))+geom_line()+theme_classic() +ggplot(modelOutputRS,aes(x="t",y="N"))+geom_line()+theme_classic() +ggplot(modelOutputNS,aes(x="t",y="N"))+geom_line()+theme_classic() + + +#---Part 2--- +import numpy +import pandas +import scipy +import scipy.integrate as spint +from plotnine import * +#SIR model in epidemiology +#N=S+I+R (susceptible-->B-->infected-->g-->resistant); unidirectional +#R0=((B*(S+I+R))/g) ; R0 is how bad the disease is in the population; see YY's notes for where max daily incidence, max daily prevalence, percent affected, and basic reproduction number are found in a graph + +#state vairables=I/S=state parameters=r/B +#basic reproduction only initial values + +#equations are dS/dt=-BIS, dI/dt=BIS-RI and dR/dt=RI +def epiSim(y,t,B,r): + I=y[0] + S=y[1] + dSdt= (-B*I*S) + dIdt= (B*I*S)-(r*I) + dRdt= (r*I) + return[dSdt,dIdt,dRdt] + +#generate a dataframe containing a list of the B and r parameters +parB=[0.0005,0.005,0.0001,0.00005,0.0001,0.0002,0.0001] +parr=[0.05,0.5,0.1,0.1,0.05,0.05,0.06] +#create list for parameters +plotnames=["Parameter 1", "Parameter 2","Parameter 3","Parameter 4","Parameter 5","Parameter 6", "Parameter 7"] +#create list to hold plots +plots=["Plot1", "Plot2","Plot3","Plot4","Plot5","Plot6","Plot7"] +DataOut=pandas.DataFrame(numpy.zeros([7,6]),columns=['MaxIncidence','MaxPrevalence','PercentAffected','R0','Beta','Gamma']) +times=range(0,500) +#extract params from dataframe +for i in range(0,6): + I0=[1] + S0=[999,1,0] + parms=(parB[i],parr[i]) + times=range(0,500) + #simuate model + sim=spint.odeint(func=epiSim,y0=S0,t=times,args=parms) + #put output into dataframe + sim_df=pandas.DataFrame({"t": times, "dSdt": sim[:,0], "dIdt": sim[:,1], "dRdt": sim[:,2]}) + #plot S, I, and R against time + ggplot(sim_df, aes(x="t"))+geom_line(aes(y="dSdt"))+theme_classic()+xlab("Time")+ylab("N")+geom_line(aes(y="dIdt"), color = "red")+geom_line(aes(y="dRdt"), color = "blue")+ggtitle(plotnames[i]) + #calculate maximum incidence + incidence=range(500) + #creates a list to store incidence values + for j in range(0, 500): + #calculate incidence at every time step + incidence[j]=sim_df.iloc[j,0]-sim_df.iloc[j-1,0] + DataOut.iloc[i, 1]=max(incidence) #get the max + #calculate max prevalence + prevalence=sim_df.iloc[:,0]/1000 + #calculate prevalence at every time step + DataOut.iloc[i, 2]=max(prevalence) + #calculate percent affected + DataOut.iloc[i, 3]=(sim_df.iloc[499,0]+sim_df.iloc[499,1])/100 + #calculate R0 + DataOut.iloc[i,4]=(parB[i]*1000)/parr[i] + +print(plots) +print(DataOut) + +#Gleaned information: +#increases in the gamma value appear to have a inverse-correlation with the calculated values +#increases in the beta values appear to have a direct-correlation with the calculated values +#Thus...increase beta to increase the spread of the disease