Coverage for pesummary/gw/conversions/__init__.py: 58.9%

1210 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 

3error_msg = ( 

4 "Unable to install '{}'. You will not be able to use some of the inbuilt " 

5 "functions." 

6) 

7import copy 

8import numpy as np 

9from pathlib import Path 

10 

11from pesummary import conf 

12from pesummary.utils.decorators import set_docstring 

13from pesummary.utils.exceptions import EvolveSpinError 

14from pesummary.utils.utils import logger 

15 

16try: 

17 import lalsimulation 

18except ImportError: 

19 logger.warning(error_msg.format("lalsimulation")) 

20try: 

21 import astropy 

22except ImportError: 

23 logger.warning(error_msg.format("astropy")) 

24 

25from .angles import * 

26from .cosmology import * 

27from .cosmology import _source_from_detector 

28from .mass import * 

29from .remnant import * 

30from .remnant import _final_from_initial_BBH 

31from .snr import * 

32from .snr import _ifo_snr 

33from .spins import * 

34from .tidal import * 

35from .tidal import _check_NSBH_approximant 

36from .time import * 

37 

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

39_conversion_doc = """ 

40 Class to calculate all possible derived quantities 

41 

42 Parameters 

43 ---------- 

44 data: dict, list 

45 either a dictionary or samples or a list of parameters and a list of 

46 samples. See the examples below for details 

47 extra_kwargs: dict, optional 

48 dictionary of kwargs associated with this set of posterior samples. 

49 f_low: float, optional 

50 the low frequency cut-off to use when evolving the spins 

51 f_start: float, optional 

52 the starting frequency of the waveform 

53 f_ref: float, optional 

54 the reference frequency when spins are defined 

55 f_final: float, optional 

56 the final frequency to use when integrating over frequencies 

57 approximant: str, optional 

58 the approximant to use when evolving the spins 

59 approximant_flags: dict, optional 

60 dictionary of flags to control the variant of waveform model used 

61 evolve_spins_forwards: float/str, optional 

62 the final velocity to evolve the spins up to. 

63 evolve_spins_backwards: str, optional 

64 method to use when evolving the spins backwards to an infinite separation 

65 return_kwargs: Bool, optional 

66 if True, return a modified dictionary of kwargs containing information 

67 about the conversion 

68 NRSur_fits: float/str, optional 

69 the NRSurrogate model to use to calculate the remnant fits. If nothing 

70 passed, the average NR fits are used instead 

71 multipole_snr: Bool, optional 

72 if True, the SNR in the (l, m) = [(2, 1), (3, 3), (4, 4)] multipoles are 

73 calculated from the posterior samples. 

74 samples. 

75 precessing_snr: Bool, optional 

76 if True, the precessing SNR is calculated from the posterior samples. 

77 psd: dict, optional 

78 dictionary containing a psd frequency series for each detector you wish 

79 to include in calculations 

80 waveform_fits: Bool, optional 

81 if True, the approximant is used to calculate the remnant fits. Default 

82 is False which means that the average NR fits are used 

83 multi_process: int, optional 

84 number of cores to use to parallelize the computationally expensive 

85 conversions 

86 redshift_method: str, optional 

87 method you wish to use when calculating the redshift given luminosity 

88 distance samples. If redshift samples already exist, this method is not 

89 used. Default is 'approx' meaning that interpolation is used to calculate 

90 the redshift given N luminosity distance points. 

91 cosmology: str, optional 

92 cosmology you wish to use when calculating the redshift given luminosity 

93 distance samples. 

94 force_non_evolved: Bool, optional 

95 force non evolved remnant quantities to be calculated when evolved quantities 

96 already exist in the input. Default False 

97 force_BBH_remnant_computation: Bool, optional 

98 force BBH remnant quantities to be calculated for systems that include 

99 tidal deformability parameters where BBH fits may not be applicable. 

100 Default False. 

101 force_BH_spin_evolution: Bool, optional 

102 force BH spin evolution methods to be applied for systems that include 

103 tidal deformability parameters where these methods may not be applicable. 

104 Default False. 

105 disable_remnant: Bool, optional 

106 disable all remnant quantities from being calculated. Default False. 

107 add_zero_spin: Bool, optional 

108 if no spins are present in the posterior table, add spins with 0 value. 

109 Default False. 

110 psd_default: str/pycbc.psd obj, optional 

111 Default PSD to use for conversions when no other PSD is provided. 

112 regenerate: list, optional 

113 list of posterior distributions that you wish to regenerate 

114 return_dict: Bool, optional 

115 if True, return a pesummary.utils.utils.SamplesDict object 

116 resume_file: str, optional 

117 path to file to use for checkpointing. If not provided, checkpointing 

118 is not used. Default None 

119 

120 Examples 

121 -------- 

122 There are two ways of passing arguments to this conversion class, either 

123 a dictionary of samples or a list of parameters and a list of samples. See 

124 the examples below: 

125 

126 >>> samples = {"mass_1": 10, "mass_2": 5} 

127 >>> converted_samples = %(function)s(samples) 

128 

129 >>> parameters = ["mass_1", "mass_2"] 

130 >>> samples = [10, 5] 

131 >>> converted_samples = %(function)s(parameters, samples) 

132 

133 >>> samples = {"mass_1": [10, 20], "mass_2": [5, 8]} 

134 >>> converted_samples = %(function)s(samples) 

135 

136 >>> parameters = ["mass_1", "mass_2"] 

137 >>> samples = [[10, 5], [20, 8]] 

138 """ 

139 

140 

141@set_docstring(_conversion_doc % {"function": "convert"}) 

142def convert(*args, restart_from_checkpoint=False, resume_file=None, **kwargs): 

143 import os 

144 if resume_file is not None: 

145 if os.path.isfile(resume_file) and restart_from_checkpoint: 

146 return _Conversion.load_current_state(resume_file) 

147 logger.info( 

148 "Unable to find resume file for conversion. Not restarting from " 

149 "checkpoint" 

150 ) 

151 return _Conversion(*args, resume_file=resume_file, **kwargs) 

152 

153 

154class _PickledConversion(object): 

155 pass 

156 

157 

158@set_docstring(_conversion_doc % {"function": "_Conversion"}) 

159class _Conversion(object): 

160 @classmethod 

161 def load_current_state(cls, resume_file): 

162 """Load current state from a resume file 

163 

164 Parameters 

165 ---------- 

166 resume_file: str 

167 path to a resume file to restart conversion 

168 """ 

169 from pesummary.io import read 

170 logger.info( 

171 "Reading checkpoint file: {}".format(resume_file) 

172 ) 

173 state = read(resume_file, checkpoint=True) 

174 return cls( 

175 state.parameters, state.samples, extra_kwargs=state.extra_kwargs, 

176 evolve_spins_forwards=state.evolve_spins_forwards, 

177 evolve_spins_backwards=state.evolve_spins_backwards, 

178 NRSur_fits=state.NRSurrogate, 

179 waveform_fits=state.waveform_fit, multi_process=state.multi_process, 

180 redshift_method=state.redshift_method, cosmology=state.cosmology, 

181 force_non_evolved=state.force_non_evolved, 

182 force_BBH_remnant_computation=state.force_remnant, 

183 disable_remnant=state.disable_remnant, 

184 add_zero_spin=state.add_zero_spin, regenerate=state.regenerate, 

185 return_kwargs=state.return_kwargs, return_dict=state.return_dict, 

186 resume_file=state.resume_file 

187 ) 

188 

189 def write_current_state(self): 

190 """Write the current state of the conversion class to file 

191 """ 

192 from pesummary.io import write 

193 state = _PickledConversion() 

194 for key, value in vars(self).items(): 

195 setattr(state, key, value) 

196 

197 _path = Path(self.resume_file) 

198 write( 

199 state, outdir=_path.parent, file_format="pickle", 

200 filename=_path.name, overwrite=True 

201 ) 

202 logger.debug( 

203 "Written checkpoint file: {}".format(self.resume_file) 

204 ) 

205 

206 def __new__(cls, *args, **kwargs): 

207 from pesummary.utils.samples_dict import SamplesDict 

208 from pesummary.utils.parameters import Parameters 

209 

210 obj = super(_Conversion, cls).__new__(cls) 

211 base_replace = ( 

212 "'{}': {} already found in the result file. Overwriting with " 

213 "the passed {}" 

214 ) 

215 if len(args) > 2: 

216 raise ValueError( 

217 "The _Conversion module only takes as arguments a dictionary " 

218 "of samples or a list of parameters and a list of samples" 

219 ) 

220 elif isinstance(args[0], dict): 

221 parameters = Parameters(args[0].keys()) 

222 samples = np.atleast_2d( 

223 np.array([args[0][i] for i in parameters]).T 

224 ).tolist() 

225 else: 

226 if not isinstance(args[0], Parameters): 

227 parameters = Parameters(args[0]) 

228 else: 

229 parameters = args[0] 

230 samples = args[1] 

231 samples = np.atleast_2d(samples).tolist() 

232 extra_kwargs = kwargs.get("extra_kwargs", {"sampler": {}, "meta_data": {}}) 

233 f_low = kwargs.get("f_low", None) 

234 f_start = kwargs.get("f_start", None) 

235 f_ref = kwargs.get("f_ref", None) 

236 f_final = kwargs.get("f_final", None) 

237 delta_f = kwargs.get("delta_f", None) 

238 

239 for param, value in {"f_final": f_final, "delta_f": delta_f}.items(): 

240 if value is not None and param in extra_kwargs["meta_data"].keys(): 

241 logger.warning( 

242 base_replace.format( 

243 param, extra_kwargs["meta_data"][param], value 

244 ) 

245 ) 

246 extra_kwargs["meta_data"][param] = value 

247 elif value is not None: 

248 extra_kwargs["meta_data"][param] = value 

249 elif value is None and param in extra_kwargs["meta_data"].keys(): 

250 logger.debug( 

251 "Found {} in input file. Using {}Hz".format( 

252 param, extra_kwargs["meta_data"][param] 

253 ) 

254 ) 

255 else: 

256 logger.warning( 

257 "Could not find {} in input file and one was not passed " 

258 "from the command line. Using {}Hz as default".format( 

259 param, getattr(conf, "default_{}".format(param)) 

260 ) 

261 ) 

262 extra_kwargs["meta_data"][param] = getattr( 

263 conf, "default_{}".format(param) 

264 ) 

265 

266 approximant = kwargs.get("approximant", None) 

267 approximant_flags = kwargs.get("approximant_flags", {}) 

268 NRSurrogate = kwargs.get("NRSur_fits", False) 

269 redshift_method = kwargs.get("redshift_method", "approx") 

270 cosmology = kwargs.get("cosmology", "Planck15") 

271 force_non_evolved = kwargs.get("force_non_evolved", False) 

272 force_remnant = kwargs.get("force_BBH_remnant_computation", False) 

273 force_evolve = kwargs.get("force_BH_spin_evolution", False) 

274 disable_remnant = kwargs.get("disable_remnant", False) 

275 if redshift_method not in ["approx", "exact"]: 

276 raise ValueError( 

277 "'redshift_method' can either be 'approx' corresponding to " 

278 "an approximant method, or 'exact' corresponding to an exact " 

279 "method of calculating the redshift" 

280 ) 

281 if isinstance(NRSurrogate, bool) and NRSurrogate: 

282 raise ValueError( 

283 "'NRSur_fits' must be a string corresponding to the " 

284 "NRSurrogate model you wish to use to calculate the remnant " 

285 "quantities" 

286 ) 

287 waveform_fits = kwargs.get("waveform_fits", False) 

288 evolve_spins_forwards = kwargs.get("evolve_spins_forwards", False) 

289 evolve_spins_backwards = kwargs.get("evolve_spins_backwards", False) 

290 if disable_remnant and ( 

291 force_non_evolved or force_remnant 

292 or NRSurrogate or waveform_fits or evolve_spins_forwards 

293 ): 

294 _disable = [] 

295 if force_non_evolved: 

296 _disable.append("force_non_evolved") 

297 force_non_evolved = False 

298 if force_remnant: 

299 _disable.append("force_BBH_remnant_computation") 

300 force_remnant = False 

301 if NRSurrogate: 

302 _disable.append("NRSur_fits") 

303 NRSurrogate = False 

304 if waveform_fits: 

305 _disable.append("waveform_fits") 

306 waveform_fits = False 

307 if evolve_spins_forwards: 

308 _disable.append("evolve_spins_forwards") 

309 evolve_spins_forwards = False 

310 logger.warning( 

311 "Unable to use 'disable_remnant' and {}. Setting " 

312 "{} and disabling all remnant quantities from being " 

313 "calculated".format( 

314 " or ".join(_disable), 

315 " and ".join(["{}=False".format(_p) for _p in _disable]) 

316 ) 

317 ) 

