-
Notifications
You must be signed in to change notification settings - Fork 0
/
rabbits-and-foxes.py
107 lines (83 loc) · 2.74 KB
/
rabbits-and-foxes.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
97
98
99
100
101
102
103
104
105
106
107
from matplotlib import pyplot as plt
import operator
# # Implementing Kinetic Monte Carlo Simulations
R = 400
F = 200
k1 = 0.015
k2 = 0.00004
k3 = 0.0004
k4 = 0.04
days = 600
import random
class Population:
'''
This class represents the population that we are modeling with the rate of rabbits and foxes
dying and the current population is stored and accessible at any point
'''
def __init__(self, R, F):
self.R = R
self.F = F
def RabbitBorn(self):
return k1*self.R
def RabbitDies(self):
return k2*self.R*self.F
def FoxBorn(self):
return k3*self.F*self.R
def FoxDies(self):
return k4*self.F
def rateOfExpectedEvents(self):
return self.RabbitDies() + self.RabbitBorn() + self.FoxDies() + self.FoxBorn()
def event(self):
totalRate = self.rateOfExpectedEvents()
pRabbitDeath = self.RabbitDies()/totalRate
pRabbitBirth = self.RabbitBorn()/totalRate
pFoxDeath = self.FoxDies()/totalRate
pFoxBirth = self.FoxBorn()/totalRate
referenceList = [pRabbitDeath, pRabbitDeath + pRabbitBirth, pRabbitDeath+pRabbitBirth+pFoxDeath, 1]
u = 1 - random.uniform(0,1)
for n, probability in enumerate(referenceList):
if u <= probability:
break
if n == 0:
self.R -= 1
elif n == 1:
self.R += 1
elif n == 2:
self.F -= 1
else:
self.F += 1
n = 100
secondPeakMax = []
secondPeakTime = []
foxPopulationDied = 0
for i in range(n):
t = 0
tList = [0]
RabbitList = [R]
FoxList = [F]
Simulation = Population(R, F)
while t < days:
timing = random.uniform(0,1)
j = random.uniform(0,1)
# t += math.log(1/(1-j))/(Simulation.rateOfExpectedEvents())
t += random.expovariate(Simulation.rateOfExpectedEvents())
tList.append(t)
Simulation.event()
RabbitList.append(Simulation.R)
FoxList.append(Simulation.F)
if Simulation.F == 0:
# print('Fox Population has Died at {} days'.format(t))
break
else:
foxMaxTime, maxFoxMC = max(enumerate(FoxList[200:]), key=operator.itemgetter(1))
secondPeakMax.append(maxFoxMC)
secondPeakTime.append(tList[foxMaxTime])
continue
foxPopulationDied += 1
print('Fox Population Died out {0} times out of {1}'.format(foxPopulationDied, n))
print('Average maximum fox population (2nd Peak) was {}'.format(sum(secondPeakMax)/len(secondPeakMax)))
movingAverage = []
for i in range(len(secondPeakMax)):
movingAverage.append(sum(secondPeakMax[:i+1])/len(secondPeakMax[:i+1]))
plt.plot(movingAverage)
plt.show()