-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmodel_pathogen_prey.py
96 lines (69 loc) · 3.71 KB
/
model_pathogen_prey.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
from pcraster import *
from pcraster.framework import *
import numpy as np
class PredPreyModel(DynamicModel):
def __init__(self):
DynamicModel.__init__(self)
# set clone directly:
# size 200x200, cell size 1, west coord: 0, north coord: 200
setclone(200, 200, 1, 0, 200)
def initial(self):
# percentage of infected prey
percInf = 0.5
# create maps with that percentage of total pred/prey present
prey = uniform(1) < percPrey
self.pred = uniform(1) < percPred
# infect some of the prey
infPrey = pcrand(uniform(1) < percInf, prey)
# create a nominal prey map, consisting of infected prey (with the value of 2,
# sound prey (1), and empty cells (0)
self.prey = ifthenelse(infPrey, nominal(2), ifthenelse(prey, nominal(1), nominal(0)))
self.report(self.prey, 'outputInfectedPrey/preyInit')
self.report(self.pred, 'outputInfectedPrey/predInit')
def dynamic(self):
# find locations of all prey, infect neighbours, define the sound ones
allPrey = self.prey != 0
infPrey = pcrand(window4total(scalar(self.prey == 2)) >=1, allPrey)
soundPrey = pcrand(allPrey, ~infPrey)
# find all locations with both predator and prey
both = pcrand(self.pred, allPrey)
self.report(both, 'outputInfectedPrey/both')
# predators reproduce at these locations
self.pred = (window4total(scalar(both)) + scalar(both)) >= 1
self.report(self.pred, 'outputInfectedPrey/pred')
# find all surviving prey
surviveInfected = pcrand(infPrey, ~both)
surviveSound = pcrand(soundPrey, ~both)
# sound prey reproduces in own cell and neighbourhood, infected only in own cell
reprSoundPrey = (window4total(scalar(surviveSound)) + scalar(surviveSound)) >= 1
reprInfPrey = (scalar(surviveInfected) >= 1)
# save and report the nominal map
self.prey = ifthenelse(reprInfPrey, nominal(2), ifthenelse(reprSoundPrey, nominal(1), nominal(0)))
self.report(self.prey, 'outputInfectedPrey/prey')
#define number of Time Steps
nrOfTimeSteps=100
# define proportion (prey and predator) step size
propstepSize=0.1
# Open file to save results
f= open("DataInfectedPrey_propstep_"+str(propstepSize)+".txt","w+")
f.write("PercPrey-Ini,PercPreyInf-Ini,PercPreySound-Ini,PercPred-Ini,PercPreyInf-Final,PercPreySound-Final,PercPred-Final"+"\n")
# set percentage of infected preys
# set percentages for prey and predator populations
for percInf in np.arange(0,1.1,propstepSize):
for percPrey in np.arange(0,1,propstepSize):
for percPred in np.arange(0,1,propstepSize):
myModel = PredPreyModel()
dynamicModel = DynamicFramework(myModel,nrOfTimeSteps)
dynamicModel.run()
# extract information from the map
preyEq = readmap("outputInfectedPrey/prey0000."+str(nrOfTimeSteps))
predEq = readmap("outputInfectedPrey/pred0000."+str(nrOfTimeSteps))
# compute calculations
PercPredFinal= maptotal(scalar(predEq==1))/(200*200)
PercInfectedPreyFinal= maptotal(scalar(preyEq==2))/(200*200)
PercSoundPreyFinal= maptotal(scalar(preyEq==1))/(200*200)
# Writting individial results into the file - Initial and final conditions are saved
f.write(str(float(percPrey)) + "," + str(float(percPrey*percInf)) + "," + str(float(percPrey-percPrey*percInf)) + "," + str(float(percPred))+","+str(float(PercInfectedPreyFinal))+","+str(float(PercSoundPreyFinal))+","+str(float(PercPredFinal))+"\n")
print(str(float(percPrey)) + "," + str(float(percPrey*percInf)) + "," + str(float(percPrey-percPrey*percInf)) + "," + str(float(percPred))+","+str(float(PercInfectedPreyFinal))+","+str(float(PercSoundPreyFinal))+","+str(float(PercPredFinal))+"\n")
# close file
f.close()