-
Notifications
You must be signed in to change notification settings - Fork 63
/
likelihood.py
92 lines (75 loc) · 3.2 KB
/
likelihood.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
from isochrones.mags import interp_mag
from .utils import fast_addmags
import numba as nb
from math import pi, log, sqrt
LOG_ONE_OVER_ROOT_2PI = log(1./sqrt(2*pi))
@nb.jit(nopython=True)
def gauss_lnprob(val, unc, model_val):
resid = val - model_val
return LOG_ONE_OVER_ROOT_2PI + log(unc) - 0.5 * resid * resid / (unc * unc)
@nb.jit(nopython=True)
def star_lnlike(pars, index_order,
spec_vals, spec_uncs,
mag_vals, mag_uncs, i_mags,
model_grid, i_Teff, i_logg, i_feh, i_Mbol,
model_ii0, model_ii1, model_ii2,
bc_grid, bc_ii0, bc_ii1, bc_ii2, bc_ii3):
n_pars = len(pars)
has_binary = False
has_triple = False
if n_pars == 5:
single_pars = [pars[0], pars[1], pars[2], pars[3], pars[4]]
elif n_pars == 6: # binary system
single_pars = [pars[0], pars[2], pars[3], pars[4], pars[5]]
binary_pars = [pars[1], pars[2], pars[3], pars[4], pars[5]]
has_binary = True
elif n_pars == 7: # triple system
single_pars = [pars[0], pars[3], pars[4], pars[5], pars[6]]
binary_pars = [pars[1], pars[3], pars[4], pars[5], pars[6]]
triple_pars = [pars[2], pars[3], pars[4], pars[5], pars[6]]
has_binary = True
has_triple = True
Teff, logg, feh, mags = interp_mag(single_pars, index_order, model_grid,
i_Teff, i_logg, i_feh, i_Mbol,
model_ii0, model_ii1, model_ii2,
bc_grid, i_mags,
bc_ii0, bc_ii1, bc_ii2, bc_ii3)
if has_binary:
_, _, _, mags_binary = interp_mag(binary_pars, index_order, model_grid,
i_Teff, i_logg, i_feh, i_Mbol,
model_ii0, model_ii1, model_ii2,
bc_grid, i_mags,
bc_ii0, bc_ii1, bc_ii2, bc_ii3)
if has_triple:
_, _, _, mags_triple = interp_mag(triple_pars, index_order, model_grid,
i_Teff, i_logg, i_feh, i_Mbol,
model_ii0, model_ii1, model_ii2,
bc_grid, i_mags,
bc_ii0, bc_ii1, bc_ii2, bc_ii3)
if n_pars == 6:
for i in range(len(mags)):
mags[i] = fast_addmags([mags[i], mags_binary[i]])
elif n_pars == 7:
for i in range(len(mags)):
mags[i] = fast_addmags([mags[i], mags_binary[i], mags_triple[i]])
lnlike = 0
# Spec_vals are Teff, logg, feh
val = spec_vals[0]
unc = spec_uncs[0]
if val == val: # Skip if nan
lnlike += gauss_lnprob(val, unc, Teff)
# logg
val = spec_vals[1]
unc = spec_uncs[1]
if val == val: # Skip if nan
lnlike += gauss_lnprob(val, unc, logg)
# feh
val = spec_vals[2]
unc = spec_uncs[2]
if val == val: # Skip if nan
lnlike += gauss_lnprob(val, unc, feh)
for i in range(len(mag_vals)):
val = mag_vals[i]
unc = mag_uncs[i]
lnlike += gauss_lnprob(val, unc, mags[i])
return lnlike