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

1# Licensed under an MIT style license -- see LICENSE.md 

2 

3import copy 

4import numpy as np 

5from pathlib import Path 

6 

7from pesummary import conf 

8from pesummary.utils.decorators import set_docstring 

9from pesummary.utils.exceptions import EvolveSpinError 

10from pesummary.utils.utils import logger, import_error_msg 

11 

12try: 

13 import lalsimulation 

14except ImportError: 

15 logger.warning(import_error_msg.format("lalsimulation")) 

16try: 

17 import astropy 

18except ImportError: 

19 logger.warning(import_error_msg.format("astropy")) 

20 

21from .angles import * 

22from .cosmology import * 

23from .cosmology import _source_from_detector 

24from .mass import * 

25from .remnant import * 

26from .remnant import _final_from_initial_BBH 

27from .snr import * 

28from .snr import _ifo_snr 

29from .spins import * 

30from .tidal import * 

31from .tidal import _check_NSBH_approximant 

32from .time import * 

33 

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

35_conversion_doc = """ 

36 Class to calculate all possible derived quantities 

37 

38 Parameters 

39 ---------- 

40 data: dict, list 

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

42 samples. See the examples below for details 

43 extra_kwargs: dict, optional 

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

45 f_low: float, optional 

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

47 f_start: float, optional 

48 the starting frequency of the waveform 

49 f_ref: float, optional 

50 the reference frequency when spins are defined 

51 f_final: float, optional 

52 the final frequency to use when integrating over frequencies 

53 approximant: str, optional 

54 the approximant to use when evolving the spins 

55 approximant_flags: dict, optional 

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

57 evolve_spins_forwards: float/str, optional 

58 the final velocity to evolve the spins up to. 

59 evolve_spins_backwards: str, optional 

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

61 return_kwargs: Bool, optional 

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

63 about the conversion 

64 NRSur_fits: float/str, optional 

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

66 passed, the average NR fits are used instead 

67 multipole_snr: Bool, optional 

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

69 calculated from the posterior samples. 

70 samples. 

71 precessing_snr: Bool, optional 

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

73 psd: dict, optional 

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

75 to include in calculations 

76 waveform_fits: Bool, optional 

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

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

79 multi_process: int, optional 

80 number of cores to use to parallelize the computationally expensive 

81 conversions 

82 redshift_method: str, optional 

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

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

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

86 the redshift given N luminosity distance points. 

87 cosmology: str, optional 

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

89 distance samples. 

90 force_non_evolved: Bool, optional 

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

92 already exist in the input. Default False 

93 force_BBH_remnant_computation: Bool, optional 

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

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

96 Default False. 

97 force_BH_spin_evolution: Bool, optional 

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

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

100 Default False. 

101 disable_remnant: Bool, optional 

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

103 add_zero_spin: Bool, optional 

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

105 Default False. 

106 psd_default: str/pycbc.psd obj, optional 

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

108 regenerate: list, optional 

109 list of posterior distributions that you wish to regenerate 

110 return_dict: Bool, optional 

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

112 resume_file: str, optional 

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

114 is not used. Default None 

115 

116 Examples 

117 -------- 

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

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

120 the examples below: 

121 

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

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

124 

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

126 >>> samples = [10, 5] 

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

128 

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

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

131 

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

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

134 """ 

135 

136 

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

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

139 import os 

140 if resume_file is not None: 

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

142 return _Conversion.load_current_state(resume_file) 

143 logger.info( 

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

145 "checkpoint" 

146 ) 

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

148 

149 

150class _PickledConversion(object): 

151 pass 

152 

153 

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

155class _Conversion(object): 

156 @classmethod 

157 def load_current_state(cls, resume_file): 

158 """Load current state from a resume file 

159 

160 Parameters 

161 ---------- 

162 resume_file: str 

163 path to a resume file to restart conversion 

164 """ 

165 from pesummary.io import read 

166 logger.info( 

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

168 ) 

169 state = read(resume_file, checkpoint=True) 

