-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathfitting.py
316 lines (223 loc) · 8.79 KB
/
fitting.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
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
from collections import OrderedDict
import json
import numpy as np
from scipy.special import gammaln
import pymultinest
import galaxy
def linear_prior(low, upp):
"""Linear prior
Maps x [0-1] uniformly between low and upp
Parameters
----------
low : float
lower limit of range
upp : float
upper limit of range
Returns
-------
y : function that takes one argument
x : float [between 0 and 1]
percentile
"""
y = lambda x: (upp-low) * x + low
return y
def logarithmic_prior(low, upp):
"""Logarithmic prior
Maps x [0-1] logarithmically between low and upp
Parameters
----------
low : float
lower limit of range
upp : float
upper limit of range
Returns
-------
y : function that takes one argument
x : float [between 0 and 1]
percentile
"""
l_low = np.log10(low)
l_upp = np.log10(upp)
y = lambda x: 10. ** ((l_upp-l_low) * x + l_low)
return y
class MultinestFitting(object):
def __init__(self, lines, flux, var, obssim, model_err,
likelihood='student', dlogZ_prior_range=(-0.3, 0.3), **kwargs):
"""Fitting object for ObsSim and data
Parameters
----------
lines : list of strings
nested list of strings identifing emission lines
e.g. lines = [['O2-3727', 'O2-3729'], ['Hg'], ['Hb'], ['O3-5007']]
here 'O2-3727', 'O2-3729' will be coadded
flux : (Nbins*Nlines) array of floats
line fluxes in [10^-20 erg/s/cm^2], Nlines = len(lines)
var : (Nbins*Nlines) array of floats
associated variances of flux
obssim : instance of BaseObsSim or subclass
model to be fitted
model_err : array of floats
systematic error to add to model fluxes
array must match flattened line list
likelihood : string [default='student']
type of likelihood function to use
(options: 'student', 'gaussian')
Keywords
--------
dof : int [default=3]
if using likelihood='student', set the Degrees of Freedom to use
"""
if flux.shape != var.shape:
raise Exception("flux and var arrays must be the same shape")
if flux.shape[1] != len(lines):
raise Exception("2nd dimension of flux must match len(lines)")
if likelihood == 'student':
self.loglikelihood = self._student_loglikelihood
try:
self.dof = float(kwargs['dof'])
except KeyError:
self.dof = 3.
elif likelihood == 'gaussian':
self.loglikelihood = self._gaussian_loglikelihood
else:
raise Exception("likelihood function '%s' unknown" % likelihood)
self.dlogZ_prior_range = np.array(dlogZ_prior_range)
self.obs_flux = flux
self.obs_var = var
self.obssim = obssim
self.params = OrderedDict()
if type(obssim.galaxy) == galaxy.GalaxyDisc:
self.init_galdisc_params()
elif type(obssim.galaxy) == galaxy.GalaxyDiscFixedrd:
self.init_galdisc_fixedrd_params()
elif type(obssim.galaxy) == galaxy.GalaxyMap:
self.init_galmap_params()
else:
raise Exception('Galaxy model %s unknown' % str(obssim.galaxy))
self.init_basegal_params()
lines, line_mapping = self.parse_lines(lines)
if len(model_err) != line_mapping.shape[1]:
raise Exception("len(model_err) must equal number of lines")
self.lines = lines
self.line_mapping_flux = line_mapping
self.line_mapping_var = (line_mapping * model_err) ** 2.
@staticmethod
def parse_lines(lines):
"""Given a nested line list create a flat list of lines and associated mapping matrix
"""
#flatten linelist and MAKE IT UNIQUE
lines_flat = OrderedDict()
for line_components in lines:
for l in line_components:
lines_flat[l] = None
lines_flat = tuple(lines_flat.keys())
mapping = np.zeros([len(lines),len(lines_flat)], dtype=float)
for i_out, line_components in enumerate(lines):
for l in line_components:
j_in = lines_flat.index(l)
mapping[i_out,j_in] = 1.
return lines_flat, mapping
def _student_loglikelihood(self, obs_flux, model_flux, obs_var, model_var):
"""Student log-likelihood function
Parameters
----------
obs_flux : array of floats
observed flux
model_flux : array of floats
model flux
obs_var : array of floats
observed variance
model_var : array of floats
model variance
Returns
-------
out : float
log-likelihood
"""
mask = np.isfinite(obs_flux) # NaNs are below snr thresh
res = obs_flux[mask] - model_flux[mask]
var = obs_var[mask] + model_var[mask]
dof = self.dof
scale_sqd = (dof-2.) * var / dof
out = ((gammaln((dof+1.)/2.) - gammaln(dof/2.)) *
np.ones_like(var, dtype=float))
out -= 0.5 * np.log(np.pi * dof * scale_sqd)
out -= (dof+1.)/2. * np.log(1. + 1./dof * res**2. / scale_sqd)
out = out.sum()
return out
def _gaussian_loglikelihood(self, obs_flux, model_flux, obs_var, model_var):
"""Gaussian log-likelihood function
Parameters
----------
obs_flux : array of floats
observed flux
model_flux : array of floats
model flux
obs_var : array of floats
observed variance
model_var : array of floats
model variance
Returns
-------
out : float
log-likelihood
"""
mask = np.isfinite(obs_flux)
res = obs_flux[mask] - model_flux[mask]
var = obs_var[mask] + model_var[mask]
n = res.size
out = -n/2. * np.log(2.*np.pi)
out -= np.log(var).sum() / 2.
out -= 1./2. * ((res**2.) / var).sum()
return out
def init_galdisc_params(self):
self.params['SFRtotal'] = logarithmic_prior(0.01, 100)
self.params['r_d'] = linear_prior(0., 2.)
def init_galdisc_fixedrd_params(self):
self.params['SFRtotal'] = logarithmic_prior(0.01, 100)
def init_galmap_params(self):
self.params['SFRtotal'] = logarithmic_prior(0.01, 100)
def init_basegal_params(self):
#get fluxgrid range in logZ and logU
fluxgrid = self.obssim.galaxy.fluxgrid
logZ_min = fluxgrid.logZ_min
logZ_max = fluxgrid.logZ_max
logU_min = fluxgrid.logU_min
logU_max = fluxgrid.logU_max
logU_0_max = logU_max + 0.8 * logZ_max
logU_0_min = logU_min + 0.8 * logZ_min
self.params['logZ_0'] = linear_prior(logZ_min, logZ_max)
dlogZ_min, dlogZ_max = self.dlogZ_prior_range
self.params['dlogZ'] = linear_prior(dlogZ_min, dlogZ_max)
self.params['logU_sol'] = linear_prior(logU_0_min, logU_0_max)
self.params['tauV_0'] = linear_prior(0., 4.)
def multinest_prior(self, cube, ndim, nparams):
for i_param, prior_func in enumerate(self.params.itervalues()):
cube[i_param] = prior_func(cube[i_param])
def cube_to_params(self, cube):
params = OrderedDict(zip(self.params, cube))
return params
def multinest_loglike(self, cube, ndim, nparams):
model_flux, model_var = self.model(cube)
loglike = self.loglikelihood(self.obs_flux, model_flux,
self.obs_var, model_var)
return loglike
def model(self, cube):
params = self.cube_to_params(cube)
flux = self.obssim(self.lines, params)
flux /= 1e-20
model_flux = np.dot(self.line_mapping_flux, flux.T).T
model_var = np.dot(self.line_mapping_var, (flux**2.).T).T
return model_flux, model_var
def multinest_run(self, basename, sampling_efficiency='parameter',
n_live_points=1000, verbose=True):
parameters = [i for i in self.params.iterkeys()]
n_params = len(parameters)
pymultinest.run(self.multinest_loglike, self.multinest_prior, n_params,
outputfiles_basename=basename,
importance_nested_sampling=False, multimodal=True,
resume=False, verbose=verbose,
n_live_points=n_live_points,
sampling_efficiency=sampling_efficiency)
# save parameter names
json.dump(parameters, open(basename+'params.json', 'w'))