diff --git a/bilby/core/prior/analytical.py b/bilby/core/prior/analytical.py
index 9443fc4c5a1a9fedbf9ddf4bd3688e8e1821932d..dd503c98911ac4abfc8264d0a710116b35b275ce 100644
--- a/bilby/core/prior/analytical.py
+++ b/bilby/core/prior/analytical.py
@@ -1517,3 +1517,107 @@ class Categorical(Prior):
idxs = np.isin(val, self.categories)
probs[idxs] = self.lnp
return probs
+
+
+class Triangular(Prior):
+ """
+ Define a new prior class which draws from a triangular distribution.
+
+ For distribution details see: wikipedia.org/wiki/Triangular_distribution
+
+ Here, minimum <= mode <= maximum,
+ where the mode has the highest pdf value.
+
+ """
+ def __init__(self, mode, minimum, maximum, name=None, latex_label=None, unit=None):
+ super(Triangular, self).__init__(
+ name=name,
+ latex_label=latex_label,
+ unit=unit,
+ minimum=minimum,
+ maximum=maximum,
+ )
+ self.mode = mode
+ self.fractional_mode = (self.mode - self.minimum) / (
+ self.maximum - self.minimum
+ )
+ self.scale = self.maximum - self.minimum
+ self.rescaled_minimum = self.minimum - (self.minimum == self.mode) * self.scale
+ self.rescaled_maximum = self.maximum + (self.maximum == self.mode) * self.scale
+
+ def rescale(self, val):
+ """
+ 'Rescale' a sample from standard uniform to a triangular distribution.
+
+ This maps to the inverse CDF. This has been analytically solved for this case.
+
+ Parameters
+ ==========
+ val: Union[float, int, array_like]
+ Uniform probability
+
+ Returns
+ =======
+ Union[float, array_like]: Rescaled probability
+
+ """
+ below_mode = (val * self.scale * (self.mode - self.minimum)) ** 0.5
+ above_mode = ((1 - val) * self.scale * (self.maximum - self.mode)) ** 0.5
+ return (self.minimum + below_mode) * (val < self.fractional_mode) + (
+ self.maximum - above_mode
+ ) * (val >= self.fractional_mode)
+
+ def prob(self, val):
+ """
+ Return the prior probability of val
+
+ Parameters
+ ==========
+ val: Union[float, int, array_like]
+
+ Returns
+ =======
+ float: Prior probability of val
+
+ """
+ between_minimum_and_mode = (
+ (val < self.mode)
+ * (self.minimum <= val)
+ * (val - self.rescaled_minimum)
+ / (self.mode - self.rescaled_minimum)
+ )
+ between_mode_and_maximum = (
+ (self.mode <= val)
+ * (val <= self.maximum)
+ * (self.rescaled_maximum - val)
+ / (self.rescaled_maximum - self.mode)
+ )
+ return 2.0 * (between_minimum_and_mode + between_mode_and_maximum) / self.scale
+
+ def cdf(self, val):
+ """
+ Return the prior cumulative probability at val
+
+ Parameters
+ ==========
+ val: Union[float, int, array_like]
+
+ Returns
+ =======
+ float: prior cumulative probability at val
+
+ """
+ return (
+ + (val > self.mode)
+ + (val > self.minimum)
+ * (val <= self.maximum)
+ / (self.scale)
+ * (
+ (val > self.mode)
+ * (self.rescaled_maximum - val) ** 2.0
+ / (self.mode - self.rescaled_maximum)
+ + (val <= self.mode)
+ * (val - self.rescaled_minimum) ** 2.0
+ / (self.mode - self.rescaled_minimum)
+ )
+ )
diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py
index 2a06a228412235e7ce6ebd270f868ee1cb923fb5..091727c5d1f6fcda39449788ce31c86265371f7a 100644
--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
@@ -12,11 +12,25 @@ import numpy as np
from pandas import DataFrame, Series
from scipy.stats import norm
+from .utils import (lalsim_SimNeutronStarEOS4ParamSDGammaCheck,
+ lalsim_SimNeutronStarEOS4ParameterSpectralDecomposition,
+ lalsim_SimNeutronStarEOS4ParamSDViableFamilyCheck,
+ lalsim_SimNeutronStarEOS3PieceDynamicPolytrope,
+ lalsim_SimNeutronStarEOS3PieceCausalAnalytic,
+ lalsim_SimNeutronStarEOS3PDViableFamilyCheck,
+ lalsim_CreateSimNeutronStarFamily,
+ lalsim_SimNeutronStarEOSMaxPseudoEnthalpy,
+ lalsim_SimNeutronStarEOSSpeedOfSoundGeometerized,
+ lalsim_SimNeutronStarFamMinimumMass,
+ lalsim_SimNeutronStarMaximumMass,
+ lalsim_SimNeutronStarRadius,
+ lalsim_SimNeutronStarLoveNumberK2)
+
from ..core.likelihood import MarginalizedLikelihoodReconstructionError
-from ..core.utils import logger, solar_mass, command_line_args, safe_file_dump
+from ..core.utils import logger, solar_mass, gravitational_constant, speed_of_light, command_line_args, safe_file_dump
from ..core.prior import DeltaFunction
from .utils import lalsim_SimInspiralTransformPrecessingNewInitialConditions
-from .eos.eos import SpectralDecompositionEOS, EOSFamily, IntegrateTOV
+from .eos.eos import IntegrateTOV
from .cosmology import get_cosmology, z_at_value
@@ -261,9 +275,7 @@ def convert_to_lal_binary_black_hole_parameters(parameters):
converted_parameters["delta_phase"]
- np.sign(np.cos(converted_parameters["theta_jn"]))
* converted_parameters["psi"],
- 2 * np.pi
- )
-
+ 2 * np.pi)
added_keys = [key for key in converted_parameters.keys()
if key not in original_keys]
@@ -280,9 +292,11 @@ def convert_to_lal_binary_neutron_star_parameters(parameters):
Mass: mass_1, mass_2
Spin: a_1, a_2, tilt_1, tilt_2, phi_12, phi_jl
Extrinsic: luminosity_distance, theta_jn, phase, ra, dec, geocent_time, psi
+ Tidal: lambda_1, lamda_2, lambda_tilde, delta_lambda_tilde
This involves popping a lot of things from parameters.
The keys in added_keys should be popped after evaluating the waveform.
+ For details on tidal parameters see https://arxiv.org/pdf/1402.5156.pdf.
Parameters
==========
@@ -302,8 +316,9 @@ def convert_to_lal_binary_neutron_star_parameters(parameters):
convert_to_lal_binary_black_hole_parameters(converted_parameters)
if not any([key in converted_parameters for key in
- ['lambda_1', 'lambda_2', 'lambda_tilde', 'delta_lambda_tilde',
- 'eos_spectral_gamma_0', 'lambda_symmetric']]):
+ ['lambda_1', 'lambda_2',
+ 'lambda_tilde', 'delta_lambda_tilde', 'lambda_symmetric',
+ 'eos_polytrope_gamma_0', 'eos_spectral_pca_gamma_0', 'eos_v1']]):
converted_parameters['lambda_1'] = 0
converted_parameters['lambda_2'] = 0
added_keys = added_keys + ['lambda_1', 'lambda_2']
@@ -330,31 +345,215 @@ def convert_to_lal_binary_neutron_star_parameters(parameters):
converted_parameters['lambda_1']\
* converted_parameters['mass_1']**5\
/ converted_parameters['mass_2']**5
- elif 'eos_spectral_gamma_0' in converted_parameters.keys(): # FIXME: This is a clunky way to do this
- # Pick out the eos parameters from dict of parameters and sort them
- eos_parameter_keys = sorted([key for key in original_keys if 'eos_spectral_gamma_' in key])
- gammas = [converted_parameters[key] for key in eos_parameter_keys]
-
- eos = SpectralDecompositionEOS(gammas, sampling_flag=True, e0=1.2856e14, p0=5.3716e32)
- if eos.warning_flag:
- converted_parameters['lambda_1'] = 0.0
- converted_parameters['lambda_2'] = 0.0
- converted_parameters['eos_check'] = False
- elif eos_family_physical_check(eos):
- converted_parameters['lambda_1'] = 0.0
- converted_parameters['lambda_2'] = 0.0
- converted_parameters['eos_check'] = False
- else:
- fam = EOSFamily(eos)
-
- if (converted_parameters['mass_1'] <= fam.maximum_mass and
- converted_parameters['mass_2'] <= fam.maximum_mass):
- converted_parameters['lambda_1'] = fam.lambda_from_mass(converted_parameters['mass_1'])
- converted_parameters['lambda_2'] = fam.lambda_from_mass(converted_parameters['mass_2'])
- else:
- converted_parameters['lambda_1'] = 0.0
- converted_parameters['lambda_2'] = 0.0
- converted_parameters['eos_check'] = False
+ elif 'eos_spectral_pca_gamma_0' in converted_parameters.keys(): # FIXME: This is a clunky way to do this
+ converted_parameters = generate_source_frame_parameters(converted_parameters)
+ float_eos_params = {}
+ max_len = 1
+ eos_keys = ['eos_spectral_pca_gamma_0',
+ 'eos_spectral_pca_gamma_1',
+ 'eos_spectral_pca_gamma_2',
+ 'eos_spectral_pca_gamma_3',
+ 'mass_1_source', 'mass_2_source']
+ for key in eos_keys:
+ try:
+ if (len(converted_parameters[key]) > max_len):
+ max_len = len(converted_parameters[key])
+ except TypeError:
+ float_eos_params[key] = converted_parameters[key]
+ if len(float_eos_params) == len(eos_keys): # case where all eos params are floats (pinned)
+ g0, g1, g2, g3 = spectral_pca_to_spectral(
+ converted_parameters['eos_spectral_pca_gamma_0'],
+ converted_parameters['eos_spectral_pca_gamma_1'],
+ converted_parameters['eos_spectral_pca_gamma_2'],
+ converted_parameters['eos_spectral_pca_gamma_3'])
+ converted_parameters['lambda_1'], converted_parameters['lambda_2'], converted_parameters['eos_check'] = \
+ spectral_params_to_lambda_1_lambda_2(
+ g0, g1, g2, g3, converted_parameters['mass_1_source'], converted_parameters['mass_2_source'])
+ elif len(float_eos_params) < len(eos_keys): # case where some or none of the eos params are floats (pinned)
+ for key in float_eos_params.keys():
+ converted_parameters[key] = np.ones(max_len) * converted_parameters[key]
+ g0pca = converted_parameters['eos_spectral_pca_gamma_0']
+ g1pca = converted_parameters['eos_spectral_pca_gamma_1']
+ g2pca = converted_parameters['eos_spectral_pca_gamma_2']
+ g3pca = converted_parameters['eos_spectral_pca_gamma_3']
+ m1s = converted_parameters['mass_1_source']
+ m2s = converted_parameters['mass_2_source']
+ all_lambda_1 = np.empty(0)
+ all_lambda_2 = np.empty(0)
+ all_eos_check = np.empty(0, dtype=bool)
+ for (g_0pca, g_1pca, g_2pca, g_3pca, m1_s, m2_s) in zip(g0pca, g1pca, g2pca, g3pca, m1s, m2s):
+ g_0, g_1, g_2, g_3 = spectral_pca_to_spectral(g_0pca, g_1pca, g_2pca, g_3pca)
+ lambda_1, lambda_2, eos_check = \
+ spectral_params_to_lambda_1_lambda_2(g_0, g_1, g_2, g_3, m1_s, m2_s)
+ all_lambda_1 = np.append(all_lambda_1, lambda_1)
+ all_lambda_2 = np.append(all_lambda_2, lambda_2)
+ all_eos_check = np.append(all_eos_check, eos_check)
+ converted_parameters['lambda_1'] = all_lambda_1
+ converted_parameters['lambda_2'] = all_lambda_2
+ converted_parameters['eos_check'] = all_eos_check
+ for key in float_eos_params.keys():
+ converted_parameters[key] = float_eos_params[key]
+ elif 'eos_polytrope_gamma_0' and 'eos_polytrope_log10_pressure_1' in converted_parameters.keys():
+ converted_parameters = generate_source_frame_parameters(converted_parameters)
+ float_eos_params = {}
+ max_len = 1
+ eos_keys = ['eos_polytrope_gamma_0',
+ 'eos_polytrope_gamma_1',
+ 'eos_polytrope_gamma_2',
+ 'eos_polytrope_log10_pressure_1',
+ 'eos_polytrope_log10_pressure_2',
+ 'mass_1_source', 'mass_2_source']
+ for key in eos_keys:
+ try:
+ if (len(converted_parameters[key]) > max_len):
+ max_len = len(converted_parameters[key])
+ except TypeError:
+ float_eos_params[key] = converted_parameters[key]
+ if len(float_eos_params) == len(eos_keys): # case where all eos params are floats (pinned)
+ converted_parameters['lambda_1'], converted_parameters['lambda_2'], converted_parameters['eos_check'] = \
+ polytrope_or_causal_params_to_lambda_1_lambda_2(
+ converted_parameters['eos_polytrope_gamma_0'],
+ converted_parameters['eos_polytrope_log10_pressure_1'],
+ converted_parameters['eos_polytrope_gamma_1'],
+ converted_parameters['eos_polytrope_log10_pressure_2'],
+ converted_parameters['eos_polytrope_gamma_2'],
+ converted_parameters['mass_1_source'],
+ converted_parameters['mass_2_source'],
+ causal=0)
+ elif len(float_eos_params) < len(eos_keys): # case where some or none are floats (pinned)
+ for key in float_eos_params.keys():
+ converted_parameters[key] = np.ones(max_len) * converted_parameters[key]
+ pg0 = converted_parameters['eos_polytrope_gamma_0']
+ pg1 = converted_parameters['eos_polytrope_gamma_1']
+ pg2 = converted_parameters['eos_polytrope_gamma_2']
+ logp1 = converted_parameters['eos_polytrope_log10_pressure_1']
+ logp2 = converted_parameters['eos_polytrope_log10_pressure_2']
+ m1s = converted_parameters['mass_1_source']
+ m2s = converted_parameters['mass_2_source']
+ all_lambda_1 = np.empty(0)
+ all_lambda_2 = np.empty(0)
+ all_eos_check = np.empty(0, dtype=bool)
+ for (pg_0, pg_1, pg_2, logp_1, logp_2, m1_s, m2_s) in zip(pg0, pg1, pg2, logp1, logp2, m1s, m2s):
+ lambda_1, lambda_2, eos_check = \
+ polytrope_or_causal_params_to_lambda_1_lambda_2(
+ pg_0, logp_1, pg_1, logp_2, pg_2, m1_s, m2_s, causal=0)
+ all_lambda_1 = np.append(all_lambda_1, lambda_1)
+ all_lambda_2 = np.append(all_lambda_2, lambda_2)
+ all_eos_check = np.append(all_eos_check, eos_check)
+ converted_parameters['lambda_1'] = all_lambda_1
+ converted_parameters['lambda_2'] = all_lambda_2
+ converted_parameters['eos_check'] = all_eos_check
+ for key in float_eos_params.keys():
+ converted_parameters[key] = float_eos_params[key]
+ elif 'eos_polytrope_gamma_0' and 'eos_polytrope_scaled_pressure_ratio' in converted_parameters.keys():
+ converted_parameters = generate_source_frame_parameters(converted_parameters)
+ float_eos_params = {}
+ max_len = 1
+ eos_keys = ['eos_polytrope_gamma_0',
+ 'eos_polytrope_gamma_1',
+ 'eos_polytrope_gamma_2',
+ 'eos_polytrope_scaled_pressure_ratio',
+ 'eos_polytrope_scaled_pressure_2',
+ 'mass_1_source', 'mass_2_source']
+ for key in eos_keys:
+ try:
+ if (len(converted_parameters[key]) > max_len):
+ max_len = len(converted_parameters[key])
+ except TypeError:
+ float_eos_params[key] = converted_parameters[key]
+ if len(float_eos_params) == len(eos_keys): # case where all eos params are floats (pinned)
+ logp1, logp2 = log_pressure_reparameterization_conversion(
+ converted_parameters['eos_polytrope_scaled_pressure_ratio'],
+ converted_parameters['eos_polytrope_scaled_pressure_2'])
+ converted_parameters['lambda_1'], converted_parameters['lambda_2'], converted_parameters['eos_check'] = \
+ polytrope_or_causal_params_to_lambda_1_lambda_2(
+ converted_parameters['eos_polytrope_gamma_0'],
+ logp1,
+ converted_parameters['eos_polytrope_gamma_1'],
+ logp2,
+ converted_parameters['eos_polytrope_gamma_2'],
+ converted_parameters['mass_1_source'],
+ converted_parameters['mass_2_source'],
+ causal=0)
+ elif len(float_eos_params) < len(eos_keys): # case where some or none are floats (pinned)
+ for key in float_eos_params.keys():
+ converted_parameters[key] = np.ones(max_len) * converted_parameters[key]
+ pg0 = converted_parameters['eos_polytrope_gamma_0']
+ pg1 = converted_parameters['eos_polytrope_gamma_1']
+ pg2 = converted_parameters['eos_polytrope_gamma_2']
+ scaledratio = converted_parameters['eos_polytrope_scaled_pressure_ratio']
+ scaled_p2 = converted_parameters['eos_polytrope_scaled_pressure_2']
+ m1s = converted_parameters['mass_1_source']
+ m2s = converted_parameters['mass_2_source']
+ all_lambda_1 = np.empty(0)
+ all_lambda_2 = np.empty(0)
+ all_eos_check = np.empty(0, dtype=bool)
+ for (pg_0, pg_1, pg_2, scaled_ratio, scaled_p_2, m1_s,
+ m2_s) in zip(pg0, pg1, pg2, scaledratio, scaled_p2, m1s, m2s):
+ logp_1, logp_2 = log_pressure_reparameterization_conversion(scaled_ratio, scaled_p_2)
+ lambda_1, lambda_2, eos_check = \
+ polytrope_or_causal_params_to_lambda_1_lambda_2(
+ pg_0, logp_1, pg_1, logp_2, pg_2, m1_s, m2_s, causal=0)
+ all_lambda_1 = np.append(all_lambda_1, lambda_1)
+ all_lambda_2 = np.append(all_lambda_2, lambda_2)
+ all_eos_check = np.append(all_eos_check, eos_check)
+ converted_parameters['lambda_1'] = all_lambda_1
+ converted_parameters['lambda_2'] = all_lambda_2
+ converted_parameters['eos_check'] = all_eos_check
+ for key in float_eos_params.keys():
+ converted_parameters[key] = float_eos_params[key]
+ elif 'eos_v1' in converted_parameters.keys():
+ converted_parameters = generate_source_frame_parameters(converted_parameters)
+ float_eos_params = {}
+ max_len = 1
+ eos_keys = ['eos_v1',
+ 'eos_v2',
+ 'eos_v3',
+ 'eos_log10_pressure1_cgs',
+ 'eos_log10_pressure2_cgs',
+ 'mass_1_source', 'mass_2_source']
+ for key in eos_keys:
+ try:
+ if (len(converted_parameters[key]) > max_len):
+ max_len = len(converted_parameters[key])
+ except TypeError:
+ float_eos_params[key] = converted_parameters[key]
+ if len(float_eos_params) == len(eos_keys): # case where all eos params are floats (pinned)
+ converted_parameters['lambda_1'], converted_parameters['lambda_2'], converted_parameters['eos_check'] = \
+ polytrope_or_causal_params_to_lambda_1_lambda_2(
+ converted_parameters['eos_v1'],
+ converted_parameters['eos_log10_pressure1_cgs'],
+ converted_parameters['eos_v2'],
+ converted_parameters['eos_log10_pressure2_cgs'],
+ converted_parameters['eos_v3'],
+ converted_parameters['mass_1_source'],
+ converted_parameters['mass_2_source'],
+ causal=1)
+ elif len(float_eos_params) < len(eos_keys): # case where some or none are floats (pinned)
+ for key in float_eos_params.keys():
+ converted_parameters[key] = np.ones(max_len) * converted_parameters[key]
+ v1 = converted_parameters['eos_v1']
+ v2 = converted_parameters['eos_v2']
+ v3 = converted_parameters['eos_v3']
+ logp1 = converted_parameters['eos_log10_pressure1_cgs']
+ logp2 = converted_parameters['eos_log10_pressure2_cgs']
+ m1s = converted_parameters['mass_1_source']
+ m2s = converted_parameters['mass_2_source']
+ all_lambda_1 = np.empty(0)
+ all_lambda_2 = np.empty(0)
+ all_eos_check = np.empty(0, dtype=bool)
+ for (v_1, v_2, v_3, logp_1, logp_2, m1_s, m2_s) in zip(v1, v2, v3, logp1, logp2, m1s, m2s):
+ lambda_1, lambda_2, eos_check = \
+ polytrope_or_causal_params_to_lambda_1_lambda_2(
+ v_1, logp_1, v_2, logp_2, v_3, m1_s, m2_s, causal=1)
+ all_lambda_1 = np.append(all_lambda_1, lambda_1)
+ all_lambda_2 = np.append(all_lambda_2, lambda_2)
+ all_eos_check = np.append(all_eos_check, eos_check)
+ converted_parameters['lambda_1'] = all_lambda_1
+ converted_parameters['lambda_2'] = all_lambda_2
+ converted_parameters['eos_check'] = all_eos_check
+ for key in float_eos_params.keys():
+ converted_parameters[key] = float_eos_params[key]
elif 'lambda_symmetric' in converted_parameters.keys():
if 'lambda_antisymmetric' in converted_parameters.keys():
converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\
@@ -380,6 +579,233 @@ def convert_to_lal_binary_neutron_star_parameters(parameters):
return converted_parameters, added_keys
+def log_pressure_reparameterization_conversion(scaled_pressure_ratio, scaled_pressure_2):
+ '''
+ Converts the reparameterization joining pressures from
+ (scaled_pressure_ratio,scaled_pressure_2) to (log10_pressure_1,log10_pressure_2).
+ This reparameterization with a triangular prior (with mode = max) on scaled_pressure_2
+ and a uniform prior on scaled_pressure_ratio
+ mimics identical uniform priors on log10_pressure_1 and log10_pressure_2
+ where samples with log10_pressure_2 > log10_pressure_1 are rejected.
+ This reparameterization allows for a faster initialization.
+ A minimum log10_pressure of 33 (in cgs units) is chosen to be slightly higher than the low-density crust EOS
+ that is stitched to the dynamic polytrope EOS model in LALSimulation.
+
+ Parameters
+ ----------
+ scaled_pressure_ratio, scaled_pressure_2: float
+ reparameterizations of the dividing pressures
+
+ Returns
+ -------
+ log10_pressure_1, log10_pressure_2: float
+ joining pressures in the original parameterization
+
+ '''
+ minimum_pressure = 33.
+ log10_pressure_1 = (scaled_pressure_ratio * scaled_pressure_2) + minimum_pressure
+ log10_pressure_2 = minimum_pressure + scaled_pressure_2
+
+ return log10_pressure_1, log10_pressure_2
+
+
+def spectral_pca_to_spectral(gamma_pca_0, gamma_pca_1, gamma_pca_2, gamma_pca_3):
+ '''
+ Change of basis on parameter space
+ from an efficient space to sample in (sample space)
+ to the space used in spectral eos model (model space).
+ Uses principle component analysis (PCA) to sample spectral gamma parameters more efficiently.
+ See arxiv:2001.01747 for the PCA conversion transformation matrix, mean, and standard deviation.
+ (Note that this transformation matrix is the inverse of the one referenced
+ in order to perform the inverse transformation.)
+
+ Parameters
+ ----------
+ gamma_pca_0, gamma_pca_1, gamma_pca_2, gamma_pca_3: float
+ Spectral gamma PCA parameters to be converted to spectral gamma parameters.
+
+ Returns
+ -------
+ converted_gamma_parameters: np.array()
+ array of gamma_0, gamma_1, gamma_2, gamma_3 in model space
+
+ '''
+ sampled_pca_gammas = np.array([gamma_pca_0, gamma_pca_1, gamma_pca_2, gamma_pca_3])
+ transformation_matrix = np.array(
+ [
+ [0.43801, -0.76705, 0.45143, 0.12646],
+ [-0.53573, 0.17169, 0.67968, 0.47070],
+ [0.52660, 0.31255, -0.19454, 0.76626],
+ [-0.49379, -0.53336, -0.54444, 0.41868]
+ ]
+ )
+
+ model_space_mean = np.array([0.89421, 0.33878, -0.07894, 0.00393])
+ model_space_standard_deviation = np.array([0.35700, 0.25769, 0.05452, 0.00312])
+ converted_gamma_parameters = \
+ model_space_mean + model_space_standard_deviation * np.dot(transformation_matrix, sampled_pca_gammas)
+
+ return converted_gamma_parameters
+
+
+def spectral_params_to_lambda_1_lambda_2(gamma_0, gamma_1, gamma_2, gamma_3, mass_1_source, mass_2_source):
+ '''
+ Converts from the 4 spectral decomposition parameters and the source masses
+ to the tidal deformability parameters.
+
+ Parameters
+ ----------
+ gamma_0, gamma_1, gamma_2, gamma_3: float
+ sampled spectral decomposition parameters
+ mass_1_source, mass_2_source: float
+ sampled component mass parameters converted to source frame in solar masses
+
+ Returns
+ -------
+ lambda_1, lambda_2: float
+ component tidal deformability parameters
+ eos_check: bool
+ whether or not the equation of state is viable /
+ if eos_check = False, lambdas are 0 and the sample is rejected.
+
+ '''
+ eos_check = True
+ if lalsim_SimNeutronStarEOS4ParamSDGammaCheck(gamma_0, gamma_1, gamma_2, gamma_3) != 0:
+ lambda_1 = 0.0
+ lambda_2 = 0.0
+ eos_check = False
+ else:
+ eos = lalsim_SimNeutronStarEOS4ParameterSpectralDecomposition(gamma_0, gamma_1, gamma_2, gamma_3)
+ if lalsim_SimNeutronStarEOS4ParamSDViableFamilyCheck(gamma_0, gamma_1, gamma_2, gamma_3) != 0:
+ lambda_1 = 0.0
+ lambda_2 = 0.0
+ eos_check = False
+ else:
+ lambda_1, lambda_2, eos_check = neutron_star_family_physical_check(eos, mass_1_source, mass_2_source)
+
+ return lambda_1, lambda_2, eos_check
+
+
+def polytrope_or_causal_params_to_lambda_1_lambda_2(
+ param1, log10_pressure1_cgs, param2, log10_pressure2_cgs, param3, mass_1_source, mass_2_source, causal):
+ """
+ Converts parameters from sampled dynamic piecewise polytrope parameters
+ to component tidal deformablity parameters.
+ Enforces log10_pressure1_cgs < log10_pressure2_cgs.
+ Checks number of points in the equation of state for viability.
+ Note that subtracting 1 from the log10 pressure in cgs converts it to
+ log10 pressure in si units.
+
+ Parameters
+ ----------
+ param1, param2, param3: float
+ either the sampled adiabatic indices in piecewise polytrope model
+ or the sampled causal model params v1, v2, v3
+ log10_pressure1_cgs, log10_pressure2_cgs: float
+ dividing pressures in piecewise polytrope model or causal model
+ mass_1_source, mass_2_source: float
+ source frame component mass parameters in Msuns
+ causal: bool
+ whether or not to use causal polytrope model
+ 1 - causal; 0 - not causal
+
+ Returns
+ -------
+ lambda_1: float
+ tidal deformability parameter associated with mass 1
+ lambda_2: float
+ tidal deformability parameter associated with mass 2
+ eos_check: bool
+ whether eos is valid or not
+
+ """
+ eos_check = True
+ if log10_pressure1_cgs >= log10_pressure2_cgs:
+ lambda_1 = 0.0
+ lambda_2 = 0.0
+ eos_check = False
+ else:
+ if causal == 0:
+ eos = lalsim_SimNeutronStarEOS3PieceDynamicPolytrope(
+ param1, log10_pressure1_cgs - 1., param2, log10_pressure2_cgs - 1., param3)
+ else:
+ eos = lalsim_SimNeutronStarEOS3PieceCausalAnalytic(
+ param1, log10_pressure1_cgs - 1., param2, log10_pressure2_cgs - 1., param3)
+ if lalsim_SimNeutronStarEOS3PDViableFamilyCheck(
+ param1, log10_pressure1_cgs - 1., param2, log10_pressure2_cgs - 1., param3, causal) != 0:
+ lambda_1 = 0.0
+ lambda_2 = 0.0
+ eos_check = False
+ else:
+ lambda_1, lambda_2, eos_check = neutron_star_family_physical_check(eos, mass_1_source, mass_2_source)
+
+ return lambda_1, lambda_2, eos_check
+
+
+def neutron_star_family_physical_check(eos, mass_1_source, mass_2_source):
+ """
+ Takes in a lalsim eos object. Performs causal and max/min mass eos checks.
+ Calculates component lambdas if eos object passes causality.
+ Returns lambda = 0 if not.
+
+ Parameters
+ ----------
+ eos: lalsim swig-wrapped eos object
+ the neutron star equation of state
+ mass_1_source, mass_2_source: float
+ source frame component masses 1 and 2 in solar masses
+
+ Returns
+ -------
+ lambda_1, lambda_2: float
+ component tidal deformability parameters
+ eos_check: bool
+ whether or not the equation of state is physically allowed
+
+ """
+ eos_check = True
+ family = lalsim_CreateSimNeutronStarFamily(eos)
+ max_pseudo_enthalpy = lalsim_SimNeutronStarEOSMaxPseudoEnthalpy(eos)
+ max_speed_of_sound = lalsim_SimNeutronStarEOSSpeedOfSoundGeometerized(max_pseudo_enthalpy, eos)
+ min_mass = lalsim_SimNeutronStarFamMinimumMass(family) / solar_mass
+ max_mass = lalsim_SimNeutronStarMaximumMass(family) / solar_mass
+ if max_speed_of_sound <= 1.1 and min_mass <= mass_1_source <= max_mass and min_mass <= mass_2_source <= max_mass:
+ lambda_1 = lambda_from_mass_and_family(mass_1_source, family)
+ lambda_2 = lambda_from_mass_and_family(mass_2_source, family)
+ else:
+ lambda_1 = 0.0
+ lambda_2 = 0.0
+ eos_check = False
+
+ return lambda_1, lambda_2, eos_check
+
+
+def lambda_from_mass_and_family(mass_i, family):
+ """
+ Convert from equation of state model parameters to
+ component tidal parameters.
+
+ Parameters
+ ----------
+ family: lalsim family object
+ EOS family of type lalsimulation.SimNeutronStarFamily.
+ mass_i: Component mass of neutron star in solar masses.
+
+ Returns
+ -------
+ lambda_1: float
+ component tidal deformability parameter
+
+ """
+ radius = lalsim_SimNeutronStarRadius(mass_i * solar_mass, family)
+ love_number_k2 = lalsim_SimNeutronStarLoveNumberK2(mass_i * solar_mass, family)
+ mass_geometrized = mass_i * solar_mass * gravitational_constant / speed_of_light ** 2.
+ compactness = mass_geometrized / radius
+ lambda_i = (2. / 3.) * love_number_k2 / compactness ** 5.
+
+ return lambda_i
+
+
def eos_family_physical_check(eos):
"""
Function that determines if the EoS family contains
@@ -712,6 +1138,7 @@ def lambda_1_lambda_2_to_lambda_tilde(lambda_1, lambda_2, mass_1, mass_2):
=======
lambda_tilde: float
Dominant tidal term.
+
"""
eta = component_masses_to_symmetric_mass_ratio(mass_1, mass_2)
lambda_plus = lambda_1 + lambda_2
@@ -750,9 +1177,8 @@ def lambda_1_lambda_2_to_delta_lambda_tilde(lambda_1, lambda_2, mass_1, mass_2):
lambda_minus = lambda_1 - lambda_2
delta_lambda_tilde = 1 / 2 * (
(1 - 4 * eta) ** 0.5 * (1 - 13272 / 1319 * eta + 8944 / 1319 * eta**2) *
- lambda_plus + (1 - 15910 / 1319 * eta + 32850 / 1319 * eta**2 +
- 3380 / 1319 * eta**3) * lambda_minus)
-
+ lambda_plus + (1 - 15910 / 1319 * eta + 32850 / 1319 * eta ** 2 +
+ 3380 / 1319 * eta ** 3) * lambda_minus)
return delta_lambda_tilde
@@ -780,6 +1206,7 @@ def lambda_tilde_delta_lambda_tilde_to_lambda_1_lambda_2(
Tidal parameter of more massive neutron star.
lambda_2: float
Tidal parameter of less massive neutron star.
+
"""
eta = component_masses_to_symmetric_mass_ratio(mass_1, mass_2)
coefficient_1 = (1 + 7 * eta - 31 * eta**2)
@@ -798,6 +1225,7 @@ def lambda_tilde_delta_lambda_tilde_to_lambda_1_lambda_2(
2 * delta_lambda_tilde * (coefficient_1 + coefficient_2)) \
/ ((coefficient_1 - coefficient_2) * (coefficient_3 + coefficient_4) -
(coefficient_1 + coefficient_2) * (coefficient_3 - coefficient_4))
+
return lambda_1, lambda_2
diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index 11b249f75d0c42257c6a5860091fb56347558603..a825fa5e8d7764f28d34a0c7e2e93b6523374b01 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -791,6 +791,160 @@ def lalsim_SimInspiralWaveformParamsInsertTidalLambda2(
waveform_dictionary, lambda_2)
+def lalsim_SimNeutronStarEOS4ParamSDGammaCheck(g0, g1, g2, g3):
+ from lalsimulation import SimNeutronStarEOS4ParamSDGammaCheck
+ try:
+ g0 = float(g0)
+ g1 = float(g1)
+ g2 = float(g2)
+ g3 = float(g3)
+ except ValueError:
+ raise ValueError("Unable to convert EOS spectral parameters to floats")
+ except TypeError:
+ raise TypeError("Unable to convert EOS spectral parameters to floats")
+
+ return SimNeutronStarEOS4ParamSDGammaCheck(g0, g1, g2, g3)
+
+
+def lalsim_SimNeutronStarEOS4ParameterSpectralDecomposition(g0, g1, g2, g3):
+ from lalsimulation import SimNeutronStarEOS4ParameterSpectralDecomposition
+ try:
+ g0 = float(g0)
+ g1 = float(g1)
+ g2 = float(g2)
+ g3 = float(g3)
+ except ValueError:
+ raise ValueError("Unable to convert EOS spectral parameters to floats")
+ except TypeError:
+ raise TypeError("Unable to convert EOS spectral parameters to floats")
+
+ return SimNeutronStarEOS4ParameterSpectralDecomposition(g0, g1, g2, g3)
+
+
+def lalsim_SimNeutronStarEOS4ParamSDViableFamilyCheck(g0, g1, g2, g3):
+ from lalsimulation import SimNeutronStarEOS4ParamSDViableFamilyCheck
+ try:
+ g0 = float(g0)
+ g1 = float(g1)
+ g2 = float(g2)
+ g3 = float(g3)
+ except ValueError:
+ raise ValueError("Unable to convert EOS spectral parameters to floats")
+ except TypeError:
+ raise TypeError("Unable to convert EOS spectral parameters to floats")
+
+ return SimNeutronStarEOS4ParamSDViableFamilyCheck(g0, g1, g2, g3)
+
+
+def lalsim_SimNeutronStarEOS3PieceDynamicPolytrope(g0, log10p1_si, g1, log10p2_si, g2):
+ from lalsimulation import SimNeutronStarEOS3PieceDynamicPolytrope
+ try:
+ g0 = float(g0)
+ g1 = float(g1)
+ g2 = float(g2)
+ log10p1_si = float(log10p1_si)
+ log10p2_si = float(log10p2_si)
+ except ValueError:
+ raise ValueError("Unable to convert EOS polytrope parameters to floats")
+ except TypeError:
+ raise TypeError("Unable to convert EOS polytrope parameters to floats")
+
+ return SimNeutronStarEOS3PieceDynamicPolytrope(g0, log10p1_si, g1, log10p2_si, g2)
+
+
+def lalsim_SimNeutronStarEOS3PieceCausalAnalytic(v1, log10p1_si, v2, log10p2_si, v3):
+ from lalsimulation import SimNeutronStarEOS3PieceCausalAnalytic
+ try:
+ v1 = float(v1)
+ v2 = float(v2)
+ v3 = float(v3)
+ log10p1_si = float(log10p1_si)
+ log10p2_si = float(log10p2_si)
+ except ValueError:
+ raise ValueError("Unable to convert EOS causal parameters to floats")
+ except TypeError:
+ raise TypeError("Unable to convert EOS causal parameters to floats")
+
+ return SimNeutronStarEOS3PieceCausalAnalytic(v1, log10p1_si, v2, log10p2_si, v3)
+
+
+def lalsim_SimNeutronStarEOS3PDViableFamilyCheck(p0, log10p1_si, p1, log10p2_si, p2, causal):
+ from lalsimulation import SimNeutronStarEOS3PDViableFamilyCheck
+ try:
+ p0 = float(p0)
+ p1 = float(p1)
+ p2 = float(p2)
+ log10p1_si = float(log10p1_si)
+ log10p2_si = float(log10p2_si)
+ causal = int(causal)
+ except ValueError:
+ raise ValueError("Unable to convert EOS parameters to floats or int")
+ except TypeError:
+ raise TypeError("Unable to convert EOS parameters to floats or int")
+
+ return SimNeutronStarEOS3PDViableFamilyCheck(p0, log10p1_si, p1, log10p2_si, p2, causal)
+
+
+def lalsim_CreateSimNeutronStarFamily(eos):
+ from lalsimulation import CreateSimNeutronStarFamily
+
+ return CreateSimNeutronStarFamily(eos)
+
+
+def lalsim_SimNeutronStarEOSMaxPseudoEnthalpy(eos):
+ from lalsimulation import SimNeutronStarEOSMaxPseudoEnthalpy
+
+ return SimNeutronStarEOSMaxPseudoEnthalpy(eos)
+
+
+def lalsim_SimNeutronStarEOSSpeedOfSoundGeometerized(max_pseudo_enthalpy, eos):
+ from lalsimulation import SimNeutronStarEOSSpeedOfSoundGeometerized
+ try:
+ max_pseudo_enthalpy = float(max_pseudo_enthalpy)
+ except ValueError:
+ raise ValueError("Unable to convert max_pseudo_enthalpy to float.")
+ except TypeError:
+ raise TypeError("Unable to convert max_pseudo_enthalpy to float.")
+
+ return SimNeutronStarEOSSpeedOfSoundGeometerized(max_pseudo_enthalpy, eos)
+
+
+def lalsim_SimNeutronStarFamMinimumMass(fam):
+ from lalsimulation import SimNeutronStarFamMinimumMass
+
+ return SimNeutronStarFamMinimumMass(fam)
+
+
+def lalsim_SimNeutronStarMaximumMass(fam):
+ from lalsimulation import SimNeutronStarMaximumMass
+
+ return SimNeutronStarMaximumMass(fam)
+
+
+def lalsim_SimNeutronStarRadius(mass_in_SI, fam):
+ from lalsimulation import SimNeutronStarRadius
+ try:
+ mass_in_SI = float(mass_in_SI)
+ except ValueError:
+ raise ValueError("Unable to convert mass_in_SI to float.")
+ except TypeError:
+ raise TypeError("Unable to convert mass_in_SI to float.")
+
+ return SimNeutronStarRadius(mass_in_SI, fam)
+
+
+def lalsim_SimNeutronStarLoveNumberK2(mass_in_SI, fam):
+ from lalsimulation import SimNeutronStarLoveNumberK2
+ try:
+ mass_in_SI = float(mass_in_SI)
+ except ValueError:
+ raise ValueError("Unable to convert mass_in_SI to float.")
+ except TypeError:
+ raise TypeError("Unable to convert mass_in_SI to float.")
+
+ return SimNeutronStarLoveNumberK2(mass_in_SI, fam)
+
+
def spline_angle_xform(delta_psi):
"""
Returns the angle in degrees corresponding to the spline
diff --git a/examples/gw_examples/injection_examples/bns_polytrope_eos_example.py b/examples/gw_examples/injection_examples/bns_polytrope_eos_example.py
new file mode 100644
index 0000000000000000000000000000000000000000..e25119d436d20622281c63fb25b16f36ab169f25
--- /dev/null
+++ b/examples/gw_examples/injection_examples/bns_polytrope_eos_example.py
@@ -0,0 +1,171 @@
+#!/usr/bin/env python
+"""
+Tutorial to demonstrate running parameter estimation on a binary neutron star
+with equation of state inference using a 3 piece dynamic polytrope model.
+
+This example estimates the masses and polytrope model parameters.
+
+This polytrope model is generically introduced in https://arxiv.org/pdf/1804.04072.pdf.
+LALSimulation has an implementation of this generic model with 3 polytrope pieces attached
+to a crust EOS, and that can be found in LALSimNeutronStarEOSDynamicPolytrope.c.
+
+Details on a 4 piece static polytrope model can be found here https://arxiv.org/pdf/0812.2163.pdf.
+The reparameterization here uses a 3-piece dynamic polytrope model implemented in lalsimulation
+LALSimNeutronStarEOSDynamicPolytrope.c.
+"""
+
+
+import bilby
+import numpy as np
+
+# Specify the output directory and the name of the simulation.
+outdir = "outdir"
+label = "bns_polytrope_eos_simulation"
+bilby.core.utils.setup_logger(outdir=outdir, label=label)
+# Set up a random seed for result reproducibility. This is optional!
+np.random.seed(88170235)
+
+# We are going to inject a binary neutron star waveform with component mass and lambda
+# values consistent with the MPA1 equation of state.
+# The lambda values are generated from LALSimNeutronStarEOS_MPA1.dat
+# lalsim-ns-params -n MPA1 -m 1.523194
+# lalsim-ns-params -n MPA1 -m 1.5221469
+# Note that we injection with tidal params (lambdas) but recover in eos model params.
+injection_parameters = dict(
+ mass_1=1.523194,
+ mass_2=1.5221469,
+ lambda_1=311.418,
+ lambda_2=312.717,
+ chi_1=0.0001,
+ chi_2=0.0001,
+ luminosity_distance=57.628867,
+ theta_jn=0.66246341,
+ psi=0.18407784,
+ phase=5.4800181,
+ geocent_time=10.0,
+ ra=3.6309322,
+ dec=-0.30355747,
+)
+
+# Set the duration and sampling frequency of the data segment that we're going
+# to inject the signal into. For the
+# TaylorF2 waveform, we cut the signal close to the isco frequency
+duration = 64
+sampling_frequency = 2 * 2048
+start_time = injection_parameters["geocent_time"] + 2 - duration
+
+# Fixed arguments passed into the source model. The analysis starts at 40 Hz.
+waveform_arguments = dict(
+ waveform_approximant="TaylorF2",
+ reference_frequency=30.0,
+ minimum_frequency=30.0,
+)
+
+# Create the waveform_generator using a LAL Binary Neutron Star source function
+waveform_generator = bilby.gw.WaveformGenerator(
+ duration=duration,
+ sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=bilby.gw.source.lal_binary_neutron_star,
+ parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_neutron_star_parameters,
+ waveform_arguments=waveform_arguments,
+)
+
+# Set up interferometers. In this case we'll use three interferometers
+# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)).
+# These default to their design sensitivity and start at 40 Hz.
+interferometers = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
+for interferometer in interferometers:
+ interferometer.minimum_frequency = 30.0
+interferometers.set_strain_data_from_zero_noise(
+ sampling_frequency=sampling_frequency, duration=duration, start_time=start_time
+)
+interferometers.inject_signal(
+ parameters=injection_parameters, waveform_generator=waveform_generator
+)
+
+# Load the default prior for binary neutron stars.
+# We're going to sample in chirp_mass, mass_ratio, and model parameters
+# rather than mass_1, mass_2, or lamba_tilde, delta_lambda_tilde.
+# BNS have aligned spins by default, if you want to allow precessing spins
+# pass aligned_spin=False to the BNSPriorDict
+priors = bilby.gw.prior.BNSPriorDict()
+
+priors["chirp_mass"] = bilby.gw.prior.Uniform(0.4, 4.4)
+priors["mass_ratio"] = bilby.gw.prior.Uniform(0.125, 1)
+
+# The following are dynamic polytrope model priors
+# They are required for EOS inference
+priors["eos_polytrope_gamma_0"] = bilby.core.prior.Uniform(
+ 1.0, 5.0, name="Gamma0", latex_label="$\\Gamma_0$"
+)
+priors["eos_polytrope_gamma_1"] = bilby.core.prior.Uniform(
+ 1.0, 5.0, name="Gamma1", latex_label="$\\Gamma_1$"
+)
+priors["eos_polytrope_gamma_2"] = bilby.core.prior.Uniform(
+ 1.0, 5.0, name="Gamma2", latex_label="$\\Gamma_2$"
+)
+"""
+One can run this model without the reparameterization using the following
+priors in place of the scaled pressure priors. The reparameterization approximates
+the flat priors in pressure however the reparameterization initializes twice as fast.
+
+priors["eos_polytrope_log10_pressure_1"] = bilby.core.prior.Uniform(
+ 34.0, 37.0, name="plogp1", latex_label="$log(p_1)$"
+)
+priors["eos_polytrope_log10_pressure_2"] = bilby.core.prior.Uniform(
+ 34.0, 37.0, name="plogp2", latex_label="$log(p_2)$"
+)
+"""
+priors["eos_polytrope_scaled_pressure_ratio"] = bilby.core.prior.Uniform(
+ minimum=0.0, maximum=1, name="Pressure_Ratio", latex_label="$\\rho$"
+)
+priors["eos_polytrope_scaled_pressure_2"] = bilby.core.prior.Triangular(
+ minimum=0.0, maximum=4.0, mode=4, name="Pressure_offset", latex_label="$p_0$"
+)
+priors["eos_check"] = bilby.gw.prior.EOSCheck()
+# The eos_check sets up a rejection prior for sampling.
+# It should be enabled for all eos runs.
+
+# Pinning other parameters to their injected values.
+for key in [
+ "luminosity_distance",
+ "geocent_time",
+ "chi_1",
+ "chi_2",
+ "psi",
+ "phase",
+ "ra",
+ "dec",
+ "theta_jn",
+]:
+ priors[key] = injection_parameters[key]
+
+# Remove tidal parameters from our priors.
+priors.pop("lambda_1")
+priors.pop("lambda_2")
+
+# Initialise the likelihood by passing in the interferometer data (IFOs)
+# and the waveform generator
+likelihood = bilby.gw.GravitationalWaveTransient(
+ interferometers=interferometers,
+ waveform_generator=waveform_generator,
+ time_marginalization=False,
+ phase_marginalization=False,
+ distance_marginalization=False,
+ priors=priors,
+)
+
+# Run sampler. This model has mostly been testing with the `dynesty` sampler.
+result = bilby.run_sampler(
+ conversion_function=bilby.gw.conversion.generate_all_bns_parameters,
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ nlive=1000,
+ dlogz=0.1,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+)
+
+result.plot_corner()
diff --git a/examples/gw_examples/injection_examples/bns_spectral_pca_eos_example.py b/examples/gw_examples/injection_examples/bns_spectral_pca_eos_example.py
new file mode 100644
index 0000000000000000000000000000000000000000..280c5c2013dcf7c1be58451eed9144157fabc377
--- /dev/null
+++ b/examples/gw_examples/injection_examples/bns_spectral_pca_eos_example.py
@@ -0,0 +1,151 @@
+#!/usr/bin/env python
+"""
+Tutorial to demonstrate running parameter estimation on a binary neutron star
+with equation of state inference using a spectral decomposition model.
+
+This example estimates the masses and spectral model parameters.
+
+Details on the spectral model can be found in https://arxiv.org/pdf/1805.11217.pdf.
+The principal component analysis (PCA) used for the spectral model here can be found
+in the appendix of https://arxiv.org/pdf/2001.01747.pdf.
+"""
+
+
+import bilby
+import numpy as np
+
+# Specify the output directory and the name of the simulation.
+outdir = "outdir"
+label = "bns_spectral_eos_simulation"
+bilby.core.utils.setup_logger(outdir=outdir, label=label)
+# Set up a random seed for result reproducibility. This is optional!
+np.random.seed(88170235)
+
+# We are going to inject a binary neutron star waveform with component mass and lambda
+# values consistent with the MPA1 equation of state.
+# The lambda values are generated from LALSimNeutronStarEOS_MPA1.dat
+# lalsim-ns-params -n MPA1 -m 1.523194
+# lalsim-ns-params -n MPA1 -m 1.5221469
+# Note that we injection with tidal params (lambdas) but recover in eos model params.
+injection_parameters = dict(
+ mass_1=1.523194,
+ mass_2=1.5221469,
+ lambda_1=311.418,
+ lambda_2=312.717,
+ chi_1=0.0001,
+ chi_2=0.0001,
+ luminosity_distance=57.628867,
+ theta_jn=0.66246341,
+ psi=0.18407784,
+ phase=5.4800181,
+ geocent_time=10.0,
+ ra=3.6309322,
+ dec=-0.30355747,
+)
+
+# Set the duration and sampling frequency of the data segment that we're going
+# to inject the signal into. For the
+# TaylorF2 waveform, we cut the signal close to the isco frequency
+duration = 64
+sampling_frequency = 2 * 2048
+start_time = injection_parameters["geocent_time"] + 2 - duration
+
+# Fixed arguments passed into the source model. The analysis starts at 40 Hz.
+waveform_arguments = dict(
+ waveform_approximant="TaylorF2",
+ reference_frequency=30.0,
+ minimum_frequency=30.0,
+)
+
+# Create the waveform_generator using a LAL Binary Neutron Star source function
+waveform_generator = bilby.gw.WaveformGenerator(
+ duration=duration,
+ sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=bilby.gw.source.lal_binary_neutron_star,
+ parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_neutron_star_parameters,
+ waveform_arguments=waveform_arguments,
+)
+
+# Set up interferometers. In this case we'll use three interferometers
+# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)).
+# These default to their design sensitivity and start at 30 Hz.
+interferometers = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
+for interferometer in interferometers:
+ interferometer.minimum_frequency = 30.0
+interferometers.set_strain_data_from_zero_noise(
+ sampling_frequency=sampling_frequency, duration=duration, start_time=start_time
+)
+interferometers.inject_signal(
+ parameters=injection_parameters, waveform_generator=waveform_generator
+)
+
+# Load the default prior for binary neutron stars.
+# We're going to sample in chirp_mass, mass_ratio, and model parameters
+# rather than mass_1, mass_2, or lamba_tilde, delta_lambda_tilde.
+# BNS have aligned spins by default, if you want to allow precessing spins
+# pass aligned_spin=False to the BNSPriorDict
+priors = bilby.gw.prior.BNSPriorDict()
+# The following are spectral decomposition model priors
+# Note that we sample in 'eos_spectral_pca_gamma_*', which is a
+# prior space that enables more efficient. These parameters are
+# subsequently converted to the spectral decomposition model space.
+# These priors are required for EOS inference
+priors["eos_spectral_pca_gamma_0"] = bilby.core.prior.Uniform(
+ -4.37722, 4.91227, name="pca_gamma0", latex_label="$\\gamma^{pca}_0$"
+)
+priors["eos_spectral_pca_gamma_1"] = bilby.core.prior.Uniform(
+ -1.82240, 2.06387, name="pca_gamma1", latex_label="$\\gamma^{pca}_1$"
+)
+priors["eos_spectral_pca_gamma_2"] = bilby.core.prior.Uniform(
+ -0.32445, 0.36469, name="pca_gamma2", latex_label="$\\gamma^{pca}_2$"
+)
+priors["eos_spectral_pca_gamma_3"] = bilby.core.prior.Uniform(
+ -0.09529, 0.11046, name="pca_gamma3", latex_label="$\\gamma^{pca}_3$"
+)
+priors["eos_check"] = bilby.gw.prior.EOSCheck()
+# The eos_check sets up a rejection prior for sampling.
+# It should be enabled for all eos runs.
+
+# Pinning other parameters to their injected values.
+for key in [
+ "luminosity_distance",
+ "geocent_time",
+ "chi_1",
+ "chi_2",
+ "psi",
+ "phase",
+ "ra",
+ "dec",
+ "theta_jn",
+]:
+ priors[key] = injection_parameters[key]
+
+# Remove tidal parameters from our priors.
+priors.pop("lambda_1")
+priors.pop("lambda_2")
+
+# Initialise the likelihood by passing in the interferometer data (IFOs)
+# and the waveform generator
+likelihood = bilby.gw.GravitationalWaveTransient(
+ interferometers=interferometers,
+ waveform_generator=waveform_generator,
+ time_marginalization=False,
+ phase_marginalization=False,
+ distance_marginalization=False,
+ priors=priors,
+)
+
+# Run sampler. This model has mostly been testing with the `dynesty` sampler.
+result = bilby.run_sampler(
+ conversion_function=bilby.gw.conversion.generate_all_bns_parameters,
+ likelihood=likelihood,
+ priors=priors,
+ sampler="dynesty",
+ nlive=1000,
+ dlogz=0.1,
+ injection_parameters=injection_parameters,
+ outdir=outdir,
+ label=label,
+)
+
+result.plot_corner()
diff --git a/test/core/prior/prior_test.py b/test/core/prior/prior_test.py
index 3fd77af938a4abf8fad832a669e89f797ea87f5f..364d48e9e0302ba8d97ba92de1643d77f0fc9e4c 100644
--- a/test/core/prior/prior_test.py
+++ b/test/core/prior/prior_test.py
@@ -135,6 +135,27 @@ class TestPriorClasses(unittest.TestCase):
minimum=5e0,
maximum=1e2,
),
+ bilby.core.prior.Triangular(
+ name="test",
+ unit="unit",
+ minimum=-1.1,
+ maximum=3.14,
+ mode=0.,
+ ),
+ bilby.core.prior.Triangular(
+ name="test",
+ unit="unit",
+ minimum=0.,
+ maximum=4.,
+ mode=4.,
+ ),
+ bilby.core.prior.Triangular(
+ name="test",
+ unit="unit",
+ minimum=2.,
+ maximum=5.,
+ mode=2.,
+ ),
bilby.gw.prior.ConditionalUniformComovingVolume(
condition_func=condition_func, name="redshift", minimum=0.1, maximum=1.0
),
@@ -702,6 +723,7 @@ class TestPriorClasses(unittest.TestCase):
bilby.core.prior.Gamma,
bilby.core.prior.MultivariateGaussian,
bilby.core.prior.FermiDirac,
+ bilby.core.prior.Triangular,
bilby.gw.prior.HealPixPrior,
),
):
@@ -725,6 +747,7 @@ class TestPriorClasses(unittest.TestCase):
bilby.core.prior.Gamma,
bilby.core.prior.MultivariateGaussian,
bilby.core.prior.FermiDirac,
+ bilby.core.prior.Triangular,
bilby.gw.prior.HealPixPrior,
),
):
diff --git a/test/gw/conversion_test.py b/test/gw/conversion_test.py
index 6cf05de66275debe90bb6069c04947507f7e1f3a..f484019c194c6c9c9f9e1f3a1196f2f1bc2a2f2e 100644
--- a/test/gw/conversion_test.py
+++ b/test/gw/conversion_test.py
@@ -3,7 +3,6 @@ import unittest
import numpy as np
import pandas as pd
-
import bilby
from bilby.gw import conversion
@@ -658,5 +657,217 @@ class TestGenerateMassParameters(unittest.TestCase):
self.expected_values, source=True)
+class TestEquationOfStateConversions(unittest.TestCase):
+ '''
+ Class to test equation of state conversions.
+ The test points were generated from a simulation independent of bilby using the original lalsimulation calls.
+ Specific cases tested are described within each function.
+
+ '''
+ def setUp(self):
+ self.mass_1_source_spectral = [
+ 4.922542724434885,
+ 4.350626907771598,
+ 4.206155335439082,
+ 1.7822696459661311,
+ 1.3091740103047926
+ ]
+ self.mass_2_source_spectral = [
+ 3.459974694590303,
+ 1.2276461777181447,
+ 3.7287707089639976,
+ 0.3724016563531846,
+ 1.055042934805801
+ ]
+ self.spectral_pca_gamma_0 = [
+ 0.7074873121348357,
+ 0.05855931126849878,
+ 0.7795329261793462,
+ 1.467907561566463,
+ 2.9066488405635624
+ ]
+ self.spectral_pca_gamma_1 = [
+ -0.29807111670823816,
+ 2.027708558522935,
+ -1.4415775226512115,
+ -0.7104870098896858,
+ -0.4913817181089619
+ ]
+ self.spectral_pca_gamma_2 = [
+ 0.25625095371021156,
+ -0.19574096643220049,
+ -0.2710238103460012,
+ 0.22815820981582358,
+ -0.1543413205016374
+ ]
+ self.spectral_pca_gamma_3 = [
+ -0.04030365100175101,
+ 0.05698030777919032,
+ -0.045595911403040264,
+ -0.023480394227900117,
+ -0.07114492992285618
+ ]
+ self.spectral_gamma_0 = [
+ 1.1259406796075457,
+ 0.3191335618787259,
+ 1.3651245109783452,
+ 1.3540140238735314,
+ 1.4551949842961993
+ ]
+ self.spectral_gamma_1 = [
+ 0.26791504475282835,
+ 0.3930374252139248,
+ 0.11438399886108475,
+ 0.14181113477953,
+ -0.11989033256620368
+ ]
+ self.spectral_gamma_2 = [
+ -0.06810849354463173,
+ -0.038250139296677754,
+ -0.0801540229444505,
+ -0.05230330841791625,
+ -0.005197303281460286
+ ]
+ self.spectral_gamma_3 = [
+ 0.002848121360389597,
+ 0.000872447754855139,
+ 0.005528747386660879,
+ 0.0024325946344566484,
+ 0.00043890906202786106
+ ]
+ self.mass_1_source_polytrope = [
+ 2.2466565877822573,
+ 2.869741556013239,
+ 4.123897187899834,
+ 2.014160764697004,
+ 1.414796714032148,
+ 2.0919349759766614
+ ]
+ self.mass_2_source_polytrope = [
+ 0.36696047254774256,
+ 0.8580637120326807,
+ 1.650477659961306,
+ 1.310399737462001,
+ 0.5470843356210495,
+ 1.2311162283818198
+ ]
+ self.polytrope_log10_pressure_1 = [
+ 34.05849276958394,
+ 33.06962096113144,
+ 33.07579629429792,
+ 33.93412833210738,
+ 34.24096323517809,
+ 35.293288373856534
+ ]
+ self.polytrope_log10_pressure_2 = [
+ 33.82891829901602,
+ 35.14230456819543,
+ 34.940095188881976,
+ 34.72710820593933,
+ 35.42780071717415,
+ 35.648689969687915
+ ]
+ self.polytrope_gamma_0 = [
+ 2.359580734009537,
+ 2.3111471709796945,
+ 4.784129809424835,
+ 1.4900432021657437,
+ 1.0037220431922798,
+ 4.183994058757201
+ ]
+ self.polytrope_gamma_1 = [
+ 1.9497583698697314,
+ 1.0141111305083874,
+ 2.8228335336587826,
+ 4.032519623275465,
+ 1.10894361284508,
+ 3.168076721819637
+ ]
+ self.polytrope_gamma_2 = [
+ 4.6001755196585385,
+ 4.424090418206996,
+ 4.429607300132092,
+ 1.8176338276795763,
+ 2.9938859949129797,
+ 1.300271383168368
+ ]
+ self.lambda_1_spectral = [0., 0., 0., 0., 1275.7253186286332]
+ self.lambda_2_spectral = [0., 0., 0., 0., 4504.897675043909]
+ self.lambda_1_polytrope = [0., 0., 0., 0., 0., 234.66424898184766]
+ self.lambda_2_polytrope = [0., 0., 0., 0., 0., 3710.931378294547]
+ self.eos_check_spectral = [0, 0, 0, 0, 1]
+ self.eos_check_polytrope = [0, 0, 0, 0, 0, 1]
+
+ def test_spectral_pca_to_spectral(self):
+ for i in range(len(self.mass_1_source_spectral)):
+ spectral_gamma_0, spectral_gamma_1, spectral_gamma_2, spectral_gamma_3 = \
+ conversion.spectral_pca_to_spectral(
+ self.spectral_pca_gamma_0[i],
+ self.spectral_pca_gamma_1[i],
+ self.spectral_pca_gamma_2[i],
+ self.spectral_pca_gamma_3[i]
+ )
+ self.assertAlmostEqual(spectral_gamma_0, self.spectral_gamma_0[i], places=5)
+ self.assertAlmostEqual(spectral_gamma_1, self.spectral_gamma_1[i], places=5)
+ self.assertAlmostEqual(spectral_gamma_2, self.spectral_gamma_2[i], places=5)
+ self.assertAlmostEqual(spectral_gamma_3, self.spectral_gamma_3[i], places=5)
+
+ def test_spectral_params_to_lambda_1_lambda_2(self):
+ '''
+ The points cover 5 test cases:
+ - Fail SimNeutronStarEOS4ParamSDGammaCheck()
+ - Fail max_speed_of_sound_ <=1.1
+ - Fail mass_1_source <= max_mass
+ - Fail mass_2_source >= min_mass
+ - Passes all and produces accurate lambda_1, lambda_2, eos_check values
+ '''
+ for i in range(len(self.mass_1_source_spectral)):
+ spectral_gamma_0, spectral_gamma_1, spectral_gamma_2, spectral_gamma_3 = \
+ conversion.spectral_pca_to_spectral(
+ self.spectral_pca_gamma_0[i],
+ self.spectral_pca_gamma_1[i],
+ self.spectral_pca_gamma_2[i],
+ self.spectral_pca_gamma_3[i]
+ )
+ lambda_1, lambda_2, eos_check = \
+ conversion.spectral_params_to_lambda_1_lambda_2(
+ spectral_gamma_0,
+ spectral_gamma_1,
+ spectral_gamma_2,
+ spectral_gamma_3,
+ self.mass_1_source_spectral[i],
+ self.mass_2_source_spectral[i]
+ )
+ self.assertAlmostEqual(self.lambda_1_spectral[i], lambda_1, places=0)
+ self.assertAlmostEqual(self.lambda_2_spectral[i], lambda_2, places=0)
+ self.assertAlmostEqual(self.eos_check_spectral[i], eos_check)
+
+ def test_polytrope_or_causal_params_to_lambda_1_lambda_2_causal(self):
+ '''
+ The points cover 6 test cases:
+ - Fail log10_pressure1 >= log10_pressure2
+ - Fail SimNeutronStarEOS3PDViableFamilyCheck()
+ - Fail max_speed_of_sound_ <= 1.1
+ - Fail mass_1_source <= max_mass
+ - Fail mass_2_source >= min_mass
+ - Passes all and produces accurate lambda_1, lambda_2, eos_check values
+ '''
+ for i in range(len(self.mass_1_source_polytrope)):
+ lambda_1, lambda_2, eos_check = \
+ conversion.polytrope_or_causal_params_to_lambda_1_lambda_2(
+ self.polytrope_gamma_0[i],
+ self.polytrope_log10_pressure_1[i],
+ self.polytrope_gamma_1[i],
+ self.polytrope_log10_pressure_2[i],
+ self.polytrope_gamma_2[i],
+ self.mass_1_source_polytrope[i],
+ self.mass_2_source_polytrope[i],
+ 0
+ )
+ self.assertAlmostEqual(self.lambda_1_polytrope[i], lambda_1, places=3)
+ self.assertAlmostEqual(self.lambda_2_polytrope[i], lambda_2, places=3)
+ self.assertAlmostEqual(self.eos_check_polytrope[i], eos_check)
+
+
if __name__ == "__main__":
unittest.main()