# Epidemic simulation over a large population

This is a follow-up to my question about epidemic simulation. The accepted answer works, but also mentions that I should look at modifying the external_function_call. Here's the complete working code, but it still takes quite a long time to run - in fact, it's slower than what I originally posted.

Here are the first few lines from the profiler, assuming a population size 100 (it's actually 625, which makes the code extremely slow):

     11185402 function calls (11184367 primitive calls) in 25.228 seconds

Ordered by: standard name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
1    0.108    0.108   25.229   25.229 20Dec_for_codereview.py:1(<module>)
1    0.000    0.000    0.530    0.530 20Dec_for_codereview.py:10(create_pop)
6585    0.019    0.000    0.019    0.000 20Dec_for_codereview.py:26(checkEqual)
748   10.670    0.014   16.705    0.022 20Dec_for_codereview.py:29(powerlaw_epidemic)
58    0.002    0.000    0.013    0.000 20Dec_for_codereview.py:66(inf_per_count_time)
6321    0.005    0.000    0.005    0.000 20Dec_for_codereview.py:70(<genexpr>)
1    0.008    0.008   16.731   16.731 20Dec_for_codereview.py:74(write_powerlaw_epidemics)


I'm very sure I'm not using lists as I should, which may be a part of the problem - checkEqual is checking whether all elements in the list are equal to 0.

And here's the code. The purpose of changing the list element to 0 was to keep the spot in the list, but exclude that element from the list of susceptible individuals at time t. I'm also not convinced I should define dist1 as I do, in the loop, but I can't see another way to do it. This code is far too slow, even with the improvements suggested in the epidemic simulation question. How can I make it faster?

import numpy as np
import scipy
from scipy import spatial
import json
import itertools
import collections

MAX_FAILED_ATTEMPS = 10

def create_pop(meanval, sigmaval, pop):
x_pos, y_pos = np.random.multivariate_normal(meanval, sigmaval, pop).T
dist_mat = np.asarray(zip(x_pos, y_pos))
pdistance = scipy.spatial.distance.pdist(dist_mat)
full_mat = scipy.spatial.distance.squareform(pdistance)

#Save the population - this happens, but commented out to not save files on someone else's machine
#with open('20Dec_population.txt', 'w') as f:
#    xyarray = np.array([x_pos, y_pos])
#    xyarray = xyarray.T
#    np.savetxt(f, xyarray, fmt=['%f', '%f'])

return full_mat

#http://stackoverflow.com/questions/3844801/check-if-all-elements-in-a-list-are-identical
def checkEqual(lst):
return lst[1:] == lst[:-1]

def powerlaw_epidemic(pop, susc, trans, inf_period, eps, full_mat):
#Purpose:
#Create a list of susceptible individuals based on the population
#From the susceptible individuals, choose one to be infectious at time 1
#Then the susceptible at that position in the list becomes 0
#And the infectious dictionart gets updated to be {susceptible[chosen_position]:1}
#While there are still susceptible individuals, allow the epidemic to continue

susceptible = list(range(1, pop + 1))
infectious = {}
t = 1
while not checkEqual(susceptible): #while there are still susceptible individuals (while susceptible is not all 0s)
if t == 1:
inf_1 = int(np.random.randint(1, pop, size=1))
infectious[inf_1] = t
inf_times = t
susceptible[inf_1-1] = 0
else:
for ind in susceptible:
if ind != 0:
inf_times = [i for i in infectious.values() if ((t - inf_period <= i and (i < t)))]
t_times = set(range(1, t)) #from 1 to t-1
infset = set(inf_times)
if not t_times.intersection(infset):
return infectious
l1 = [infectious.values().index(i) for i in inf_times]
inf_indivs = [infectious.keys()[i] for i in l1]
dist1 = [(full_mat[ind-1, j-1])**(-trans) for j in inf_indivs]
probinf = 1 - np.exp(-susc*sum(dist1) + eps)
if probinf >= np.random.uniform(0, 1):
infectious[ind] = t
susceptible[ind-1] = 0
else:
continue
t += 1
return infectious

def inf_per_count_time(infectious):
inf_count = {}
c_1 = 0
for t in range(1, np.max(infectious.values())+1): #to the max time
inf_count[t] = np.sum(1 for i in infectious.values() if i==t)
return inf_count #e.g. {2:3} means 3 people became infectious at time 2

#This was the accepted answer from the epidemic simulation question
def write_powerlaw_epidemics(susc, trans, inf_period, eps, repetitions, full_mat):
epi_list = []
count_list = []
new_susc = []
new_trans = []
new_inf_period = []
new_eps = []
new_count_list = []
count = 0

parameters_product = itertools.product(trans, inf_period, eps)
for transmissibility, infectious_period, epsilon in parameters_product:
while True:
for rep in range(repetitions):
for _ in range(MAX_FAILED_ATTEMPS):
g1 = powerlaw_epidemic(
pop, susc, transmissibility,
infectious_period, epsilon, full_mat)
if len(g1) >= 10 and max(g1.values()) >= 10:
g2 = inf_per_count_time(g1)
print 'parameters', transmissibility, infectious_period, epsilon
count += 1
epi_list.append(g1)
count_list.append(g2)
new_susc.append(susc)
new_trans.append(transmissibility)
new_inf_period.append(infectious_period)
new_eps.append(epsilon)
break #epi was successful, go to next rep
else:
transmissibility += 1
# Cleanup because we failed too many times
# Note: new_susc[-0:] is not always an empty list, so change if location
if rep > 0:

del epi_list[-rep:]
del count_list[-rep:]
del new_susc[-rep:]
del new_trans[-rep:]
del new_inf_period[-rep:]
del new_eps[-rep:]

count -=1
break #epi failed, break out of reps and start new para set
else:
# do not restart if we made it through the whole repetitions
print 'wrote', transmissibility, infectious_period, epsilon
break

paras =  np.array([
np.asarray(new_susc),
np.asarray(new_trans),
np.asarray(new_inf_period),
np.asarray(new_eps)
]).T
#This actually happens, but is commented out so no kind person answering this question gets files written to their machine.
#print 'number of parameter rows', paras[:,0].shape
#with open('Dec20_powerlaw_parameters.txt', 'w') as newfile1:
#    np.savetxt(newfile1, paras, fmt = ['%f', '%f', '%f', '%f'])

print 'final count', count
if __name__ == "__main__":
pop = 100
susc = 0.3
powerlaw_trans = [1.5, 2.5, 3]
inf_period = [2, 3]
eps = [0, 0.01, 0.02, 0.05]
meanval = np.array([0., 0.])
sigmaval = np.array([[1., 0.], [0., 1.]])
full_mat = create_pop(meanval, sigmaval, pop)
reps = 2

write_powerlaw_epidemics(susc, powerlaw_trans, inf_period, eps, reps, full_mat)

• Out of curiosity, did you try to implement the advices of the second solution and see if it makes some differences? – 301_Moved_Permanently Dec 20 '16 at 20:21
• Hi @MathiasEttinger, I tried, but couldn't get it to run, so I accepted the answer I could run. – StatsSorceress Dec 20 '16 at 20:27
• Ok. But at least it taught you about cProfile ;). We now have real figures to reason about. And as I said, your bottleneck is power_law_epidemic as it more than 99% of the running time of write_powerlaw_epidemics. Total time of create_pop seems also unecessaryly large. – 301_Moved_Permanently Dec 20 '16 at 21:01
• Yes, you taught me lots of useful stuff! Thanks again. I'm fiddling with power_law_epidemic but it seems storing the matrix outside the loop and just accessing inside the loop is making it worse....Ugh, n00b! – StatsSorceress Dec 20 '16 at 21:07
• That last edit looks like it's a comment on the answer. I rolled it back; please avoid editing a question in response to answers. See help/someone-answers for all the details. – Mathieu Guindon Dec 21 '16 at 17:05

Indeed, you seem to not use lists as you should. You seems to not use dictionary as you shoud either. In fact, you picked a too complicated data structure for the task at hand.

The aim of your dictionary is to store the date at which people got infected based on the index of the person. The aim of your list is to store a list of uninfected indexes. This is way too contrived as you could associate indexes to a value using just a list (and enumerate).

If you use only a single list, storing for each index the date at which the person (for this index) got infected, you can simplify your computation a lot. On top of that, when computing the people infected at time t, you include contagious people that were computed in the range $[t - \text{inf_period}; t[$, and to do so, you explicitly exclude people contaminated at time t. Instead, you could compute the list of infectious people before iterating over susceptible to avoid to:

1. Exclude those contaminated at time t;
2. Recomputing the same list over and over for each people in the population.

In the following rewrite, I use a list of infection dates, the value 0 meaning the person never got infected yet:

def powerlaw_epidemic(pop, susc, trans, inf_period, eps, full_mat):
"""Compute a list of infection date for each people in the population.
A value of 0 meaning the person never got infected during the simulation.
"""

infection_dates = np.zeros(pop, dtype=int)
initial_patient = np.random.randint(1, pop)
infection_dates[initial_patient] = 1  # First got infected on day 1

for day in itertools.count(2):
if infection_dates.all():
# Stop the simulation if everyone was infected
return infection_dates

still_virulent = max(day - inf_period, 1)
infectious_people = infection_dates >= still_virulent

if not infectious_people.any():
# Stop the simulation if noone is infectious anymore
return infection_dates

# Compute newly infected persons
for person, infected_at in enumerate(infection_dates):
if not infected_at:
dist = full_mat[person, infectious_people] ** (-trans)
probinf = 1 - np.exp(-susc * dist.sum() + eps)
if probinf >= np.random.uniform(0, 1):
infection_dates[person] = day


A few things to note:

• You don't need to use a test to differentiate the first day from the other ones, just do the special case before the loop to avoid loosing time in a useless test in the loop.
• Since you use numpy in various places, you can improve the computation a bit, especially extracting out specific indexes out of full_mat or computing the indexes of infectious_people.

You may be able to use numpy features to improve the update of infection_dates (instead of using a for loop), but I don't know numpy enough to provide them.

Now you may notice that I return an array where you returned a dictionnary, this means you have to adapt your calling code to handle the change. inf_per_count_time is easy to adapt as it only counts the number of people infected on a given day: this calls for collections.Counter. But you will also need to adapt your condition to know if a simulation was "virulent enough". Better write a function for that:

def is_virulent_enough(simulation, infected=10, duration=10):
return (simulation > 0).sum() >= infected and simulation.max() >= duration


create_pop could also be simplified as zip and .T cancel each other; making some re-creation of data superfluous.

Using various tips that I already gave in my previous answers, you can end up with the following code:

import csv
import itertools
import collections
import argparse

import numpy as np
import scipy.spatial

MAX_FAILED_ATTEMPTS = 50
EpidemyStatistics = collections.namedtuple(
'EpidemyStatistics',
'susceptibility transmissibility infectious epsilon count infected')

def create_population(pop, meanval=None, sigmaval=None):
if meanval is None:
meanval = np.array([0., 0.])
if sigmaval is None:
sigmaval = np.array([[1., 0.], [0., 1.]])

xyarray = np.random.multivariate_normal(meanval, sigmaval, pop)
pdistance = scipy.spatial.distance.pdist(xyarray)
full_mat = scipy.spatial.distance.squareform(pdistance)

with open('20Dec_population.txt', 'w') as f:
np.savetxt(f, xyarray, fmt=['%f', '%f'])

return full_mat

def powerlaw_epidemic(pop, susc, trans, inf_period, eps, full_mat):
"""Compute a list of infection date for each people in the population.
A value of 0 meaning the person never got infected during the simulation.
"""

infection_dates = np.zeros(pop, dtype=int)
initial_patient = np.random.randint(1, pop)
infection_dates[initial_patient] = 1  # First got infected on day 1

for day in itertools.count(2):
if infection_dates.all():
# Stop the simulation if everyone was infected
return infection_dates

still_virulent = max(day - inf_period, 1)
infectious_people = infection_dates >= still_virulent

if not infectious_people.any():
# Stop the simulation if noone is infectious anymore
return infection_dates

# Compute newly infected persons
for person, infected_at in enumerate(infection_dates):
if not infected_at:
dist = full_mat[person, infectious_people] ** (-trans)
probinf = 1 - np.exp(-susc * dist.sum() + eps)
if probinf >= np.random.uniform(0, 1):
infection_dates[person] = day

def is_virulent_enough(simulation, infected=10, duration=10):
return (simulation > 0).sum() >= infected and simulation.max() >= duration

def write_powerlaw_epidemics(population, susc, trans, inf_period, eps, repetitions, full_mat):
epidemies = []
count = 0

parameters_product = itertools.product(susc, trans, inf_period, eps)
for susceptibility, transmissibility, infectious_period, epsilon in parameters_product:
while True:
for rep in range(repetitions):
for _ in range(MAX_FAILED_ATTEMPTS):
infection_dates = powerlaw_epidemic(
population, susceptibility, transmissibility,
infectious_period, epsilon, full_mat)

if is_virulent_enough(infection_dates):
epidemies.append(EpidemyStatistics(
susceptibility,
transmissibility,
infectious_period,
epsilon,
collections.Counter(infection_dates),
infection_dates.tolist()))
count += 1
break
else:
transmissibility += 1

# Cleanup because we failed too many times
# Beware of del epidemies[-0:] which delete the whole list
if rep:
del epidemies[-rep:]
count -=1
break
else:
# do not restart if we made it through the whole repetitions
break

print 'number of parameter rows', len(epidemies)
with open('parameters.txt', 'w') as newfile1:
writer = csv.writer(newfile1, delimiter=' ')
writer.writerows(epidemies)

print count

if __name__ == "__main__":
parser = argparse.ArgumentParser(description='Some infos here')
parser.add_argument('-t', '--transmissibility', type=float, nargs='+', default=[1.5, 2.5, 3])
parser.add_argument('-i', '--infectious-period', type=int, nargs='+', default=[2, 3])
parser.add_argument('-e', '--epsilon', type=float, nargs='+', default=[0, 0.01, 0.02, 0.05])
args = parser.parse_args()

full_mat = create_population(args.population)
write_powerlaw_epidemics(
args.population,
args.susceptibility,
args.transmissibility,
args.infectious_period,
args.epsilon,
args.repetitions,
full_mat)

• Thank you @Mathias Ettinger, but this throws an overflow error: 21Dec_so.py:59: RuntimeWarning: overflow encountered in power dist = full_mat[person, infectious_people] ** (-trans) Any idea why the error would be thrown now, if it wasn't before? I looked it up, and it means the result is too big to store. – StatsSorceress Dec 21 '16 at 15:15
• Cool, never heard of that! [Here you go] (pastebin.com/HvsKar9L). A brief description: fullmat is a matrix of distances between individuals. The diagonal is 0, because the distance from a person to themself is 0. I only really need to store the upper (or lower) triangular matrix, but I had trouble getting the positioning of individuals to work in that case. That is, I needed the distance from either individual 1 to individual 2 OR individual 2 to individual 1, and I couldn't get that to work without using the full matrix. – StatsSorceress Dec 21 '16 at 15:25
• Ah - the link only has a population of 5 individuals to show the concept. – StatsSorceress Dec 21 '16 at 15:27
• @StatsSorceress I got no error running a full 9600 simulations for both 100 (33 seconds) and 625 (4 min 16) population. The only thing I changed was using duration=2 when checking if the simulation was successful instead of duration=10 as I saw nobody getting infected past day 6 or 7. So in your case, it must mean the transmissibility parameter kept growing and growing until generating the RuntimeWarning. I’ll make some further testing comparing your original code and the rewrite as I don't quite understand why, even with large transmissibility the number of days doesn't evolve much. – 301_Moved_Permanently Dec 21 '16 at 22:03
• @StatsSorceress I tried with your original code, the simulations doesn't seem to reach day 10 easily either. You may want to either adjust your formulas, reduce the step of increase of transmissibility or reduce your upper bound of required duration. – 301_Moved_Permanently Dec 21 '16 at 22:56