170 return cls( 

171 state.parameters, state.samples, extra_kwargs=state.extra_kwargs, 

172 evolve_spins_forwards=state.evolve_spins_forwards, 

173 evolve_spins_backwards=state.evolve_spins_backwards, 

174 NRSur_fits=state.NRSurrogate, 

175 waveform_fits=state.waveform_fit, multi_process=state.multi_process, 

176 redshift_method=state.redshift_method, cosmology=state.cosmology, 

177 force_non_evolved=state.force_non_evolved, 

178 force_BBH_remnant_computation=state.force_remnant, 

179 disable_remnant=state.disable_remnant, 

180 add_zero_spin=state.add_zero_spin, regenerate=state.regenerate, 

181 return_kwargs=state.return_kwargs, return_dict=state.return_dict, 

182 resume_file=state.resume_file 

183 ) 

184 

185 def write_current_state(self): 

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

187 """ 

188 from pesummary.io import write 

189 state = _PickledConversion() 

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

191 setattr(state, key, value) 

192 

193 _path = Path(self.resume_file) 

194 write( 

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

196 filename=_path.name, overwrite=True 

197 ) 

198 logger.debug( 

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

200 ) 

201 

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

203 from pesummary.utils.samples_dict import SamplesDict 

204 from pesummary.utils.parameters import Parameters 

205 

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

207 base_replace = ( 

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

209 "the passed {}" 

210 ) 

211 if len(args) > 2: 

212 raise ValueError( 

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

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

215 ) 

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

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

218 samples = np.atleast_2d( 

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

220 ).tolist() 

221 else: 

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

223 parameters = Parameters(args[0]) 

224 else: 

225 parameters = args[0] 

226 samples = args[1] 

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

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

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

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

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

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

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

234 

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

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

237 logger.warning( 

238 base_replace.format( 

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

240 ) 

241 ) 

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

243 elif value is not None: 

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

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

246 logger.debug( 

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

248 param, extra_kwargs["meta_data"][param] 

249 ) 

250 ) 

251 else: 

252 logger.warning( 

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

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

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

256 ) 

257 ) 

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

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

260 ) 

261 

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

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

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

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

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

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

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

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

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

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

272 raise ValueError( 

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

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

275 "method of calculating the redshift" 

276 ) 

277 if isinstance(NRSurrogate, bool) and NRSurrogate: 

278 raise ValueError( 

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

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

281 "quantities" 

282 ) 

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

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

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

286 if disable_remnant and ( 

287 force_non_evolved or force_remnant 

288 or NRSurrogate or waveform_fits or evolve_spins_forwards 

289 ): 

290 _disable = [] 

291 if force_non_evolved: 

292 _disable.append("force_non_evolved") 

293 force_non_evolved = False 

294 if force_remnant: 

295 _disable.append("force_BBH_remnant_computation") 

296 force_remnant = False 

297 if NRSurrogate: 

298 _disable.append("NRSur_fits") 

299 NRSurrogate = False 

300 if waveform_fits: 

301 _disable.append("waveform_fits") 

302 waveform_fits = False 

303 if evolve_spins_forwards: 

304 _disable.append("evolve_spins_forwards") 

305 evolve_spins_forwards = False 

306 logger.warning( 

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

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

309 "calculated".format( 

310 " or ".join(_disable), 

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

312 ) 

313 ) 

314 if NRSurrogate and waveform_fits: 

315 raise ValueError( 

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

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

318 approximant 

319 ) 

320 ) 

321 if isinstance(evolve_spins_forwards, bool) and evolve_spins_forwards: 

322 raise ValueError( 

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

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

325 "evolve the spins up to the ISCO frequency" 

326 ) 

327 if not evolve_spins_forwards and (NRSurrogate or waveform_fits): 

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

329 logger.warning( 

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

331 "{} fits.".format( 

332 "NRSurrogate" if NRSurrogate else approximant 

333 ) 

334 ) 

335 elif evolve_spins_forwards and (NRSurrogate or waveform_fits): 

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

337 logger.warning( 

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

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

340 "NRSurrogate" if NRSurrogate else approximant 

341 ) 

342 ) 

343 else: 

344 logger.warning( 

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

346 approximant 

347 ) 

348 ) 

349 evolve_spins_forwards = False 

350 

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

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

353 

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

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

356 logger.warning( 

357 base_replace.format( 

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

359 ) 

360 ) 

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

362 elif freq is not None: 

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

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

365 logger.debug( 

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

367 name, extra_kwargs["meta_data"][name] 

368 ) 

369 ) 

370 else: 

371 logger.warning( 

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

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

374 name, conf.default_flow 

375 ) 

376 ) 

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

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

379 logger.warning( 

380 base_replace.format( 

381 "approximant_flags", extra_kwargs["meta_data"]["approximant_flags"], 

382 approximant_flags 

383 ) 

384 ) 

385 extra_kwargs["meta_data"]["approximant_flags"] = approximant_flags 

386 elif len(approximant_flags): 

387 extra_kwargs["meta_data"]["approximant_flags"] = approximant_flags 

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

389 logger.warning( 

390 base_replace.format( 

391 "approximant", extra_kwargs["meta_data"]["approximant"], 

392 approximant 

393 ) 

394 ) 

395 extra_kwargs["meta_data"]["approximant"] = approximant 

396 elif approximant is not None: 

397 extra_kwargs["meta_data"]["approximant"] = approximant 

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

399 logger.warning( 

400 base_replace.format( 

401 "f_ref", extra_kwargs["meta_data"]["f_ref"], f_ref 

402 ) 

403 ) 

404 extra_kwargs["meta_data"]["f_ref"] = f_ref 

405 elif f_ref is not None: 

406 extra_kwargs["meta_data"]["f_ref"] = f_ref 

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

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

409 if multi_process is not None: 

410 multi_process = int(multi_process) 

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

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

413 if psd is None: 

414 psd = {} 

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

416 raise ValueError( 

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

418 ) 

419 ifos = list(psd.keys()) 

420 pycbc_psd = copy.deepcopy(psd) 

421 if psd != {}: 

422 from pesummary.gw.file.psd import PSD 

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

424 for ifo in ifos: 

425 try: 

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

427 extra_kwargs["meta_data"]["f_low"], 

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

429 f_high_override=True 

430 ) 

431 except (ImportError, IndexError, ValueError): 

432 pass 

433 obj.__init__( 

434 parameters, samples, extra_kwargs, evolve_spins_forwards, NRSurrogate, 

435 waveform_fits, multi_process, regenerate, redshift_method, 

436 cosmology, force_non_evolved, force_remnant, 

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

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

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

440 pycbc_psd, psd_default, evolve_spins_backwards, force_evolve 

441 ) 

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

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

444 return [ 

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

446 obj.extra_kwargs 

447 ] 

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

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

450 elif return_kwargs: 

451 return obj.parameters, obj.samples, obj.extra_kwargs 

452 else: 

453 return obj.parameters, obj.samples 

454 

455 def __init__( 

456 self, parameters, samples, extra_kwargs, evolve_spins_forwards, NRSurrogate, 

457 waveform_fits, multi_process, regenerate, redshift_method, 

458 cosmology, force_non_evolved, force_remnant, add_zero_spin, 

459 disable_remnant, return_kwargs, return_dict, resume_file, multipole_snr, 

460 precessing_snr, psd, psd_default, evolve_spins_backwards, force_evolve 

461 ): 

462 self.parameters = parameters 

463 self.samples = samples 

464 self.extra_kwargs = extra_kwargs 

465 self.evolve_spins_forwards = evolve_spins_forwards 

466 self.evolve_spins_backwards = evolve_spins_backwards 

467 self.NRSurrogate = NRSurrogate 

468 self.waveform_fit = waveform_fits 

469 self.multi_process = multi_process 

470 self.regenerate = regenerate 

471 self.redshift_method = redshift_method 

472 self.cosmology = cosmology 

473 self.force_non_evolved = force_non_evolved 

474 self.force_remnant = force_remnant 

475 self.force_evolve = force_evolve 

476 self.disable_remnant = disable_remnant 

477 self.return_kwargs = return_kwargs 

478 self.return_dict = return_dict 

479 self.resume_file = resume_file 

480 self.multipole_snr = multipole_snr 

481 self.precessing_snr = precessing_snr 

482 self.psd = psd 

483 self.psd_default = psd_default 

484 self.non_precessing = False 

485 cond1 = any( 

486 param in self.parameters for param in 

487 conf.precessing_angles + conf.precessing_spins 

488 ) 

489 if not cond1: 

490 self.non_precessing = True 

491 if "chi_p" in self.parameters: 

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

493 if not np.any(_chi_p): 

494 logger.info( 

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

496 "non-precessing system" 

497 ) 

498 self.non_precessing = True 

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

500 samples = self.specific_parameter_samples(conf.precessing_spins) 

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

502 logger.info( 

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

504 "a non-precessing system" 

505 ) 

506 cond1 = self.non_precessing and evolve_spins_forwards 

507 cond2 = self.non_precessing and evolve_spins_backwards 

508 if cond1 or cond2: 

509 logger.info( 

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

511 "transformation required." 

512 ) 

513 self.evolve_spins_forwards = False 

514 self.evolve_spins_backwards = False 

515 if not self.non_precessing and multipole_snr: 

516 logger.warning( 

517 "The calculation for computing the SNR in subdominant " 

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

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

520 "results" 

521 ) 

522 if self.non_precessing and precessing_snr: 

523 logger.info( 

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

525 "conversion required." 

526 ) 

527 self.precessing_snr = False 

528 self.has_tidal = self._check_for_tidal_parameters() 

529 self.NSBH = self._check_for_NSBH_system() 

530 self.compute_remnant = not self.disable_remnant 

531 if self.has_tidal: 

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

533 logger.warning( 

534 "Posterior samples for tidal deformability found in the " 

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

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

537 "sensible results" 

538 ) 

539 elif self.evolve_spins_forwards or self.evolve_spins_backwards: 

540 logger.warning( 

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

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

543 "for BHs." 

544 ) 

545 self.evolve_spins_forwards = False 

546 self.evolve_spins_backwards = False 

547 

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

549 logger.warning( 

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

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

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

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

554 ) 

555 elif self.NSBH and self.compute_remnant: 

556 logger.warning( 

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

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

559 "fits to this system." 

560 ) 

561 self.waveform_fit = True 

562 elif force_remnant and self.compute_remnant: 

563 logger.warning( 

564 "Posterior samples for tidal deformability found in the " 

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

566 "This may not give sensible results." 

567 ) 

568 elif self.compute_remnant: 

569 logger.info( 

570 "Skipping remnant calculations as tidal deformability " 

571 "parameters found in the posterior table." 

572 ) 

573 self.compute_remnant = False 

574 if self.regenerate is not None: 

575 for param in self.regenerate: 

576 self.remove_posterior(param) 

577 self.add_zero_spin = add_zero_spin 

578 self.generate_all_posterior_samples() 

579 

580 def _check_for_tidal_parameters(self): 

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

582 """ 

