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