318 if NRSurrogate and waveform_fits: 

319 raise ValueError( 

320 "Unable to use both the NRSurrogate and {} to calculate " 

321 "remnant quantities. Please select only one option".format( 

322 approximant 

323 ) 

324 ) 

325 if isinstance(evolve_spins_forwards, bool) and evolve_spins_forwards: 

326 raise ValueError( 

327 "'evolve_spins_forwards' must be a float, the final velocity to " 

328 "evolve the spins up to, or a string, 'ISCO', meaning " 

329 "evolve the spins up to the ISCO frequency" 

330 ) 

331 if not evolve_spins_forwards and (NRSurrogate or waveform_fits): 

332 if (approximant is not None and "eob" in approximant) or NRSurrogate: 

333 logger.warning( 

334 "Only evolved spin remnant quantities are returned by the " 

335 "{} fits.".format( 

336 "NRSurrogate" if NRSurrogate else approximant 

337 ) 

338 ) 

339 elif evolve_spins_forwards and (NRSurrogate or waveform_fits): 

340 if (approximant is not None and "eob" in approximant) or NRSurrogate: 

341 logger.warning( 

342 "The {} fits already evolve the spins. Therefore " 

343 "additional spin evolution will not be performed.".format( 

344 "NRSurrogate" if NRSurrogate else approximant 

345 ) 

346 ) 

347 else: 

348 logger.warning( 

349 "The {} fits are not applied with spin evolution.".format( 

350 approximant 

351 ) 

352 ) 

353 evolve_spins_forwards = False 

354 

355 multipole_snr = kwargs.get("multipole_snr", False) 

356 precessing_snr = kwargs.get("precessing_snr", False) 

357 

358 for freq, name in zip([f_start, f_low], ["f_start", "f_low"]): 

359 if freq is not None and name in extra_kwargs["meta_data"].keys(): 

360 logger.warning( 

361 base_replace.format( 

362 name, extra_kwargs["meta_data"][name], freq 

363 ) 

364 ) 

365 extra_kwargs["meta_data"][name] = freq 

366 elif freq is not None: 

367 extra_kwargs["meta_data"][name] = freq 

368 elif freq is None and name in extra_kwargs["meta_data"].keys(): 

369 logger.debug( 

370 "Found {} in input file. Using {}Hz".format( 

371 name, extra_kwargs["meta_data"][name] 

372 ) 

373 ) 

374 else: 

375 logger.warning( 

376 "Could not find {} in input file and one was not passed from " 

377 "the command line. Using {}Hz as default".format( 

378 name, conf.default_flow 

379 ) 

380 ) 

381 extra_kwargs["meta_data"][name] = conf.default_flow 

382 if len(approximant_flags) and "approximant_flags" in extra_kwargs["meta_data"].keys(): 

383 logger.warning( 

384 base_replace.format( 

385 "approximant_flags", extra_kwargs["meta_data"]["approximant_flags"], 

386 approximant_flags 

387 ) 

388 ) 

389 extra_kwargs["meta_data"]["approximant_flags"] = approximant_flags 

390 elif len(approximant_flags): 

391 extra_kwargs["meta_data"]["approximant_flags"] = approximant_flags 

392 if approximant is not None and "approximant" in extra_kwargs["meta_data"].keys(): 

393 logger.warning( 

394 base_replace.format( 

395 "approximant", extra_kwargs["meta_data"]["approximant"], 

396 approximant 

397 ) 

398 ) 

399 extra_kwargs["meta_data"]["approximant"] = approximant 

400 elif approximant is not None: 

401 extra_kwargs["meta_data"]["approximant"] = approximant 

402 if f_ref is not None and "f_ref" in extra_kwargs["meta_data"].keys(): 

403 logger.warning( 

404 base_replace.format( 

405 "f_ref", extra_kwargs["meta_data"]["f_ref"], f_ref 

406 ) 

407 ) 

408 extra_kwargs["meta_data"]["f_ref"] = f_ref 

409 elif f_ref is not None: 

410 extra_kwargs["meta_data"]["f_ref"] = f_ref 

411 regenerate = kwargs.get("regenerate", None) 

412 multi_process = kwargs.get("multi_process", None) 

413 if multi_process is not None: 

414 multi_process = int(multi_process) 

415 psd_default = kwargs.get("psd_default", "aLIGOZeroDetHighPower") 

416 psd = kwargs.get("psd", {}) 

417 if psd is None: 

418 psd = {} 

419 elif psd is not None and not isinstance(psd, dict): 

420 raise ValueError( 

421 "'psd' must be a dictionary of frequency series for each detector" 

422 ) 

423 ifos = list(psd.keys()) 

424 pycbc_psd = copy.deepcopy(psd) 

425 if psd != {}: 

426 from pesummary.gw.file.psd import PSD 

427 if isinstance(psd[ifos[0]], PSD): 

428 for ifo in ifos: 

429 try: 

430 pycbc_psd[ifo] = pycbc_psd[ifo].to_pycbc( 

431 extra_kwargs["meta_data"]["f_low"], 

432 f_high=extra_kwargs["meta_data"]["f_final"], 

433 f_high_override=True 

434 ) 

435 except (ImportError, IndexError, ValueError): 

436 pass 

437 obj.__init__( 

438 parameters, samples, extra_kwargs, evolve_spins_forwards, NRSurrogate, 

439 waveform_fits, multi_process, regenerate, redshift_method, 

440 cosmology, force_non_evolved, force_remnant, 

441 kwargs.get("add_zero_spin", False), disable_remnant, 

442 kwargs.get("return_kwargs", False), kwargs.get("return_dict", True), 

443 kwargs.get("resume_file", None), multipole_snr, precessing_snr, 

444 pycbc_psd, psd_default, evolve_spins_backwards, force_evolve 

445 ) 

446 return_kwargs = kwargs.get("return_kwargs", False) 

447 if kwargs.get("return_dict", True) and return_kwargs: 

448 return [ 

449 SamplesDict(obj.parameters, np.array(obj.samples).T), 

450 obj.extra_kwargs 

451 ] 

452 elif kwargs.get("return_dict", True): 

453 return SamplesDict(obj.parameters, np.array(obj.samples).T) 

454 elif return_kwargs: 

455 return obj.parameters, obj.samples, obj.extra_kwargs 

456 else: 

457 return obj.parameters, obj.samples 

458 

459 def __init__( 

460 self, parameters, samples, extra_kwargs, evolve_spins_forwards, NRSurrogate, 

461 waveform_fits, multi_process, regenerate, redshift_method, 

462 cosmology, force_non_evolved, force_remnant, add_zero_spin, 

463 disable_remnant, return_kwargs, return_dict, resume_file, multipole_snr, 

464 precessing_snr, psd, psd_default, evolve_spins_backwards, force_evolve 

465 ): 

466 self.parameters = parameters 

467 self.samples = samples 

468 self.extra_kwargs = extra_kwargs 

469 self.evolve_spins_forwards = evolve_spins_forwards 

470 self.evolve_spins_backwards = evolve_spins_backwards 

471 self.NRSurrogate = NRSurrogate 

472 self.waveform_fit = waveform_fits 

473 self.multi_process = multi_process 

474 self.regenerate = regenerate 

475 self.redshift_method = redshift_method 

476 self.cosmology = cosmology 

477 self.force_non_evolved = force_non_evolved 

478 self.force_remnant = force_remnant 

479 self.force_evolve = force_evolve 

480 self.disable_remnant = disable_remnant 

481 self.return_kwargs = return_kwargs 

482 self.return_dict = return_dict 

483 self.resume_file = resume_file 

484 self.multipole_snr = multipole_snr 

485 self.precessing_snr = precessing_snr 

486 self.psd = psd 

487 self.psd_default = psd_default 

488 self.non_precessing = False 

489 cond1 = any( 

490 param in self.parameters for param in 

491 conf.precessing_angles + conf.precessing_spins 

492 ) 

493 if not cond1: 

494 self.non_precessing = True 

495 if "chi_p" in self.parameters: 

496 _chi_p = self.specific_parameter_samples(["chi_p"]) 

497 if not np.any(_chi_p): 

498 logger.info( 

499 "chi_p = 0 for all samples. Treating this as a " 

500 "non-precessing system" 

501 ) 

502 self.non_precessing = True 

503 elif all(param in self.parameters for param in conf.precessing_spins): 

504 samples = self.specific_parameter_samples(conf.precessing_spins) 

505 if not any(np.array(samples).flatten()): 

506 logger.info( 

507 "in-plane spins are zero for all samples. Treating this as " 

508 "a non-precessing system" 

509 ) 

510 cond1 = self.non_precessing and evolve_spins_forwards 

511 cond2 = self.non_precessing and evolve_spins_backwards 

512 if cond1 or cond2: 

513 logger.info( 

514 "Spin evolution is trivial for a non-precessing system. No additional " 

515 "transformation required." 

516 ) 

517 self.evolve_spins_forwards = False 

518 self.evolve_spins_backwards = False 

519 if not self.non_precessing and multipole_snr: 

520 logger.warning( 

521 "The calculation for computing the SNR in subdominant " 

522 "multipoles assumes that the system is non-precessing. " 

523 "Since precessing samples are provided, this may give incorrect " 

524 "results" 

525 ) 

526 if self.non_precessing and precessing_snr: 

527 logger.info( 

528 "Precessing SNR is 0 for a non-precessing system. No additional " 

529 "conversion required." 

530 ) 

531 self.precessing_snr = False 

532 self.has_tidal = self._check_for_tidal_parameters() 

533 self.NSBH = self._check_for_NSBH_system() 

534 self.compute_remnant = not self.disable_remnant 

535 if self.has_tidal: 

536 if force_evolve and (self.evolve_spins_forwards or self.evolve_spins_backwards): 

537 logger.warning( 

538 "Posterior samples for tidal deformability found in the " 

539 "posterior table. 'force_evolve' provided so using BH spin " 

540 "evolution methods for this system. This may not give " 

541 "sensible results" 

542 ) 

543 elif self.evolve_spins_forwards or self.evolve_spins_backwards: 

544 logger.warning( 

545 "Tidal deformability parameters found in the posterior table. " 

546 "Skipping spin evolution as current methods are only valid " 

547 "for BHs." 

548 ) 

549 self.evolve_spins_forwards = False 

550 self.evolve_spins_backwards = False 

551 

552 if force_remnant and self.NSBH and self.compute_remnant: 

553 logger.warning( 

554 "Posterior samples for lambda_2 found in the posterior table " 

555 "but unable to find samples for lambda_1. Assuming this " 

556 "is an NSBH system. 'force_remnant' provided so using BBH remnant " 

557 "fits for this system. This may not give sensible results" 

558 ) 

559 elif self.NSBH and self.compute_remnant: 

560 logger.warning( 

561 "Posterior samples for lambda_2 found in the posterior table " 

562 "but unable to find samples for lambda_1. Applying NSBH " 

563 "fits to this system." 

564 ) 

565 self.waveform_fit = True 

566 elif force_remnant and self.compute_remnant: 

567 logger.warning( 

568 "Posterior samples for tidal deformability found in the " 

569 "posterior table. Applying BBH remnant fits to this system. " 

570 "This may not give sensible results." 

571 ) 

572 elif self.compute_remnant: 

573 logger.info( 

574 "Skipping remnant calculations as tidal deformability " 

575 "parameters found in the posterior table." 

576 ) 

577 self.compute_remnant = False 

578 if self.regenerate is not None: 

579 for param in self.regenerate: 

580 self.remove_posterior(param) 

581 self.add_zero_spin = add_zero_spin 

582 self.generate_all_posterior_samples() 

583 

584 def _check_for_tidal_parameters(self): 

585 """Check to see if any tidal parameters are stored in the table 

586 """ 

587 from pesummary.gw.file.standard_names import tidal_params 

588 

589 if any(param in self.parameters for param in tidal_params): 

590 if not all(_ in self.parameters for _ in ["lambda_1", "lambda_2"]): 

591 return True 

592 _tidal_posterior = self.specific_parameter_samples( 

593 ["lambda_1", "lambda_2"] 

594 ) 

595 if not all(np.any(_) for _ in _tidal_posterior): 

596 logger.warning( 

597 "Tidal deformability parameters found in the posterior " 

598 "table but they are all exactly 0. Assuming this is a BBH " 

599 "system." 

600 ) 

601 return False 

602 return True 

603 return False 

604 

605 def _check_for_NSBH_system(self): 

606 """Check to see if the posterior samples correspond to an NSBH 

607 system 

608 """ 

609 if "lambda_2" in self.parameters and "lambda_1" not in self.parameters: 

610 _lambda_2 = self.specific_parameter_samples(["lambda_2"]) 

