Skip to content

Instantly share code, notes, and snippets.

@ccritchfield

ccritchfield/_coronavirus_model_v2.readme

Last active July 5, 2020 20:48
Show Gist options
  • Save ccritchfield/5fe46dc5d99820434aa421692f4f4adf to your computer and use it in GitHub Desktop.
Save ccritchfield/5fe46dc5d99820434aa421692f4f4adf to your computer and use it in GitHub Desktop.
Download ZIP
Python - Coronavirus Model v2 - Dynamic People Lists
Raw
_coronavirus_model_v2.readme
v2 of the model adds and removes "attributes" from a person's list
as-needed to try to cut down on memory use.
EG: so a person won't be a list[0,0,0,0] tracking all 4 attributes
(status, incubation days, symptoms, transmission rate) all the time.
Instead, if they're uninfected, they'll simply be a 1 element list
with their status. As they get infected and incubate, then an
incubation attribute is added to count down. As they finish
incubating, the incubation attr will then get re-used to see
what symptom they have. As they become contagious, they get
a transmission attribute, and as they're no longer contagious,
that attribute is removed.
So, the idea is to spare RAM by sacrificing CPU cycles to
dynamically resize the person list.
you can read more about this version of the model here...
https://www.linkedin.com/pulse/covid-19-modeling-viral-funnel-pt-2-craig-critchfield-mis/
Raw
model_constants.py
##########################################################
# Viral Model - Global Constants
##########################################################
#
# Copyright Craig Critchfield 2020
#
##########################################################
"""
Shifted a lot of the constants over here,
b/c I had to write a lot of comments to explain
what they were doing, where I was getting them
from (internet links), and / or how they are
calculated. Gets annoying scrolling through
miles of code, so offloading them here to
reduce "main's" file size to make it easier
to organize and sift through.
-----------------------------
ABOUT THE MODEL
-----------------------------
Using Case Fatality Rate, we can pretty easily estimate
how likely someone will die from the virus if they get
infected.
The problem, though, is trying to determine how likely
someone will get infected. It's a constantly moving
number as more people become infected, thus increasing
chances of others getting infected if exposed.
It's also impacted by counter measures, like social
distancing, wearing masks, etc.
There's probably some grand algorithm that experts
use to determine all of this.
But, I just wanted to create a simple simulation
model to see how things go. The simulation will
loop through each person in a population to determine
if they come in contact with someone infected,
get infected, how long they incubate, chance to
infect others, etc.
It will do this day-by-day, and at the end of each
day it will rack up some KPI (key performance indicators)
to see how things are going and to chart out the
curves.
EG KPI...
* population mix of uninfected, infected, recovered & dead.
* generic % chance to get infected (this was my main goal with this)
* etc
I'm considering this a "black box" scenario, beacuse the
simulated people in it are just bouncing around w/o
knowing what they're doing or with any goal in mind.
We just want to see how they interact given some variables.
I pulled a lot of metrics from worldometers web site to
try to keep consistent in the data source.
Caveats to model....
1) assuming people can't get reinfected .. we're seeing a few
cases pop up where previously infected & recovered are
getting re-infected. But, it seems to be more that the
virus in them is flaring up again, not that they're getting
the virus again from a different source. Theoretically the
virus can mutate, and a new strain is formed that previuosly
infected people won't have antibodies to counter. But, for
this simple model, we'll just assume people get infected
(and recover or pass on) only once.
2) people won't be infectious forever .. we're using R0
(transmission rate) to determine how many people an
infected carrier can infect. This will be a % chance
to infect that many people, because we'll also use
counter measures % to reduce the infection chance.
So, while there are articles saying how long it takes
people to recover, that quarantining for 14 days is
prudent until you know if you have the virus or not,
all we're doing in this simulation is flagging someone
as getting the virus, they go through non-contagious
incubation period, then contagious incubation period,
and once they're contagious, each day after that they
will "try" to infect one other random person in population.
If no counter measures, then they will have 100% chance.
If counter measures, then they will have < 100% chance.
But, each "chance" adds up to the transmission rate.
Once they reach that transmission rate, they stop
"trying to infect" others. (This isn't coded to have
them act like they're purposefully trying to infect
others. They're just pebbles bumping around randomly
in a jar. They just will bump into someone every day,
and have a chance to infect them if infectious. We're
assuming no malicious intent to spread infection in
the model.)
"""
BORDER = "-" * 50 # border split to segregate terminal print feedback
#---------------------------------
# Transmission Rate (R0)
#---------------------------------
"""
Also known as Infection Rate, R0 or Reproduction rate.
R0 values are how likely something is to let another
thing through. In video game shaders, materials use
R0 values to determine how likely light is to transmit
through an opaque, translucent or transparent surface.
EG: rock would have an R0 of 0, since it blocks all light.
Clear Water would have an R0 of 1, because it lets all light
through.
So, R0 is a mathmatical term to determine transmission
of something through something else.
BBC has a nice article that explains R0 applied to
the virus perspective:
https://www.bbc.com/news/health-52473523
Basically, it's how many people, on average, a person
can infect if they have the virus. One person might
infect more. Another person might not infect any. But,
it averages out to a value where 1 person infected can
transmit it to R0 other people.
An R0 of 1 means 1:1, where 1 infected typically transmits
to another person. An R0 of 3 (which has been the suggested
value for coronavirus) means 1:3, 1 person transmits
it to 3 on average.
To determine how many people end up infected from a
"patient zero", you use the R0 in an exponential formula,
like so...
Infected = R0 Infection Rate ^ Power to Raise By
( ^ symbol in Python is the "Raise to Power of" operator )
Power to Raise by is how many degrees of separation
the virus has transmitted. EG: 1 person giving it to
3 others is 1 degree of separation, or raised to
power of 1. 1 person giving it to 3 others, then
those 3 giving it to 3 others (9 total) would
be 2 degrees of separation, or raising to power
of 2.
So, we can re-write the formula as such...
Infected = R0 Infection Rate ^ Degrees of Separation
In MS Excel you'd write it as...
Infected = power(R0, Degrees of Separation)
So, the higher the R0, the scarier it gets quickly
due to exponential increase.
Infected = 3^1 = 3
Infected = 3^2 = 9
Infected = 3^3 = 27
Infected = 3^4 = 81
Infected = 3^5 = 243
Infected = 3^6 = 729
Infected = 3^7 = 2,187
Infected = 3^8 = 6,561
Infected = 3^9 = 19,683
Infected = 3^10 = 59,049
So, 1 person infected could branch it out to 59,000
people if no counter measures (masks, distancing, etc)
are implemented.
CDC looked at studies from Wuhan China suggesting R0
of 2.2 to 2.7. Then they did their own modelling, and
it suggested an R0 of 5.7.
https://wwwnc.cdc.gov/eid/article/26/7/20-0282_article
I'm just going to go with R0 3 I've heard experts suggest.
"""
# corona virus
TRANSMISSION_RATE = 3
# flu
#TRANSMISSION_RATE = 1.35
#---------------------------------
# Counter Measures
#---------------------------------
"""
Counter measures means we're implementing proactive prevention
to try to slow the spread: wearing masks, social distancing,
isolation / quarantine, etc.
While we know counter measures work, because we can see the curve
flatten when implemented, we don't know for sure how effective
each one specifically is. IE: there's no definitive research
showing % effectiveness of wearing mask (for one person or both),
or social distancing.
There's a graphic going around saying that wearing masks means
the other person is 70% less likely to get the virus, and if
both people are wearing masks then 90% less likely. But,
Reuters did some fact-checking and found out that graphic isn't
valid. It's just some made-up numbers. I think someone felt if
they put numbers to the mask wearing, then people could latch
on to the idea of wearing masks more due to tangible / concrete
effectiveness. But, we simply don't have concrete / tangible metrics
for it yet.
https://www.reuters.com/article/uk-factcheck-coronavirus-mask-efficacy/partly-false-claim-wear-a-face-mask-covid-19-risk-reduced-by-up-to-98-5-idUSKCN2252T6
Likewise for social distancing, hand-washing, etc.
So, we're just going to distill this down into a single
"counter measures" metric we can use to skew the chance
someone can transmit the virus to others.
Counter measures will equal how effective they are, on a
scale of 0 to 100% effectiveness.
But, it acts as reduction to infection_chance. So, we invert it.
Infection Chance % * ( 100% - Counter Measures % )
So, as CM goes up, it lets less infection pass through.
Also, if we have a TRANSMISSION_RATE of 3 (1 person can
give it to 3 others), we'd need a TRANSMISSION_RATE of 1
to at least flatten the curve.
We can easily do the math to determine what kind of
counter measure % we'd need to just flatten the curve by
setting counter measures to X and solving for it...
TRANSMISSION_RATE * ( INFECTION_CHANCE - COUNTER_MEASURES ) = NEW_TRANSMISSION_RATE
3 * (100% - x%)= 1
100% - x% = 1/3 # divided both sides by 3
100% - x% = 33% # convert 1/3 to % rounded down
- x% = 33% - 100% # subtract 100% from both sides
- x% =-67% # reduce right-side equation
x% = 67% # divided both sides by -1
So, we'd need at least 67% counter measures to get TRANSMISSION_RATE
down to 1 or so. More if we want the new cases to lower via
R0 < 1
"""
"""
!!! DEPRECATED !!!
Counter measures are now lists in model_main.py
that let us control the days certain CM's start/
stop. This lets us simulate communities locking
down and opening up over time. Keeping this here
mostly for the notes above, so folks will understand
where I'm going with counter measures in the model.
"""
#COUNTER_MEASURES = 0.00 # no counter measures .. TRANSMISSION_RATE of 3.0
#COUNTER_MEASURES = 0.25 # 25% counter measures .. new TRANSMISSION_RATE of 2.25
#COUNTER_MEASURES = 0.35 # 35% counter measures .. new TRANSMISSION_RATE of 1.95
#COUNTER_MEASURES = 0.50 # 50% counter measures .. new TRANSMISSION_RATE of 1.5
#COUNTER_MEASURES = 0.66 # 66% counter measures .. new TRANSMISSION_RATE about 1.0
#COUNTER_MEASURES = 0.75 # 75% counter measures .. new TRANSMISSION_RATE about 0.75
#---------------------------------
# Vaccine Immunity
#---------------------------------
"""
A vaccine will stop the virus, but after it's
developed it takes time to manufacture,
distribute, and get masses of people injected.
We can simulate this by kicking off the day it
arrives, and then start incrementing immunity
as an extra counter measure to simulate more
and more people becoming immune over time.
"""
VACCINE_DAY = 0 # day vaccine is available (0 = no vaccine)
VACCINE_IMMUNITY_INCREMENT = 0.01 # counter measure amount vaccine adds every day after vaccine available
#---------------------------------
# Case Fatality Rate (CFR)
#---------------------------------
"""
Helps us determine liklihood of death if someone
gets the virus (and is confirmed to have it via
testing / diagnosis). As cases show up, they go
through the "pipe" until some kind of outcome
occurs, either death or recovery. So...
Outcomes = deaths + recovered
CFR = deaths / outcomes = deaths / (deaths + recovered)
EG:
worldometer US stats for May 18th, 2020...
https://www.worldometers.info/coronavirus/country/us/
Cases = 1,543,842
Deaths = 91,683
Recovered = 351,936
Outcomes = Deaths + Recovered = 91,683 + 351,936 = 443,619
CFR = Deaths / Outcomes = 91,683 / 443,619 = 0.2066705889513299
CFR = 0.21 ... 21%
I've monitored CFR for May, and for US it hangs around 20%.
So, for this simple model I'm going to just hard-code 20% as CFR
to determine death if someone gets the virus.
"""
CASE_FATALITY_RATE = 0.2
#---------------------------------
# Incubation Period
#---------------------------------
"""
https://www.acpjournals.org/doi/10.7326/M20-0504
Virus takes 3 to 14 days to incubate, 5 on-average.
When a person gets infected, we can use this function
to determine how many days they'll silently spread the
virus until symptoms show up.
Since it's 3 to 14, 5 on-average, it's a right-skewed
bell curve, meaning cases stack more towards the 3 to 5
range then the 6 to 14 range. IE: people seem to generally
show symptoms sooner rather then later.
To represent this, we'll use a pre-canned function
that lets us plug in a range (3 to 14) and average (5)
and it'll create a distribution curve from it and select
a day amount that's statistically representative of that
curve.
IE: if we took all population and gave them incubation
days, if we avg'ed them all together we'd get 5.
Python has a "normal" distribution function that we can
feed an average and standard deviation to and have it
spit out a statistically relevant value.
We know the average, it's 5 days. But, we don't know
the standard deviation. However, we can ad-hoc it by
taking the (max day - min days) / 4, which is the
standard deviation estimation formula in statistics.
"""
# INCUBATION_MIN = 3
# INCUBATION_AVG = 5
# INCUBATION_MAX = 14
# INCUBATION_STD = ( INCUBATION_MAX - INCUBATION_MIN ) / 4
# # debug ... check standard deviation
# #print("incubation std_dev = " + INCUBATION_STD)
# def IncubationPeriod():
# incubation = triangular( INCUBATION_MIN, INCUBATION_AVG, INCUBATION_MAX )
# # incubation = normal( INCUBATION_AVG, INCUBATION_STD )
# incubation = round( incubation ) # round incubation to nearest integer ( < 0.5 rounds down, 0.5+ rounds up )
# incubation = int( incubation ) # convert to integer and return value
# return incubation
# # debug ... test incubation function
# if 1:
# for i in range(1000):
# inc = IncubationPeriod()
# print(inc)
"""
!!!
normal() is giving numbers < 3 (and even negative incubation numbers).
triangular() isn't weighting the center point more towards 5 avg days.
so, until I can figure out how to create a more robust incubation curve
I'm just going to hard set incubation time to the avg 5 days
Also, note that I'm doing 5 days inclusive to incubate, which means
day 5 is also considered incubation. So, once we hit 6+, that's when
symptoms will show up.
!!!
"""
INCUBATION_DAYS = 5
#---------------------------------
# Contagious Period
#---------------------------------
"""
We need to know how long someone will be contagious, so we can track
how long they have to spread the virus.
https://www.health.harvard.edu/diseases-and-conditions/if-youve-been-exposed-to-the-coronavirus
'We know that a person with COVID-19 may be contagious 48 to 72 hours
before starting to experience symptoms. Emerging research suggests that
people may actually be most likely to spread the virus to others during
the 48 hours before they start to experience symptoms.'
'Most people with coronavirus who have symptoms will no longer be contagious
by 10 days after symptoms resolve. People who test positive for the virus
but never develop symptoms over the following 10 days after testing are
probably no longer contagious, but again there are documented exceptions.
So some experts are still recommending 14 days of isolation.'
So, what this tells us is that "it's complicated", b/c it depends on if
someone shows symptoms or not.
They're contagious 2 to 3 days before they show symptoms. With an avg
incubation time of 5 days, that means they will become contagious within
2 to 3 days after they get the virus (5 - 2 or 3 = 2 or 3). If we just
go with 2 days before symptoms, then we can get:
Incubation Days - Contagious Lead = 5 - 2 = 3 days to incubate w/o contagiousness
total days 1 2 3 4 5 6 7 8 9 10 11 12 13
incubation days 1 2 3 4 5 - - - - - - - -
contagious days - - - 1 2 3 - - - - - - -
The reason they only have 3 contagious days is because model will
have them try to infect only 1 person per-day, and after 3 days
they have hit their Transmission Rate (R0) of 3.
But, theoretically, in real life a person may be incubating and
then contagious for up to 14 days. this is why they suggest people
quarantine for 14 days, b/c if you don't have symptoms after 14
days (to err on the side of caution), then you will most likely
not be contagious after 14 days.
Symtomatic people are contagious during their recovery time + 10 days after
recovery. And, it depends on mild vs. severe symptoms...
'It depends on how sick you get. Those with mild cases appear to recover
within one to two weeks. With severe cases, recovery can take six weeks
or more. According to the most recent estimates, about 1% of infected
persons will succumb to the disease.'
mild = 2 wks * 7 days = 14 days w/symptoms + 10 days after symptoms = 24 days contagious
severe = 6 wks * 7 days = 42 days w/symptoms + 10 days after symptoms = 52 days contagious
so...
contagious days = incubating but contagious overlap days + recovery days + 10 days after recovery
If we go with 2 days incubating and contagious, then we get the following breakdown
asymptomatic = 2 incubating but contagious + 0 recovery + 10 after recovery = 12 days contagiuos
mild = 2 incubating but contagious + 14 recovery + 10 after recovery = 26 days contagious
severe = 2 incubating but contagious + 42 recovery + 10 after recovery = 54 days contagious
Since we're going to use the number as a countdown, we'll have to
kick it off after "incubation days - incubation but not contagious days"
In this case 5 - 2 = 3 days.. so 3 days a person will incubate while
not being contagious.. and then next day they will become contagious
and their contagious days will get added to them and start counting down.
During that time they will (obviously) be contagious, and we will
then run chances of them getting others infected.
"""
"""
!!!
Decided to simplify this down and skip the 'recovery days.'
Once a person hits R0, they move to outcome, which eliminates
them from processing in order to save cpu cycles.
!!!
"""
INCUBATION_AND_CONTAGIOUS_DAYS = 2
INCUBATION_NOT_CONTAGIOUS_DAYS = INCUBATION_DAYS - INCUBATION_AND_CONTAGIOUS_DAYS
RECOVERY_DAYS_ASYMPTOMATIC = 0
RECOVERY_DAYS_MILD = 14 # 2 wks * 7
RECOVERY_DAYS_SEVERE = 42 # 6 wks * 7
RECOVERY_DAYS_POST = 10 # 10 days after recovery needed to no longer be contagious
DAYS_CONTAGIOUS_ASYMPTOMATIC = INCUBATION_AND_CONTAGIOUS_DAYS + RECOVERY_DAYS_ASYMPTOMATIC + RECOVERY_DAYS_POST
DAYS_CONTAGIOUS_MILD = INCUBATION_AND_CONTAGIOUS_DAYS + RECOVERY_DAYS_MILD + RECOVERY_DAYS_POST
DAYS_CONTAGIOUS_SEVERE = INCUBATION_AND_CONTAGIOUS_DAYS + RECOVERY_DAYS_SEVERE + RECOVERY_DAYS_POST
# make a list out of it, so we can use
# person.symptoms value as index to look up value
DAYS_CONTAGIOUS = [ DAYS_CONTAGIOUS_ASYMPTOMATIC, DAYS_CONTAGIOUS_MILD, DAYS_CONTAGIOUS_SEVERE ]
# every day that goes by, a person's days_incubated goes down
# if it is <= INCUBATION_NOT_CONTAGIOUS_DAYS
#---------------------------------
# Person Attributes
#---------------------------------
"""
We use a person "object" to track people
in population bumping around. We use int
variables to track certain things, so
let's setup some constants to make it
easier to refer to them in english instead
of having to remember what their numbers are.
Originally created an actual person object,
but it ran very slow. So, then I switched to
tracking each person as a dict, but that was
also slower then just tracking them as lists.
Also, tried using enums to track person
attributes, but that ran 3x slower then just
one-off constant variables.
The irony of python is if you code properly,
eg: create objects/classes, use enums, etc like
you should, everything runs a lot slower. That's
why so many people find hacks and shortcuts to
improve performance. But, it can make code look
like a hot mess that's hard to understand.
Since I need brute force performance to process
millions of people in the model, I'm ignoring
"form over function" and instead going with
basic, brute force function over elegant form.
However, I'm still opting for code that's easy
to read. Hence I'm using if/else's and such
instead of trying to find hackish ways to
get a bit more speed by making the code more
obtuse.
STATUS
0 = uninfected
1 = infected (today, skip incubation today)
2 = incubating
3 = contagious (incubated to point of contagious, but might not be showing symptoms yet)
4 = recovered
5 = dead
INCUBATE
how many days they've incubated, up to INCUBATION_DAYS
SYMPTOMS
0 = asymotpmatic (or never infected if status = 0 uninfected)
1 = mild
2 = severe (needs hospitalization, possible death on outcome)
TRANSMIT
TRANSMISSION_RATE people person can transmit to.
Represents a percentage, wtih each infection attempt
dropping the transmit value by 100% if it's 1+,
or using the 0.x leftover decimal if less then 1.
EG: Coronavirus has transmission rate of 3. which
means a contagious person passes the virus on to
3 other people (on average). We can consider that
to be 3 people * 100% chance to infect. In our
model, we have a contagious person attempt to infect
only 1 person per-day. So, we just go 1 person * 100%
and then reduce the TRANSMIT down from 3 to 2.
If we were dealing with decimal transmission rates,
like the 1.3 for the flu, then a contagious
person could 100% try to infect 1 person, and then
30% try to infect another person. Sum 100% + 30%
and you get the 1.3 transmission rate.
This gives us a base chance to infect another person,
thought. It gets hindered by counter measures,
vaccine immunity, and the infected person's own
extra counter measures they might take via symptoms.
(EG: a severely infected, hospitalized person won't
be out and about much, so they have extra counter
measures that lower their risk of infecting others.)
"""
# person attribute constants to easily
# refer to the array columns in the code
STATUS = 0
INCUBATE = 1
TRANSMIT = 2
# person.status
UNINFECTED = 0 # has not been infected yet
INFECTED = 1 # infected today, skip incubating until tomorrow
INCUBATE_NOT_CONTAGIOUS = 2 # incubating, but not contagious yet
INCUBATE_CONTAGIOUS = 3 # incubating, and contagious (silent transmission)
SYMPTOMS_CONTAGIOUS = 4 # done incubating, contagious w/ symptoms (symptom counter measures kick in)
RECOVERED = 5 # no longer contagious, outcome = recovered
DEAD = 6 # no longer contagious, outcome = death
# we need to check if someone is incubating
# or contagious. So, create some tuples
# that will make it easier to "if x in y"
STATUS_INCUBATING = ( INCUBATE_NOT_CONTAGIOUS, INCUBATE_CONTAGIOUS )
STATUS_CONTAGIOUS = ( INCUBATE_CONTAGIOUS, SYMPTOMS_CONTAGIOUS )
# person.symptoms
# we're blowing out system memory with
# very large populations, so we're going to
# get rid of the dedicated "SYMPTOMS" attr
# on people, and just re-use "INCUBATE"
# attr. Initially, "INCUBATE" will track
# # of days incubating. After person has
# finished incubating, then attr will
# shift over to a value for their SYMPTOM.
# Incubation days will be positive integers,
# and symptoms will be negative. We can just
# do an abs() or ( x * -1) to convert it to
# positive during processing. We can
# determine if person is incubating or
# symptomatic via STATUS which reflects
# such, so INCUBATE is just a tracking
# sub-variable
ASYM = 0 # no symptoms
MILD = -1 # mild symptoms
SEVR = -2 # severe symptoms, needs hospitalization, CASE_FATALITY_RATE can cause death when determing outcome
#---------------------------------
# Symptom Ratios
#---------------------------------
"""
It's still unclear how many people that get the virus can become
asymptomatic, and whether they can spread the virus. Because,
there seems to be conflicting studies.
The issue is viral load in the human body is not an "on/off" switch.
A person may have the virus, but not show symptoms (personal immunity
and true asymptomatic). Or, they may show such mild symptoms that they
don't think they have it. Or, they think it's something else
(eg: seasonal allergies, common migraines they usually get, etc).
I stared at a few sources, but overall this article sums things up:
https://www.sciencealert.com/a-physician-answers-5-questions-about-asymptomatic-covid-19
From the variuos studies it references, we have a range of asymptomatic
in sample groups from 10% to 50% high.
And Icelandich study suggests 50% of people tested were asymptomatic:
https://www.healthline.com/health-news/50-percent-of-people-with-covid19-not-aware-have-virus
The Diamond Cruise study suggests 17% asymptomatic:
https://med.stanford.edu/content/dam/sm/id/documents/COVID/AsymptCOVID_TransmissionShip.pdf
CDC in this article suggests 1 in 4 people in US could be asymptomatic:
https://www.livescience.com/undetected-infections-coronavirus-widespread.html
The problem we're facing is that most people that get tested are ones
that present themselves with symptoms. Folks that don't have symptoms
are not likely to go get tested, b/c they think they're fine.
So, studies that focus on a proactive testing of a sample, even if some
in the sample are not showing symptoms, seem to be the best to go off
of. So, the Icelandic study of mass testing of a sample size is
a pretty decent one to take into account.
I would also take what the CDC says into account, b/c they're
aggregating all kinds of study and research into models.
So, I would put asymptomatic between the 25% the CDC says, and the
50% Icelandic study says.
The other caveat to this is that there could be factors in each
population that cause more or less transmission, infection, asymptomatics,
etc.
EG: The Icelandic study may later find out there's some segment of
population genetically predisposed to dealing with the virus easier.
This may attribute to why they had 50% asymptomatic. There could
be a larger concentration of those genetically predisposed asymptomatics
in their population. Meanwhile, in the US, we might have more folks
genetically predisposed to get symptoms.
The other factor is sample bias. If Iceland truly did a simple random
sample, then 50% would be a decent number to go with.. for their
population. But, their population is not the same mix as, say,
the United States. (Again, might be enough of a genetic difference
in population to create differences, hence why each country, state,
region, etc would need to do mass testing on their own).
The problem the CDC may be facing is they could be drawing conclusions
from their own biased samples, which would be all the reactive testing
going on in the states from people presenting themselves. This would
then make sense why the asymptomatic is lower in the US then Iceland,
b/c, again, if someone doesn't feel sick, they most likely won't go
get a test. So, we could be missing some 25% of asymptomatic people
in the US that haven't tested.
But, for now, we'll just go with a 25% of infected become asymptomatic
per CDC, since we're doing US model.
For the other 75% that get infected, this article suggests
80% of them end up with mild symptoms, which means 20% end
up with severe symptoms:
https://www.healthline.com/health-news/50-percent-of-people-with-covid19-not-aware-have-virus
Mild = 0.75 * 0.8 = 0.60 = 60% of total infected
Severe = 0.75 * 0.2 = 0.15 = 15% of total infected
So, we have the following breakdown:
--------------------------
Infected Symptoms (United States)
--------------------------
25% asymptomatic
60% mild
15% severe (need hospitalization)
--------------------------
The next issue is what is the chance an asymptomatic person can
transmit the virus. There seems to be conflicting sources.
This article by UC Davis suggests asymptomatic can transmit,
but in very limited fashion:
https://health.ucdavis.edu/coronavirus/resources/covid-19-faqs-for-health-professionals.html
This article has the CDC director saying the 1 in 4 asymptomatic
could be leading widespread transmission of the virus:
https://www.sciencealert.com/here-s-what-we-know-so-far-about-those-who-can-pass-corona-without-symptoms
Plenty of articles are saying that people are most contagious
during the 1 to 2 day incubation period before they start
showing symptoms. (in our model that would be the INCUBATION_AND_CONTAGIOUS_DAYS).
We're not going to get into a bunch of % chances to transmit,
and have them decrease over time.
Instead, we're going to go with the TRANSMISSION_RATE used
up evenly over time until it's gone and person has outcome.
TRANSMISSION_RATE is basically an average number of people
someone can infect.
our model handles this by right-skewing by front-loading
the infection chances earlier rather then later.
But, doing it this way for now means:
1) the model remains simple
2) it front-loads infection chance to others by having them try
to give it to 3 others over 3 days if possible which goes
along with the notion that a person is most infectious
during their pre-symptomatic incubation days
(IE: INCUBATION_AND_CONTAGIOUS_DAYS)
"""
ASYM_RATIO = 0.25
MILD_RATIO = 0.60
SEVR_RATIO = 0.15
# these are <= in model, EG: if random() <= ASYM_THRESHOLD
ASYM_THRESHOLD = ASYM_RATIO # range 0 to 25%
MILD_THRESHOLD = ASYM_RATIO + MILD_RATIO # range 25% to 85%
#SEVR_THRESHOLD = 1 - MILD_THRESHOLD # range 85%+ (not used, b/c it's the "else" part of the if statement)
#---------------------------------
# Symptom-specific Counter Measures
#---------------------------------
"""
when people get sick, they tend to convalesce,
more so as they get sick more. We don't
have any hard metrics for how much people with
symptoms convalesce, but we can make some
reasonable assumptions.
asymptomatic ... wouldn't know they are
ill (no symptoms), so they won't change their
behavior. So, they won't take additional counter
measures.
mild (doesn't need hospitalization) will
most likely stay home more, possibly trying to
self-quarantine just in case. So, they take more
counter measures.
severe is considered hospitalized in this model
(admitted, and either on oxygen or moved to a
ventilator). They will obviously be put under
heavier counter measures by the hospital to keep
others from getting sick. But, they are not
completely isolated; there is still a chance
someone could get infected from them. (we have
nurses and doctors getting infected from working
with covid patients). So, their counter measures
will be greatly increased, but still pose a
chance of infection.
So, we're going to have each symptom type
apply a multiplier to the infection chance to
possibly cut it down more to simulate more
counter measures taken by someone that's ill.
By using a multiplier, we can ensure that
we reduce infection chance, but it will never
go to 0% by just convalescing.
Also note that since a patient is infectious
after 3 days in this model, and then show
symptoms after 5 days, they will do silent
transmission for 2 days before showing symptoms
and having addt'l counter measures kick in.
So like this...
day 0 .. infected, but not incubating
day 1 .. incubating, not contagious
day 2 .. incubating, not contagious
day 3 .. incubating, not contagious
day 4 .. incubating, contagious, silent transmission (most dangerous)
day 5 .. incubating, contagious, silent transmission (most dangerous)
day 6 .. symptoms, contagious, addt'l counter measures (if any)
day 7 .. already hit transmission rate, apply outcome and skip going fwd
"""
ASYM_COUNTER_MEASURES_MULTIPLIER = 0.0 # asymp don't do anything different
MILD_COUNTER_MEASURES_MULTIPLIER = 0.5 # mild quarantine more and limit movement
SEVR_COUNTER_MEASURES_MULTIPLIER = 0.9 # severe are hospitalized under greater counter measures
# make a tuple out of these, so we can
# use the person.symptom int value as an
# index to look up the symptom counter
# measure value.
SYMPTOM_COUNTER_MEASURES = \
(
ASYM_COUNTER_MEASURES_MULTIPLIER,
MILD_COUNTER_MEASURES_MULTIPLIER,
SEVR_COUNTER_MEASURES_MULTIPLIER
)
# we're tracking counter measures as
# "more is better", but using them in
# "less blocks more infection", so we
# invert them when we use them in
# 1 - CM fashion. So, we can pre-calc
# that here to save some processing
ASYM_COUNTER_MEASURES_MULTIPLIER_INV = 1.0 - ASYM_COUNTER_MEASURES_MULTIPLIER
MILD_COUNTER_MEASURES_MULTIPLIER_INV = 1.0 - MILD_COUNTER_MEASURES_MULTIPLIER
SEVR_COUNTER_MEASURES_MULTIPLIER_INV = 1.0 - SEVR_COUNTER_MEASURES_MULTIPLIER
# tuple for inverted symptom cm's
SYMPTOM_COUNTER_MEASURES_INVERSE = \
(
ASYM_COUNTER_MEASURES_MULTIPLIER_INV,
MILD_COUNTER_MEASURES_MULTIPLIER_INV,
SEVR_COUNTER_MEASURES_MULTIPLIER_INV
)
#---------------------------------
# Community Simulation
#---------------------------------
"""
The problem with creating a massive population
is that we treat everyone like they can instantly
bump into anyone else in the model. This isn't
true in real life.
EG: if we're simulating US population, a person
in Florida that's contagious can't instantly
infect someone that's in California.
To work around this, we create many smaller
populations to run them at the same time, but
use a TRAVEL_CHANCE to have infectors choose
to infect locally vs. remotely.
The TRAVEL_CHANCE will act as another counter
measure, because it can keep an outbreak
contained in a population if there's a low
travel chance.
This would simulate cities or states deciding
to close borders and preven travel in order
to isolate and contain the outbreak.
The pop sizes we create depend on what we're
trying to simualte.
EG: Instead of doing a single population for
the US, we could do a population per-state,
and see how they spread.
I live in Dallas / Fort-Worth TX metroplex,
so instead of doing DFW as a single population,
we can break it out into the many suburbs and
small towns that make up the 6M+ people.
And, we don't have to be specificially accurate
here.
For instance, for DFW, we can just use an
assumed "community size" to divide the population
by.
EG: 6M divided by, say, 20 townships... that
would give us ~300K per township in the metroplex.
That would be good enough to get a ballpark
idea of how the virus might spread in the
metroplex if we implement various counter
measures and travel restrictions.
TRAVEL_CHANCE
How this will work is when an infector is
up for infecting someone, they will first
roll a dice to see if they fall within (<=)
TRAVEL_CHANCE. If they do, then they went
outside their community / town to remotely
infect. Otherwise, they stayed local.
So, TRAVEL_CHANCE of 0%, it means complete
lockdown & no travel.
TRAVEL_CHANCE is dependant on the populations
we're trying to simulate. For instance,
it's far more likely that people within
DFW will travel amongst townships, so if
we were modelling that we'd use a higher
TRAVEL_CHANCE. But, if we were modelling
the US States, we'd use a lower travel
chance, since many people are not travelling
outside their state on a regular basis.
If we really wanted to start getting fancy,
we'd apply a separate TRAVEL_CHANCE to each
population, and go even further to do
network analysis by assigning certain
populations connecting more then others.
EG: Certain cities live along highways
that act as key arteries / supply lines
in the US. So, they would be much more
likely to have people travel back-n-forth
then, say, some remote rural area hardly
anyone visits.
That is going beyond this simulation,
though. For now, we just want a basic
model we can tweak some easily-understood
values on.
"""
# single-digit travel chance for state-to-state
TRAVEL_CHANCE = 0.05 # 5% chance person was remote instead of local when infecting
# higher travel chance for city-to-city
#TRAVEL_CHANCE = 0.5
Raw
model_functions.py
##########################################################
# Viral Model - Functions / Methods
##########################################################
#
# Copyright Craig Critchfield 2020
#
##########################################################
import model_constants as mc
import model_populations as mp
from random import random
#------------------------------
# Infect
#------------------------------
"""
Helper functions to kick off
infection in people.
"""
# infect person, and rack up
# new infection counter
# newly infected will skip incubation today
def infect( person, dailystats ):
dailystats["new_infected"] += 1
person[mc.STATUS] = mc.INFECTED
# find random initial person
# and infect them to kick off first wave
def infectpatientzero( totalpops, populations ):
if totalpops > 1:
select = mp.randompopulation( totalpops )
population = populations[ select ]
else:
population = populations[ 0 ]
people = population["people"]
random_person = mp.randomperson( population["size"] )
patient_zero = people[ random_person ]
infect( patient_zero, population["dailystats"] )
# find someone uninfected and infect them
# to kick off new wave of virus
# this represents virus coming back
# into population from outside source
# and will then play out inside pop
# once again
def infectnextwave( totalpops, populations ):
if totalpops > 1:
select = mp.randompopulation( totalpops )
population = populations[ select ]
else:
population = populations[ 0 ]
people = population["people"]
size = population["size"]
while True:
select = mp.randomperson( size )
person = people[ select ]
if person[mc.STATUS] == mc.UNINFECTED:
infect( person, population["dailystats"] )
break
#------------------------------
# Determine Symptoms
#------------------------------
"""
When person past incubation period,
determine if Asymp, Mild or Severe
symptom breaks down as follows
======================
0-25% ... 25% asymptomatic
25-85%... 60% mild
85%+ ... 15% severe (need hospitalization)
======================
100% total
we can generate a random 0.0 to 1.0 number,
then if/else it to add a symptom based on
cumulative category values
"""
def symptoms( person, dailystats ):
category = random() # random number from 0.0 to 1.0
if category <= mc.ASYM_THRESHOLD: # range 0.00 to 0.25 asymp
dailystats["new_asym"] += 1
dailystats["asym_contagious"] += 1
person[mc.INCUBATE] = mc.ASYM
elif category <= mc.MILD_THRESHOLD: # range 0.25 asymp to (0.25 asymp + 0.60 mild) = range 0.25 to 0.85
dailystats["new_mild"] += 1
dailystats["mild_contagious"] += 1
person[mc.INCUBATE] = mc.MILD
else: # range (0.25 asymp + 0.60 mild) to 1.00 = range 0.85 to 1.00
dailystats["new_sevr"] += 1
dailystats["sevr_contagious"] += 1
person[mc.INCUBATE] = mc.SEVR
#------------------------------
# Incubation Cycle
#------------------------------
"""
if someone infected previous day, kick off incubation now.
otherwise, incubate people that have been incubating.
once they hit incubation days where they're contagious,
they'll start possibly infecting others (silent transmission).
when they hit max incubation days, they're assigned symptoms.
We do an initial incubation pass, so if someone up the chain
infects someone down the chain, we don't want that person
down the chain incubating the same day they got infected.
We want them to spend a day getting infected, then next
day+ they're incubating.
"""
# if incubating, determine when
# contagiuos and when symptoms hit
def incubate( population ):
dailystats = population["dailystats"]
people = population["people"]
for person in people:
# if person was infected yesterday
# flag them to start incubating today
# and add INCUBATE attr to person list
if person[mc.STATUS] == mc.INFECTED:
person[mc.STATUS] = mc.INCUBATE_NOT_CONTAGIOUS
person.append( 0 )
# if person is in an incubating status
# rack up incubation days
if person[mc.STATUS] in mc.STATUS_INCUBATING:
person[mc.INCUBATE] += 1
# incubated to point of being contagious
# flip status, rack up metrics, add TRANSMIT attribute
if person[mc.STATUS] == mc.INCUBATE_NOT_CONTAGIOUS \
and person[mc.INCUBATE] > mc.INCUBATION_NOT_CONTAGIOUS_DAYS:
person[mc.STATUS] = mc.INCUBATE_CONTAGIOUS
dailystats["incu_contagious"] += 1
person.append( mc.TRANSMISSION_RATE )
# incubated to point of showing symptoms
# flip status, rack up metrics, flip INCUBATE to track symptom
elif person[mc.STATUS] == mc.INCUBATE_CONTAGIOUS \
and person[mc.INCUBATE] > mc.INCUBATION_DAYS:
person[mc.STATUS] = mc.SYMPTOMS_CONTAGIOUS
dailystats["incu_contagious"] -= 1
symptoms( person, dailystats )
#------------------------------
# Determine Outcomes
#------------------------------
"""
Once person has possibly infected up to
TRANSMISSION_RATE of other people, we'll
jump right to an outcome for them to take
them out of processing.
In real life, mild symptoms would take
up to 2 weeks or so to recover while
severe cases would take up to 5 weeks.
Since we're modelling transmission based
on transmission rate, not how long someone
is recovering, then after they infect
up to TRANSMISSION_RATE, we can just
take them out of play.
People with severe symptoms (hospitalized)
have CASE_FAtALITY_RATE chance to die
upon outcome. Otherwise, they'll recover.
And, folks that are Asymptomatic or suffering
from mild symptoms are considered to automatically
recover as well.
"""
def outcomes( population ):
dailystats = population["dailystats"]
people = population["people"]
for person in people:
# if person in a contagious status
# and done transmitting to others
if person[mc.STATUS] in mc.STATUS_CONTAGIOUS \
and person[mc.TRANSMIT] == 0:
# done with TRANSMIT field, so delete it
# asymptomatic & mild auto-recover.
# severe has chance of death per CFR
del person[mc.TRANSMIT]
symptom = person[mc.INCUBATE]
recovered = True
if symptom == mc.ASYM:
dailystats["asym_contagious"] -= 1
elif symptom == mc.MILD:
dailystats["mild_contagious"] -= 1
elif symptom == mc.SEVR:
dailystats["sevr_contagious"] -= 1
if random() <= mc.CASE_FATALITY_RATE:
recovered = False
if recovered:
dailystats["new_recovered"] += 1 # rack up recovered count
person[mc.STATUS] = mc.RECOVERED
else:
dailystats["new_dead"] += 1 # rack up death toll
person[mc.STATUS] = mc.DEAD
#------------------------------
# Counter Measures & Vaccine Immunity
#------------------------------
"""
Counter Measures in model might change
as days go by, and Vaccine Immunity
can kick in on a certain day and need
to rack up over time. So, we need
to re-calc them on a daily basis.
"""
# if population is doing vaccine (VACCINE_DAY > 0),
# and we've hit day it kicks in
def vaccineimmunity( dailystats ):
# if mc.VACCINE_DAY > 0 \
# and dailystats["sum_days"] >= mc.VACCINE_DAY:
if dailystats["sum_days"] >= mc.VACCINE_DAY > 0:
# add more immunity, but cap at 100%
dailystats["vaccine_immunity"] += mc.VACCINE_IMMUNITY_INCREMENT
dailystats["vaccine_immunity"] = max( 1.0, dailystats["vaccine_immunity"] )
# want to have counter measures
# dynamically change over the days
# https://covid19.healthdata.org/united-states-of-america/texas
# worldometer tracks when counter measures started
# impacting TX, and when they drop off
# we can use that to determine some counter
# measure time lines via days that pass in sim
# https://www.worldometers.info/coronavirus/usa/texas/
# texas virus deaths to compare by
"""
Populations can have different counter
measures kick in on different days.
So, we need to loop through the CM
day ranges to see what CM's apply for
the current day to the population.
Since we're doing a loop using i,
eventually we'll reach a 'IndexError:
list index out of range' error by
going beyond the range of days we
have entries for. When we do, we'll
just error trap that and return the
last counter measure value from
the population as an on-going,
final CM value.
We track the CM to use by looking
at the current day the sim is on
and comparing to the beginning and
ending date values, with ending
day being exclusive...
begin day <= current day < end day
So, if we had the following...
cm beg = 0, 30, 60, 90
cm end = 30, 60, 90
cm val = 0.0, 0.1, 0.2, 0.1
The function would do the following...
if 0 <= current day < 30
return 0.0 CM
if 30 <= current day < 60
return 0.1 CM
if 60 <= current day < 90
return 0.2
if 90 <= current day < (IndexError)
go to error trap
return 0.1 as on-going final CM
"""
def countermeasures( population ):
dailystats = population["dailystats"]
daynow = dailystats["sum_days"]
daybeg = population["countermeasurebeg"]
dayend = population["countermeasureend"]
dayval = population["countermeasureval"]
lenbeg = len(daybeg)
try:
for i in range( lenbeg ):
if daybeg[i] <= daynow < dayend[i]:
dailystats["counter_measures"] = dayval[i]
except IndexError:
lenval = len(dayval)
dailystats["counter_measures"] = dayval[ lenval - 1 ]
"""
one-stop function that re-calc's
daily counter measures and vaccinne
immunity to create an overall CM
for the day for the population.
"""
def immunity( population ):
dailystats = population["dailystats"]
countermeasures( population )
vaccineimmunity( dailystats )
# counter measures are tracked in
# "more is better" % fashion, but
# are used to reduce infection %,
# so we need to 1-% to invert them
# into reduction multipliers. We
# can prep these values now to
# avoid doing it every time a
# contagious person calculates
# an infection chance.
inv_counter_measures = 1.0 - dailystats["counter_measures"]
inv_vaccine_immunity = 1.0 - dailystats["vaccine_immunity"]
# since infection chance, counter measures,
# vaccine immunity, etc are multiplying
# together, we can pre-mul the cm & vi
# into a single inverse metric now
# for later use.
dailystats["inverse_counter_measures"] = inv_counter_measures \
* inv_vaccine_immunity
# to save cpu cycles, we can pre-calculate
# the cm's each contagious group will use,
# because folks with symptoms might get
# extra counter measures
asym = mc.SYMPTOM_COUNTER_MEASURES_INVERSE[ mc.ASYM ]
mild = mc.SYMPTOM_COUNTER_MEASURES_INVERSE[ mc.MILD ]
sevr = mc.SYMPTOM_COUNTER_MEASURES_INVERSE[ mc.SEVR ]
# incubating contagious don't receive any addtl
# counter measures
incu = dailystats["inverse_counter_measures"]
asym *= dailystats["inverse_counter_measures"]
mild *= dailystats["inverse_counter_measures"]
sevr *= dailystats["inverse_counter_measures"]
# plug them into our daily metrics
dailystats["inverse_counter_measures_incu"] = incu
dailystats["inverse_counter_measures_asym"] = asym
dailystats["inverse_counter_measures_mild"] = mild
dailystats["inverse_counter_measures_sevr"] = sevr
#------------------------------
# Infection Chance
#------------------------------
"""
Calculate chance contagious person has of
infecting another. This takes into account
their "transmit" value, counter measures,
vaccine immunity, and symptomatic counter
measures.
Since vaccine immunity & counter measures
are tracked in a "more is better" % value,
but are used to decrease the infection
chance %, we have to invert them via
( 1 - % ) to get the % amount they
decrease something.
EG: 75% counter measures will let
( 1 - 75% ) = 25% of the virus through
"""
def infectionchancebase( person ):
# if infector can infect 1 whole person
# use 100% infection chance
# and subtract 100% from transmit
if person[mc.TRANSMIT] >= 1:
basechance = 1.0
person[mc.TRANSMIT] -= 1
# otherwise, use whatever % they have left
# and zero out their transmit chances
# (they'll stop being contagious and
# have an outcome next day)
else:
basechance = person[mc.TRANSMIT]
person[mc.TRANSMIT] = 0
return basechance
# apply global counter measures &
# vaccine immunity (using inverted CM).
# also, determine person's contagious status
# which helps us determine any extra
# symptom counter measures. Also see how
# much of the daily infection chance they
# contribute (infectchance * 1 person in pop ratio)
def infectionchance( person, dailystats, oneperson ):
if person[mc.STATUS] == mc.INCUBATE_CONTAGIOUS:
countermeasures = "inverse_counter_measures_incu"
onepersoninfect = "new_infect_chance_incu"
else:
symptom = person[mc.INCUBATE]
if symptom == mc.ASYM:
countermeasures = "inverse_counter_measures_asym"
onepersoninfect = "new_infect_chance_asym"
elif symptom == mc.MILD:
countermeasures = "inverse_counter_measures_mild"
onepersoninfect = "new_infect_chance_mild"
else: # if symptom == mc.SEVR:
countermeasures = "inverse_counter_measures_sevr"
onepersoninfect = "new_infect_chance_sevr"
# determine their infection chance
# and rack up their daily inf chance contribution
infectchance = infectionchancebase( person )
infectchance *= dailystats[ countermeasures ]
dailystats[onepersoninfect] += ( infectchance * oneperson )
return infectchance
#------------------------------
# Transmission Cycle
#------------------------------
"""
personA can transmit virus to another if...
* infected
* incubated to contagious (shedding virus cells)
* still have uninfected to transmit to
If able to transmit, find random personB, and
attempt to infect. It's an attempt, because
it might be blocked:
* counter measures lower infection chance, so chance to not transmit
* personB may already be infected, dead or recovered, so can't get infected again
"""
def transmit( population, population_index, totalpops, populations ):
travelchance = population["travelchance"]
dailystats = population["dailystats"]
people = population["people"]
person_index = 0
for personA in people:
# if person is contagious
# see if they infect someone
if personA[mc.STATUS] in mc.STATUS_CONTAGIOUS:
infectchance = infectionchance( personA, dailystats, population["oneperson"] )
# if random 0.0 to 1.0 within threshold
# infect someone, otherwise it's been averted
if random() <= infectchance:
# if doing 2+ pops
# and travelling
if totalpops > 1 \
and random() <= travelchance:
# infect remotely
infectotherpopulation( totalpops, population_index, populations )
else:
# infect locally
infectsamepopulation( population, person_index )
# increment index to next person
person_index += 1
# when infecting a different population
# exclude infector's population,
# but pick anyone in that other population
def infectotherpopulation( totalpops, population_index, populations ):
select = mp.randomexcludepopulation( totalpops, population_index )
otherpop = populations[ select ]
select = mp.randomperson( otherpop["size"] )
personB = otherpop["people"][ select ]
if personB[mc.STATUS] == mc.UNINFECTED:
infect( personB, otherpop["dailystats"] )
# when infecting same population
# exclude infector from random person selection
def infectsamepopulation( samepop, person_index ):
select = mp.randomexcludeperson( samepop["size"], person_index )
personB = samepop["people"][ select ]
# if they're uninfected, infect them
if personB[mc.STATUS] == mc.UNINFECTED:
infect( personB, samepop["dailystats"] )
Raw
model_logging.py
##########################################################
# Viral Model - Debug Logging & Summary Outputs
##########################################################
#
# Copyright Craig Critchfield 2020
#
##########################################################
import model_constants as mc
import model_utility as mu
#------------------------------
# Compile Stats
#------------------------------
"""
Take our daily stats, add the new
stuff to the sums, then copy it
all to the total stat tracker.
"""
def compilestats( population ):
dailystats = population["dailystats"]
# add pop name to the stats
dailystats["popname"] = population["name"]
# infected
dailystats["sum_infected"] += dailystats["new_infected"]
dailystats["sum_uninfected"] = population["size"] \
- dailystats["sum_infected"]
# outcomes
dailystats["sum_recovered"] += dailystats["new_recovered"]
dailystats["sum_dead"] += dailystats["new_dead"]
dailystats["sum_outcome"] = dailystats["sum_recovered"] \
+ dailystats["sum_dead"]
# symptoms
dailystats["sum_asym"] += dailystats["new_asym"]
dailystats["sum_mild"] += dailystats["new_mild"]
dailystats["sum_sevr"] += dailystats["new_sevr"]
# contagious
dailystats["sum_contagious"] = dailystats["incu_contagious"] \
+ dailystats["asym_contagious"] \
+ dailystats["mild_contagious"] \
+ dailystats["sevr_contagious"]
# daily infection chance
dailystats["new_infect_chance"] = dailystats["new_infect_chance_incu"] \
+ dailystats["new_infect_chance_asym"] \
+ dailystats["new_infect_chance_mild"] \
+ dailystats["new_infect_chance_sevr"]
# get max contagious & daily infection chance
# to easily track peak of outbreak
# mostly for summary output to terminal
dailystats["max_contagious"] = max( dailystats["max_contagious"], dailystats["sum_contagious"] )
dailystats["new_infect_chance_max"] = max( dailystats["new_infect_chance"], dailystats["new_infect_chance_max"] )
# reset certain values in
# dailystats for next day
# accumulation
def restartstats( population):
dailystats = population["dailystats"]
# reset new counts
dailystats["new_infected"] = 0
dailystats["new_recovered"] = 0
dailystats["new_dead"] = 0
dailystats["new_asym"] = 0
dailystats["new_mild"] = 0
dailystats["new_sevr"] = 0
dailystats["new_infect_chance_incu"] = 0
dailystats["new_infect_chance_asym"] = 0
dailystats["new_infect_chance_mild"] = 0
dailystats["new_infect_chance_sevr"] = 0
dailystats["new_infect_chance"]
#------------------------------
# Debug Terminal Stats Output
#------------------------------
# how long we want the
# "value = " part
#before showing the actual value
PADLEN = 20
PADVAL = " "
# create string "name = value" with single var value
# padded "name =" to look more uniform in output
def modelsummaryvalue( varname, varval ):
strvar = varname.ljust( PADLEN, PADVAL )
strval = str( varval )
msgval = strvar + " = " + strval
return msgval
# create string "name = value" with dict value
# padded "name =" to look more uniform in output
def modelsummarystats( varname, dailystats ):
strvar = varname.ljust( PADLEN, PADVAL )
strval = str( dailystats[varname] )
msgval = strvar + " = " + strval
return msgval
# create string of "name = value (%)"
# padded "name =" to look more uniform in output
def modelsummaryratio( varname, dailystats, percent ):
strvar = varname.ljust( PADLEN, PADVAL )
strval = str( dailystats[varname] )
strper = str( percent )
msgval = strvar + " = " + strval + " (" + strper + ")"
return msgval
# debug terminal output
def modelsummary( population, decimals ):
# from sys import getsizeof
# see how badly increasing pop
# sizes are blowing out my RAM.
# "We're gonna need a bigger boat!"
# mem_population = getsizeof( population )
# mem_totalstats = getsizeof( totalstats )
dailystats = population["dailystats"]
# total folks with symptoms assigned
sum_symptom = dailystats["sum_asym"] \
+ dailystats["sum_mild"] \
+ dailystats["sum_sevr"]
# symptom percents
percent_asym = mu.formatcalculatepercent( dailystats["sum_asym"], sum_symptom, decimals)
percent_mild = mu.formatcalculatepercent( dailystats["sum_mild"], sum_symptom, decimals)
percent_sevr = mu.formatcalculatepercent( dailystats["sum_sevr"], sum_symptom, decimals)
# severe cases still alive
sum_sevr_alive = dailystats["sum_sevr"] \
- dailystats["sum_dead"]
# severe case alive / dead ratios
percent_sevr_alive = mu.formatcalculatepercent( sum_sevr_alive, dailystats["sum_sevr"], decimals )
percent_sevr_dead = mu.formatcalculatepercent( dailystats["sum_dead"], dailystats["sum_sevr"], decimals )
# population percentages
percent_pop_infected = mu.formatcalculatepercent( dailystats["sum_infected"], population["size"], decimals)
percent_pop_uninfected = mu.formatcalculatepercent( dailystats["sum_uninfected"], population["size"], decimals)
percent_pop_contagious = mu.formatcalculatepercent( dailystats["max_contagious"], population["size"], decimals)
# immunity ratios
percent_vaccine_immunity = mu.formatpercent( dailystats["vaccine_immunity"], decimals )
percent_counter_measures = mu.formatpercent( dailystats["counter_measures"], decimals )
percent_daily_inf_max = mu.formatpercent( dailystats["new_infect_chance_max"], decimals )
# death ratios
percent_mortality_rate = mu.formatcalculatepercent( dailystats["sum_dead"], population["size"], decimals )
from datetime import datetime
msg = []
msg.append( mc.BORDER )
msg.append( str(datetime.now()) )
msg.append( modelsummaryvalue( "popname", population["name"] ) )
msg.append( modelsummarystats( "wave", dailystats ) )
msg.append( modelsummarystats( "sum_days", dailystats ) )
msg.append( mc.BORDER )
msg.append( modelsummaryvalue( "total_pop", population["size"] ) )
msg.append( modelsummaryratio( "sum_infected", dailystats, percent_pop_infected ) )
msg.append( modelsummaryratio( "sum_uninfected", dailystats, percent_pop_uninfected ) )
msg.append( modelsummaryratio( "max_contagious", dailystats, percent_pop_contagious ) )
msg.append( mc.BORDER )
msg.append( " Herd Immunity @ " + str(percent_pop_infected) + " population" )
msg.append( " Vaccine Immunity @ " + str(percent_vaccine_immunity) )
msg.append( " Counter Measures @ " + str(percent_counter_measures) )
msg.append( " Max Daily Inf Chance @ " + str(percent_daily_inf_max) )
msg.append( mc.BORDER )
msg.append( modelsummarystats( "sum_outcome", dailystats ) )
msg.append( modelsummarystats( "sum_recovered", dailystats ) )
msg.append( modelsummaryratio( "sum_dead", dailystats, percent_mortality_rate ) + " mortality rate" )
msg.append( mc.BORDER )
msg.append( modelsummaryvalue( "sum_symptom", sum_symptom ) )
msg.append( modelsummaryratio( "sum_asym", dailystats, percent_asym ) )
msg.append( modelsummaryratio( "sum_mild", dailystats, percent_mild ) )
msg.append( modelsummaryratio( "sum_sevr", dailystats, percent_sevr ) )
msg.append( mc.BORDER )
msg.append( "sum_sevr_alive = " + str(sum_sevr_alive) + " (" + percent_sevr_alive + ")" )
msg.append( "sum_sevr_dead = " + str(dailystats["sum_dead"]) + " (" + percent_sevr_dead + ")" )
# print( mc.BORDER )
# print("population RAM = " + mu.bytestring(mem_population) )
# print("totalstats RAM = " + mu.bytestring(mem_totalstats) )
msg.append( mc.BORDER )
return msg
Raw
model_main.py
##########################################################
# COVID-19 Infection Model
##########################################################
#
# Copyright Craig Critchfield 2020
#
##########################################################
"""
Main processing loop
"""
#------------------------------
# TO DO - Possible Improvements
#------------------------------
"""
* Resource .. to monitor RAM / CPU use
* use triangular() to give better range of
incubation days for each person
(would need to track incubate days
as a countdown instead of count-up
then, triggering contagious when it
hits < incubate contagious days)
* error trap memory error and revert to I/O
on disk as fallback instead of crashing
(disk I/O adds tons of time to sim run,
so not going to implement this.)
* travel chance per-population, and needs
to be two-way, IE: population won't let
folks leave or enter, so have to check
both ways. Or, take both pops' travel
chances and multiply together to get a
single chance for a person A to go
to another community and community B
to let them in.
* population interaction ratios to have
certain pops interact more then others
(probably way beyond this sim, b/c
that gets into advanced social
network stuff that other programs
are better at, and we'd need a
lot more research to see which
cities interact more then others.)
"""
#------------------------------
# Imports
#------------------------------
"""
chucked a lot of code into separate modules
to organize and tidy up, b/c this was getting
super-long.
"""
import model_constants as mc # import model constants
import model_functions as mf # import model functions
import model_logging as ml # import model logging functions
import model_utility as mu # import model utility functions
import model_populations as mp # import model population funcs and pre-cans
from process_timer import procTimer
#from os import getcwd
from memory_profiler import profile
#----------------------------------
# Constants / Variables
#----------------------------------
# get working dir, b/c we'll shove
# out logs and csv's to it
# nix ... kept using directory where code
# was called from. using a .bat file
# to run the code from pypy folder would
# cause it to dump things to pypy folder.
# don't want that, so just hard-coding the
# output folder.
#FILEPATH = getcwd() + "\\"
FILEPATH = "D:\\python\\projects\\viral model\\"
#----------------------------------
# Main
#----------------------------------
@profile
def main( dolog = True, docsv = True ):
percent_decimal_places = 1
days_without_contagious = 0
totalstats = []
populations = []
# first US case hit Jan 20, 2020
# it's easier to track the model in 30 day
# month-to-month chunks, so we're going to
# assume model starts Feb+. We'll create
# counter measures for each month for US.
# it's important to note that with Transmission
# Rate (R0), what matters is if it's > 1, = 1
# or < 1...
# > 1 = cases rising, contagious person infecting more folks then resolving
# = 1 = cases steady, contagious person infecting 1 other person then resolving
# < 1 = cases lowering, contagious person infecting < 1 other person.
# Coronavirus has a Transmission Rate of 3, and
# we consider that "100% chance to infect 3 people."
# to get it to = or < 1, we have to drop that 100% down
# per-person we use it on. We can easily figure out
# what we need to do to flatten the curve (R0 = 1)
# by dividing 100% by our Transmission Rate.
# 100% / 3 = 1 / 3 = 33%
# so, we need to multiply 100% by 33% to get 33%
# infection through... and that would be ...
# 3 people * 33% = 0.99 = 99% .. about 100%
# which would mean our R0 has flattened.
# since we're doing (1 - CM) to invert counter
# measure amount, we can do (1 - Infect Chance)
# to see what CM's we need to set.
# 1 - 33% = 66.6% = 67% counter measures to flatten curve
# CM's > 67% will lower R0 curve
# with that in mind, we can model out how the
# United States responded to the virus each month...
#----------------------------------
# days month response
#----------------------------------
# Jan .. late January, first cases
# 0 - 30 Feb .. few counter measures ("locked down China", no masks or distancing)
# 30 - 60 Mar .. some counter measures ("wait n see" stance, isolation, distancing, but no masks)
# 60 - 90 Apr .. lockdown (isolation, distancing & masks), curve flattening
# 90 - 120 May .. places opening up, people not wearing masks, curve lowering
# 120 - 150 Jun .. protests .. lots of congregation w/o masks
# 150+ Jul+ we're just going to go with Jun CM's going fwd
# we use June's CM going fwd, b/c a lot of states
# are now open and acting like the virus is done.
# (even though it was transmission of 1 at best in
# some places, lowering a bit on some places, and
# just a slower-rising curve in others.) Overall,
# the US showed lowering curve in March, but June
# is seeing rising cases, people congregating at
# protests and jobs, people not wearing masks,
# people protesting having to wear masks, etc.
# so, our model will just assume we're going fwd
# with low counter measures.. that might change
# when they see another spike hit August or so.
# Feb Mar Apr May Jun Jul Aug
countermeasurebeg = [ 0, 30, 60, 90, 120, 150, 180 ]
countermeasureend = [ 30, 60, 90, 120, 150, 180 ]
# no counter measures at all.. to model out worst-case scenario
# countermeasureval = [0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 ]
# counter measures based on how US responded each month
# since symptoms also provide counter measures, we're not
# setting exactly 67% to get flattened curve.. as more folks
# go mild or severe, we automatically start isolating certain
# parts of population
countermeasureval = [0.00, 0.25, 0.60, 0.63, 0.45, 0.45, 0.45 ]
# if we had better proactive counter measures early on
# countermeasureval = [0.65, 0.65, 0.65, 0.65, 0.65, 0.65, 0.65 ]
#-------------------------
# multi-population experiment
#-------------------------
# generate population from pre-canned
# population = mp.loadpopulation( mp.test100, countermeasurebeg, countermeasureend, countermeasureval )
# populations.append( population.copy() )
# population = mp.loadpopulation( mp.test1K, countermeasurebeg, countermeasureend, countermeasureval )
# populations.append( population.copy() )
# population = mp.loadpopulation( mp.test10K, countermeasurebeg, countermeasureend, countermeasureval )
# populations.append( population.copy() )
# population = mp.loadpopulation( mp.test1M, countermeasurebeg, countermeasureend, countermeasureval )
# populations.append( population.copy() )
#-------------------------
# dallas / fort-worth pop as whole
#-------------------------
# population = mp.loadpopulation( mp.tx_dfw, countermeasurebeg, countermeasureend, countermeasureval )
# populations.append( population.copy() )
#-------------------------
# dallas / fort-worth multi-pop
#-------------------------
# load up 20 sub-communities (0 to 19)
# for Dallas / Fort-Worth, TX.
# rename them to keep track of each
"""
for i in range(20):
population = mp.loadpopulation( mp.tx_dfw_sub, countermeasurebeg, countermeasureend, countermeasureval )
population["name"] = population["name"] + "_" + str(i)
population["travelchance"] = 0.05 # limited travel during pandemic
populations.append( population.copy() )
"""
#-------------------------
# US multi-pop
#-------------------------
# US population is too big for my comp's RAM.
# stalls out while running with Mem Error.
# so need to use HDD I/O to read/write temp files.
# already loaded the state pop sizes into dicts
# and then loaded those into a list, so we can
# just iterate through the list to easily load the pops.
# !!!!!!!!!!!!!!!!!!!
# started updating code to do this, then once
# I got the read/write functions done decided to
# experiment with just California population
# of 39M. Making people for it, and saving
# to csv then reading back in just one time
# took 2.5 min... and that was just the file
# read / write; didn't include time it took
# to generate the population.
# for the model run, it'd have to handle
# read/writes of pops for each state, and
# do so for every "for pop in pops", so that's
# a boat-load of read/writes, disk thrashing,
# and tons of time.
# My computer can't handle this massive of
# population analysis, so I'm going to keep
# it more local, like TX or DFW, for now.
# since I can run that in RAM.
# !!!!!!!!!!!!!!!!!!!
diskIO = False
for state in mp.pop_US:
population = mp.loadpopulation( state, countermeasurebeg, countermeasureend, countermeasureval )
if diskIO:
population["filepath"] = FILEPATH
population["filename"] = "temp_" + population["name"] + ".csv"
# dump out the people list to temp file
mu.outputcsv( population["people"], population["filename"], population["filepath"] )
# purge people list
population["people"] = []
# now we can append it to pops
population["travelchance"] = 0.05 # limited travel during pandemic
populations.append( population.copy() )
#-----------------------------
# Main
#-----------------------------
print(mc.BORDER)
procTimer("run sim")
totalpops = len(populations)
wavenow = 1 # which wave we're on
wavemax = 2 # how many waves to do
# infect patient zero.
# if 2+ populations, pick a person from one at random
# otherwise pick a person from our only population at random
mf.infectpatientzero( totalpops, populations )
#-----------------------------
# loop until out of contagious people
#-----------------------------
while True:
# model runs in 3 phase...
# phase 1: maintenance
for population in populations:
mf.immunity( population ) # re-calc counter measures for day
mf.outcomes( population ) # determine outcomes
mf.incubate( population ) # incubate folks a day
# phase 2: activity
# now that daily statuses are sorted
# we can run around and infect people
population_index = 0
for population in populations:
mf.transmit( population, population_index, totalpops, populations )
population_index += 1
# phase 3: stats
# compile the daily numbers
# into totalstats tracker for output
total_contagious = 0
for population in populations:
ml.compilestats( population )
dailystats = population["dailystats"]
totalstats.append( dailystats.copy() )
dailystats["sum_days"] += 1
total_contagious += dailystats["sum_contagious"]
#-------------------------------------
# rack up how many days we've gone
# without someone being contagious
if total_contagious > 0:
days_without_contagious = 0
else:
days_without_contagious += 1
# if we've gone INCUBATION_DAYS w/o
# more contagious, then virus has
# snuffed out ... end sim
if days_without_contagious > mc.INCUBATION_DAYS:
# if we haven't hit our wavemax yet
# try to kick off another wave
if wavenow < wavemax:
wavenow += 1
for population in populations:
population["dailystats"]["wave"] = wavenow
# infect next wave patient zero.
# if 2+ populations, pick a person from one at random
# otherwise pick a person from our only population at random
mf.infectnextwave( totalpops, populations )
days_without_contagious = 0
# otherwise, end sim
# and dump log file
else:
# output model summary to terminal
if dolog:
# dump summary to log file
import datetime
nowstamp = datetime.datetime.now()
nowstamp = nowstamp.strftime( "%Y-%m-%d_%H-%M-%S" )
filename = "model_logging." + nowstamp + ".txt"
pathname = FILEPATH + filename
with open( pathname, 'a', newline = '' ) as txtfile:
for population in populations:
msg = ml.modelsummary( population, percent_decimal_places )
# dump summary to debug terminal
for i in range( 0, len(msg) ):
print( msg[i] )
txtfile.write( msg[i] + "\n" )
# exit loop
# which ends sim
break
# if we're still looping,
# reset daily stat trackers
for population in populations:
ml.restartstats( population )
#-----------------------------
procTimer("run sim")
# output stats to CSV
if docsv:
filename = "model_output.csv"
mu.outputcsv( totalstats, filename, FILEPATH )
#---------------------------
# Run Main
#---------------------------
profile = False
if profile:
# profiling performance for bottlenecks
fileprof = "model_profile.prof"
filetext = "model_profile.txt"
mu.profilemain( FILEPATH, fileprof, filetext )
else:
# not profiling, just run main
main()
#---------------------------
# End
#---------------------------
print(mc.BORDER)
Raw
model_populations.py
##########################################################
# Viral Model - Preset Populations
##########################################################
#
# Copyright Craig Critchfield 2020
#
##########################################################
"""
This contains variables and functions
used to manage people and populations
for the model.
It also contains pre-canned population
data to make it easy to make copies
of populations into the model to work
with on-the-fly.
"""
#------------------------------
# Imports
#------------------------------
import model_constants as mc
from random import randrange
"""
NOTE - I originally used randint(start, end),
but profiling showed it added overhead to the
program as an intermediary to just call randrange.
So, I switched over to use randrange( start, end ).
But, when trying to pick between 2 populations in
a populations list of 2 (randrange( 0, 1)), it
would only return 0 and get stuck in a loop if
the 1st pop was the one to call it.
With some experimentation, I realized that
randint is end-inclusive while randrange is
end-exclusive.
EG: randint(0, 2) will spit out 0, 1, but
randrange(0, 2) will only spit out 0 & 1.
I guess that's because typically if you're
calling randrange it's because you're giving
it the len() of the range you're dealing with,
so an object with 100 elements would have a
len of 100. But, a range of 0 to 99. So, if
you wanted to pick a random value from the
range you'd want 0 to 99... and returning
100 would give you an error since you don't
have item 100 in your list.
Anyways, I retooled the code to use randrange.
It's better this way, b/c I don't have to
pre-calc an "end" for randint of "len()-1"
to keep from blowing past the max elements
in population & people lists.
But, this is something to note in case I
decide to change the code later.
"""
#------------------------------
# Variables
#------------------------------
"""
Daily Stats is a template for
tracking the day's metrics for
a population. When we create
a population, we make a copy
of it, so each population can
track their own metrics. At
the end of each day, it will
get copied to the pop's totalstats
list to create a "spreadsheet"
of what's happened.
"""
dailystats = {
"popname": "",
"wave": 1,
"sum_days": 0,
"sum_uninfected": 0, # total pop - sum_infected
"sum_infected": 0,
"new_infected": 0,
"sum_recovered": 0,
"new_recovered": 0,
"sum_dead": 0,
"new_dead": 0,
"sum_asym": 0,
"new_asym": 0,
"sum_mild": 0,
"new_mild": 0,
"sum_sevr": 0,
"new_sevr": 0,
"sum_outcome": 0, # sum_dead + sum_recovered
"counter_measures": 0, # "more is better" value
"vaccine_immunity": 0, # "more is better" value
"inverse_counter_measures": 0, # 1 - counter_measures, so we can reduce infection chance
"inverse_counter_measures_incu": 0, # pre-calc'ed for incu contagious
"inverse_counter_measures_asym": 0, # pre-calc'ed for asym symptoms
"inverse_counter_measures_mild": 0, # pre-calc'ed for mild symptoms
"inverse_counter_measures_sevr": 0, # pre-calc'ed for sevr symptoms
"inverse_vaccine_immunity": 0, # 1 - vaccine_immunity, so we can reduce infection chance
"new_infect_chance_incu": 0,
"new_infect_chance_asym": 0,
"new_infect_chance_mild": 0,
"new_infect_chance_sevr": 0,
"new_infect_chance": 0, # sum of incu + asym + mild + sevr infect chances
"new_infect_chance_max": 0, # max chance at peak of outbreak
"incu_contagious": 0,
"asym_contagious": 0,
"mild_contagious": 0,
"sevr_contagious": 0,
"sum_contagious": 0, # incu + asym + mild + sevr contagious
"max_contagious": 0 # peak contagious people during outbreak
}
#------------------------------
# Load People into Population
#------------------------------
"""
load uninfected people into population.
give them TRANSMISSION_RATE to individually
work down through as they possibly
infect others.
"""
# make a person list object out of person attributes.
# since our attr vals are 0 to 255, we can run
# it as a bytearray instead of a list which speeds
# up processing some
# !!!
# don't use bytearray... can't do floats,
# and if we use model for something with a
# transmission rate that isn't an integer
# eg: flu with R0 of 1.3 to 1.4, then bytearray
# blows up with error. So, just keep it a normal
# list
# !!!
# !!!
# just kick off uninfected folks with STATUS
# attribute.. we'll append more attributes to
# person as they get infected, but we want to
# min / max memory usage.
# !!!
def makeperson():
status = mc.UNINFECTED
person = [ status ]
return person
# make a person, and use them to
# create the people list in population
def loadpeople( size ):
person = makeperson()
people = []
for i in range( size ):
people.append( person.copy() )
return people
#------------------------------
# Generate Population
#------------------------------
"""
We store precanned population names
and sizes, because we know what those
values are now. But, we load them into
a population dict when we run the model,
that way we can adjust transmission rate,
counter measures for the pop, etc.
"oneperson"
we rack up each person's infection chance
(chance when possibly infecting someone)
to see what the daily infection chance is.
(IE: to see "what the weather's like").
to do that, we have to multiply the chance
they calculate when possibly infecting by
1 / POP_TOTAL, then we add it to their
contagious group's infection sum to see
how much infection chance each group is
giving, which we can then sum to a total
infection chance. So, we need to pre-calc
the % ratio of 1 person in the population.
"travelchance"
how likely someone in a population is to
go to a remote population to infect instead
of infecting locally. A population will
use "travelchance" to head out, but a
population will "defend" itself from
remote visitors via counter measures.
It stands to reason that a population
taking a pandemic seriuosly will both
limit travel from it's occupants and
increase counter measures to protect
them. But, we can model out odd pops
that decide to have high counter measures
to protect themselves, but also have
high travel chance letting them head
out and potentially spread. This would
be a place where they have a hard time
spreading in their local population,
but people are easily spreading it to
other populations (if those populations
don't have good counter measures
themselves).
"""
def loadpopulation( poptoload, countermeasurebeg, countermeasureend, countermeasureval ):
name = poptoload["name"]
size = poptoload["size"]
population = {
"name": name,
"size": size,
"filename": "", # if doing HDD I/O, this will be temp file location
"filepath": "",
"oneperson": 1.0 / size, # % of 1 person in pop for daily infection chance sums
"travelchance": 0, # chance folks in pop will go to other pops to infect
"countermeasure": 0, # day's amount infect chance is reduced
"countermeasurebeg": countermeasurebeg.copy(),
"countermeasureend": countermeasureend.copy(),
"countermeasureval": countermeasureval.copy(),
"dailystats": dailystats.copy(),
"people": loadpeople( size )
}
return population
#------------------------------
# Pick Random
#------------------------------
"""
need to pick simple random samples
of various things. randint calls
randrange, but creates overhead,
so got rid of it and just went
with randrange. Lists use zero-
based counting,
"""
# pick a random population to infect first
# if dealing with 2 populations, randrange(0, 1)
# will only return 0, not 0 or 1.
def randompopulation( totalpops ):
return randrange( 0, totalpops )
# pick random person to infect first
def randomperson( popsize ):
return randrange( 0, popsize )
"""
When someone is infecting others, we
need to exclude infector when picking
random, and exclude their population
when they travel to infect another pop.
Python doesn't have a "random excluding"
pre-canned func, so we make our own.
"""
def randomexclude( end, exclusion ):
while True:
randomvalue = randrange( 0, end )
if randomvalue != exclusion:
break
return randomvalue
def randomexcludeperson( popsize, this_person ):
return randomexclude( popsize, this_person )
def randomexcludepopulation( totalpops, this_pop ):
return randomexclude( totalpops, this_pop )
#------------------------------
# Pre-Canned Populations
#------------------------------
"""
US population 2020
https://www.worldometers.info/world-population/us-population/
Experts think herd immunity starts kicking in when 70% of the
population is infected.
https://www.jhsph.edu/covid-19/articles/achieving-herd-immunity-with-covid19.html
Our model just keeps running until it
hits INCUBATION_DAYS w/o having any more
contagious people. This lets us see how
much of the pop creates herd immunity
"naturally" in the simulation instead
of forcing a pre-conveived end-limit on
it. Because the whole point of this model
is to let all the marbles bump around and
see how they interact instead of just
using the preconveived end-result assumptions.
Setting up pre-canned populations to
import and work with.
"""
test100 = {
"name": "test100",
"size": 100
}
test1K = {
"name": "test1K",
"size": 1000
}
test10K = {
"name": "test10K",
"size": 10000
}
test100K = {
"name": "test100K",
"size": 100000
}
test1M = {
"name": "test1M",
"size": 1000000
}
test10M = {
"name": "test10M",
"size": 10000000
}
# Dallas / Fort-Worth, TX population
tx_dfw = {
"name": "tx_dfw",
"size": 6301000
}
# some places are better modelled as a lot of
# sub-pops or communities, EG: DFW is like 20+
# smaller towns all converged into a metroplex.
# to model it out better, we can take the
# population, and divide it by a general
# sub-community count, then iterate and copy
# that sub-community pre-canned pop to generate
# sub-pops totalling the main pop.
# EG: for DFW we can take 6,301,000 metro pop
# divided by 20 "major" communities to get
# 315,050‬ in 20 sub-pops.. we can load those
# 20 sub-pops and then use travel restrictions
# to see how individual communities restricting
# travel might help curtail viral spread.
# Dallas / Fort-Worth, TX generic sub-community
tx_dfw_sub = {
"name": "tx_dfw_sub",
"size": 315050
}
#---------------------------
# US state populations
#---------------------------
"""
Want to model the US, but it's
better to do each state as a population
to let them interact a bit organically
instead of just using US pop as a whole.
!!! trying to do US pop blows out my 16gb of ram !!!
"""
pop_AK = {"name": "pop_AK", "size": 734002}
pop_AL = {"name": "pop_AL", "size": 4908621}
pop_AR = {"name": "pop_AR", "size": 3038999}
pop_AZ = {"name": "pop_AZ", "size": 7378494}
pop_CA = {"name": "pop_CA", "size": 39937489}
pop_CO = {"name": "pop_CO", "size": 5845526}
pop_CT = {"name": "pop_CT", "size": 3563077}
pop_DC = {"name": "pop_DC", "size": 720687}
pop_DE = {"name": "pop_DE", "size": 982895}
pop_FL = {"name": "pop_FL", "size": 21992985}
pop_GA = {"name": "pop_GA", "size": 10736059}
pop_HI = {"name": "pop_HI", "size": 1412687}
pop_IA = {"name": "pop_IA", "size": 3179849}
pop_ID = {"name": "pop_ID", "size": 1826156}
pop_IL = {"name": "pop_IL", "size": 12659682}
pop_IN = {"name": "pop_IN", "size": 6745354}
pop_KS = {"name": "pop_KS", "size": 2910357}
pop_KY = {"name": "pop_KY", "size": 4499692}
pop_LA = {"name": "pop_LA", "size": 4645184}
pop_MA = {"name": "pop_MA", "size": 6976597}
pop_MD = {"name": "pop_MD", "size": 6083116}
pop_ME = {"name": "pop_ME", "size": 1345790}
pop_MI = {"name": "pop_MI", "size": 10045029}
pop_MN = {"name": "pop_MN", "size": 5700671}
pop_MO = {"name": "pop_MO", "size": 6169270}
pop_MS = {"name": "pop_MS", "size": 2989260}
pop_MT = {"name": "pop_MT", "size": 1086759}
pop_NC = {"name": "pop_NC", "size": 10611862}
pop_ND = {"name": "pop_ND", "size": 761723}
pop_NE = {"name": "pop_NE", "size": 1952570}
pop_NH = {"name": "pop_NH", "size": 8936574}
pop_NH = {"name": "pop_NH", "size": 1371246}
pop_NM = {"name": "pop_NM", "size": 2096640}
pop_NV = {"name": "pop_NV", "size": 3139658}
pop_NY = {"name": "pop_NY", "size": 19440469}
pop_OH = {"name": "pop_OH", "size": 11747694}
pop_OK = {"name": "pop_OK", "size": 3954821}
pop_OR = {"name": "pop_OR", "size": 4301089}
pop_PA = {"name": "pop_PA", "size": 12820878}
pop_PR = {"name": "pop_PR", "size": 3032165}
pop_RI = {"name": "pop_RI", "size": 1056161}
pop_SC = {"name": "pop_SC", "size": 5210095}
pop_SD = {"name": "pop_SD", "size": 903027}
pop_TN = {"name": "pop_TN", "size": 6897576}
pop_TX = {"name": "pop_TX", "size": 29472295}
pop_UT = {"name": "pop_UT", "size": 3282115}
pop_VA = {"name": "pop_VA", "size": 8626207}
pop_VT = {"name": "pop_VT", "size": 628061}
pop_WA = {"name": "pop_WA", "size": 7797095}
pop_WI = {"name": "pop_WI", "size": 5851754}
pop_WV = {"name": "pop_WV", "size": 1778070}
pop_WY = {"name": "pop_WY", "size": 567025}
# load all US state pops into a list
# to make it easier to iterate through
# and build population objects from
pop_US = [
pop_AK,
pop_AL,
pop_AR,
pop_AZ,
pop_CA,
pop_CO,
pop_CT,
pop_DC,
pop_DE,
pop_FL,
pop_GA,
pop_HI,
pop_IA,
pop_ID,
pop_IL,
pop_IN,
pop_KS,
pop_KY,
pop_LA,
pop_MA,
pop_MD,
pop_ME,
pop_MI,
pop_MN,
pop_MO,
pop_MS,
pop_MT,
pop_NC,
pop_ND,
pop_NE,
pop_NH,
pop_NH,
pop_NM,
pop_NV,
pop_NY,
pop_OH,
pop_OK,
pop_OR,
pop_PA,
pop_PR,
pop_RI,
pop_SC,
pop_SD,
pop_TN,
pop_TX,
pop_UT,
pop_VA,
pop_VT,
pop_WA,
pop_WI,
pop_WV,
pop_WY
]
Raw
model_utility.py
##########################################################
# Viral Model - General Utility Functions
##########################################################
#
# Copyright Craig Critchfield 2020
#
##########################################################
"""
These are util functions I created
to help with certain things, but
they're more general purpose then
model-specific. So, chucking them
into a util module to help trim
down the size of the model_functions
module to make it easier to scroll
and work on things there.
"""
#------------------------------
# Calculate & Format Percents
#------------------------------
# catch "divide by 0" errors and return 0
def calculatepercent( numerator, denominator ):
try:
return numerator / denominator
except ZeroDivisionError:
return 0
# format % as string to decimal places
def formatpercent( percentage, decimals ):
percent = "{:." + str(decimals) + "%}"
return percent.format( percentage )
# both above as one function call
def formatcalculatepercent( numerator, denominator, decimals ):
result = calculatepercent( numerator, denominator )
return formatpercent( result, decimals )
#------------------------------
# Format Bytes to String
#------------------------------
"""
sys.getsizeof() returns bytes. I want to
reduce huge byte sizes down to more
reasonable numbers with proper size
suffix on end.
EG: 100000 bytes / 1000 = 100 KB
To get ballpark sizes, we can just
divide by 1000 until we get a num len < 4
"""
"""
!!!
nixed this, b/c some folks feel getsizeof()
isn't a very accurate profile of resource.
So, think I'll try implementing resource
module instead later on.
!!!
"""
"""
def bytestring( bytesize ):
suffix = ( "Bytes", "KB", "MB", "GB", "TB", "PB", "EB", "ZB", "YB" )
result = bytesize
append = 0
# if we have even 1000,
# bump it to the next
# unit size as 1.0
while result > 999:
result /= 1000
append += 1
if 0: # ... debug
print("result = " + str(result))
print("append = " + str(append))
# round the result to 1 decimal place
result = round( result, 1 )
return str( result ) + " " + suffix[ append ]
"""
#------------------------------
# Output CSV
#------------------------------
"""
output list of lists or
list of dicts to csv file.
started with stat tracker
being a list of lists, then
switched to being a list of
dicts. Decided to keep the
code to do either or in
the future.
"""
# pass a list of lists of list of dicts
# will determine which by testing 1st
# object in list to see which csv func
# to call below for parsing data to csv
def outputcsv( output, filename, filepath ):
pathname = filepath + filename
if isinstance( output[0], list ):
listcsv( output, pathname )
elif isinstance( output[0], dict ):
dictcsv( output, pathname )
# in csv outputs below...
# 'w' ... opens file to overwrite
# newline = '' ... keeps from inserting blank rows between items
# list of lists to csv
def listcsv( output, pathname ):
from csv import writer#, QUOTE_ALL
with open( pathname, 'w', newline='' ) as newfile:
wr = writer( newfile )#, quoting = QUOTE_ALL )
wr.writerows( output )
# list of dicts to csv
def dictcsv( output, pathname ):
from csv import DictWriter
keys = output[0].keys()
with open( pathname, 'w', newline='' ) as newfile:
dictwriter = DictWriter( newfile, keys )
dictwriter.writeheader()
dictwriter.writerows( output )
def importcsvlist( filename, filepath ):
from csv import reader
pathname = filepath + filename
intlist = []
with open( pathname, 'r' ) as csvfile:
r = reader ( csvfile, quotechar = ' ' )
rowlist = list( r )
for row in rowlist:
rowlist = []
for value in row:
rowlist.append( int(value) )
intlist.append( rowlist )
return intlist
"""
!!!!!!!!!!!!!!!!
pickle serializes IO faster then just
dumping out csv's. So, might expand on this
later to do large file IO of populations
to from disk to keep from blowing out RAM.
Might find another way to do this. Main
issue is speed at this point. What would
be the fastest way to do disk IO w/o
thrashing HDD in the process?
!!!!!!!!!!!!!!!!
filepath = "D:\\python\\projects\\viral model\\"
filename = "model_test_pickle"
pathname = filepath + filename
# wb = write binary
pickle.dump( pop, open(pathname, "wb") )
# rb = read binary
pop = pickle.load( open(pathname, "rb") )
"""
#------------------------------
# Profile Code for Bottlenecks
#------------------------------
"""
cProfile dumps out a binary file
that has to get loaded into pstats
to then dump out as a txt file.
Really convoluted to have to write
so many lines of code to just dump
out profile stats to txt file.
This just profiles time constraints
in code (what's using up time), not
RAM or CPU or other resources.
"""
def profilemain( filepath, fileprof, filetext ):
# lazy import
from cProfile import run
from pstats import Stats
# generate path + name locations
outputprof = filepath + fileprof
outputtext = filepath + filetext
# profile function call
run( "main()", outputprof )
# open binary file,
# sort by cumulative times
# then dump to txt file
with open( outputtext, "w") as f:
ps = Stats( outputprof, stream = f)
ps.sort_stats('cumulative')
ps.print_stats()
Raw
process_timer.py
########################################
"""
Process Timer
As you do more complex programming, you might
find yourself in need of a process timer...
something that lets you time your overall
code run, or snippets of it.
This is where an ad-hoc process timer comes in.
---------------------------------------
Main Process Timer functions are:
procTimer(name: string) -> (returns nothing, but prints/logs results if ending a timer)
Send it a timer name as a string.
If it can't find the timer name in the timer dictionary, it creates
the timer and snapshots the current time it added it.
If it can find the timer name, then it removes the timer from
the dictionary and returns the elapsed time (stop time - start time)
procLog(activate: bool, filename: string) -> (returns nothing)
Use this function to activate text file logging of timers.
The text file dumps out to the same folder that the
"process_timer.py" code file is in (eg: if it's in
c:\temp\ then your log fille will end up in c:\temp\ too.
activate = boolean True/False ...
True = activate and log timers to log file
False = stop logging timers
filename = filename to use for appending timer string messages to.
default = "timer_log.txt", so if you dont' provide a value,
it will use "timer_log.txt".
procInfo() -> (returns nothing, but prints summary of proc timer stuff)
Call it to get a summary of whether logging is on/off,
what the log file name currently is, and a list of timers
still active in the process timer dictionary.
This is mainly useful if you're running code from the
console one line at a time, b/c the code and variables
there remain running after each line of input. So,
you can run procInfo() to get a quick check of timers
still active.
---------------------------------------
If you want to import process_timer.py stuff from any
code you're doing in the IDE without having to worry
about copy/paste'ing process_timer.py file to the
same folder as your code file, then you can toss
process_timer.py into ...
C:\\Users\\(username)\\Anaconda3\\Lib
(have to use double \\, otherwise Python thinks
the single slash is an escape character, even in comment )
That's the Windows directory that houses all of the
extra boiler-plate Python code libraries, eg: random.py.
(Not sure what directory to use for Linux or Mac).
The problem with doing this is that if you turn on
timer logging, it will spit out a log file to that
lib directory still. (I could have imported os.py
and used that to come up with directory creation
and such, but logging is an extracurricular feature
and I didn't want to open up the os.py can of worms.
This code file is already hitting scope-creep levels
of bloat as it is.)
I may expand on this later if I think of more
stuff to add to it, or optimize it or whatever.
But, for now, I feel it's good enough.
- Craig Critchfield
"""
########################################
# CODE LIBRARY IMPORTS
########################################
# datetime.datetime.now() function
# tracks time in milliseconds
from datetime import datetime
########################################
# GLOBAL VARIABLES
########################################
proc_timers = {} # dictionary stores timer names & start times
log_timers = False # flag if we're logging timers to txt log or not (default = no)
log_file = "timer_log.txt" # default text file name to log to
########################################
# FUNCTIONS - TIMER
########################################
# single function to call to start/stop timers
def procTimer( name ):
if name in proc_timers:
_endTimer(name)
else:
_addTimer(name)
# add new timer to dict
# with snapshot of current time
def _addTimer( name ):
proc_timers[name] = datetime.now()
# pop timer from dict,
# printing results
def _endTimer( name ):
stop = datetime.now()
start = proc_timers.pop(name)
procTime = stop - start
msg = '{:^20}'.format(name) # center name within 20 char padding
logtime = datetime.now().strftime("%Y-%m-%d %H:%M:%S")
msg = logtime + " ... " + "timer ... " + msg + " ... " + str(procTime) + " (h:mm:ss.ms)"
print(msg)
if log_timers:
logTimer( msg )
########################################
# FUNCTIONS - LOG FILE
########################################
# allows user to activate / deactivate log file
# and change the filename that's used
# Note that we have these vars setup as globals
# so we use the "global" keyword to tell the
# function to refer to the global vars, and not
# simply kick off local variables with the same names.
def procLog( activate = log_timers, filename = log_file):
global log_timers
global log_file
log_timers = activate
log_file = filename
# append timer to output log file
def logTimer( msg ):
file = open(log_file, "a") # create / open file to append
file.write(msg + "\n") # add line break after input
file.close # close file to free object
########################################
# FUNCTIONS - TIMER INFO
########################################
def procInfo():
border = "-" * 50
print(border)
print( "logging timers = " + str(log_timers))
print( "log file name = " + log_file)
print(border)
print( "active timers...")
for name in proc_timers:
print("\t* " + name)
print(border)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment