Hier is het simpele python programmaatje dat die simulatie doet.
(het forum heeft niet echt een goeie format om code te tonen heb ik de indruk...)
import numpy as np
import matplotlib.pyplot as plt
import random as rn
mutationrate = 0.005
naiveR = 2.5;
immuneR = 1.6;
previousR = 1.6;
previous2R = 1.6;
oldR = 2.0;
trials = 50.0;
class MyPatient:
def __init__(self,time,virus):
self.time = time
self.virus = virus
self.immune = 0
def vaccinate(self,level):
self.immune = level
def evolve(self):
if (self.time > 1):
self.time -= 1
elif (self.time == 1):
self.time = 0
self.immune = self.virus
self.virus = 0
def contaminate(self,virus1):
if (virus1 > 0):
# virus mutation
probi = rn.random()
if (probi < mutationrate):
virus = virus1 + 1
else:
virus = virus1
else:
virus = 0
probi = rn.random()
if (virus > 0) and (self.immune == 0):
# naive patient
if (probi < naiveR/trials):
self.time = 10
self.virus = virus
elif (self.immune >= virus):
if (probi < immuneR/trials) and (virus > 0):
# contamination even if immune
self.time = 10
self.virus = virus
elif (self.immune == virus - 1):
# cross immunity of previous virus
if (probi < previousR / trials) and (virus > 0):
self.time = 10
self.virus = virus
elif (self.immune == virus - 2):
# cross immunity of previous previous virus
if (probi < previous2R / trials) and (virus > 0):
self.time = 10
self.virus = virus
else:
# distant cross immunity
if (probi < oldR / trials) and (virus > 0):
self.time = 10
self.virus = virus
nofpatients = 10000
noftime = 1000
nofcontacts = 10
nofillpatients = np.zeros([noftime,80])
patients = []
for patient in range(nofpatients):
# generate initial ill or native patients
probi = rn.random()
if (probi < 0.01):
patients.append(MyPatient(10,1))
else:
patients.append(MyPatient(0,0))
for time in range(noftime):
# counting number of ill patients
for p in patients:
nofillpatients[time,p.virus] += 1
# do contaminate others
for p in patients:
for k in range(nofcontacts):
ct = rn.randint(0,nofpatients-1)
q = patients[ct]
if (p.time < 6):
# after incubation time, contaminate
q.contaminate(p.virus)
# evolve patients
for p in patients:
p.evolve()
allill = np.sum(nofillpatients[:,1:],1)
plt.imshow(nofillpatients[::20,:])
plt.show()
plt.plot(allill)
plt.show()