611 if not np.any(_lambda_2): 

612 logger.warning( 

613 "Posterior samples for lambda_2 found in the posterior " 

614 "table but they are all exactly 0. Assuming this is a BBH " 

615 "system." 

616 ) 

617 return False 

618 return True 

619 elif "lambda_2" in self.parameters and "lambda_1" in self.parameters: 

620 _lambda_1, _lambda_2 = self.specific_parameter_samples( 

621 ["lambda_1", "lambda_2"] 

622 ) 

623 if not np.any(_lambda_1) and not np.any(_lambda_2): 

624 return False 

625 elif not np.any(_lambda_1): 

626 logger.warning( 

627 "Posterior samples for lambda_1 and lambda_2 found in the " 

628 "posterior table but lambda_1 is always exactly 0. " 

629 "Assuming this is an NSBH system." 

630 ) 

631 return True 

632 return False 

633 

634 def remove_posterior(self, parameter): 

635 if parameter in self.parameters: 

636 logger.info( 

637 "Removing the posterior samples for '{}'".format(parameter) 

638 ) 

639 ind = self.parameters.index(parameter) 

640 self.parameters.remove(self.parameters[ind]) 

641 for i in self.samples: 

642 del i[ind] 

643 else: 

644 logger.info( 

645 "'{}' is not in the table of posterior samples. Unable to " 

646 "remove".format(parameter) 

647 ) 

648 

649 def _specific_parameter_samples(self, param): 

650 """Return the samples for a specific parameter 

651 

652 Parameters 

653 ---------- 

654 param: str 

655 the parameter that you would like to return the samples for 

656 """ 

657 if param == "empty": 

658 return np.array(np.zeros(len(self.samples))) 

659 ind = self.parameters.index(param) 

660 samples = np.array([i[ind] for i in self.samples]) 

661 return samples 

662 

663 def specific_parameter_samples(self, param): 

664 """Return the samples for either a list or a single parameter 

665 

666 Parameters 

667 ---------- 

668 param: list/str 

669 the parameter/parameters that you would like to return the samples 

670 for 

671 """ 

672 if type(param) == list: 

673 samples = [self._specific_parameter_samples(i) for i in param] 

674 else: 

675 samples = self._specific_parameter_samples(param) 

676 return samples 

677 

678 def append_data(self, parameter, samples): 

679 """Add a list of samples to the existing samples data object 

680 

681 Parameters 

682 ---------- 

683 parameter: str 

684 the name of the parameter you would like to append 

685 samples: list 

686 the list of samples that you would like to append 

687 """ 

688 if parameter not in self.parameters: 

689 self.parameters.append(parameter) 

690 for num, i in enumerate(self.samples): 

691 self.samples[num].append(samples[num]) 

692 if self.resume_file is not None: 

693 self.write_current_state() 

694 

695 def _mchirp_from_mchirp_source_z(self): 

696 samples = self.specific_parameter_samples(["chirp_mass_source", "redshift"]) 

697 chirp_mass = mchirp_from_mchirp_source_z(samples[0], samples[1]) 

698 self.append_data("chirp_mass", chirp_mass) 

699 

700 def _q_from_eta(self): 

701 samples = self.specific_parameter_samples("symmetric_mass_ratio") 

702 mass_ratio = q_from_eta(samples) 

703 self.append_data("mass_ratio", mass_ratio) 

704 

705 def _q_from_m1_m2(self): 

706 samples = self.specific_parameter_samples(["mass_1", "mass_2"]) 

707 mass_ratio = q_from_m1_m2(samples[0], samples[1]) 

708 self.append_data("mass_ratio", mass_ratio) 

709 

710 def _invert_q(self): 

711 ind = self.parameters.index("mass_ratio") 

712 for num, i in enumerate(self.samples): 

713 self.samples[num][ind] = 1. / self.samples[num][ind] 

714 

715 def _invq_from_q(self): 

716 samples = self.specific_parameter_samples("mass_ratio") 

717 inverted_mass_ratio = invq_from_q(samples) 

718 self.append_data("inverted_mass_ratio", inverted_mass_ratio) 

719 

720 def _mchirp_from_mtotal_q(self): 

721 samples = self.specific_parameter_samples(["total_mass", "mass_ratio"]) 

722 chirp_mass = mchirp_from_mtotal_q(samples[0], samples[1]) 

723 self.append_data("chirp_mass", chirp_mass) 

724 

725 def _m1_from_mchirp_q(self): 

726 samples = self.specific_parameter_samples(["chirp_mass", "mass_ratio"]) 

727 mass_1 = m1_from_mchirp_q(samples[0], samples[1]) 

728 self.append_data("mass_1", mass_1) 

729 

730 def _m2_from_mchirp_q(self): 

731 samples = self.specific_parameter_samples(["chirp_mass", "mass_ratio"]) 

732 mass_2 = m2_from_mchirp_q(samples[0], samples[1]) 

733 self.append_data("mass_2", mass_2) 

734 

735 def _m1_from_mtotal_q(self): 

736 samples = self.specific_parameter_samples(["total_mass", "mass_ratio"]) 

737 mass_1 = m1_from_mtotal_q(samples[0], samples[1]) 

738 self.append_data("mass_1", mass_1) 

739 

740 def _m2_from_mtotal_q(self): 

741 samples = self.specific_parameter_samples(["total_mass", "mass_ratio"]) 

742 mass_2 = m2_from_mtotal_q(samples[0], samples[1]) 

743 self.append_data("mass_2", mass_2) 

744 

745 def _reference_frequency(self): 

746 nsamples = len(self.samples) 

747 extra_kwargs = self.extra_kwargs["meta_data"] 

748 if extra_kwargs != {} and "f_ref" in list(extra_kwargs.keys()): 

749 self.append_data( 

750 "reference_frequency", [float(extra_kwargs["f_ref"])] * nsamples 

751 ) 

752 else: 

753 logger.warning( 

754 "Could not find reference_frequency in input file. Using 20Hz " 

755 "as default") 

756 self.append_data("reference_frequency", [20.] * nsamples) 

757 

758 def _mtotal_from_m1_m2(self): 

759 samples = self.specific_parameter_samples(["mass_1", "mass_2"]) 

760 m_total = m_total_from_m1_m2(samples[0], samples[1]) 

761 self.append_data("total_mass", m_total) 

762 

763 def _mchirp_from_m1_m2(self): 

764 samples = self.specific_parameter_samples(["mass_1", "mass_2"]) 

765 chirp_mass = mchirp_from_m1_m2(samples[0], samples[1]) 

766 self.append_data("chirp_mass", chirp_mass) 

767 

768 def _eta_from_m1_m2(self): 

769 samples = self.specific_parameter_samples(["mass_1", "mass_2"]) 

770 eta = eta_from_m1_m2(samples[0], samples[1]) 

771 self.append_data("symmetric_mass_ratio", eta) 

772 

773 def _phi_12_from_phi1_phi2(self): 

774 samples = self.specific_parameter_samples(["phi_1", "phi_2"]) 

775 phi_12 = phi_12_from_phi1_phi2(samples[0], samples[1]) 

776 self.append_data("phi_12", phi_12) 

777 

778 def _phi1_from_spins(self): 

779 samples = self.specific_parameter_samples(["spin_1x", "spin_1y"]) 

780 phi_1 = phi1_from_spins(samples[0], samples[1]) 

781 self.append_data("phi_1", phi_1) 

782 

783 def _phi2_from_spins(self): 

784 samples = self.specific_parameter_samples(["spin_2x", "spin_2y"]) 

785 phi_2 = phi2_from_spins(samples[0], samples[1]) 

786 self.append_data("phi_2", phi_2) 

787 

788 def _spin_angles(self): 

789 _spin_angles = ["theta_jn", "phi_jl", "tilt_1", "tilt_2", "phi_12", 

790 "a_1", "a_2"] 

791 spin_angles_to_calculate = [ 

792 i for i in _spin_angles if i not in self.parameters] 

793 spin_components = [ 

794 "mass_1", "mass_2", "iota", "spin_1x", "spin_1y", "spin_1z", 

795 "spin_2x", "spin_2y", "spin_2z", "reference_frequency"] 

796 samples = self.specific_parameter_samples(spin_components) 

797 if "phase" in self.parameters: 

798 spin_components.append("phase") 

799 samples.append(self.specific_parameter_samples("phase")) 

800 else: 

801 logger.warning( 

802 "Phase it not given, we will be assuming that a " 

803 "reference phase of 0 to calculate all the spin angles" 

804 ) 

805 samples.append([0] * len(samples[0])) 

806 angles = spin_angles( 

807 samples[0], samples[1], samples[2], samples[3], samples[4], 

808 samples[5], samples[6], samples[7], samples[8], samples[9], 

809 samples[10]) 

810 

811 for i in spin_angles_to_calculate: 

812 ind = _spin_angles.index(i) 

813 data = np.array([i[ind] for i in angles]) 

814 self.append_data(i, data) 

815 

816 def _non_precessing_component_spins(self): 

817 spins = ["iota", "spin_1x", "spin_1y", "spin_1z", "spin_2x", "spin_2y", 

818 "spin_2z"] 

819 angles = ["a_1", "a_2", "theta_jn", "tilt_1", "tilt_2"] 

820 if all(i in self.parameters for i in angles): 

821 samples = self.specific_parameter_samples(angles) 

822 cond1 = all(i in [0, np.pi] for i in samples[3]) 

823 cond2 = all(i in [0, np.pi] for i in samples[4]) 

824 spins_to_calculate = [ 

825 i for i in spins if i not in self.parameters] 

826 if cond1 and cond1: 

827 spin_1x = np.array([0.] * len(samples[0])) 

828 spin_1y = np.array([0.] * len(samples[0])) 

829 spin_1z = samples[0] * np.cos(samples[3]) 

830 spin_2x = np.array([0.] * len(samples[0])) 

831 spin_2y = np.array([0.] * len(samples[0])) 

832 spin_2z = samples[1] * np.cos(samples[4]) 

833 iota = np.array(samples[2]) 

834 spin_components = [ 

835 iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z] 

836 

837 for i in spins_to_calculate: 

838 ind = spins.index(i) 

839 data = spin_components[ind] 

840 self.append_data(i, data) 

841 

842 def _component_spins(self): 

843 spins = ["iota", "spin_1x", "spin_1y", "spin_1z", "spin_2x", "spin_2y", 

844 "spin_2z"] 

845 spins_to_calculate = [ 

846 i for i in spins if i not in self.parameters] 

847 angles = [ 

848 "theta_jn", "phi_jl", "tilt_1", "tilt_2", "phi_12", "a_1", "a_2", 

849 "mass_1", "mass_2", "reference_frequency"] 

850 samples = self.specific_parameter_samples(angles) 

851 if "phase" in self.parameters: 

852 angles.append("phase") 

853 samples.append(self.specific_parameter_samples("phase")) 

854 else: 

855 logger.warning( 

856 "Phase it not given, we will be assuming that a " 

857 "reference phase of 0 to calculate all the spin angles" 

858 ) 

859 samples.append([0] * len(samples[0])) 

860 spin_components = component_spins( 

861 samples[0], samples[1], samples[2], samples[3], samples[4], 

862 samples[5], samples[6], samples[7], samples[8], samples[9], 

863 samples[10]) 

864 

865 for i in spins_to_calculate: 

866 ind = spins.index(i) 

867 data = np.array([i[ind] for i in spin_components]) 

868 self.append_data(i, data) 

869 

870 def _component_spins_from_azimuthal_and_polar_angles(self): 

871 spins = ["spin_1x", "spin_1y", "spin_1z", "spin_2x", "spin_2y", 

872 "spin_2z"] 

873 spins_to_calculate = [ 

874 i for i in spins if i not in self.parameters] 

875 angles = [ 

876 "a_1", "a_2", "a_1_azimuthal", "a_1_polar", "a_2_azimuthal", 

877 "a_2_polar"] 

878 samples = self.specific_parameter_samples(angles) 

879 spin_components = spin_angles_from_azimuthal_and_polar_angles( 

880 samples[0], samples[1], samples[2], samples[3], samples[4], 

881 samples[5]) 

882 for i in spins_to_calculate: 

883 ind = spins.index(i) 

884 data = np.array([i[ind] for i in spin_components]) 

885 self.append_data(i, data) 

886 

887 def _chi_p(self): 

888 parameters = [ 

889 "mass_1", "mass_2", "spin_1x", "spin_1y", "spin_2x", "spin_2y"] 

890 samples = self.specific_parameter_samples(parameters) 

891 chi_p_samples = chi_p( 

892 samples[0], samples[1], samples[2], samples[3], samples[4], 

893 samples[5]) 

