-
Notifications
You must be signed in to change notification settings - Fork 0
/
cliques_MCR.py
140 lines (93 loc) · 3.9 KB
/
cliques_MCR.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
# Written 2/7/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
import galaxy as gal
import params
import numpy as np
import matplotlib.pyplot as plt
from sys import argv
# 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
# 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
#print "Counts: ", myGalaxy.groupcount
#print "Sizes: ", myGalaxy.groupsizes
# 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'
# 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)
fig1 = plt.figure(1)
ax = fig1.add_subplot(111)
ax.set_ylabel('Relative Frequency')
ax.set_xlabel('Group Population')
ax.hist(mcr_groupcounts, bins=100, normed=True)
plt.show()
fig1 = plt.figure(1)
ax = fig1.add_subplot(111)
ax.set_ylabel('Relative Frequency')
ax.set_xlabel('Group Size (kpc)')
ax.hist(mcr_groupsizes, bins=100)
plt.show()
fig1 = plt.figure(1)
ax = fig1.add_subplot(111)
ax.set_ylabel('Relative Frequency')
ax.set_xlabel('Group Arrival Time (Myr)')
ax.hist(mcr_grouparrive, bins=100)
plt.show()
# Write MCR data to output file
outputfile = 'output.MCR'
print 'Outputting to 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)
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)
f_obj.close()