Coverage for pesummary/tests/conversion_test.py: 16.4%
572 statements
« 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
3import os
4import socket
5import shutil
7import numpy as np
8import h5py
10from pesummary.gw.conversions import *
11from pesummary.gw.conversions.nrutils import *
12from pycbc import conversions
13import pytest
14import tempfile
16tmpdir = tempfile.TemporaryDirectory(prefix=".", dir=".").name
18__author__ = ["Charlie Hoy <charlie.hoy@ligo.org>"]
21def conversion_check(
22 pesummary_function, pesummary_args, other_function, other_args,
23 dp=8
24):
25 """Check that the conversions made by PESummary are the same as those
26 from pycbc
27 """
28 _pesummary = pesummary_function(*pesummary_args)
29 _other = other_function(*other_args)
30 assert np.testing.assert_almost_equal(_pesummary, _other, dp) is None
33class TestConversions(object):
34 def setup_method(self):
35 class Arguments(object):
36 mass1 = 10.
37 mass2 = 5.
38 mtotal = 30.
39 mchirp = 10.
40 q = 1. / 4.
41 eta = 0.214
42 iota = 0.5
43 spin1x = 0.75
44 spin1y = 0.
45 spin1z = 0.5
46 spin2x = 0.
47 spin2y = 0.
48 spin2z = 0.5
49 lambda1 = 500.
50 lambda2 = 500.
51 lambda_tilde = 1000.
53 theta_jn = [0.5, 0.5]
54 phi_jl = [0.3, 0.3]
55 tilt_1 = [0.5, 0.]
56 tilt_2 = [0., 0.]
57 phi_12 = [0., 0.]
58 a_1 = [0.5, 0.5]
59 a_2 = [0., 0.]
60 f_ref = [20., 20.]
61 phase = [0., 0.]
63 redshift = 0.5
64 l_distance = 500.
66 self.opts = Arguments()
68 def test_z_from_dL(self):
69 from bilby.gw.conversion import luminosity_distance_to_redshift
71 l_distance = np.random.randint(100, 5000, 20)
72 bilby_function = luminosity_distance_to_redshift
73 pesummary_function = z_from_dL_exact
74 conversion_check(
75 pesummary_function, [l_distance], bilby_function, [l_distance]
76 )
77 pesummary_function = z_from_dL_approx
78 conversion_check(
79 pesummary_function, [l_distance], bilby_function, [l_distance],
80 dp=4
81 )
83 def test_change_of_cosmology_for_z_from_dL(self):
84 from lalinference.bayespputils import calculate_redshift
86 l_distance = np.random.randint(100, 5000, 20)
87 lal_redshift = calculate_redshift(
88 np.atleast_2d(l_distance).T
89 ).T[0]
90 redshift = z_from_dL_exact(
91 l_distance, cosmology="Planck15_lal"
92 )
93 np.testing.assert_almost_equal(lal_redshift, redshift, 8)
95 def test_z_from_comoving_volume(self):
96 from pycbc.cosmology import redshift_from_comoving_volume
98 volume = np.random.randint(100**3, 5000**3, 20)
99 pycbc_function = redshift_from_comoving_volume
100 pesummary_function = z_from_comoving_volume_exact
101 conversion_check(
102 pesummary_function, [volume], pycbc_function, [volume, False]
103 )
104 pesummary_function = z_from_comoving_volume_approx
105 conversion_check(
106 pesummary_function, [volume], pycbc_function, [volume, False],
107 dp=4
108 )
110 def test_dL_from_z(self):
111 from bilby.gw.conversion import redshift_to_luminosity_distance
113 redshift = np.random.uniform(1, 5, 100)
114 bilby_function = redshift_to_luminosity_distance
115 pesummary_function = dL_from_z
116 conversion_check(
117 pesummary_function, [redshift], bilby_function, [redshift]
118 )
120 def test_comoving_volume_from_z(self):
121 from pycbc.cosmology import cosmological_quantity_from_redshift
123 redshift = np.random.uniform(1, 5, 100)
124 pycbc_function = cosmological_quantity_from_redshift
125 pesummary_function = comoving_volume_from_z
126 conversion_check(
127 pesummary_function, [redshift], pycbc_function,
128 [redshift, "comoving_volume"]
129 )
130 def test_comoving_distance_from_z(self):
131 from bilby.gw.conversion import redshift_to_comoving_distance
133 redshift = np.random.randint(1, 5, 100)
134 bilby_function = redshift_to_comoving_distance
135 pesummary_function = comoving_distance_from_z
136 conversion_check(
137 pesummary_function, [redshift], bilby_function, [redshift]
138 )
140 def test_m1_source_from_m1_z(self):
141 from bilby.gw.conversion import generate_source_frame_parameters
143 mass_1 = np.random.randint(5, 100, 100)
144 mass_2 = np.random.randint(2, mass_1, 100)
145 luminosity_distance = np.random.randint(100, 500, 100)
146 sample = generate_source_frame_parameters(
147 {"mass_1": mass_1, "mass_2": mass_2,
148 "luminosity_distance": luminosity_distance}
149 )
150 source_frame = generate_source_frame_parameters(sample)
151 assert np.testing.assert_almost_equal(
152 m1_source_from_m1_z(mass_1, sample["redshift"]), sample["mass_1_source"],
153 8
154 ) is None
156 def test_m1_m2_from_m1_source_m2_source_z(self):
157 from bilby.gw.conversion import generate_source_frame_parameters
159 mass_1_source = np.random.randint(5, 100, 100)
160 mass_2_source = np.random.randint(2, mass_1_source, 100)
161 redshift = np.random.randint(1, 10, 100)
162 luminosity_distance = dL_from_z(redshift)
164 # calculate mass_1 and mass_2 using pesummary
165 mass_1 = m1_from_m1_source_z(mass_1_source, redshift)
166 mass_2 = m2_from_m2_source_z(mass_2_source, redshift)
167 # use calculated mass_1/mass_2 to calculate mass_1_source/mass_2_source using
168 # bilby
169 sample = generate_source_frame_parameters(
170 {"mass_1": mass_1, "mass_2": mass_2,
171 "luminosity_distance": luminosity_distance}
172 )
173 source_frame = generate_source_frame_parameters(sample)
174 # confirm that bilby's mass_1_source/mass_2_source is the same as
175 # mass_1_source/mass_2_source that pesummary used
176 np.testing.assert_almost_equal(sample["mass_1_source"], mass_1_source, 6)
177 np.testing.assert_almost_equal(sample["mass_2_source"], mass_2_source, 6)
179 def test_m2_source_from_m2_z(self):
180 from bilby.gw.conversion import generate_source_frame_parameters
182 mass_1 = np.random.randint(5, 100, 100)
183 mass_2 = np.random.randint(2, mass_1, 100)
184 luminosity_distance = np.random.randint(100, 500, 100)
185 sample = generate_source_frame_parameters(
186 {"mass_1": mass_1, "mass_2": mass_2,
187 "luminosity_distance": luminosity_distance}
188 )
189 source_frame = generate_source_frame_parameters(sample)
190 assert np.testing.assert_almost_equal(
191 m2_source_from_m2_z(mass_1, sample["redshift"]), sample["mass_1_source"],
192 8
193 ) is None
195 def test_m_total_source_from_mtotal_z(self):
196 from bilby.gw.conversion import generate_source_frame_parameters
198 total_mass = np.random.randint(5, 100, 100)
199 luminosity_distance = np.random.randint(100, 500, 100)
200 sample = generate_source_frame_parameters(
201 {"total_mass": total_mass, "luminosity_distance": luminosity_distance}
202 )
203 source_frame = generate_source_frame_parameters(sample)
204 assert np.testing.assert_almost_equal(
205 m_total_source_from_mtotal_z(total_mass, sample["redshift"]),
206 sample["total_mass_source"], 8
207 ) is None
209 def test_m_total_from_mtotal_source_z(self):
210 total_mass_source, redshift = 20., self.opts.redshift
211 m_total = mtotal_from_mtotal_source_z(total_mass_source, redshift)
212 assert np.round(m_total, 4) == self.opts.mtotal
214 def test_mchirp_source_from_mchirp_z(self):
215 from bilby.gw.conversion import generate_source_frame_parameters
217 chirp_mass = np.random.randint(5, 100, 100)
218 luminosity_distance = np.random.randint(100, 500, 100)
219 sample = generate_source_frame_parameters(
220 {"chirp_mass": chirp_mass,
221 "luminosity_distance": luminosity_distance}
222 )
223 source_frame = generate_source_frame_parameters(sample)
224 assert np.testing.assert_almost_equal(
225 mchirp_source_from_mchirp_z(chirp_mass, sample["redshift"]),
226 sample["chirp_mass_source"],
227 8
228 ) is None
230 def test_mchirp_from_mchirp_source_z(self):
231 mchirp_source, redshift = 20./3., self.opts.redshift
232 mchirp = mchirp_from_mchirp_source_z(mchirp_source, redshift)
233 assert np.round(mchirp, 4) == self.opts.mchirp
235 def test_mchirp_from_m1_m2(self):
236 mass1 = np.random.randint(5, 100, 100)
237 mass2 = np.random.randint(2, mass1, 100)
238 pycbc_function = conversions.mchirp_from_mass1_mass2
239 pesummary_function = mchirp_from_m1_m2
240 conversion_check(
241 pesummary_function, [mass1, mass2], pycbc_function, [mass1, mass2]
242 )
244 def test_m_total_from_m1_m2(self):
245 mass1 = np.random.randint(5, 100, 100)
246 mass2 = np.random.randint(2, mass1, 100)
247 pycbc_function = conversions.mtotal_from_mass1_mass2
248 pesummary_function = m_total_from_m1_m2
249 conversion_check(
250 pesummary_function, [mass1, mass2], pycbc_function, [mass1, mass2]
251 )
253 def test_m1_from_mchirp_q(self):
254 mchirp = np.random.randint(5, 100, 100)
255 q = np.random.random(100)
256 mchirp, q = self.opts.mchirp, self.opts.q
257 pycbc_function = conversions.mass1_from_mchirp_q
258 pesummary_function = m1_from_mchirp_q
259 conversion_check(
260 pesummary_function, [mchirp, q], pycbc_function, [mchirp, 1./q]
261 )
263 def test_m2_from_mchirp_q(self):
264 mchirp = np.random.randint(5, 100, 100)
265 q = np.random.random(100)
266 pycbc_function = conversions.mass2_from_mchirp_q
267 pesummary_function = m2_from_mchirp_q
268 conversion_check(
269 pesummary_function, [mchirp, q], pycbc_function, [mchirp, 1./q]
270 )
272 def test_m1_from_mtotal_q(self):
273 mtotal = np.random.randint(5, 100, 100)
274 q = np.random.random(100)
275 pycbc_function = conversions.mass1_from_mtotal_q
276 pesummary_function = m1_from_mtotal_q
277 conversion_check(
278 pesummary_function, [mtotal, q], pycbc_function, [mtotal, 1./q]
279 )
281 def test_m2_from_mtotal_q(self):
282 mtotal = np.random.randint(5, 100, 100)
283 q = np.random.random(100)
284 pycbc_function = conversions.mass2_from_mtotal_q
285 pesummary_function = m2_from_mtotal_q
286 conversion_check(
287 pesummary_function, [mtotal, q], pycbc_function, [mtotal, 1./q]
288 )
290 def test_eta_from_m1_m2(self):
291 mass1 = np.random.randint(5, 100, 100)
292 mass2 = np.random.randint(2, mass1, 100)
293 pycbc_function = conversions.eta_from_mass1_mass2
294 pesummary_function = eta_from_m1_m2
295 conversion_check(
296 pesummary_function, [mass1, mass2], pycbc_function, [mass1, mass2]
297 )
299 def test_q_from_m1_m2(self):
300 mass1 = np.random.randint(5, 100, 100)
301 mass2 = np.random.randint(2, mass1, 100)
302 pycbc_function = conversions.invq_from_mass1_mass2
303 pesummary_function = q_from_m1_m2
304 conversion_check(
305 pesummary_function, [mass1, mass2], pycbc_function, [mass1, mass2]
306 )
308 def test_q_from_eta(self):
309 from bilby.gw.conversion import symmetric_mass_ratio_to_mass_ratio
311 eta = np.random.uniform(0, 0.25, 100)
312 bilby_function = symmetric_mass_ratio_to_mass_ratio
313 pesummary_function = q_from_eta
314 conversion_check(
315 pesummary_function, [eta], bilby_function, [eta]
316 )
318 def test_mchirp_from_mtotal_q(self):
319 mtotal, q = self.opts.mtotal, self.opts.q
320 mchirp = mchirp_from_mtotal_q(mtotal, q)
321 assert np.round(mchirp, 4) == 9.9906
323 def test_chi_p(self):
324 mass1 = np.random.randint(5, 100, 100)
325 mass2 = np.random.randint(2, mass1, 100)
326 spin1_mag = np.random.random(100)
327 spin1_ang = np.random.random(100)
328 spin2_mag = np.random.random(100)
329 spin2_ang = np.random.random(100)
330 spin1x = spin1_mag * np.cos(spin1_ang)
331 spin1y = spin1_mag * np.sin(spin1_ang)
332 spin2x = spin2_mag * np.cos(spin2_ang)
333 spin2y = spin2_mag * np.sin(spin2_ang)
334 pycbc_function = conversions.chi_p
335 pesummary_function = chi_p
336 conversion_check(
337 pesummary_function, [mass1, mass2, spin1x, spin1y, spin2x, spin2y],
338 pycbc_function, [mass1, mass2, spin1x, spin1y, spin2x, spin2y]
339 )
341 from lalsimulation import SimPhenomUtilsChiP
343 mass1, mass2 = self.opts.mass1, self.opts.mass2
344 spin1x, spin1y = self.opts.spin1x, self.opts.spin1y
345 spin1z, spin2x = self.opts.spin1z, self.opts.spin2x
346 spin2y, spin2z = self.opts.spin2y, self.opts.spin2z
347 chi_p_value = chi_p(mass1, mass2, spin1x, spin1y, spin2y, spin2y)
348 assert chi_p_value == 0.75
349 for i in range(100):
350 mass_1 = np.random.randint(10, 100)
351 mass_2 = np.random.randint(5, mass_1)
352 spin1 = np.random.random(3)
353 norm = np.sqrt(np.sum(np.square(spin1)))
354 spin1 /= norm
355 spin2 = np.random.random(3)
356 norm = np.sqrt(np.sum(np.square(spin2)))
357 spin2 /= norm
358 chi_p_value = chi_p(
359 mass_1, mass_2, spin1[0], spin1[1], spin2[0], spin2[1]
360 )
361 lal_value = SimPhenomUtilsChiP(
362 mass_1, mass_2, spin1[0], spin1[1], spin2[0], spin2[1]
363 )
364 assert np.testing.assert_almost_equal(chi_p_value, lal_value, 9) is None
366 def test_chi_eff(self):
367 mass1 = np.random.randint(5, 100, 100)
368 mass2 = np.random.randint(2, mass1, 100)
369 spin1z = np.random.uniform(-1, 1, 100)
370 spin2z = np.random.uniform(-1, 1, 100)
371 pycbc_function = conversions.chi_eff
372 pesummary_function = chi_eff
373 conversion_check(
374 pesummary_function, [mass1, mass2, spin1z, spin2z], pycbc_function,
375 [mass1, mass2, spin1z, spin2z]
376 )
378 def test_phi_12_from_phi1_phi2(self):
379 data = phi_12_from_phi1_phi2(0.2, 0.5)
380 assert data == 0.3
381 data = phi_12_from_phi1_phi2(0.5, 0.2)
382 rounded_data = np.round(data, 2)
383 assert rounded_data == 5.98
384 data = phi_12_from_phi1_phi2(np.array([0.5, 0.2]), np.array([0.3, 0.7]))
385 rounded_data = np.round(data, 2)
386 assert all(i == j for i,j in zip(rounded_data, [6.08, 0.5]))
388 def test_phi_from_spins(self):
389 def cart2sph(x, y):
390 return np.fmod(2 * np.pi + np.arctan2(y,x), 2 * np.pi)
391 assert phi1_from_spins(0.5, 0.2) == cart2sph(0.5, 0.2)
392 assert phi2_from_spins(0.1, 0.5) == cart2sph(0.1, 0.5)
394 def test_spin_angles(self):
395 from lalsimulation import SimInspiralTransformPrecessingWvf2PE
397 mass1 = np.random.uniform(5., 100., 100)
398 mass2 = np.random.uniform(2., mass1, 100)
399 inc = np.random.uniform(0, np.pi, 100)
400 spin1_mag = np.random.random(100)
401 spin1_ang = np.random.random(100)
402 spin2_mag = np.random.random(100)
403 spin2_ang = np.random.random(100)
404 spin1x = spin1_mag * np.cos(spin1_ang)
405 spin1y = spin1_mag * np.sin(spin1_ang)
406 spin2x = spin2_mag * np.cos(spin2_ang)
407 spin2y = spin2_mag * np.sin(spin2_ang)
408 spin1z = np.random.random(100) - (spin1x**2 + spin1y**2)**0.5
409 spin2z = np.random.random(100) - (spin1x**2 + spin1y**2)**0.5
410 f_ref = [20.0] * len(mass1)
411 phase = [0.4] * len(mass1)
412 lalsimulation_function = SimInspiralTransformPrecessingWvf2PE
413 pesummary_function = spin_angles
414 for ii in np.arange(len(mass1)):
415 conversion_check(
416 pesummary_function,
417 [mass1[ii], mass2[ii], inc[ii], spin1x[ii], spin1y[ii],
418 spin1z[ii], spin2x[ii], spin2y[ii], spin2z[ii], f_ref[ii],
419 phase[ii]],
420 lalsimulation_function,
421 [inc[ii], spin1x[ii], spin1y[ii], spin1z[ii], spin2x[ii],
422 spin2y[ii], spin2z[ii], mass1[ii], mass2[ii], f_ref[ii],
423 phase[ii]]
424 )
426 def test_component_spins(self):
427 from bilby.gw.conversion import bilby_to_lalsimulation_spins
428 from lal import MSUN_SI
430 mass1 = np.random.uniform(5., 100., 100)
431 mass2 = np.random.uniform(2., mass1, 100)
432 theta_jn = np.random.uniform(0, np.pi, 100)
433 phi_jl = np.random.uniform(0, np.pi, 100)
434 phi_12 = np.random.uniform(0, np.pi, 100)
435 a_1 = np.random.uniform(0, 1, 100)
436 a_2 = np.random.uniform(0, 1, 100)
437 tilt_1 = np.random.uniform(0, np.pi, 100)
438 tilt_2 = np.random.uniform(0, np.pi, 100)
439 f_ref = [20.] * len(mass1)
440 phase = [0.5] * len(mass2)
442 bilby_function = bilby_to_lalsimulation_spins
443 pesummary_function = component_spins
444 for ii in np.arange(len(mass1)):
445 conversion_check(
446 pesummary_function,
447 [theta_jn[ii], phi_jl[ii], tilt_1[ii], tilt_1[ii], phi_12[ii],
448 a_1[ii], a_2[ii], mass1[ii], mass2[ii], f_ref[ii], phase[ii]],
449 bilby_function,
450 [theta_jn[ii], phi_jl[ii], tilt_1[ii], tilt_1[ii], phi_12[ii],
451 a_1[ii], a_2[ii], mass1[ii]*MSUN_SI, mass2[ii]*MSUN_SI,
452 f_ref[ii], phase[ii]]
453 )
455 def test_time_in_each_ifo(self):
456 from pycbc.detector import Detector
457 from lal import TimeDelayFromEarthCenter
458 from lalsimulation import DetectorPrefixToLALDetector
460 optimal_ra, optimal_dec = -0.2559168059473027, 0.81079526383
461 time = time_in_each_ifo("H1", optimal_ra, optimal_dec, 0)
462 light_time = 6371*10**3 / (3.0*10**8)
463 assert -np.round(light_time, 4) == np.round(time, 4)
465 ra = np.random.uniform(-np.pi, np.pi, 100)
466 dec = np.random.uniform(0, 2*np.pi, 100)
467 time = np.random.uniform(10000, 20000, 100)
468 H1_time = time_in_each_ifo("H1", ra, dec, time)
469 pycbc_time = time + Detector("H1").time_delay_from_earth_center(
470 ra, dec, time
471 )
472 lal_time = time + np.array([TimeDelayFromEarthCenter(
473 DetectorPrefixToLALDetector('H1').location, ra[ii], dec[ii],
474 time[ii]) for ii in range(len(ra))
475 ])
476 difference = np.abs(lal_time - pycbc_time)
477 dp = np.floor(np.abs(np.log10(np.max(difference))))
478 assert np.testing.assert_almost_equal(H1_time, pycbc_time, dp) is None
479 assert np.testing.assert_almost_equal(H1_time, lal_time, dp) is None
482 def test_lambda_tilde_from_lambda1_lambda2(self):
483 lambda1 = np.random.uniform(0, 5000, 100)
484 lambda2 = np.random.uniform(0, 5000, 100)
485 mass1 = np.random.randint(5, 100, 100)
486 mass2 = np.random.randint(2, mass1, 100)
487 pycbc_function = conversions.lambda_tilde
488 pesummary_function = lambda_tilde_from_lambda1_lambda2
489 conversion_check(
490 pesummary_function, [lambda1, lambda2, mass1, mass2],
491 pycbc_function, [mass1, mass2, lambda1, lambda2]
492 )
494 def test_delta_lambda_from_lambda1_lambda2(self):
495 from bilby.gw.conversion import lambda_1_lambda_2_to_delta_lambda_tilde
497 lambda1 = np.random.uniform(0, 5000, 100)
498 lambda2 = np.random.uniform(0, 5000, 100)
499 mass1 = np.random.randint(5, 100, 100)
500 mass2 = np.random.randint(2, mass1, 100)
501 conversion_check(
502 delta_lambda_from_lambda1_lambda2,
503 [lambda1, lambda2, mass1, mass2],
504 lambda_1_lambda_2_to_delta_lambda_tilde,
505 [lambda1, lambda2, mass1, mass2]
506 )
508 def test_lambda1_from_lambda_tilde(self):
509 from bilby.gw.conversion import lambda_tilde_to_lambda_1_lambda_2
511 mass1 = np.random.randint(5, 100, 100)
512 mass2 = np.random.randint(2, mass1, 100)
513 lambda_tilde = np.random.uniform(-100, 100, 100)
514 lambda_1 = lambda_tilde_to_lambda_1_lambda_2(
515 lambda_tilde, mass1, mass2
516 )[0]
517 lambda1 = lambda1_from_lambda_tilde(lambda_tilde, mass1, mass2)
518 assert np.testing.assert_almost_equal(lambda1, lambda_1) is None
519 #assert np.round(lambda1, 4) == 192.8101
521 def test_lambda2_from_lambda1(self):
522 from bilby.gw.conversion import convert_to_lal_binary_neutron_star_parameters
524 mass1 = np.random.uniform(5, 50, 100)
525 mass2 = np.random.uniform(2, mass1, 100)
526 lambda_1 = np.random.uniform(0, 5000, 100)
527 sample = {"mass_1": mass1, "mass_2": mass2, "lambda_1": lambda_1}
528 convert = convert_to_lal_binary_neutron_star_parameters(sample)[0]
529 lambda2 = lambda2_from_lambda1(lambda_1, mass1, mass2)
530 diff = np.abs(lambda2 - convert["lambda_2"])
531 ind = np.argmax(diff)
532 assert np.testing.assert_almost_equal(lambda2, convert["lambda_2"], 5) is None
534 def test_network_snr(self):
535 snr_H1 = snr_L1 = snr_V1 = np.array([2., 3.])
536 assert network_snr([snr_H1[0], snr_L1[0], snr_V1[0]]) == np.sqrt(3) * 2
537 print(snr_H1)
538 network = network_snr([snr_H1, snr_L1, snr_V1])
539 print(network)
540 assert network[0] == np.sqrt(3) * 2
541 assert network[1] == np.sqrt(3) * 3
543 def test_network_matched_filter_snr(self):
544 """Samples taken from a lalinference result file
545 """
546 snr_mf_H1 = 7.950207935574794
547 snr_mf_L1 = 19.232672412819483
548 snr_mf_V1 = 3.666438738845737
549 snr_opt_H1 = 9.668043620320788
550 snr_opt_L1 = 19.0826463504282
551 snr_opt_V1 = 3.5578582036515236
552 network = network_matched_filter_snr(
553 [snr_mf_H1, snr_mf_L1, snr_mf_V1],
554 [snr_opt_H1, snr_opt_L1, snr_opt_V1]
555 )
556 np.testing.assert_almost_equal(network, 21.06984787727566)
557 network = network_matched_filter_snr(
558 [[snr_mf_H1] * 2, [snr_mf_L1] * 2, [snr_mf_V1] * 2],
559 [[snr_opt_H1] * 2, [snr_opt_L1] * 2, [snr_opt_V1] * 2]
560 )
561 np.testing.assert_almost_equal(
562 network, [21.06984787727566] * 2
563 )
564 snr_mf_H1 = 7.950207935574794 - 1.004962343498161 * 1j
565 snr_mf_L1 = 19.232672412819483 - 0.4646531569951501 * 1j
566 snr_mf_V1 = 3.666438738845737 - 0.08177741915398137 * 1j
567 network = network_matched_filter_snr(
568 [snr_mf_H1, snr_mf_L1, snr_mf_V1],
569 [snr_opt_H1, snr_opt_L1, snr_opt_V1]
570 )
571 np.testing.assert_almost_equal(network, 21.06984787727566)
573 def test_full_conversion(self):
574 from pesummary.utils.array import Array
575 from pesummary.gw.conversions import convert
576 from bilby.gw.conversion import convert_to_lal_binary_black_hole_parameters
577 from pandas import DataFrame
579 dictionary = {
580 "mass_1": [10.],
581 "mass_2": [2.],
582 "a_1": [0.2],
583 "a_2": [0.2],
584 "cos_theta_jn": [0.5],
585 "cos_tilt_1": [0.2],
586 "cos_tilt_2": [0.2],
587 "luminosity_distance": [0.2],
588 "geocent_time": [0.2],
589 "ra": [0.2],
590 "dec": [0.2],
591 "psi": [0.2],
592 "phase": [0.5],
593 "phi_12": [0.7],
594 "phi_jl": [0.25],
595 "H1_matched_filter_abs_snr": [10.0],
596 "H1_matched_filter_snr_angle": [0.],
597 "H1_matched_filter_snr": [10.0],
598 "H1_optimal_snr": [10.0]
599 }
600 data = convert(dictionary)
601 true_params = [
602 'mass_1', 'mass_2', 'a_1', 'a_2', 'cos_theta_jn', 'cos_tilt_1',
603 'cos_tilt_2', 'luminosity_distance', 'geocent_time', 'ra', 'dec',
604 'psi', 'phase', 'phi_12', 'phi_jl', 'H1_matched_filter_abs_snr',
605 'H1_matched_filter_snr_angle', 'H1_matched_filter_snr', 'theta_jn',
606 'tilt_1', 'tilt_2', 'mass_ratio', 'inverted_mass_ratio',
607 'total_mass', 'chirp_mass', 'symmetric_mass_ratio', 'iota',
608 'spin_1x', 'spin_1y', 'spin_1z', 'spin_2x', 'spin_2y', 'spin_2z',
609 'phi_1', 'phi_2', 'chi_eff', 'chi_p', 'final_spin_non_evolved',
610 'peak_luminosity_non_evolved', 'final_mass_non_evolved', 'redshift',
611 'comoving_distance', 'mass_1_source', 'mass_2_source',
612 'total_mass_source', 'chirp_mass_source',
613 'final_mass_source_non_evolved', 'radiated_energy_non_evolved',
614 'network_matched_filter_snr', 'cos_iota'
615 ]
616 assert all(i in data.parameters for i in true_params)
617 assert len(data.parameters) == len(set(data.parameters))
618 true = {
619 'mass_1': Array(dictionary["mass_1"]),
620 'mass_2': Array(dictionary["mass_2"]),
621 'a_1': Array(dictionary["a_1"]),
622 'a_2': Array(dictionary["a_2"]),
623 'cos_theta_jn': Array(dictionary["cos_theta_jn"]),
624 'cos_tilt_1': Array(dictionary["cos_tilt_1"]),
625 'cos_tilt_2': Array(dictionary["cos_tilt_2"]),
626 'luminosity_distance': Array(dictionary["luminosity_distance"]),
627 'geocent_time': Array(dictionary["geocent_time"]),
628 'ra': Array(dictionary["ra"]), 'dec': Array(dictionary["dec"]),
629 'psi': Array(dictionary["psi"]), 'phase': Array(dictionary["phase"]),
630 'phi_12': Array(dictionary["phi_12"]),
631 'phi_jl': Array(dictionary["phi_jl"]),
632 'H1_matched_filter_abs_snr': Array(dictionary["H1_matched_filter_abs_snr"]),
633 'H1_matched_filter_snr_angle': Array(dictionary["H1_matched_filter_snr_angle"]),
634 'H1_matched_filter_snr': Array(dictionary["H1_matched_filter_snr"]),
635 'H1_optimal_snr': Array(dictionary["H1_optimal_snr"]),
636 'theta_jn': Array(np.arccos(dictionary["cos_theta_jn"])),
637 'tilt_1': Array(np.arccos(dictionary["cos_tilt_1"])),
638 'tilt_2': Array(np.arccos(dictionary["cos_tilt_2"])),
639 'mass_ratio': Array([dictionary["mass_2"][0] / dictionary["mass_1"][0]]),
640 'inverted_mass_ratio': Array([dictionary["mass_1"][0] / dictionary["mass_2"][0]]),
641 'total_mass': Array([dictionary["mass_1"][0] + dictionary["mass_2"][0]]),
642 'chirp_mass': Array([3.67097772]),
643 'symmetric_mass_ratio': Array([0.13888889]),
644 'iota': Array([1.01719087]), 'spin_1x': Array([-0.18338713]),
645 'spin_1y': Array([0.06905911]), 'spin_1z': Array([0.04]),
646 'spin_2x': Array([-0.18475131]), 'spin_2y': Array([-0.06532192]),
647 'spin_2z': Array([0.04]), 'phi_1': Array([2.78144141]),
648 'phi_2': Array([3.48144141]), 'chi_eff': Array([0.04]),
649 'chi_p': Array([0.19595918]),
650 'final_spin_non_evolved': Array([0.46198316]),
651 'peak_luminosity_non_evolved': Array([1.01239394]),
652 'final_mass_non_evolved': Array([11.78221674]),
653 'redshift': Array([4.5191183e-05]),
654 'comoving_distance': Array([0.19999755]),
655 'mass_1_source': Array([9.99954811]),
656 'mass_2_source': Array([1.99990962]),
657 'total_mass_source': Array([11.99945773]),
658 'chirp_mass_source': Array([3.67081183]),
659 'final_mass_source_non_evolved': Array([11.78168431]),
660 'radiated_energy_non_evolved': Array([0.21777342]),
661 'network_matched_filter_snr': Array([10.0]),
662 'cos_iota': Array([0.52575756])
663 }
664 for key in true.keys():
665 np.testing.assert_almost_equal(
666 true[key][0], data[key][0], 5
667 )
668 convert = convert_to_lal_binary_black_hole_parameters(dictionary)[0]
669 for key, item in convert.items():
670 assert np.testing.assert_almost_equal(
671 item, true[key], 5
672 ) is None
674 def test_remove_parameter(self):
675 from pesummary.gw.conversions import convert
677 dictionary = {
678 "mass_1": np.random.uniform(5, 100, 100),
679 "mass_ratio": [0.1] * 100
680 }
681 dictionary["mass_2"] = np.random.uniform(2, dictionary["mass_1"], 100)
682 incorrect_mass_ratio = convert(dictionary)
683 data = convert(dictionary, regenerate=["mass_ratio"])
684 assert all(i != j for i, j in zip(
685 incorrect_mass_ratio["mass_ratio"], data["mass_ratio"]
686 ))
687 assert np.testing.assert_almost_equal(
688 data["mass_ratio"],
689 q_from_m1_m2(dictionary["mass_1"], dictionary["mass_2"]), 8
690 ) is None
693class TestPrecessingSNR(object):
694 """Test the precessing_snr conversion
695 """
696 def setup_method(self):
697 """Setup the testing class
698 """
699 np.random.seed(1234)
700 self.n_samples = 20
701 self.approx = "IMRPhenomPv2"
702 self.mass_1 = np.random.uniform(20, 100, self.n_samples)
703 self.mass_2 = np.random.uniform(5, self.mass_1, self.n_samples)
704 self.a_1 = np.random.uniform(0, 1, self.n_samples)
705 self.a_2 = np.random.uniform(0, 1, self.n_samples)
706 self.tilt_1 = np.arccos(np.random.uniform(-1, 1, self.n_samples))
707 self.tilt_2 = np.arccos(np.random.uniform(-1, 1, self.n_samples))
708 self.phi_12 = np.random.uniform(0, 1, self.n_samples)
709 self.theta_jn = np.arccos(np.random.uniform(-1, 1, self.n_samples))
710 self.phi_jl = np.random.uniform(0, 1, self.n_samples)
711 self.f_low = [20.] * self.n_samples
712 self.f_final = [1024.] * self.n_samples
713 self.phase = np.random.uniform(0, 1, self.n_samples)
714 self.distance = np.random.uniform(100, 500, self.n_samples)
715 self.ra = np.random.uniform(0, np.pi, self.n_samples)
716 self.dec = np.arccos(np.random.uniform(-1, 1, self.n_samples))
717 self.psi_l = np.random.uniform(0, 1, self.n_samples)
718 self.time = [10000.] * self.n_samples
719 self.beta = opening_angle(
720 self.mass_1, self.mass_2, self.phi_jl, self.tilt_1,
721 self.tilt_2, self.phi_12, self.a_1, self.a_2,
722 [20.] * self.n_samples, self.phase
723 )
724 self.spin_1z = self.a_1 * np.cos(self.tilt_1)
725 self.spin_2z = self.a_2 * np.cos(self.tilt_2)
726 self.psi_J = psi_J(self.psi_l, self.theta_jn, self.phi_jl, self.beta)
728 def test_harmonic_overlap(self):
729 """Test that the sum of 5 precessing harmonics matches exactly to the
730 original precessing waveform.
731 """
732 from pycbc import pnutils
733 from pycbc.psd import aLIGOZeroDetHighPower
734 from pycbc.detector import Detector
735 from pesummary.gw.conversions.snr import (
736 _calculate_precessing_harmonics, _dphi,
737 _make_waveform_from_precessing_harmonics
738 )
739 from pesummary.gw.waveform import fd_waveform
741 for i in range(self.n_samples):
742 duration = pnutils.get_imr_duration(
743 self.mass_1[i], self.mass_2[i], self.spin_1z[i],
744 self.spin_2z[i], self.f_low[i], "IMRPhenomD"
745 )
746 t_len = 2**np.ceil(np.log2(duration) + 1)
747 df = 1./t_len
748 flen = int(self.f_final[i] / df) + 1
749 aLIGOpsd = aLIGOZeroDetHighPower(flen, df, self.f_low[i])
750 psd = aLIGOpsd
751 h = fd_waveform(
752 {
753 "theta_jn": self.theta_jn[i], "phase": self.phase[i],
754 "phi_jl": self.phi_jl[i], "psi": self.psi_l[i],
755 "mass_1": self.mass_1[i], "mass_2": self.mass_2[i],
756 "tilt_1": self.tilt_1[i], "tilt_2": self.tilt_2[i],
757 "phi_12": self.phi_12[i], "a_1": self.a_1[i],
758 "a_2": self.a_2[i], "luminosity_distance": self.distance[i],
759 "ra": self.ra[i], "dec": self.dec[i],
760 "geocent_time": self.time[i]
761 }, self.approx, df, self.f_low[i], self.f_final[i],
762 f_ref=self.f_low[i], flen=flen, pycbc=True, project="L1"
763 )
764 harmonics = _calculate_precessing_harmonics(
765 self.mass_1[i], self.mass_2[i], self.a_1[i],
766 self.a_2[i], self.tilt_1[i], self.tilt_2[i],
767 self.phi_12[i], self.beta[i], self.distance[i],
768 approx=self.approx, f_final=self.f_final[i],
769 flen=flen, f_ref=self.f_low[i], f_low=self.f_low[i],
770 df=df, harmonics=[0, 1, 2, 3, 4]
771 )
772 dphi = _dphi(
773 self.theta_jn[i], self.phi_jl[i], self.beta[i]
774 )
775 f_plus_j, f_cross_j = Detector("L1").antenna_pattern(
776 self.ra[i], self.dec[i], self.psi_J[i], self.time[i]
777 )
778 h_all = _make_waveform_from_precessing_harmonics(
779 harmonics, self.theta_jn[i], self.phi_jl[i],
780 self.phase[i] - dphi, f_plus_j, f_cross_j
781 )
782 overlap = compute_the_overlap(
783 h, h_all, psd, low_frequency_cutoff=self.f_low[i],
784 high_frequency_cutoff=self.f_final[i], normalized=True
785 )
786 np.testing.assert_almost_equal(np.abs(overlap), 1.0, decimal=6)
787 np.testing.assert_almost_equal(np.angle(overlap), 0.0, decimal=6)
789 def test_precessing_snr(self):
790 """Test the pesummary.gw.conversions.snr.precessing_snr function
791 """
792 # Use default PSD
793 rho_p = precessing_snr(
794 self.mass_1, self.mass_2, self.beta, self.psi_J, self.a_1, self.a_2,
795 self.tilt_1, self.tilt_2, self.phi_12, self.theta_jn,
796 self.ra, self.dec, self.time, self.phi_jl, self.distance,
797 self.phase, f_low=self.f_low, spin_1z=self.spin_1z,
798 spin_2z=self.spin_2z, multi_process=1, debug=False, df=1./8
799 )
800 assert len(rho_p) == len(self.mass_1)
801 np.testing.assert_almost_equal(
802 rho_p, [
803 0.68388587, 4.44970478, 1.50271423, 16.38278551, 8.36573961,
804 10.76045275, 5.56389142, 38.75092545, 19.04936648,
805 46.08235298, 11.04476225, 22.15809274, 21.83931449,
806 0.52940244, 18.51671749, 60.36654189, 20.90566175, 7.59963947,
807 28.81494433, 4.80448464
808 ]
809 )
812class TestMultipoleSNR(TestPrecessingSNR):
813 """Test the multipole_snr conversion
814 """
815 def setup_method(self):
816 super(TestMultipoleSNR, self).setup_method()
818 @pytest.mark.skip(reason="Inherited test")
819 def test_harmonic_overlap(self):
820 pass
822 @pytest.mark.skip(reason="Inherited test")
823 def test_precessing_snr(self):
824 pass
826 def test_multipole_snr(self):
827 """Test the pesummary.gw.conversions.multipole_snr function
828 """
829 rho = multipole_snr(
830 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z, self.psi_l,
831 self.theta_jn, self.ra, self.dec, self.time, self.distance,
832 self.phase, multi_process=1, df=1./8, multipole=[21, 33, 44]
833 )
834 assert rho.shape[0] == 3
835 np.testing.assert_almost_equal(
836 rho[0], [
837 0.65751988, 2.1414448, 1.92966995, 1.36583827, 0.91202824,
838 3.76485092, 7.65392604, 1.16167184, 10.25460489, 9.33561953,
839 0.21747679, 4.44893768, 1.4904289, 4.70318423, 2.76921055,
840 1.90083035, 3.42330456, 1.12028809, 15.47861815, 0.15595173]
841 )
842 np.testing.assert_almost_equal(
843 rho[1], [
844 2.20813922, 5.19271057, 7.4048816, 3.81037071, 1.62958146,
845 4.28125682, 22.54228317, 14.46282899, 21.05195085, 26.81466547,
846 2.33247145, 11.60509714, 6.4064751, 16.00204021, 6.21510341,
847 14.83850082, 12.07577107, 0.7150659, 14.36457899, 3.02642992
848 ]
849 )
850 np.testing.assert_almost_equal(
851 rho[2], [
852 0.18102321, 2.20770255, 0.93261726, 1.7573492, 4.60563054,
853 3.2601658, 9.10279199, 13.04963478, 9.7026812, 15.3289052,
854 1.58697352, 5.36571405, 8.15103903, 4.51177119, 4.90159264,
855 11.23723456, 4.3787358, 1.02301209, 7.38082088, 9.24657289
856 ]
857 )
858 with pytest.raises(ValueError):
859 rho = multipole_snr(
860 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z, self.psi_l,
861 self.theta_jn, self.ra, self.dec, self.time, self.distance,
862 self.phase, multi_process=1, df=1./8, multipole=[21, 33, 44, 55]
863 )
866class TestNRutils(object):
868 def setup_method(self):
869 self.mass_1 = 100
870 self.mass_2 = 5
871 self.total_mass = m_total_from_m1_m2(self.mass_1, self.mass_2)
872 self.eta = eta_from_m1_m2(self.mass_1, self.mass_2)
873 self.spin_1z = 0.3
874 self.spin_2z = 0.5
875 self.chi_eff = chi_eff(
876 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z
877 )
878 self.final_spin = bbh_final_spin_non_precessing_Healyetal(
879 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z,
880 version="2016"
881 )
883 def test_bbh_peak_luminosity_non_precessing_Healyetal(self):
884 from lalinference.imrtgr.nrutils import \
885 bbh_peak_luminosity_non_precessing_Healyetal as lal_Healyetal
887 assert np.round(
888 bbh_peak_luminosity_non_precessing_Healyetal(
889 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z
890 ), 8
891 ) == np.round(
892 lal_Healyetal(
893 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z
894 ), 8
895 )
897 def test_bbh_peak_luminosity_non_precessing_T1600018(self):
898 from lalinference.imrtgr.nrutils import \
899 bbh_peak_luminosity_non_precessing_T1600018 as lal_T1600018
901 assert np.round(
902 bbh_peak_luminosity_non_precessing_T1600018(
903 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z
904 ), 8
905 ) == np.round(
906 lal_T1600018(
907 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z
908 ), 8
909 )
911 def test_bbh_peak_luminosity_non_precessing_UIB2016(self):
912 from lalinference.imrtgr.nrutils import \
913 bbh_peak_luminosity_non_precessing_UIB2016 as lal_UIB2016
915 assert np.round(
916 bbh_peak_luminosity_non_precessing_UIB2016(
917 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z
918 ), 8
919 ) == np.round(
920 lal_UIB2016(
921 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z
922 ), 8
923 )
925 def test_bbh_final_spin_non_precessing_Healyetal(self):
926 from lalinference.imrtgr.nrutils import \
927 bbh_final_spin_non_precessing_Healyetal as lal_Healyetal
929 assert np.round(self.final_spin, 8) == np.round(
930 lal_Healyetal(
931 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z,
932 version="2016"
933 ), 8
934 )
936 def test_bbh_final_mass_non_precessing_Healyetal(self):
937 from lalinference.imrtgr.nrutils import \
938 bbh_final_mass_non_precessing_Healyetal as lal_Healyetal
940 assert np.round(
941 bbh_final_mass_non_precessing_Healyetal(
942 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z,
943 final_spin=self.final_spin, version="2016"
944 ), 8
945 ) == np.round(
946 lal_Healyetal(
947 self.mass_1, self.mass_2, self.spin_1z, self.spin_2z,
948 chif=self.final_spin, version="2016"
949 ), 8
950 )
953class TestConvert(object):
954 """Test the pesummary.gw.conversions._Conversion class
955 """
956 def setup_method(self):
957 """Setup the TestConvert class
958 """
959 self.dirs = [tmpdir]
960 for dd in self.dirs:
961 if not os.path.isdir(dd):
962 os.mkdir(dd)
964 def teardown_method(self):
965 """Remove the files and directories created from this class
966 """
967 for dd in self.dirs:
968 if os.path.isdir(dd):
969 shutil.rmtree(dd)
971 @staticmethod
972 def _convert(resume_file=None, **kwargs):
973 """
974 """
975 np.random.seed(100)
976 data = {
977 "mass_1": np.random.uniform(10, 1000, 10000),
978 "mass_2": np.random.uniform(2, 10, 10000),
979 "spin_1z": np.random.uniform(0, 1, 10000),
980 "spin_2z": np.random.uniform(0, 1, 10000)
981 }
982 converted = convert(data, resume_file=resume_file, **kwargs)
983 return converted
985 def test_from_checkpoint(self):
986 """Check that when restarted from checkpoint, the output is the same
987 """
988 import time
989 import multiprocessing
990 from pesummary.io import read
992 t0 = time.time()
993 no_checkpoint = self._convert()
994 t1 = time.time()
995 # check that when interrupted and restarted, the results are the same
996 process = multiprocessing.Process(
997 target=self._convert, args=["{}/checkpoint.pickle".format(tmpdir)]
998 )
999 process.start()
1000 time.sleep(5)
1001 process.terminate()
1002 # check that not all samples have been made
1003 _checkpoint = read("{}/checkpoint.pickle".format(tmpdir), checkpoint=True)
1004 assert os.path.isfile("{}/checkpoint.pickle".format(tmpdir))
1005 assert not all(
1006 param in _checkpoint.parameters for param in no_checkpoint.keys()
1007 )
1008 # restart from checkpoint
1009 checkpoint = self._convert(resume_file="{}/checkpoint.pickle".format(tmpdir))
1010 for param, value in no_checkpoint.items():
1011 np.testing.assert_almost_equal(value, checkpoint[param])
1014def test_evolve_angles_forwards():
1015 """Check that the pesummary.gw.conversions.evolve.evolve_angles_forwards
1016 function works as expected
1017 """
1018 from pesummary.gw.conversions.evolve import (
1019 evolve_spins, _wrapper_for_evolve_angles_forwards
1020 )
1021 from lalsimulation import (
1022 SimInspiralSpinTaylorPNEvolveOrbit, SIM_INSPIRAL_SPINS_FLOW
1023 )
1024 from lal import MTSUN_SI, MSUN_SI
1025 import lalsimulation
1026 input_data = [
1027 (
1028 146.41677334700145, 123.72830811599943, 0.4788630899579355,
1029 0.4486656607260327, 2.121165143160338, 0.5194005460918241,
1030 3.0144238369366736, 11., 11., "IMRPhenomPv3HM"
1031 ),
1032 (
1033 145.60646217583144, 97.70678957171464, 0.418373390477266,
1034 0.22414039975174402, 2.2857994587400494, 2.730311388309907,
1035 3.47318438014925, 11., 11., "IMRPhenomPv3HM"
1036 ),
1037 (
1038 144.27815417432444, 99.32482850107179, 0.3413842190485782,
1039 0.12003981467617035, 2.4429395586527884, 0.9993057630904596,
1040 1.422769967575501, 11., 11., "IMRPhenomPv3HM"
1041 ),
1042 ]
1043 for num, sample in enumerate(input_data):
1044 tilt_1_evol, tilt_2_evol, phi_12_evol = evolve_spins(
1045 *sample, multi_process=1, evolve_limit="ISCO", dt=0.1,
1046 evolution_approximant="SpinTaylorT5"
1047 )
1048 _tilt_1_evol, _tilt_2_evol, _phi_12_evol = _wrapper_for_evolve_angles_forwards(
1049 list(sample)[:7] + [11., 6 ** -0.5, 1e-3, 0.1, "SpinTaylorT5"]
1050 )
1051 np.testing.assert_almost_equal(tilt_1_evol, _tilt_1_evol)
1052 np.testing.assert_almost_equal(tilt_2_evol, _tilt_2_evol)
1053 np.testing.assert_almost_equal(phi_12_evol, _phi_12_evol)
1054 total_mass = (sample[0] + sample[1]) * MTSUN_SI
1055 spinfreq_enum = SimInspiralGetSpinFreqFromApproximant(
1056 getattr(lalsimulation, sample[-1])
1057 )
1058 f_start = float(np.where(
1059 np.array(spinfreq_enum == SIM_INSPIRAL_SPINS_FLOW), sample[-3],
1060 sample[-2]
1061 ))
1062 f_final = (6 ** -0.5) ** 3 / (total_mass * np.pi)
1063 _approx = lalsimulation.SpinTaylorT5
1064 data = SimInspiralSpinTaylorPNEvolveOrbit(
1065 deltaT=0.1 * total_mass, m1=sample[0] * MSUN_SI,
1066 m2=sample[1] * MSUN_SI, fStart=f_start, fEnd=f_final,
1067 s1x=sample[2] * np.sin(sample[4]), s1y=0.,
1068 s1z=sample[2] * np.cos(sample[4]),
1069 s2x=sample[3] * np.sin(sample[5]) * np.cos(sample[6]),
1070 s2y=sample[3] * np.sin(sample[5]) * np.sin(sample[6]),
1071 s2z=sample[3] * np.cos(sample[5]), lnhatx=0., lnhaty=0., lnhatz=1.,
1072 e1x=1., e1y=0., e1z=0., lambda1=0., lambda2=0., quadparam1=1.,
1073 quadparam2=1., spinO=6, tideO=0, phaseO=7, lscorr=0,
1074 approx=_approx
1075 )
1076 a_1_evolve = np.array(
1077 [
1078 data[2].data.data[-1], data[3].data.data[-1],
1079 data[4].data.data[-1]
1080 ]
1081 )
1082 a_2_evolve = np.array(
1083 [
1084 data[5].data.data[-1], data[6].data.data[-1],
1085 data[7].data.data[-1]
1086 ]
1087 )
1088 Ln_evolve = np.array(
1089 [
1090 data[8].data.data[-1], data[9].data.data[-1],
1091 data[10].data.data[-1]
1092 ]
1093 )
1094 # code taken from https://git.ligo.org/lscsoft/lalsuite/-/blob/master/
1095 # lalinference/bin/lalinference_evolve_spins_and_append_samples.py#L26-53
1096 # to test
1097 # pesummary.gw.conversions.tilt_angles_and_phi_12_from_spin_vectors_and_L
1098 chi1_v_norm = np.linalg.norm(a_1_evolve)
1099 chi2_v_norm = np.linalg.norm(a_2_evolve)
1100 Ln_evolve /= np.linalg.norm(Ln_evolve)
1101 chi1dL_v = np.dot(a_1_evolve, Ln_evolve)
1102 chi2dL_v = np.dot(a_2_evolve, Ln_evolve)
1103 chi1inplane = a_1_evolve - chi1dL_v * Ln_evolve
1104 chi2inplane = a_2_evolve - chi2dL_v * Ln_evolve
1105 cos_tilt1 = chi1dL_v / chi1_v_norm
1106 cos_tilt2 = chi2dL_v / chi2_v_norm
1107 cos_phi12 = np.dot(chi1inplane, chi2inplane) / (
1108 np.linalg.norm(chi1inplane) * np.linalg.norm(chi2inplane)
1109 )
1110 phi12_evol_i = np.arccos(cos_phi12)
1111 if np.sign(np.dot(Ln_evolve, np.cross(a_1_evolve, a_2_evolve))) < 0:
1112 phi12_evol_i = 2. * np.pi - phi12_evol_i
1113 tilt_1_evol_true = np.arccos(cos_tilt1)
1114 tilt_2_evol_true = np.arccos(cos_tilt2)
1115 phi_12_evol_true = phi12_evol_i
1116 np.testing.assert_almost_equal(tilt_1_evol, tilt_1_evol_true)
1117 np.testing.assert_almost_equal(tilt_2_evol, tilt_2_evol_true)
1118 np.testing.assert_almost_equal(phi_12_evol, phi_12_evol_true)
1121def test_evolve_angles_backwards():
1122 """Check that the pesummary.gw.conversions.evolve.evolve_angles_backwards
1123 function works as expected
1124 """
1125 from pesummary.gw.conversions.evolve import evolve_spins
1126 from lal import MSUN_SI
1127 from packaging import version
1128 import lalsimulation
1129 input_data = [
1130 (
1131 1.3862687342652575e+32, 1.5853186050191907e+31, 0.8768912154180827,
1132 0.9635416612042661, 2.8861591668037119, 2.7423707262813442,
1133 4.750537251642867, 8.000000000000000, 8.000000000000000, "IMRPhenomXPHM",
1134 2.8861301463160993, 2.7425208816155378, 2.8861542118177956,
1135 2.7426347054696230, 'v1'
1136 ),
1137 (
1138 4.0380177255695994e+31, 2.1111685497317552e+31, 0.9442047756726544,
1139 0.2197148251155545, 2.7060072810080551, 0.8920951236808333,
1140 1.7330264974887994, 14.0000000000000, 14.0000000000000, "IMRPhenomXPHM",
1141 2.7082295817672812, 0.8821772303625787, 2.7084623305332132,
1142 0.8811491121799003, 'v1'
1143 ),
1144 (
1145 1.4778236544770486e+32, 2.6197742077777032e+31, 0.4650532384488123,
1146 0.4135203147241133, 2.5477872046486589, 1.3374887745402186,
1147 5.8300235171959054, 15.0000000000000, 15.0000000000000, "IMRPhenomXPHM",
1148 2.5307614455226255, 1.3999636874375283, 2.5310639813744329,
1149 1.4020896531141123, 'v1'
1150 )
1151 ]
1152 for num, method in enumerate(["precession_averaged", "hybrid_orbit_averaged"]):
1153 for sample in input_data:
1154 _sample = np.array(sample[:9])
1155 _sample[:2] = _sample[:2] / MSUN_SI
1156 tilt_1_inf, tilt_2_inf, phi_12 = evolve_spins(
1157 *_sample, sample[9], method=method, multi_process=1,
1158 evolve_limit="infinite_separation", version=sample[-1],
1159 approx="SpinTaylorT5"
1160 )
1161 if num == 0:
1162 np.testing.assert_almost_equal(tilt_1_inf, sample[10], 5)
1163 np.testing.assert_almost_equal(tilt_2_inf, sample[11], 5)
1164 else:
1165 np.testing.assert_almost_equal(tilt_1_inf, sample[12], 5)
1166 np.testing.assert_almost_equal(tilt_2_inf, sample[13], 5)