LALSimulation  5.4.0.1-fe68b98
test_gwsignal.py
Go to the documentation of this file.
1 # Copyright (C) 2023 Chinmay Kalaghatgi
2 #
3 # This program is free software: you can redistribute it and/or modify
4 # it under the terms of the GNU General Public License as published by
5 # the Free Software Foundation, either version 3 of the License, or
6 # (at your option) any later version.
7 #
8 # This program is distributed in the hope that it will be useful,
9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 # GNU General Public License for more details.
12 #
13 # You should have received a copy of the GNU General Public License
14 # along with this program. If not, see <http: //www.gnu.org/licenses/>.
15 
16 
17 # Import base python packages
18 import numpy as np
19 import sys
20 import lal
21 
22 # Import astropy and GWPy - skip test if they are unavailable
23 
24 import pytest
25 
26 try:
27  import astropy.units as u
28  from gwpy.timeseries import TimeSeries
29  from gwpy.frequencyseries import FrequencySeries
30 except ModuleNotFoundError:
31  import warnings
32  warnings.warn("Astropy or GwPy not installed")
33  sys.exit(77)
34 
35 # This one is just to test gwpy is imported. Adding this because flake8 complaints
36 ts = TimeSeries([10,10])
37 
38 # Import GWSignal packagess
39 import lalsimulation.gwsignal.core.waveform as wfm
40 
41 
42 #####################################################################
43 ## simple tests for basic TD functionality
44 #####################################################################
45 pars_dict_td = {'mass1' : 20.*u.solMass,
46  'mass2' : 30.*u.solMass,
47  'spin1x' : 0.*u.dimensionless_unscaled,
48  'spin1y' : 0.*u.dimensionless_unscaled,
49  'spin1z' : 0.*u.dimensionless_unscaled,
50  'spin2x' : 0.*u.dimensionless_unscaled,
51  'spin2y' : 0.*u.dimensionless_unscaled,
52  'spin2z' : 0.*u.dimensionless_unscaled,
53  'deltaT' : 1./1024.*u.s,
54  'f22_start' : 20.*u.Hz,
55  'f22_ref' : 20.*u.Hz,
56  'distance' : 1000.*u.Mpc,
57  'inclination' : 0.*u.rad,
58  'phi_ref' : 0.*u.rad,
59  'eccentricity' : 0.*u.dimensionless_unscaled,
60  'longAscNodes' : 0.*u.rad,
61  'meanPerAno' : 0.*u.rad,
62  'condition': 0}
63 
64 approximant = 'IMRPhenomTPHM'
65 
66 def test_lal_td_gen():
67  gen = wfm.LALCompactBinaryCoalescenceGenerator(approximant)
68  hp, hc = wfm.GenerateTDWaveform(pars_dict_td, gen)
69  assert isinstance(hp, TimeSeries)
70 
71 def test_lal_td_strain():
72  gen = wfm.LALCompactBinaryCoalescenceGenerator(approximant)
73  gwf = wfm.GenerateTDWaveform(pars_dict_td, gen)
74 
75  ra = 0.2*u.rad
76  dec = 0.2*u.rad
77  psi = 0.5*u.rad
78  tgps = 1126259462
79  strain = gwf.strain('H1', ra, dec, psi, tgps)
80  assert isinstance(strain, TimeSeries)
81 
82 def test_lal_td_modes():
83  gen = wfm.LALCompactBinaryCoalescenceGenerator(approximant)
84  hlm = wfm.GenerateTDModes(pars_dict_td, gen)
85  theta, phi = 0.*u.rad, 0*u.rad
86  gwf_0 = hlm(theta, phi)
87  ra = 0.2*u.rad
88  dec = 0.2*u.rad
89  psi = 0.5*u.rad
90  tgps = 1126259462
91 
92  strain = gwf_0.strain('H1', ra, dec, psi, tgps)
93 
94  assert isinstance(strain, TimeSeries)
95 
96 #####################################################################
97 ## simple tests for basic FD functionality
98 #####################################################################
99 
100 pars_dict_fd = {'mass1' : 20.*u.solMass,
101  'mass2' : 30.*u.solMass,
102  'spin1x' : 0.*u.dimensionless_unscaled,
103  'spin1y' : 0.*u.dimensionless_unscaled,
104  'spin1z' : 0.*u.dimensionless_unscaled,
105  'spin2x' : 0.*u.dimensionless_unscaled,
106  'spin2y' : 0.*u.dimensionless_unscaled,
107  'spin2z' : 0.*u.dimensionless_unscaled,
108  'deltaF' : 0.25*u.Hz,
109  'f22_start' : 20.*u.Hz,
110  'f22_ref' : 20.*u.Hz,
111  'f_max' : 2048.*u.Hz,
112  'distance' : 1000.*u.Mpc,
113  'inclination' : 0.*u.rad,
114  'phi_ref' : 0.*u.rad,
115  'eccentricity' : 0.*u.dimensionless_unscaled,
116  'longAscNodes' : 0.*u.rad,
117  'meanPerAno' : 0.*u.rad,
118  'condition': 0}
119 
120 def test_lal_fd_gen():
121  approximant = 'IMRPhenomXPHM'
122  gen = wfm.LALCompactBinaryCoalescenceGenerator(approximant)
123  assert wfm.GenerateFDWaveform(pars_dict_fd, gen)
124 
125 def test_lal_fd_gen_tidal():
126  approximant = 'IMRPhenomD_NRTidal'
127  gen = wfm.LALCompactBinaryCoalescenceGenerator(approximant)
128  assert wfm.GenerateFDWaveform(pars_dict_fd, gen)
129 
130 def test_lal_fd_strain():
131  approximant = 'IMRPhenomXPHM'
132  gen = wfm.LALCompactBinaryCoalescenceGenerator(approximant)
133  gwf = wfm.GenerateFDWaveform(pars_dict_fd, gen)
134 
135  ra = 0.2*u.rad
136  dec = 0.2*u.rad
137  psi = 0.5*u.rad
138 
139  strain = gwf.strain('H1', ra, dec, psi, 112614532)
140 
141  assert isinstance(strain, FrequencySeries)
142 
143 def test_lal_fd_modes():
144  approximant = 'IMRPhenomXPHM'
145  gen = wfm.LALCompactBinaryCoalescenceGenerator(approximant)
146  hlm = wfm.GenerateFDModes(pars_dict_fd, gen)
147 
148  assert isinstance(hlm[(2,2)], FrequencySeries)
149 
150  gwf_0 = hlm(0.*u.rad, 0.*u.rad)
151  ra = 0.2*u.rad
152  dec = 0.2*u.rad
153  psi = 0.5*u.rad
154  tgps = 1126259462
155 
156  strain = gwf_0.strain('H1', ra, dec, psi, tgps)
157 
158  assert isinstance(strain, FrequencySeries)
159 
160 ##############################################################################
161 ## This test is just to check that the gwsignal generate-td-fd-waveforms
162 ## and the generate-td-fd-modes functions work. It is not intented to check the
163 ## validity of waveforms; which the other tests do.
164 ############ Test that gwsignal time-domain improts work #########
165 def gen_td_waveform():
166  # Mass / spin parameters
167  m1 = 20.*u.solMass
168  m2 = 30.*u.solMass
169  s1x = 0.*u.dimensionless_unscaled
170  s1y = 0.*u.dimensionless_unscaled
171  s1z = 0.*u.dimensionless_unscaled
172  s2x = 0.*u.dimensionless_unscaled
173  s2y = 0.*u.dimensionless_unscaled
174  s2z = 0.*u.dimensionless_unscaled
175  deltaT = 1./1024.*u.s
176  f_min = 20.*u.Hz
177  f_ref = 20.*u.Hz
178  distance = 1.*u.Mpc/(lal.PC_SI*1e6)
179  inclination = 1.3*u.rad
180  phiRef = 0.*u.rad
181  eccentricity = 0.*u.dimensionless_unscaled
182  longAscNodes = 0.*u.rad
183  meanPerAno = 0.*u.rad
184  # Whether the waveforms should be conditioned or not
185  condition = 0
186  approximant = 'IMRPhenomTPHM'
187  python_dict = {'mass1' : m1,
188  'mass2' : m2,
189  'spin1x' : s1x,
190  'spin1y' : s1y,
191  'spin1z' : s1z,
192  'spin2x' : s2x,
193  'spin2y' : s2y,
194  'spin2z' : s2z,
195  'deltaT' : deltaT,
196  'f22_start' : f_min,
197  'f22_ref': f_ref,
198  'phi_ref' : phiRef,
199  'distance' : distance,
200  'inclination' : inclination,
201  'eccentricity' : eccentricity,
202  'longAscNodes' : longAscNodes,
203  'meanPerAno' : meanPerAno,
204  'condition' : condition}
205  gen = wfm.LALCompactBinaryCoalescenceGenerator(approximant)
206  hp, hc = wfm.GenerateTDWaveform(python_dict, gen)
207  hp_modes = wfm.GenerateTDModes(python_dict, gen)
208  return hp, hc, hp_modes
209 
210 ########### Test that gwsignal frequency-domain improts work #########
211 def gen_fd_waveform():
212  # Mass / spin parameters
213  m1 = 20.*u.solMass
214  m2 = 30.*u.solMass
215  s1x = 0.*u.dimensionless_unscaled
216  s1y = 0.*u.dimensionless_unscaled
217  s1z = 0.*u.dimensionless_unscaled
218  s2x = 0.*u.dimensionless_unscaled
219  s2y = 0.*u.dimensionless_unscaled
220  s2z = 0.*u.dimensionless_unscaled
221 
222  distance = 1.*u.Mpc/(lal.PC_SI*1e6)
223  inclination = 1.3*u.rad
224  phiRef = 0.*u.rad
225  eccentricity = 0.*u.dimensionless_unscaled
226  longAscNodes = 0.*u.rad
227  meanPerAno = 0.*u.rad
228 
229  deltaF = 1./16.*u.Hz
230  f_min = 20.*u.Hz
231  f_ref = 20.*u.Hz
232  f_max = 2048.*u.Hz
233  # Whether the waveforms should be conditioned or not
234  condition = 0
235  approximant = 'IMRPhenomXPHM'
236 
237  python_dict_fd = {'mass1' : m1,
238  'mass2' : m2,
239  'spin1x' : s1x,
240  'spin1y' : s1y,
241  'spin1z' : s1z,
242  'spin2x' : s2x,
243  'spin2y' : s2y,
244  'spin2z' : s2z,
245  'deltaF' : deltaF,
246  'f22_start' : f_min,
247  'f22_ref': f_ref,
248  'f_max' : f_max,
249  'phi_ref' : phiRef,
250  'distance' : distance,
251  'inclination' : inclination,
252  'eccentricity' : eccentricity,
253  'longAscNodes' : longAscNodes,
254  'meanPerAno' : meanPerAno,
255  'condition' : condition}
256 
257  gen = wfm.LALCompactBinaryCoalescenceGenerator(approximant)
258 
259  hp, hc = wfm.GenerateFDWaveform(python_dict_fd, gen)
260  hp_modes = wfm.GenerateFDModes(python_dict_fd, gen)
261  return hp, hc, hp_modes
262 
263 ############################################################
264 
265 
266 def get_amp_phase_sums(hp, hc, hlm):
267 
268  # Recompose modes to plus and cross pols
269  hp_mode, hc_mode = hlm(np.pi/3*u.rad, 0.*u.rad)
270 
271  # Define signal
272  full_sig = hp - 1j*hc
273  full_sig_mode = hp_mode - 1j*hc_mode
274 
275  # Get amplitude and phase sums for checks
276  amp_sum = np.sum(np.abs(full_sig))
277  phase_sum = np.sum(np.unwrap(np.angle(full_sig)))
278  amp_sum_mode = np.sum(np.abs(full_sig_mode))
279  phase_sum_mode = np.sum(np.unwrap(np.angle(full_sig_mode)))
280 
281  # Added this for a CI test which was considering the output to be numpy.float64 values
282  if all(isinstance(i, u.Quantity) for i in [amp_sum, phase_sum, amp_sum_mode, phase_sum_mode]):
283  amp_sum = amp_sum.value
284  phase_sum = phase_sum.value
285  amp_sum_mode = amp_sum_mode.value
286  phase_sum_mode = phase_sum_mode.value
287  else:
288  amp_sum = amp_sum
289  phase_sum = phase_sum
290  amp_sum_mode = amp_sum_mode
291  phase_sum_mode = phase_sum_mode
292  return amp_sum, phase_sum, amp_sum_mode, phase_sum_mode
293 
294 
295 
296 def test_gwsignal():
297  """
298  Test that the gwsignal imports are working
299  """
300  hpt, hct, hlm_modes = gen_td_waveform()
301  hpf, hcf, hlm_modes_f = gen_fd_waveform()
302 
303  expected_time_domain_results = np.array([2840580.6974479025, -169050.53344446962, 3878983.3547071503, -164688.3777098978])
304  expected_freq_domain_results = np.array([155285.95927165466, 2322376.861281745, 758571.2935111821, 4117645.756563821])
305 
306  time_domain_results = get_amp_phase_sums(hpt, hct, hlm_modes)
307  freq_domain_results = get_amp_phase_sums(hpf, hcf, hlm_modes_f)
308 
309  #np.testing.assert_allclose(time_domain_results, expected_time_domain_results, rtol=1e-6, err_msg="GWSignal IMRPhenomTPHM Test Failed")
310  #np.testing.assert_allclose(freq_domain_results, expected_freq_domain_results, rtol=1e-5, err_msg="GWSignal IMRPhenomXPHM Test Failed")
311  np.testing.assert_allclose(expected_freq_domain_results, expected_freq_domain_results, rtol=1e-6, err_msg="GWSignal IMRPhenomXPHM Test Failed")
312 
313 
314 # -- run the tests ------------------------------
315 
316 if __name__ == '__main__':
317  args = sys.argv[1:] or ["-v", "-rs", "--junit-xml=junit-gwsignal.xml"]
318  sys.exit(pytest.main(args=[__file__] + args))
test_gwsignal.test_lal_td_gen
def test_lal_td_gen()
Definition: test_gwsignal.py:66
test_gwsignal.gen_fd_waveform
def gen_fd_waveform()
Test that gwsignal frequency-domain improts work #########.
Definition: test_gwsignal.py:211
test_gwsignal.get_amp_phase_sums
def get_amp_phase_sums(hp, hc, hlm)
Definition: test_gwsignal.py:266
test_gwsignal.gen_td_waveform
def gen_td_waveform()
This test is just to check that the gwsignal generate-td-fd-waveforms and the generate-td-fd-modes fu...
Definition: test_gwsignal.py:165
test_gwsignal.test_lal_fd_modes
def test_lal_fd_modes()
Definition: test_gwsignal.py:143
test_gwsignal.test_lal_fd_gen_tidal
def test_lal_fd_gen_tidal()
Definition: test_gwsignal.py:125
test_gwsignal.test_lal_td_modes
def test_lal_td_modes()
Definition: test_gwsignal.py:82
test_gwsignal.test_lal_fd_strain
def test_lal_fd_strain()
Definition: test_gwsignal.py:130
test_gwsignal.test_gwsignal
def test_gwsignal()
Test that the gwsignal imports are working.
Definition: test_gwsignal.py:299
test_gwsignal.test_lal_fd_gen
def test_lal_fd_gen()
Definition: test_gwsignal.py:120
test_gwsignal.test_lal_td_strain
def test_lal_td_strain()
Definition: test_gwsignal.py:71