Coverage for pesummary/tests/waveform_test.py: 100.0%

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

2 

3import numpy as np 

4from pesummary.gw.waveform import fd_waveform, td_waveform 

5from pycbc.waveform import get_fd_waveform, get_td_waveform 

6from lal import MSUN_SI 

7from lalsimulation import SimInspiralTransformPrecessingNewInitialConditions 

8 

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

10np.random.seed(1234) 

11 

12 

13class TestWaveformModule(): 

14 """Class to test pesummary.gw.waveform module 

15 """ 

16 def setup_method(self): 

17 """Setup the testing class 

18 """ 

19 self.n_samples = 20 

20 self.approx = "IMRPhenomPv2" 

21 self.mass_1 = np.random.uniform(20, 100, self.n_samples) 

22 self.mass_2 = np.random.uniform(5, self.mass_1, self.n_samples) 

23 self.a_1 = np.random.uniform(0, 1, self.n_samples) 

24 self.a_2 = np.random.uniform(0, 1, self.n_samples) 

25 self.tilt_1 = np.arccos(np.random.uniform(-1, 1, self.n_samples)) 

26 self.tilt_2 = np.arccos(np.random.uniform(-1, 1, self.n_samples)) 

27 self.phi_12 = np.random.uniform(0, 1, self.n_samples) 

28 self.theta_jn = np.arccos(np.random.uniform(-1, 1, self.n_samples)) 

29 self.phi_jl = np.random.uniform(0, 1, self.n_samples) 

30 self.f_low = [10.] * self.n_samples 

31 self.f_final = [1024.] * self.n_samples 

32 self.phase = np.random.uniform(0, 1, self.n_samples) 

33 self.distance = np.random.uniform(100, 500, self.n_samples) 

34 self.ra = np.random.uniform(0, np.pi, self.n_samples) 

35 self.dec = np.arccos(np.random.uniform(-1, 1, self.n_samples)) 

36 self.psi_l = np.random.uniform(0, 1, self.n_samples) 

37 self.time = [10000.] * self.n_samples 

38 self.spin_1z = self.a_1 * np.cos(self.tilt_1) 

39 self.spin_2z = self.a_2 * np.cos(self.tilt_2) 

40 

41 def _make_waveform(self, index, df, func=fd_waveform, approx=None, **kwargs): 

42 """Make a waveform for testing 

43 

44 Parameters 

45 ---------- 

46 index: int 

47 index of mass_1/mass_2/... array to use for waveform generation 

48 df: float 

49 difference in frequency samples to use when generating waveform. 

50 This should match the difference in frequency samples used for the 

51 PSD 

52 """ 

53 if approx is None: 

54 approx = self.approx 

55 base = dict( 

56 theta_jn=self.theta_jn[index], 

57 phi_jl=self.phi_jl[index], phase=self.phase[index], 

58 mass_1=self.mass_1[index], mass_2=self.mass_2[index], 

59 tilt_1=self.tilt_1[index], tilt_2=self.tilt_2[index], 

60 phi_12=self.phi_12[index], a_1=self.a_1[index], 

61 a_2=self.a_2[index], luminosity_distance=self.distance[index] 

62 ) 

63 try: 

64 wvfs = func( 

65 base, approx, df, self.f_low[index], self.f_final[index] / 2, 

66 f_ref=self.f_low[index], pycbc=True, **kwargs 

67 ) 

68 except Exception: 

69 wvfs = func( 

70 base, approx, df, self.f_low[index], f_ref=self.f_low[index], 

71 pycbc=True, **kwargs 

72 ) 

73 hp, hc = wvfs["h_plus"], wvfs["h_cross"] 

74 return hp, hc 

75 

76 def test_fd_waveform(self): 

77 """Test the pesummary.gw.waveform.fd_waveform function 

78 """ 

79 df = 1./128 

80 for i in range(self.n_samples): 

81 iota, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z = \ 

82 SimInspiralTransformPrecessingNewInitialConditions( 

83 self.theta_jn[i], self.phi_jl[i], self.tilt_1[i], 

84 self.tilt_2[i], self.phi_12[i], self.a_1[i], 

85 self.a_2[i], self.mass_1[i] * MSUN_SI, 

86 self.mass_2[i] * MSUN_SI, self.f_low[i], 

87 self.phase[i] 

88 ) 

89 pycbc_hp, pycbc_hc = get_fd_waveform( 

90 approximant=self.approx, mass1=self.mass_1[i], 

91 mass2=self.mass_2[i], spin1x=spin1x, 

92 spin1y=spin1y, spin1z=spin1z, spin2x=spin2x, 

93 spin2y=spin2y, spin2z=spin2z, 

94 inclination=iota, distance=self.distance[i], 

95 coa_phase=self.phase[i], f_lower=self.f_low[i], 

96 f_final=self.f_final[i] / 2, delta_f=df, 

97 f_ref=self.f_low[i] 

98 ) 

99 pesummary_hp, pesummary_hc = self._make_waveform(i, df) 

