-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgaussian_likelihood_sym_params.py
150 lines (125 loc) · 6.2 KB
/
gaussian_likelihood_sym_params.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
import numpy as np
import scipy.stats
import emcee
import faulthandler
faulthandler.enable()
import os
import sys
nseospy_tk = os.getenv('NSEOSPY_TK')
sys.path.insert(0, nseospy_tk)
from multiprocessing import Pool
from nseospy import setup_den_tot
from nseospy import setup_den_tov
from nseospy import tovsolver
from nseospy import build_eos
np.random.seed(10)
# Array of central densities (baryons fm^-3) for EOS (logarithmic spacing)
# Solves for pressure at 200 density points and interpolates
Den = setup_den_tot(model="n-d-log", var2=0.0, n_min=2.e-7, n_max=1.7, npt=200)
# List of EoS params passed to TOV solver
list_param_sly4 = np.array([-15.97, 0.1595, 230.0, -225.0, -443.0, 32.01, 46.00, -120.0, 350.0, -690.0, 1.0000, 0.000, 0.0, 0.0, 0.0, 0.0, 0.0, 6.90, 0.00], dtype=np.float64)
list_param_qyc = np.array([250, 0.3])
list_param_pair = np.array([1.573e-03, 3.105, 8.551e-02, 1.386, 0.0, 0.0, 0.0, 0.0], dtype=np.float64)
list_ns_massRef = np.array([8.67, 68.0, 30.0], dtype=np.float64) # SLy5
list_aFsFlag = np.array([4, 1, 0], dtype=np.int32)
list_aInit = np.array([60.0, 0.167, 0.4, 0.0], dtype=np.float64)
list_flags_ns = np.array([0, 1], dtype=np.int32)
# list_fs_param = [0.967, 1.039, 1.257, 1.116, 0.980, 1.000, 1.000, 1.000]
list_fs_param = 'SLy5'
list_ns_outFileFormat = np.array([0, 0, 0, 0], dtype=np.int32)
# Define physically valid ranges of parameters
n = 20
sly4_0 = np.linspace(-16.35, -15.31, n) # E_sat
sly4_1 = np.linspace(0.145, 0.1746, n) # n_sat
sly4_2 = np.linspace(201, 338, n) # K_sat
sly4_3 = np.linspace(-570, 950, n) # Q_sat
sly4_4 = np.linspace(-903, 9997, n) # Z_sat
sly4_5 = np.linspace(26.83, 38.71, n) # E_sym
sly4_6 = np.linspace(9.9, 64, n) # L_sym
sly4_7 = np.linspace(-235, 213, n) # K_sym
sly4_8 = np.linspace(-86, 846, n) # Q_sym
sly4_9 = np.linspace(-1450, -5, n) # Z_sym
param_ranges = np.array([np.max(sly4_5)-np.min(sly4_5), np.max(sly4_6)-np.min(sly4_6), np.max(sly4_7)-np.min(sly4_7)])
# =============================================================================
# param_ranges = np.array([np.max(sly4_0)-np.min(sly4_0), np.max(sly4_1)-np.min(sly4_1)]) #,
# # np.max(sly4_2)-np.min(sly4_2), np.max(sly4_3)-np.min(sly4_3),
# # np.max(sly4_4)-np.min(sly4_4), np.max(sly4_5)-np.min(sly4_5),
# # np.max(sly4_6)-np.min(sly4_6), np.max(sly4_7)-np.min(sly4_7),
# # np.max(sly4_8)-np.min(sly4_8), np.max(sly4_9)-np.min(sly4_9)])
# =============================================================================
# Define a likelihood function that takes nuc params generated by MCMC, uses them to solve the TOV
# eqns, produces the radius of a 1.4 solar mass star and calculates the likelihood of generating the
# observed data
def ln_likelihood(params):
# Append the params not being varied
params = np.array([-15.97, 0.1595, 230.0, -225.0, -443.0, params[0], params[1], params[2], 350.0, -690.0, 1.0000, 0.000, 0.0, 0.0, 0.0, 0.0, 0.0, 6.90, 0.00], dtype=np.float64)
# Build the EoS
Eos = build_eos(
Den,
form="tov",
mf_model="mm_", # mm_ is nucleonic, qyc is quarkyonic
pair_model="no_",
ffg_type="nr",
mm_type="mmnr",
# mm_param="SLy5opt2",
mm_param=params,
qyc_type="qycis__",
qyc_param="250-3",
# qyc_param=list_param_qyc,
pair_type="pair1",
pair_bcs="bcs",
pair_param=list_param_pair,
fmuon="y",
fneutrino="n",
fbnuc="y",
fblep="y",
ns_massRef=list_ns_massRef,
aFsFlag=list_aFsFlag,
force="new",
aInit=list_aInit,
flags_ns=list_flags_ns,
fs_param=list_fs_param,
ns_outFileFormat=list_ns_outFileFormat
)
# Central density array for NS
NSnbc = setup_den_tov(model="test")
# Solves TOV eqns for each central density and solves to give masses and radii
atov = tovsolver(NSnbc, Eos)
print("Radii:", atov.rad, "in", atov.rad_unit)
print("Densities:", atov.aNSnbc)
print("Masses:", atov.m, "in", atov.m_unit)
if np.all(atov.rad):
# Use linear interpolation to determine the radius of a 1.4 solar mass star
difference_array = 1.4 - atov.m
nearest_above = np.where(difference_array < 0, difference_array, -np.inf).argmax()
nearest_below = np.where(difference_array > 0, difference_array, np.inf).argmin()
rad_at_1_point_4_msun = atov.rad[nearest_below] + (1.4-atov.m[nearest_below])*((atov.rad[nearest_above]-atov.rad[nearest_below])/(atov.m[nearest_above]-atov.m[nearest_below]))
return [np.log(scipy.stats.norm(11, 0.5).pdf(rad_at_1_point_4_msun)), np.max(atov.m)]
return -np.inf
# Prior distribution, just uniform for now
def ln_prior(params):
if np.min(sly4_5) < params[0] < np.max(sly4_5) and np.min(sly4_6) < params[1] < np.max(sly4_6) and np.min(sly4_7) < params[2] < np.max(sly4_7): # and \
# np.min(sly4_2) < params[2] < np.max(sly4_2) and np.min(sly4_3) < params[3] < np.max(sly4_3) and \
# np.min(sly4_4) < params[4] < np.max(sly4_4) and np.min(sly4_5) < params[5] < np.max(sly4_5) and \
# np.min(sly4_6) < params[6] < np.max(sly4_6) and np.min(sly4_7) < params[7] < np.max(sly4_7) and \
# np.min(sly4_8) < params[8] < np.max(sly4_8) and np.min(sly4_9) < params[9] < np.max(sly4_9):
return 0.5
return -np.inf
# As long as the prior is finite,
def ln_posterior(params):
ln_prior_val = ln_prior(params)
if not np.isfinite(ln_prior_val) or not np.isfinite(ln_likelihood(params)[0]):
return -np.inf, ln_likelihood(params)[1]
return ln_prior_val + ln_likelihood(params)[0], ln_likelihood(params)[1]
ndim = 3
nwalkers = 32
nsteps = 1000
initial_positions = [list_param_sly4[5:8] + 1e-2 * np.random.randn(ndim) * param_ranges for i in range(nwalkers)]
# Setup the backend to ensure data is consistently logged
# Good backup to have in case the code crashes halfway through for example
filename = "gaussian_likelihood_sym_params.h5"
backend = emcee.backends.HDFBackend(filename)
backend.reset(nwalkers, ndim)
#with Pool() as pool:
sampler = emcee.EnsembleSampler(nwalkers, ndim, ln_posterior, backend=backend) #, pool=pool)
sampler.run_mcmc(initial_positions, nsteps, **{'skip_initial_state_check':True})