From 40b47d52f25d6d2f5f264226fe728308b32d8d19 Mon Sep 17 00:00:00 2001 From: ayamasaki2011 Date: Fri, 3 Nov 2017 10:39:23 -0400 Subject: [PATCH 01/12] AEY: Push initial file for Question 1 --- Question1.py | 4 ++++ 1 file changed, 4 insertions(+) create mode 100755 Question1.py diff --git a/Question1.py b/Question1.py new file mode 100755 index 0000000..ff1979a --- /dev/null +++ b/Question1.py @@ -0,0 +1,4 @@ +import pandas +import scipy +import scipy.integrate as spint +from plotnine import * \ No newline at end of file From e8a02adb79f1c578c4a852114b09773a89cd9925 Mon Sep 17 00:00:00 2001 From: ayamasaki2011 Date: Fri, 3 Nov 2017 11:00:06 -0400 Subject: [PATCH 02/12] AEY: Add for loop for varying growth rates and store in data frame --- Question1.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/Question1.py b/Question1.py index ff1979a..f4b8a9a 100755 --- a/Question1.py +++ b/Question1.py @@ -1,4 +1,22 @@ import pandas import scipy import scipy.integrate as spint -from plotnine import * \ No newline at end of file +from plotnine import * + +def ddSim (y,t0,r,K): + N=y[0] + dNdt=r*(1-N/K)*N + return [dNdt] + +K=100 +N0=[10] +times=range(0,500) + +rs=[-0.1,0.1,0.4,0.8,1.0] +store_rs=pandas.DataFrame({"time":times,"r1":0,"r2":0,"r3":0,"r4":0,"r5":0}) + +for i in range(0,len(rs)): + params=(rs[i],K) + sim=spint.odeint(func=ddSim,y0=N0,t=times,args=params) + store_rs.iloc[:,i]=sim[:,0] + From b3b6f076bb085954497098950c29a1bacfa1855b Mon Sep 17 00:00:00 2001 From: ayamasaki2011 Date: Fri, 3 Nov 2017 11:10:39 -0400 Subject: [PATCH 03/12] AEY: Complete for loop and plot for varying growth rates - should probably figure out how to plot with multiple colors and a legend --- Question1.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/Question1.py b/Question1.py index f4b8a9a..53959eb 100755 --- a/Question1.py +++ b/Question1.py @@ -20,3 +20,11 @@ def ddSim (y,t0,r,K): sim=spint.odeint(func=ddSim,y0=N0,t=times,args=params) store_rs.iloc[:,i]=sim[:,0] +rates=ggplot(store_rs,aes(x="time",y="N")) +rates=rates+geom_line(store_rs,aes(y="r1")) +rates=rates+geom_line(store_rs,aes(y="r2")) +rates=rates+geom_line(store_rs,aes(y="r3")) +rates=rates+geom_line(store_rs,aes(y="r4")) +rates=rates+geom_line(store_rs,aes(y="r5")) + +rates \ No newline at end of file From 0e23030781c95906d2a877731542cf2409d66bdf Mon Sep 17 00:00:00 2001 From: ayamasaki2011 Date: Fri, 3 Nov 2017 11:22:35 -0400 Subject: [PATCH 04/12] AEY: Complete for loop and plot for varying carrying capacities - see previous point about aesthetics of graph --- Question1.py | 41 ++++++++++++++++++++++++++++++----------- 1 file changed, 30 insertions(+), 11 deletions(-) diff --git a/Question1.py b/Question1.py index 53959eb..286a89c 100755 --- a/Question1.py +++ b/Question1.py @@ -3,28 +3,47 @@ import scipy.integrate as spint from plotnine import * -def ddSim (y,t0,r,K): +def ddSim (y,t0,r,K): #Define basic model for density-dependent growth N=y[0] dNdt=r*(1-N/K)*N return [dNdt] -K=100 +K=100 #Define a set carrying capacity & initial population size N0=[10] times=range(0,500) -rs=[-0.1,0.1,0.4,0.8,1.0] -store_rs=pandas.DataFrame({"time":times,"r1":0,"r2":0,"r3":0,"r4":0,"r5":0}) - -for i in range(0,len(rs)): - params=(rs[i],K) +rs=[-0.1,0.1,0.4,0.8,1.0] #Define a list of growth rates for each individual population +store_rs=pandas.DataFrame({"time":times,"r1":0,"r2":0,"r3":0,"r4":0,"r5":0}) #Initialize a data frame to store population size, with + #columns labeled for times, and each growth rate +for i in range(0,len(rs)): #Set up a for loop that runs through each growth rate in the list, and models population growth using + params=(rs[i],K) #established parameters sim=spint.odeint(func=ddSim,y0=N0,t=times,args=params) - store_rs.iloc[:,i]=sim[:,0] + store_rs.iloc[:,i]=sim[:,0] #Store population values in the column appropriate for the current growth rate being modeled -rates=ggplot(store_rs,aes(x="time",y="N")) -rates=rates+geom_line(store_rs,aes(y="r1")) +rates=ggplot(store_rs,aes(x="time",y="N"))+theme_classic() #Set up a base plot for the population curves +rates=rates+geom_line(store_rs,aes(y="r1")) #Add each population curve to the graph rates=rates+geom_line(store_rs,aes(y="r2")) rates=rates+geom_line(store_rs,aes(y="r3")) rates=rates+geom_line(store_rs,aes(y="r4")) rates=rates+geom_line(store_rs,aes(y="r5")) -rates \ No newline at end of file +rates #Display the graph + +r=0.2 #Define a set growth rate and iitial population size +N0=[1] +times=range(0,500) + +Ks=[10,50,100] #Define a list of carrying capacities for each individual population +store_Ks=pandas.DataFrame({"time":times,"K1":0,"K2":0,"K3":0}) #Initialize a data frame to store population size with columns for each + #carrying capacity, and the times modeled +for i in range(0,len(Ks)): #Set up a for loop that runs through each carrying capacity in the list and models population growth using + params=(r,Ks[i]) #other established parameters + sim=spint.odeint(func=ddSim,y0=N0,t=times,args=params) + store_Ks.iloc[:,i]=sim[:,0] #Store population values in the column appropriate for the current carrying capacity + +capacity=ggplot(store_Ks,aes(x="time",y="N"))+theme_classic()#Set up base plot for population curves +capacity=capacity+geom_line(store_Ks,aes(y="K1")) +capacity=capacity+geom_line(store_Ks,aes(y="K2")) +capacity=capacity+geom_line(store_Ks,aes(y="K3")) + +capacity #Display graph \ No newline at end of file From b1d94762e1597018bf3de3c7156301510df26916 Mon Sep 17 00:00:00 2001 From: Joe C Date: Wed, 8 Nov 2017 10:06:27 -0500 Subject: [PATCH 05/12] I tried answering part 3 of question 1, but I am clearly missing some part as it says it is 2D and the graph needs to be 1D --- Ex10_Answers.py | 70 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) create mode 100755 Ex10_Answers.py diff --git a/Ex10_Answers.py b/Ex10_Answers.py new file mode 100755 index 0000000..963ea85 --- /dev/null +++ b/Ex10_Answers.py @@ -0,0 +1,70 @@ +import pandas +import scipy +import scipy.integrate as spint +from plotnine import * + +def ddSim (y,t0,r,K): #Define basic model for density-dependent growth + N=y[0] + dNdt=r*(1-N/K)*N + return [dNdt] + +K=100 #Define a set carrying capacity & initial population size +N0=[10] +times=range(0,500) + +rs=[-0.1,0.1,0.4,0.8,1.0] #Define a list of growth rates for each individual population +store_rs=pandas.DataFrame({"time":times,"r1":0,"r2":0,"r3":0,"r4":0,"r5":0}) #Initialize a data frame to store population size, with + #columns labeled for times, and each growth rate +for i in range(0,len(rs)): #Set up a for loop that runs through each growth rate in the list, and models population growth using + params=(rs[i],K) #established parameters + sim=spint.odeint(func=ddSim,y0=N0,t=times,args=params) + store_rs.iloc[:,i]=sim[:,0] #Store population values in the column appropriate for the current growth rate being modeled + +rates=ggplot(store_rs,aes(x="time",y="N"))+theme_classic() #Set up a base plot for the population curves +rates=rates+geom_line(store_rs,aes(y="r1")) #Add each population curve to the graph +rates=rates+geom_line(store_rs,aes(y="r2")) +rates=rates+geom_line(store_rs,aes(y="r3")) +rates=rates+geom_line(store_rs,aes(y="r4")) +rates=rates+geom_line(store_rs,aes(y="r5")) + +rates #Display the graph + +r=0.2 #Define a set growth rate and iitial population size +N0=[1] +times=range(0,500) + +Ks=[10,50,100] #Define a list of carrying capacities for each individual population +store_Ks=pandas.DataFrame({"time":times,"K1":0,"K2":0,"K3":0}) #Initialize a data frame to store population size with columns for each + #carrying capacity, and the times modeled +for i in range(0,len(Ks)): #Set up a for loop that runs through each carrying capacity in the list and models population growth using + params=(r,Ks[i]) #other established parameters + sim=spint.odeint(func=ddSim,y0=N0,t=times,args=params) + store_Ks.iloc[:,i]=sim[:,0] #Store population values in the column appropriate for the current carrying capacity + +capacity=ggplot(store_Ks,aes(x="time",y="N"))+theme_classic()#Set up base plot for population curves +capacity=capacity+geom_line(store_Ks,aes(y="K1")) +capacity=capacity+geom_line(store_Ks,aes(y="K2")) +capacity=capacity+geom_line(store_Ks,aes(y="K3")) + +capacity #Display graph +## number 1 part 3 + +r=0.1 #Define a set growth rate and carrying capacity +K=[50] +times=range(0,500) + +N0s=[1,50,100] #Define a list of initial sizes for each individual population +store_N0s=pandas.DataFrame({"time":times,"N01":0,"N02":0,"N03":0}) #Initialize a data frame to store initial sizes with columns for each + #initial size, and the times modeled +for i in range(0,len(N0s)): #Set up a for loop that runs through each N0 in the list and models population growth using + params=(r,K) #other established parameters + sim=spint.odeint(func=ddSim,y0=N0,t=times,args=params) + store_N0s.iloc[:,i]=sim[:,0] #Store population values in the column appropriate for the current carrying capacity + +capacity=ggplot(store_N0s,aes(x="time",y="N"))+theme_classic()#Set up base plot for population curves +capacity=capacity+geom_line(store_N0s,aes(y="N01")) +capacity=capacity+geom_line(store_N0s,aes(y="N02")) +capacity=capacity+geom_line(store_N0s,aes(y="N03")) + +capacity #Display graph + From bed4719dd1955144bdf6c0660b8c6f860152ed4d Mon Sep 17 00:00:00 2001 From: Joe C Date: Wed, 8 Nov 2017 21:37:39 -0500 Subject: [PATCH 06/12] attempts made on question 2 but keep getting the same error at odeint step --- Ex10_Answers.py | 97 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 96 insertions(+), 1 deletion(-) diff --git a/Ex10_Answers.py b/Ex10_Answers.py index 963ea85..293c8c4 100755 --- a/Ex10_Answers.py +++ b/Ex10_Answers.py @@ -54,7 +54,7 @@ def ddSim (y,t0,r,K): #Define basic model for density-dependent growth times=range(0,500) N0s=[1,50,100] #Define a list of initial sizes for each individual population -store_N0s=pandas.DataFrame({"time":times,"N01":0,"N02":0,"N03":0}) #Initialize a data frame to store initial sizes with columns for each +store_N0s=pandas.DataFrame({"time":times,"N01":1,"N02":50,"N03":1500}) #Initialize a data frame to store initial sizes with columns for each #initial size, and the times modeled for i in range(0,len(N0s)): #Set up a for loop that runs through each N0 in the list and models population growth using params=(r,K) #other established parameters @@ -68,3 +68,98 @@ def ddSim (y,t0,r,K): #Define basic model for density-dependent growth capacity #Display graph +# Question 2 + + +import pandas +import scipy +import scipy.integrate as spint +from plotnine import * +import numpy as np +from scipy.integrate import odeint + + +N=1000 +S=999 +I=1 +R=0 +T=range(0,500) + +B=[.0005] +gamma=[.05] + +params= (N, B, gamma) + +def ddSim(y, t, N, B, gamma): + y = S, I, R + dSdt=(-B*I*S) + dIdt=(B*I*S-(gamma*I)) + dRdt=(gamma*I) + return (dSdt, dIdt, dRdt) + +y0= S, I, R +ret=odeint(func=ddSim,y0=N,t=T,args=params) +S, I, R =ret.T + +fig = plt.figure(facecolor='w') +ax = fig.add_subplot(111, axis_bgcolor='#dddddd', axisbelow=True) +ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible') +ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected') +ax.plot(t, R/1000, 'g', alpha=0.5, lw=2, label='Recovered with immunity') +ax.set_xlabel('Time /days') +ax.set_ylabel('Number (1000s)') +ax.set_ylim(0,1.2) +ax.yaxis.set_tick_params(length=0) +ax.xaxis.set_tick_params(length=0) +ax.grid(b=True, which='major', c='w', lw=2, ls='-') +legend = ax.legend() +legend.get_frame().set_alpha(0.5) +for spine in ('top', 'right', 'bottom', 'left'): + ax.spines[spine].set_visible(False) +plt.show() + + + +incidence=(It0-It1) +prevalence=(I/(S+I+R)) +percentAffected=((I+R)/(S+I+R)) +basicReproductionNumber=(((B(S+I+R))/gamma)) + + + +######################################################### +import scipy +import scipy.integrate as spint +import numpy +import matplotlib.pyplot as plt + +def SIR_model(y, t, beta, gamma): + S, I, R = y + + dSdt=(-beta*I*S) + dIdt=(beta*I*S-(gamma*I)) + dRdt=(gamma*I) + + return (dSdt, dIdt, dRdt) + +N=1000 +S0=999 +I0=1 +R0=0 +beta=[.0005] +gamma0=[.05] +t=numpy.linspace(0,500) + +solution=scipy.integrate.odeint(SIR_model, [S0, I0, R0], t, args=(beta, gamma)) +solution= numpy.array(solution) + +plt.figure(figsize=[6,4]) +plt.plot(t, solution[:, 0], label="S(t)") +plt.plot(t, solution[:, 1], label="I(t)") +plt.plot(t, solution[:, 2], label="R(t)") +plt.show() + + + + + From 031e1d50a95f591f624ac48cda81ed0c127a2ac3 Mon Sep 17 00:00:00 2001 From: Joe C Date: Thu, 9 Nov 2017 14:56:35 -0500 Subject: [PATCH 07/12] got question 2 to work but only for one at a time. still cant get 1c. --- Ex10_Answers.py | 66 ++++++------------------------------------------- 1 file changed, 7 insertions(+), 59 deletions(-) diff --git a/Ex10_Answers.py b/Ex10_Answers.py index 293c8c4..86420ac 100755 --- a/Ex10_Answers.py +++ b/Ex10_Answers.py @@ -70,64 +70,10 @@ def ddSim (y,t0,r,K): #Define basic model for density-dependent growth # Question 2 +######################################################### -import pandas -import scipy -import scipy.integrate as spint -from plotnine import * -import numpy as np -from scipy.integrate import odeint - - -N=1000 -S=999 -I=1 -R=0 -T=range(0,500) - -B=[.0005] -gamma=[.05] - -params= (N, B, gamma) - -def ddSim(y, t, N, B, gamma): - y = S, I, R - dSdt=(-B*I*S) - dIdt=(B*I*S-(gamma*I)) - dRdt=(gamma*I) - return (dSdt, dIdt, dRdt) - -y0= S, I, R -ret=odeint(func=ddSim,y0=N,t=T,args=params) -S, I, R =ret.T - -fig = plt.figure(facecolor='w') -ax = fig.add_subplot(111, axis_bgcolor='#dddddd', axisbelow=True) -ax.plot(t, S/1000, 'b', alpha=0.5, lw=2, label='Susceptible') -ax.plot(t, I/1000, 'r', alpha=0.5, lw=2, label='Infected') -ax.plot(t, R/1000, 'g', alpha=0.5, lw=2, label='Recovered with immunity') -ax.set_xlabel('Time /days') -ax.set_ylabel('Number (1000s)') -ax.set_ylim(0,1.2) -ax.yaxis.set_tick_params(length=0) -ax.xaxis.set_tick_params(length=0) -ax.grid(b=True, which='major', c='w', lw=2, ls='-') -legend = ax.legend() -legend.get_frame().set_alpha(0.5) -for spine in ('top', 'right', 'bottom', 'left'): - ax.spines[spine].set_visible(False) -plt.show() - - - -incidence=(It0-It1) -prevalence=(I/(S+I+R)) -percentAffected=((I+R)/(S+I+R)) -basicReproductionNumber=(((B(S+I+R))/gamma)) - - +# this works but only for one beta and gamma at a time -######################################################### import scipy import scipy.integrate as spint import numpy @@ -146,13 +92,15 @@ def SIR_model(y, t, beta, gamma): S0=999 I0=1 R0=0 -beta=[.0005] -gamma0=[.05] +beta=.0005 +gamma=.05 t=numpy.linspace(0,500) -solution=scipy.integrate.odeint(SIR_model, [S0, I0, R0], t, args=(beta, gamma)) +solution=scipy.integrate.odeint(SIR_model, [S0, I0, R0], t, args=(.005, .5)) solution= numpy.array(solution) + + plt.figure(figsize=[6,4]) plt.plot(t, solution[:, 0], label="S(t)") plt.plot(t, solution[:, 1], label="I(t)") From 65db75f60b55d3cce9ae88ea3dae4cdb2a9e162b Mon Sep 17 00:00:00 2001 From: ayamasaki2011 Date: Thu, 9 Nov 2017 18:01:35 -0500 Subject: [PATCH 08/12] AEY: Got Question 1 Part 3 working --- Question1.py | 28 +++++++++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/Question1.py b/Question1.py index 286a89c..f48946d 100755 --- a/Question1.py +++ b/Question1.py @@ -29,9 +29,8 @@ def ddSim (y,t0,r,K): #Define basic model for density-dependent growth rates #Display the graph -r=0.2 #Define a set growth rate and iitial population size +r=0.2 #Define a set growth rate and initial population size N0=[1] -times=range(0,500) Ks=[10,50,100] #Define a list of carrying capacities for each individual population store_Ks=pandas.DataFrame({"time":times,"K1":0,"K2":0,"K3":0}) #Initialize a data frame to store population size with columns for each @@ -46,4 +45,27 @@ def ddSim (y,t0,r,K): #Define basic model for density-dependent growth capacity=capacity+geom_line(store_Ks,aes(y="K2")) capacity=capacity+geom_line(store_Ks,aes(y="K3")) -capacity #Display graph \ No newline at end of file +capacity #Display graph + +r=0.1 #Define a set growth rate and carrying capacity +K=50 +times=range(0,500) + +N0s=[1,50,100] #Define a list of initial sizes for each individual population +store_N0s=pandas.DataFrame({"time":times,"N01":1,"N02":50,"N03":1500}) #Initialize a data frame to store initial sizes with columns for each + #initial size, and the times modeled + +for i in range(0,len(N0s)): #Set up a for loop that runs through each N0 in the list and models population growth using + params=(r,K)#other established parameters + N0=N0s[i] #Assign N0 a single value from the array of initial population sizes + sim=spint.odeint(func=ddSim,y0=N0,t=times,args=params) + store_N0s.iloc[:,i]=sim[:,0] #Store population values in the column appropriate for the current carrying capacity + +initpop=ggplot(store_N0s,aes(x="time",y="N"))+theme_classic()#Set up base plot for population curves +initpop=initpop+geom_line(store_N0s,aes(y="N01")) +initpop=initpop+geom_line(store_N0s,aes(y="N02")) +initpop=initpop+geom_line(store_N0s,aes(y="N03")) + +initpop #Display graph + + From b21910c7a684b055c708f5c09ff3452a80168f2c Mon Sep 17 00:00:00 2001 From: ayamasaki2011 Date: Thu, 9 Nov 2017 18:03:17 -0500 Subject: [PATCH 09/12] AEY: Add completed Q1 Part 3 to final answers file --- Ex10_Answers.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/Ex10_Answers.py b/Ex10_Answers.py index 86420ac..d7f6448 100755 --- a/Ex10_Answers.py +++ b/Ex10_Answers.py @@ -50,23 +50,25 @@ def ddSim (y,t0,r,K): #Define basic model for density-dependent growth ## number 1 part 3 r=0.1 #Define a set growth rate and carrying capacity -K=[50] +K=50 times=range(0,500) N0s=[1,50,100] #Define a list of initial sizes for each individual population store_N0s=pandas.DataFrame({"time":times,"N01":1,"N02":50,"N03":1500}) #Initialize a data frame to store initial sizes with columns for each #initial size, and the times modeled + for i in range(0,len(N0s)): #Set up a for loop that runs through each N0 in the list and models population growth using - params=(r,K) #other established parameters + params=(r,K)#other established parameters + N0=N0s[i] #Assign N0 a single value from the array of initial population sizes sim=spint.odeint(func=ddSim,y0=N0,t=times,args=params) store_N0s.iloc[:,i]=sim[:,0] #Store population values in the column appropriate for the current carrying capacity -capacity=ggplot(store_N0s,aes(x="time",y="N"))+theme_classic()#Set up base plot for population curves -capacity=capacity+geom_line(store_N0s,aes(y="N01")) -capacity=capacity+geom_line(store_N0s,aes(y="N02")) -capacity=capacity+geom_line(store_N0s,aes(y="N03")) +initpop=ggplot(store_N0s,aes(x="time",y="N"))+theme_classic()#Set up base plot for population curves +initpop=initpop+geom_line(store_N0s,aes(y="N01")) +initpop=initpop+geom_line(store_N0s,aes(y="N02")) +initpop=initpop+geom_line(store_N0s,aes(y="N03")) -capacity #Display graph +initpop #Display graph # Question 2 From ad27c5c28c2184752c905b64b82ce509d62bcec5 Mon Sep 17 00:00:00 2001 From: Joe C Date: Thu, 9 Nov 2017 20:05:41 -0500 Subject: [PATCH 10/12] question 2 completed --- Ex10_Answers.py | 72 ++++++++++++++++++++++++++++++++++--------------- 1 file changed, 51 insertions(+), 21 deletions(-) diff --git a/Ex10_Answers.py b/Ex10_Answers.py index 86420ac..c75359b 100755 --- a/Ex10_Answers.py +++ b/Ex10_Answers.py @@ -10,7 +10,7 @@ def ddSim (y,t0,r,K): #Define basic model for density-dependent growth K=100 #Define a set carrying capacity & initial population size N0=[10] -times=range(0,500) +times=range(0,100) rs=[-0.1,0.1,0.4,0.8,1.0] #Define a list of growth rates for each individual population store_rs=pandas.DataFrame({"time":times,"r1":0,"r2":0,"r3":0,"r4":0,"r5":0}) #Initialize a data frame to store population size, with @@ -31,7 +31,7 @@ def ddSim (y,t0,r,K): #Define basic model for density-dependent growth r=0.2 #Define a set growth rate and iitial population size N0=[1] -times=range(0,500) +times=range(0,50) Ks=[10,50,100] #Define a list of carrying capacities for each individual population store_Ks=pandas.DataFrame({"time":times,"K1":0,"K2":0,"K3":0}) #Initialize a data frame to store population size with columns for each @@ -50,37 +50,37 @@ def ddSim (y,t0,r,K): #Define basic model for density-dependent growth ## number 1 part 3 r=0.1 #Define a set growth rate and carrying capacity -K=[50] -times=range(0,500) +K=50 +times=range(0,100) N0s=[1,50,100] #Define a list of initial sizes for each individual population -store_N0s=pandas.DataFrame({"time":times,"N01":1,"N02":50,"N03":1500}) #Initialize a data frame to store initial sizes with columns for each +store_N0s=pandas.DataFrame({"time":times,"N01":0,"N02":0,"N03":0}) #Initialize a data frame to store initial sizes with columns for each #initial size, and the times modeled for i in range(0,len(N0s)): #Set up a for loop that runs through each N0 in the list and models population growth using params=(r,K) #other established parameters - sim=spint.odeint(func=ddSim,y0=N0,t=times,args=params) + sim=spint.odeint(func=ddSim,y0=N0s[i],t=times,args=params) store_N0s.iloc[:,i]=sim[:,0] #Store population values in the column appropriate for the current carrying capacity -capacity=ggplot(store_N0s,aes(x="time",y="N"))+theme_classic()#Set up base plot for population curves -capacity=capacity+geom_line(store_N0s,aes(y="N01")) -capacity=capacity+geom_line(store_N0s,aes(y="N02")) -capacity=capacity+geom_line(store_N0s,aes(y="N03")) +popsize=ggplot(store_N0s,aes(x="time",y="N"))+theme_classic()#Set up base plot for population curves +popsize=popsize+geom_line(store_N0s,aes(y="N01")) +popsize=popsize+geom_line(store_N0s,aes(y="N02")) +popsize=popsize+geom_line(store_N0s,aes(y="N03")) -capacity #Display graph +popsize #Display graph # Question 2 ######################################################### -# this works but only for one beta and gamma at a time - import scipy import scipy.integrate as spint import numpy import matplotlib.pyplot as plt def SIR_model(y, t, beta, gamma): - S, I, R = y + S = y[0] + I = y[1] + R = y[2] dSdt=(-beta*I*S) dIdt=(beta*I*S-(gamma*I)) @@ -88,17 +88,25 @@ def SIR_model(y, t, beta, gamma): return (dSdt, dIdt, dRdt) -N=1000 +N = [999, 1, 0] S0=999 I0=1 R0=0 -beta=.0005 -gamma=.05 -t=numpy.linspace(0,500) +beta=numpy.array([.0005, .005, .0001, .00005, .0001, .0002, .0001]) +gamma=numpy.array([.05, .5, .1, .1, .05, .05, .06]) + +store_S=pandas.DataFrame({"time":times,"beta1":0,"beta2":0,"beta3":0,"beta4":0,"beta5":0,"beta6":0,"beta7":0}) +store_I=pandas.DataFrame({"time":times,"beta1":0,"beta2":0,"beta3":0,"beta4":0,"beta5":0,"beta6":0,"beta7":0}) +store_R=pandas.DataFrame({"time":times,"beta1":0,"beta2":0,"beta3":0,"beta4":0,"beta5":0,"beta6":0,"beta7":0}) -solution=scipy.integrate.odeint(SIR_model, [S0, I0, R0], t, args=(.005, .5)) -solution= numpy.array(solution) +t=numpy.linspace(0,100) +for i in range(0,len(beta)): + params=(beta[i],gamma[i]) + sim=spint.odeint(func=SIR_model,y0=N,t=times,args=params) + store_S.iloc[:,i]=sim[:,0] + store_I.iloc[:,i]=sim[:,1] + store_R.iloc[:,i]=sim[:,2] plt.figure(figsize=[6,4]) @@ -107,7 +115,29 @@ def SIR_model(y, t, beta, gamma): plt.plot(t, solution[:, 2], label="R(t)") plt.show() +#max daily incidence +mdi=np.diff(store_I.iloc[:,1]) +Imax=numpy.max(mdi) +print("max daily incidence") +print(Imax) +#max daily prevelance +mdp=Imax/1000 +print("max daily prevalence") +print(mdp) +#percent affected over simulation +smallS=numpy.min(store_S.iloc[:,0]) +pas=(1/smallS) +print("percent affected over simulation") +print(pas) - +#basic reproduction number +S0=999 +I0=1 +R0=0 +beta=numpy.array([.0005, .005, .0001, .00005, .0001, .0002, .0001]) +gamma=numpy.array([.05, .5, .1, .1, .05, .05, .06]) +basicReproductionNumber=(((beta*(S0+I0+R0))/gamma)) +print("BasicReproductionNumber") +print(basicReproductionNumber) From a5a6e99e0bfcd00ea18d8ad115ac4f54f3a38ba7 Mon Sep 17 00:00:00 2001 From: Joe C Date: Thu, 9 Nov 2017 20:17:59 -0500 Subject: [PATCH 11/12] Completed 1 and 2 with both codes working. --- Ex10_Answers.py | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/Ex10_Answers.py b/Ex10_Answers.py index b8ff547..afa2b70 100755 --- a/Ex10_Answers.py +++ b/Ex10_Answers.py @@ -1,8 +1,15 @@ +##Ex 10 + +## Question 1 + +#import packages + import pandas import scipy import scipy.integrate as spint from plotnine import * + def ddSim (y,t0,r,K): #Define basic model for density-dependent growth N=y[0] dNdt=r*(1-N/K)*N @@ -29,6 +36,8 @@ def ddSim (y,t0,r,K): #Define basic model for density-dependent growth rates #Display the graph +####### part 2 of quesiton 1 + r=0.2 #Define a set growth rate and iitial population size N0=[1] times=range(0,50) @@ -47,32 +56,20 @@ def ddSim (y,t0,r,K): #Define basic model for density-dependent growth capacity=capacity+geom_line(store_Ks,aes(y="K3")) capacity #Display graph -## number 1 part 3 + +######## number 1 part 3 r=0.1 #Define a set growth rate and carrying capacity K=50 times=range(0,100) -times=range(0,500) - N0s=[1,50,100] #Define a list of initial sizes for each individual population store_N0s=pandas.DataFrame({"time":times,"N01":0,"N02":0,"N03":0}) #Initialize a data frame to store initial sizes with columns for each #initial size, and the times modeled for i in range(0,len(N0s)): #Set up a for loop that runs through each N0 in the list and models population growth using - params=(r,K) #other established parameters - sim=spint.odeint(func=ddSim,y0=N0s[i],t=times,args=params) - store_N0s.iloc[:,i]=sim[:,0] #Store population values in the column appropriate for the current carrying capacity - -popsize=ggplot(store_N0s,aes(x="time",y="N"))+theme_classic()#Set up base plot for population curves -popsize=popsize+geom_line(store_N0s,aes(y="N01")) -popsize=popsize+geom_line(store_N0s,aes(y="N02")) -popsize=popsize+geom_line(store_N0s,aes(y="N03")) - -popsize #Display graph - params=(r,K)#other established parameters N0=N0s[i] #Assign N0 a single value from the array of initial population sizes sim=spint.odeint(func=ddSim,y0=N0,t=times,args=params) @@ -158,3 +155,5 @@ def SIR_model(y, t, beta, gamma): basicReproductionNumber=(((beta*(S0+I0+R0))/gamma)) print("BasicReproductionNumber") print(basicReproductionNumber) + +print("!!! please see file named README for analysis completing question 2 !!!") From 3e25b5307f6a31f97e0a932e029022c2c993e2bf Mon Sep 17 00:00:00 2001 From: Joe C Date: Thu, 9 Nov 2017 20:49:31 -0500 Subject: [PATCH 12/12] updated the file to include a summary of observations --- README.md | 23 ++++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 897565a..6330f26 100644 --- a/README.md +++ b/README.md @@ -1 +1,22 @@ -# Intro_Biocom_ND_319_Tutorial10 \ No newline at end of file +# Intro_Biocom_ND_319_Tutorial10 + +# Answer for Question 2 + +## Please see Ex10_Answers.py for the code used to determine the following paragraph + + + +### The patterns I found while answering question 2 of the exercise are determined from observing the change in S, I, and R +### depending on which beta and gamma combination being observed --- for example to compare the effect of changing B I looked at +### numbers 3 and 4 in the chart given to us (starting with 1 not 0 in the list), and to determine the effect of gamma I +### looked at numbers 5 and 7 as gamma changed while beta stayed constant at 0.0001 --- + +### I determined that the greater the value of beta (assuming constant gamma) the faster S changed. +### Inversely, S changed faster when gamma was .05 compared to .06. +### therefore, in terms of change in S, beta and gamma are inversely related. + +### The change in I was faster when both beta and gamma are lower - therefore they are directly correlated (in a negative manner) + +### the change in R was greater when both beta and gamma are higher - therefore they are directly correlated in a positive manner. + +### Overall, the basic reproduction number Ro is nothing more than a ratio of beta/gamma. \ No newline at end of file