-
Notifications
You must be signed in to change notification settings - Fork 0
/
chi2.py
350 lines (303 loc) · 11.8 KB
/
chi2.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
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
""" This module defines the chi-squared and related functions
Module author: Chen Sun
Year: 2020
Email: chensun@mail.tau.ac.il
"""
import numpy as np
import model
def chi2_no_soliton(c, Rs, ups_disk, ups_bulg, gal, DM_profile="NFW"):
"""chi2 for an NFW fit (c, Rs; ups_disk, ups_bulg). Runs over a single galaxy
:param c: concentration param of the NFW profile, and delc for Burkert
:param Rs: critical radius of the NFW profile
:param ups_disk: disk surface brightness
:param ups_bulg: bulge surface brightness
:param gal: a Galaxy instance
"""
chi2 = 0
Vth2_arr = model.v2_rot(gal, c, Rs, ups_bulg, ups_disk, DM_profile)
for i, r in enumerate(gal.R):
# treat the i-th bin of the rot curve
#
# # TODO: move this part out to a dedicated function
# # V thory due to DM
# if DM_profile == "NFW":
# M_enclosed = model.M_NFW(r, Rs, c)
# elif DM_profile == "Burkert":
# M_enclosed = model.M_Burkert(
# r, delc=c, Rs=Rs)
# else:
# raise Exception(
# "Only NFW and Burkert are implemented at the moment.")
# VDM2 = model._G_Msun_over_kpc * model._c**2 * (M_enclosed/1.) * (1./r)
# # combine DM with the baryon mass model (from SB data)
# Vb2 = (ups_bulg*np.abs(gal.Vbul[i])*gal.Vbul[i]
# + ups_disk*np.abs(gal.Vdisk[i])*gal.Vdisk[i]
# + np.abs(gal.Vgas[i])*gal.Vgas[i]
# )
# Vth2 = VDM2 + Vb2
# ODOT
Vth2 = Vth2_arr[i]
if Vth2 > 0:
Vth = np.sqrt(Vth2)
else:
Vth = 0.
Vobs = gal.Vobs[i]
dVobs = gal.dVobs[i]
# compute chi2 for this bin
chi2 += (Vth - Vobs)**2/dVobs**2
# construct Vtot for visual/sanity checks
# construct Vtot
return chi2
def fit_rot_curve(gal, DM_profile='NFW', gridsize=50):
"""A quick fit with NFW/Burkert only
"""
rs_arr = np.linspace(1, 80, gridsize)
c_arr = np.logspace(0, 1.5, gridsize)
rs_mesh, c_mesh = np.meshgrid(rs_arr, c_arr, indexing='ij')
rs_flat, c_flat = rs_mesh.reshape(-1), c_mesh.reshape(-1)
chi2_flat = []
chi2_val_rec = 1e9
rec_idx = None
for i in range(len(rs_flat)):
rs = rs_flat[i]
c = c_flat[i]
chi2_val = chi2_no_soliton(c=c, Rs=rs, ups_disk=0.5, ups_bulg=0.5,
gal=gal, DM_profile=DM_profile)
chi2_flat.append(chi2_val)
if chi2_val < chi2_val_rec:
rec_idx = i
chi2_val_rec = chi2_val
# print(chi2_val)
chi2_flat = np.array(chi2_flat)
# best fit
rs = rs_flat[rec_idx]
c = c_flat[rec_idx]
# output
gal.rs = rs
gal.c = c
return (rs_mesh, c_mesh, chi2_flat)
def chi2_single_gal(m, M, c, Rs, ups_disk, ups_bulg, gal, flg_Vtot=False, DM_profile="NFW", flg_overshoot=False, combine_mode=None):
"""chi2 for a fixed theory point (m, M, c, Rs; ups_disk, ups_bulg). Runs over a single galaxy, gal
up
:param m: scalar mass [eV]
:param M: soliton mass [Msun]
:param c: concentration param of the NFW profile, and delc for Burkert
:param Rs: critical radius of the NFW profile
:param ups_disk: disk surface brightness
:param ups_bulg: bulge surface brightness
:param gal: a Galaxy instance
:param flg_Vtot: flg to signal returning of the total rot velocity
defined as Vtot**2 = Vb**2 + VNFW**2 + Vsol**2 [km/s]
:param DM_profile: flag to choose between NFW and Burkert. (Default: NFW)
:returns: chi2 value
"""
chi2 = 0
if flg_Vtot:
Vtot = np.array([])
Vth2_arr = model.v2_rot(gal, c, Rs, ups_bulg,
ups_disk, DM_profile,
m=m, M=M, combine_mode=combine_mode)
for i, r in enumerate(gal.R):
# treat the i-th bin of the rot curve
#
# # TODO: move this part out to a dedicated function
# # V thory due to DM
# if DM_profile == "NFW":
# M_enclosed = model.M_NFW(r, Rs, c) + model.M_sol(r, m, M)
# elif DM_profile == "Burkert":
# M_enclosed = model.M_Burkert(
# r, delc=c, Rs=Rs) + model.M_sol(r, m, M)
# else:
# raise Exception(
# "Only NFW and Burkert are implemented at the moment.")
# VDM2 = model._G_Msun_over_kpc * model._c**2 * (M_enclosed/1.) * (1./r)
# # combine DM with the baryon mass model (from SB data)
# Vb2 = (ups_bulg*np.abs(gal.Vbul[i])*gal.Vbul[i]
# + ups_disk*np.abs(gal.Vdisk[i])*gal.Vdisk[i]
# + np.abs(gal.Vgas[i])*gal.Vgas[i]
# )
# Vth2 = VDM2 + Vb2
# ODOT
Vth2 = Vth2_arr[i]
if Vth2 > 0:
Vth = np.sqrt(Vth2)
else:
Vth = 0.
Vobs = gal.Vobs[i]
dVobs = gal.dVobs[i]
# compute chi2 for this bin
if not flg_overshoot:
chi2 += (Vth - Vobs)**2/dVobs**2
else:
if Vth - Vobs > 0:
chi2 += (Vth - Vobs)**2/dVobs**2
# construct Vtot for visual/sanity checks
# construct Vtot
if flg_Vtot:
Vtot = np.append(Vtot, Vth)
if flg_Vtot:
return (chi2, Vtot)
else:
return chi2
def chi2_single_gal_overshooting(m, M, ups_disk, ups_bulg, gal, flg_Vtot=False):
"""chi2 for a fixed theory point (m, M; ups_disk, ups_bulg)
run over a single galaxy, gal. It only has soliton core in it and only counts
the overshooting of the data
:param m: scalar mass [eV]
:param M: soliton mass [Msun]
:param ups_disk: disk surface brightness
:param ups_bulg: bulge surface brightness
:param gal: a Galaxy instance
:param flg_Vtot: flg to signal returning of the total rot velocity
defined as Vtot**2 = Vb**2 + VNFW**2 + Vsol**2 [km/s]
:returns: chi2 value
"""
chi2 = 0
if flg_Vtot:
Vtot = np.array([])
for i, r in enumerate(gal.R):
# treat the i-th bin of the rot curve
#
# TODO: move this part out to a dedicated function
# V thory due to DM
M_enclosed = model.M_sol(r, m, M) # now only with soliton, no NFW
VDM2 = model._G_Msun_over_kpc * model._c**2 * (M_enclosed/1.) * (1./r)
# combine DM with the baryon mass model (from SB data)
Vb2 = (ups_bulg*np.abs(gal.Vbul[i])*gal.Vbul[i]
+ ups_disk*np.abs(gal.Vdisk[i])*gal.Vdisk[i]
+ np.abs(gal.Vgas[i])*gal.Vgas[i]
)
Vth2 = VDM2 + Vb2
if Vth2 > 0:
Vth = np.sqrt(Vth2)
else:
Vth = 0.
# ODOT
Vobs = gal.Vobs[i]
dVobs = gal.dVobs[i]
# compute chi2 for this bin
if Vth > Vobs:
chi2 += (Vth - Vobs)**2/dVobs**2
# construct Vtot for visual/sanity checks
# construct Vtot
if flg_Vtot:
Vtot = np.append(Vtot, Vth)
if flg_Vtot:
return (chi2, Vtot)
else:
return chi2
def chi2_gals(x, keys, keys_fixed, data, params, verbose=0):
"""This function mainly does the parsing of x from lnprob(x) and feeds each galaxy with parameters relevant to it.
:param x: The theory point whose chi2 is being computed. n-tuple
:param keys: The description of each entry of x. list of length n
:param keys_fixed: The description of each fixed entry
:param data: the galaxies to be run over
:param params: param card
:returns: chi2 value
"""
chi2_tot = 0
chi2_comp = [] # the list to save all the chi2 components
for gal in data:
# for each galaxy, we pack a dict
param_for_gal = {}
# deal with free params
for i, key in enumerate(keys):
words = key.split()
if (len(words) == 2) and (words[1] == gal.name):
# dealing with galaxy-specific variables
param_for_gal[words[0]] = x[i]
if len(words) == 1:
# dealing with universal variables
param_for_gal[words[0]] = x[i]
# now deal with fixed params
for i, key in enumerate(keys_fixed):
param_for_gal[key] = params[key+' fixed']
# from now on only read param_for_gal param card
m = 10**param_for_gal['logm']
# combining of DM components
# catch the combine mode
# try:
# combine_mode = params['combine_mode']
# except KeyError:
# combine_mode = 'matching'
combine_mode = params['combine_mode']
# DM profile
try:
DM_profile = params['DM_profile']
except KeyError:
if verbose > 0:
print('chi2.py:DM_profile is not specified, default to NFW')
DM_profile = "NFW"
try:
Rs = param_for_gal['Rs']
except Exception:
if verbose > 5:
print('-----param_for_gal= %s' % param_for_gal)
print('-----Rs is not set, so soliton is the only component')
if DM_profile != "Soliton":
raise Exception(
"You specify the DM profile to be either NFW or Burkert, but didn't supply Rs")
flg_catch_c = False
flg_catch_logc = False
try:
# in NFW use c for concentration
c = param_for_gal['c']
flg_catch_c = True
except:
pass
try:
# in Burkert profile try to catch logc for log(del_c)
c = 10**param_for_gal['logc']
flg_catch_logc = True
except:
pass
# I'm being sloppy here by postponing the raise until here when it is actually used
# let's call it a just-in-time warning :-)
if flg_catch_c and flg_catch_logc:
raise Exception('You gave both c and logc. Only one can be kept.')
if (not flg_catch_c) and (not flg_catch_logc) and (DM_profile != "Soliton"):
raise Exception(
"You specify the DM profile to be either NFW or Burkert, but didn't supply c or logc")
try:
ups_disk = param_for_gal['ups_disk']
except KeyError:
raise Exception('no ups_disk is specified.')
try:
ups_bulg = param_for_gal['ups_bulg']
except KeyError:
raise Exception('no ups_bulg is specified')
# try:
# h = param_for_gal['h']
# except KeyError:
# raise Exception('h is not specified.')
# now deal with derived params, if any, such as M
try:
M = 10**param_for_gal['logM']
except KeyError:
raise Exception('logM is not specified.')
# sum over
if (DM_profile == "NFW") or (DM_profile == "Burkert"):
chi2_i = chi2_single_gal(m=m,
M=M,
c=c,
Rs=Rs,
ups_disk=ups_disk,
ups_bulg=ups_bulg,
gal=gal,
DM_profile=DM_profile,
combine_mode=combine_mode)
elif DM_profile == "Soliton":
chi2_i = chi2_single_gal_overshooting(m=m,
M=M,
ups_disk=ups_disk,
ups_bulg=ups_bulg,
gal=gal,
flg_Vtot=False)
chi2_tot += chi2_i
# chi2_comp.append(chi2_i)
chi2_comp.append(M) # use chi2_comp to contain M temporarily
chi2_comp.insert(0, chi2_tot)
if verbose > 5:
print('-----chi2.py:chi2_comp=%s' % chi2_comp)
# return chi2_tot
return np.array(chi2_comp)