Coverage for pesummary/gw/pycbc.py: 27.3%
44 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-12-09 22:34 +0000
« 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)
7from pesummary.utils.utils import logger
8import numpy as np
9try:
10 import pycbc.psd
11 from pycbc.filter import sigma, overlap_cplx
12 from pycbc import pnutils
13except ImportError:
14 logger.warning(error_msg.format("pycbc"))
16__author__ = ["Charlie Hoy <charlie.hoy@ligo.org>"]
19def optimal_snr(
20 template, psd, low_frequency_cutoff=20., high_frequency_cutoff=1024.
21):
22 """Calculate the loudness of the waveform. See Duncan Brown’s thesis for
23 details
25 Parameters
26 ----------
27 template: pycbc.type.frequencyseries.FrequencySeries
28 waveform you wish to calculate the loudness for
29 psd: pycbc.type.frequencyseries.FrequencySeries
30 psd to use when calculating the integral
31 low_frequency_cutoff: float, optional
32 low frequency cut-off to start calculating the integral
33 """
34 return sigma(
35 template, psd, low_frequency_cutoff=low_frequency_cutoff,
36 high_frequency_cutoff=high_frequency_cutoff
37 )
40def compute_the_overlap(template, data, psd, **kwargs):
41 """Wrapper for pycbc.filter.overlap_cplx
43 Parameters
44 ----------
45 template: pycbc.types.frequencyseries.FrequencySeries
46 frequency domain template
47 data: pycbc.types.frequencyseries.FrequencySeries
48 frequency domain data
49 psd: pycbc.types.frequencyseries.FrequencySeries
50 the PSD to use when computing the overlap
51 **kwargs: dict, optional
52 all additional kwargs passed to pycbc.filter.overlap_cplx
53 """
54 return overlap_cplx(template, data, psd=psd, **kwargs)
57def pycbc_default_psd(
58 mass_1=None, mass_2=None, chi_eff=None, spin_1z=None, spin_2z=None, psd=None,
59 duration=None, df=None, detectors=["H1", "L1"], f_low=20., f_final=1024.
60):
61 """Return a dictionary of pycbc psds for a given detector network
63 Parameters
64 ----------
65 mass_1: np.ndarray, optional
66 array of primary masses to use when working out the duration and
67 consequently delta_f for the PSD
68 mass_2: np.ndarray, optional
69 array of secondary masses to use when working out the duration and
70 consequently delta_f for the PSD
71 chi_eff: np.ndarray, optional
72 array of effective spin samples to use when working out the duration and
73 consequently delta_f for the psd. Used as an approximant when spin_1z
74 and spin_2z samples are not provided
75 spin_1z: np.ndarray, optional
76 array of samples for the primary spin parallel to the orbital angular
77 momentum to use when computing the duration and consequently delta_f for
78 the psd
79 spin_2z: np.ndarray, optional
80 array of samples for the secondary spin parallel to the orbital angular
81 momentum to use when computing the duration and consequently delta_f for
82 the psd
83 psd: func, optional
84 pycbc function to use when calculating the default psd. Default is
85 aLIGOZeroDetHighPower
86 duration: np.npdarray
87 array containing the durations for a given set of waveform. The maximum
88 duration is then used to calculate the delta_f for the PSD
89 detectors: list, optional
90 the detector network you wish to calculate psds for
91 f_low: float, optional
92 the low frequency to start computing the psd from. Default 20Hz
93 f_final: float, optional
94 the highest frequency to finish computing the psd. Default 1024Hz
95 """
96 if psd is None:
97 psd = pycbc.psd.aLIGOZeroDetHighPower
98 elif isinstance(psd, str):
99 psd = getattr(pycbc.psd, psd)
101 logger.warning("No PSD provided. Using '{}' for the psd".format(psd.__name__))
102 _required = [mass_1, mass_2]
103 if df is None:
104 cond1 = all(i is not None for i in _required + [spin_1z, spin_2z])
105 cond2 = all(i is not None for i in _required + [chi_eff])
106 if cond1 and duration is None:
107 duration = pnutils.get_imr_duration(
108 mass_1, mass_2, spin_1z, spin_2z, f_low
109 )
110 elif cond2 and duration is None:
111 logger.warning(
112 "Could not find samples for spin_1z and spin_2z. We will "
113 "assume that spin_1z = spin_2z = chi_eff to estimate the "
114 "maximum IMR duration. This is used to estimate delta_f for "
115 "the PSD."
116 )
117 duration = pnutils.get_imr_duration(
118 mass_1, mass_2, chi_eff, chi_eff, f_low
119 )
120 elif duration is None:
121 raise ValueError(
122 "Please provide either 'spin_1z' and 'spin_2z' or 'chi_eff' "
123 "samples to estimate the maximum IMR duration. This is used to "
124 "estimate delta_f for the PSD."
125 )
126 duration = np.max(duration)
127 t_len = 2**np.ceil(np.log2(duration) + 1)
128 df = 1. / t_len
129 flen = int(f_final / df) + 1
130 _psd = psd(flen, df, f_low)
131 setattr(_psd, "__name__", psd.__name__)
132 return {ifo: _psd for ifo in detectors}
135def interpolate_psd(psd, low_freq_cutoff, delta_f):
136 """Interpolate PSD to a new delta_f
138 Parameters
139 ----------
140 low_freq_cutoff: float
141 Frequencies below this value are set to zero.
142 delta_f : float, optional
143 Frequency resolution of the frequency series in Hertz.
144 """
145 from pesummary.gw.file.psd import PSD
146 from pycbc.psd import estimate
147 if isinstance(psd, PSD):
148 psd = psd.to_pycbc(low_freq_cutoff)
149 return estimate.interpolate(psd, delta_f)