Source code for neost.Likelihood

import numpy as np
from scipy.interpolate import UnivariateSpline
import warnings

from neost.Star import Star
from neost import global_imports
from neost.utils import m1_from_mc_m2

c = global_imports._c
G = global_imports._G
Msun = global_imports._M_s
pi = global_imports._pi
rho_ns = global_imports._rhons
n_ns = global_imports._n_ns


[docs] class Likelihood(): def __init__(self, prior, likelihood_functions, likelihood_params, chirp_masses): self.prior = prior self.likelihood_functions = likelihood_functions self.likelihood_params = likelihood_params self.chirp_masses = chirp_masses
[docs] def call(self, pr): likelihoods = [] pr_dict = self.prior.pr constraints = self.prior.EOS.check_constraints() if constraints is False: return -1e101 for i in range(self.prior.number_stars): if self.prior.EOS.adm_type == 'Bosonic': epsdm_cent = self.prior.EOS.find_epsdm_cent(ADM_fraction=pr_dict['adm_fraction'], epscent = 10**(pr_dict['rhoc_' + str(i + 1)])) if (self.prior.EOS.reach_fraction == False): return -1e101 elif (self.prior.EOS.reach_fraction == True): star = Star(10**(pr_dict['rhoc_' + str(i + 1)]), epsdm_cent) star.solve_structure(self.prior.EOS.energydensities, self.prior.EOS.pressures, self.prior.EOS.energydensities_dm, self.prior.EOS.pressures_dm, self.prior.EOS.dm_halo, self.prior.EOS.two_fluid_tidal) Mgrav = star.Mrot Req = star.Req Rdm_halo = star.Rdm_halo tidal = star.tidal #print('adm: ', Mgrav, Req, Rdm_halo,np.log10(pr_dict['mchi']), np.log10(pr_dict['gchi_over_mphi']), pr_dict['adm_fraction']) if self.prior.EOS.adm_type == 'Fermionic': #Hard cut-off imposed as all stars within these boxes have masses well below 1 Msun [~0.4 Msun down to ~0.001 Msun], thus this will save computation time if the code doesn't even have to compute them. if (pr_dict['mchi'] >= pow(10,6) and pr_dict['gchi_over_mphi'] <= pow(10,-3.5) and pr_dict['adm_fraction'] >= 0.01): return -1e101 elif (pr_dict['mchi'] >= pow(10,7) and pr_dict['gchi_over_mphi'] <= pow(10,-3) and pr_dict['adm_fraction'] >= 0.01): return -1e101 elif (pr_dict['mchi'] >= pow(10,8) and pr_dict['gchi_over_mphi'] <= pow(10,-2) and pr_dict['adm_fraction'] >= 0.01): return -1e101 elif (pr_dict['mchi'] >= pow(10,8.5) and pr_dict['gchi_over_mphi'] <= pow(10,-1) and pr_dict['adm_fraction'] >= 0.01): return -1e101 else: epsdm_cent = self.prior.EOS.find_epsdm_cent(ADM_fraction=pr_dict['adm_fraction'],epscent = 10**(pr_dict['rhoc_' + str(i + 1)])) if (self.prior.EOS.reach_fraction == False): return -1e101 elif (self.prior.EOS.reach_fraction == True): star = Star(10**(pr_dict['rhoc_' + str(i + 1)]), epsdm_cent) star.solve_structure(self.prior.EOS.energydensities, self.prior.EOS.pressures, self.prior.EOS.energydensities_dm, self.prior.EOS.pressures_dm, self.prior.EOS.dm_halo, self.prior.EOS.two_fluid_tidal) Mgrav = star.Mrot Req = star.Req Rdm_halo = star.Rdm_halo tidal = star.tidal if self.prior.EOS.adm_type == 'None': star = Star(10**(pr_dict['rhoc_' + str(i + 1)]), 0.0) star.solve_structure(self.prior.EOS.energydensities,self.prior.EOS.pressures) Mgrav = star.Mrot Req = star.Req tidal = star.tidal Rdm_halo = 0.0 if (Mgrav > 3. or Mgrav < 1. or Req > 16. or Req < Mgrav / Msun * 2. * G / c**2. / 1e5): return -1e101 if self.prior.EOS.dm_halo == False and Rdm_halo > 0: return -1e101 if self.chirp_masses[i] is None: MassRadiusTidal = {'Mass':Mgrav, 'Radius':Req, 'Lambda':tidal} tmp = list(map(MassRadiusTidal.get, self.likelihood_params[i])) like = self.likelihood_functions[i](tmp) likelihoods.append(self.array_to_scalar(like)) else: if self.prior.EOS.adm_type != 'None': warnings.warn("Two-fluid tidal deformability only implemented in python. Performance will slow!!") M2 = Mgrav M1 = m1_from_mc_m2(self.chirp_masses[i], M2) if M1 < M2 or M1 > self.prior.MRT[:,0][-1]: return -1e101 MTspline = UnivariateSpline(self.prior.MRT[:,0], self.prior.MRT[:,2], k=1, s=0) point = np.array([self.chirp_masses[i], M2 / M1, MTspline(M1), tidal]) like = self.likelihood_functions[i](point) likelihoods.append(self.array_to_scalar(like)) like_total = np.prod(np.array(likelihoods)) # print('lnlike is', like_total) if like_total == 0.0: return -1e101 return np.log(like_total)
[docs] def array_to_scalar(self, var): # np no longer accepts ragged arrays, # check that if var is an array it only has one element and return it if hasattr(var, '__len__'): try: assert(len(var) == 1) except AssertionError: raise ValueError('Malformed likelihood') var = var[0] return var
[docs] def loglike_prior(self,pr): pr_dict = self.prior.pr constraints = self.prior.EOS.check_constraints() if constraints is False: return -1e101 if self.prior.EOS.adm_type == 'None': star = Star(self.prior.EOS.max_edsc) star.solve_structure(self.prior.EOS.energydensities, self.prior.EOS.pressures) if(star.Mrot < 1): return -1e101 for i in range(self.prior.number_stars): star = Star(10**(pr_dict['rhoc_' + str(i + 1)]), 0.0) star.solve_structure(self.prior.EOS.energydensities, self.prior.EOS.pressures) if(star.Mrot < 1.): return -1e101 if self.prior.EOS.adm_type == 'Bosonic' or self.prior.EOS.adm_type == 'Fermionic': if (pr_dict['mchi'] >= pow(10,6) and pr_dict['gchi_over_mphi'] <= pow(10,-3.5) and pr_dict['adm_fraction'] >= 0.01): return -1e101 elif (pr_dict['mchi'] >= pow(10,7) and pr_dict['gchi_over_mphi'] <= pow(10,-3) and pr_dict['adm_fraction'] >= 0.01): return -1e101 elif (pr_dict['mchi'] >= pow(10,8) and pr_dict['gchi_over_mphi'] <= pow(10,-2) and pr_dict['adm_fraction'] >= 0.01): return -1e101 elif (pr_dict['mchi'] >= pow(10,8.5) and pr_dict['gchi_over_mphi'] <= pow(10,-1) and pr_dict['adm_fraction'] >= 0.01): return -1e101 else: epsdm_cent = self.prior.EOS.find_epsdm_cent(ADM_fraction=pr_dict['adm_fraction'], epscent = self.prior.EOS.max_edsc) if self.prior.EOS.reach_fraction == False: return -1e101 star = Star(self.prior.EOS.max_edsc,epsdm_cent) star.solve_structure(self.prior.EOS.energydensities, self.prior.EOS.pressures, self.prior.EOS.energydensities_dm, self.prior.EOS.pressures_dm, self.prior.EOS.dm_halo, self.prior.EOS.two_fluid_tidal) if(star.Mrot < 1.): return -1e101 if self.prior.EOS.dm_halo == False and star.Rdm_halo > 0: return -1e101 for i in range(self.prior.number_stars): epsdm_cent = self.prior.EOS.find_epsdm_cent(ADM_fraction=pr_dict['adm_fraction'], epscent = 10**(pr_dict['rhoc_' + str(i + 1)])) if self.prior.EOS.reach_fraction == False: return -1e101 star = Star(10**(pr_dict['rhoc_' + str(i + 1)]), epsdm_cent) star.solve_structure(self.prior.EOS.energydensities, self.prior.EOS.pressures, self.prior.EOS.energydensities_dm, self.prior.EOS.pressures_dm, self.prior.EOS.dm_halo, self.prior.EOS.two_fluid_tidal) if(star.Mrot < 1.): return -1e101 if self.prior.EOS.dm_halo == False and star.Rdm_halo > 0: return -1e101 loglike_const = 1. return loglike_const