exercise_10 Bruzzese_Inskeep_Nemera

exercise_10.py

@@ -0,0 +1,78 @@
+## exercise 10 ##
+
+# 1
+#Load packages
+import pandas
+import scipy
+import scipy.integrate as si
+from plotnine import *
+
+
+#Define custom function
+def popGrowth(y,t0,r,K,N):
+    N=y[0]
+    dNdt=r*(1-N/K)*N
+    return [dNdt]
+times=range(0,100)
+NO=[.01]
+#Set a pool of values for growth rate
+growthRates=[-.1,.1,.4,.8,1.0]
+#Dataframe for storing model output
+store_growthRates=pandas.DataFrame({"time":times,"r1":0,"r2":0,"r3":0,"r4":0,"r5":0})
+
+#Using a for loop to make my life easier
+for i in range(0,len(growthRates)):
+    pars=(growthRates[i],100,10)
+    sim=si.odeint(func=popGrowth,y0=NO,t=times,args=pars)
+    store_growthRates.iloc[:,i]=sim[:,0]
+print(store_growthRates)    
+
+#Plot 1- pop size as function of time
+g=(ggplot(data=store_growthRates)
+    +geom_line(store_growthRates,aes(x="time",y="r1"))+theme_classic()
+    +ylab("population size")
+    +geom_line(store_growthRates,aes(x="time",y="r2"),color="green")
+    +geom_line(store_growthRates,aes(x="time",y="r3"),color="red")
+    +geom_line(store_growthRates,aes(x="time",y="r4"),color="orange")
+    +geom_line(store_growthRates,aes(x="time",y="r5"),color="purple"))
+print(g)
+
+#set a pool of values for carrying capacity
+carryCap=[10,50,100]
+#Dataframe for storing model output
+store_carryCap=pandas.DataFrame({"time":times,"K1":0,"K2":0,"K3":0})
+NO = 1
+#for loop for Plot 2
+for i in range(0,len(carryCap)):
+    pars=(.2,carryCap[i],1)
+    sim=si.odeint(func=popGrowth,y0=NO,t=times,args=pars)
+    store_carryCap.iloc[:,i]=sim[:,0]
+print(store_carryCap) 
+
+#Plot 2- 
+k=(ggplot(data=store_carryCap)
+    +geom_line(store_carryCap,aes(x="time",y="K1"),color="blue")+theme_classic()
+    +ylab("population size")
+    +geom_line(store_carryCap,aes(x="time",y="K2"),color="green")
+    +geom_line(store_carryCap,aes(x="time",y="K3"),color="purple"))
+print(k)
+
+#Set a pool of values for init pop size
+initPop=[1,50,100]
+#Dataframe for storing model output
+store_initPop=pandas.DataFrame({"time":times,"N1":0,"N2":0,"N3":0})
+
+#Using a for loop to make my life easier
+for i in range(0,len(initPop)):
+    pars=(.1,50,initPop[i])
+    sim=si.odeint(func=popGrowth,y0=initPop[i],t=times,args=pars)
+    store_initPop.iloc[:,i]=sim[:,0]
+print(store_initPop)    
+
+#Plot 3- effect of initial Pop size differences
+p=(ggplot(data=store_initPop)
+    +geom_line(store_initPop,aes(x="time",y="N1"),color="orange")+theme_classic()
+    +ylab("population size")
+    +geom_line(store_initPop,aes(x="time",y="N2"),color="yellow")
+    +geom_line(store_initPop,aes(x="time",y="N3"),color="red"))
+print(p)

exercise_10_B.py

