Coverage for bilby/gw/conversion.py: 63%
898 statements
« prev ^ index » next coverage.py v7.6.1, created at 2025-05-06 04:57 +0000
1"""
2A collection of functions to convert between parameters describing
3gravitational-wave sources.
4"""
6import os
7import sys
8import multiprocessing
9import pickle
11import numpy as np
12from pandas import DataFrame, Series
13from scipy.stats import norm
15from .utils import (lalsim_SimNeutronStarEOS4ParamSDGammaCheck,
16 lalsim_SimNeutronStarEOS4ParameterSpectralDecomposition,
17 lalsim_SimNeutronStarEOS4ParamSDViableFamilyCheck,
18 lalsim_SimNeutronStarEOS3PieceDynamicPolytrope,
19 lalsim_SimNeutronStarEOS3PieceCausalAnalytic,
20 lalsim_SimNeutronStarEOS3PDViableFamilyCheck,
21 lalsim_CreateSimNeutronStarFamily,
22 lalsim_SimNeutronStarEOSMaxPseudoEnthalpy,
23 lalsim_SimNeutronStarEOSSpeedOfSoundGeometerized,
24 lalsim_SimNeutronStarFamMinimumMass,
25 lalsim_SimNeutronStarMaximumMass,
26 lalsim_SimNeutronStarRadius,
27 lalsim_SimNeutronStarLoveNumberK2)
29from ..core.likelihood import MarginalizedLikelihoodReconstructionError
30from ..core.utils import logger, solar_mass, gravitational_constant, speed_of_light, command_line_args, safe_file_dump
31from ..core.prior import DeltaFunction
32from .utils import lalsim_SimInspiralTransformPrecessingNewInitialConditions
33from .eos.eos import IntegrateTOV
34from .cosmology import get_cosmology, z_at_value
37def redshift_to_luminosity_distance(redshift, cosmology=None):
38 cosmology = get_cosmology(cosmology)
39 return cosmology.luminosity_distance(redshift).value
42def redshift_to_comoving_distance(redshift, cosmology=None):
43 cosmology = get_cosmology(cosmology)
44 return cosmology.comoving_distance(redshift).value
47def luminosity_distance_to_redshift(distance, cosmology=None):
48 from astropy import units
49 cosmology = get_cosmology(cosmology)
50 if isinstance(distance, Series):
51 distance = distance.values
52 return z_at_value(cosmology.luminosity_distance, distance * units.Mpc)
55def comoving_distance_to_redshift(distance, cosmology=None):
56 from astropy import units
57 cosmology = get_cosmology(cosmology)
58 if isinstance(distance, Series):
59 distance = distance.values
60 return z_at_value(cosmology.comoving_distance, distance * units.Mpc)
63def comoving_distance_to_luminosity_distance(distance, cosmology=None):
64 cosmology = get_cosmology(cosmology)
65 redshift = comoving_distance_to_redshift(distance, cosmology)
66 return redshift_to_luminosity_distance(redshift, cosmology)
69def luminosity_distance_to_comoving_distance(distance, cosmology=None):
70 cosmology = get_cosmology(cosmology)
71 redshift = luminosity_distance_to_redshift(distance, cosmology)
72 return redshift_to_comoving_distance(redshift, cosmology)
75_cosmology_docstring = """
76Convert from {input} to {output}
78Parameters
79----------
80{input}: float
81 The {input} to convert.
82cosmology: astropy.cosmology.Cosmology
83 The cosmology to use for the transformation.
84 See :code:`bilby.gw.cosmology.get_cosmology` for details of how to
85 specify this.
87Returns
88-------
89float
90 The {output} corresponding to the provided {input}.
91"""
93for _func in [
94 comoving_distance_to_luminosity_distance,
95 comoving_distance_to_redshift,
96 luminosity_distance_to_comoving_distance,
97 luminosity_distance_to_redshift,
98 redshift_to_comoving_distance,
99 redshift_to_luminosity_distance,
100]:
101 input, output = _func.__name__.split("_to_")
102 _func.__doc__ = _cosmology_docstring.format(input=input, output=output)
105def bilby_to_lalsimulation_spins(
106 theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2,
107 reference_frequency, phase
108):
109 """
110 Convert from Bilby spin parameters to lalsimulation ones.
112 All parameters are defined at the reference frequency and in SI units.
114 Parameters
115 ==========
116 theta_jn: float
117 Inclination angle
118 phi_jl: float
119 Spin phase angle
120 tilt_1: float
121 Primary object tilt
122 tilt_2: float
123 Secondary object tilt
124 phi_12: float
125 Relative spin azimuthal angle
126 a_1: float
127 Primary dimensionless spin magnitude
128 a_2: float
129 Secondary dimensionless spin magnitude
130 mass_1: float
131 Primary mass in SI units
132 mass_2: float
133 Secondary mass in SI units
134 reference_frequency: float
135 phase: float
136 Orbital phase
138 Returns
139 =======
140 iota: float
141 Transformed inclination
142 spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z: float
143 Cartesian spin components
144 """
145 if (a_1 == 0 or tilt_1 in [0, np.pi]) and (a_2 == 0 or tilt_2 in [0, np.pi]):
146 spin_1x = 0
147 spin_1y = 0
148 spin_1z = a_1 * np.cos(tilt_1)
149 spin_2x = 0
150 spin_2y = 0
151 spin_2z = a_2 * np.cos(tilt_2)
152 iota = theta_jn
153 else:
154 from numbers import Number
155 args = (
156 theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
157 mass_2, reference_frequency, phase
158 )
159 float_inputs = all([isinstance(arg, Number) for arg in args])
160 if float_inputs:
161 func = lalsim_SimInspiralTransformPrecessingNewInitialConditions
162 else:
163 func = transform_precessing_spins
164 iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = func(*args)
165 return iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z
168@np.vectorize
169def transform_precessing_spins(*args):
170 """
171 Vectorized wrapper for
172 lalsimulation.SimInspiralTransformPrecessingNewInitialConditions
174 For detailed documentation see
175 :code:`bilby.gw.conversion.bilby_to_lalsimulation_spins`.
176 This will be removed from the public API in a future release.
177 """
178 return lalsim_SimInspiralTransformPrecessingNewInitialConditions(*args)
181def convert_to_lal_binary_black_hole_parameters(parameters):
182 """
183 Convert parameters we have into parameters we need.
185 This is defined by the parameters of bilby.source.lal_binary_black_hole()
188 Mass: mass_1, mass_2
189 Spin: a_1, a_2, tilt_1, tilt_2, phi_12, phi_jl
190 Extrinsic: luminosity_distance, theta_jn, phase, ra, dec, geocent_time, psi
192 This involves popping a lot of things from parameters.
193 The keys in added_keys should be popped after evaluating the waveform.
195 Parameters
196 ==========
197 parameters: dict
198 dictionary of parameter values to convert into the required parameters
200 Returns
201 =======
202 converted_parameters: dict
203 dict of the required parameters
204 added_keys: list
205 keys which are added to parameters during function call
206 """
208 converted_parameters = parameters.copy()
209 original_keys = list(converted_parameters.keys())
210 if 'luminosity_distance' not in original_keys:
211 if 'redshift' in converted_parameters.keys():
212 converted_parameters['luminosity_distance'] = \
213 redshift_to_luminosity_distance(parameters['redshift'])
214 elif 'comoving_distance' in converted_parameters.keys():
215 converted_parameters['luminosity_distance'] = \
216 comoving_distance_to_luminosity_distance(
217 parameters['comoving_distance'])
219 for key in original_keys:
220 if key[-7:] == '_source':
221 if 'redshift' not in converted_parameters.keys():
222 converted_parameters['redshift'] =\
223 luminosity_distance_to_redshift(
224 parameters['luminosity_distance'])
225 converted_parameters[key[:-7]] = converted_parameters[key] * (
226 1 + converted_parameters['redshift'])
228 # we do not require the component masses be added if no mass parameters are present
229 converted_parameters = generate_component_masses(converted_parameters, require_add=False)
231 for idx in ['1', '2']:
232 key = 'chi_{}'.format(idx)
233 if key in original_keys:
234 if "chi_{}_in_plane".format(idx) in original_keys:
235 converted_parameters["a_{}".format(idx)] = (
236 converted_parameters[f"chi_{idx}"] ** 2
237 + converted_parameters[f"chi_{idx}_in_plane"] ** 2
238 ) ** 0.5
239 converted_parameters[f"cos_tilt_{idx}"] = (
240 converted_parameters[f"chi_{idx}"]
241 / converted_parameters[f"a_{idx}"]
242 )
243 elif "a_{}".format(idx) not in original_keys:
244 converted_parameters['a_{}'.format(idx)] = abs(
245 converted_parameters[key])
246 converted_parameters['cos_tilt_{}'.format(idx)] = \
247 np.sign(converted_parameters[key])
248 else:
249 with np.errstate(invalid="raise"):
250 try:
251 converted_parameters[f"cos_tilt_{idx}"] = (
252 converted_parameters[key] / converted_parameters[f"a_{idx}"]
253 )
254 except (FloatingPointError, ZeroDivisionError):
255 logger.debug(
256 "Error in conversion to spherical spin tilt. "
257 "This is often due to the spin parameters being zero. "
258 f"Setting cos_tilt_{idx} = 1."
259 )
260 converted_parameters[f"cos_tilt_{idx}"] = 1.0
262 for key in ["phi_jl", "phi_12"]:
263 if key not in converted_parameters:
264 converted_parameters[key] = 0.0
266 for angle in ['tilt_1', 'tilt_2', 'theta_jn']:
267 cos_angle = str('cos_' + angle)
268 if cos_angle in converted_parameters.keys():
269 with np.errstate(invalid="ignore"):
270 converted_parameters[angle] = np.arccos(converted_parameters[cos_angle])
272 if "delta_phase" in original_keys:
273 with np.errstate(invalid="ignore"):
274 converted_parameters["phase"] = np.mod(
275 converted_parameters["delta_phase"]
276 - np.sign(np.cos(converted_parameters["theta_jn"]))
277 * converted_parameters["psi"],
278 2 * np.pi)
279 added_keys = [key for key in converted_parameters.keys()
280 if key not in original_keys]
282 return converted_parameters, added_keys
285def convert_to_lal_binary_neutron_star_parameters(parameters):
286 """
287 Convert parameters we have into parameters we need.
289 This is defined by the parameters of bilby.source.lal_binary_black_hole()
292 Mass: mass_1, mass_2
293 Spin: a_1, a_2, tilt_1, tilt_2, phi_12, phi_jl
294 Extrinsic: luminosity_distance, theta_jn, phase, ra, dec, geocent_time, psi
295 Tidal: lambda_1, lamda_2, lambda_tilde, delta_lambda_tilde
297 This involves popping a lot of things from parameters.
298 The keys in added_keys should be popped after evaluating the waveform.
299 For details on tidal parameters see https://arxiv.org/pdf/1402.5156.pdf.
301 Parameters
302 ==========
303 parameters: dict
304 dictionary of parameter values to convert into the required parameters
306 Returns
307 =======
308 converted_parameters: dict
309 dict of the required parameters
310 added_keys: list
311 keys which are added to parameters during function call
312 """
313 converted_parameters = parameters.copy()
314 original_keys = list(converted_parameters.keys())
315 converted_parameters, added_keys =\
316 convert_to_lal_binary_black_hole_parameters(converted_parameters)
318 if not any([key in converted_parameters for key in
319 ['lambda_1', 'lambda_2',
320 'lambda_tilde', 'delta_lambda_tilde', 'lambda_symmetric',
321 'eos_polytrope_gamma_0', 'eos_spectral_pca_gamma_0', 'eos_v1']]):
322 converted_parameters['lambda_1'] = 0
323 converted_parameters['lambda_2'] = 0
324 added_keys = added_keys + ['lambda_1', 'lambda_2']
325 return converted_parameters, added_keys
327 if 'delta_lambda_tilde' in converted_parameters.keys():
328 converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\
329 lambda_tilde_delta_lambda_tilde_to_lambda_1_lambda_2(
330 converted_parameters['lambda_tilde'],
331 parameters['delta_lambda_tilde'], converted_parameters['mass_1'],
332 converted_parameters['mass_2'])
333 elif 'lambda_tilde' in converted_parameters.keys():
334 converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\
335 lambda_tilde_to_lambda_1_lambda_2(
336 converted_parameters['lambda_tilde'],
337 converted_parameters['mass_1'], converted_parameters['mass_2'])
338 if 'lambda_2' not in converted_parameters.keys() and 'lambda_1' in converted_parameters.keys():
339 converted_parameters['lambda_2'] =\
340 converted_parameters['lambda_1']\
341 * converted_parameters['mass_1']**5\
342 / converted_parameters['mass_2']**5
343 elif 'lambda_2' in converted_parameters.keys() and converted_parameters['lambda_2'] is None:
344 converted_parameters['lambda_2'] =\
345 converted_parameters['lambda_1']\
346 * converted_parameters['mass_1']**5\
347 / converted_parameters['mass_2']**5
348 elif 'eos_spectral_pca_gamma_0' in converted_parameters.keys(): # FIXME: This is a clunky way to do this
349 converted_parameters = generate_source_frame_parameters(converted_parameters)
350 float_eos_params = {}
351 max_len = 1
352 eos_keys = ['eos_spectral_pca_gamma_0',
353 'eos_spectral_pca_gamma_1',
354 'eos_spectral_pca_gamma_2',
355 'eos_spectral_pca_gamma_3',
356 'mass_1_source', 'mass_2_source']
357 for key in eos_keys:
358 try:
359 if (len(converted_parameters[key]) > max_len):
360 max_len = len(converted_parameters[key])
361 except TypeError:
362 float_eos_params[key] = converted_parameters[key]
363 if len(float_eos_params) == len(eos_keys): # case where all eos params are floats (pinned)
364 g0, g1, g2, g3 = spectral_pca_to_spectral(
365 converted_parameters['eos_spectral_pca_gamma_0'],
366 converted_parameters['eos_spectral_pca_gamma_1'],
367 converted_parameters['eos_spectral_pca_gamma_2'],
368 converted_parameters['eos_spectral_pca_gamma_3'])
369 converted_parameters['lambda_1'], converted_parameters['lambda_2'], converted_parameters['eos_check'] = \
370 spectral_params_to_lambda_1_lambda_2(
371 g0, g1, g2, g3, converted_parameters['mass_1_source'], converted_parameters['mass_2_source'])
372 elif len(float_eos_params) < len(eos_keys): # case where some or none of the eos params are floats (pinned)
373 for key in float_eos_params.keys():
374 converted_parameters[key] = np.ones(max_len) * converted_parameters[key]
375 g0pca = converted_parameters['eos_spectral_pca_gamma_0']
376 g1pca = converted_parameters['eos_spectral_pca_gamma_1']
377 g2pca = converted_parameters['eos_spectral_pca_gamma_2']
378 g3pca = converted_parameters['eos_spectral_pca_gamma_3']
379 m1s = converted_parameters['mass_1_source']
380 m2s = converted_parameters['mass_2_source']
381 all_lambda_1 = np.empty(0)
382 all_lambda_2 = np.empty(0)
383 all_eos_check = np.empty(0, dtype=bool)
384 for (g_0pca, g_1pca, g_2pca, g_3pca, m1_s, m2_s) in zip(g0pca, g1pca, g2pca, g3pca, m1s, m2s):
385 g_0, g_1, g_2, g_3 = spectral_pca_to_spectral(g_0pca, g_1pca, g_2pca, g_3pca)
386 lambda_1, lambda_2, eos_check = \
387 spectral_params_to_lambda_1_lambda_2(g_0, g_1, g_2, g_3, m1_s, m2_s)
388 all_lambda_1 = np.append(all_lambda_1, lambda_1)
389 all_lambda_2 = np.append(all_lambda_2, lambda_2)
390 all_eos_check = np.append(all_eos_check, eos_check)
391 converted_parameters['lambda_1'] = all_lambda_1
392 converted_parameters['lambda_2'] = all_lambda_2
393 converted_parameters['eos_check'] = all_eos_check
394 for key in float_eos_params.keys():
395 converted_parameters[key] = float_eos_params[key]
396 elif 'eos_polytrope_gamma_0' and 'eos_polytrope_log10_pressure_1' in converted_parameters.keys():
397 converted_parameters = generate_source_frame_parameters(converted_parameters)
398 float_eos_params = {}
399 max_len = 1
400 eos_keys = ['eos_polytrope_gamma_0',
401 'eos_polytrope_gamma_1',
402 'eos_polytrope_gamma_2',
403 'eos_polytrope_log10_pressure_1',
404 'eos_polytrope_log10_pressure_2',
405 'mass_1_source', 'mass_2_source']
406 for key in eos_keys:
407 try:
408 if (len(converted_parameters[key]) > max_len):
409 max_len = len(converted_parameters[key])
410 except TypeError:
411 float_eos_params[key] = converted_parameters[key]
412 if len(float_eos_params) == len(eos_keys): # case where all eos params are floats (pinned)
413 converted_parameters['lambda_1'], converted_parameters['lambda_2'], converted_parameters['eos_check'] = \
414 polytrope_or_causal_params_to_lambda_1_lambda_2(
415 converted_parameters['eos_polytrope_gamma_0'],
416 converted_parameters['eos_polytrope_log10_pressure_1'],
417 converted_parameters['eos_polytrope_gamma_1'],
418 converted_parameters['eos_polytrope_log10_pressure_2'],
419 converted_parameters['eos_polytrope_gamma_2'],
420 converted_parameters['mass_1_source'],
421 converted_parameters['mass_2_source'],
422 causal=0)
423 elif len(float_eos_params) < len(eos_keys): # case where some or none are floats (pinned)
424 for key in float_eos_params.keys():
425 converted_parameters[key] = np.ones(max_len) * converted_parameters[key]
426 pg0 = converted_parameters['eos_polytrope_gamma_0']
427 pg1 = converted_parameters['eos_polytrope_gamma_1']
428 pg2 = converted_parameters['eos_polytrope_gamma_2']
429 logp1 = converted_parameters['eos_polytrope_log10_pressure_1']
430 logp2 = converted_parameters['eos_polytrope_log10_pressure_2']
431 m1s = converted_parameters['mass_1_source']
432 m2s = converted_parameters['mass_2_source']
433 all_lambda_1 = np.empty(0)
434 all_lambda_2 = np.empty(0)
435 all_eos_check = np.empty(0, dtype=bool)
436 for (pg_0, pg_1, pg_2, logp_1, logp_2, m1_s, m2_s) in zip(pg0, pg1, pg2, logp1, logp2, m1s, m2s):
437 lambda_1, lambda_2, eos_check = \
438 polytrope_or_causal_params_to_lambda_1_lambda_2(
439 pg_0, logp_1, pg_1, logp_2, pg_2, m1_s, m2_s, causal=0)
440 all_lambda_1 = np.append(all_lambda_1, lambda_1)
441 all_lambda_2 = np.append(all_lambda_2, lambda_2)
442 all_eos_check = np.append(all_eos_check, eos_check)
443 converted_parameters['lambda_1'] = all_lambda_1
444 converted_parameters['lambda_2'] = all_lambda_2
445 converted_parameters['eos_check'] = all_eos_check
446 for key in float_eos_params.keys():
447 converted_parameters[key] = float_eos_params[key]
448 elif 'eos_polytrope_gamma_0' and 'eos_polytrope_scaled_pressure_ratio' in converted_parameters.keys():
449 converted_parameters = generate_source_frame_parameters(converted_parameters)
450 float_eos_params = {}
451 max_len = 1
452 eos_keys = ['eos_polytrope_gamma_0',
453 'eos_polytrope_gamma_1',
454 'eos_polytrope_gamma_2',
455 'eos_polytrope_scaled_pressure_ratio',
456 'eos_polytrope_scaled_pressure_2',
457 'mass_1_source', 'mass_2_source']
458 for key in eos_keys:
459 try:
460 if (len(converted_parameters[key]) > max_len):
461 max_len = len(converted_parameters[key])
462 except TypeError:
463 float_eos_params[key] = converted_parameters[key]
464 if len(float_eos_params) == len(eos_keys): # case where all eos params are floats (pinned)
465 logp1, logp2 = log_pressure_reparameterization_conversion(
466 converted_parameters['eos_polytrope_scaled_pressure_ratio'],
467 converted_parameters['eos_polytrope_scaled_pressure_2'])
468 converted_parameters['lambda_1'], converted_parameters['lambda_2'], converted_parameters['eos_check'] = \
469 polytrope_or_causal_params_to_lambda_1_lambda_2(
470 converted_parameters['eos_polytrope_gamma_0'],
471 logp1,
472 converted_parameters['eos_polytrope_gamma_1'],
473 logp2,
474 converted_parameters['eos_polytrope_gamma_2'],
475 converted_parameters['mass_1_source'],
476 converted_parameters['mass_2_source'],
477 causal=0)
478 elif len(float_eos_params) < len(eos_keys): # case where some or none are floats (pinned)
479 for key in float_eos_params.keys():
480 converted_parameters[key] = np.ones(max_len) * converted_parameters[key]
481 pg0 = converted_parameters['eos_polytrope_gamma_0']
482 pg1 = converted_parameters['eos_polytrope_gamma_1']
483 pg2 = converted_parameters['eos_polytrope_gamma_2']
484 scaledratio = converted_parameters['eos_polytrope_scaled_pressure_ratio']
485 scaled_p2 = converted_parameters['eos_polytrope_scaled_pressure_2']
486 m1s = converted_parameters['mass_1_source']
487 m2s = converted_parameters['mass_2_source']
488 all_lambda_1 = np.empty(0)
489 all_lambda_2 = np.empty(0)
490 all_eos_check = np.empty(0, dtype=bool)
491 for (pg_0, pg_1, pg_2, scaled_ratio, scaled_p_2, m1_s,
492 m2_s) in zip(pg0, pg1, pg2, scaledratio, scaled_p2, m1s, m2s):
493 logp_1, logp_2 = log_pressure_reparameterization_conversion(scaled_ratio, scaled_p_2)
494 lambda_1, lambda_2, eos_check = \
495 polytrope_or_causal_params_to_lambda_1_lambda_2(
496 pg_0, logp_1, pg_1, logp_2, pg_2, m1_s, m2_s, causal=0)
497 all_lambda_1 = np.append(all_lambda_1, lambda_1)
498 all_lambda_2 = np.append(all_lambda_2, lambda_2)
499 all_eos_check = np.append(all_eos_check, eos_check)
500 converted_parameters['lambda_1'] = all_lambda_1
501 converted_parameters['lambda_2'] = all_lambda_2
502 converted_parameters['eos_check'] = all_eos_check
503 for key in float_eos_params.keys():
504 converted_parameters[key] = float_eos_params[key]
505 elif 'eos_v1' in converted_parameters.keys():
506 converted_parameters = generate_source_frame_parameters(converted_parameters)
507 float_eos_params = {}
508 max_len = 1
509 eos_keys = ['eos_v1',
510 'eos_v2',
511 'eos_v3',
512 'eos_log10_pressure1_cgs',
513 'eos_log10_pressure2_cgs',
514 'mass_1_source', 'mass_2_source']
515 for key in eos_keys:
516 try:
517 if (len(converted_parameters[key]) > max_len):
518 max_len = len(converted_parameters[key])
519 except TypeError:
520 float_eos_params[key] = converted_parameters[key]
521 if len(float_eos_params) == len(eos_keys): # case where all eos params are floats (pinned)
522 converted_parameters['lambda_1'], converted_parameters['lambda_2'], converted_parameters['eos_check'] = \
523 polytrope_or_causal_params_to_lambda_1_lambda_2(
524 converted_parameters['eos_v1'],
525 converted_parameters['eos_log10_pressure1_cgs'],
526 converted_parameters['eos_v2'],
527 converted_parameters['eos_log10_pressure2_cgs'],
528 converted_parameters['eos_v3'],
529 converted_parameters['mass_1_source'],
530 converted_parameters['mass_2_source'],
531 causal=1)
532 elif len(float_eos_params) < len(eos_keys): # case where some or none are floats (pinned)
533 for key in float_eos_params.keys():
534 converted_parameters[key] = np.ones(max_len) * converted_parameters[key]
535 v1 = converted_parameters['eos_v1']
536 v2 = converted_parameters['eos_v2']
537 v3 = converted_parameters['eos_v3']
538 logp1 = converted_parameters['eos_log10_pressure1_cgs']
539 logp2 = converted_parameters['eos_log10_pressure2_cgs']
540 m1s = converted_parameters['mass_1_source']
541 m2s = converted_parameters['mass_2_source']
542 all_lambda_1 = np.empty(0)
543 all_lambda_2 = np.empty(0)
544 all_eos_check = np.empty(0, dtype=bool)
545 for (v_1, v_2, v_3, logp_1, logp_2, m1_s, m2_s) in zip(v1, v2, v3, logp1, logp2, m1s, m2s):
546 lambda_1, lambda_2, eos_check = \
547 polytrope_or_causal_params_to_lambda_1_lambda_2(
548 v_1, logp_1, v_2, logp_2, v_3, m1_s, m2_s, causal=1)
549 all_lambda_1 = np.append(all_lambda_1, lambda_1)
550 all_lambda_2 = np.append(all_lambda_2, lambda_2)
551 all_eos_check = np.append(all_eos_check, eos_check)
552 converted_parameters['lambda_1'] = all_lambda_1
553 converted_parameters['lambda_2'] = all_lambda_2
554 converted_parameters['eos_check'] = all_eos_check
555 for key in float_eos_params.keys():
556 converted_parameters[key] = float_eos_params[key]
557 elif 'lambda_symmetric' in converted_parameters.keys():
558 if 'lambda_antisymmetric' in converted_parameters.keys():
559 converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\
560 lambda_symmetric_lambda_antisymmetric_to_lambda_1_lambda_2(
561 converted_parameters['lambda_symmetric'],
562 converted_parameters['lambda_antisymmetric'])
563 elif 'mass_ratio' in converted_parameters.keys():
564 if 'binary_love_uniform' in converted_parameters.keys():
565 converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\
566 binary_love_lambda_symmetric_to_lambda_1_lambda_2_manual_marginalisation(
567 converted_parameters['binary_love_uniform'],
568 converted_parameters['lambda_symmetric'],
569 converted_parameters['mass_ratio'])
570 else:
571 converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\
572 binary_love_lambda_symmetric_to_lambda_1_lambda_2_automatic_marginalisation(
573 converted_parameters['lambda_symmetric'],
574 converted_parameters['mass_ratio'])
576 added_keys = [key for key in converted_parameters.keys()
577 if key not in original_keys]
579 return converted_parameters, added_keys
582def log_pressure_reparameterization_conversion(scaled_pressure_ratio, scaled_pressure_2):
583 '''
584 Converts the reparameterization joining pressures from
585 (scaled_pressure_ratio,scaled_pressure_2) to (log10_pressure_1,log10_pressure_2).
586 This reparameterization with a triangular prior (with mode = max) on scaled_pressure_2
587 and a uniform prior on scaled_pressure_ratio
588 mimics identical uniform priors on log10_pressure_1 and log10_pressure_2
589 where samples with log10_pressure_2 > log10_pressure_1 are rejected.
590 This reparameterization allows for a faster initialization.
591 A minimum log10_pressure of 33 (in cgs units) is chosen to be slightly higher than the low-density crust EOS
592 that is stitched to the dynamic polytrope EOS model in LALSimulation.
594 Parameters
595 ----------
596 scaled_pressure_ratio, scaled_pressure_2: float
597 reparameterizations of the dividing pressures
599 Returns
600 -------
601 log10_pressure_1, log10_pressure_2: float
602 joining pressures in the original parameterization
604 '''
605 minimum_pressure = 33.
606 log10_pressure_1 = (scaled_pressure_ratio * scaled_pressure_2) + minimum_pressure
607 log10_pressure_2 = minimum_pressure + scaled_pressure_2
609 return log10_pressure_1, log10_pressure_2
612def spectral_pca_to_spectral(gamma_pca_0, gamma_pca_1, gamma_pca_2, gamma_pca_3):
613 '''
614 Change of basis on parameter space
615 from an efficient space to sample in (sample space)
616 to the space used in spectral eos model (model space).
617 Uses principle component analysis (PCA) to sample spectral gamma parameters more efficiently.
618 See arxiv:2001.01747 for the PCA conversion transformation matrix, mean, and standard deviation.
619 (Note that this transformation matrix is the inverse of the one referenced
620 in order to perform the inverse transformation.)
622 Parameters
623 ----------
624 gamma_pca_0, gamma_pca_1, gamma_pca_2, gamma_pca_3: float
625 Spectral gamma PCA parameters to be converted to spectral gamma parameters.
627 Returns
628 -------
629 converted_gamma_parameters: np.array()
630 array of gamma_0, gamma_1, gamma_2, gamma_3 in model space
632 '''
633 sampled_pca_gammas = np.array([gamma_pca_0, gamma_pca_1, gamma_pca_2, gamma_pca_3])
634 transformation_matrix = np.array(
635 [
636 [0.43801, -0.76705, 0.45143, 0.12646],
637 [-0.53573, 0.17169, 0.67968, 0.47070],
638 [0.52660, 0.31255, -0.19454, 0.76626],
639 [-0.49379, -0.53336, -0.54444, 0.41868]
640 ]
641 )
643 model_space_mean = np.array([0.89421, 0.33878, -0.07894, 0.00393])
644 model_space_standard_deviation = np.array([0.35700, 0.25769, 0.05452, 0.00312])
645 converted_gamma_parameters = \
646 model_space_mean + model_space_standard_deviation * np.dot(transformation_matrix, sampled_pca_gammas)
648 return converted_gamma_parameters
651def spectral_params_to_lambda_1_lambda_2(gamma_0, gamma_1, gamma_2, gamma_3, mass_1_source, mass_2_source):
652 '''
653 Converts from the 4 spectral decomposition parameters and the source masses
654 to the tidal deformability parameters.
656 Parameters
657 ----------
658 gamma_0, gamma_1, gamma_2, gamma_3: float
659 sampled spectral decomposition parameters
660 mass_1_source, mass_2_source: float
661 sampled component mass parameters converted to source frame in solar masses
663 Returns
664 -------
665 lambda_1, lambda_2: float
666 component tidal deformability parameters
667 eos_check: bool
668 whether or not the equation of state is viable /
669 if eos_check = False, lambdas are 0 and the sample is rejected.
671 '''
672 eos_check = True
673 if lalsim_SimNeutronStarEOS4ParamSDGammaCheck(gamma_0, gamma_1, gamma_2, gamma_3) != 0:
674 lambda_1 = 0.0
675 lambda_2 = 0.0
676 eos_check = False
677 else:
678 eos = lalsim_SimNeutronStarEOS4ParameterSpectralDecomposition(gamma_0, gamma_1, gamma_2, gamma_3)
679 if lalsim_SimNeutronStarEOS4ParamSDViableFamilyCheck(gamma_0, gamma_1, gamma_2, gamma_3) != 0:
680 lambda_1 = 0.0
681 lambda_2 = 0.0
682 eos_check = False
683 else:
684 lambda_1, lambda_2, eos_check = neutron_star_family_physical_check(eos, mass_1_source, mass_2_source)
686 return lambda_1, lambda_2, eos_check
689def polytrope_or_causal_params_to_lambda_1_lambda_2(
690 param1, log10_pressure1_cgs, param2, log10_pressure2_cgs, param3, mass_1_source, mass_2_source, causal):
691 """
692 Converts parameters from sampled dynamic piecewise polytrope parameters
693 to component tidal deformablity parameters.
694 Enforces log10_pressure1_cgs < log10_pressure2_cgs.
695 Checks number of points in the equation of state for viability.
696 Note that subtracting 1 from the log10 pressure in cgs converts it to
697 log10 pressure in si units.
699 Parameters
700 ----------
701 param1, param2, param3: float
702 either the sampled adiabatic indices in piecewise polytrope model
703 or the sampled causal model params v1, v2, v3
704 log10_pressure1_cgs, log10_pressure2_cgs: float
705 dividing pressures in piecewise polytrope model or causal model
706 mass_1_source, mass_2_source: float
707 source frame component mass parameters in Msuns
708 causal: bool
709 whether or not to use causal polytrope model
710 1 - causal; 0 - not causal
712 Returns
713 -------
714 lambda_1: float
715 tidal deformability parameter associated with mass 1
716 lambda_2: float
717 tidal deformability parameter associated with mass 2
718 eos_check: bool
719 whether eos is valid or not
721 """
722 eos_check = True
723 if log10_pressure1_cgs >= log10_pressure2_cgs:
724 lambda_1 = 0.0
725 lambda_2 = 0.0
726 eos_check = False
727 else:
728 if causal == 0:
729 eos = lalsim_SimNeutronStarEOS3PieceDynamicPolytrope(
730 param1, log10_pressure1_cgs - 1., param2, log10_pressure2_cgs - 1., param3)
731 else:
732 eos = lalsim_SimNeutronStarEOS3PieceCausalAnalytic(
733 param1, log10_pressure1_cgs - 1., param2, log10_pressure2_cgs - 1., param3)
734 if lalsim_SimNeutronStarEOS3PDViableFamilyCheck(
735 param1, log10_pressure1_cgs - 1., param2, log10_pressure2_cgs - 1., param3, causal) != 0:
736 lambda_1 = 0.0
737 lambda_2 = 0.0
738 eos_check = False
739 else:
740 lambda_1, lambda_2, eos_check = neutron_star_family_physical_check(eos, mass_1_source, mass_2_source)
742 return lambda_1, lambda_2, eos_check
745def neutron_star_family_physical_check(eos, mass_1_source, mass_2_source):
746 """
747 Takes in a lalsim eos object. Performs causal and max/min mass eos checks.
748 Calculates component lambdas if eos object passes causality.
749 Returns lambda = 0 if not.
751 Parameters
752 ----------
753 eos: lalsim swig-wrapped eos object
754 the neutron star equation of state
755 mass_1_source, mass_2_source: float
756 source frame component masses 1 and 2 in solar masses
758 Returns
759 -------
760 lambda_1, lambda_2: float
761 component tidal deformability parameters
762 eos_check: bool
763 whether or not the equation of state is physically allowed
765 """
766 eos_check = True
767 family = lalsim_CreateSimNeutronStarFamily(eos)
768 max_pseudo_enthalpy = lalsim_SimNeutronStarEOSMaxPseudoEnthalpy(eos)
769 max_speed_of_sound = lalsim_SimNeutronStarEOSSpeedOfSoundGeometerized(max_pseudo_enthalpy, eos)
770 min_mass = lalsim_SimNeutronStarFamMinimumMass(family) / solar_mass
771 max_mass = lalsim_SimNeutronStarMaximumMass(family) / solar_mass
772 if max_speed_of_sound <= 1.1 and min_mass <= mass_1_source <= max_mass and min_mass <= mass_2_source <= max_mass:
773 lambda_1 = lambda_from_mass_and_family(mass_1_source, family)
774 lambda_2 = lambda_from_mass_and_family(mass_2_source, family)
775 else:
776 lambda_1 = 0.0
777 lambda_2 = 0.0
778 eos_check = False
780 return lambda_1, lambda_2, eos_check
783def lambda_from_mass_and_family(mass_i, family):
784 """
785 Convert from equation of state model parameters to
786 component tidal parameters.
788 Parameters
789 ----------
790 family: lalsim family object
791 EOS family of type lalsimulation.SimNeutronStarFamily.
792 mass_i: Component mass of neutron star in solar masses.
794 Returns
795 -------
796 lambda_1: float
797 component tidal deformability parameter
799 """
800 radius = lalsim_SimNeutronStarRadius(mass_i * solar_mass, family)
801 love_number_k2 = lalsim_SimNeutronStarLoveNumberK2(mass_i * solar_mass, family)
802 mass_geometrized = mass_i * solar_mass * gravitational_constant / speed_of_light ** 2.
803 compactness = mass_geometrized / radius
804 lambda_i = (2. / 3.) * love_number_k2 / compactness ** 5.
806 return lambda_i
809def eos_family_physical_check(eos):
810 """
811 Function that determines if the EoS family contains
812 sufficient number of points before maximum mass is reached.
814 e_min is chosen to be sufficiently small so that the entire
815 EoS is captured when converting to mass-radius space.
817 Returns True if family is valid, False if not.
818 """
819 e_min = 1.6e-10
820 e_central = eos.e_pdat[-1, 1]
821 loge_min = np.log(e_min)
822 loge_central = np.log(e_central)
823 logedat = np.linspace(loge_min, loge_central, num=eos.npts)
824 edat = np.exp(logedat)
826 # Generate m, r, and k2 lists
827 mdat = []
828 rdat = []
829 k2dat = []
830 for i in range(8):
831 tov_solver = IntegrateTOV(eos, edat[i])
832 m, r, k2 = tov_solver.integrate_TOV()
833 mdat.append(m)
834 rdat.append(r)
835 k2dat.append(k2)
837 # Check if maximum mass has been found
838 if i > 0 and mdat[i] <= mdat[i - 1]:
839 break
841 if len(mdat) < 8:
842 return False
843 else:
844 return True
847def total_mass_and_mass_ratio_to_component_masses(mass_ratio, total_mass):
848 """
849 Convert total mass and mass ratio of a binary to its component masses.
851 Parameters
852 ==========
853 mass_ratio: float
854 Mass ratio (mass_2/mass_1) of the binary
855 total_mass: float
856 Total mass of the binary
858 Returns
859 =======
860 mass_1: float
861 Mass of the heavier object
862 mass_2: float
863 Mass of the lighter object
864 """
866 mass_1 = total_mass / (1 + mass_ratio)
867 mass_2 = mass_1 * mass_ratio
868 return mass_1, mass_2
871def chirp_mass_and_mass_ratio_to_component_masses(chirp_mass, mass_ratio):
872 """
873 Convert total mass and mass ratio of a binary to its component masses.
875 Parameters
876 ==========
877 chirp_mass: float
878 Chirp mass of the binary
879 mass_ratio: float
880 Mass ratio (mass_2/mass_1) of the binary
882 Returns
883 =======
884 mass_1: float
885 Mass of the heavier object
886 mass_2: float
887 Mass of the lighter object
888 """
889 total_mass = chirp_mass_and_mass_ratio_to_total_mass(chirp_mass=chirp_mass,
890 mass_ratio=mass_ratio)
891 mass_1, mass_2 = (
892 total_mass_and_mass_ratio_to_component_masses(
893 total_mass=total_mass, mass_ratio=mass_ratio)
894 )
895 return mass_1, mass_2
898def symmetric_mass_ratio_to_mass_ratio(symmetric_mass_ratio):
899 """
900 Convert the symmetric mass ratio to the normal mass ratio.
902 Parameters
903 ==========
904 symmetric_mass_ratio: float
905 Symmetric mass ratio of the binary
907 Returns
908 =======
909 mass_ratio: float
910 Mass ratio of the binary
911 """
913 temp = (1 / symmetric_mass_ratio / 2 - 1)
914 return temp - (temp ** 2 - 1) ** 0.5
917def chirp_mass_and_total_mass_to_symmetric_mass_ratio(chirp_mass, total_mass):
918 """
919 Convert chirp mass and total mass of a binary to its symmetric mass ratio.
921 Parameters
922 ==========
923 chirp_mass: float
924 Chirp mass of the binary
925 total_mass: float
926 Total mass of the binary
928 Returns
929 =======
930 symmetric_mass_ratio: float
931 Symmetric mass ratio of the binary
932 """
934 return (chirp_mass / total_mass) ** (5 / 3)
937def chirp_mass_and_primary_mass_to_mass_ratio(chirp_mass, mass_1):
938 """
939 Convert chirp mass and mass ratio of a binary to its total mass.
941 Rearranging the relation for chirp mass (as a function of mass_1 and
942 mass_2) and q = mass_2 / mass_1, it can be shown that
944 (chirp_mass/mass_1)^5 = q^3 / (1 + q)
946 Solving for q, we find the relation expressed in python below for q.
948 Parameters
949 ==========
950 chirp_mass: float
951 Chirp mass of the binary
952 mass_1: float
953 The primary mass
955 Returns
956 =======
957 mass_ratio: float
958 Mass ratio (mass_2/mass_1) of the binary
959 """
960 a = (chirp_mass / mass_1) ** 5
961 t0 = np.cbrt(9 * a + np.sqrt(3) * np.sqrt(27 * a ** 2 - 4 * a ** 3))
962 t1 = np.cbrt(2) * 3 ** (2 / 3)
963 t2 = np.cbrt(2 / 3) * a
964 return t2 / t0 + t0 / t1
967def chirp_mass_and_mass_ratio_to_total_mass(chirp_mass, mass_ratio):
968 """
969 Convert chirp mass and mass ratio of a binary to its total mass.
971 Parameters
972 ==========
973 chirp_mass: float
974 Chirp mass of the binary
975 mass_ratio: float
976 Mass ratio (mass_2/mass_1) of the binary
978 Returns
979 =======
980 mass_1: float
981 Mass of the heavier object
982 mass_2: float
983 Mass of the lighter object
984 """
986 with np.errstate(invalid="ignore"):
987 return chirp_mass * (1 + mass_ratio) ** 1.2 / mass_ratio ** 0.6
990def component_masses_to_chirp_mass(mass_1, mass_2):
991 """
992 Convert the component masses of a binary to its chirp mass.
994 Parameters
995 ==========
996 mass_1: float
997 Mass of the heavier object
998 mass_2: float
999 Mass of the lighter object
1001 Returns
1002 =======
1003 chirp_mass: float
1004 Chirp mass of the binary
1005 """
1007 return (mass_1 * mass_2) ** 0.6 / (mass_1 + mass_2) ** 0.2
1010def component_masses_to_total_mass(mass_1, mass_2):
1011 """
1012 Convert the component masses of a binary to its total mass.
1014 Parameters
1015 ==========
1016 mass_1: float
1017 Mass of the heavier object
1018 mass_2: float
1019 Mass of the lighter object
1021 Returns
1022 =======
1023 total_mass: float
1024 Total mass of the binary
1025 """
1027 return mass_1 + mass_2
1030def component_masses_to_symmetric_mass_ratio(mass_1, mass_2):
1031 """
1032 Convert the component masses of a binary to its symmetric mass ratio.
1034 Parameters
1035 ==========
1036 mass_1: float
1037 Mass of the heavier object
1038 mass_2: float
1039 Mass of the lighter object
1041 Returns
1042 =======
1043 symmetric_mass_ratio: float
1044 Symmetric mass ratio of the binary
1045 """
1047 return np.minimum((mass_1 * mass_2) / (mass_1 + mass_2) ** 2, 1 / 4)
1050def component_masses_to_mass_ratio(mass_1, mass_2):
1051 """
1052 Convert the component masses of a binary to its chirp mass.
1054 Parameters
1055 ==========
1056 mass_1: float
1057 Mass of the heavier object
1058 mass_2: float
1059 Mass of the lighter object
1061 Returns
1062 =======
1063 mass_ratio: float
1064 Mass ratio of the binary
1065 """
1067 return mass_2 / mass_1
1070def mass_1_and_chirp_mass_to_mass_ratio(mass_1, chirp_mass):
1071 """
1072 Calculate mass ratio from mass_1 and chirp_mass.
1074 This involves solving mc = m1 * q**(3/5) / (1 + q)**(1/5).
1076 Parameters
1077 ==========
1078 mass_1: float
1079 Mass of the heavier object
1080 chirp_mass: float
1081 Mass of the lighter object
1083 Returns
1084 =======
1085 mass_ratio: float
1086 Mass ratio of the binary
1087 """
1088 temp = (chirp_mass / mass_1) ** 5
1089 mass_ratio = (2 / 3 / (3 ** 0.5 * (27 * temp ** 2 - 4 * temp ** 3) ** 0.5 +
1090 9 * temp)) ** (1 / 3) * temp + \
1091 ((3 ** 0.5 * (27 * temp ** 2 - 4 * temp ** 3) ** 0.5 +
1092 9 * temp) / (2 * 3 ** 2)) ** (1 / 3)
1093 return mass_ratio
1096def mass_2_and_chirp_mass_to_mass_ratio(mass_2, chirp_mass):
1097 """
1098 Calculate mass ratio from mass_1 and chirp_mass.
1100 This involves solving mc = m2 * (1/q)**(3/5) / (1 + (1/q))**(1/5).
1102 Parameters
1103 ==========
1104 mass_2: float
1105 Mass of the lighter object
1106 chirp_mass: float
1107 Chirp mass of the binary
1109 Returns
1110 =======
1111 mass_ratio: float
1112 Mass ratio of the binary
1113 """
1114 # Passing mass_2, the expression from the function above
1115 # returns 1/q (because chirp mass is invariant under
1116 # mass_1 <-> mass_2)
1117 return 1 / mass_1_and_chirp_mass_to_mass_ratio(mass_2, chirp_mass)
1120def lambda_1_lambda_2_to_lambda_tilde(lambda_1, lambda_2, mass_1, mass_2):
1121 """
1122 Convert from individual tidal parameters to domainant tidal term.
1124 See, e.g., Wade et al., https://arxiv.org/pdf/1402.5156.pdf.
1126 Parameters
1127 ==========
1128 lambda_1: float
1129 Tidal parameter of more massive neutron star.
1130 lambda_2: float
1131 Tidal parameter of less massive neutron star.
1132 mass_1: float
1133 Mass of more massive neutron star.
1134 mass_2: float
1135 Mass of less massive neutron star.
1137 Returns
1138 =======
1139 lambda_tilde: float
1140 Dominant tidal term.
1142 """
1143 eta = component_masses_to_symmetric_mass_ratio(mass_1, mass_2)
1144 lambda_plus = lambda_1 + lambda_2
1145 lambda_minus = lambda_1 - lambda_2
1146 lambda_tilde = 8 / 13 * (
1147 (1 + 7 * eta - 31 * eta**2) * lambda_plus +
1148 (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta**2) * lambda_minus)
1150 return lambda_tilde
1153def lambda_1_lambda_2_to_delta_lambda_tilde(lambda_1, lambda_2, mass_1, mass_2):
1154 """
1155 Convert from individual tidal parameters to second domainant tidal term.
1157 See, e.g., Wade et al., https://arxiv.org/pdf/1402.5156.pdf.
1159 Parameters
1160 ==========
1161 lambda_1: float
1162 Tidal parameter of more massive neutron star.
1163 lambda_2: float
1164 Tidal parameter of less massive neutron star.
1165 mass_1: float
1166 Mass of more massive neutron star.
1167 mass_2: float
1168 Mass of less massive neutron star.
1170 Returns
1171 =======
1172 delta_lambda_tilde: float
1173 Second dominant tidal term.
1174 """
1175 eta = component_masses_to_symmetric_mass_ratio(mass_1, mass_2)
1176 lambda_plus = lambda_1 + lambda_2
1177 lambda_minus = lambda_1 - lambda_2
1178 delta_lambda_tilde = 1 / 2 * (
1179 (1 - 4 * eta) ** 0.5 * (1 - 13272 / 1319 * eta + 8944 / 1319 * eta**2) *
1180 lambda_plus + (1 - 15910 / 1319 * eta + 32850 / 1319 * eta ** 2 +
1181 3380 / 1319 * eta ** 3) * lambda_minus)
1182 return delta_lambda_tilde
1185def lambda_tilde_delta_lambda_tilde_to_lambda_1_lambda_2(
1186 lambda_tilde, delta_lambda_tilde, mass_1, mass_2):
1187 """
1188 Convert from dominant tidal terms to individual tidal parameters.
1190 See, e.g., Wade et al., https://arxiv.org/pdf/1402.5156.pdf.
1192 Parameters
1193 ==========
1194 lambda_tilde: float
1195 Dominant tidal term.
1196 delta_lambda_tilde: float
1197 Secondary tidal term.
1198 mass_1: float
1199 Mass of more massive neutron star.
1200 mass_2: float
1201 Mass of less massive neutron star.
1203 Returns
1204 =======
1205 lambda_1: float
1206 Tidal parameter of more massive neutron star.
1207 lambda_2: float
1208 Tidal parameter of less massive neutron star.
1210 """
1211 eta = component_masses_to_symmetric_mass_ratio(mass_1, mass_2)
1212 coefficient_1 = (1 + 7 * eta - 31 * eta**2)
1213 coefficient_2 = (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta**2)
1214 coefficient_3 = (1 - 4 * eta)**0.5 *\
1215 (1 - 13272 / 1319 * eta + 8944 / 1319 * eta**2)
1216 coefficient_4 = (1 - 15910 / 1319 * eta + 32850 / 1319 * eta**2 +
1217 3380 / 1319 * eta**3)
1218 lambda_1 =\
1219 (13 * lambda_tilde / 8 * (coefficient_3 - coefficient_4) -
1220 2 * delta_lambda_tilde * (coefficient_1 - coefficient_2))\
1221 / ((coefficient_1 + coefficient_2) * (coefficient_3 - coefficient_4) -
1222 (coefficient_1 - coefficient_2) * (coefficient_3 + coefficient_4))
1223 lambda_2 =\
1224 (13 * lambda_tilde / 8 * (coefficient_3 + coefficient_4) -
1225 2 * delta_lambda_tilde * (coefficient_1 + coefficient_2)) \
1226 / ((coefficient_1 - coefficient_2) * (coefficient_3 + coefficient_4) -
1227 (coefficient_1 + coefficient_2) * (coefficient_3 - coefficient_4))
1229 return lambda_1, lambda_2
1232def lambda_tilde_to_lambda_1_lambda_2(
1233 lambda_tilde, mass_1, mass_2):
1234 """
1235 Convert from dominant tidal term to individual tidal parameters
1236 assuming lambda_1 * mass_1**5 = lambda_2 * mass_2**5.
1238 See, e.g., Wade et al., https://arxiv.org/pdf/1402.5156.pdf.
1240 Parameters
1241 ==========
1242 lambda_tilde: float
1243 Dominant tidal term.
1244 mass_1: float
1245 Mass of more massive neutron star.
1246 mass_2: float
1247 Mass of less massive neutron star.
1249 Returns
1250 =======
1251 lambda_1: float
1252 Tidal parameter of more massive neutron star.
1253 lambda_2: float
1254 Tidal parameter of less massive neutron star.
1255 """
1256 eta = component_masses_to_symmetric_mass_ratio(mass_1, mass_2)
1257 q = mass_2 / mass_1
1258 lambda_1 = 13 / 8 * lambda_tilde / (
1259 (1 + 7 * eta - 31 * eta**2) * (1 + q**-5) +
1260 (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta**2) * (1 - q**-5))
1261 lambda_2 = lambda_1 / q**5
1262 return lambda_1, lambda_2
1265def lambda_1_lambda_2_to_lambda_symmetric(lambda_1, lambda_2):
1266 """
1267 Convert from individual tidal parameters to symmetric tidal term.
1269 See, e.g., Yagi & Yunes, https://arxiv.org/abs/1512.02639
1271 Parameters
1272 ==========
1273 lambda_1: float
1274 Tidal parameter of more massive neutron star.
1275 lambda_2: float
1276 Tidal parameter of less massive neutron star.
1278 Returns
1279 ======
1280 lambda_symmetric: float
1281 Symmetric tidal parameter.
1282 """
1283 lambda_symmetric = (lambda_2 + lambda_1) / 2.
1284 return lambda_symmetric
1287def lambda_1_lambda_2_to_lambda_antisymmetric(lambda_1, lambda_2):
1288 """
1289 Convert from individual tidal parameters to antisymmetric tidal term.
1291 See, e.g., Yagi & Yunes, https://arxiv.org/abs/1512.02639
1293 Parameters
1294 ==========
1295 lambda_1: float
1296 Tidal parameter of more massive neutron star.
1297 lambda_2: float
1298 Tidal parameter of less massive neutron star.
1300 Returns
1301 ======
1302 lambda_antisymmetric: float
1303 Antisymmetric tidal parameter.
1304 """
1305 lambda_antisymmetric = (lambda_2 - lambda_1) / 2.
1306 return lambda_antisymmetric
1309def lambda_symmetric_lambda_antisymmetric_to_lambda_1(lambda_symmetric, lambda_antisymmetric):
1310 """
1311 Convert from symmetric and antisymmetric tidal terms to lambda_1
1313 See, e.g., Yagi & Yunes, https://arxiv.org/abs/1512.02639
1315 Parameters
1316 ==========
1317 lambda_symmetric: float
1318 Symmetric tidal parameter.
1319 lambda_antisymmetric: float
1320 Antisymmetric tidal parameter.
1322 Returns
1323 ======
1324 lambda_1: float
1325 Tidal parameter of more massive neutron star.
1326 """
1327 lambda_1 = lambda_symmetric - lambda_antisymmetric
1328 return lambda_1
1331def lambda_symmetric_lambda_antisymmetric_to_lambda_2(lambda_symmetric, lambda_antisymmetric):
1332 """
1333 Convert from symmetric and antisymmetric tidal terms to lambda_2
1335 See, e.g., Yagi & Yunes, https://arxiv.org/abs/1512.02639
1337 Parameters
1338 ==========
1339 lambda_symmetric: float
1340 Symmetric tidal parameter.
1341 lambda_antisymmetric: float
1342 Antisymmetric tidal parameter.
1344 Returns
1345 ======
1346 lambda_2: float
1347 Tidal parameter of less massive neutron star.
1348 """
1349 lambda_2 = lambda_symmetric + lambda_antisymmetric
1350 return lambda_2
1353def lambda_symmetric_lambda_antisymmetric_to_lambda_1_lambda_2(lambda_symmetric, lambda_antisymmetric):
1354 """
1355 Convert from symmetric and antisymmetric tidal terms to lambda_1 and lambda_2
1357 See, e.g., Yagi & Yunes, https://arxiv.org/abs/1512.02639
1359 Parameters
1360 ==========
1361 lambda_symmetric: float
1362 Symmetric tidal parameter.
1363 lambda_antisymmetric: float
1364 Antisymmetric tidal parameter.
1366 Returns
1367 ======
1368 lambda_1: float
1369 Tidal parameter of more massive neutron star.
1370 lambda_2: float
1371 Tidal parameter of less massive neutron star.
1372 """
1373 lambda_1 = lambda_symmetric_lambda_antisymmetric_to_lambda_1(lambda_symmetric, lambda_antisymmetric)
1374 lambda_2 = lambda_symmetric_lambda_antisymmetric_to_lambda_2(lambda_symmetric, lambda_antisymmetric)
1375 return lambda_1, lambda_2
1378def binary_love_fit_lambda_symmetric_mass_ratio_to_lambda_antisymmetric(lambda_symmetric, mass_ratio):
1380 """
1381 Convert from symmetric tidal terms and mass ratio to antisymmetric tidal terms
1382 using BinaryLove relations.
1383 This function does only the fit itself, and doesn't account for marginalisation
1384 over the uncertanties in the fit
1386 See Yagi & Yunes, https://arxiv.org/abs/1512.02639
1387 and Chatziioannou, Haster, Zimmerman, https://arxiv.org/abs/1804.03221
1389 This function will use the implementation from the CHZ paper, which
1390 marginalises over the uncertainties in the BinaryLove fits.
1391 This is also implemented in LALInference through the function
1392 LALInferenceBinaryLove.
1394 Parameters
1395 ==========
1396 lambda_symmetric: float
1397 Symmetric tidal parameter.
1398 mass_ratio: float
1399 Mass ratio (mass_2/mass_1, with mass_2 < mass_1) of the binary
1401 Returns
1402 ======
1403 lambda_antisymmetric: float
1404 Antisymmetric tidal parameter.
1405 """
1406 lambda_symmetric_m1o5 = np.power(lambda_symmetric, -1. / 5.)
1407 lambda_symmetric_m2o5 = lambda_symmetric_m1o5 * lambda_symmetric_m1o5
1408 lambda_symmetric_m3o5 = lambda_symmetric_m2o5 * lambda_symmetric_m1o5
1410 q = mass_ratio
1411 q2 = np.square(mass_ratio)
1413 # Eqn.2 from CHZ, incorporating the dependence on mass ratio
1415 n_polytropic = 0.743 # average polytropic index for the EoSs included in the fit
1416 q_for_Fnofq = np.power(q, 10. / (3. - n_polytropic))
1417 Fnofq = (1. - q_for_Fnofq) / (1. + q_for_Fnofq)
1419 # b_ij and c_ij coefficients are given in Table I of CHZ
1421 b11 = -27.7408
1422 b12 = 8.42358
1423 b21 = 122.686
1424 b22 = -19.7551
1425 b31 = -175.496
1426 b32 = 133.708
1428 c11 = -25.5593
1429 c12 = 5.58527
1430 c21 = 92.0337
1431 c22 = 26.8586
1432 c31 = -70.247
1433 c32 = -56.3076
1435 # Eqn 1 from CHZ, giving the lambda_antisymmetric_fitOnly
1436 # not yet accounting for the uncertainty in the fit
1438 numerator = 1.0 + \
1439 (b11 * q * lambda_symmetric_m1o5) + (b12 * q2 * lambda_symmetric_m1o5) + \
1440 (b21 * q * lambda_symmetric_m2o5) + (b22 * q2 * lambda_symmetric_m2o5) + \
1441 (b31 * q * lambda_symmetric_m3o5) + (b32 * q2 * lambda_symmetric_m3o5)
1443 denominator = 1.0 + \
1444 (c11 * q * lambda_symmetric_m1o5) + (c12 * q2 * lambda_symmetric_m1o5) + \
1445 (c21 * q * lambda_symmetric_m2o5) + (c22 * q2 * lambda_symmetric_m2o5) + \
1446 (c31 * q * lambda_symmetric_m3o5) + (c32 * q2 * lambda_symmetric_m3o5)
1448 lambda_antisymmetric_fitOnly = Fnofq * lambda_symmetric * numerator / denominator
1450 return lambda_antisymmetric_fitOnly
1453def binary_love_lambda_symmetric_to_lambda_1_lambda_2_manual_marginalisation(binary_love_uniform,
1454 lambda_symmetric, mass_ratio):
1455 """
1456 Convert from symmetric tidal terms to lambda_1 and lambda_2
1457 using BinaryLove relations
1459 See Yagi & Yunes, https://arxiv.org/abs/1512.02639
1460 and Chatziioannou, Haster, Zimmerman, https://arxiv.org/abs/1804.03221
1462 This function will use the implementation from the CHZ paper, which
1463 marginalises over the uncertainties in the BinaryLove fits.
1464 This is also implemented in LALInference through the function
1465 LALInferenceBinaryLove.
1467 Parameters
1468 ==========
1469 binary_love_uniform: float (defined in the range [0,1])
1470 Uniformly distributed variable used in BinaryLove uncertainty marginalisation
1471 lambda_symmetric: float
1472 Symmetric tidal parameter.
1473 mass_ratio: float
1474 Mass ratio (mass_2/mass_1, with mass_2 < mass_1) of the binary
1476 Returns
1477 ======
1478 lambda_1: float
1479 Tidal parameter of more massive neutron star.
1480 lambda_2: float
1481 Tidal parameter of less massive neutron star.
1482 """
1483 lambda_antisymmetric_fitOnly = binary_love_fit_lambda_symmetric_mass_ratio_to_lambda_antisymmetric(lambda_symmetric,
1484 mass_ratio)
1486 lambda_symmetric_sqrt = np.sqrt(lambda_symmetric)
1488 q = mass_ratio
1489 q2 = np.square(mass_ratio)
1491 # mu_i and sigma_i coefficients are given in Table II of CHZ
1493 mu_1 = 137.1252739
1494 mu_2 = -32.8026613
1495 mu_3 = 0.5168637
1496 mu_4 = -11.2765281
1497 mu_5 = 14.9499544
1498 mu_6 = - 4.6638851
1500 sigma_1 = -0.0000739
1501 sigma_2 = 0.0103778
1502 sigma_3 = 0.4581717
1503 sigma_4 = -0.8341913
1504 sigma_5 = -201.4323962
1505 sigma_6 = 273.9268276
1506 sigma_7 = -71.2342246
1508 # Eqn 6 from CHZ, correction on fit for lambdaA caused by
1509 # uncertainty in the mean of the lambdaS residual fit,
1510 # using coefficients mu_1, mu_2 and mu_3 from Table II of CHZ
1512 lambda_antisymmetric_lambda_symmetric_meanCorr = \
1513 (mu_1 / (lambda_symmetric * lambda_symmetric)) + \
1514 (mu_2 / lambda_symmetric) + mu_3
1516 # Eqn 8 from CHZ, correction on fit for lambdaA caused by
1517 # uncertainty in the standard deviation of lambdaS residual fit,
1518 # using coefficients sigma_1, sigma_2, sigma_3 and sigma_4 from Table II
1520 lambda_antisymmetric_lambda_symmetric_stdCorr = \
1521 (sigma_1 * lambda_symmetric * lambda_symmetric_sqrt) + \
1522 (sigma_2 * lambda_symmetric) + \
1523 (sigma_3 * lambda_symmetric_sqrt) + sigma_4
1525 # Eqn 7, correction on fit for lambdaA caused by
1526 # uncertainty in the mean of the q residual fit,
1527 # using coefficients mu_4, mu_5 and mu_6 from Table II
1529 lambda_antisymmetric_mass_ratio_meanCorr = \
1530 (mu_4 * q2) + (mu_5 * q) + mu_6
1532 # Eqn 9 from CHZ, correction on fit for lambdaA caused by
1533 # uncertainty in the standard deviation of the q residual fit,
1534 # using coefficients sigma_5, sigma_6 and sigma_7 from Table II
1536 lambda_antisymmetric_mass_ratio_stdCorr = \
1537 (sigma_5 * q2) + (sigma_6 * q) + sigma_7
1539 # Eqn 4 from CHZ, averaging the corrections from the
1540 # mean of the residual fits
1542 lambda_antisymmetric_meanCorr = \
1543 (lambda_antisymmetric_lambda_symmetric_meanCorr +
1544 lambda_antisymmetric_mass_ratio_meanCorr) / 2.
1546 # Eqn 5 from CHZ, averaging the corrections from the
1547 # standard deviations of the residual fits
1549 lambda_antisymmetric_stdCorr = \
1550 np.sqrt(np.square(lambda_antisymmetric_lambda_symmetric_stdCorr) +
1551 np.square(lambda_antisymmetric_mass_ratio_stdCorr))
1553 # Draw a correction on the fit from a
1554 # Gaussian distribution with width lambda_antisymmetric_stdCorr
1555 # this is done by sampling a percent point function (inverse cdf)
1556 # through a U{0,1} variable called binary_love_uniform
1558 lambda_antisymmetric_scatter = norm.ppf(binary_love_uniform, loc=0.,
1559 scale=lambda_antisymmetric_stdCorr)
1561 # Add the correction of the residual mean
1562 # and the Gaussian scatter to the lambda_antisymmetric_fitOnly value
1564 lambda_antisymmetric = lambda_antisymmetric_fitOnly + \
1565 (lambda_antisymmetric_meanCorr + lambda_antisymmetric_scatter)
1567 lambda_1 = lambda_symmetric_lambda_antisymmetric_to_lambda_1(lambda_symmetric, lambda_antisymmetric)
1568 lambda_2 = lambda_symmetric_lambda_antisymmetric_to_lambda_2(lambda_symmetric, lambda_antisymmetric)
1570 # The BinaryLove model is only physically valid where
1571 # lambda_2 > lambda_1 as it assumes mass_1 > mass_2
1572 # It also assumes both lambda1 and lambda2 to be positive
1573 # This is an explicit feature of the "raw" model fit,
1574 # but since this implementation also incorporates
1575 # marginalisation over the fit uncertainty, there can
1576 # be instances where those assumptions are randomly broken.
1577 # For those cases, set lambda_1 and lambda_2 to negative
1578 # values, which in turn will cause (effectively all) the
1579 # waveform models in LALSimulation to fail, thus setting
1580 # logL = -infinity
1582 # if np.greater(lambda_1, lambda_2) or np.less(lambda_1, 0) or np.less(lambda_2, 0):
1583 # lambda_1 = -np.inf
1584 # lambda_2 = -np.inf
1586 # For now set this through an explicit constraint prior instead
1588 return lambda_1, lambda_2
1591def binary_love_lambda_symmetric_to_lambda_1_lambda_2_automatic_marginalisation(lambda_symmetric, mass_ratio):
1592 """
1593 Convert from symmetric tidal terms to lambda_1 and lambda_2
1594 using BinaryLove relations
1596 See Yagi & Yunes, https://arxiv.org/abs/1512.02639
1597 and Chatziioannou, Haster, Zimmerman, https://arxiv.org/abs/1804.03221
1599 This function will use the implementation from the CHZ paper, which
1600 marginalises over the uncertainties in the BinaryLove fits.
1601 This is also implemented in LALInference through the function
1602 LALInferenceBinaryLove.
1604 This function should be used when the BinaryLove marginalisation wasn't
1605 explicitly active in the sampling stage. It will draw a random number in U[0,1]
1606 here instead.
1608 Parameters
1609 ==========
1610 lambda_symmetric: float
1611 Symmetric tidal parameter.
1612 mass_ratio: float
1613 Mass ratio (mass_2/mass_1, with mass_2 < mass_1) of the binary
1615 Returns
1616 ======
1617 lambda_1: float
1618 Tidal parameter of more massive neutron star.
1619 lambda_2: float
1620 Tidal parameter of less massive neutron star.
1621 """
1622 from ..core.utils.random import rng
1624 binary_love_uniform = rng.uniform(0, 1, len(lambda_symmetric))
1626 lambda_1, lambda_2 = binary_love_lambda_symmetric_to_lambda_1_lambda_2_manual_marginalisation(
1627 binary_love_uniform, lambda_symmetric, mass_ratio)
1629 return lambda_1, lambda_2
1632def _generate_all_cbc_parameters(sample, defaults, base_conversion,
1633 likelihood=None, priors=None, npool=1):
1634 """Generate all cbc parameters, helper function for BBH/BNS"""
1635 output_sample = sample.copy()
1636 waveform_defaults = defaults
1637 for key in waveform_defaults:
1638 try:
1639 output_sample[key] = \
1640 likelihood.waveform_generator.waveform_arguments[key]
1641 except (KeyError, AttributeError):
1642 default = waveform_defaults[key]
1643 output_sample[key] = default
1644 logger.debug('Assuming {} = {}'.format(key, default))
1646 output_sample = fill_from_fixed_priors(output_sample, priors)
1647 output_sample, _ = base_conversion(output_sample)
1648 if likelihood is not None:
1649 compute_per_detector_log_likelihoods(
1650 samples=output_sample, likelihood=likelihood, npool=npool)
1652 marginalized_parameters = getattr(likelihood, "_marginalized_parameters", list())
1653 if len(marginalized_parameters) > 0:
1654 try:
1655 generate_posterior_samples_from_marginalized_likelihood(
1656 samples=output_sample, likelihood=likelihood, npool=npool)
1657 except MarginalizedLikelihoodReconstructionError as e:
1658 logger.warning(
1659 "Marginalised parameter reconstruction failed with message "
1660 "{}. Some parameters may not have the intended "
1661 "interpretation.".format(e)
1662 )
1663 if priors is not None:
1664 misnamed_marginalizations = dict(
1665 luminosity_distance="distance",
1666 geocent_time="time",
1667 recalib_index="calibration",
1668 )
1669 for par in marginalized_parameters:
1670 name = misnamed_marginalizations.get(par, par)
1671 if (
1672 getattr(likelihood, f'{name}_marginalization', False)
1673 and par in likelihood.priors
1674 ):
1675 priors[par] = likelihood.priors[par]
1677 if (
1678 not getattr(likelihood, "reference_frame", "sky") == "sky"
1679 or not getattr(likelihood, "time_reference", "geocenter") == "geocenter"
1680 ):
1681 try:
1682 generate_sky_frame_parameters(
1683 samples=output_sample, likelihood=likelihood
1684 )
1685 except TypeError:
1686 logger.info(
1687 "Failed to generate sky frame parameters for type {}"
1688 .format(type(output_sample))
1689 )
1690 if likelihood is not None:
1691 compute_snrs(output_sample, likelihood, npool=npool)
1692 for key, func in zip(["mass", "spin", "source frame"], [
1693 generate_mass_parameters, generate_spin_parameters,
1694 generate_source_frame_parameters]):
1695 try:
1696 output_sample = func(output_sample)
1697 except KeyError as e:
1698 logger.info(
1699 "Generation of {} parameters failed with message {}".format(
1700 key, e))
1701 return output_sample
1704def generate_all_bbh_parameters(sample, likelihood=None, priors=None, npool=1):
1705 """
1706 From either a single sample or a set of samples fill in all missing
1707 BBH parameters, in place.
1709 Parameters
1710 ==========
1711 sample: dict or pandas.DataFrame
1712 Samples to fill in with extra parameters, this may be either an
1713 injection or posterior samples.
1714 likelihood: bilby.gw.likelihood.GravitationalWaveTransient, optional
1715 GravitationalWaveTransient used for sampling, used for waveform and
1716 likelihood.interferometers.
1717 priors: dict, optional
1718 Dictionary of prior objects, used to fill in non-sampled parameters.
1719 """
1720 waveform_defaults = {
1721 'reference_frequency': 50.0, 'waveform_approximant': 'IMRPhenomPv2',
1722 'minimum_frequency': 20.0}
1723 output_sample = _generate_all_cbc_parameters(
1724 sample, defaults=waveform_defaults,
1725 base_conversion=convert_to_lal_binary_black_hole_parameters,
1726 likelihood=likelihood, priors=priors, npool=npool)
1727 return output_sample
1730def generate_all_bns_parameters(sample, likelihood=None, priors=None, npool=1):
1731 """
1732 From either a single sample or a set of samples fill in all missing
1733 BNS parameters, in place.
1735 Since we assume BNS waveforms are aligned, component spins won't be
1736 calculated.
1738 Parameters
1739 ==========
1740 sample: dict or pandas.DataFrame
1741 Samples to fill in with extra parameters, this may be either an
1742 injection or posterior samples.
1743 likelihood: bilby.gw.likelihood.GravitationalWaveTransient, optional
1744 GravitationalWaveTransient used for sampling, used for waveform and
1745 likelihood.interferometers.
1746 priors: dict, optional
1747 Dictionary of prior objects, used to fill in non-sampled parameters.
1748 npool: int, (default=1)
1749 If given, perform generation (where possible) using a multiprocessing pool
1751 """
1752 waveform_defaults = {
1753 'reference_frequency': 50.0, 'waveform_approximant': 'TaylorF2',
1754 'minimum_frequency': 20.0}
1755 output_sample = _generate_all_cbc_parameters(
1756 sample, defaults=waveform_defaults,
1757 base_conversion=convert_to_lal_binary_neutron_star_parameters,
1758 likelihood=likelihood, priors=priors, npool=npool)
1759 try:
1760 output_sample = generate_tidal_parameters(output_sample)
1761 except KeyError as e:
1762 logger.debug(
1763 "Generation of tidal parameters failed with message {}".format(e))
1764 return output_sample
1767def generate_specific_parameters(sample, parameters):
1768 """
1769 Generate a specific subset of parameters that can be generated.
1771 Parameters
1772 ----------
1773 sample: dict
1774 The input sample to be converted.
1775 parameters: list
1776 The list of parameters to return.
1778 Returns
1779 -------
1780 output_sample: dict
1781 The converted parameters
1783 Notes
1784 -----
1785 This is _not_ an optimized function. Under the hood, it generates all
1786 possible parameters and then downselects.
1788 If the passed :code:`parameters` do not include the input parameters,
1789 those will not be returned.
1790 """
1791 updated_sample = generate_all_bns_parameters(sample=sample.copy())
1792 output_sample = sample.__class__()
1793 for key in parameters:
1794 if key in updated_sample:
1795 output_sample[key] = updated_sample[key]
1796 else:
1797 raise KeyError("{} not in converted sample.".format(key))
1798 return output_sample
1801def fill_from_fixed_priors(sample, priors):
1802 """Add parameters with delta function prior to the data frame/dictionary.
1804 Parameters
1805 ==========
1806 sample: dict
1807 A dictionary or data frame
1808 priors: dict
1809 A dictionary of priors
1811 Returns
1812 =======
1813 dict:
1814 """
1815 output_sample = sample.copy()
1816 if priors is not None:
1817 for name in priors:
1818 if isinstance(priors[name], DeltaFunction):
1819 output_sample[name] = priors[name].peak
1820 return output_sample
1823def generate_component_masses(sample, require_add=False, source=False):
1824 """"
1825 Add the component masses to the dataframe/dictionary
1826 We add:
1827 mass_1, mass_2
1828 Or if source=True
1829 mass_1_source, mass_2_source
1830 We also add any other masses which may be necessary for
1831 intermediate steps, i.e. typically the total mass is necessary, along
1832 with the mass ratio, so these will usually be added to the dictionary
1834 If `require_add` is True, then having an incomplete set of mass
1835 parameters (so that the component mass parameters cannot be added)
1836 will throw an error, otherwise it will quietly add nothing to the
1837 dictionary.
1839 Parameters
1840 =========
1841 sample : dict
1842 The input dictionary with at least one
1843 component with overall mass scaling (i.e.
1844 chirp_mass, mass_1, mass_2, total_mass) and
1845 then any other mass parameter.
1846 source : bool, default False
1847 If True, then perform the conversions for source mass parameters
1848 i.e. mass_1_source instead of mass_1
1850 Returns
1851 dict : the updated dictionary
1852 """
1853 def check_and_return_quietly(require_add, sample):
1854 if require_add:
1855 raise KeyError("Insufficient mass parameters in input dictionary")
1856 else:
1857 return sample
1858 output_sample = sample.copy()
1860 if source:
1861 mass_1_key = "mass_1_source"
1862 mass_2_key = "mass_2_source"
1863 total_mass_key = "total_mass_source"
1864 chirp_mass_key = "chirp_mass_source"
1865 else:
1866 mass_1_key = "mass_1"
1867 mass_2_key = "mass_2"
1868 total_mass_key = "total_mass"
1869 chirp_mass_key = "chirp_mass"
1871 if mass_1_key in sample.keys():
1872 if mass_2_key in sample.keys():
1873 return output_sample
1874 if total_mass_key in sample.keys():
1875 output_sample[mass_2_key] = output_sample[total_mass_key] - (
1876 output_sample[mass_1_key]
1877 )
1878 return output_sample
1880 elif "mass_ratio" in sample.keys():
1881 pass
1882 elif "symmetric_mass_ratio" in sample.keys():
1883 output_sample["mass_ratio"] = (
1884 symmetric_mass_ratio_to_mass_ratio(
1885 output_sample["symmetric_mass_ratio"])
1886 )
1887 elif chirp_mass_key in sample.keys():
1888 output_sample["mass_ratio"] = (
1889 mass_1_and_chirp_mass_to_mass_ratio(
1890 mass_1=output_sample[mass_1_key],
1891 chirp_mass=output_sample[chirp_mass_key])
1892 )
1893 else:
1894 return check_and_return_quietly(require_add, sample)
1896 output_sample[mass_2_key] = (
1897 output_sample["mass_ratio"] * output_sample[mass_1_key]
1898 )
1900 return output_sample
1902 elif mass_2_key in sample.keys():
1903 # mass_1 is not in the dict
1904 if total_mass_key in sample.keys():
1905 output_sample[mass_1_key] = (
1906 output_sample[total_mass_key] - output_sample[mass_2_key]
1907 )
1908 return output_sample
1909 elif "mass_ratio" in sample.keys():
1910 pass
1911 elif "symmetric_mass_ratio" in sample.keys():
1912 output_sample["mass_ratio"] = (
1913 symmetric_mass_ratio_to_mass_ratio(
1914 output_sample["symmetric_mass_ratio"])
1915 )
1916 elif chirp_mass_key in sample.keys():
1917 output_sample["mass_ratio"] = (
1918 mass_2_and_chirp_mass_to_mass_ratio(
1919 mass_2=output_sample[mass_2_key],
1920 chirp_mass=output_sample[chirp_mass_key])
1921 )
1922 else:
1923 check_and_return_quietly(require_add, sample)
1925 output_sample[mass_1_key] = 1 / output_sample["mass_ratio"] * (
1926 output_sample[mass_2_key]
1927 )
1929 return output_sample
1931 # Only if neither mass_1 or mass_2 is in the input sample
1932 if total_mass_key in sample.keys():
1933 if "mass_ratio" in sample.keys():
1934 pass # We have everything we need already
1935 elif "symmetric_mass_ratio" in sample.keys():
1936 output_sample["mass_ratio"] = (
1937 symmetric_mass_ratio_to_mass_ratio(
1938 output_sample["symmetric_mass_ratio"])
1939 )
1940 elif chirp_mass_key in sample.keys():
1941 output_sample["symmetric_mass_ratio"] = (
1942 chirp_mass_and_total_mass_to_symmetric_mass_ratio(
1943 chirp_mass=output_sample[chirp_mass_key],
1944 total_mass=output_sample[total_mass_key])
1945 )
1946 output_sample["mass_ratio"] = (
1947 symmetric_mass_ratio_to_mass_ratio(
1948 output_sample["symmetric_mass_ratio"])
1949 )
1950 else:
1951 return check_and_return_quietly(require_add, sample)
1953 elif chirp_mass_key in sample.keys():
1954 if "mass_ratio" in sample.keys():
1955 pass
1956 elif "symmetric_mass_ratio" in sample.keys():
1957 output_sample["mass_ratio"] = (
1958 symmetric_mass_ratio_to_mass_ratio(
1959 sample["symmetric_mass_ratio"])
1960 )
1961 else:
1962 return check_and_return_quietly(require_add, sample)
1964 output_sample[total_mass_key] = (
1965 chirp_mass_and_mass_ratio_to_total_mass(
1966 chirp_mass=output_sample[chirp_mass_key],
1967 mass_ratio=output_sample["mass_ratio"])
1968 )
1970 # We haven't matched any of the criteria
1971 if total_mass_key not in output_sample.keys() or (
1972 "mass_ratio" not in output_sample.keys()):
1973 return check_and_return_quietly(require_add, sample)
1974 mass_1, mass_2 = (
1975 total_mass_and_mass_ratio_to_component_masses(
1976 total_mass=output_sample[total_mass_key],
1977 mass_ratio=output_sample["mass_ratio"])
1978 )
1979 output_sample[mass_1_key] = mass_1
1980 output_sample[mass_2_key] = mass_2
1982 return output_sample
1985def generate_mass_parameters(sample, source=False):
1986 """
1987 Add the known mass parameters to the data frame/dictionary. We do
1988 not recompute keys already present in the dictionary
1990 We add, potentially:
1991 chirp_mass, total_mass, symmetric_mass_ratio, mass_ratio, mass_1, mass_2
1992 Or if source=True:
1993 chirp_mass_source, total_mass_source, symmetric_mass_ratio, mass_ratio, mass_1_source, mass_2_source
1995 Parameters
1996 ==========
1997 sample : dict
1998 The input dictionary with two "spanning" mass parameters
1999 e.g. (mass_1, mass_2), or (chirp_mass, mass_ratio), but not e.g. only
2000 (mass_ratio, symmetric_mass_ratio)
2001 source : bool, default False
2002 If True, then perform the conversions for source mass parameters
2003 i.e. mass_1_source instead of mass_1
2005 Returns
2006 =======
2007 dict: The updated dictionary
2009 """
2010 # Only add the parameters if they're not already present
2011 intermediate_sample = generate_component_masses(sample, source=source)
2012 output_sample = intermediate_sample.copy()
2014 if source:
2015 mass_1_key = 'mass_1_source'
2016 mass_2_key = 'mass_2_source'
2017 total_mass_key = 'total_mass_source'
2018 chirp_mass_key = 'chirp_mass_source'
2019 else:
2020 mass_1_key = 'mass_1'
2021 mass_2_key = 'mass_2'
2022 total_mass_key = 'total_mass'
2023 chirp_mass_key = 'chirp_mass'
2025 if chirp_mass_key not in output_sample.keys():
2026 output_sample[chirp_mass_key] = (
2027 component_masses_to_chirp_mass(output_sample[mass_1_key],
2028 output_sample[mass_2_key])
2029 )
2030 if total_mass_key not in output_sample.keys():
2031 output_sample[total_mass_key] = (
2032 component_masses_to_total_mass(output_sample[mass_1_key],
2033 output_sample[mass_2_key])
2034 )
2035 if 'symmetric_mass_ratio' not in output_sample.keys():
2036 output_sample['symmetric_mass_ratio'] = (
2037 component_masses_to_symmetric_mass_ratio(output_sample[mass_1_key],
2038 output_sample[mass_2_key])
2039 )
2040 if 'mass_ratio' not in output_sample.keys():
2041 output_sample['mass_ratio'] = (
2042 component_masses_to_mass_ratio(output_sample[mass_1_key],
2043 output_sample[mass_2_key])
2044 )
2046 return output_sample
2049def generate_spin_parameters(sample):
2050 """
2051 Add all spin parameters to the data frame/dictionary.
2053 We add:
2054 cartesian spin components, chi_eff, chi_p cos tilt 1, cos tilt 2
2056 Parameters
2057 ==========
2058 sample : dict, pandas.DataFrame
2059 The input dictionary with some spin parameters
2061 Returns
2062 =======
2063 dict: The updated dictionary
2065 """
2066 output_sample = sample.copy()
2068 output_sample = generate_component_spins(output_sample)
2070 output_sample['chi_eff'] = (output_sample['spin_1z'] +
2071 output_sample['spin_2z'] *
2072 output_sample['mass_ratio']) /\
2073 (1 + output_sample['mass_ratio'])
2075 output_sample['chi_1_in_plane'] = np.sqrt(
2076 output_sample['spin_1x'] ** 2 + output_sample['spin_1y'] ** 2
2077 )
2078 output_sample['chi_2_in_plane'] = np.sqrt(
2079 output_sample['spin_2x'] ** 2 + output_sample['spin_2y'] ** 2
2080 )
2082 output_sample['chi_p'] = np.maximum(
2083 output_sample['chi_1_in_plane'],
2084 (4 * output_sample['mass_ratio'] + 3) /
2085 (3 * output_sample['mass_ratio'] + 4) * output_sample['mass_ratio'] *
2086 output_sample['chi_2_in_plane'])
2088 try:
2089 output_sample['cos_tilt_1'] = np.cos(output_sample['tilt_1'])
2090 output_sample['cos_tilt_2'] = np.cos(output_sample['tilt_2'])
2091 except KeyError:
2092 pass
2094 return output_sample
2097def generate_component_spins(sample):
2098 """
2099 Add the component spins to the data frame/dictionary.
2101 This function uses a lalsimulation function to transform the spins.
2103 Parameters
2104 ==========
2105 sample: A dictionary with the necessary spin conversion parameters:
2106 'theta_jn', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1',
2107 'mass_2', 'reference_frequency', 'phase'
2109 Returns
2110 =======
2111 dict: The updated dictionary
2113 """
2114 output_sample = sample.copy()
2115 spin_conversion_parameters =\
2116 ['theta_jn', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2',
2117 'mass_1', 'mass_2', 'reference_frequency', 'phase']
2118 if all(key in output_sample.keys() for key in spin_conversion_parameters):
2119 (
2120 output_sample['iota'], output_sample['spin_1x'],
2121 output_sample['spin_1y'], output_sample['spin_1z'],
2122 output_sample['spin_2x'], output_sample['spin_2y'],
2123 output_sample['spin_2z']
2124 ) = np.vectorize(bilby_to_lalsimulation_spins)(
2125 output_sample['theta_jn'], output_sample['phi_jl'],
2126 output_sample['tilt_1'], output_sample['tilt_2'],
2127 output_sample['phi_12'], output_sample['a_1'], output_sample['a_2'],
2128 output_sample['mass_1'] * solar_mass,
2129 output_sample['mass_2'] * solar_mass,
2130 output_sample['reference_frequency'], output_sample['phase']
2131 )
2133 output_sample['phi_1'] =\
2134 np.fmod(2 * np.pi + np.arctan2(
2135 output_sample['spin_1y'], output_sample['spin_1x']), 2 * np.pi)
2136 output_sample['phi_2'] =\
2137 np.fmod(2 * np.pi + np.arctan2(
2138 output_sample['spin_2y'], output_sample['spin_2x']), 2 * np.pi)
2140 elif 'chi_1' in output_sample and 'chi_2' in output_sample:
2141 output_sample['spin_1x'] = 0
2142 output_sample['spin_1y'] = 0
2143 output_sample['spin_1z'] = output_sample['chi_1']
2144 output_sample['spin_2x'] = 0
2145 output_sample['spin_2y'] = 0
2146 output_sample['spin_2z'] = output_sample['chi_2']
2147 else:
2148 logger.debug("Component spin extraction failed.")
2150 return output_sample
2153def generate_tidal_parameters(sample):
2154 """
2155 Generate all tidal parameters
2157 lambda_tilde, delta_lambda_tilde
2159 Parameters
2160 ==========
2161 sample: dict, pandas.DataFrame
2162 Should contain lambda_1, lambda_2
2164 Returns
2165 =======
2166 output_sample: dict, pandas.DataFrame
2167 Updated sample
2168 """
2169 output_sample = sample.copy()
2171 output_sample['lambda_tilde'] =\
2172 lambda_1_lambda_2_to_lambda_tilde(
2173 output_sample['lambda_1'], output_sample['lambda_2'],
2174 output_sample['mass_1'], output_sample['mass_2'])
2175 output_sample['delta_lambda_tilde'] = \
2176 lambda_1_lambda_2_to_delta_lambda_tilde(
2177 output_sample['lambda_1'], output_sample['lambda_2'],
2178 output_sample['mass_1'], output_sample['mass_2'])
2180 return output_sample
2183def generate_source_frame_parameters(sample):
2184 """
2185 Generate source frame masses along with redshifts and comoving distance.
2187 Parameters
2188 ==========
2189 sample: dict, pandas.DataFrame
2191 Returns
2192 =======
2193 output_sample: dict, pandas.DataFrame
2194 """
2195 output_sample = sample.copy()
2197 output_sample['redshift'] =\
2198 luminosity_distance_to_redshift(output_sample['luminosity_distance'])
2199 output_sample['comoving_distance'] =\
2200 redshift_to_comoving_distance(output_sample['redshift'])
2202 for key in ['mass_1', 'mass_2', 'chirp_mass', 'total_mass']:
2203 if key in output_sample:
2204 output_sample['{}_source'.format(key)] =\
2205 output_sample[key] / (1 + output_sample['redshift'])
2207 return output_sample
2210def compute_snrs(sample, likelihood, npool=1):
2211 """
2212 Compute the optimal and matched filter snrs of all posterior samples
2213 and print it out.
2215 Parameters
2216 ==========
2217 sample: dict or array_like
2219 likelihood: bilby.gw.likelihood.GravitationalWaveTransient
2220 Likelihood function to be applied on the posterior
2222 """
2223 if likelihood is not None:
2224 if isinstance(sample, dict):
2225 likelihood.parameters.update(sample)
2226 signal_polarizations = likelihood.waveform_generator.frequency_domain_strain(likelihood.parameters.copy())
2227 for ifo in likelihood.interferometers:
2228 per_detector_snr = likelihood.calculate_snrs(signal_polarizations, ifo)
2229 sample['{}_matched_filter_snr'.format(ifo.name)] =\
2230 per_detector_snr.complex_matched_filter_snr
2231 sample['{}_optimal_snr'.format(ifo.name)] = \
2232 per_detector_snr.optimal_snr_squared.real ** 0.5
2233 else:
2234 from tqdm.auto import tqdm
2235 logger.info('Computing SNRs for every sample.')
2237 fill_args = [(ii, row) for ii, row in sample.iterrows()]
2238 if npool > 1:
2239 from ..core.sampler.base_sampler import _initialize_global_variables
2240 pool = multiprocessing.Pool(
2241 processes=npool,
2242 initializer=_initialize_global_variables,
2243 initargs=(likelihood, None, None, False),
2244 )
2245 logger.info(
2246 "Using a pool with size {} for nsamples={}".format(npool, len(sample))
2247 )
2248 new_samples = pool.map(_compute_snrs, tqdm(fill_args, file=sys.stdout))
2249 pool.close()
2250 pool.join()
2251 else:
2252 from ..core.sampler.base_sampler import _sampling_convenience_dump
2253 _sampling_convenience_dump.likelihood = likelihood
2254 new_samples = [_compute_snrs(xx) for xx in tqdm(fill_args, file=sys.stdout)]
2256 for ii, ifo in enumerate(likelihood.interferometers):
2257 snr_updates = dict()
2258 for key in new_samples[0][ii].snrs_as_sample.keys():
2259 snr_updates[f"{ifo.name}_{key}"] = list()
2260 for new_sample in new_samples:
2261 snr_update = new_sample[ii].snrs_as_sample
2262 for key, val in snr_update.items():
2263 snr_updates[f"{ifo.name}_{key}"].append(val)
2264 for k, v in snr_updates.items():
2265 sample[k] = v
2266 else:
2267 logger.debug('Not computing SNRs.')
2270def _compute_snrs(args):
2271 """A wrapper of computing the SNRs to enable multiprocessing"""
2272 from ..core.sampler.base_sampler import _sampling_convenience_dump
2273 likelihood = _sampling_convenience_dump.likelihood
2274 ii, sample = args
2275 sample = dict(sample).copy()
2276 likelihood.parameters.update(sample)
2277 signal_polarizations = likelihood.waveform_generator.frequency_domain_strain(
2278 likelihood.parameters.copy()
2279 )
2280 snrs = list()
2281 for ifo in likelihood.interferometers:
2282 snrs.append(likelihood.calculate_snrs(signal_polarizations, ifo, return_array=False))
2283 return snrs
2286def compute_per_detector_log_likelihoods(samples, likelihood, npool=1, block=10):
2287 """
2288 Calculate the log likelihoods in each detector.
2290 Parameters
2291 ==========
2292 samples: DataFrame
2293 Posterior from run with a marginalised likelihood.
2294 likelihood: bilby.gw.likelihood.GravitationalWaveTransient
2295 Likelihood used during sampling.
2296 npool: int, (default=1)
2297 If given, perform generation (where possible) using a multiprocessing pool
2298 block: int, (default=10)
2299 Size of the blocks to use in multiprocessing
2301 Returns
2302 =======
2303 sample: DataFrame
2304 Returns the posterior with new samples.
2305 """
2306 if likelihood is not None:
2307 if not callable(likelihood.compute_per_detector_log_likelihood):
2308 logger.debug('Not computing per-detector log likelihoods.')
2309 return samples
2311 if isinstance(samples, dict):
2312 likelihood.parameters.update(samples)
2313 samples = likelihood.compute_per_detector_log_likelihood()
2314 return samples
2316 elif not isinstance(samples, DataFrame):
2317 raise ValueError("Unable to handle input samples of type {}".format(type(samples)))
2318 from tqdm.auto import tqdm
2320 logger.info('Computing per-detector log likelihoods.')
2322 # Initialize cache dict
2323 cached_samples_dict = dict()
2325 # Store samples to convert for checking
2326 cached_samples_dict["_samples"] = samples
2328 # Set up the multiprocessing
2329 if npool > 1:
2330 from ..core.sampler.base_sampler import _initialize_global_variables
2331 pool = multiprocessing.Pool(
2332 processes=npool,
2333 initializer=_initialize_global_variables,
2334 initargs=(likelihood, None, None, False),
2335 )
2336 logger.info(
2337 "Using a pool with size {} for nsamples={}"
2338 .format(npool, len(samples))
2339 )
2340 else:
2341 from ..core.sampler.base_sampler import _sampling_convenience_dump
2342 _sampling_convenience_dump.likelihood = likelihood
2343 pool = None
2345 fill_args = [(ii, row) for ii, row in samples.iterrows()]
2346 ii = 0
2347 pbar = tqdm(total=len(samples), file=sys.stdout)
2348 while ii < len(samples):
2349 if ii in cached_samples_dict:
2350 ii += block
2351 pbar.update(block)
2352 continue
2354 if pool is not None:
2355 subset_samples = pool.map(_compute_per_detector_log_likelihoods,
2356 fill_args[ii: ii + block])
2357 else:
2358 subset_samples = [list(_compute_per_detector_log_likelihoods(xx))
2359 for xx in fill_args[ii: ii + block]]
2361 cached_samples_dict[ii] = subset_samples
2363 ii += block
2364 pbar.update(len(subset_samples))
2365 pbar.close()
2367 if pool is not None:
2368 pool.close()
2369 pool.join()
2371 new_samples = np.concatenate(
2372 [np.array(val) for key, val in cached_samples_dict.items() if key != "_samples"]
2373 )
2375 for ii, key in \
2376 enumerate([f'{ifo.name}_log_likelihood' for ifo in likelihood.interferometers]):
2377 samples[key] = new_samples[:, ii]
2379 return samples
2381 else:
2382 logger.debug('Not computing per-detector log likelihoods.')
2385def _compute_per_detector_log_likelihoods(args):
2386 """A wrapper of computing the per-detector log likelihoods to enable multiprocessing"""
2387 from ..core.sampler.base_sampler import _sampling_convenience_dump
2388 likelihood = _sampling_convenience_dump.likelihood
2389 ii, sample = args
2390 sample = dict(sample).copy()
2391 likelihood.parameters.update(dict(sample).copy())
2392 new_sample = likelihood.compute_per_detector_log_likelihood()
2393 return tuple((new_sample[key] for key in
2394 [f'{ifo.name}_log_likelihood' for ifo in likelihood.interferometers]))
2397def generate_posterior_samples_from_marginalized_likelihood(
2398 samples, likelihood, npool=1, block=10, use_cache=True):
2399 """
2400 Reconstruct the distance posterior from a run which used a likelihood which
2401 explicitly marginalised over time/distance/phase.
2403 See Eq. (C29-C32) of https://arxiv.org/abs/1809.02293
2405 Parameters
2406 ==========
2407 samples: DataFrame
2408 Posterior from run with a marginalised likelihood.
2409 likelihood: bilby.gw.likelihood.GravitationalWaveTransient
2410 Likelihood used during sampling.
2411 npool: int, (default=1)
2412 If given, perform generation (where possible) using a multiprocessing pool
2413 block: int, (default=10)
2414 Size of the blocks to use in multiprocessing
2415 use_cache: bool, (default=True)
2416 If true, cache the generation so that reconstuction can begin from the
2417 cache on restart.
2419 Returns
2420 =======
2421 sample: DataFrame
2422 Returns the posterior with new samples.
2423 """
2424 from ..core.utils.random import generate_seeds
2426 marginalized_parameters = getattr(likelihood, "_marginalized_parameters", list())
2427 if len(marginalized_parameters) == 0:
2428 return samples
2430 # pass through a dictionary
2431 if isinstance(samples, dict):
2432 return samples
2433 elif not isinstance(samples, DataFrame):
2434 raise ValueError("Unable to handle input samples of type {}".format(type(samples)))
2435 from tqdm.auto import tqdm
2437 logger.info('Reconstructing marginalised parameters.')
2439 try:
2440 cache_filename = f"{likelihood.outdir}/.{likelihood.label}_generate_posterior_cache.pickle"
2441 except AttributeError:
2442 logger.warning("Likelihood has no outdir and label attribute: caching disabled")
2443 use_cache = False
2445 if use_cache and os.path.exists(cache_filename) and not command_line_args.clean:
2446 try:
2447 with open(cache_filename, "rb") as f:
2448 cached_samples_dict = pickle.load(f)
2449 except EOFError:
2450 logger.warning("Cache file is empty")
2451 cached_samples_dict = None
2453 # Check the samples are identical between the cache and current
2454 if (cached_samples_dict is not None) and (cached_samples_dict["_samples"].equals(samples)):
2455 # Calculate reconstruction percentage and print a log message
2456 nsamples_converted = np.sum(
2457 [len(val) for key, val in cached_samples_dict.items() if key != "_samples"]
2458 )
2459 perc = 100 * nsamples_converted / len(cached_samples_dict["_samples"])
2460 logger.info(f'Using cached reconstruction with {perc:0.1f}% converted.')
2461 else:
2462 logger.info("Cached samples dict out of date, ignoring")
2463 cached_samples_dict = dict(_samples=samples)
2465 else:
2466 # Initialize cache dict
2467 cached_samples_dict = dict()
2469 # Store samples to convert for checking
2470 cached_samples_dict["_samples"] = samples
2472 # Set up the multiprocessing
2473 if npool > 1:
2474 from ..core.sampler.base_sampler import _initialize_global_variables
2475 pool = multiprocessing.Pool(
2476 processes=npool,
2477 initializer=_initialize_global_variables,
2478 initargs=(likelihood, None, None, False),
2479 )
2480 logger.info(
2481 "Using a pool with size {} for nsamples={}"
2482 .format(npool, len(samples))
2483 )
2484 else:
2485 from ..core.sampler.base_sampler import _sampling_convenience_dump
2486 _sampling_convenience_dump.likelihood = likelihood
2487 pool = None
2489 seeds = generate_seeds(len(samples))
2490 fill_args = [(ii, row, seed) for (ii, row), seed in zip(samples.iterrows(), seeds)]
2491 ii = 0
2492 pbar = tqdm(total=len(samples), file=sys.stdout)
2493 while ii < len(samples):
2494 if ii in cached_samples_dict:
2495 ii += block
2496 pbar.update(block)
2497 continue
2499 if pool is not None:
2500 subset_samples = pool.map(fill_sample, fill_args[ii: ii + block])
2501 else:
2502 subset_samples = [list(fill_sample(xx)) for xx in fill_args[ii: ii + block]]
2504 cached_samples_dict[ii] = subset_samples
2506 if use_cache:
2507 safe_file_dump(cached_samples_dict, cache_filename, "pickle")
2509 ii += block
2510 pbar.update(len(subset_samples))
2511 pbar.close()
2513 if pool is not None:
2514 pool.close()
2515 pool.join()
2517 new_samples = np.concatenate(
2518 [np.array(val) for key, val in cached_samples_dict.items() if key != "_samples"]
2519 )
2521 for ii, key in enumerate(marginalized_parameters):
2522 samples[key] = new_samples[:, ii]
2524 return samples
2527def generate_sky_frame_parameters(samples, likelihood):
2528 if isinstance(samples, dict):
2529 likelihood.parameters.update(samples)
2530 samples.update(likelihood.get_sky_frame_parameters())
2531 return
2532 elif not isinstance(samples, DataFrame):
2533 raise ValueError
2534 from tqdm.auto import tqdm
2536 logger.info('Generating sky frame parameters.')
2537 new_samples = list()
2538 for ii in tqdm(range(len(samples)), file=sys.stdout):
2539 sample = dict(samples.iloc[ii]).copy()
2540 likelihood.parameters.update(sample)
2541 new_samples.append(likelihood.get_sky_frame_parameters())
2542 new_samples = DataFrame(new_samples)
2543 for key in new_samples:
2544 samples[key] = new_samples[key]
2547def fill_sample(args):
2548 from ..core.sampler.base_sampler import _sampling_convenience_dump
2549 from ..core.utils.random import seed
2551 ii, sample, rseed = args
2552 seed(rseed)
2553 likelihood = _sampling_convenience_dump.likelihood
2554 marginalized_parameters = getattr(likelihood, "_marginalized_parameters", list())
2555 sample = dict(sample).copy()
2556 likelihood.parameters.update(dict(sample).copy())
2557 new_sample = likelihood.generate_posterior_sample_from_marginalized_likelihood()
2558 return tuple((new_sample[key] for key in marginalized_parameters))
2561def identity_map_conversion(parameters):
2562 """An identity map conversion function that makes no changes to the parameters,
2563 but returns the correct signature expected by other conversion functions
2564 (e.g. convert_to_lal_binary_black_hole_parameters)"""
2565 return parameters, []
2568def identity_map_generation(sample, likelihood=None, priors=None, npool=1):
2569 """An identity map generation function that handles marginalizations, SNRs, etc. correctly,
2570 but does not attempt e.g. conversions in mass or spins
2572 Parameters
2573 ==========
2574 sample: dict or pandas.DataFrame
2575 Samples to fill in with extra parameters, this may be either an
2576 injection or posterior samples.
2577 likelihood: bilby.gw.likelihood.GravitationalWaveTransient, optional
2578 GravitationalWaveTransient used for sampling, used for waveform and
2579 likelihood.interferometers.
2580 priors: dict, optional
2581 Dictionary of prior objects, used to fill in non-sampled parameters.
2583 Returns
2584 =======
2586 """
2587 output_sample = sample.copy()
2589 output_sample = fill_from_fixed_priors(output_sample, priors)
2591 if likelihood is not None:
2592 compute_per_detector_log_likelihoods(
2593 samples=output_sample, likelihood=likelihood, npool=npool)
2595 marginalized_parameters = getattr(likelihood, "_marginalized_parameters", list())
2596 if len(marginalized_parameters) > 0:
2597 try:
2598 generate_posterior_samples_from_marginalized_likelihood(
2599 samples=output_sample, likelihood=likelihood, npool=npool)
2600 except MarginalizedLikelihoodReconstructionError as e:
2601 logger.warning(
2602 "Marginalised parameter reconstruction failed with message "
2603 "{}. Some parameters may not have the intended "
2604 "interpretation.".format(e)
2605 )
2607 if ("ra" in output_sample.keys() and "dec" in output_sample.keys() and "psi" in output_sample.keys()):
2608 compute_snrs(output_sample, likelihood, npool=npool)
2609 else:
2610 logger.info(
2611 "Skipping SNR computation since samples have insufficient sky location information"
2612 )
2614 return output_sample