Coverage for pesummary/gw/pycbc.py: 41.9%

1# Licensed under an MIT style license -- see LICENSE.md 

2 

3from pesummary.utils.utils import logger, import_error_msg 

4import numpy as np 

5try: 

6 import pycbc.psd 

7 from pycbc.filter import sigma, overlap_cplx 

8 from pycbc import pnutils 

9except ImportError: 

10 logger.warning(import_error_msg.format("pycbc")) 

11 

12__author__ = ["Charlie Hoy <charlie.hoy@ligo.org>"] 

13 

14 

15def optimal_snr( 

16 template, psd, low_frequency_cutoff=20., high_frequency_cutoff=1024. 

17): 

18 """Calculate the loudness of the waveform. See Duncan Brown’s thesis for 

19 details 

20 

21 Parameters 

22 ---------- 

23 template: pycbc.type.frequencyseries.FrequencySeries 

24 waveform you wish to calculate the loudness for 

25 psd: pycbc.type.frequencyseries.FrequencySeries 

26 psd to use when calculating the integral 

27 low_frequency_cutoff: float, optional 

28 low frequency cut-off to start calculating the integral 

29 """ 

30 return sigma( 

31 template, psd, low_frequency_cutoff=low_frequency_cutoff, 

32 high_frequency_cutoff=high_frequency_cutoff 

33 ) 

34 

35 

36def compute_the_overlap(template, data, psd, **kwargs): 

37 """Wrapper for pycbc.filter.overlap_cplx 

38 

39 Parameters 

40 ---------- 

41 template: pycbc.types.frequencyseries.FrequencySeries 

42 frequency domain template 

43 data: pycbc.types.frequencyseries.FrequencySeries 

44 frequency domain data 

45 psd: pycbc.types.frequencyseries.FrequencySeries 

46 the PSD to use when computing the overlap 

47 **kwargs: dict, optional 

48 all additional kwargs passed to pycbc.filter.overlap_cplx 

49 """ 

50 return overlap_cplx(template, data, psd=psd, **kwargs) 

51 

52 

53def pycbc_default_psd( 

54 mass_1=None, mass_2=None, chi_eff=None, spin_1z=None, spin_2z=None, psd=None, 

55 duration=None, df=None, detectors=["H1", "L1"], f_low=20., f_final=1024. 

56): 

57 """Return a dictionary of pycbc psds for a given detector network 

58 

59 Parameters 

60 ---------- 

61 mass_1: np.ndarray, optional 

62 array of primary masses to use when working out the duration and 

63 consequently delta_f for the PSD 

64 mass_2: np.ndarray, optional 

65 array of secondary masses to use when working out the duration and 

66 consequently delta_f for the PSD 

67 chi_eff: np.ndarray, optional 

68 array of effective spin samples to use when working out the duration and 

69 consequently delta_f for the psd. Used as an approximant when spin_1z 

70 and spin_2z samples are not provided 

71 spin_1z: np.ndarray, optional 

72 array of samples for the primary spin parallel to the orbital angular 

73 momentum to use when computing the duration and consequently delta_f for 

74 the psd 

75 spin_2z: np.ndarray, optional 

76 array of samples for the secondary spin parallel to the orbital angular 

77 momentum to use when computing the duration and consequently delta_f for 

78 the psd 

79 psd: func, optional 

80 pycbc function to use when calculating the default psd. Default is 

81 aLIGOZeroDetHighPower 

82 duration: np.npdarray 

83 array containing the durations for a given set of waveform. The maximum 

84 duration is then used to calculate the delta_f for the PSD 

85 detectors: list, optional 

86 the detector network you wish to calculate psds for 

87 f_low: float, optional 

88 the low frequency to start computing the psd from. Default 20Hz 

89 f_final: float, optional 

90 the highest frequency to finish computing the psd. Default 1024Hz 

91 """ 

92 if psd is None: 

93 psd = pycbc.psd.aLIGOZeroDetHighPower 

94 elif isinstance(psd, str): 

95 psd = getattr(pycbc.psd, psd) 

96 

97 logger.warning("No PSD provided. Using '{}' for the psd".format(psd.__name__)) 

98 _required = [mass_1, mass_2] 

99 if df is None: 

100 cond1 = all(i is not None for i in _required + [spin_1z, spin_2z]) 

101 cond2 = all(i is not None for i in _required + [chi_eff]) 

102 if cond1 and duration is None: 

103 duration = pnutils.get_imr_duration( 

104 mass_1, mass_2, spin_1z, spin_2z, f_low 

105 ) 

106 elif cond2 and duration is None: 

107 logger.warning( 

108 "Could not find samples for spin_1z and spin_2z. We will " 

109 "assume that spin_1z = spin_2z = chi_eff to estimate the " 

110 "maximum IMR duration. This is used to estimate delta_f for " 

111 "the PSD." 

112 ) 

113 duration = pnutils.get_imr_duration( 

114 mass_1, mass_2, chi_eff, chi_eff, f_low 

115 ) 

116 elif duration is None: 

117 raise ValueError( 

118 "Please provide either 'spin_1z' and 'spin_2z' or 'chi_eff' " 

119 "samples to estimate the maximum IMR duration. This is used to " 

120 "estimate delta_f for the PSD." 

121 ) 

122 duration = np.max(duration) 

123 t_len = 2**np.ceil(np.log2(duration) + 1) 

124 df = 1. / t_len 

125 flen = int(f_final / df) + 1 

126 _psd = psd(flen, df, f_low) 

127 setattr(_psd, "__name__", psd.__name__) 

128 return {ifo: _psd for ifo in detectors} 

129 

130 

131def interpolate_psd(psd, low_freq_cutoff, delta_f): 

132 """Interpolate PSD to a new delta_f 

133 

134 Parameters 

135 ---------- 

136 low_freq_cutoff: float 

137 Frequencies below this value are set to zero. 

138 delta_f : float, optional 

139 Frequency resolution of the frequency series in Hertz. 

140 """ 

141 from pesummary.gw.file.psd import PSD 

142 from pycbc.psd import estimate 

143 if isinstance(psd, PSD): 

144 psd = psd.to_pycbc(low_freq_cutoff) 

145 return estimate.interpolate(psd, delta_f)