-
Notifications
You must be signed in to change notification settings - Fork 1
/
foxes_in_parallel.py
181 lines (145 loc) · 6.04 KB
/
foxes_in_parallel.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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
"""
parallel program
To run in your computer write in the terminal
mpiexec -n 4 python [file_name.py]
where 4 is the arbitrary number of processors.
Must install mpi4py using (pip install mpi4py) in terminal
"""
import numpy as np
from matplotlib import pyplot as plt
import random
#random.seed(1) # so results don't change every time I execute
from mpi4py import MPI
from mpi4py.MPI import ANY_SOURCE
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()
wt = MPI.Wtime()
k1 = 0.015
k2 = 0.00004
k3 = 0.0004
k4 = 0.04
end_time = 600
def get_rates(rabbits, foxes):
"""
Return the rates (expected events per day) as a tuple:
(rabbit_birth, rabbit_death, fox_birth, fox_death)
"""
rabbit_birth = k1 * rabbits
rabbit_death = k2 * rabbits * foxes
fox_birth = k3 * rabbits * foxes
fox_death = k4 * foxes
return (rabbit_birth, rabbit_death, fox_birth, fox_death)
dead_foxes = 0
dead_everything = 0
runs = int(10000/size)
second_peak_times = []
second_peak_foxes = []
mean_times = np.zeros(runs)
mean_foxes = np.zeros(runs)
upper_quartile_times = np.zeros(runs)
lower_quartile_times = np.zeros(runs)
upper_quartile_foxes = np.zeros(runs)
lower_quartile_foxes = np.zeros(runs)
for run in range(runs):
time = 0
rabbit = 400
fox = 200
# we don't know how long these will be so start as lists and convert to arrays later
times = []
rabbits = []
foxes = []
while time < end_time:
times.append(time)
rabbits.append(rabbit)
foxes.append(fox)
(rabbit_birth, rabbit_death, fox_birth, fox_death) = rates = get_rates(rabbit, fox)
sum_rates = sum(rates)
if sum_rates == 0:
# print("everything dead at t=",time)
dead_everything += 1
times.append(end_time)
rabbits.append(rabbit)
foxes.append(fox)
break
wait_time = random.expovariate( sum_rates )
time += wait_time
choice = random.uniform(0, sum_rates)
# Imagine we threw a dart at a number line with span (0, sum_rates) and it hit at "choice"
# Foxes change more often than rabbits, so we'll be faster if we check them first!
choice -= fox_birth
if choice < 0:
fox += 1 # fox born
continue
choice -= fox_death
if choice < 0:
fox -= 1 # fox died
if fox == 0:
#print("Foxes all died at t=",time)
dead_foxes += 1
## Break here to speed things up (and not track the growing rabbit population)
continue
if choice < rabbit_birth:
rabbit += 1 # rabbit born
continue
rabbit -= 1 # rabbit died
times = np.array(times)
rabbits = np.array(rabbits)
foxes = np.array(foxes)
index_of_second_peak = np.argmax(foxes*(times>200)*(foxes>100))
if index_of_second_peak:
second_peak_times.append(times[index_of_second_peak])
second_peak_foxes.append(foxes[index_of_second_peak])
if len(second_peak_times)>0:
mean_times[run] = np.mean(second_peak_times)
mean_foxes[run] = np.mean(second_peak_foxes)
upper_quartile_times[run] = np.percentile(second_peak_times,75)
lower_quartile_times[run] = np.percentile(second_peak_times,25)
upper_quartile_foxes[run] = np.percentile(second_peak_foxes,75)
lower_quartile_foxes[run] = np.percentile(second_peak_foxes,25)
# We don't want to plot too many lines, but would be fun to see a few
if run < 20 :
plt.plot(times, rabbits, 'b')
plt.plot(times, foxes, 'g')
comm.Barrier()
if rank == 0:
wt = MPI.Wtime() - wt
print ("run time {:.2f} with number of runs {} and number of processors {}".format(wt, runs*size, size))
comm.Barrier()
plt.legend(['rabbits','foxes'],loc="best") # put the legend at the best location to avoid overlapping things
plt.ylim(0,3000)
plt.title("This plot from processor {}".format(rank))
plt.show()
total_mean_times = np.zeros(1)
comm.Reduce(mean_times[-1] , total_mean_times , op=MPI.SUM , root=0)
total_mean_times = total_mean_times /size
total_mean_foxes = np.zeros(1)
comm.Reduce(mean_foxes[-1] , total_mean_foxes , op=MPI.SUM , root=0)
total_mean_foxes = total_mean_foxes /size
total_lower_quartile_times = np.zeros(1)
comm.Reduce(lower_quartile_times[-1] , total_lower_quartile_times , op=MPI.SUM , root=0)
total_lower_quartile_times = total_lower_quartile_times /size
total_upper_quartile_times = np.zeros(1)
comm.Reduce(upper_quartile_times[-1] , total_upper_quartile_times , op=MPI.SUM , root=0)
total_upper_quartile_times = total_upper_quartile_times /size
total_lower_quartile_foxes = np.zeros(1)
comm.Reduce(lower_quartile_foxes[-1] , total_lower_quartile_foxes , op=MPI.SUM , root=0)
total_lower_quartile_foxes = total_lower_quartile_foxes /size
total_upper_quartile_foxes = np.zeros(1)
comm.Reduce(upper_quartile_foxes[-1] , total_upper_quartile_foxes , op=MPI.SUM , root=0)
total_upper_quartile_foxes = total_upper_quartile_foxes /size
total_runs = np.zeros(1)
runs_MPI = np.ones(1)* runs
comm.Reduce(runs_MPI , total_runs , op=MPI.SUM , root=0)
total_dead_foxes = np.zeros(1)
dead_foxes_MPI = np.ones(1)* dead_foxes
comm.Reduce(dead_foxes_MPI , total_dead_foxes , op=MPI.SUM , root=0)
total_dead_everything = np.zeros(1)
dead_everything_MPI = np.ones(1)* dead_everything
comm.Reduce(dead_everything_MPI , total_dead_everything , op=MPI.SUM , root=0)
if rank ==0:
print("Number of total runs {}".format(total_runs))
print("Everything died {} times out of {} or {:.1f}%".format(total_dead_everything, total_runs, 100*total_dead_everything[0]/total_runs[0]))
print("Foxes died {} times out of {} or {:.1f}%".format(total_dead_foxes, total_runs, 100*total_dead_foxes[0]/total_runs[0]))
print("Second peak (days) is {:.1f} with IQR [{:.1f}-{:.1f}] ".format(total_mean_times[0], total_lower_quartile_times[0], total_upper_quartile_times[0]))
print("Second peak (foxes) is {:.1f} with IQR [{:.1f}-{:.1f}] ".format(total_mean_foxes[0], total_lower_quartile_foxes[0], total_upper_quartile_foxes[0]))