-
Notifications
You must be signed in to change notification settings - Fork 0
/
chiplotter.py
192 lines (174 loc) · 8.78 KB
/
chiplotter.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
from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter, MaxNLocator
'''
This handy program turns your ELC output file party into something more useful.
--> Makes a plot of chi^2 vs. fit parameters from a markovELC / geneticELC run
(Only parameters in the generation files with the lowest 10k chi2 values are plotted)
--> Spits out the fit parameters with a measure of uncertainty (+ and -) for all
the parameters in both the generation and ELCparm files
Required file: 'generation.all', created by 'cat generation.1* > generation.all'
Required file: 'ELCparm.all', created by 'cat ELCparm.1* > ELCparm.all'
Required file: 'key.ELCparm' (this is automatically generated by ELC)
Required file: 'gridloop.opt' (you need this to run ELC to start with)
NOTE that in the case of demcmc-ELC, substitute demcmc_fitparm.1* for generation.1*,
substitute demcmc_starparm.1* for ELCparm.1*, and also make a chi.all file from chi.1*.
IF THERE ARE TOO MANY FILES, you may need to use these commands:
find . -name "demcmc_fitparm.1*" -print0 | xargs -0 cat > fitparm.all
find . -name "demcmc_starparm.1*" -print0 | xargs -0 cat > starparm.all
find . -name "chi.1*" -print0 | xargs -0 cat > chi.all
The fit parameters in gridloop.opt are reported in the generation files.
Other parameters of interest which follow from these are reported in the ELCparm files.
The plot will be a 4 x 5 grid. If you are fitting more than 20 parameters in gridloop.opt,
the last ones will be omitted.
'''
# Important filename definitions2
gridloopfile = '../../RG_ELCmodeling/9246715/demcmc001_jerry/gridloop.opt'
generationfile = '../../RG_ELCmodeling/9246715/demcmc001_jerry/fitparm.all'
parmfile = '../../RG_ELCmodeling/9246715/demcmc001_jerry/starparm.all'
parmkeyfile = '../../RG_ELCmodeling/9246715/demcmc001_jerry/key.ELCparm'
chi2file = '../../RG_ELCmodeling/9246715/demcmc001_jerry/chi.all' # OMIT FOR MARKOVELC
outfile = '../../RG_ELCmodeling/9246715/demcmc001_jerry/chiplotout_WTF.txt'
out = open(outfile, 'w')
gridloop = [line.rstrip('\n') for line in open(gridloopfile)]
nvars = int(gridloop[10]) # reads the number of fit variables from gridloop file
print('Working, please be patient...')
print('Results will be written to {0}'.format(outfile))
# Read in names and limits for fit variables from gridloop file
varnames = []; varlower = []; varupper = []; varlist = []
for i in range(0, nvars):
varnames.append(gridloop[10+i+1].rstrip())
varlimits = gridloop[10+nvars+i+1]
values = varlimits.split()
varlower.append(float(values[0]))
varupper.append(float(values[1]))
# manually include 2 systemic velocity columns and 4 final columns (t0, tconj, ecc, argper)
varnames.append('gamma1'); varnames.append('gamma2'); varnames.append('t0v2')
varnames.append('tconjv2'); varnames.append('ecc'); varnames.append('argper')
varlower.append(-100); varlower.append(-100); varlower.append(0); varlower.append(0)
varlower.append(0); varlower.append(0)
varupper.append(100); varupper.append(100); varupper.append(1000); varupper.append(1000)
varupper.append(1); varupper.append(360)
# Read in chi^2 and parameter values from generation/fitparm file
# The number of cols in the generation/fitparm file varies for markovELC/demcmcELC
try:
varlist_gen = np.loadtxt(generationfile, usecols=(range(1,nvars+8)), dtype=np.float64, unpack=True)
chi2_gen = varlist_gen[0]
varlist_gen = np.delete(varlist_gen, 0, 0)
demcmc = False
except:
varlist_gen = np.loadtxt(generationfile, usecols=(range(0,nvars+5)), dtype=np.float64, unpack=True)
demcmc = True
print('Read in generation/fitparm file')
# Read in chi^2 and parameter values from ELCparm/starparm file
# The number of cols in the ELCparm/starparm file varies for markovELC/demcmcELC
parmkeys = np.loadtxt(parmkeyfile, comments='#', usecols=(1,), dtype={'names':('parmkeys',),'formats':('|S11',)}, unpack=True)
for idx, entry in enumerate(parmkeys): # remove 'index' from parmkeyfile
entry = str(entry)
if ('index' in entry):
parmkeys = np.delete(parmkeys, idx, axis=0) # remove 'chi^2' from parmkeyfile
for idx, entry in enumerate(parmkeys):
entry = str(entry)
if ('chi^2' in entry):
parmkeys = np.delete(parmkeys, idx, axis=0)
nparms = len(parmkeys)
if demcmc == False:
varlist_par = np.loadtxt(parmfile, usecols=(range(1,nparms)), dtype=np.float64, unpack=True)
chi2_par = varlist_par[0]
varlist_par = np.delete(varlist_par, 0, 0)
chi2s = chi2_par # higher precision than chi2_gen, but otherwise identical
elif demcmc == True:
varlist_par = np.loadtxt(parmfile, usecols=(range(0,nparms-1)), dtype=np.float64, unpack=True)
chi2s = np.loadtxt(chi2file, usecols=(1,), dtype=np.float64, unpack=True)
print('Read in ELCparm/starparm file (and chi2 file)')
# Ensure the gen/fitparm, ELC/starparm, and chi2 arrays are all the same length
newlength = min(varlist_gen.shape[1], varlist_par.shape[1], chi2s.shape[0]) # find length of shortest array
varlist_gen = varlist_gen[:,0:newlength]
varlist_par = varlist_par[:,0:newlength]
chi2s = chi2s[0:newlength]
# Sort parameter arrays by chi2, and only keep values with the lowest chi2 values
sorted_varlist_gen = []; sorted_varlist_par = []
chi2cutval = 10000
for array in varlist_gen:
sorted_varlist_gen.append(array[np.argsort(chi2s)][:chi2cutval])
for array in varlist_par:
sorted_varlist_par.append(array[np.argsort(chi2s)][:chi2cutval])
sorted_chi2s = chi2s[np.argsort(chi2s)][:chi2cutval]
deltachi = np.min(sorted_chi2s) + 1
print('Sorted by chi2 and kept only the {0} values with lowest chi2 values'.format(chi2cutval))
# Loop over generation file things (fit parameters) to define x-variables for plot
for i in range (0, nvars+5):
if demcmc == False:
i = i + 1
xvalues = sorted_varlist_gen[i-1]
elif demcmc == True:
xvalues = sorted_varlist_gen[i]
##############################
### PLOTTING STUFF IS HERE ###
if i <= 5:
ax = plt.subplot2grid((4,5),(0,i-1))
elif i > 5 and i <= 10:
ax = plt.subplot2grid((4,5),(1,i-5-1))
elif i > 10 and i <= 15:
ax = plt.subplot2grid((4,5),(2,i-10-1))
elif i > 15 and i <= 20:
ax = plt.subplot2grid((4,5),(3,i-15-1))
#elif i > 20 and i <= 25:
# ax = plt.subplot2grid((5,5),(4,i-20-1))
ax.xaxis.set_major_locator(MaxNLocator(3)) # no more than 2 ticks per axis
#ax.yaxis.set_major_formatter(FormatStrFormatter('%.3f'))
ax.xaxis.set_major_formatter(FormatStrFormatter('%.2f'))
plt.subplots_adjust(wspace=0, hspace=0.4)
if i % 5 != 1: ax.set_yticklabels(())
plt.plot(xvalues, sorted_chi2s, marker='o', color='k', mec=None, ms=2, ls='None', markevery=10)
plt.xlabel(varnames[i])
#plt.text(xmin*1.1, ymin*1.1, varnames[i-1])
#xmin = np.min(xfilter)
#xmax = np.max(xfilter)
#plt.axis([xmin, xmax, ymin, ymax])
#plt.ylim((plotymin, plotymax))
#plt.axis([np.min(errorxs), np.max(errorxs), plotymin, plotymax])
### PLOTTING STUFF IS HERE ###
##############################
# Print out fit variables from generation file with uncertainties
bestx = xvalues[0]
errorx = []
for idx, (value, chi2) in enumerate(zip(xvalues, chi2s)):
if (chi2 - chi2s[0]) <= deltachi:
errorx.append(value)
errorxplus = np.max(errorx)
errorxminus = np.min(errorx)
print(varnames[i], '\t =', bestx, '+', errorxplus-bestx, '\ -', bestx-errorxminus, file=out)
# Loop over ELCparm file things
# Print out fit variables from ELCparm file with uncertainties
# Two options for different number of cols, as before
if demcmc == False:
for i in range(1, nparms+1):
xvalues = sorted_varlist_par[i-1]
bestx = xvalues[0]
errorx = []
for idx, (value, chi2) in enumerate(zip(xvalues, chi2s)):
if (chi2 - chi2s[0]) <= deltachi:
errorx.append(value)
errorxplus = np.max(errorx)
errorxminus = np.min(errorx)
print(parmkeys[i-1][0], '\t =', bestx, '+', errorxplus-bestx, '\ -', bestx-errorxminus, file=out)
elif demcmc == True:
for i in range(0, nparms-1):
xvalues = sorted_varlist_par[i]
bestx = xvalues[0]
errorx = []
for idx, (value, chi2) in enumerate(zip(xvalues, chi2s)):
if (chi2 - chi2s[0]) <= deltachi:
errorx.append(value)
errorxplus = np.max(errorx)
errorxminus = np.min(errorx)
print(parmkeys[i][0], '\t =', bestx, '+', errorxplus-bestx, '\ -', bestx-errorxminus, file=out)
out.close()
print('Chi2 ranges from {0} to {1}'.format(np.min(chi2s), np.max(chi2s)))
print('There are {0} parameters explicitly being fit in gridloop.'.format(nvars))
print('Error bars assume a delta-chi2 threshold of {0}.'.format(deltachi))
print('Skipping a plot because that is ugly and slow.')
#print('Here comes a plot...')
#plt.show()