@@ -0,0 +1,141 @@
+#Load packages
+import pandas
+import scipy
+import scipy.integrate as si
+from plotnine import *
+
+# function
+def SIR (y,t0,beta,gamma):
+    S = y[0]
+    I = y[1]
+    R = y[2]
+    dS = -1*(beta*I*S)
+    dI = (beta*I*S)-(gamma*I)
+    dR = (gamma*I)
+    return dS, dI, dR
+
+# initial conditions
+times = range(0,500)
+NO = [999, 1, 0]
+
+# dataframe of gamma and beta values
+data = [{'beta' : .0005, 'gamma' : .05},
+        {'beta': .005, 'gamma': .5},
+        {'beta': .0001, 'gamma': .1},
+        {'beta': .00005, 'gamma': .1},
+        {'beta': .0001, 'gamma': .05},
+        {'beta': .0002, 'gamma': .05},
+        {'beta': .0001, 'gamma': .06},
+        {'beta': .0001, 'gamma': .9},
+        {'beta': .0001, 'gamma': .5},
+        {'beta': .0001, 'gamma': .25},
+        {'beta': .0001, 'gamma': .1},
+        {'beta': .0001, 'gamma': .05},
+        {'beta': .0001, 'gamma': .01},
+        {'beta': .0001, 'gamma': .001},
+        {'beta': .0001, 'gamma': .0001},
+        {'beta': .9, 'gamma': .0001},
+        {'beta': .5, 'gamma': .0001},
+        {'beta': .25, 'gamma': .0001},
+        {'beta': .1, 'gamma': .0001},
+        {'beta': .05, 'gamma': .0001},
+        {'beta': .01, 'gamma': .0001},
+        {'beta': .001, 'gamma': .0001},
+        {'beta': .0001, 'gamma': .0001},
+        {'beta': .00001, 'gamma': .0001},
+        ]
+my_data = pandas.DataFrame(data)
+
+# make lists to hold the results
+mdi = []
+mdp = []
+pa = []
+ro = []
+b = []
+g = []
+
+# start big  loop here
+for line in range(0,len(my_data),):
+    q = my_data.iloc[line]['beta']
+    p = my_data.iloc[line]['gamma']
+    params = (q, p)  # make tuple
+
+    b.append(params[0])  # append list
+    g.append(params[1])
+
+    # make dataframe
+    infection = pandas.DataFrame({"time":times,"S":0,"I":0,"R":0})
+
+    # sim shite
+    sim = si.odeint(func=SIR, y0=NO, t=times, args=params)
+
+    # fill dataframe
+    infection.iloc[:,2]=sim[:,0]
+    infection.iloc[:,0]=sim[:,1]
+    infection.iloc[:,1]=sim[:,2]
+
+    # calc max daily incidence
+    daily_incidence = []
+    for i in range(0,len(infection),):
+        if infection.time[i]==0:
+            continue
+        else:
+            I = infection.iloc[i]['I']
+            Iold = infection.iloc[i-1]['I']
+            incidence = I-Iold
+            daily_incidence.append(incidence)
+    max_daily_incidence = max(daily_incidence)
+    mdi.append(max_daily_incidence)
+
+    # calc max daily prevalence
+    daily_prev = []
+    for i in range(0,len(infection),):
+        I = infection.iloc[i]['I']
+        R = infection.iloc[i]['R']
+        S = infection.iloc[i]['S']
+        prev = I/(S+I+R)
+        daily_prev.append(prev)
+    max_daily_prev = max(daily_prev)
+    mdp.append(max_daily_prev)
+
+    #calc percent affected over simulation- use last time step (499)
+    I= infection.iloc[499]['I']
+    R= infection.iloc[499]['R']
+    S= infection.iloc[499]['S']
+    percent_affected = (I+R)/(S+I+R)
+    pa.append(percent_affected)
+
+    # basic reproduction number initial SIR
+    beta = params[0]
+    gamma = params[1]
+    I= infection.iloc[0]['I']
+    R= infection.iloc[0]['R']
+    S= infection.iloc[0]['S']
+    repo_number = (beta*(S+I+R))/gamma
+    ro.append(repo_number)
+
+# make a dataframe for results from all the lists
+results = pandas.DataFrame(
+    {'beta' : b,
+     'gamma' : g,
+     'max_daily_incide' : mdi,
+     'max_daily_prev' : mdp,
+     'percent_affect' : pa,
+     'repo_num' : ro})
+
+print results
+#need these to fill into a list or a dataframe
+# need to put all this intoa bigger loop
+
+'''
+* observations *
+We noticed that if beta is held constant while gamma gets smaller the max incidence, daily prevalence, percent infection, 
+and reproductive number all rise.  Thus a bigger gamma causes, things to get smaller If we hold gamma constant, as  beta
+gets bigger the max incidence, daily prevalence, percent affected, and  reproductive number all rise. Thus a small beta, 
+causes things to get smaller To summarize, a high beta and a small gamma cause the disease to have a higher rate and 
+srpead of infection in the population.
+
+'''
+
+
+