894 self.append_data("chi_p", chi_p_samples) 

895 

896 def _chi_p_from_tilts(self, suffix=""): 

897 parameters = [ 

898 "mass_1", "mass_2", "a_1", "tilt_1{}".format(suffix), "a_2", 

899 "tilt_2{}".format(suffix) 

900 ] 

901 samples = self.specific_parameter_samples(parameters) 

902 chi_p_samples = chi_p_from_tilts( 

903 samples[0], samples[1], samples[2], samples[3], samples[4], 

904 samples[5] 

905 ) 

906 self.append_data("chi_p{}".format(suffix), chi_p_samples) 

907 

908 def _chi_p_2spin(self): 

909 parameters = [ 

910 "mass_1", "mass_2", "spin_1x", "spin_1y", "spin_2x", "spin_2y"] 

911 samples = self.specific_parameter_samples(parameters) 

912 chi_p_2spin_samples = chi_p_2spin( 

913 samples[0], samples[1], samples[2], samples[3], samples[4], 

914 samples[5]) 

915 self.append_data("chi_p_2spin", chi_p_2spin_samples) 

916 

917 def _chi_eff(self, suffix=""): 

918 parameters = [ 

919 "mass_1", "mass_2", "spin_1z{}".format(suffix), 

920 "spin_2z{}".format(suffix) 

921 ] 

922 samples = self.specific_parameter_samples(parameters) 

923 chi_eff_samples = chi_eff( 

924 samples[0], samples[1], samples[2], samples[3]) 

925 self.append_data("chi_eff{}".format(suffix), chi_eff_samples) 

926 

927 def _aligned_spin_from_magnitude_tilts( 

928 self, primary=False, secondary=False, suffix="" 

929 ): 

930 if primary: 

931 parameters = ["a_1", "tilt_1{}".format(suffix)] 

932 param_to_add = "spin_1z{}".format(suffix) 

933 elif secondary: 

934 parameters = ["a_2", "tilt_2{}".format(suffix)] 

935 param_to_add = "spin_2z{}".format(suffix) 

936 samples = self.specific_parameter_samples(parameters) 

937 spin_samples = samples[0] * np.cos(samples[1]) 

938 self.append_data(param_to_add, spin_samples) 

939 

940 def _cos_tilt_1_from_tilt_1(self): 

941 samples = self.specific_parameter_samples("tilt_1") 

942 cos_tilt_1 = np.cos(samples) 

943 self.append_data("cos_tilt_1", cos_tilt_1) 

944 

945 def _cos_tilt_2_from_tilt_2(self): 

946 samples = self.specific_parameter_samples("tilt_2") 

947 cos_tilt_2 = np.cos(samples) 

948 self.append_data("cos_tilt_2", cos_tilt_2) 

949 

950 def _viewing_angle(self): 

951 samples = self.specific_parameter_samples("theta_jn") 

952 viewing_angle = viewing_angle_from_inclination(samples) 

953 self.append_data("viewing_angle", viewing_angle) 

954 

955 def _dL_from_z(self): 

956 samples = self.specific_parameter_samples("redshift") 

957 distance = dL_from_z(samples, cosmology=self.cosmology) 

958 self.extra_kwargs["meta_data"]["cosmology"] = self.cosmology 

959 self.append_data("luminosity_distance", distance) 

960 

961 def _z_from_dL(self): 

962 samples = self.specific_parameter_samples("luminosity_distance") 

963 func = getattr(Redshift.Distance, self.redshift_method) 

964 redshift = func( 

965 samples, cosmology=self.cosmology, multi_process=self.multi_process 

966 ) 

967 self.extra_kwargs["meta_data"]["cosmology"] = self.cosmology 

968 self.append_data("redshift", redshift) 

969 

970 def _z_from_comoving_volume(self): 

971 samples = self.specific_parameter_samples("comoving_volume") 

972 func = getattr(Redshift.ComovingVolume, self.redshift_method) 

973 redshift = func( 

974 samples, cosmology=self.cosmology, multi_process=self.multi_process 

975 ) 

976 self.extra_kwargs["meta_data"]["cosmology"] = self.cosmology 

977 self.append_data("redshift", redshift) 

978 

979 def _comoving_distance_from_z(self): 

980 samples = self.specific_parameter_samples("redshift") 

981 distance = comoving_distance_from_z(samples, cosmology=self.cosmology) 

982 self.extra_kwargs["meta_data"]["cosmology"] = self.cosmology 

983 self.append_data("comoving_distance", distance) 

984 

985 def _comoving_volume_from_z(self): 

986 samples = self.specific_parameter_samples("redshift") 

987 volume = comoving_volume_from_z(samples, cosmology=self.cosmology) 

988 self.extra_kwargs["meta_data"]["cosmology"] = self.cosmology 

989 self.append_data("comoving_volume", volume) 

990 

991 def _m1_source_from_m1_z(self): 

992 samples = self.specific_parameter_samples(["mass_1", "redshift"]) 

993 mass_1_source = m1_source_from_m1_z(samples[0], samples[1]) 

994 self.append_data("mass_1_source", mass_1_source) 

995 

996 def _m1_from_m1_source_z(self): 

997 samples = self.specific_parameter_samples(["mass_1_source", "redshift"]) 

998 mass_1 = m1_from_m1_source_z(samples[0], samples[1]) 

999 self.append_data("mass_1", mass_1) 

1000 

1001 def _m2_source_from_m2_z(self): 

1002 samples = self.specific_parameter_samples(["mass_2", "redshift"]) 

1003 mass_2_source = m2_source_from_m2_z(samples[0], samples[1]) 

1004 self.append_data("mass_2_source", mass_2_source) 

1005 

1006 def _m2_from_m2_source_z(self): 

1007 samples = self.specific_parameter_samples(["mass_2_source", "redshift"]) 

1008 mass_2 = m2_from_m2_source_z(samples[0], samples[1]) 

1009 self.append_data("mass_2", mass_2) 

1010 

1011 def _mtotal_source_from_mtotal_z(self): 

1012 samples = self.specific_parameter_samples(["total_mass", "redshift"]) 

1013 total_mass_source = m_total_source_from_mtotal_z(samples[0], samples[1]) 

1014 self.append_data("total_mass_source", total_mass_source) 

1015 

1016 def _mtotal_from_mtotal_source_z(self): 

1017 samples = self.specific_parameter_samples(["total_mass_source", "redshift"]) 

1018 total_mass = mtotal_from_mtotal_source_z(samples[0], samples[1]) 

1019 self.append_data("total_mass", total_mass) 

1020 

1021 def _mchirp_source_from_mchirp_z(self): 

1022 samples = self.specific_parameter_samples(["chirp_mass", "redshift"]) 

1023 chirp_mass_source = mchirp_source_from_mchirp_z(samples[0], samples[1]) 

1024 self.append_data("chirp_mass_source", chirp_mass_source) 

1025 

1026 def _beta(self): 

1027 samples = self.specific_parameter_samples([ 

1028 "mass_1", "mass_2", "phi_jl", "tilt_1", "tilt_2", "phi_12", 

1029 "a_1", "a_2", "reference_frequency", "phase" 

1030 ]) 

1031 beta = opening_angle( 

1032 samples[0], samples[1], samples[2], samples[3], samples[4], 

1033 samples[5], samples[6], samples[7], samples[8], samples[9] 

1034 ) 

1035 self.append_data("beta", beta) 

1036 

1037 def _psi_J(self): 

1038 samples = self.specific_parameter_samples([ 

1039 "psi", "theta_jn", "phi_jl", "beta" 

1040 ]) 

1041 psi = psi_J(samples[0], samples[1], samples[2], samples[3]) 

1042 self.append_data("psi_J", psi) 

1043 

1044 def _retrieve_detectors(self): 

1045 detectors = [] 

1046 if "IFOs" in list(self.extra_kwargs["meta_data"].keys()): 

1047 detectors = self.extra_kwargs["meta_data"]["IFOs"].split(" ") 

1048 else: 

1049 for i in self.parameters: 

1050 if "optimal_snr" in i and i != "network_optimal_snr": 

1051 det = i.split("_optimal_snr")[0] 

1052 detectors.append(det) 

1053 return detectors 

1054 

1055 def _time_in_each_ifo(self): 

1056 detectors = self._retrieve_detectors() 

1057 samples = self.specific_parameter_samples(["ra", "dec", "geocent_time"]) 

1058 for i in detectors: 

1059 time = time_in_each_ifo(i, samples[0], samples[1], samples[2]) 

1060 self.append_data("%s_time" % (i), time) 

1061 

1062 def _time_delay_between_ifos(self): 

1063 from itertools import combinations 

1064 detectors = sorted(self._retrieve_detectors()) 

1065 for ifos in combinations(detectors, 2): 

1066 samples = self.specific_parameter_samples( 

1067 ["{}_time".format(_ifo) for _ifo in ifos] 

1068 ) 

1069 self.append_data( 

1070 "%s_%s_time_delay" % (ifos[0], ifos[1]), 

1071 samples[1] - samples[0] 

1072 ) 

1073 

1074 def _lambda1_from_lambda_tilde(self): 

1075 samples = self.specific_parameter_samples([ 

1076 "lambda_tilde", "mass_1", "mass_2"]) 

1077 lambda_1 = lambda1_from_lambda_tilde(samples[0], samples[1], samples[2]) 

1078 self.append_data("lambda_1", lambda_1) 

1079 

1080 def _lambda2_from_lambda1(self): 

1081 samples = self.specific_parameter_samples([ 

1082 "lambda_1", "mass_1", "mass_2"]) 

1083 lambda_2 = lambda2_from_lambda1(samples[0], samples[1], samples[2]) 

1084 self.append_data("lambda_2", lambda_2) 

1085 

1086 def _lambda_tilde_from_lambda1_lambda2(self): 

1087 samples = self.specific_parameter_samples([ 

1088 "lambda_1", "lambda_2", "mass_1", "mass_2"]) 

1089 lambda_tilde = lambda_tilde_from_lambda1_lambda2( 

1090 samples[0], samples[1], samples[2], samples[3]) 

1091 self.append_data("lambda_tilde", lambda_tilde) 

1092 

1093 def _delta_lambda_from_lambda1_lambda2(self): 

1094 samples = self.specific_parameter_samples([ 

1095 "lambda_1", "lambda_2", "mass_1", "mass_2"]) 

1096 delta_lambda = delta_lambda_from_lambda1_lambda2( 

1097 samples[0], samples[1], samples[2], samples[3]) 

1098 self.append_data("delta_lambda", delta_lambda) 

1099 

1100 def _NS_compactness_from_lambda(self, parameter="lambda_1"): 

1101 if parameter not in ["lambda_1", "lambda_2"]: 

1102 logger.warning( 

1103 "Can only use Love-compactness relation for 'lambda_1' and/or " 

1104 "'lambda_2'. Skipping conversion" 

1105 ) 

1106 return 

1107 ind = parameter.split("lambda_")[1] 

1108 samples = self.specific_parameter_samples([parameter]) 

1109 compactness = NS_compactness_from_lambda(samples[0]) 

1110 self.append_data("compactness_{}".format(ind), compactness) 

1111 self.extra_kwargs["meta_data"]["compactness_fits"] = ( 

1112 "YagiYunes2017_with_BBHlimit" 

1113 ) 

1114 

1115 def _NS_baryonic_mass(self, primary=True): 

1116 if primary: 

1117 params = ["compactness_1", "mass_1"] 

1118 else: 

1119 params = ["compactness_2", "mass_2"] 

1120 samples = self.specific_parameter_samples(params) 

1121 mass = NS_baryonic_mass(samples[0], samples[1]) 

1122 if primary: 

1123 self.append_data("baryonic_mass_1", mass) 

1124 else: 

1125 self.append_data("baryonic_mass_2", mass) 

1126 self.extra_kwargs["meta_data"]["baryonic_mass_fits"] = "Breu2016" 

1127 

1128 def _lambda1_lambda2_from_polytrope_EOS(self): 

1129 samples = self.specific_parameter_samples([ 

1130 "log_pressure", "gamma_1", "gamma_2", "gamma_3", "mass_1", "mass_2" 

1131 ]) 

1132 lambda_1, lambda_2 = \ 

1133 lambda1_lambda2_from_4_parameter_piecewise_polytrope_equation_of_state( 

1134 *samples, multi_process=self.multi_process 

1135 ) 

1136 if "lambda_1" not in self.parameters: 

1137 self.append_data("lambda_1", lambda_1) 

1138 if "lambda_2" not in self.parameters: 

1139 self.append_data("lambda_2", lambda_2) 

1140 

1141 def _lambda1_lambda2_from_spectral_decomposition_EOS(self): 

