Coverage for pesummary/gw/conversions/__init__.py: 58.9%
1210 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-12-09 22:34 +0000
1# Licensed under an MIT style license -- see LICENSE.md
3error_msg = (
4 "Unable to install '{}'. You will not be able to use some of the inbuilt "
5 "functions."
6)
7import copy
8import numpy as np
9from pathlib import Path
11from pesummary import conf
12from pesummary.utils.decorators import set_docstring
13from pesummary.utils.exceptions import EvolveSpinError
14from pesummary.utils.utils import logger
16try:
17 import lalsimulation
18except ImportError:
19 logger.warning(error_msg.format("lalsimulation"))
20try:
21 import astropy
22except ImportError:
23 logger.warning(error_msg.format("astropy"))
25from .angles import *
26from .cosmology import *
27from .cosmology import _source_from_detector
28from .mass import *
29from .remnant import *
30from .remnant import _final_from_initial_BBH
31from .snr import *
32from .snr import _ifo_snr
33from .spins import *
34from .tidal import *
35from .tidal import _check_NSBH_approximant
36from .time import *
38__author__ = ["Charlie Hoy <charlie.hoy@ligo.org>"]
39_conversion_doc = """
40 Class to calculate all possible derived quantities
42 Parameters
43 ----------
44 data: dict, list
45 either a dictionary or samples or a list of parameters and a list of
46 samples. See the examples below for details
47 extra_kwargs: dict, optional
48 dictionary of kwargs associated with this set of posterior samples.
49 f_low: float, optional
50 the low frequency cut-off to use when evolving the spins
51 f_start: float, optional
52 the starting frequency of the waveform
53 f_ref: float, optional
54 the reference frequency when spins are defined
55 f_final: float, optional
56 the final frequency to use when integrating over frequencies
57 approximant: str, optional
58 the approximant to use when evolving the spins
59 approximant_flags: dict, optional
60 dictionary of flags to control the variant of waveform model used
61 evolve_spins_forwards: float/str, optional
62 the final velocity to evolve the spins up to.
63 evolve_spins_backwards: str, optional
64 method to use when evolving the spins backwards to an infinite separation
65 return_kwargs: Bool, optional
66 if True, return a modified dictionary of kwargs containing information
67 about the conversion
68 NRSur_fits: float/str, optional
69 the NRSurrogate model to use to calculate the remnant fits. If nothing
70 passed, the average NR fits are used instead
71 multipole_snr: Bool, optional
72 if True, the SNR in the (l, m) = [(2, 1), (3, 3), (4, 4)] multipoles are
73 calculated from the posterior samples.
74 samples.
75 precessing_snr: Bool, optional
76 if True, the precessing SNR is calculated from the posterior samples.
77 psd: dict, optional
78 dictionary containing a psd frequency series for each detector you wish
79 to include in calculations
80 waveform_fits: Bool, optional
81 if True, the approximant is used to calculate the remnant fits. Default
82 is False which means that the average NR fits are used
83 multi_process: int, optional
84 number of cores to use to parallelize the computationally expensive
85 conversions
86 redshift_method: str, optional
87 method you wish to use when calculating the redshift given luminosity
88 distance samples. If redshift samples already exist, this method is not
89 used. Default is 'approx' meaning that interpolation is used to calculate
90 the redshift given N luminosity distance points.
91 cosmology: str, optional
92 cosmology you wish to use when calculating the redshift given luminosity
93 distance samples.
94 force_non_evolved: Bool, optional
95 force non evolved remnant quantities to be calculated when evolved quantities
96 already exist in the input. Default False
97 force_BBH_remnant_computation: Bool, optional
98 force BBH remnant quantities to be calculated for systems that include
99 tidal deformability parameters where BBH fits may not be applicable.
100 Default False.
101 force_BH_spin_evolution: Bool, optional
102 force BH spin evolution methods to be applied for systems that include
103 tidal deformability parameters where these methods may not be applicable.
104 Default False.
105 disable_remnant: Bool, optional
106 disable all remnant quantities from being calculated. Default False.
107 add_zero_spin: Bool, optional
108 if no spins are present in the posterior table, add spins with 0 value.
109 Default False.
110 psd_default: str/pycbc.psd obj, optional
111 Default PSD to use for conversions when no other PSD is provided.
112 regenerate: list, optional
113 list of posterior distributions that you wish to regenerate
114 return_dict: Bool, optional
115 if True, return a pesummary.utils.utils.SamplesDict object
116 resume_file: str, optional
117 path to file to use for checkpointing. If not provided, checkpointing
118 is not used. Default None
120 Examples
121 --------
122 There are two ways of passing arguments to this conversion class, either
123 a dictionary of samples or a list of parameters and a list of samples. See
124 the examples below:
126 >>> samples = {"mass_1": 10, "mass_2": 5}
127 >>> converted_samples = %(function)s(samples)
129 >>> parameters = ["mass_1", "mass_2"]
130 >>> samples = [10, 5]
131 >>> converted_samples = %(function)s(parameters, samples)
133 >>> samples = {"mass_1": [10, 20], "mass_2": [5, 8]}
134 >>> converted_samples = %(function)s(samples)
136 >>> parameters = ["mass_1", "mass_2"]
137 >>> samples = [[10, 5], [20, 8]]
138 """
141@set_docstring(_conversion_doc % {"function": "convert"})
142def convert(*args, restart_from_checkpoint=False, resume_file=None, **kwargs):
143 import os
144 if resume_file is not None:
145 if os.path.isfile(resume_file) and restart_from_checkpoint:
146 return _Conversion.load_current_state(resume_file)
147 logger.info(
148 "Unable to find resume file for conversion. Not restarting from "
149 "checkpoint"
150 )
151 return _Conversion(*args, resume_file=resume_file, **kwargs)
154class _PickledConversion(object):
155 pass
158@set_docstring(_conversion_doc % {"function": "_Conversion"})
159class _Conversion(object):
160 @classmethod
161 def load_current_state(cls, resume_file):
162 """Load current state from a resume file
164 Parameters
165 ----------
166 resume_file: str
167 path to a resume file to restart conversion
168 """
169 from pesummary.io import read
170 logger.info(
171 "Reading checkpoint file: {}".format(resume_file)
172 )
173 state = read(resume_file, checkpoint=True)
174 return cls(
175 state.parameters, state.samples, extra_kwargs=state.extra_kwargs,
176 evolve_spins_forwards=state.evolve_spins_forwards,
177 evolve_spins_backwards=state.evolve_spins_backwards,
178 NRSur_fits=state.NRSurrogate,
179 waveform_fits=state.waveform_fit, multi_process=state.multi_process,
180 redshift_method=state.redshift_method, cosmology=state.cosmology,
181 force_non_evolved=state.force_non_evolved,
182 force_BBH_remnant_computation=state.force_remnant,
183 disable_remnant=state.disable_remnant,
184 add_zero_spin=state.add_zero_spin, regenerate=state.regenerate,
185 return_kwargs=state.return_kwargs, return_dict=state.return_dict,
186 resume_file=state.resume_file
187 )
189 def write_current_state(self):
190 """Write the current state of the conversion class to file
191 """
192 from pesummary.io import write
193 state = _PickledConversion()
194 for key, value in vars(self).items():
195 setattr(state, key, value)
197 _path = Path(self.resume_file)
198 write(
199 state, outdir=_path.parent, file_format="pickle",
200 filename=_path.name, overwrite=True
201 )
202 logger.debug(
203 "Written checkpoint file: {}".format(self.resume_file)
204 )
206 def __new__(cls, *args, **kwargs):
207 from pesummary.utils.samples_dict import SamplesDict
208 from pesummary.utils.parameters import Parameters
210 obj = super(_Conversion, cls).__new__(cls)
211 base_replace = (
212 "'{}': {} already found in the result file. Overwriting with "
213 "the passed {}"
214 )
215 if len(args) > 2:
216 raise ValueError(
217 "The _Conversion module only takes as arguments a dictionary "
218 "of samples or a list of parameters and a list of samples"
219 )
220 elif isinstance(args[0], dict):
221 parameters = Parameters(args[0].keys())
222 samples = np.atleast_2d(
223 np.array([args[0][i] for i in parameters]).T
224 ).tolist()
225 else:
226 if not isinstance(args[0], Parameters):
227 parameters = Parameters(args[0])
228 else:
229 parameters = args[0]
230 samples = args[1]
231 samples = np.atleast_2d(samples).tolist()
232 extra_kwargs = kwargs.get("extra_kwargs", {"sampler": {}, "meta_data": {}})
233 f_low = kwargs.get("f_low", None)
234 f_start = kwargs.get("f_start", None)
235 f_ref = kwargs.get("f_ref", None)
236 f_final = kwargs.get("f_final", None)
237 delta_f = kwargs.get("delta_f", None)
239 for param, value in {"f_final": f_final, "delta_f": delta_f}.items():
240 if value is not None and param in extra_kwargs["meta_data"].keys():
241 logger.warning(
242 base_replace.format(
243 param, extra_kwargs["meta_data"][param], value
244 )
245 )
246 extra_kwargs["meta_data"][param] = value
247 elif value is not None:
248 extra_kwargs["meta_data"][param] = value
249 elif value is None and param in extra_kwargs["meta_data"].keys():
250 logger.debug(
251 "Found {} in input file. Using {}Hz".format(
252 param, extra_kwargs["meta_data"][param]
253 )
254 )
255 else:
256 logger.warning(
257 "Could not find {} in input file and one was not passed "
258 "from the command line. Using {}Hz as default".format(
259 param, getattr(conf, "default_{}".format(param))
260 )
261 )
262 extra_kwargs["meta_data"][param] = getattr(
263 conf, "default_{}".format(param)
264 )
266 approximant = kwargs.get("approximant", None)
267 approximant_flags = kwargs.get("approximant_flags", {})
268 NRSurrogate = kwargs.get("NRSur_fits", False)
269 redshift_method = kwargs.get("redshift_method", "approx")
270 cosmology = kwargs.get("cosmology", "Planck15")
271 force_non_evolved = kwargs.get("force_non_evolved", False)
272 force_remnant = kwargs.get("force_BBH_remnant_computation", False)
273 force_evolve = kwargs.get("force_BH_spin_evolution", False)
274 disable_remnant = kwargs.get("disable_remnant", False)
275 if redshift_method not in ["approx", "exact"]:
276 raise ValueError(
277 "'redshift_method' can either be 'approx' corresponding to "
278 "an approximant method, or 'exact' corresponding to an exact "
279 "method of calculating the redshift"
280 )
281 if isinstance(NRSurrogate, bool) and NRSurrogate:
282 raise ValueError(
283 "'NRSur_fits' must be a string corresponding to the "
284 "NRSurrogate model you wish to use to calculate the remnant "
285 "quantities"
286 )
287 waveform_fits = kwargs.get("waveform_fits", False)
288 evolve_spins_forwards = kwargs.get("evolve_spins_forwards", False)
289 evolve_spins_backwards = kwargs.get("evolve_spins_backwards", False)
290 if disable_remnant and (
291 force_non_evolved or force_remnant
292 or NRSurrogate or waveform_fits or evolve_spins_forwards
293 ):
294 _disable = []
295 if force_non_evolved:
296 _disable.append("force_non_evolved")
297 force_non_evolved = False
298 if force_remnant:
299 _disable.append("force_BBH_remnant_computation")
300 force_remnant = False
301 if NRSurrogate:
302 _disable.append("NRSur_fits")
303 NRSurrogate = False
304 if waveform_fits:
305 _disable.append("waveform_fits")
306 waveform_fits = False
307 if evolve_spins_forwards:
308 _disable.append("evolve_spins_forwards")
309 evolve_spins_forwards = False
310 logger.warning(
311 "Unable to use 'disable_remnant' and {}. Setting "
312 "{} and disabling all remnant quantities from being "
313 "calculated".format(
314 " or ".join(_disable),
315 " and ".join(["{}=False".format(_p) for _p in _disable])
316 )
317 )
318 if NRSurrogate and waveform_fits:
319 raise ValueError(
320 "Unable to use both the NRSurrogate and {} to calculate "
321 "remnant quantities. Please select only one option".format(
322 approximant
323 )
324 )
325 if isinstance(evolve_spins_forwards, bool) and evolve_spins_forwards:
326 raise ValueError(
327 "'evolve_spins_forwards' must be a float, the final velocity to "
328 "evolve the spins up to, or a string, 'ISCO', meaning "
329 "evolve the spins up to the ISCO frequency"
330 )
331 if not evolve_spins_forwards and (NRSurrogate or waveform_fits):
332 if (approximant is not None and "eob" in approximant) or NRSurrogate:
333 logger.warning(
334 "Only evolved spin remnant quantities are returned by the "
335 "{} fits.".format(
336 "NRSurrogate" if NRSurrogate else approximant
337 )
338 )
339 elif evolve_spins_forwards and (NRSurrogate or waveform_fits):
340 if (approximant is not None and "eob" in approximant) or NRSurrogate:
341 logger.warning(
342 "The {} fits already evolve the spins. Therefore "
343 "additional spin evolution will not be performed.".format(
344 "NRSurrogate" if NRSurrogate else approximant
345 )
346 )
347 else:
348 logger.warning(
349 "The {} fits are not applied with spin evolution.".format(
350 approximant
351 )
352 )
353 evolve_spins_forwards = False
355 multipole_snr = kwargs.get("multipole_snr", False)
356 precessing_snr = kwargs.get("precessing_snr", False)
358 for freq, name in zip([f_start, f_low], ["f_start", "f_low"]):
359 if freq is not None and name in extra_kwargs["meta_data"].keys():
360 logger.warning(
361 base_replace.format(
362 name, extra_kwargs["meta_data"][name], freq
363 )
364 )
365 extra_kwargs["meta_data"][name] = freq
366 elif freq is not None:
367 extra_kwargs["meta_data"][name] = freq
368 elif freq is None and name in extra_kwargs["meta_data"].keys():
369 logger.debug(
370 "Found {} in input file. Using {}Hz".format(
371 name, extra_kwargs["meta_data"][name]
372 )
373 )
374 else:
375 logger.warning(
376 "Could not find {} in input file and one was not passed from "
377 "the command line. Using {}Hz as default".format(
378 name, conf.default_flow
379 )
380 )
381 extra_kwargs["meta_data"][name] = conf.default_flow
382 if len(approximant_flags) and "approximant_flags" in extra_kwargs["meta_data"].keys():
383 logger.warning(
384 base_replace.format(
385 "approximant_flags", extra_kwargs["meta_data"]["approximant_flags"],
386 approximant_flags
387 )
388 )
389 extra_kwargs["meta_data"]["approximant_flags"] = approximant_flags
390 elif len(approximant_flags):
391 extra_kwargs["meta_data"]["approximant_flags"] = approximant_flags
392 if approximant is not None and "approximant" in extra_kwargs["meta_data"].keys():
393 logger.warning(
394 base_replace.format(
395 "approximant", extra_kwargs["meta_data"]["approximant"],
396 approximant
397 )
398 )
399 extra_kwargs["meta_data"]["approximant"] = approximant
400 elif approximant is not None:
401 extra_kwargs["meta_data"]["approximant"] = approximant
402 if f_ref is not None and "f_ref" in extra_kwargs["meta_data"].keys():
403 logger.warning(
404 base_replace.format(
405 "f_ref", extra_kwargs["meta_data"]["f_ref"], f_ref
406 )
407 )
408 extra_kwargs["meta_data"]["f_ref"] = f_ref
409 elif f_ref is not None:
410 extra_kwargs["meta_data"]["f_ref"] = f_ref
411 regenerate = kwargs.get("regenerate", None)
412 multi_process = kwargs.get("multi_process", None)
413 if multi_process is not None:
414 multi_process = int(multi_process)
415 psd_default = kwargs.get("psd_default", "aLIGOZeroDetHighPower")
416 psd = kwargs.get("psd", {})
417 if psd is None:
418 psd = {}
419 elif psd is not None and not isinstance(psd, dict):
420 raise ValueError(
421 "'psd' must be a dictionary of frequency series for each detector"
422 )
423 ifos = list(psd.keys())
424 pycbc_psd = copy.deepcopy(psd)
425 if psd != {}:
426 from pesummary.gw.file.psd import PSD
427 if isinstance(psd[ifos[0]], PSD):
428 for ifo in ifos:
429 try:
430 pycbc_psd[ifo] = pycbc_psd[ifo].to_pycbc(
431 extra_kwargs["meta_data"]["f_low"],
432 f_high=extra_kwargs["meta_data"]["f_final"],
433 f_high_override=True
434 )
435 except (ImportError, IndexError, ValueError):
436 pass
437 obj.__init__(
438 parameters, samples, extra_kwargs, evolve_spins_forwards, NRSurrogate,
439 waveform_fits, multi_process, regenerate, redshift_method,
440 cosmology, force_non_evolved, force_remnant,
441 kwargs.get("add_zero_spin", False), disable_remnant,
442 kwargs.get("return_kwargs", False), kwargs.get("return_dict", True),
443 kwargs.get("resume_file", None), multipole_snr, precessing_snr,
444 pycbc_psd, psd_default, evolve_spins_backwards, force_evolve
445 )
446 return_kwargs = kwargs.get("return_kwargs", False)
447 if kwargs.get("return_dict", True) and return_kwargs:
448 return [
449 SamplesDict(obj.parameters, np.array(obj.samples).T),
450 obj.extra_kwargs
451 ]
452 elif kwargs.get("return_dict", True):
453 return SamplesDict(obj.parameters, np.array(obj.samples).T)
454 elif return_kwargs:
455 return obj.parameters, obj.samples, obj.extra_kwargs
456 else:
457 return obj.parameters, obj.samples
459 def __init__(
460 self, parameters, samples, extra_kwargs, evolve_spins_forwards, NRSurrogate,
461 waveform_fits, multi_process, regenerate, redshift_method,
462 cosmology, force_non_evolved, force_remnant, add_zero_spin,
463 disable_remnant, return_kwargs, return_dict, resume_file, multipole_snr,
464 precessing_snr, psd, psd_default, evolve_spins_backwards, force_evolve
465 ):
466 self.parameters = parameters
467 self.samples = samples
468 self.extra_kwargs = extra_kwargs
469 self.evolve_spins_forwards = evolve_spins_forwards
470 self.evolve_spins_backwards = evolve_spins_backwards
471 self.NRSurrogate = NRSurrogate
472 self.waveform_fit = waveform_fits
473 self.multi_process = multi_process
474 self.regenerate = regenerate
475 self.redshift_method = redshift_method
476 self.cosmology = cosmology
477 self.force_non_evolved = force_non_evolved
478 self.force_remnant = force_remnant
479 self.force_evolve = force_evolve
480 self.disable_remnant = disable_remnant
481 self.return_kwargs = return_kwargs
482 self.return_dict = return_dict
483 self.resume_file = resume_file
484 self.multipole_snr = multipole_snr
485 self.precessing_snr = precessing_snr
486 self.psd = psd
487 self.psd_default = psd_default
488 self.non_precessing = False
489 cond1 = any(
490 param in self.parameters for param in
491 conf.precessing_angles + conf.precessing_spins
492 )
493 if not cond1:
494 self.non_precessing = True
495 if "chi_p" in self.parameters:
496 _chi_p = self.specific_parameter_samples(["chi_p"])
497 if not np.any(_chi_p):
498 logger.info(
499 "chi_p = 0 for all samples. Treating this as a "
500 "non-precessing system"
501 )
502 self.non_precessing = True
503 elif all(param in self.parameters for param in conf.precessing_spins):
504 samples = self.specific_parameter_samples(conf.precessing_spins)
505 if not any(np.array(samples).flatten()):
506 logger.info(
507 "in-plane spins are zero for all samples. Treating this as "
508 "a non-precessing system"
509 )
510 cond1 = self.non_precessing and evolve_spins_forwards
511 cond2 = self.non_precessing and evolve_spins_backwards
512 if cond1 or cond2:
513 logger.info(
514 "Spin evolution is trivial for a non-precessing system. No additional "
515 "transformation required."
516 )
517 self.evolve_spins_forwards = False
518 self.evolve_spins_backwards = False
519 if not self.non_precessing and multipole_snr:
520 logger.warning(
521 "The calculation for computing the SNR in subdominant "
522 "multipoles assumes that the system is non-precessing. "
523 "Since precessing samples are provided, this may give incorrect "
524 "results"
525 )
526 if self.non_precessing and precessing_snr:
527 logger.info(
528 "Precessing SNR is 0 for a non-precessing system. No additional "
529 "conversion required."
530 )
531 self.precessing_snr = False
532 self.has_tidal = self._check_for_tidal_parameters()
533 self.NSBH = self._check_for_NSBH_system()
534 self.compute_remnant = not self.disable_remnant
535 if self.has_tidal:
536 if force_evolve and (self.evolve_spins_forwards or self.evolve_spins_backwards):
537 logger.warning(
538 "Posterior samples for tidal deformability found in the "
539 "posterior table. 'force_evolve' provided so using BH spin "
540 "evolution methods for this system. This may not give "
541 "sensible results"
542 )
543 elif self.evolve_spins_forwards or self.evolve_spins_backwards:
544 logger.warning(
545 "Tidal deformability parameters found in the posterior table. "
546 "Skipping spin evolution as current methods are only valid "
547 "for BHs."
548 )
549 self.evolve_spins_forwards = False
550 self.evolve_spins_backwards = False
552 if force_remnant and self.NSBH and self.compute_remnant:
553 logger.warning(
554 "Posterior samples for lambda_2 found in the posterior table "
555 "but unable to find samples for lambda_1. Assuming this "
556 "is an NSBH system. 'force_remnant' provided so using BBH remnant "
557 "fits for this system. This may not give sensible results"
558 )
559 elif self.NSBH and self.compute_remnant:
560 logger.warning(
561 "Posterior samples for lambda_2 found in the posterior table "
562 "but unable to find samples for lambda_1. Applying NSBH "
563 "fits to this system."
564 )
565 self.waveform_fit = True
566 elif force_remnant and self.compute_remnant:
567 logger.warning(
568 "Posterior samples for tidal deformability found in the "
569 "posterior table. Applying BBH remnant fits to this system. "
570 "This may not give sensible results."
571 )
572 elif self.compute_remnant:
573 logger.info(
574 "Skipping remnant calculations as tidal deformability "
575 "parameters found in the posterior table."
576 )
577 self.compute_remnant = False
578 if self.regenerate is not None:
579 for param in self.regenerate:
580 self.remove_posterior(param)
581 self.add_zero_spin = add_zero_spin
582 self.generate_all_posterior_samples()
584 def _check_for_tidal_parameters(self):
585 """Check to see if any tidal parameters are stored in the table
586 """
587 from pesummary.gw.file.standard_names import tidal_params
589 if any(param in self.parameters for param in tidal_params):
590 if not all(_ in self.parameters for _ in ["lambda_1", "lambda_2"]):
591 return True
592 _tidal_posterior = self.specific_parameter_samples(
593 ["lambda_1", "lambda_2"]
594 )
595 if not all(np.any(_) for _ in _tidal_posterior):
596 logger.warning(
597 "Tidal deformability parameters found in the posterior "
598 "table but they are all exactly 0. Assuming this is a BBH "
599 "system."
600 )
601 return False
602 return True
603 return False
605 def _check_for_NSBH_system(self):
606 """Check to see if the posterior samples correspond to an NSBH
607 system
608 """
609 if "lambda_2" in self.parameters and "lambda_1" not in self.parameters:
610 _lambda_2 = self.specific_parameter_samples(["lambda_2"])
611 if not np.any(_lambda_2):
612 logger.warning(
613 "Posterior samples for lambda_2 found in the posterior "
614 "table but they are all exactly 0. Assuming this is a BBH "
615 "system."
616 )
617 return False
618 return True
619 elif "lambda_2" in self.parameters and "lambda_1" in self.parameters:
620 _lambda_1, _lambda_2 = self.specific_parameter_samples(
621 ["lambda_1", "lambda_2"]
622 )
623 if not np.any(_lambda_1) and not np.any(_lambda_2):
624 return False
625 elif not np.any(_lambda_1):
626 logger.warning(
627 "Posterior samples for lambda_1 and lambda_2 found in the "
628 "posterior table but lambda_1 is always exactly 0. "
629 "Assuming this is an NSBH system."
630 )
631 return True
632 return False
634 def remove_posterior(self, parameter):
635 if parameter in self.parameters:
636 logger.info(
637 "Removing the posterior samples for '{}'".format(parameter)
638 )
639 ind = self.parameters.index(parameter)
640 self.parameters.remove(self.parameters[ind])
641 for i in self.samples:
642 del i[ind]
643 else:
644 logger.info(
645 "'{}' is not in the table of posterior samples. Unable to "
646 "remove".format(parameter)
647 )
649 def _specific_parameter_samples(self, param):
650 """Return the samples for a specific parameter
652 Parameters
653 ----------
654 param: str
655 the parameter that you would like to return the samples for
656 """
657 if param == "empty":
658 return np.array(np.zeros(len(self.samples)))
659 ind = self.parameters.index(param)
660 samples = np.array([i[ind] for i in self.samples])
661 return samples
663 def specific_parameter_samples(self, param):
664 """Return the samples for either a list or a single parameter
666 Parameters
667 ----------
668 param: list/str
669 the parameter/parameters that you would like to return the samples
670 for
671 """
672 if type(param) == list:
673 samples = [self._specific_parameter_samples(i) for i in param]
674 else:
675 samples = self._specific_parameter_samples(param)
676 return samples
678 def append_data(self, parameter, samples):
679 """Add a list of samples to the existing samples data object
681 Parameters
682 ----------
683 parameter: str
684 the name of the parameter you would like to append
685 samples: list
686 the list of samples that you would like to append
687 """
688 if parameter not in self.parameters:
689 self.parameters.append(parameter)
690 for num, i in enumerate(self.samples):
691 self.samples[num].append(samples[num])
692 if self.resume_file is not None:
693 self.write_current_state()
695 def _mchirp_from_mchirp_source_z(self):
696 samples = self.specific_parameter_samples(["chirp_mass_source", "redshift"])
697 chirp_mass = mchirp_from_mchirp_source_z(samples[0], samples[1])
698 self.append_data("chirp_mass", chirp_mass)
700 def _q_from_eta(self):
701 samples = self.specific_parameter_samples("symmetric_mass_ratio")
702 mass_ratio = q_from_eta(samples)
703 self.append_data("mass_ratio", mass_ratio)
705 def _q_from_m1_m2(self):
706 samples = self.specific_parameter_samples(["mass_1", "mass_2"])
707 mass_ratio = q_from_m1_m2(samples[0], samples[1])
708 self.append_data("mass_ratio", mass_ratio)
710 def _invert_q(self):
711 ind = self.parameters.index("mass_ratio")
712 for num, i in enumerate(self.samples):
713 self.samples[num][ind] = 1. / self.samples[num][ind]
715 def _invq_from_q(self):
716 samples = self.specific_parameter_samples("mass_ratio")
717 inverted_mass_ratio = invq_from_q(samples)
718 self.append_data("inverted_mass_ratio", inverted_mass_ratio)
720 def _mchirp_from_mtotal_q(self):
721 samples = self.specific_parameter_samples(["total_mass", "mass_ratio"])
722 chirp_mass = mchirp_from_mtotal_q(samples[0], samples[1])
723 self.append_data("chirp_mass", chirp_mass)
725 def _m1_from_mchirp_q(self):
726 samples = self.specific_parameter_samples(["chirp_mass", "mass_ratio"])
727 mass_1 = m1_from_mchirp_q(samples[0], samples[1])
728 self.append_data("mass_1", mass_1)
730 def _m2_from_mchirp_q(self):
731 samples = self.specific_parameter_samples(["chirp_mass", "mass_ratio"])
732 mass_2 = m2_from_mchirp_q(samples[0], samples[1])
733 self.append_data("mass_2", mass_2)
735 def _m1_from_mtotal_q(self):
736 samples = self.specific_parameter_samples(["total_mass", "mass_ratio"])
737 mass_1 = m1_from_mtotal_q(samples[0], samples[1])
738 self.append_data("mass_1", mass_1)
740 def _m2_from_mtotal_q(self):
741 samples = self.specific_parameter_samples(["total_mass", "mass_ratio"])
742 mass_2 = m2_from_mtotal_q(samples[0], samples[1])
743 self.append_data("mass_2", mass_2)
745 def _reference_frequency(self):
746 nsamples = len(self.samples)
747 extra_kwargs = self.extra_kwargs["meta_data"]
748 if extra_kwargs != {} and "f_ref" in list(extra_kwargs.keys()):
749 self.append_data(
750 "reference_frequency", [float(extra_kwargs["f_ref"])] * nsamples
751 )
752 else:
753 logger.warning(
754 "Could not find reference_frequency in input file. Using 20Hz "
755 "as default")
756 self.append_data("reference_frequency", [20.] * nsamples)
758 def _mtotal_from_m1_m2(self):
759 samples = self.specific_parameter_samples(["mass_1", "mass_2"])
760 m_total = m_total_from_m1_m2(samples[0], samples[1])
761 self.append_data("total_mass", m_total)
763 def _mchirp_from_m1_m2(self):
764 samples = self.specific_parameter_samples(["mass_1", "mass_2"])
765 chirp_mass = mchirp_from_m1_m2(samples[0], samples[1])
766 self.append_data("chirp_mass", chirp_mass)
768 def _eta_from_m1_m2(self):
769 samples = self.specific_parameter_samples(["mass_1", "mass_2"])
770 eta = eta_from_m1_m2(samples[0], samples[1])
771 self.append_data("symmetric_mass_ratio", eta)
773 def _phi_12_from_phi1_phi2(self):
774 samples = self.specific_parameter_samples(["phi_1", "phi_2"])
775 phi_12 = phi_12_from_phi1_phi2(samples[0], samples[1])
776 self.append_data("phi_12", phi_12)
778 def _phi1_from_spins(self):
779 samples = self.specific_parameter_samples(["spin_1x", "spin_1y"])
780 phi_1 = phi1_from_spins(samples[0], samples[1])
781 self.append_data("phi_1", phi_1)
783 def _phi2_from_spins(self):
784 samples = self.specific_parameter_samples(["spin_2x", "spin_2y"])
785 phi_2 = phi2_from_spins(samples[0], samples[1])
786 self.append_data("phi_2", phi_2)
788 def _spin_angles(self):
789 _spin_angles = ["theta_jn", "phi_jl", "tilt_1", "tilt_2", "phi_12",
790 "a_1", "a_2"]
791 spin_angles_to_calculate = [
792 i for i in _spin_angles if i not in self.parameters]
793 spin_components = [
794 "mass_1", "mass_2", "iota", "spin_1x", "spin_1y", "spin_1z",
795 "spin_2x", "spin_2y", "spin_2z", "reference_frequency"]
796 samples = self.specific_parameter_samples(spin_components)
797 if "phase" in self.parameters:
798 spin_components.append("phase")
799 samples.append(self.specific_parameter_samples("phase"))
800 else:
801 logger.warning(
802 "Phase it not given, we will be assuming that a "
803 "reference phase of 0 to calculate all the spin angles"
804 )
805 samples.append([0] * len(samples[0]))
806 angles = spin_angles(
807 samples[0], samples[1], samples[2], samples[3], samples[4],
808 samples[5], samples[6], samples[7], samples[8], samples[9],
809 samples[10])
811 for i in spin_angles_to_calculate:
812 ind = _spin_angles.index(i)
813 data = np.array([i[ind] for i in angles])
814 self.append_data(i, data)
816 def _non_precessing_component_spins(self):
817 spins = ["iota", "spin_1x", "spin_1y", "spin_1z", "spin_2x", "spin_2y",
818 "spin_2z"]
819 angles = ["a_1", "a_2", "theta_jn", "tilt_1", "tilt_2"]
820 if all(i in self.parameters for i in angles):
821 samples = self.specific_parameter_samples(angles)
822 cond1 = all(i in [0, np.pi] for i in samples[3])
823 cond2 = all(i in [0, np.pi] for i in samples[4])
824 spins_to_calculate = [
825 i for i in spins if i not in self.parameters]
826 if cond1 and cond1:
827 spin_1x = np.array([0.] * len(samples[0]))
828 spin_1y = np.array([0.] * len(samples[0]))
829 spin_1z = samples[0] * np.cos(samples[3])
830 spin_2x = np.array([0.] * len(samples[0]))
831 spin_2y = np.array([0.] * len(samples[0]))
832 spin_2z = samples[1] * np.cos(samples[4])
833 iota = np.array(samples[2])
834 spin_components = [
835 iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z]
837 for i in spins_to_calculate:
838 ind = spins.index(i)
839 data = spin_components[ind]
840 self.append_data(i, data)
842 def _component_spins(self):
843 spins = ["iota", "spin_1x", "spin_1y", "spin_1z", "spin_2x", "spin_2y",
844 "spin_2z"]
845 spins_to_calculate = [
846 i for i in spins if i not in self.parameters]
847 angles = [
848 "theta_jn", "phi_jl", "tilt_1", "tilt_2", "phi_12", "a_1", "a_2",
849 "mass_1", "mass_2", "reference_frequency"]
850 samples = self.specific_parameter_samples(angles)
851 if "phase" in self.parameters:
852 angles.append("phase")
853 samples.append(self.specific_parameter_samples("phase"))
854 else:
855 logger.warning(
856 "Phase it not given, we will be assuming that a "
857 "reference phase of 0 to calculate all the spin angles"
858 )
859 samples.append([0] * len(samples[0]))
860 spin_components = component_spins(
861 samples[0], samples[1], samples[2], samples[3], samples[4],
862 samples[5], samples[6], samples[7], samples[8], samples[9],
863 samples[10])
865 for i in spins_to_calculate:
866 ind = spins.index(i)
867 data = np.array([i[ind] for i in spin_components])
868 self.append_data(i, data)
870 def _component_spins_from_azimuthal_and_polar_angles(self):
871 spins = ["spin_1x", "spin_1y", "spin_1z", "spin_2x", "spin_2y",
872 "spin_2z"]
873 spins_to_calculate = [
874 i for i in spins if i not in self.parameters]
875 angles = [
876 "a_1", "a_2", "a_1_azimuthal", "a_1_polar", "a_2_azimuthal",
877 "a_2_polar"]
878 samples = self.specific_parameter_samples(angles)
879 spin_components = spin_angles_from_azimuthal_and_polar_angles(
880 samples[0], samples[1], samples[2], samples[3], samples[4],
881 samples[5])
882 for i in spins_to_calculate:
883 ind = spins.index(i)
884 data = np.array([i[ind] for i in spin_components])
885 self.append_data(i, data)
887 def _chi_p(self):
888 parameters = [
889 "mass_1", "mass_2", "spin_1x", "spin_1y", "spin_2x", "spin_2y"]
890 samples = self.specific_parameter_samples(parameters)
891 chi_p_samples = chi_p(
892 samples[0], samples[1], samples[2], samples[3], samples[4],
893 samples[5])
894 self.append_data("chi_p", chi_p_samples)
896 def _chi_p_from_tilts(self, suffix=""):
897 parameters = [
898 "mass_1", "mass_2", "a_1", "tilt_1{}".format(suffix), "a_2",
899 "tilt_2{}".format(suffix)
900 ]
901 samples = self.specific_parameter_samples(parameters)
902 chi_p_samples = chi_p_from_tilts(
903 samples[0], samples[1], samples[2], samples[3], samples[4],
904 samples[5]
905 )
906 self.append_data("chi_p{}".format(suffix), chi_p_samples)
908 def _chi_p_2spin(self):
909 parameters = [
910 "mass_1", "mass_2", "spin_1x", "spin_1y", "spin_2x", "spin_2y"]
911 samples = self.specific_parameter_samples(parameters)
912 chi_p_2spin_samples = chi_p_2spin(
913 samples[0], samples[1], samples[2], samples[3], samples[4],
914 samples[5])
915 self.append_data("chi_p_2spin", chi_p_2spin_samples)
917 def _chi_eff(self, suffix=""):
918 parameters = [
919 "mass_1", "mass_2", "spin_1z{}".format(suffix),
920 "spin_2z{}".format(suffix)
921 ]
922 samples = self.specific_parameter_samples(parameters)
923 chi_eff_samples = chi_eff(
924 samples[0], samples[1], samples[2], samples[3])
925 self.append_data("chi_eff{}".format(suffix), chi_eff_samples)
927 def _aligned_spin_from_magnitude_tilts(
928 self, primary=False, secondary=False, suffix=""
929 ):
930 if primary:
931 parameters = ["a_1", "tilt_1{}".format(suffix)]
932 param_to_add = "spin_1z{}".format(suffix)
933 elif secondary:
934 parameters = ["a_2", "tilt_2{}".format(suffix)]
935 param_to_add = "spin_2z{}".format(suffix)
936 samples = self.specific_parameter_samples(parameters)
937 spin_samples = samples[0] * np.cos(samples[1])
938 self.append_data(param_to_add, spin_samples)
940 def _cos_tilt_1_from_tilt_1(self):
941 samples = self.specific_parameter_samples("tilt_1")
942 cos_tilt_1 = np.cos(samples)
943 self.append_data("cos_tilt_1", cos_tilt_1)
945 def _cos_tilt_2_from_tilt_2(self):
946 samples = self.specific_parameter_samples("tilt_2")
947 cos_tilt_2 = np.cos(samples)
948 self.append_data("cos_tilt_2", cos_tilt_2)
950 def _viewing_angle(self):
951 samples = self.specific_parameter_samples("theta_jn")
952 viewing_angle = viewing_angle_from_inclination(samples)
953 self.append_data("viewing_angle", viewing_angle)
955 def _dL_from_z(self):
956 samples = self.specific_parameter_samples("redshift")
957 distance = dL_from_z(samples, cosmology=self.cosmology)
958 self.extra_kwargs["meta_data"]["cosmology"] = self.cosmology
959 self.append_data("luminosity_distance", distance)
961 def _z_from_dL(self):
962 samples = self.specific_parameter_samples("luminosity_distance")
963 func = getattr(Redshift.Distance, self.redshift_method)
964 redshift = func(
965 samples, cosmology=self.cosmology, multi_process=self.multi_process
966 )
967 self.extra_kwargs["meta_data"]["cosmology"] = self.cosmology
968 self.append_data("redshift", redshift)
970 def _z_from_comoving_volume(self):
971 samples = self.specific_parameter_samples("comoving_volume")
972 func = getattr(Redshift.ComovingVolume, self.redshift_method)
973 redshift = func(
974 samples, cosmology=self.cosmology, multi_process=self.multi_process
975 )
976 self.extra_kwargs["meta_data"]["cosmology"] = self.cosmology
977 self.append_data("redshift", redshift)
979 def _comoving_distance_from_z(self):
980 samples = self.specific_parameter_samples("redshift")
981 distance = comoving_distance_from_z(samples, cosmology=self.cosmology)
982 self.extra_kwargs["meta_data"]["cosmology"] = self.cosmology
983 self.append_data("comoving_distance", distance)
985 def _comoving_volume_from_z(self):
986 samples = self.specific_parameter_samples("redshift")
987 volume = comoving_volume_from_z(samples, cosmology=self.cosmology)
988 self.extra_kwargs["meta_data"]["cosmology"] = self.cosmology
989 self.append_data("comoving_volume", volume)
991 def _m1_source_from_m1_z(self):
992 samples = self.specific_parameter_samples(["mass_1", "redshift"])
993 mass_1_source = m1_source_from_m1_z(samples[0], samples[1])
994 self.append_data("mass_1_source", mass_1_source)
996 def _m1_from_m1_source_z(self):
997 samples = self.specific_parameter_samples(["mass_1_source", "redshift"])
998 mass_1 = m1_from_m1_source_z(samples[0], samples[1])
999 self.append_data("mass_1", mass_1)
1001 def _m2_source_from_m2_z(self):
1002 samples = self.specific_parameter_samples(["mass_2", "redshift"])
1003 mass_2_source = m2_source_from_m2_z(samples[0], samples[1])
1004 self.append_data("mass_2_source", mass_2_source)
1006 def _m2_from_m2_source_z(self):
1007 samples = self.specific_parameter_samples(["mass_2_source", "redshift"])
1008 mass_2 = m2_from_m2_source_z(samples[0], samples[1])
1009 self.append_data("mass_2", mass_2)
1011 def _mtotal_source_from_mtotal_z(self):
1012 samples = self.specific_parameter_samples(["total_mass", "redshift"])
1013 total_mass_source = m_total_source_from_mtotal_z(samples[0], samples[1])
1014 self.append_data("total_mass_source", total_mass_source)
1016 def _mtotal_from_mtotal_source_z(self):
1017 samples = self.specific_parameter_samples(["total_mass_source", "redshift"])
1018 total_mass = mtotal_from_mtotal_source_z(samples[0], samples[1])
1019 self.append_data("total_mass", total_mass)
1021 def _mchirp_source_from_mchirp_z(self):
1022 samples = self.specific_parameter_samples(["chirp_mass", "redshift"])
1023 chirp_mass_source = mchirp_source_from_mchirp_z(samples[0], samples[1])
1024 self.append_data("chirp_mass_source", chirp_mass_source)
1026 def _beta(self):
1027 samples = self.specific_parameter_samples([
1028 "mass_1", "mass_2", "phi_jl", "tilt_1", "tilt_2", "phi_12",
1029 "a_1", "a_2", "reference_frequency", "phase"
1030 ])
1031 beta = opening_angle(
1032 samples[0], samples[1], samples[2], samples[3], samples[4],
1033 samples[5], samples[6], samples[7], samples[8], samples[9]
1034 )
1035 self.append_data("beta", beta)
1037 def _psi_J(self):
1038 samples = self.specific_parameter_samples([
1039 "psi", "theta_jn", "phi_jl", "beta"
1040 ])
1041 psi = psi_J(samples[0], samples[1], samples[2], samples[3])
1042 self.append_data("psi_J", psi)
1044 def _retrieve_detectors(self):
1045 detectors = []
1046 if "IFOs" in list(self.extra_kwargs["meta_data"].keys()):
1047 detectors = self.extra_kwargs["meta_data"]["IFOs"].split(" ")
1048 else:
1049 for i in self.parameters:
1050 if "optimal_snr" in i and i != "network_optimal_snr":
1051 det = i.split("_optimal_snr")[0]
1052 detectors.append(det)
1053 return detectors
1055 def _time_in_each_ifo(self):
1056 detectors = self._retrieve_detectors()
1057 samples = self.specific_parameter_samples(["ra", "dec", "geocent_time"])
1058 for i in detectors:
1059 time = time_in_each_ifo(i, samples[0], samples[1], samples[2])
1060 self.append_data("%s_time" % (i), time)
1062 def _time_delay_between_ifos(self):
1063 from itertools import combinations
1064 detectors = sorted(self._retrieve_detectors())
1065 for ifos in combinations(detectors, 2):
1066 samples = self.specific_parameter_samples(
1067 ["{}_time".format(_ifo) for _ifo in ifos]
1068 )
1069 self.append_data(
1070 "%s_%s_time_delay" % (ifos[0], ifos[1]),
1071 samples[1] - samples[0]
1072 )
1074 def _lambda1_from_lambda_tilde(self):
1075 samples = self.specific_parameter_samples([
1076 "lambda_tilde", "mass_1", "mass_2"])
1077 lambda_1 = lambda1_from_lambda_tilde(samples[0], samples[1], samples[2])
1078 self.append_data("lambda_1", lambda_1)
1080 def _lambda2_from_lambda1(self):
1081 samples = self.specific_parameter_samples([
1082 "lambda_1", "mass_1", "mass_2"])
1083 lambda_2 = lambda2_from_lambda1(samples[0], samples[1], samples[2])
1084 self.append_data("lambda_2", lambda_2)
1086 def _lambda_tilde_from_lambda1_lambda2(self):
1087 samples = self.specific_parameter_samples([
1088 "lambda_1", "lambda_2", "mass_1", "mass_2"])
1089 lambda_tilde = lambda_tilde_from_lambda1_lambda2(
1090 samples[0], samples[1], samples[2], samples[3])
1091 self.append_data("lambda_tilde", lambda_tilde)
1093 def _delta_lambda_from_lambda1_lambda2(self):
1094 samples = self.specific_parameter_samples([
1095 "lambda_1", "lambda_2", "mass_1", "mass_2"])
1096 delta_lambda = delta_lambda_from_lambda1_lambda2(
1097 samples[0], samples[1], samples[2], samples[3])
1098 self.append_data("delta_lambda", delta_lambda)
1100 def _NS_compactness_from_lambda(self, parameter="lambda_1"):
1101 if parameter not in ["lambda_1", "lambda_2"]:
1102 logger.warning(
1103 "Can only use Love-compactness relation for 'lambda_1' and/or "
1104 "'lambda_2'. Skipping conversion"
1105 )
1106 return
1107 ind = parameter.split("lambda_")[1]
1108 samples = self.specific_parameter_samples([parameter])
1109 compactness = NS_compactness_from_lambda(samples[0])
1110 self.append_data("compactness_{}".format(ind), compactness)
1111 self.extra_kwargs["meta_data"]["compactness_fits"] = (
1112 "YagiYunes2017_with_BBHlimit"
1113 )
1115 def _NS_baryonic_mass(self, primary=True):
1116 if primary:
1117 params = ["compactness_1", "mass_1"]
1118 else:
1119 params = ["compactness_2", "mass_2"]
1120 samples = self.specific_parameter_samples(params)
1121 mass = NS_baryonic_mass(samples[0], samples[1])
1122 if primary:
1123 self.append_data("baryonic_mass_1", mass)
1124 else:
1125 self.append_data("baryonic_mass_2", mass)
1126 self.extra_kwargs["meta_data"]["baryonic_mass_fits"] = "Breu2016"
1128 def _lambda1_lambda2_from_polytrope_EOS(self):
1129 samples = self.specific_parameter_samples([
1130 "log_pressure", "gamma_1", "gamma_2", "gamma_3", "mass_1", "mass_2"
1131 ])
1132 lambda_1, lambda_2 = \
1133 lambda1_lambda2_from_4_parameter_piecewise_polytrope_equation_of_state(
1134 *samples, multi_process=self.multi_process
1135 )
1136 if "lambda_1" not in self.parameters:
1137 self.append_data("lambda_1", lambda_1)
1138 if "lambda_2" not in self.parameters:
1139 self.append_data("lambda_2", lambda_2)
1141 def _lambda1_lambda2_from_spectral_decomposition_EOS(self):
1142 samples = self.specific_parameter_samples([
1143 "spectral_decomposition_gamma_0", "spectral_decomposition_gamma_1",
1144 "spectral_decomposition_gamma_2", "spectral_decomposition_gamma_3",
1145 "mass_1", "mass_2"
1146 ])
1147 lambda_1, lambda_2 = lambda1_lambda2_from_spectral_decomposition(
1148 *samples, multi_process=self.multi_process
1149 )
1150 if "lambda_1" not in self.parameters:
1151 self.append_data("lambda_1", lambda_1)
1152 if "lambda_2" not in self.parameters:
1153 self.append_data("lambda_2", lambda_2)
1155 def _ifo_snr(self):
1156 abs_snrs = [
1157 i for i in self.parameters if "_matched_filter_abs_snr" in i
1158 ]
1159 angle_snrs = [
1160 i for i in self.parameters if "_matched_filter_snr_angle" in i
1161 ]
1162 for ifo in [snr.split("_matched_filter_abs_snr")[0] for snr in abs_snrs]:
1163 if "{}_matched_filter_snr".format(ifo) not in self.parameters:
1164 samples = self.specific_parameter_samples(
1165 [
1166 "{}_matched_filter_abs_snr".format(ifo),
1167 "{}_matched_filter_snr_angle".format(ifo)
1168 ]
1169 )
1170 snr = _ifo_snr(samples[0], samples[1])
1171 self.append_data("{}_matched_filter_snr".format(ifo), snr)
1173 def _optimal_network_snr(self):
1174 snrs = [i for i in self.parameters if "_optimal_snr" in i]
1175 samples = self.specific_parameter_samples(snrs)
1176 snr = network_snr(samples)
1177 self.append_data("network_optimal_snr", snr)
1179 def _matched_filter_network_snr(self):
1180 mf_snrs = sorted([
1181 i for i in self.parameters if "_matched_filter_snr" in i
1182 and "_angle" not in i and "_abs" not in i
1183 ])
1184 opt_snrs = sorted([
1185 i for i in self.parameters if "_optimal_snr" in i and "network" not
1186 in i
1187 ])
1188 _mf_detectors = [
1189 param.split("_matched_filter_snr")[0] for param in mf_snrs
1190 ]
1191 _opt_detectors = [
1192 param.split("_optimal_snr")[0] for param in opt_snrs
1193 ]
1194 if _mf_detectors == _opt_detectors:
1195 mf_samples = self.specific_parameter_samples(mf_snrs)
1196 opt_samples = self.specific_parameter_samples(opt_snrs)
1197 snr = network_matched_filter_snr(mf_samples, opt_samples)
1198 self.append_data("network_matched_filter_snr", snr)
1199 else:
1200 logger.warning(
1201 "Unable to generate 'network_matched_filter_snr' as "
1202 "there is an inconsistency in the detector network based on "
1203 "the 'optimal_snrs' and the 'matched_filter_snrs'. We find "
1204 "that from the 'optimal_snrs', the detector network is: {} "
1205 "while we find from the 'matched_filter_snrs', the detector "
1206 "network is: {}".format(_opt_detectors, _mf_detectors)
1207 )
1209 def _rho_hm(self, multipoles):
1210 import copy
1211 required = [
1212 "mass_1", "mass_2", "spin_1z", "spin_2z", "psi", "iota", "ra",
1213 "dec", "geocent_time", "luminosity_distance", "phase",
1214 "reference_frequency"
1215 ]
1216 samples = self.specific_parameter_samples(required)
1217 _f_low = self._retrieve_stored_frequency("f_low")
1218 if isinstance(_f_low, (np.ndarray)):
1219 f_low = _f_low() * len(samples[0])
1220 else:
1221 f_low = [_f_low] * len(samples[0])
1222 original_list = copy.deepcopy(multipoles)
1223 rho, data_used = multipole_snr(
1224 *samples[:-1], f_low=f_low, psd=self.psd, f_ref=samples[-1],
1225 f_final=self.extra_kwargs["meta_data"]["f_final"],
1226 return_data_used=True, multi_process=self.multi_process,
1227 psd_default=self.psd_default, multipole=multipoles,
1228 df=self.extra_kwargs["meta_data"]["delta_f"]
1229 )
1230 for num, mm in enumerate(original_list):
1231 self.append_data("network_{}_multipole_snr".format(mm), rho[num])
1232 self.extra_kwargs["meta_data"]["multipole_snr"] = data_used
1234 def _rho_p(self):
1235 required = [
1236 "mass_1", "mass_2", "beta", "psi_J", "a_1", "a_2", "tilt_1",
1237 "tilt_2", "phi_12", "theta_jn", "ra", "dec", "geocent_time",
1238 "phi_jl", "reference_frequency", "luminosity_distance", "phase"
1239 ]
1240 samples = self.specific_parameter_samples(required)
1241 try:
1242 spins = self.specific_parameter_samples(["spin_1z", "spin_2z"])
1243 except ValueError:
1244 spins = [None, None]
1245 _f_low = self._retrieve_stored_frequency("f_low")
1246 if isinstance(_f_low, (np.ndarray)):
1247 f_low = _f_low() * len(samples[0])
1248 else:
1249 f_low = [_f_low] * len(samples[0])
1250 [rho_p, b_bar, overlap, snrs], data_used = precessing_snr(
1251 samples[0], samples[1], samples[2], samples[3], samples[4],
1252 samples[5], samples[6], samples[7], samples[8], samples[9], samples[10],
1253 samples[11], samples[12], samples[13], samples[15], samples[16], f_low=f_low,
1254 spin_1z=spins[0], spin_2z=spins[1], psd=self.psd, return_data_used=True,
1255 f_final=self.extra_kwargs["meta_data"]["f_final"], f_ref=samples[14],
1256 multi_process=self.multi_process, psd_default=self.psd_default,
1257 df=self.extra_kwargs["meta_data"]["delta_f"], debug=True
1258 )
1259 self.append_data("network_precessing_snr", rho_p)
1260 self.append_data("_b_bar", b_bar)
1261 self.append_data("_precessing_harmonics_overlap", overlap)
1262 nbreakdown = len(np.argwhere(b_bar > 0.3))
1263 if nbreakdown > 0:
1264 logger.warning(
1265 "{}/{} ({}%) samples have b_bar greater than 0.3. For these "
1266 "samples, the two-harmonic approximation used to calculate "
1267 "the precession SNR may not be valid".format(
1268 nbreakdown, len(b_bar),
1269 np.round((nbreakdown / len(b_bar)) * 100, 2)
1270 )
1271 )
1272 try:
1273 _samples = self.specific_parameter_samples("network_optimal_snr")
1274 if np.logical_or(
1275 np.median(snrs) > 1.1 * np.median(_samples),
1276 np.median(snrs) < 0.9 * np.median(_samples)
1277 ):
1278 logger.warning(
1279 "The two-harmonic SNR is different from the stored SNR. "
1280 "This indicates that the provided PSD may be different "
1281 "from the one used in the sampling."
1282 )
1283 except Exception:
1284 pass
1285 self.extra_kwargs["meta_data"]["precessing_snr"] = data_used
1287 def _retrieve_stored_frequency(self, name):
1288 extra_kwargs = self.extra_kwargs["meta_data"]
1289 if extra_kwargs != {} and name in list(extra_kwargs.keys()):
1290 freq = extra_kwargs[name]
1291 else:
1292 raise ValueError(
1293 "Could not find f_low in input file. Please either modify the "
1294 "input file or pass it from the command line"
1295 )
1296 return freq
1298 def _retrieve_approximant(self):
1299 extra_kwargs = self.extra_kwargs["meta_data"]
1300 if extra_kwargs != {} and "approximant" in list(extra_kwargs.keys()):
1301 approximant = extra_kwargs["approximant"]
1302 else:
1303 raise ValueError(
1304 "Unable to find the approximant used to generate the posterior "
1305 "samples in the result file."
1306 )
1307 return approximant
1309 def _evolve_spins(self, final_velocity="ISCO", forward=True):
1310 from .evolve import evolve_spins
1311 from ..waveform import _check_approximant_from_string
1313 samples = self.specific_parameter_samples(
1314 ["mass_1", "mass_2", "a_1", "a_2", "tilt_1", "tilt_2",
1315 "phi_12", "reference_frequency"]
1316 )
1317 approximant = self._retrieve_approximant()
1318 if not _check_approximant_from_string(approximant):
1319 _msg = (
1320 'Not evolving spins: approximant {0} unknown to '
1321 'lalsimulation and gwsignal'.format(approximant)
1322 )
1323 logger.warning(_msg)
1324 raise EvolveSpinError(_msg)
1325 f_low = self._retrieve_stored_frequency("f_low")
1326 if not forward:
1327 [tilt_1_evolved, tilt_2_evolved, phi_12_evolved], fits_used = evolve_spins(
1328 samples[0], samples[1], samples[2], samples[3], samples[4],
1329 samples[5], samples[6], f_low, samples[7][0],
1330 approximant, evolve_limit="infinite_separation",
1331 multi_process=self.multi_process, return_fits_used=True,
1332 method=self.evolve_spins_backwards
1333 )
1334 suffix = ""
1335 if self.evolve_spins_backwards.lower() == "precession_averaged":
1336 suffix = "_only_prec_avg"
1337 self.append_data("tilt_1_infinity{}".format(suffix), tilt_1_evolved)
1338 self.append_data("tilt_2_infinity{}".format(suffix), tilt_2_evolved)
1339 self.extra_kwargs["meta_data"]["backward_spin_evolution"] = fits_used
1340 return
1341 else:
1342 tilt_1_evolved, tilt_2_evolved, phi_12_evolved = evolve_spins(
1343 samples[0], samples[1], samples[2], samples[3], samples[4],
1344 samples[5], samples[6], f_low, samples[7][0],
1345 approximant, final_velocity=final_velocity,
1346 multi_process=self.multi_process
1347 )
1348 self.extra_kwargs["meta_data"]["forward_spin_evolution"] = final_velocity
1349 spin_1z_evolved = samples[2] * np.cos(tilt_1_evolved)
1350 spin_2z_evolved = samples[3] * np.cos(tilt_2_evolved)
1351 self.append_data("tilt_1_evolved", tilt_1_evolved)
1352 self.append_data("tilt_2_evolved", tilt_2_evolved)
1353 self.append_data("phi_12_evolved", phi_12_evolved)
1354 self.append_data("spin_1z_evolved", spin_1z_evolved)
1355 self.append_data("spin_2z_evolved", spin_2z_evolved)
1357 @staticmethod
1358 def _evolved_vs_non_evolved_parameter(
1359 parameter, evolved=False, core_param=False, non_precessing=False
1360 ):
1361 if non_precessing:
1362 base_string = ""
1363 elif evolved and core_param:
1364 base_string = "_evolved"
1365 elif evolved:
1366 base_string = ""
1367 elif core_param:
1368 base_string = ""
1369 else:
1370 base_string = "_non_evolved"
1371 return "{}{}".format(parameter, base_string)
1373 def _precessing_vs_non_precessing_parameters(
1374 self, non_precessing=False, evolved=False
1375 ):
1376 if not non_precessing:
1377 tilt_1 = self._evolved_vs_non_evolved_parameter(
1378 "tilt_1", evolved=evolved, core_param=True
1379 )
1380 tilt_2 = self._evolved_vs_non_evolved_parameter(
1381 "tilt_2", evolved=evolved, core_param=True
1382 )
1383 samples = self.specific_parameter_samples([
1384 "mass_1", "mass_2", "a_1", "a_2", tilt_1, tilt_2
1385 ])
1386 if "phi_12" in self.parameters and evolved:
1387 phi_12_samples = self.specific_parameter_samples([
1388 self._evolved_vs_non_evolved_parameter(
1389 "phi_12", evolved=True, core_param=True
1390 )
1391 ])[0]
1392 elif "phi_12" in self.parameters:
1393 phi_12_samples = self.specific_parameter_samples(["phi_12"])[0]
1394 else:
1395 phi_12_samples = np.zeros_like(samples[0])
1396 samples.append(phi_12_samples)
1397 if self.NRSurrogate:
1398 NRSurrogate_samples = self.specific_parameter_samples([
1399 "phi_jl", "theta_jn", "phase"
1400 ])
1401 for ss in NRSurrogate_samples:
1402 samples.append(ss)
1403 else:
1404 spin_1z = self._evolved_vs_non_evolved_parameter(
1405 "spin_1z", evolved=evolved, core_param=True, non_precessing=True
1406 )
1407 spin_2z = self._evolved_vs_non_evolved_parameter(
1408 "spin_2z", evolved=evolved, core_param=True, non_precessing=True
1409 )
1410 samples = self.specific_parameter_samples([
1411 "mass_1", "mass_2", spin_1z, spin_2z
1412 ])
1413 samples = [
1414 samples[0], samples[1], np.abs(samples[2]), np.abs(samples[3]),
1415 0.5 * np.pi * (1 - np.sign(samples[2])),
1416 0.5 * np.pi * (1 - np.sign(samples[3])),
1417 np.zeros_like(samples[0])
1418 ]
1419 return samples
1421 def _peak_luminosity_of_merger(self, evolved=False):
1422 param = self._evolved_vs_non_evolved_parameter(
1423 "peak_luminosity", evolved=evolved, non_precessing=self.non_precessing
1424 )
1425 spin_1z_param = self._evolved_vs_non_evolved_parameter(
1426 "spin_1z", evolved=evolved, core_param=True,
1427 non_precessing=self.non_precessing
1428 )
1429 spin_2z_param = self._evolved_vs_non_evolved_parameter(
1430 "spin_2z", evolved=evolved, core_param=True,
1431 non_precessing=self.non_precessing
1432 )
1434 samples = self.specific_parameter_samples([
1435 "mass_1", "mass_2", spin_1z_param, spin_2z_param
1436 ])
1437 peak_luminosity, fits = peak_luminosity_of_merger(
1438 samples[0], samples[1], samples[2], samples[3],
1439 return_fits_used=True
1440 )
1441 self.append_data(param, peak_luminosity)
1442 self.extra_kwargs["meta_data"]["peak_luminosity_NR_fits"] = fits
1444 def _final_remnant_properties_from_NRSurrogate(
1445 self, non_precessing=False,
1446 parameters=["final_mass", "final_spin", "final_kick"]
1447 ):
1448 f_low = self._retrieve_stored_frequency("f_start")
1449 approximant = self._retrieve_approximant()
1450 samples = self._precessing_vs_non_precessing_parameters(
1451 non_precessing=non_precessing, evolved=False
1452 )
1453 frequency_samples = self.specific_parameter_samples([
1454 "reference_frequency"
1455 ])
1456 data, fits = final_remnant_properties_from_NRSurrogate(
1457 *samples, f_low=f_low, f_ref=frequency_samples[0],
1458 properties=parameters, return_fits_used=True,
1459 approximant=approximant
1460 )
1461 for param in parameters:
1462 self.append_data(param, data[param])
1463 self.extra_kwargs["meta_data"]["{}_NR_fits".format(param)] = fits
1465 def _final_remnant_properties_from_NSBH_waveform(
1466 self, source=False, parameters=[
1467 "baryonic_torus_mass", "final_mass", "final_spin"
1468 ]
1469 ):
1470 approximant = self._retrieve_approximant()
1471 if source:
1472 sample_params = [
1473 "mass_1_source", "mass_2_source", "spin_1z", "lambda_2"
1474 ]
1475 else:
1476 sample_params = ["mass_1", "mass_2", "spin_1z", "lambda_2"]
1477 samples = self.specific_parameter_samples(sample_params)
1478 _data = _check_NSBH_approximant(
1479 approximant, samples[0], samples[1], samples[2], samples[3],
1480 _raise=False
1481 )
1482 if _data is None:
1483 return
1484 _mapping = {
1485 "220_quasinormal_mode_frequency": 0, "tidal_disruption_frequency": 1,
1486 "baryonic_torus_mass": 2, "compactness_2": 3,
1487 "final_mass": 4, "final_spin": 5
1488 }
1489 for param in parameters:
1490 self.append_data(param, _data[_mapping[param]])
1491 if "final_mass" in parameters:
1492 self.extra_kwargs["meta_data"]["final_mass_NR_fits"] = "Zappa2019"
1493 if "final_spin" in parameters:
1494 self.extra_kwargs["meta_data"]["final_spin_NR_fits"] = "Zappa2019"
1495 if "baryonic_torus_mass" in parameters:
1496 self.extra_kwargs["meta_data"]["baryonic_torus_mass_fits"] = (
1497 "Foucart2012"
1498 )
1499 if "220_quasinormal_mode_frequency" in parameters:
1500 self.extra_kwargs["meta_data"]["quasinormal_mode_fits"] = (
1501 "London2019"
1502 )
1503 if "tidal_disruption_frequency" in parameters:
1504 probabilities = NSBH_merger_type(
1505 samples[0], samples[1], samples[2], samples[3],
1506 approximant=approximant,
1507 _ringdown=_data[_mapping["220_quasinormal_mode_frequency"]],
1508 _disruption=_data[_mapping["tidal_disruption_frequency"]],
1509 _torus=_data[_mapping["baryonic_torus_mass"]], percentages=True
1510 )
1511 self.extra_kwargs["meta_data"]["NSBH_merger_type_probabilities"] = (
1512 probabilities
1513 )
1514 self.extra_kwargs["meta_data"]["tidal_disruption_frequency_fits"] = (
1515 "Pannarale2018"
1516 )
1517 ratio = (
1518 _data[_mapping["tidal_disruption_frequency"]]
1519 / _data[_mapping["220_quasinormal_mode_frequency"]]
1520 )
1521 self.append_data(
1522 "tidal_disruption_frequency_ratio", ratio
1523 )
1525 def _final_remnant_properties_from_waveform(
1526 self, non_precessing=False, parameters=["final_mass", "final_spin"],
1527 ):
1528 f_low = self._retrieve_stored_frequency("f_start")
1529 approximant = self._retrieve_approximant()
1530 if "delta_t" in self.extra_kwargs["meta_data"].keys():
1531 delta_t = self.extra_kwargs["meta_data"]["delta_t"]
1532 else:
1533 delta_t = 1. / 4096
1534 if "seob" in approximant.lower():
1535 logger.warning(
1536 "Could not find 'delta_t' in the meta data. Using {} as "
1537 "default.".format(delta_t)
1538 )
1539 if non_precessing:
1540 sample_params = [
1541 "mass_1", "mass_2", "empty", "empty", "spin_1z", "empty",
1542 "empty", "spin_2z", "iota", "luminosity_distance",
1543 "phase"
1544 ]
1545 else:
1546 sample_params = [
1547 "mass_1", "mass_2", "spin_1x", "spin_1y", "spin_1z",
1548 "spin_2x", "spin_2y", "spin_2z", "iota", "luminosity_distance",
1549 "phase", "reference_frequency"
1550 ]
1551 samples = self.specific_parameter_samples(sample_params)
1552 ind = self.parameters.index("spin_1x")
1553 _data, fits = _final_from_initial_BBH(
1554 *samples[:8], iota=samples[8], luminosity_distance=samples[9],
1555 f_low=[f_low] * len(samples[0]), f_ref=samples[-1],
1556 phi_ref=samples[10], delta_t=1. / 4096, approximant=approximant,
1557 xphm_flags=self.extra_kwargs["meta_data"].get("approximant_flags", {}),
1558 return_fits_used=True, multi_process=self.multi_process
1559 )
1560 data = {"final_mass": _data[0], "final_spin": _data[1]}
1561 for param in parameters:
1562 self.append_data(param, data[param])
1563 self.extra_kwargs["meta_data"]["{}_NR_fits".format(param)] = fits
1565 def _final_mass_of_merger(self, evolved=False):
1566 param = self._evolved_vs_non_evolved_parameter(
1567 "final_mass", evolved=evolved, non_precessing=self.non_precessing
1568 )
1569 spin_1z_param = self._evolved_vs_non_evolved_parameter(
1570 "spin_1z", evolved=evolved, core_param=True,
1571 non_precessing=self.non_precessing
1572 )
1573 spin_2z_param = self._evolved_vs_non_evolved_parameter(
1574 "spin_2z", evolved=evolved, core_param=True,
1575 non_precessing=self.non_precessing
1576 )
1577 samples = self.specific_parameter_samples([
1578 "mass_1", "mass_2", spin_1z_param, spin_2z_param
1579 ])
1580 final_mass, fits = final_mass_of_merger(
1581 *samples, return_fits_used=True,
1582 )
1583 self.append_data(param, final_mass)
1584 self.extra_kwargs["meta_data"]["final_mass_NR_fits"] = fits
1586 def _final_mass_source(self, evolved=False):
1587 param = self._evolved_vs_non_evolved_parameter(
1588 "final_mass", evolved=evolved, non_precessing=self.non_precessing
1589 )
1590 samples = self.specific_parameter_samples([param, "redshift"])
1591 final_mass_source = _source_from_detector(
1592 samples[0], samples[1]
1593 )
1594 self.append_data(param.replace("mass", "mass_source"), final_mass_source)
1596 def _final_spin_of_merger(self, non_precessing=False, evolved=False):
1597 param = self._evolved_vs_non_evolved_parameter(
1598 "final_spin", evolved=evolved, non_precessing=self.non_precessing
1599 )
1600 samples = self._precessing_vs_non_precessing_parameters(
1601 non_precessing=non_precessing, evolved=evolved
1602 )
1603 final_spin, fits = final_spin_of_merger(
1604 *samples, return_fits_used=True,
1605 )
1606 self.append_data(param, final_spin)
1607 self.extra_kwargs["meta_data"]["final_spin_NR_fits"] = fits
1609 def _radiated_energy(self, evolved=False):
1610 param = self._evolved_vs_non_evolved_parameter(
1611 "radiated_energy", evolved=evolved, non_precessing=self.non_precessing
1612 )
1613 final_mass_param = self._evolved_vs_non_evolved_parameter(
1614 "final_mass_source", evolved=evolved, non_precessing=self.non_precessing
1615 )
1616 samples = self.specific_parameter_samples([
1617 "total_mass_source", final_mass_param
1618 ])
1619 radiated_energy = samples[0] - samples[1]
1620 self.append_data(param, radiated_energy)
1622 def _cos_angle(self, parameter_to_add, reverse=False):
1623 if reverse:
1624 samples = self.specific_parameter_samples(
1625 ["cos_" + parameter_to_add])
1626 cos_samples = np.arccos(samples[0])
1627 else:
1628 samples = self.specific_parameter_samples(
1629 [parameter_to_add.split("cos_")[1]]
1630 )
1631 cos_samples = np.cos(samples[0])
1632 self.append_data(parameter_to_add, cos_samples)
1634 def source_frame_from_detector_frame(self, detector_frame_parameter):
1635 samples = self.specific_parameter_samples(
1636 [detector_frame_parameter, "redshift"]
1637 )
1638 source_frame = _source_from_detector(samples[0], samples[1])
1639 self.append_data(
1640 "{}_source".format(detector_frame_parameter), source_frame
1641 )
1643 def _check_parameters(self):
1644 params = ["mass_1", "mass_2", "a_1", "a_2", "mass_1_source", "mass_2_source",
1645 "mass_ratio", "total_mass", "chirp_mass"]
1646 for i in params:
1647 if i in self.parameters:
1648 samples = self.specific_parameter_samples([i])
1649 if "mass" in i:
1650 cond = any(np.array(samples[0]) <= 0.)
1651 else:
1652 cond = any(np.array(samples[0]) < 0.)
1653 if cond:
1654 if "mass" in i:
1655 ind = np.argwhere(np.array(samples[0]) <= 0.)
1656 else:
1657 ind = np.argwhere(np.array(samples[0]) < 0.)
1658 logger.warning(
1659 "Removing %s samples because they have unphysical "
1660 "values (%s < 0)" % (len(ind), i)
1661 )
1662 for i in np.arange(len(ind) - 1, -1, -1):
1663 self.samples.remove(list(np.array(self.samples)[ind[i][0]]))
1665 def generate_all_posterior_samples(self):
1666 logger.debug("Starting to generate all derived posteriors")
1667 evolve_condition = (
1668 True if self.evolve_spins_forwards and self.compute_remnant else False
1669 )
1670 if "cos_theta_jn" in self.parameters and "theta_jn" not in self.parameters:
1671 self._cos_angle("theta_jn", reverse=True)
1672 if "cos_iota" in self.parameters and "iota" not in self.parameters:
1673 self._cos_angle("iota", reverse=True)
1674 if "cos_tilt_1" in self.parameters and "tilt_1" not in self.parameters:
1675 self._cos_angle("tilt_1", reverse=True)
1676 if "cos_tilt_2" in self.parameters and "tilt_2" not in self.parameters:
1677 self._cos_angle("tilt_2", reverse=True)
1678 angles = [
1679 "a_1", "a_2", "a_1_azimuthal", "a_1_polar", "a_2_azimuthal",
1680 "a_2_polar"
1681 ]
1682 if all(i in self.parameters for i in angles):
1683 self._component_spins_from_azimuthal_and_polar_angles()
1684 spin_magnitudes = ["a_1", "a_2"]
1685 angles = ["phi_jl", "tilt_1", "tilt_2", "phi_12"]
1686 cartesian = ["spin_1x", "spin_1y", "spin_1z", "spin_2x", "spin_2y", "spin_2z"]
1687 cond1 = all(i in self.parameters for i in spin_magnitudes)
1688 cond2 = all(i in self.parameters for i in angles)
1689 cond3 = all(i in self.parameters for i in cartesian)
1690 for _param in spin_magnitudes:
1691 if _param in self.parameters and not cond2 and not cond3:
1692 _index = _param.split("a_")[1]
1693 _spin = self.specific_parameter_samples(_param)
1694 _tilt = np.arccos(np.sign(_spin))
1695 self.append_data("tilt_{}".format(_index), _tilt)
1696 _spin_ind = self.parameters.index(_param)
1697 for num, i in enumerate(self.samples):
1698 self.samples[num][_spin_ind] = abs(self.samples[num][_spin_ind])
1700 if not cond2 and not cond3 and self.add_zero_spin:
1701 for _param in spin_magnitudes:
1702 if _param not in self.parameters:
1703 _spin = np.zeros(len(self.samples))
1704 self.append_data(_param, _spin)
1705 _index = _param.split("a_")[1]
1706 self.append_data("spin_{}z".format(_index), _spin)
1707 self._check_parameters()
1708 if "cos_theta_jn" in self.parameters and "theta_jn" not in self.parameters:
1709 self._cos_angle("theta_jn", reverse=True)
1710 if "cos_iota" in self.parameters and "iota" not in self.parameters:
1711 self._cos_angle("iota", reverse=True)
1712 if "cos_tilt_1" in self.parameters and "tilt_1" not in self.parameters:
1713 self._cos_angle("tilt_1", reverse=True)
1714 if "cos_tilt_2" in self.parameters and "tilt_2" not in self.parameters:
1715 self._cos_angle("tilt_2", reverse=True)
1716 if "luminosity_distance" not in self.parameters:
1717 if "redshift" in self.parameters:
1718 self._dL_from_z()
1719 if "redshift" not in self.parameters:
1720 if "luminosity_distance" in self.parameters:
1721 self._z_from_dL()
1722 if "comoving_volume" in self.parameters:
1723 self._z_from_comoving_volume()
1724 if "comoving_distance" not in self.parameters:
1725 if "redshift" in self.parameters:
1726 self._comoving_distance_from_z()
1727 if "comoving_volume" not in self.parameters:
1728 if "redshift" in self.parameters:
1729 self._comoving_volume_from_z()
1731 if "mass_ratio" not in self.parameters and "symmetric_mass_ratio" in \
1732 self.parameters:
1733 self._q_from_eta()
1734 if "mass_ratio" not in self.parameters and "mass_1" in self.parameters \
1735 and "mass_2" in self.parameters:
1736 self._q_from_m1_m2()
1737 if "mass_ratio" in self.parameters:
1738 ind = self.parameters.index("mass_ratio")
1739 median = np.median([i[ind] for i in self.samples])
1740 if median > 1.:
1741 self._invert_q()
1742 if "inverted_mass_ratio" not in self.parameters and "mass_ratio" in \
1743 self.parameters:
1744 self._invq_from_q()
1745 if "chirp_mass" not in self.parameters and "total_mass" in self.parameters:
1746 self._mchirp_from_mtotal_q()
1747 if "mass_1" not in self.parameters and "chirp_mass" in self.parameters:
1748 self._m1_from_mchirp_q()
1749 if "mass_2" not in self.parameters and "chirp_mass" in self.parameters:
1750 self._m2_from_mchirp_q()
1751 if "mass_1" not in self.parameters and "total_mass" in self.parameters:
1752 self._m1_from_mtotal_q()
1753 if "mass_2" not in self.parameters and "total_mass" in self.parameters:
1754 self._m2_from_mtotal_q()
1755 if "mass_1" in self.parameters and "mass_2" in self.parameters:
1756 if "total_mass" not in self.parameters:
1757 self._mtotal_from_m1_m2()
1758 if "chirp_mass" not in self.parameters:
1759 self._mchirp_from_m1_m2()
1760 if "symmetric_mass_ratio" not in self.parameters:
1761 self._eta_from_m1_m2()
1762 if "redshift" in self.parameters:
1763 if "mass_1_source" not in self.parameters:
1764 if "mass_1" in self.parameters:
1765 self._m1_source_from_m1_z()
1766 if "mass_1_source" in self.parameters:
1767 if "mass_1" not in self.parameters:
1768 self._m1_from_m1_source_z()
1769 if "mass_2_source" not in self.parameters:
1770 if "mass_2" in self.parameters:
1771 self._m2_source_from_m2_z()
1772 if "mass_2_source" in self.parameters:
1773 if "mass_2" not in self.parameters:
1774 self._m2_from_m2_source_z()
1775 if "total_mass_source" not in self.parameters:
1776 if "total_mass" in self.parameters:
1777 self._mtotal_source_from_mtotal_z()
1778 if "total_mass_source" in self.parameters:
1779 if "total_mass" not in self.parameters:
1780 self._mtotal_from_mtotal_source_z()
1781 if "chirp_mass_source" not in self.parameters:
1782 if "chirp_mass" in self.parameters:
1783 self._mchirp_source_from_mchirp_z()
1784 if "chirp_mass_source" in self.parameters:
1785 if "chirp_mass" not in self.parameters:
1786 self._mchirp_from_mchirp_source_z()
1788 if "reference_frequency" not in self.parameters:
1789 self._reference_frequency()
1790 condition1 = "phi_12" not in self.parameters
1791 condition2 = "phi_1" in self.parameters and "phi_2" in self.parameters
1792 if condition1 and condition2:
1793 self._phi_12_from_phi1_phi2()
1795 check_for_evolved_parameter = lambda suffix, param, params: (
1796 param not in params and param + suffix not in params if
1797 len(suffix) else param not in params
1798 )
1800 if "mass_1" in self.parameters and "mass_2" in self.parameters:
1801 spin_components = [
1802 "spin_1x", "spin_1y", "spin_1z", "spin_2x", "spin_2y", "spin_2z",
1803 "iota"
1804 ]
1805 angles = ["a_1", "a_2", "tilt_1", "tilt_2", "theta_jn"]
1806 if all(i in self.parameters for i in spin_components):
1807 self._spin_angles()
1808 if all(i in self.parameters for i in angles):
1809 samples = self.specific_parameter_samples(["tilt_1", "tilt_2"])
1810 cond1 = all(i in [0, np.pi] for i in samples[0])
1811 cond2 = all(i in [0, np.pi] for i in samples[1])
1812 if cond1 and cond1:
1813 self._non_precessing_component_spins()
1814 else:
1815 angles = [
1816 "phi_jl", "phi_12", "reference_frequency"]
1817 if all(i in self.parameters for i in angles):
1818 self._component_spins()
1819 cond1 = "spin_1x" in self.parameters and "spin_1y" in self.parameters
1820 if "phi_1" not in self.parameters and cond1:
1821 self._phi1_from_spins()
1822 cond1 = "spin_2x" in self.parameters and "spin_2y" in self.parameters
1823 if "phi_2" not in self.parameters and cond1:
1824 self._phi2_from_spins()
1825 evolve_spins_params = ["tilt_1", "tilt_2", "phi_12"]
1826 if self.evolve_spins_backwards:
1827 if all(i in self.parameters for i in evolve_spins_params):
1828 try:
1829 self._evolve_spins(forward=False)
1830 except EvolveSpinError:
1831 # Raised when approximant is unknown to lalsimulation or
1832 # lalsimulation.SimInspiralGetSpinFreqFromApproximant is
1833 # equal to lalsimulation.SIM_INSPIRAL_SPINS_CASEBYCASE
1834 pass
1835 for suffix in ["_infinity", "_infinity_only_prec_avg", ""]:
1836 if "spin_1z{}".format(suffix) not in self.parameters:
1837 _params = ["a_1", "tilt_1{}".format(suffix)]
1838 if all(i in self.parameters for i in _params):
1839 self._aligned_spin_from_magnitude_tilts(
1840 primary=True, suffix=suffix
1841 )
1842 if "spin_2z{}".format(suffix) not in self.parameters:
1843 _params = ["a_2", "tilt_2{}".format(suffix)]
1844 if all(i in self.parameters for i in _params):
1845 self._aligned_spin_from_magnitude_tilts(
1846 secondary=True, suffix=suffix
1847 )
1848 if "chi_eff{}".format(suffix) not in self.parameters:
1849 _params = ["spin_1z{}".format(suffix), "spin_2z{}".format(suffix)]
1850 if all(i in self.parameters for i in _params):
1851 self._chi_eff(suffix=suffix)
1852 if any(
1853 _p.format(suffix) not in self.parameters for _p in
1854 ["chi_p{}", "chi_p_2spin"]
1855 ):
1856 _params = [
1857 "a_1", "tilt_1{}".format(suffix), "a_2",
1858 "tilt_2{}".format(suffix)
1859 ]
1860 _cartesian_params = ["spin_1x", "spin_1y", "spin_2x", "spin_2y"]
1861 if "chi_p{}".format(suffix) not in self.parameters:
1862 if all(i in self.parameters for i in _params):
1863 self._chi_p_from_tilts(suffix=suffix)
1864 elif all(i in self.parameters for i in _cartesian_params):
1865 self._chi_p()
1866 if "chi_p_2spin" not in self.parameters:
1867 if all(i in self.parameters for i in _cartesian_params):
1868 self._chi_p_2spin()
1869 if "beta" not in self.parameters:
1870 beta_components = [
1871 "mass_1", "mass_2", "phi_jl", "tilt_1", "tilt_2", "phi_12",
1872 "a_1", "a_2", "reference_frequency", "phase"
1873 ]
1874 if all(i in self.parameters for i in beta_components):
1875 self._beta()
1876 polytrope_params = ["log_pressure", "gamma_1", "gamma_2", "gamma_3"]
1877 if all(param in self.parameters for param in polytrope_params):
1878 if "lambda_1" not in self.parameters or "lambda_2" not in self.parameters:
1879 self._lambda1_lambda2_from_polytrope_EOS()
1880 spectral_params = [
1881 "spectral_decomposition_gamma_{}".format(num) for num in
1882 np.arange(4)
1883 ]
1884 if all(param in self.parameters for param in spectral_params):
1885 if "lambda_1" not in self.parameters or "lambda_2" not in self.parameters:
1886 self._lambda1_lambda2_from_spectral_decomposition_EOS()
1887 if "lambda_tilde" in self.parameters and "lambda_1" not in self.parameters:
1888 self._lambda1_from_lambda_tilde()
1889 if "lambda_2" not in self.parameters and "lambda_1" in self.parameters:
1890 self._lambda2_from_lambda1()
1891 if "lambda_1" in self.parameters and "lambda_2" in self.parameters:
1892 if "lambda_tilde" not in self.parameters:
1893 self._lambda_tilde_from_lambda1_lambda2()
1894 if "delta_lambda" not in self.parameters:
1895 self._delta_lambda_from_lambda1_lambda2()
1896 if "psi" in self.parameters:
1897 dpsi_parameters = ["theta_jn", "phi_jl", "beta"]
1898 if all(i in self.parameters for i in dpsi_parameters):
1899 if "psi_J" not in self.parameters:
1900 self._psi_J()
1902 evolve_suffix = "_non_evolved"
1903 final_spin_params = ["a_1", "a_2"]
1904 non_precessing_NR_params = ["spin_1z", "spin_2z"]
1905 if evolve_condition:
1906 final_spin_params += [
1907 "tilt_1_evolved", "tilt_2_evolved", "phi_12_evolved"
1908 ]
1909 non_precessing_NR_params = [
1910 "{}_evolved".format(i) for i in non_precessing_NR_params
1911 ]
1912 evolve_suffix = "_evolved"
1913 if all(i in self.parameters for i in evolve_spins_params):
1914 try:
1915 self._evolve_spins(final_velocity=self.evolve_spins_forwards)
1916 except EvolveSpinError:
1917 # Raised when approximant is unknown to lalsimulation or
1918 # lalsimulation.SimInspiralGetSpinFreqFromApproximant is
1919 # equal to lalsimulation.SIM_INSPIRAL_SPINS_CASEBYCASE
1920 evolve_condition = False
1921 else:
1922 evolve_condition = False
1923 else:
1924 final_spin_params += ["tilt_1", "tilt_2", "phi_12"]
1926 condition_peak_luminosity = check_for_evolved_parameter(
1927 evolve_suffix, "peak_luminosity", self.parameters
1928 )
1929 condition_final_spin = check_for_evolved_parameter(
1930 evolve_suffix, "final_spin", self.parameters
1931 )
1932 condition_final_mass = check_for_evolved_parameter(
1933 evolve_suffix, "final_mass", self.parameters
1934 )
1935 if (self.NRSurrogate or self.waveform_fit) and self.compute_remnant:
1936 parameters = []
1937 _default = ["final_mass", "final_spin"]
1938 if self.NRSurrogate:
1939 _default.append("final_kick")
1940 function = self._final_remnant_properties_from_NRSurrogate
1941 else:
1942 final_spin_params = [
1943 "spin_1x", "spin_1y", "spin_1z", "spin_2x",
1944 "spin_2y", "spin_2z"
1945 ]
1946 function = self._final_remnant_properties_from_waveform
1948 for param in _default:
1949 if param not in self.parameters:
1950 parameters.append(param)
1951 # We already know that lambda_2 is in the posterior table if
1952 # self.NSBH = True
1953 if self.NSBH and "spin_1z" in self.parameters:
1954 self._final_remnant_properties_from_NSBH_waveform()
1955 elif all(i in self.parameters for i in final_spin_params):
1956 function(non_precessing=False, parameters=parameters)
1957 elif all(i in self.parameters for i in non_precessing_NR_params):
1958 function(non_precessing=True, parameters=parameters)
1959 if all(i in self.parameters for i in non_precessing_NR_params):
1960 if condition_peak_luminosity or self.force_non_evolved:
1961 if not self.NSBH:
1962 self._peak_luminosity_of_merger(evolved=evolve_condition)
1963 elif self.compute_remnant:
1964 if all(i in self.parameters for i in final_spin_params):
1965 if condition_final_spin or self.force_non_evolved:
1966 self._final_spin_of_merger(evolved=evolve_condition)
1967 elif all(i in self.parameters for i in non_precessing_NR_params):
1968 if condition_final_spin or self.force_non_evolved:
1969 self._final_spin_of_merger(
1970 non_precessing=True, evolved=False
1971 )
1972 if all(i in self.parameters for i in non_precessing_NR_params):
1973 if condition_peak_luminosity or self.force_non_evolved:
1974 self._peak_luminosity_of_merger(evolved=evolve_condition)
1975 if condition_final_mass or self.force_non_evolved:
1976 self._final_mass_of_merger(evolved=evolve_condition)
1978 # if NSBH system and self.compute_remnant = False and/or BBH fits
1979 # fits used, only calculate baryonic_torus_mass
1980 if self.NSBH and "spin_1z" in self.parameters:
1981 if "baryonic_torus_mass" not in self.parameters:
1982 self._final_remnant_properties_from_NSBH_waveform(
1983 parameters=["baryonic_torus_mass"]
1984 )
1985 # calculate compactness from Love-compactness relation
1986 if "lambda_1" in self.parameters and "compactness_1" not in self.parameters:
1987 self._NS_compactness_from_lambda(parameter="lambda_1")
1988 if "mass_1" in self.parameters and "baryonic_mass_1" not in self.parameters:
1989 self._NS_baryonic_mass(primary=True)
1990 if "lambda_2" in self.parameters and "compactness_2" not in self.parameters:
1991 self._NS_compactness_from_lambda(parameter="lambda_2")
1992 if "mass_2" in self.parameters and "baryonic_mass_2" not in self.parameters:
1993 self._NS_baryonic_mass(primary=False)
1994 for suffix in ["_infinity", "_infinity_only_prec_avg", ""]:
1995 for tilt in ["tilt_1", "tilt_2"]:
1996 cond1 = "cos_{}{}".format(tilt, suffix) not in self.parameters
1997 cond2 = "{}{}".format(tilt, suffix) in self.parameters
1998 if cond1 and cond2:
1999 self._cos_angle("cos_{}{}".format(tilt, suffix))
2000 evolve_suffix = "_non_evolved"
2001 if evolve_condition or self.NRSurrogate or self.waveform_fit or self.non_precessing:
2002 evolve_suffix = ""
2003 evolve_condition = True
2004 if "redshift" in self.parameters:
2005 condition_final_mass_source = check_for_evolved_parameter(
2006 evolve_suffix, "final_mass_source", self.parameters
2007 )
2008 if condition_final_mass_source or self.force_non_evolved:
2009 if "final_mass{}".format(evolve_suffix) in self.parameters:
2010 self._final_mass_source(evolved=evolve_condition)
2011 if "baryonic_torus_mass" in self.parameters:
2012 if "baryonic_torus_mass_source" not in self.parameters:
2013 self.source_frame_from_detector_frame(
2014 "baryonic_torus_mass"
2015 )
2016 if "baryonic_mass_1" in self.parameters:
2017 if "baryonic_mass_1_source" not in self.parameters:
2018 self.source_frame_from_detector_frame(
2019 "baryonic_mass_1"
2020 )
2021 if "baryonic_mass_2" in self.parameters:
2022 if "baryonic_mass_2_source" not in self.parameters:
2023 self.source_frame_from_detector_frame(
2024 "baryonic_mass_2"
2025 )
2026 if "total_mass_source" in self.parameters:
2027 if "final_mass_source{}".format(evolve_suffix) in self.parameters:
2028 condition_radiated_energy = check_for_evolved_parameter(
2029 evolve_suffix, "radiated_energy", self.parameters
2030 )
2031 if condition_radiated_energy or self.force_non_evolved:
2032 self._radiated_energy(evolved=evolve_condition)
2033 if self.NSBH and "spin_1z" in self.parameters:
2034 if all(_p in self.parameters for _p in ["mass_1_source", "mass_2_source"]):
2035 _NSBH_parameters = []
2036 if "tidal_disruption_frequency" not in self.parameters:
2037 _NSBH_parameters.append("tidal_disruption_frequency")
2038 if "220_quasinormal_mode_frequency" not in self.parameters:
2039 _NSBH_parameters.append("220_quasinormal_mode_frequency")
2040 if len(_NSBH_parameters):
2041 self._final_remnant_properties_from_NSBH_waveform(
2042 parameters=_NSBH_parameters, source=True
2043 )
2044 location = ["geocent_time", "ra", "dec"]
2045 if all(i in self.parameters for i in location):
2046 try:
2047 self._time_in_each_ifo()
2048 self._time_delay_between_ifos()
2049 except Exception as e:
2050 logger.warning(
2051 "Failed to generate posterior samples for the time in each "
2052 "detector because %s" % (e)
2053 )
2054 if any("_matched_filter_snr_angle" in i for i in self.parameters):
2055 if any("_matched_filter_abs_snr" in i for i in self.parameters):
2056 self._ifo_snr()
2057 if any("_optimal_snr" in i for i in self.parameters):
2058 if "network_optimal_snr" not in self.parameters:
2059 self._optimal_network_snr()
2060 if any("_matched_filter_snr" in i for i in self.parameters):
2061 if "network_matched_filter_snr" not in self.parameters:
2062 self._matched_filter_network_snr()
2063 if self.multipole_snr:
2064 rho_hm_parameters = [
2065 "mass_1", "mass_2", "spin_1z", "spin_2z", "psi", "iota", "ra",
2066 "dec", "geocent_time", "luminosity_distance", "phase"
2067 ]
2068 cond = [
2069 int(mm) for mm in ['21', '33', '44'] if
2070 "network_{}_multipole_snr".format(mm) not in self.parameters
2071 ]
2072 if all(i in self.parameters for i in rho_hm_parameters) and len(cond):
2073 try:
2074 logger.warning(
2075 "Starting to calculate the SNR in the {} multipole{}. "
2076 "This may take some time".format(
2077 " and ".join([str(mm) for mm in cond]),
2078 "s" if len(cond) > 1 else ""
2079 )
2080 )
2081 self._rho_hm(cond)
2082 except ImportError as e:
2083 logger.warning(e)
2084 elif len(cond):
2085 logger.warning(
2086 "Unable to calculate the multipole SNR because it requires "
2087 "samples for {}".format(
2088 ", ".join(
2089 [i for i in rho_hm_parameters if i not in self.parameters]
2090 )
2091 )
2092 )
2093 if "network_precessing_snr" not in self.parameters and self.precessing_snr:
2094 rho_p_parameters = [
2095 "mass_1", "mass_2", "beta", "psi_J", "a_1", "a_2", "tilt_1",
2096 "tilt_2", "phi_12", "theta_jn", "phi_jl", "ra", "dec", "geocent_time",
2097 "phi_jl"
2098 ]
2099 if all(i in self.parameters for i in rho_p_parameters):
2100 try:
2101 logger.warning(
2102 "Starting to calculate the precessing SNR. This may take "
2103 "some time"
2104 )
2105 self._rho_p()
2106 except ImportError as e:
2107 logger.warning(e)
2108 else:
2109 logger.warning(
2110 "Unable to calculate the precessing SNR because requires "
2111 "samples for {}".format(
2112 ", ".join(
2113 [i for i in rho_p_parameters if i not in self.parameters]
2114 )
2115 )
2116 )
2117 if "theta_jn" in self.parameters and "cos_theta_jn" not in self.parameters:
2118 self._cos_angle("cos_theta_jn")
2119 if "theta_jn" in self.parameters and "viewing_angle" not in self.parameters:
2120 self._viewing_angle()
2121 if "iota" in self.parameters and "cos_iota" not in self.parameters:
2122 self._cos_angle("cos_iota")
2123 remove_parameters = [
2124 "tilt_1_evolved", "tilt_2_evolved", "phi_12_evolved",
2125 "spin_1z_evolved", "spin_2z_evolved", "reference_frequency",
2126 "minimum_frequency"
2127 ]
2128 for param in remove_parameters:
2129 if param in self.parameters:
2130 ind = self.parameters.index(param)
2131 self.parameters.remove(self.parameters[ind])
2132 for i in self.samples:
2133 del i[ind]