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

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

2 

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")) 

15 

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

17 

18 

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 

24 

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 ) 

38 

39 

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

41 """Wrapper for pycbc.filter.overlap_cplx 

42 

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) 

55 

56 

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 

62 

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) 

100 

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} 

133 

134 

135def interpolate_psd(psd, low_freq_cutoff, delta_f): 

136 """Interpolate PSD to a new delta_f 

137 

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)