Post by David Drake on Mar 17, 2020 11:03:22 GMT -5
Just as an interesting exercise, I created the code below to simulate the effects of different factors on the spread of an infectious disease. This isn't scientific and it wasn't designed by an epidemiologist, but it is interesting to observe what happens as you vary population density, how contagious it is, fraction of the population that stays put, the length of time someone can spread the disease, and so forth. I think it helps emphasize the value of not moving around.
It isn't pretty or efficient, but it works. Enjoy!
It isn't pretty or efficient, but it works. Enjoy!
population=100 'population
seeds=5 'number of starting sick
wellCount=population-seeds
transferRate=0.50'probability of contracting on contact
recoCount=0'number of sims that were sick but recovered
fractionStationary=0.75
daysInfectious=21
prompt "Population: ";population
prompt "Seeds: ";seeds
prompt "Stationary fraction: ";fractionStationary
prompt "Transfer on contact fraction: ";transferRate
prompt "Days infectious: ";daysInfectious
sickCount=seeds
'nomainwin
UpperLeftX=1
UpperLeftY=1
WindowWidth=800
WindowHeight=600
graphicbox #1.g, 1, 1, 500, 500
graphicbox #1.t1, 1, 501, 800, 100
graphicbox #1.t2, 501, 1, 300, 500
open "Sim" for graphics_nf_nsb as #1
#1.g "fill white;flush"
#1.t1 "fill darkblue;flush"
#1.t2 "fill black;flush"
#1 "trapclose [q]"
dim sim(500,10)'500 sims, 10 parameters each
'1=x location
'2=y location
'3=x velocity
'4=y velocity
'5=wellness parameter (0=well, 1=sick, 2=recovered)
'6=time since initial infection
'set initial paramaters
'random location
for a =1 to population
sim(a,1)=int(rnd(1)*450)+20
sim(a,2)=int(rnd(1)*450)+20
sim(a,3)=(0.5-rnd(rnd(1)))*5
if sim(a,3)<1 and sim(a,3)>-1 then
if sim(a,3)<1 then sim(a,3)=sim(a,3)-2
if sim(a,3)>1 then sim(a,3)=sim(a,3)+2
end if
sim(a,4)=(0.5-rnd(rnd(1)))*5
if sim(a,4)<1 and sim(a,4)>-1 then
if sim(a,4)<1 then sim(a,4)=sim(a,4)-2
if sim(a,4)>1 then sim(a,4)=sim(a,4)+2
end if
if a/population<fractionStationary then sim(a,3)=0:sim(a,4)=0
next a
'make sprites
#1.t2 "down;place 1 1;color white;backcolor white;boxfilled 10 10;place 5 5;backcolor black;color black;circlefilled 5;place 5 15;color green;backcolor green;circlefilled 5"
#1.t2 "getbmp well 1 1 8 20"
#1.t2 "color red;backcolor red;circlefilled 4"
#1.t2 "getbmp sick 1 1 8 20"
#1.t2 "color blue;backcolor blue;circlefilled 4"
#1.t2 "getbmp reco 1 1 8 20"
#1.g "getbmp bg 1 1 500 500"
#1.g "background bg"
#1.g "drawbmp well 10 10;drawbmp sick 40 40;drawbmp reco 80 80"
for a = 1 to population
#1.g "addsprite sim";a;" well sick reco"
#1.g "spritexy sim";a;" ";sim(a,1);" ";sim(a,2)
#1.g "spritemovexy sim";a;" ";sim(a,3);" ";sim(a,4)
next a
for a = 1 to seeds ' seed sick people with random times since infaction
sim(a,5)=1
sim(a,6)=rnd(1)*14*24
#1.g "spriteimage sim";a;" sick"
next a
#1.g "drawsprites"
timer 25, [animate]
wait
[animate]
scan
if sickCount<=0 then
waitCount=waitCount+1
if waitCount>50 then
timer 0
goto [endSim]
end if
end if
hours=hours+1
#1.g "drawsprites"
for a = 1 to population
if sim(a,6)>0 then 'decrement time sick
sim(a,6)=sim(a,6)-1
if sim(a,6)<=0 then
sim(a,5)=2 'move to reco
#1.g "spriteimage sim";a;" reco"
recoCount=recoCount+1
sickCount=sickCount-1
end if
end if
#1.g "spritecollides sim";a;" col$"
#1.g "spritexy? sim";a;" x y"
if col$<>"" then
if rnd(1)<0.5 then sim(a,3)=sim(a,3)*(-1) 'change direction
if rnd(1)<0.5 then sim(a,4)=sim(a,4)*(-1) 'change direction
col=val(left$(after$(col$,"sim"),2)) 'get the colliding sim
if sim(a,5)=1 and sim(col,5)=0 and rnd(1)<transferRate then 'move from well to sick
sim(col,5)=1:sim(col,6)=daysInfectious*24
#1.g "spriteimage sim";col;" sick"
wellCount=wellCount-1
sickCount=sickCount+1
else
if sim(a,5)=0 and sim(col,5)=1 and rnd(1)<transferRate then
sim(a,5)=1:sim(a,6)=daysInfectious*24 'move from well to sick
#1.g "spriteimage sim";a;" sick"
wellCount=wellCount-1
sickCount=sickCount+1
end if
end if
end if
if x>=480 or x<=10 then sim(a,3)=sim(a,3)*(-1) 'change direction at perimeter
if y>=480 or y<=10 then sim(a,4)=sim(a,4)*(-1) 'change direction at perimeter
#1.g "spritemovexy sim";a;" ";sim(a,3);" ";sim(a,4)
next a
if sickCount>maxSickCount then maxSickCount=sickCount
#1.t1, "place 10 40;\Time: ";hours;" Well: ";wellCount;" Sick: ";sickCount;" Reco: ";recoCount;" "
#1.t2 "discard;cls"
#1.t2, "place 10 100;backcolor white;color black;boxfilled 40 0"
#1.t2, "place 40 100;backcolor white;color black;boxfilled 70 0"
#1.t2, "place 70 100;backcolor white;color black;boxfilled 100 0"
#1.t2, "place 40 100;backcolor pink;color black;boxfilled 70 ";100-(100*maxSickCount/population)
#1.t2, "place 10 100;backcolor green;color black;boxfilled 40 ";100-(100*wellCount/population)
#1.t2, "place 40 100;backcolor red;color black;boxfilled 70 ";100-(100*sickCount/population)
#1.t2, "place 70 100;backcolor blue;color black;boxfilled 100 ";100-(100*recoCount/population)
wait
[endSim]
timer 0
print "Simulation complete"
print "Population: ";population
print "Seed sick: ";seeds
print "Stationary: ";fractionStationary
print "Tran rate: ";transferRate
print "Days infec: ";daysInfectious
print
print "Hrs of end: ";hours
print "% infected: ";int(recoCount/population*100);"%"
print "Max % sick: ";int(maxSickCount/population*100);"%"
print "------------------------"
wait
[q]
timer 0
unloadbmp "well"
unloadbmp "sick"
unloadbmp "reco"
#1.g "cls;discard"
#1.t1 "cls;discard"
#1.t2 "cls;discard"
close #1
end