583 from pesummary.gw.file.standard_names import tidal_params 

584 

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

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

587 return True 

588 _tidal_posterior = self.specific_parameter_samples( 

589 ["lambda_1", "lambda_2"] 

590 ) 

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

592 logger.warning( 

593 "Tidal deformability parameters found in the posterior " 

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

595 "system." 

596 ) 

597 return False 

598 return True 

599 return False 

600 

601 def _check_for_NSBH_system(self): 

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

603 system 

604 """ 

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

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

607 if not np.any(_lambda_2): 

608 logger.warning( 

609 "Posterior samples for lambda_2 found in the posterior " 

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

611 "system." 

612 ) 

613 return False 

614 return True 

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

616 _lambda_1, _lambda_2 = self.specific_parameter_samples( 

617 ["lambda_1", "lambda_2"] 

618 ) 

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

620 return False 

621 elif not np.any(_lambda_1): 

622 logger.warning( 

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

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

625 "Assuming this is an NSBH system." 

626 ) 

627 return True 

628 return False 

629 

630 def remove_posterior(self, parameter): 

631 if parameter in self.parameters: 

632 logger.info( 

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

634 ) 

635 ind = self.parameters.index(parameter) 

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

637 for i in self.samples: 

638 del i[ind] 

639 else: 

640 logger.info( 

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

642 "remove".format(parameter) 

643 ) 

644 

645 def _specific_parameter_samples(self, param): 

646 """Return the samples for a specific parameter 