1142 samples = self.specific_parameter_samples([ 

1143 "spectral_decomposition_gamma_0", "spectral_decomposition_gamma_1", 

1144 "spectral_decomposition_gamma_2", "spectral_decomposition_gamma_3", 

1145 "mass_1", "mass_2" 

1146 ]) 

1147 lambda_1, lambda_2 = lambda1_lambda2_from_spectral_decomposition( 

1148 *samples, multi_process=self.multi_process 

1149 ) 

1150 if "lambda_1" not in self.parameters: 

1151 self.append_data("lambda_1", lambda_1) 

1152 if "lambda_2" not in self.parameters: 

1153 self.append_data("lambda_2", lambda_2) 

1154 

1155 def _ifo_snr(self): 

1156 abs_snrs = [ 

1157 i for i in self.parameters if "_matched_filter_abs_snr" in i 

1158 ] 

1159 angle_snrs = [ 

1160 i for i in self.parameters if "_matched_filter_snr_angle" in i 

1161 ] 

1162 for ifo in [snr.split("_matched_filter_abs_snr")[0] for snr in abs_snrs]: 

1163 if "{}_matched_filter_snr".format(ifo) not in self.parameters: 

1164 samples = self.specific_parameter_samples( 

1165 [ 

1166 "{}_matched_filter_abs_snr".format(ifo), 

1167 "{}_matched_filter_snr_angle".format(ifo) 

1168 ] 

1169 ) 

1170 snr = _ifo_snr(samples[0], samples[1]) 

1171 self.append_data("{}_matched_filter_snr".format(ifo), snr) 

1172 

1173 def _optimal_network_snr(self): 

1174 snrs = [i for i in self.parameters if "_optimal_snr" in i] 

1175 samples = self.specific_parameter_samples(snrs) 

1176 snr = network_snr(samples) 

1177 self.append_data("network_optimal_snr", snr) 

1178 

1179 def _matched_filter_network_snr(self): 

1180 mf_snrs = sorted([ 

1181 i for i in self.parameters if "_matched_filter_snr" in i 

1182 and "_angle" not in i and "_abs" not in i 

1183 ]) 

1184 opt_snrs = sorted([ 

1185 i for i in self.parameters if "_optimal_snr" in i and "network" not 

1186 in i 

1187 ]) 

1188 _mf_detectors = [ 

1189 param.split("_matched_filter_snr")[0] for param in mf_snrs 

1190 ] 

1191 _opt_detectors = [ 

1192 param.split("_optimal_snr")[0] for param in opt_snrs 

1193 ] 

1194 if _mf_detectors == _opt_detectors: 

1195 mf_samples = self.specific_parameter_samples(mf_snrs) 

1196 opt_samples = self.specific_parameter_samples(opt_snrs) 

1197 snr = network_matched_filter_snr(mf_samples, opt_samples) 

1198 self.append_data("network_matched_filter_snr", snr) 

1199 else: 

1200 logger.warning( 

1201 "Unable to generate 'network_matched_filter_snr' as " 

1202 "there is an inconsistency in the detector network based on " 

1203 "the 'optimal_snrs' and the 'matched_filter_snrs'. We find " 

1204 "that from the 'optimal_snrs', the detector network is: {} " 

1205 "while we find from the 'matched_filter_snrs', the detector " 

1206 "network is: {}".format(_opt_detectors, _mf_detectors) 

1207 ) 

1208 

1209 def _rho_hm(self, multipoles): 

1210 import copy 

1211 required = [ 

1212 "mass_1", "mass_2", "spin_1z", "spin_2z", "psi", "iota", "ra", 

1213 "dec", "geocent_time", "luminosity_distance", "phase", 

1214 "reference_frequency" 

1215 ] 

1216 samples = self.specific_parameter_samples(required) 

1217 _f_low = self._retrieve_stored_frequency("f_low") 

1218 if isinstance(_f_low, (np.ndarray)): 

1219 f_low = _f_low() * len(samples[0]) 

1220 else: 

1221 f_low = [_f_low] * len(samples[0]) 

1222 original_list = copy.deepcopy(multipoles) 

1223 rho, data_used = multipole_snr( 

1224 *samples[:-1], f_low=f_low, psd=self.psd, f_ref=samples[-1], 

1225 f_final=self.extra_kwargs["meta_data"]["f_final"], 

1226 return_data_used=True, multi_process=self.multi_process, 

1227 psd_default=self.psd_default, multipole=multipoles, 

1228 df=self.extra_kwargs["meta_data"]["delta_f"] 

1229 ) 

1230 for num, mm in enumerate(original_list): 

1231 self.append_data("network_{}_multipole_snr".format(mm), rho[num]) 

1232 self.extra_kwargs["meta_data"]["multipole_snr"] = data_used 

1233 

1234 def _rho_p(self): 

1235 required = [ 

1236 "mass_1", "mass_2", "beta", "psi_J", "a_1", "a_2", "tilt_1", 

1237 "tilt_2", "phi_12", "theta_jn", "ra", "dec", "geocent_time", 

1238 "phi_jl", "reference_frequency", "luminosity_distance", "phase" 

1239 ] 

1240 samples = self.specific_parameter_samples(required) 

1241 try: 

1242 spins = self.specific_parameter_samples(["spin_1z", "spin_2z"]) 

1243 except ValueError: 

1244 spins = [None, None] 

1245 _f_low = self._retrieve_stored_frequency("f_low") 

1246 if isinstance(_f_low, (np.ndarray)): 

1247 f_low = _f_low() * len(samples[0]) 

1248 else: 

1249 f_low = [_f_low] * len(samples[0]) 

1250 [rho_p, b_bar, overlap, snrs], data_used = precessing_snr( 

1251 samples[0], samples[1], samples[2], samples[3], samples[4], 

1252 samples[5], samples[6], samples[7], samples[8], samples[9], samples[10], 

1253 samples[11], samples[12], samples[13], samples[15], samples[16], f_low=f_low, 

1254 spin_1z=spins[0], spin_2z=spins[1], psd=self.psd, return_data_used=True, 

1255 f_final=self.extra_kwargs["meta_data"]["f_final"], f_ref=samples[14], 

1256 multi_process=self.multi_process, psd_default=self.psd_default, 

1257 df=self.extra_kwargs["meta_data"]["delta_f"], debug=True 

1258 ) 

1259 self.append_data("network_precessing_snr", rho_p) 

1260 self.append_data("_b_bar", b_bar) 

1261 self.append_data("_precessing_harmonics_overlap", overlap) 

1262 nbreakdown = len(np.argwhere(b_bar > 0.3)) 

1263 if nbreakdown > 0: 

1264 logger.warning( 

1265 "{}/{} ({}%) samples have b_bar greater than 0.3. For these " 

1266 "samples, the two-harmonic approximation used to calculate " 

1267 "the precession SNR may not be valid".format( 

1268 nbreakdown, len(b_bar), 

1269 np.round((nbreakdown / len(b_bar)) * 100, 2) 

1270 ) 

1271 ) 

1272 try: 

1273 _samples = self.specific_parameter_samples("network_optimal_snr") 

1274 if np.logical_or( 

1275 np.median(snrs) > 1.1 * np.median(_samples), 

1276 np.median(snrs) < 0.9 * np.median(_samples) 

1277 ): 

1278 logger.warning( 

1279 "The two-harmonic SNR is different from the stored SNR. " 

1280 "This indicates that the provided PSD may be different " 

1281 "from the one used in the sampling." 

1282 ) 

1283 except Exception: 

1284 pass 

1285 self.extra_kwargs["meta_data"]["precessing_snr"] = data_used 

1286 

1287 def _retrieve_stored_frequency(self, name): 

1288 extra_kwargs = self.extra_kwargs["meta_data"] 

1289 if extra_kwargs != {} and name in list(extra_kwargs.keys()): 

1290 freq = extra_kwargs[name] 

1291 else: 

1292 raise ValueError( 

1293 "Could not find f_low in input file. Please either modify the " 

1294 "input file or pass it from the command line" 

1295 ) 

1296 return freq 

1297 

1298 def _retrieve_approximant(self): 

1299 extra_kwargs = self.extra_kwargs["meta_data"] 

1300 if extra_kwargs != {} and "approximant" in list(extra_kwargs.keys()): 

1301 approximant = extra_kwargs["approximant"] 

1302 else: 

1303 raise ValueError( 

1304 "Unable to find the approximant used to generate the posterior " 

1305 "samples in the result file." 

1306 ) 

1307 return approximant 

1308 

1309 def _evolve_spins(self, final_velocity="ISCO", forward=True): 

1310 from .evolve import evolve_spins 

1311 from ..waveform import _check_approximant_from_string 

1312 

1313 samples = self.specific_parameter_samples( 

1314 ["mass_1", "mass_2", "a_1", "a_2", "tilt_1", "tilt_2", 

1315 "phi_12", "reference_frequency"] 

1316 ) 

1317 approximant = self._retrieve_approximant() 

1318 if not _check_approximant_from_string(approximant): 

1319 _msg = ( 

1320 'Not evolving spins: approximant {0} unknown to ' 

1321 'lalsimulation and gwsignal'.format(approximant) 

1322 ) 

1323 logger.warning(_msg) 

1324 raise EvolveSpinError(_msg) 

1325 f_low = self._retrieve_stored_frequency("f_low") 

1326 if not forward: 

1327 [tilt_1_evolved, tilt_2_evolved, phi_12_evolved], fits_used = evolve_spins( 

1328 samples[0], samples[1], samples[2], samples[3], samples[4], 

1329 samples[5], samples[6], f_low, samples[7][0], 

1330 approximant, evolve_limit="infinite_separation", 

1331 multi_process=self.multi_process, return_fits_used=True, 

1332 method=self.evolve_spins_backwards 

1333 ) 

1334 suffix = "" 

1335 if self.evolve_spins_backwards.lower() == "precession_averaged": 

1336 suffix = "_only_prec_avg" 

1337 self.append_data("tilt_1_infinity{}".format(suffix), tilt_1_evolved) 

1338 self.append_data("tilt_2_infinity{}".format(suffix), tilt_2_evolved) 

1339 self.extra_kwargs["meta_data"]["backward_spin_evolution"] = fits_used 

1340 return 

1341 else: 

1342 tilt_1_evolved, tilt_2_evolved, phi_12_evolved = evolve_spins( 

1343 samples[0], samples[1], samples[2], samples[3], samples[4], 

1344 samples[5], samples[6], f_low, samples[7][0], 

1345 approximant, final_velocity=final_velocity, 

1346 multi_process=self.multi_process 

1347 ) 

1348 self.extra_kwargs["meta_data"]["forward_spin_evolution"] = final_velocity 

1349 spin_1z_evolved = samples[2] * np.cos(tilt_1_evolved) 

1350 spin_2z_evolved = samples[3] * np.cos(tilt_2_evolved) 

1351 self.append_data("tilt_1_evolved", tilt_1_evolved) 

1352 self.append_data("tilt_2_evolved", tilt_2_evolved) 

1353 self.append_data("phi_12_evolved", phi_12_evolved) 

1354 self.append_data("spin_1z_evolved", spin_1z_evolved) 

1355 self.append_data("spin_2z_evolved", spin_2z_evolved) 

1356 

1357 @staticmethod 

1358 def _evolved_vs_non_evolved_parameter( 

1359 parameter, evolved=False, core_param=False, non_precessing=False 

1360 ): 

1361 if non_precessing: 

1362 base_string = "" 

1363 elif evolved and core_param: 

1364 base_string = "_evolved" 

1365 elif evolved: 

1366 base_string = "" 

1367 elif core_param: 

1368 base_string = "" 

1369 else: 

1370 base_string = "_non_evolved" 

1371 return "{}{}".format(parameter, base_string) 

1372 

1373 def _precessing_vs_non_precessing_parameters( 

1374 self, non_precessing=False, evolved=False 

1375 ): 

1376 if not non_precessing: 

1377 tilt_1 = self._evolved_vs_non_evolved_parameter( 

1378 "tilt_1", evolved=evolved, core_param=True 

1379 ) 

1380 tilt_2 = self._evolved_vs_non_evolved_parameter( 

1381 "tilt_2", evolved=evolved, core_param=True 

1382 ) 

1383 samples = self.specific_parameter_samples([ 

1384 "mass_1", "mass_2", "a_1", "a_2", tilt_1, tilt_2 

1385 ]) 

1386 if "phi_12" in self.parameters and evolved: 

1387 phi_12_samples = self.specific_parameter_samples([ 

1388 self._evolved_vs_non_evolved_parameter( 

1389 "phi_12", evolved=True, core_param=True 

1390 ) 

1391 ])[0] 

1392 elif "phi_12" in self.parameters: 

1393 phi_12_samples = self.specific_parameter_samples(["phi_12"])[0] 

