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