###################################################################### # This is my solution to HW1, which is now your template for HW2. # # You will require one or more functions from the random module; here, # we’ll just import all of them. from random import * # You’ll also need functions to check if files exist and are readable. from os.path import isfile from os import access, R_OK ###################################################################### # Specification: flip(p) flips a (weighted) coin. Returns True if a

######################################################################

 

# This is my solution to HW1, which is now your template for HW2.

 

#

 

# You will require one or more functions from the random module; here,

 

# we’ll just import all of them.

 

from random import *

 

# You’ll also need functions to check if files exist and are readable.

 

from os.path import isfile

 

from os import access, R_OK

 

######################################################################

 

# Specification: flip(p) flips a (weighted) coin. Returns True if a

 

# randomly generated float between 0 and 1 (inclusive) is less than or

 

# equal to the specified probability p.

 

#

 

def flip(p=0.5):

 

return(random() <= p)

 

######################################################################

 

# Specification: newPop(N,I,vp,de,di) creates a new population

 

# agents. A population is an N-element tuple of agents, where each

 

# agent is represented by a dictionary with two entries. The ‘state’

 

# of the agent (-1 <= s <= di+de) represents that agent’s current

 

# disease condition, while the agent’s ‘vaccine’ status (a Boolean)

 

# represents whether the agent has been vaccinated or not.

 

#

 

def newPop(N, I, vp, de, di):

 

pop=tuple([ {‘state’:-1, ‘vaccine’:flip(vp)} for i in range(N) ])

 

for i in sample(range(N), I):

 

# Will be preemptively reduced on first call to update().

 

pop[i][‘state’]=di+de+1

 

return(pop)

 

######################################################################

 

# Specification: update(pop, rp) is called at the beginning of each

 

# simulation day to update each agent’s infection status. Agents with

 

# infection state 1 (i.e., at the end of their infectious period) are

 

# set to 0 (recovered state) with probability rp, and otherwise set to

 

# -1 (susceptible state). All other infectious agents have their

 

# infection state reduced by 1. Returns the number of active

 

# infections remaining.

 

#

 

def update(pop, rp):

 

infections = 0

 

for i in range(len(pop)):

 

if pop[i][‘state’] == 1:

 

# Newly susceptible with p=(1-rp), else recovered. To

 

# avoid 0’s Boolean interpretation, we’re using 1-rp so

 

# that we can return -1 and default to 0. Might be more

 

# easily interpretable using an explicit if/else here.

 

pop[i][‘state’] = (flip(1-rp) and -1) or 0

 

elif pop[i][‘state’] > 0:

 

# Decrement state if still infectious

 

pop[i][‘state’] = pop[i][‘state’]-1

 

infections = infections + 1

 

return(infections)

 

######################################################################

 

# Specification: sim(N=100, I=1, m=5, vp=0, tp=(0.02, 0.01), de=3,

 

# di=5, rp=0.5, max=100, verbose=False) returns a list of integers,

 

# where each integer represents the number of infections on the

 

# corresponding simulation day (i.e., the pandemic curve). There are a

 

# large number of parameters, all with default values:

 

# N: number of agents in the simulation

 

# I: number of infections during simulation; initially,

 

# the number of infected agents on day 1.

 

# m: mixing parameter is the max number of daily interactions initiated

 

# by a given agent: randint(0,m) generates a given day’s tally.

 

# vp: probability that an agent has been vaccinated

 

# tp: tuple of transmission probabilities; first value corresponds

 

# to the shedding rate while exposed, the second to the shedding

 

# rate while infected (generally higher).

 

# de: number of days agent remains in exposed state

 

# di: number of days agent remains in infected state

 

# rp: probability of recovery (as opposed to becoming susceptible again)

 

# max: failsafe simulation time limit, in days

 

# verbose: Boolean value toggles user feedback

 

#

 

# A few more words about how we use our simple infection model (see

 

# Carrat et al, 2008). People are exposed for de days (the exposed

 

# period), then infected for di additional days (the infected period):

 

# they are infectious (that is, they can infect others) while they are

 

# either exposed or infected. We’ll use a “daily countdown” to model

 

# each agent’s disease state.

 

#

 

# Agents generally start the simulation as susceptible, that is, a

 

# status of -1. If they are infected on day t, we assume their exposed

 