1394 else: 

1395 phi_12_samples = np.zeros_like(samples[0]) 

1396 samples.append(phi_12_samples) 

1397 if self.NRSurrogate: 

1398 NRSurrogate_samples = self.specific_parameter_samples([ 

1399 "phi_jl", "theta_jn", "phase" 

1400 ]) 

1401 for ss in NRSurrogate_samples: 

1402 samples.append(ss) 

1403 else: 

1404 spin_1z = self._evolved_vs_non_evolved_parameter( 

1405 "spin_1z", evolved=evolved, core_param=True, non_precessing=True 

1406 ) 

1407 spin_2z = self._evolved_vs_non_evolved_parameter( 

1408 "spin_2z", evolved=evolved, core_param=True, non_precessing=True 

1409 ) 

1410 samples = self.specific_parameter_samples([ 

1411 "mass_1", "mass_2", spin_1z, spin_2z 

1412 ]) 

1413 samples = [ 

1414 samples[0], samples[1], np.abs(samples[2]), np.abs(samples[3]), 

1415 0.5 * np.pi * (1 - np.sign(samples[2])), 

1416 0.5 * np.pi * (1 - np.sign(samples[3])), 

1417 np.zeros_like(samples[0]) 

1418 ] 

1419 return samples 

1420 

1421 def _peak_luminosity_of_merger(self, evolved=False): 

1422 param = self._evolved_vs_non_evolved_parameter( 

1423 "peak_luminosity", evolved=evolved, non_precessing=self.non_precessing 

1424 ) 

1425 spin_1z_param = self._evolved_vs_non_evolved_parameter( 

1426 "spin_1z", evolved=evolved, core_param=True, 

1427 non_precessing=self.non_precessing 

1428 ) 

1429 spin_2z_param = self._evolved_vs_non_evolved_parameter( 

1430 "spin_2z", evolved=evolved, core_param=True, 

1431 non_precessing=self.non_precessing 

1432 ) 

1433 

1434 samples = self.specific_parameter_samples([ 

1435 "mass_1", "mass_2", spin_1z_param, spin_2z_param 

1436 ]) 

1437 peak_luminosity, fits = peak_luminosity_of_merger( 

1438 samples[0], samples[1], samples[2], samples[3], 

1439 return_fits_used=True 

1440 ) 

1441 self.append_data(param, peak_luminosity) 

1442 self.extra_kwargs["meta_data"]["peak_luminosity_NR_fits"] = fits 

1443 

1444 def _final_remnant_properties_from_NRSurrogate( 

1445 self, non_precessing=False, 

1446 parameters=["final_mass", "final_spin", "final_kick"] 

1447 ): 

1448 f_low = self._retrieve_stored_frequency("f_start") 

1449 approximant = self._retrieve_approximant() 

1450 samples = self._precessing_vs_non_precessing_parameters( 

1451 non_precessing=non_precessing, evolved=False 

1452 ) 

1453 frequency_samples = self.specific_parameter_samples([ 

1454 "reference_frequency" 

1455 ]) 

1456 data, fits = final_remnant_properties_from_NRSurrogate( 

1457 *samples, f_low=f_low, f_ref=frequency_samples[0], 

1458 properties=parameters, return_fits_used=True, 

1459 approximant=approximant 

1460 ) 

1461 for param in parameters: 

1462 self.append_data(param, data[param]) 

1463 self.extra_kwargs["meta_data"]["{}_NR_fits".format(param)] = fits 

1464 

1465 def _final_remnant_properties_from_NSBH_waveform( 

1466 self, source=False, parameters=[ 

1467 "baryonic_torus_mass", "final_mass", "final_spin" 

1468 ] 

1469 ): 

1470 approximant = self._retrieve_approximant() 

1471 if source: 

1472 sample_params = [ 

1473 "mass_1_source", "mass_2_source", "spin_1z", "lambda_2" 

1474 ] 

1475 else: 

1476 sample_params = ["mass_1", "mass_2", "spin_1z", "lambda_2"] 

1477 samples = self.specific_parameter_samples(sample_params) 

1478 _data = _check_NSBH_approximant( 

1479 approximant, samples[0], samples[1], samples[2], samples[3], 

1480 _raise=False 

1481 ) 

1482 if _data is None: 

1483 return 

1484 _mapping = { 

1485 "220_quasinormal_mode_frequency": 0, "tidal_disruption_frequency": 1, 

1486 "baryonic_torus_mass": 2, "compactness_2": 3, 

1487 "final_mass": 4, "final_spin": 5 

1488 } 

1489 for param in parameters: 

1490 self.append_data(param, _data[_mapping[param]]) 

1491 if "final_mass" in parameters: 

1492 self.extra_kwargs["meta_data"]["final_mass_NR_fits"] = "Zappa2019" 

1493 if "final_spin" in parameters: 

1494 self.extra_kwargs["meta_data"]["final_spin_NR_fits"] = "Zappa2019" 

1495 if "baryonic_torus_mass" in parameters: 

1496 self.extra_kwargs["meta_data"]["baryonic_torus_mass_fits"] = ( 

1497 "Foucart2012" 

1498 ) 

1499 if "220_quasinormal_mode_frequency" in parameters: 

1500 self.extra_kwargs["meta_data"]["quasinormal_mode_fits"] = ( 

1501 "London2019" 

1502 ) 

1503 if "tidal_disruption_frequency" in parameters: 

1504 probabilities = NSBH_merger_type( 

1505 samples[0], samples[1], samples[2], samples[3], 

1506 approximant=approximant, 

1507 _ringdown=_data[_mapping["220_quasinormal_mode_frequency"]], 

1508 _disruption=_data[_mapping["tidal_disruption_frequency"]], 

1509 _torus=_data[_mapping["baryonic_torus_mass"]], percentages=True 

1510 ) 

1511 self.extra_kwargs["meta_data"]["NSBH_merger_type_probabilities"] = ( 

1512 probabilities 

1513 ) 

1514 self.extra_kwargs["meta_data"]["tidal_disruption_frequency_fits"] = ( 

1515 "Pannarale2018" 

1516 ) 

1517 ratio = ( 

1518 _data[_mapping["tidal_disruption_frequency"]] 

1519 / _data[_mapping["220_quasinormal_mode_frequency"]] 

1520 ) 

1521 self.append_data( 

1522 "tidal_disruption_frequency_ratio", ratio 

1523 ) 

1524 

1525 def _final_remnant_properties_from_waveform( 

1526 self, non_precessing=False, parameters=["final_mass", "final_spin"], 

1527 ): 

1528 f_low = self._retrieve_stored_frequency("f_start") 

1529 approximant = self._retrieve_approximant() 

1530 if "delta_t" in self.extra_kwargs["meta_data"].keys(): 

1531 delta_t = self.extra_kwargs["meta_data"]["delta_t"] 

1532 else: 

1533 delta_t = 1. / 4096 

1534 if "seob" in approximant.lower(): 

1535 logger.warning( 

1536 "Could not find 'delta_t' in the meta data. Using {} as " 

1537 "default.".format(delta_t) 

1538 ) 

1539 if non_precessing: 

1540 sample_params = [ 

1541 "mass_1", "mass_2", "empty", "empty", "spin_1z", "empty", 

1542 "empty", "spin_2z", "iota", "luminosity_distance", 

1543 "phase" 

1544 ] 

1545 else: 

1546 sample_params = [ 

1547 "mass_1", "mass_2", "spin_1x", "spin_1y", "spin_1z", 

1548 "spin_2x", "spin_2y", "spin_2z", "iota", "luminosity_distance", 

1549 "phase", "reference_frequency" 

1550 ] 

1551 samples = self.specific_parameter_samples(sample_params) 

1552 ind = self.parameters.index("spin_1x") 

1553 _data, fits = _final_from_initial_BBH( 

1554 *samples[:8], iota=samples[8], luminosity_distance=samples[9], 

1555 f_low=[f_low] * len(samples[0]), f_ref=samples[-1], 

1556 phi_ref=samples[10], delta_t=1. / 4096, approximant=approximant, 

1557 xphm_flags=self.extra_kwargs["meta_data"].get("approximant_flags", {}), 

1558 return_fits_used=True, multi_process=self.multi_process 

1559 ) 

1560 data = {"final_mass": _data[0], "final_spin": _data[1]} 

1561 for param in parameters: 

1562 self.append_data(param, data[param]) 

1563 self.extra_kwargs["meta_data"]["{}_NR_fits".format(param)] = fits 

1564 

1565 def _final_mass_of_merger(self, evolved=False): 

1566 param = self._evolved_vs_non_evolved_parameter( 

1567 "final_mass", evolved=evolved, non_precessing=self.non_precessing 

1568 ) 

1569 spin_1z_param = self._evolved_vs_non_evolved_parameter( 

1570 "spin_1z", evolved=evolved, core_param=True, 

1571 non_precessing=self.non_precessing 

1572 ) 

1573 spin_2z_param = self._evolved_vs_non_evolved_parameter( 

1574 "spin_2z", evolved=evolved, core_param=True, 

1575 non_precessing=self.non_precessing 

1576 ) 

1577 samples = self.specific_parameter_samples([ 

1578 "mass_1", "mass_2", spin_1z_param, spin_2z_param 

1579 ]) 

1580 final_mass, fits = final_mass_of_merger( 

1581 *samples, return_fits_used=True, 

1582 ) 

1583 self.append_data(param, final_mass) 

1584 self.extra_kwargs["meta_data"]["final_mass_NR_fits"] = fits 

1585 

1586 def _final_mass_source(self, evolved=False): 

1587 param = self._evolved_vs_non_evolved_parameter( 

1588 "final_mass", evolved=evolved, non_precessing=self.non_precessing 

1589 ) 

1590 samples = self.specific_parameter_samples([param, "redshift"]) 

1591 final_mass_source = _source_from_detector( 

1592 samples[0], samples[1] 

1593 ) 

1594 self.append_data(param.replace("mass", "mass_source"), final_mass_source) 

1595 

1596 def _final_spin_of_merger(self, non_precessing=False, evolved=False): 

1597 param = self._evolved_vs_non_evolved_parameter( 

1598 "final_spin", evolved=evolved, non_precessing=self.non_precessing 

1599 ) 

1600 samples = self._precessing_vs_non_precessing_parameters( 

1601 non_precessing=non_precessing, evolved=evolved 

1602 ) 

1603 final_spin, fits = final_spin_of_merger( 

1604 *samples, return_fits_used=True, 

1605 ) 

1606 self.append_data(param, final_spin) 

1607 self.extra_kwargs["meta_data"]["final_spin_NR_fits"] = fits 

1608 

1609 def _radiated_energy(self, evolved=False): 

1610 param = self._evolved_vs_non_evolved_parameter( 

1611 "radiated_energy", evolved=evolved, non_precessing=self.non_precessing 

1612 ) 

1613 final_mass_param = self._evolved_vs_non_evolved_parameter( 

1614 "final_mass_source", evolved=evolved, non_precessing=self.non_precessing 

1615 ) 

1616 samples = self.specific_parameter_samples([ 

1617 "total_mass_source", final_mass_param 

1618 ]) 

1619 radiated_energy = samples[0] - samples[1] 

1620 self.append_data(param, radiated_energy) 

1621 

1622 def _cos_angle(self, parameter_to_add, reverse=False): 

1623 if reverse: 

1624 samples = self.specific_parameter_samples( 

1625 ["cos_" + parameter_to_add]) 

1626 cos_samples = np.arccos(samples[0]) 

1627 else: 

1628 samples = self.specific_parameter_samples( 

1629 [parameter_to_add.split("cos_")[1]] 

1630 ) 

1631 cos_samples = np.cos(samples[0]) 

1632 self.append_data(parameter_to_add, cos_samples) 

1633 

1634 def source_frame_from_detector_frame(self, detector_frame_parameter): 

1635 samples = self.specific_parameter_samples( 

1636 [detector_frame_parameter, "redshift"] 

1637 ) 

1638 source_frame = _source_from_detector(samples[0], samples[1]) 

1639 self.append_data( 

1640 "{}_source".format(detector_frame_parameter), source_frame 

1641 ) 

1642 

1643 def _check_parameters(self): 

1644 params = ["mass_1", "mass_2", "a_1", "a_2", "mass_1_source", "mass_2_source", 

1645 "mass_ratio", "total_mass", "chirp_mass"] 

1646 for i in params: 

1647 if i in self.parameters: 

1648 samples = self.specific_parameter_samples([i]) 

1649 if "mass" in i: 

1650 cond = any(np.array(samples[0]) <= 0.) 

1651 else: 

1652 cond = any(np.array(samples[0]) < 0.) 

1653 if cond: 

1654 if "mass" in i: 

