-
Notifications
You must be signed in to change notification settings - Fork 0
/
cliques_MCR_mu_vs_sigma_life.py
260 lines (171 loc) · 8.34 KB
/
cliques_MCR_mu_vs_sigma_life.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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
# Written 1/8/14 by dh4gan
# Generates a collection of civilisations with some distribution in space and arrival time
# Then calculates causally connected groups
# Repeats this to allow MCR analysis
# Repeats MCR analysis over various values of the mean lifetime
import galaxy as gal
import params
import numpy as np
import scipy.ndimage as ndimage
import matplotlib.pyplot as plt
from sys import argv
from os.path import isfile
# Read in input parameters from file
if len(argv)==1:
paramfile = raw_input("Enter the parameter file: ")
elif len(argv)==2:
paramfile = argv[1]
nruns,nciv, rinner,router,rscale, zmin,zmax, mu_life,sigma_life,mu_t,sigma_t = params.read_parameters_mcr(paramfile)
iseed = -47
# Override definitions of mu_life, sigma_life, specify the values to be explored
nmu = 100
mu_life_values = np.linspace(0.001, 5.0, num=nmu)
nsigma = 100
sigma_life_values = np.linspace(1.0e-1,1.0e0, num=nsigma)
print "Surveying following parameter space:"
print "mu: \n",mu_life_values
print "sigma: \n",sigma_life_values
# Arrays to hold values of outputs for plotting
mean_group_values = np.zeros((nmu,nsigma))
sd_group_values = np.zeros((nmu,nsigma))
mean_size_values = np.zeros((nmu,nsigma))
sd_size_values = np.zeros((nmu,nsigma))
mean_arrive_values = np.zeros((nmu,nsigma))
sd_arrive_values =np.zeros((nmu,nsigma))
# Set up MCR output file
outputfile = 'output.MCR'
# Detect if this file exists - if it does, and user wishes, then simply read out the data and skip the calculation
plot_instead = "n"
if(isfile(outputfile)):
plot_instead = raw_input("Outputfile "+ outputfile+ " detected - do you want to skip calculation and plot this instead? ")
if plot_instead=="n":
print 'Outputting to new file ', outputfile
f_obj = open(outputfile,'w')
line = 'Nciv meanLife sdLife meanArrive sdArrive '
line += 'meanNgroup sdNgroup meanGroupSize sdGroupSize meanGroupArrive sdGroupArrive \n'
f_obj.write(line)
imu = 0
for mu_life in mu_life_values:
isigma = 0
for sigma_lifefrac in sigma_life_values:
sigma_life = sigma_lifefrac*mu_life
#sigma_life = sigma_lifefrac
# Now loop over runs
mcr_ngroups = np.empty(0)
mcr_groupcounts = np.empty(0)
mcr_groupsizes = np.empty(0)
mcr_grouparrive = np.empty(0)
for irun in range(nruns):
#print 'Beginning run ', irun+1
# Generate a galaxy of civilisations
myGalaxy = gal.galaxy(nciv,iseed)
# Spatial distribution of civilisations
#myGalaxy.generate_uniform_GHZ(rinner, router, zmin, zmax)
myGalaxy.generate_exponential_GHZ(rinner, router,rscale, zmin, zmax)
# Distribution of civilisation arrival times
#myGalaxy.generate_uniform_arrival_time(tmin, tmax)
myGalaxy.generate_gaussian_arrival_time(mu_t,sigma_t)
# Distribution of civilisation lifetimes
#myGalaxy.generate_fixed_lifetimes(mu_life)
myGalaxy.generate_gaussian_lifetimes(mu_life, sigma_life)
# Arrange civilisations by arrival time and construct groups
myGalaxy.sort_by_arrival_time()
myGalaxy.check_for_groups()
myGalaxy.count_groups()
myGalaxy.get_group_sizes()
#print "There are ", myGalaxy.ngroups, " groups, with latest group established by civilisation ", myGalaxy.groupmax
# Plot this galaxy and its groups
#myGalaxy.arrival_time_histogram()
#myGalaxy.plot_spatial2D()
# Write data to outputfile
#outputfile = 'output.'+str(irun)
#myGalaxy.output_group_statistics(outputfile)
#print 'Run ',irun+1, ' complete'
# Store the data in an array?
mcr_ngroups = np.append(mcr_ngroups,np.asarray(myGalaxy.ngroups))
mcr_groupcounts = np.append(mcr_groupcounts,np.asarray(myGalaxy.groupcount))
mcr_groupsizes = np.append(mcr_groupsizes,np.asarray(myGalaxy.groupsizes))
mcr_grouparrive = np.append(mcr_grouparrive,np.asarray(myGalaxy.grouparrive))
#print 'MCR complete - now producing means and standard deviations for ', mu_life, sigma_life
# Groups with a single member have zero size - remove them
mcr_groupsizes = mcr_groupsizes[np.nonzero(mcr_groupsizes)]
# Calculate means and standard deviations
mean_ngroups = np.mean(mcr_ngroups)
sd_ngroups = np.std(mcr_ngroups)
mean_groupsize = np.mean(mcr_groupsizes)
sd_groupsize = np.std(mcr_groupsizes)
mean_grouparrive = np.mean(mcr_grouparrive)
sd_grouparrive = np.std(mcr_grouparrive)
# Write MCR data to output file
print "Writing data for mu_life ", mu_life, " sigma_life ",sigma_life, "to file ",outputfile
line = str(nciv)+'\t'+str(mu_life)+'\t'+str(sigma_life)+'\t'+str(mu_t)+'\t'+str(sigma_t)+'\t'
line += str(mean_ngroups)+'\t'+str(sd_ngroups)+'\t'+str(mean_groupsize)+'\t'+str(sd_groupsize)+'\t'+str(mean_grouparrive)+'\t'+str(sd_grouparrive)+'\n'
f_obj.write(line)
# Add data to output arrays
mean_group_values[imu,isigma] = mean_ngroups
sd_group_values[imu,isigma] = sd_ngroups
mean_size_values[imu,isigma] = mean_groupsize
sd_size_values[imu,isigma] = sd_groupsize
mean_arrive_values[imu,isigma] = mean_grouparrive
sd_arrive_values[imu,isigma] = sd_grouparrive
isigma = isigma+1
# End of loop over sigma values
imu = imu +1
# End of loop over mu values
print "MCR complete for all values of mu_life, sigma_life "
f_obj.close()
if(plot_instead=="y"):
# Read data from outputfile
f_obj= open(outputfile, "r")
# Read header line
line = f_obj.readline()
imu = 0
isigma = 0
with open(outputfile) as f_obj:
for line in f_obj:
if '#' in line: continue
if 'Nciv' in line: continue
isigma = isigma +1
if(isigma > nsigma-1):
imu = imu+1
isigma = 0
if(imu>nmu-1):
break
#for imu in range(nmu):
# isigma = 0
# for isigma in range(nsigma):
# # Read line, break up into individual values
# line = f_obj.readline()
values = str.split(line)
print imu, isigma, values
# Add data to output arrays
#if(float(values[1])!=mu_life_values[imu]):
# print "WARNING: mu mismatch: ", values[1], mu_life_values[imu]
#if(float(values[2])!=sigma_life_values[isigma]):
# print "WARNING: sigma mismatch: ", values[2], sigma_life_values[isigma]
mean_group_values[imu,isigma] = float(values[5])
sd_group_values[imu,isigma] = float(values[6])
mean_size_values[imu,isigma] = float(values[7])
sd_size_values[imu,isigma] = float(values[8])
mean_arrive_values[imu,isigma] = float(values[9])
sd_arrive_values[imu,isigma] = float(values[10])
# End of loop over sigma values
# End of loop over mu values
f_obj.close()
# Now Plot data
mean_group_values = mean_group_values.astype(int)
print mean_group_values
contourlevels = [1,10,50,100,150,200,300,400,450,490,500]
fig1 = plt.figure()
ax = fig1.add_subplot(111)
ax.set_xlabel('$\mu_{Life}$ (Myr)')
ax.set_ylabel('$\sigma_{Life}$ (Myr)')
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlim(1.0e-1,5.0e0)
#plt.hexbin(mu_life_values,sigma_life_values,C=mean_group_values.transpose(),gridsize = int(nmu*0.01), xscale='log',yscale='log',mincnt = 1)
#cb = plt.colorbar()
#cb.set_label('Mean Group Number', fontsize=20)
cs = ax.contour(mu_life_values,sigma_life_values, mean_group_values.transpose(), levels=contourlevels, colors='k')
ax.clabel(cs,contourlevels)
fig1.savefig('contour.png', format='png')