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 

2 

3import os 

4import socket 

5import shutil 

6 

7import numpy as np 

8import h5py 

9 

10from pesummary.gw.conversions import * 

11from pesummary.gw.conversions.nrutils import * 

12from pycbc import conversions 

13import pytest 

14import tempfile 

15 

16tmpdir = tempfile.TemporaryDirectory(prefix=".", dir=".").name 

17 

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

19 

20 

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 

31 

32 

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. 

52 

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.] 

62 

63 redshift = 0.5 

64 l_distance = 500. 

65 

66 self.opts = Arguments() 

67 

68 def test_z_from_dL(self): 

69 from bilby.gw.conversion import luminosity_distance_to_redshift 

70 

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 ) 

82 

83 def test_change_of_cosmology_for_z_from_dL(self): 

84 from lalinference.bayespputils import calculate_redshift 

85 

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) 

94 

95 def test_z_from_comoving_volume(self): 

96 from pycbc.cosmology import redshift_from_comoving_volume 

97 

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 ) 

109 

110 def test_dL_from_z(self): 

111 from bilby.gw.conversion import redshift_to_luminosity_distance 

112 

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 ) 

119 

120 def test_comoving_volume_from_z(self): 

121 from pycbc.cosmology import cosmological_quantity_from_redshift 

122 

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 

132 

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 ) 

139 

140 def test_m1_source_from_m1_z(self): 

141 from bilby.gw.conversion import generate_source_frame_parameters 

142 

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 

155 

156 def test_m1_m2_from_m1_source_m2_source_z(self): 

157 from bilby.gw.conversion import generate_source_frame_parameters 

158 

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) 

163 

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) 

178 

179 def test_m2_source_from_m2_z(self): 

180 from bilby.gw.conversion import generate_source_frame_parameters 

181 

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 

194 

195 def test_m_total_source_from_mtotal_z(self): 

196 from bilby.gw.conversion import generate_source_frame_parameters 

197 

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 

208 

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 

213 

214 def test_mchirp_source_from_mchirp_z(self): 

215 from bilby.gw.conversion import generate_source_frame_parameters 

216 

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 

229 

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 

234 

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 ) 

243 

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 ) 

252 

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 ) 

262 

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 ) 

271 

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 ) 

280 

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 ) 

289 

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 ) 

298 

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 ) 

307 

308 def test_q_from_eta(self): 

309 from bilby.gw.conversion import symmetric_mass_ratio_to_mass_ratio 

310 

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 ) 

317 

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 

322 

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 ) 

340 

341 from lalsimulation import SimPhenomUtilsChiP 

342 

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 

365 

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 ) 

377 

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

387 

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) 

393 

394 def test_spin_angles(self): 

395 from lalsimulation import SimInspiralTransformPrecessingWvf2PE 

396 

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 ) 

425 

426 def test_component_spins(self): 

427 from bilby.gw.conversion import bilby_to_lalsimulation_spins 

428 from lal import MSUN_SI 

429 

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) 

441 

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 ) 

454 

455 def test_time_in_each_ifo(self): 

456 from pycbc.detector import Detector 

457 from lal import TimeDelayFromEarthCenter 

458 from lalsimulation import DetectorPrefixToLALDetector 

459 

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) 

464 

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 

480 

481 

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 ) 

493 

494 def test_delta_lambda_from_lambda1_lambda2(self): 

495 from bilby.gw.conversion import lambda_1_lambda_2_to_delta_lambda_tilde 

496 

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 ) 

507 

508 def test_lambda1_from_lambda_tilde(self): 

509 from bilby.gw.conversion import lambda_tilde_to_lambda_1_lambda_2 

510 

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 

520 

521 def test_lambda2_from_lambda1(self): 

522 from bilby.gw.conversion import convert_to_lal_binary_neutron_star_parameters 

523 

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 

533 

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 

542 

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) 

572 

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 

578 

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 

673 

674 def test_remove_parameter(self): 

675 from pesummary.gw.conversions import convert 

676 

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 

691 

692 

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) 

727 

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 

740 

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) 

788 

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 ) 

810 

811 

812class TestMultipoleSNR(TestPrecessingSNR): 

813 """Test the multipole_snr conversion 

814 """ 

815 def setup_method(self): 

816 super(TestMultipoleSNR, self).setup_method() 

817 

818 @pytest.mark.skip(reason="Inherited test") 

819 def test_harmonic_overlap(self): 

820 pass 

821 

822 @pytest.mark.skip(reason="Inherited test") 

823 def test_precessing_snr(self): 

824 pass 

825 

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 ) 

864 

865 

866class TestNRutils(object): 

867 

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 ) 

882 

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 

886 

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 ) 

896 

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 

900 

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 ) 

910 

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 

914 

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 ) 

924 

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 

928 

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 ) 

935 

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 

939 

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 ) 

951 

952 

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) 

963 

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) 

970 

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 

984 

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 

991 

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

1012 

1013 

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) 

1119 

1120 

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)