1655 ind = np.argwhere(np.array(samples[0]) <= 0.) 

1656 else: 

1657 ind = np.argwhere(np.array(samples[0]) < 0.) 

1658 logger.warning( 

1659 "Removing %s samples because they have unphysical " 

1660 "values (%s < 0)" % (len(ind), i) 

1661 ) 

1662 for i in np.arange(len(ind) - 1, -1, -1): 

1663 self.samples.remove(list(np.array(self.samples)[ind[i][0]])) 

1664 

1665 def generate_all_posterior_samples(self): 

1666 logger.debug("Starting to generate all derived posteriors") 

1667 evolve_condition = ( 

1668 True if self.evolve_spins_forwards and self.compute_remnant else False 

1669 ) 

1670 if "cos_theta_jn" in self.parameters and "theta_jn" not in self.parameters: 

1671 self._cos_angle("theta_jn", reverse=True) 

1672 if "cos_iota" in self.parameters and "iota" not in self.parameters: 

1673 self._cos_angle("iota", reverse=True) 

1674 if "cos_tilt_1" in self.parameters and "tilt_1" not in self.parameters: 

1675 self._cos_angle("tilt_1", reverse=True) 

1676 if "cos_tilt_2" in self.parameters and "tilt_2" not in self.parameters: 

1677 self._cos_angle("tilt_2", reverse=True) 

1678 angles = [ 

1679 "a_1", "a_2", "a_1_azimuthal", "a_1_polar", "a_2_azimuthal", 

1680 "a_2_polar" 

1681 ] 

1682 if all(i in self.parameters for i in angles): 

1683 self._component_spins_from_azimuthal_and_polar_angles() 

1684 spin_magnitudes = ["a_1", "a_2"] 

1685 angles = ["phi_jl", "tilt_1", "tilt_2", "phi_12"] 

1686 cartesian = ["spin_1x", "spin_1y", "spin_1z", "spin_2x", "spin_2y", "spin_2z"] 

1687 cond1 = all(i in self.parameters for i in spin_magnitudes) 

1688 cond2 = all(i in self.parameters for i in angles) 

1689 cond3 = all(i in self.parameters for i in cartesian) 

1690 for _param in spin_magnitudes: 

1691 if _param in self.parameters and not cond2 and not cond3: 

1692 _index = _param.split("a_")[1] 

1693 _spin = self.specific_parameter_samples(_param) 

1694 _tilt = np.arccos(np.sign(_spin)) 

1695 self.append_data("tilt_{}".format(_index), _tilt) 

1696 _spin_ind = self.parameters.index(_param) 

1697 for num, i in enumerate(self.samples): 

1698 self.samples[num][_spin_ind] = abs(self.samples[num][_spin_ind]) 

1699 

1700 if not cond2 and not cond3 and self.add_zero_spin: 

1701 for _param in spin_magnitudes: 

1702 if _param not in self.parameters: 

1703 _spin = np.zeros(len(self.samples)) 

1704 self.append_data(_param, _spin) 

1705 _index = _param.split("a_")[1] 

1706 self.append_data("spin_{}z".format(_index), _spin) 

1707 self._check_parameters() 

1708 if "cos_theta_jn" in self.parameters and "theta_jn" not in self.parameters: 

1709 self._cos_angle("theta_jn", reverse=True) 

1710 if "cos_iota" in self.parameters and "iota" not in self.parameters: 

1711 self._cos_angle("iota", reverse=True) 

1712 if "cos_tilt_1" in self.parameters and "tilt_1" not in self.parameters: 

1713 self._cos_angle("tilt_1", reverse=True) 

1714 if "cos_tilt_2" in self.parameters and "tilt_2" not in self.parameters: 

1715 self._cos_angle("tilt_2", reverse=True) 

1716 if "luminosity_distance" not in self.parameters: 

1717 if "redshift" in self.parameters: 

1718 self._dL_from_z() 

1719 if "redshift" not in self.parameters: 

1720 if "luminosity_distance" in self.parameters: 

1721 self._z_from_dL() 

1722 if "comoving_volume" in self.parameters: 

1723 self._z_from_comoving_volume() 

1724 if "comoving_distance" not in self.parameters: 

1725 if "redshift" in self.parameters: 

1726 self._comoving_distance_from_z() 

1727 if "comoving_volume" not in self.parameters: 

1728 if "redshift" in self.parameters: 

1729 self._comoving_volume_from_z() 

1730 

1731 if "mass_ratio" not in self.parameters and "symmetric_mass_ratio" in \ 

1732 self.parameters: 

1733 self._q_from_eta() 

1734 if "mass_ratio" not in self.parameters and "mass_1" in self.parameters \ 

1735 and "mass_2" in self.parameters: 

1736 self._q_from_m1_m2() 

1737 if "mass_ratio" in self.parameters: 

1738 ind = self.parameters.index("mass_ratio") 

1739 median = np.median([i[ind] for i in self.samples]) 

1740 if median > 1.: 

1741 self._invert_q() 

1742 if "inverted_mass_ratio" not in self.parameters and "mass_ratio" in \ 

1743 self.parameters: 

1744 self._invq_from_q() 

1745 if "chirp_mass" not in self.parameters and "total_mass" in self.parameters: 

1746 self._mchirp_from_mtotal_q() 

1747 if "mass_1" not in self.parameters and "chirp_mass" in self.parameters: 

1748 self._m1_from_mchirp_q() 

1749 if "mass_2" not in self.parameters and "chirp_mass" in self.parameters: 

1750 self._m2_from_mchirp_q() 

1751 if "mass_1" not in self.parameters and "total_mass" in self.parameters: 

1752 self._m1_from_mtotal_q() 

1753 if "mass_2" not in self.parameters and "total_mass" in self.parameters: 

1754 self._m2_from_mtotal_q() 

1755 if "mass_1" in self.parameters and "mass_2" in self.parameters: 

1756 if "total_mass" not in self.parameters: 

1757 self._mtotal_from_m1_m2() 

1758 if "chirp_mass" not in self.parameters: 

1759 self._mchirp_from_m1_m2() 

1760 if "symmetric_mass_ratio" not in self.parameters: 

1761 self._eta_from_m1_m2() 

1762 if "redshift" in self.parameters: 

1763 if "mass_1_source" not in self.parameters: 

1764 if "mass_1" in self.parameters: 

1765 self._m1_source_from_m1_z() 

1766 if "mass_1_source" in self.parameters: 

1767 if "mass_1" not in self.parameters: 

1768 self._m1_from_m1_source_z() 

1769 if "mass_2_source" not in self.parameters: 

1770 if "mass_2" in self.parameters: 

1771 self._m2_source_from_m2_z() 

1772 if "mass_2_source" in self.parameters: 

1773 if "mass_2" not in self.parameters: 

1774 self._m2_from_m2_source_z() 

1775 if "total_mass_source" not in self.parameters: 

1776 if "total_mass" in self.parameters: 

1777 self._mtotal_source_from_mtotal_z() 

1778 if "total_mass_source" in self.parameters: 

1779 if "total_mass" not in self.parameters: 

1780 self._mtotal_from_mtotal_source_z() 

1781 if "chirp_mass_source" not in self.parameters: 

1782 if "chirp_mass" in self.parameters: 

1783 self._mchirp_source_from_mchirp_z() 

1784 if "chirp_mass_source" in self.parameters: 

1785 if "chirp_mass" not in self.parameters: 

1786 self._mchirp_from_mchirp_source_z() 

1787 

1788 if "reference_frequency" not in self.parameters: 

1789 self._reference_frequency() 

1790 condition1 = "phi_12" not in self.parameters 

1791 condition2 = "phi_1" in self.parameters and "phi_2" in self.parameters 

1792 if condition1 and condition2: 

1793 self._phi_12_from_phi1_phi2() 

1794 

1795 check_for_evolved_parameter = lambda suffix, param, params: ( 

1796 param not in params and param + suffix not in params if 

1797 len(suffix) else param not in params 

1798 ) 

1799 

1800 if "mass_1" in self.parameters and "mass_2" in self.parameters: 

1801 spin_components = [ 

1802 "spin_1x", "spin_1y", "spin_1z", "spin_2x", "spin_2y", "spin_2z", 

1803 "iota" 

1804 ] 

1805 angles = ["a_1", "a_2", "tilt_1", "tilt_2", "theta_jn"] 

1806 if all(i in self.parameters for i in spin_components): 

1807 self._spin_angles() 

1808 if all(i in self.parameters for i in angles): 

1809 samples = self.specific_parameter_samples(["tilt_1", "tilt_2"]) 

1810 cond1 = all(i in [0, np.pi] for i in samples[0]) 

1811 cond2 = all(i in [0, np.pi] for i in samples[1]) 

1812 if cond1 and cond1: 

1813 self._non_precessing_component_spins() 

1814 else: 

1815 angles = [ 

1816 "phi_jl", "phi_12", "reference_frequency"] 

1817 if all(i in self.parameters for i in angles): 

1818 self._component_spins() 

1819 cond1 = "spin_1x" in self.parameters and "spin_1y" in self.parameters 

1820 if "phi_1" not in self.parameters and cond1: 

1821 self._phi1_from_spins() 

1822 cond1 = "spin_2x" in self.parameters and "spin_2y" in self.parameters 

1823 if "phi_2" not in self.parameters and cond1: 

1824 self._phi2_from_spins() 

1825 evolve_spins_params = ["tilt_1", "tilt_2", "phi_12"] 

1826 if self.evolve_spins_backwards: 

1827 if all(i in self.parameters for i in evolve_spins_params): 

1828 try: 

1829 self._evolve_spins(forward=False) 

1830 except EvolveSpinError: 

1831 # Raised when approximant is unknown to lalsimulation or 

1832 # lalsimulation.SimInspiralGetSpinFreqFromApproximant is 

1833 # equal to lalsimulation.SIM_INSPIRAL_SPINS_CASEBYCASE 

1834 pass 

1835 for suffix in ["_infinity", "_infinity_only_prec_avg", ""]: 

1836 if "spin_1z{}".format(suffix) not in self.parameters: 

1837 _params = ["a_1", "tilt_1{}".format(suffix)] 

1838 if all(i in self.parameters for i in _params): 

1839 self._aligned_spin_from_magnitude_tilts( 

1840 primary=True, suffix=suffix 

1841 ) 

1842 if "spin_2z{}".format(suffix) not in self.parameters: 

1843 _params = ["a_2", "tilt_2{}".format(suffix)] 

1844 if all(i in self.parameters for i in _params): 

1845 self._aligned_spin_from_magnitude_tilts( 

1846 secondary=True, suffix=suffix 

1847 ) 

1848 if "chi_eff{}".format(suffix) not in self.parameters: 

1849 _params = ["spin_1z{}".format(suffix), "spin_2z{}".format(suffix)] 

1850 if all(i in self.parameters for i in _params): 

1851 self._chi_eff(suffix=suffix) 

1852 if any( 

1853 _p.format(suffix) not in self.parameters for _p in 

1854 ["chi_p{}", "chi_p_2spin"] 

1855 ): 

1856 _params = [ 

1857 "a_1", "tilt_1{}".format(suffix), "a_2", 

1858 "tilt_2{}".format(suffix) 

1859 ] 

1860 _cartesian_params = ["spin_1x", "spin_1y", "spin_2x", "spin_2y"] 

1861 if "chi_p{}".format(suffix) not in self.parameters: 

1862 if all(i in self.parameters for i in _params): 

1863 self._chi_p_from_tilts(suffix=suffix) 

1864 elif all(i in self.parameters for i in _cartesian_params): 

1865 self._chi_p() 

1866 if "chi_p_2spin" not in self.parameters: 

1867 if all(i in self.parameters for i in _cartesian_params): 

1868 self._chi_p_2spin() 

1869 if "beta" not in self.parameters: 

1870 beta_components = [ 

1871 "mass_1", "mass_2", "phi_jl", "tilt_1", "tilt_2", "phi_12", 

1872 "a_1", "a_2", "reference_frequency", "phase" 

1873 ] 

1874 if all(i in self.parameters for i in beta_components): 

1875 self._beta() 

1876 polytrope_params = ["log_pressure", "gamma_1", "gamma_2", "gamma_3"] 

1877 if all(param in self.parameters for param in polytrope_params): 

1878 if "lambda_1" not in self.parameters or "lambda_2" not in self.parameters: 

1879 self._lambda1_lambda2_from_polytrope_EOS() 

1880 spectral_params = [ 

1881 "spectral_decomposition_gamma_{}".format(num) for num in 

1882 np.arange(4) 

1883 ] 

1884 if all(param in self.parameters for param in spectral_params): 

1885 if "lambda_1" not in self.parameters or "lambda_2" not in self.parameters: 