100 np.testing.assert_almost_equal( 

101 10**30 * np.array(pycbc_hp), 10**30 * np.array(pesummary_hp) 

102 ) 

103 np.testing.assert_almost_equal( 

104 np.array(pycbc_hp.sample_frequencies), 

105 np.array(pesummary_hp.sample_frequencies) 

106 ) 

107 

108 def test_td_waveform(self): 

109 """Test the pesummary.gw.waveform.td_waveform function 

110 """ 

111 dt = 1./1024 

112 for i in range(5): 

113 iota, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z = \ 

114 SimInspiralTransformPrecessingNewInitialConditions( 

115 self.theta_jn[i], self.phi_jl[i], self.tilt_1[i], 

116 self.tilt_2[i], self.phi_12[i], self.a_1[i], 

117 self.a_2[i], self.mass_1[i] * MSUN_SI, 

118 self.mass_2[i] * MSUN_SI, self.f_low[i], 

119 self.phase[i] 

120 ) 

121 pycbc_hp, pycbc_hc = get_td_waveform( 

122 approximant="SEOBNRv3", mass1=self.mass_1[i], 

123 mass2=self.mass_2[i], spin1x=spin1x, 

124 spin1y=spin1y, spin1z=spin1z, spin2x=spin2x, 

125 spin2y=spin2y, spin2z=spin2z, 

126 inclination=iota, distance=self.distance[i], 

127 coa_phase=self.phase[i], f_lower=self.f_low[i], 

128 delta_t=dt, f_ref=self.f_low[i] 

129 ) 

130 pesummary_hp, pesummary_hc = self._make_waveform( 

131 i, dt, func=td_waveform, approx="SEOBNRv3" 

132 ) 

133 np.testing.assert_almost_equal( 

134 10**30 * np.array(pycbc_hp), 10**30 * np.array(pesummary_hp) 

135 ) 

136 np.testing.assert_almost_equal( 

137 np.array(pycbc_hp.sample_times), 

138 np.array(pesummary_hp.sample_times) 

139 ) 

140 

141 def test_mode_array(self): 

142 """Test that the mode array option works for both td and fd waveforms 

143 """ 

144 df = 1./256 

145 iota, spin1x, spin1y, spin1z, spin2x, spin2y, spin2z = \ 

146 SimInspiralTransformPrecessingNewInitialConditions( 

147 self.theta_jn[0], self.phi_jl[0], self.tilt_1[0], 

148 self.tilt_2[0], self.phi_12[0], self.a_1[0], 

149 self.a_2[0], self.mass_1[0] * MSUN_SI, 

150 self.mass_2[0] * MSUN_SI, self.f_low[0], 

151 self.phase[0] 

152 ) 

153 pycbc_hp, pycbc_hc = get_fd_waveform( 

154 approximant="IMRPhenomPv3HM", mass1=self.mass_1[0], 

155 mass2=self.mass_2[0], spin1x=spin1x, 

156 spin1y=spin1y, spin1z=spin1z, spin2x=spin2x, 

157 spin2y=spin2y, spin2z=spin2z, 

158 inclination=iota, distance=self.distance[0], 

159 coa_phase=self.phase[0], f_lower=self.f_low[0], 

160 f_final=self.f_final[0] / 2, delta_f=df, 

161 f_ref=self.f_low[0], mode_array=[[2,2], [3,3]] 

162 ) 

163 pesummary_hp, pesummary_hc = self._make_waveform( 

164 0, df, approx="IMRPhenomPv3HM", mode_array=[[2,2], [3,3]] 

165 ) 

166 np.testing.assert_almost_equal( 

167 10**30 * np.array(pycbc_hp), 10**30 * np.array(pesummary_hp) 

168 ) 

169 np.testing.assert_almost_equal( 

170 np.array(pycbc_hp.sample_frequencies), 

171 np.array(pesummary_hp.sample_frequencies) 

172 ) 

173 dt = 1./2048 

174 pycbc_hp, pycbc_hc = get_td_waveform( 

175 approximant="SEOBNRv4PHM", mass1=self.mass_1[0], 

176 mass2=self.mass_2[0], spin1x=spin1x, 

177 spin1y=spin1y, spin1z=spin1z, spin2x=spin2x, 

178 spin2y=spin2y, spin2z=spin2z, 

179 inclination=iota, distance=self.distance[0], 

180 coa_phase=self.phase[0], f_lower=self.f_low[0], 

181 delta_t=dt, f_ref=self.f_low[0], mode_array=[[2,2], [3,3]] 

182 ) 

183 pesummary_hp, pesummary_hc = self._make_waveform( 

184 0, dt, approx="SEOBNRv4PHM", mode_array=[[2,2], [3,3]], 

185 func=td_waveform 

186 ) 

187 np.testing.assert_almost_equal( 

188 10**30 * np.array(pycbc_hp), 10**30 * np.array(pesummary_hp) 

189 ) 

190 np.testing.assert_almost_equal( 

191 np.array(pycbc_hp.sample_times), 

192 np.array(pesummary_hp.sample_times) 

193 )