647 

648 Parameters 

649 ---------- 

650 param: str 

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

652 """ 

653 if param == "empty": 

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

655 ind = self.parameters.index(param) 

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

657 return samples 

658 

659 def specific_parameter_samples(self, param): 

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

661 

662 Parameters 

663 ---------- 

664 param: list/str 

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

666 for 

667 """ 

668 if type(param) == list: 

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

670 else: 

671 samples = self._specific_parameter_samples(param) 

672 return samples 

673 

674 def append_data(self, parameter, samples): 

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

676 

677 Parameters 

678 ---------- 

679 parameter: str 

680 the name of the parameter you would like to append 

681 samples: list 

682 the list of samples that you would like to append 

683 """ 

684 if parameter not in self.parameters: 

685 self.parameters.append(parameter) 

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

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

688 if self.resume_file is not None: 

689 self.write_current_state() 

690 

691 def _mchirp_from_mchirp_source_z(self): 

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

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

694 self.append_data("chirp_mass", chirp_mass) 

695 

696 def _q_from_eta(self): 

697 samples = self.specific_parameter_samples("symmetric_mass_ratio") 

698 mass_ratio = q_from_eta(samples) 

699 self.append_data("mass_ratio", mass_ratio) 

700 

701 def _q_from_m1_m2(self): 

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

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

704 self.append_data("mass_ratio", mass_ratio) 

705 

706 def _invert_q(self): 

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

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

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

710 

711 def _invq_from_q(self): 

712 samples = self.specific_parameter_samples("mass_ratio") 

713 inverted_mass_ratio = invq_from_q(samples) 

714 self.append_data("inverted_mass_ratio", inverted_mass_ratio) 

715 

716 def _mchirp_from_mtotal_q(self): 

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

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

719 self.append_data("chirp_mass", chirp_mass) 

720 

721 def _m1_from_mchirp_q(self): 

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

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

724 self.append_data("mass_1", mass_1) 

725 

726 def _m2_from_mchirp_q(self): 

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

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

729 self.append_data("mass_2", mass_2) 

730 

731 def _m1_from_mtotal_q(self): 

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

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

734 self.append_data("mass_1", mass_1) 

735 

736 def _m2_from_mtotal_q(self): 

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

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

739 self.append_data("mass_2", mass_2) 

740 

741 def _reference_frequency(self): 

742 nsamples = len(self.samples) 

743 extra_kwargs = self.extra_kwargs["meta_data"] 

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

745 self.append_data( 

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

747 ) 

748 else: 

749 logger.warning( 

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

751 "as default") 

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

753 

754 def _mtotal_from_m1_m2(self): 

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

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

757 self.append_data("total_mass", m_total) 

758 

759 def _mchirp_from_m1_m2(self): 

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

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

762 self.append_data("chirp_mass", chirp_mass) 

763 

764 def _eta_from_m1_m2(self): 

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

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

767 self.append_data("symmetric_mass_ratio", eta) 

768 

769 def _phi_12_from_phi1_phi2(self): 

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

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

772 self.append_data("phi_12", phi_12) 

773 

774 def _phi1_from_spins(self): 

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

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

777 self.append_data("phi_1", phi_1) 

778 

779 def _phi2_from_spins(self): 

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

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

782 self.append_data("phi_2", phi_2) 

783 

784 def _spin_angles(self): 

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

786 "a_1", "a_2"] 

787 spin_angles_to_calculate = [ 

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

789 spin_components = [ 

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

791 "spin_2x", "spin_2y", "spin_2z", "reference_frequency"] 

792 samples = self.specific_parameter_samples(spin_components) 

793 if "phase" in self.parameters: 

794 spin_components.append("phase") 

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

796 else: 

797 logger.warning( 

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

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

800 ) 

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

802 angles = spin_angles( 

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

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

805 samples[10])