1886 self._lambda1_lambda2_from_spectral_decomposition_EOS() 

1887 if "lambda_tilde" in self.parameters and "lambda_1" not in self.parameters: 

1888 self._lambda1_from_lambda_tilde() 

1889 if "lambda_2" not in self.parameters and "lambda_1" in self.parameters: 

1890 self._lambda2_from_lambda1() 

1891 if "lambda_1" in self.parameters and "lambda_2" in self.parameters: 

1892 if "lambda_tilde" not in self.parameters: 

1893 self._lambda_tilde_from_lambda1_lambda2() 

1894 if "delta_lambda" not in self.parameters: 

1895 self._delta_lambda_from_lambda1_lambda2() 

1896 if "psi" in self.parameters: 

1897 dpsi_parameters = ["theta_jn", "phi_jl", "beta"] 

1898 if all(i in self.parameters for i in dpsi_parameters): 

1899 if "psi_J" not in self.parameters: 

1900 self._psi_J() 

1901 

1902 evolve_suffix = "_non_evolved" 

1903 final_spin_params = ["a_1", "a_2"] 

1904 non_precessing_NR_params = ["spin_1z", "spin_2z"] 

1905 if evolve_condition: 

1906 final_spin_params += [ 

1907 "tilt_1_evolved", "tilt_2_evolved", "phi_12_evolved" 

1908 ] 

1909 non_precessing_NR_params = [ 

1910 "{}_evolved".format(i) for i in non_precessing_NR_params 

1911 ] 

1912 evolve_suffix = "_evolved" 

1913 if all(i in self.parameters for i in evolve_spins_params): 

1914 try: 

1915 self._evolve_spins(final_velocity=self.evolve_spins_forwards) 

1916 except EvolveSpinError: 

1917 # Raised when approximant is unknown to lalsimulation or 

1918 # lalsimulation.SimInspiralGetSpinFreqFromApproximant is 

1919 # equal to lalsimulation.SIM_INSPIRAL_SPINS_CASEBYCASE 

1920 evolve_condition = False 

1921 else: 

1922 evolve_condition = False 

1923 else: 

1924 final_spin_params += ["tilt_1", "tilt_2", "phi_12"] 

1925 

1926 condition_peak_luminosity = check_for_evolved_parameter( 

1927 evolve_suffix, "peak_luminosity", self.parameters 

1928 ) 

1929 condition_final_spin = check_for_evolved_parameter( 

1930 evolve_suffix, "final_spin", self.parameters 

1931 ) 

1932 condition_final_mass = check_for_evolved_parameter( 

1933 evolve_suffix, "final_mass", self.parameters 

1934 ) 

1935 if (self.NRSurrogate or self.waveform_fit) and self.compute_remnant: 

1936 parameters = [] 

1937 _default = ["final_mass", "final_spin"] 

1938 if self.NRSurrogate: 

1939 _default.append("final_kick") 

1940 function = self._final_remnant_properties_from_NRSurrogate 

1941 else: 

1942 final_spin_params = [ 

1943 "spin_1x", "spin_1y", "spin_1z", "spin_2x", 

1944 "spin_2y", "spin_2z" 

1945 ] 

1946 function = self._final_remnant_properties_from_waveform 

1947 

1948 for param in _default: 

1949 if param not in self.parameters: 

1950 parameters.append(param) 

1951 # We already know that lambda_2 is in the posterior table if 

1952 # self.NSBH = True 

1953 if self.NSBH and "spin_1z" in self.parameters: 

1954 self._final_remnant_properties_from_NSBH_waveform() 

1955 elif all(i in self.parameters for i in final_spin_params): 

1956 function(non_precessing=False, parameters=parameters) 

1957 elif all(i in self.parameters for i in non_precessing_NR_params): 

1958 function(non_precessing=True, parameters=parameters) 

1959 if all(i in self.parameters for i in non_precessing_NR_params): 

1960 if condition_peak_luminosity or self.force_non_evolved: 

1961 if not self.NSBH: 

1962 self._peak_luminosity_of_merger(evolved=evolve_condition) 

1963 elif self.compute_remnant: 

1964 if all(i in self.parameters for i in final_spin_params): 

1965 if condition_final_spin or self.force_non_evolved: 

1966 self._final_spin_of_merger(evolved=evolve_condition) 

1967 elif all(i in self.parameters for i in non_precessing_NR_params): 

1968 if condition_final_spin or self.force_non_evolved: 

1969 self._final_spin_of_merger( 

1970 non_precessing=True, evolved=False 

1971 ) 

1972 if all(i in self.parameters for i in non_precessing_NR_params): 

1973 if condition_peak_luminosity or self.force_non_evolved: 

1974 self._peak_luminosity_of_merger(evolved=evolve_condition) 

1975 if condition_final_mass or self.force_non_evolved: 

1976 self._final_mass_of_merger(evolved=evolve_condition) 

1977 

1978 # if NSBH system and self.compute_remnant = False and/or BBH fits 

1979 # fits used, only calculate baryonic_torus_mass 

1980 if self.NSBH and "spin_1z" in self.parameters: 

1981 if "baryonic_torus_mass" not in self.parameters: 

1982 self._final_remnant_properties_from_NSBH_waveform( 

1983 parameters=["baryonic_torus_mass"] 

1984 ) 

1985 # calculate compactness from Love-compactness relation 

1986 if "lambda_1" in self.parameters and "compactness_1" not in self.parameters: 

1987 self._NS_compactness_from_lambda(parameter="lambda_1") 

1988 if "mass_1" in self.parameters and "baryonic_mass_1" not in self.parameters: 

1989 self._NS_baryonic_mass(primary=True) 

1990 if "lambda_2" in self.parameters and "compactness_2" not in self.parameters: 

1991 self._NS_compactness_from_lambda(parameter="lambda_2") 

1992 if "mass_2" in self.parameters and "baryonic_mass_2" not in self.parameters: 

1993 self._NS_baryonic_mass(primary=False) 

1994 for suffix in ["_infinity", "_infinity_only_prec_avg", ""]: 

1995 for tilt in ["tilt_1", "tilt_2"]: 

1996 cond1 = "cos_{}{}".format(tilt, suffix) not in self.parameters 

1997 cond2 = "{}{}".format(tilt, suffix) in self.parameters 

1998 if cond1 and cond2: 

1999 self._cos_angle("cos_{}{}".format(tilt, suffix)) 

2000 evolve_suffix = "_non_evolved" 

2001 if evolve_condition or self.NRSurrogate or self.waveform_fit or self.non_precessing: 

2002 evolve_suffix = "" 

2003 evolve_condition = True 

2004 if "redshift" in self.parameters: 

2005 condition_final_mass_source = check_for_evolved_parameter( 

2006 evolve_suffix, "final_mass_source", self.parameters 

2007 ) 

2008 if condition_final_mass_source or self.force_non_evolved: 

2009 if "final_mass{}".format(evolve_suffix) in self.parameters: 

2010 self._final_mass_source(evolved=evolve_condition) 

2011 if "baryonic_torus_mass" in self.parameters: 

2012 if "baryonic_torus_mass_source" not in self.parameters: 

2013 self.source_frame_from_detector_frame( 

2014 "baryonic_torus_mass" 

2015 ) 

2016 if "baryonic_mass_1" in self.parameters: 

2017 if "baryonic_mass_1_source" not in self.parameters: 

2018 self.source_frame_from_detector_frame( 

2019 "baryonic_mass_1" 

2020 ) 

2021 if "baryonic_mass_2" in self.parameters: 

2022 if "baryonic_mass_2_source" not in self.parameters: 

2023 self.source_frame_from_detector_frame( 

2024 "baryonic_mass_2" 

2025 ) 

2026 if "total_mass_source" in self.parameters: 

2027 if "final_mass_source{}".format(evolve_suffix) in self.parameters: 

2028 condition_radiated_energy = check_for_evolved_parameter( 

2029 evolve_suffix, "radiated_energy", self.parameters 

2030 ) 

2031 if condition_radiated_energy or self.force_non_evolved: 

2032 self._radiated_energy(evolved=evolve_condition) 

2033 if self.NSBH and "spin_1z" in self.parameters: 

2034 if all(_p in self.parameters for _p in ["mass_1_source", "mass_2_source"]): 

2035 _NSBH_parameters = [] 

2036 if "tidal_disruption_frequency" not in self.parameters: 

2037 _NSBH_parameters.append("tidal_disruption_frequency") 

2038 if "220_quasinormal_mode_frequency" not in self.parameters: 

2039 _NSBH_parameters.append("220_quasinormal_mode_frequency") 

2040 if len(_NSBH_parameters): 

2041 self._final_remnant_properties_from_NSBH_waveform( 

2042 parameters=_NSBH_parameters, source=True 

2043 ) 

2044 location = ["geocent_time", "ra", "dec"] 

2045 if all(i in self.parameters for i in location): 

2046 try: 

2047 self._time_in_each_ifo() 

2048 self._time_delay_between_ifos() 

2049 except Exception as e: 

2050 logger.warning( 

2051 "Failed to generate posterior samples for the time in each " 

2052 "detector because %s" % (e) 

2053 ) 

2054 if any("_matched_filter_snr_angle" in i for i in self.parameters): 

2055 if any("_matched_filter_abs_snr" in i for i in self.parameters): 

2056 self._ifo_snr() 

2057 if any("_optimal_snr" in i for i in self.parameters): 

2058 if "network_optimal_snr" not in self.parameters: 

2059 self._optimal_network_snr() 

2060 if any("_matched_filter_snr" in i for i in self.parameters): 

2061 if "network_matched_filter_snr" not in self.parameters: 

2062 self._matched_filter_network_snr() 

2063 if self.multipole_snr: 

2064 rho_hm_parameters = [ 

2065 "mass_1", "mass_2", "spin_1z", "spin_2z", "psi", "iota", "ra", 

2066 "dec", "geocent_time", "luminosity_distance", "phase" 

2067 ] 

2068 cond = [ 

2069 int(mm) for mm in ['21', '33', '44'] if 

2070 "network_{}_multipole_snr".format(mm) not in self.parameters 

2071 ] 

2072 if all(i in self.parameters for i in rho_hm_parameters) and len(cond): 

2073 try: 

2074 logger.warning( 

2075 "Starting to calculate the SNR in the {} multipole{}. " 

2076 "This may take some time".format( 

2077 " and ".join([str(mm) for mm in cond]), 

2078 "s" if len(cond) > 1 else "" 

2079 ) 

2080 ) 

2081 self._rho_hm(cond) 

2082 except ImportError as e: 

2083 logger.warning(e) 

2084 elif len(cond): 

2085 logger.warning( 

2086 "Unable to calculate the multipole SNR because it requires " 

2087 "samples for {}".format( 

2088 ", ".join( 

2089 [i for i in rho_hm_parameters if i not in self.parameters] 

2090 ) 

2091 ) 

2092 ) 

2093 if "network_precessing_snr" not in self.parameters and self.precessing_snr: 

2094 rho_p_parameters = [ 

2095 "mass_1", "mass_2", "beta", "psi_J", "a_1", "a_2", "tilt_1", 

2096 "tilt_2", "phi_12", "theta_jn", "phi_jl", "ra", "dec", "geocent_time", 

2097 "phi_jl" 

2098 ] 

2099 if all(i in self.parameters for i in rho_p_parameters): 

2100 try: 

2101 logger.warning( 

2102 "Starting to calculate the precessing SNR. This may take " 

2103 "some time" 

2104 ) 

2105 self._rho_p() 

2106 except ImportError as e: 

2107 logger.warning(e) 

2108 else: 

2109 logger.warning( 

2110 "Unable to calculate the precessing SNR because requires " 

2111 "samples for {}".format( 

2112 ", ".join( 

2113 [i for i in rho_p_parameters if i not in self.parameters] 

2114 ) 

2115 ) 

2116 ) 

2117 if "theta_jn" in self.parameters and "cos_theta_jn" not in self.parameters: 

2118 self._cos_angle("cos_theta_jn") 

2119 if "theta_jn" in self.parameters and "viewing_angle" not in self.parameters: 

2120 self._viewing_angle() 

2121 if "iota" in self.parameters and "cos_iota" not in self.parameters: 

2122 self._cos_angle("cos_iota") 

2123 remove_parameters = [ 

2124 "tilt_1_evolved", "tilt_2_evolved", "phi_12_evolved", 

2125 "spin_1z_evolved", "spin_2z_evolved", "reference_frequency", 

2126 "minimum_frequency" 

2127 ] 

2128 for param in remove_parameters: 

2129 if param in self.parameters: 

2130 ind = self.parameters.index(param) 

2131 self.parameters.remove(self.parameters[ind]) 

2132 for i in self.samples: 

2133 del i[ind]