# period runs for de days starting the next day, t+1, and then their

 

# infected period starts on day t+de and runs for di days. At the

 

# beginning of each day, update() updates each agent’s infection

 

# status.

 

#

 

# Carrat et al. (2008) indicate E~2 and I~7 for influenza, so let’s

 

# see a concrete example how this works. At infection on dy t, we set

 

# the agent’s state to di+de+1 (the extra 1 ensures the agent enters

 

# the infectious range on the day following infectoin, after the

 

# beginning-of-day status update). After the infection has run its

 

# course, update() flips a coin to decide if the agent is now immune

 

# (recovered, 0) or becomes suscptible once again (-1).

 

#

 

# If I=7, E=2:

 

# SUS REC I I+E I+E+1

 

# | | | | |

 

# -1 0 1 2 3 4 5 6 7 8 9 10

 

# . . |===================||====| . <= days

 

#

 

def sim(N=100, I=1, m=5, vp=0, tp=(0.01, 0.02), de=3, di=5, rp=0.5, max=100, verbose=False):

 

# Specification: susceptible(i) returns True if and only if agent

 

# i is susceptible (unvaccinated and has state counter set to -1).

 

#

 

def susceptible(i):

 

return(not pop[i][‘vaccine’] and pop[i][‘state’] == -1)

 

# Specification: exposed(i) returns True if and only if agent i is

 

# exposed (has state counter between di and di+de).

 

#

 

def exposed(i):

 

return(di < pop[i][‘state’] <= di+de)

 

# Specification: infected(i) returns True if and only if agent i

 

# is exposed (has state counter between 0 and di).

 

#

 

def infected(i):

 

return(0 < pop[i][‘state’] <= di)

 

# Specification: infectious(i) returns True if and only if agent i

 

# is infectious, that is, either infected or exposed.

 

#

 

def infectious(i):

 

return(0 < pop[i][‘state’] <= di+de)

 

# Specification: recovered(i) returns True if and only if agent i

 

# is in the recovered state.

 

#

 

def recovered(i):

 

return(pop[i][‘state’] == 0)

 

# Create agents.

 

pop = newPop(N, I, vp, de, di)

 

# Good to go, now run the simulation. We’ll use totinf to keep

 

# track of infection events.

 

totinf = I

 

# Similarly, we’ll use curve, a list, to keep track of how many

 

# infecteds there are each day, starting from day 0.

 

curve = [totinf]

 

# Run the simulation (i.e., keep looping) while there are

 

# infecteds remaining (that is, while curve[-1] > 0) or until you

 

# hit the failsafe number of rounds limit, max.

 

rounds = 0

 

while rounds < max:

 

# Beginning-of-day agent status update: returns the number of

 

# active infections on the new day, which we append to the

 

# pandemic curve.

 

curve.append(update(pop, rp))

 

# If verbose is True, show user feedback.

 

if verbose:

 

print(“Day {}: {} of {} agents infected.”.format(len(curve)-1, curve[-1], N))

 

# Quit the simulation if there are no infections left.

 

if curve[-1] == 0:

 

print(“Pandemic extinguished: {} days, {} infecteds, attack rate is {:3.1f}%.”.format(len(curve), totinf, (100*totinf)/N))

 

break

 

# For each infected agent, sample the population to find who

 

# they mix with, and then flip a coin to determine if a new

 

# infection occurs. For each new infection, update the

 

# appropriate agent’s status as well as totinf, the total

 

# infection count.

 

for i in range(N):

 

# Only worry about infectious agents

 

if infectious(i):

 

# Generate a list encounters by random sampling

 

# without replacement (see Python documentation).

 

for j in sample(range(N), randint(0,m)):

 

# Is there a new infection?

 

if susceptible(j) and flip(tp[int(infected(i))]):

 

# Update totinf, the total infection count.

 

totinf = totinf + 1

 

pop[j][‘state’] = di+de+1

 

# If verbose is True, show user feedback.

 

if verbose:

 

print(” Agent {} infected by agent {}.”.format(j, i))

 

# Update number of simulation rounds

 

rounds = rounds + 1

 

else:

 

print(“Pandemic persists: {} days, {} infecteds, attack rate is {:3.1f}%.”.format(len(curve), totinf, (100*totinf)/N))

 

# Return simulation curve.

 

return(curve)