Coverage for bilby/core/result.py: 78%

1import datetime

2import inspect

3import json

4import os

5from collections import namedtuple

6from copy import copy

7from importlib import import_module

8from itertools import product

9import multiprocessing

10from functools import partial

11import numpy as np

12import pandas as pd

13import scipy.stats

15from . import utils

16from .utils import (

17 logger, infer_parameters_from_function,

18 check_directory_exists_and_if_not_mkdir,

19 latex_plot_format, safe_save_figure,

20 BilbyJsonEncoder, load_json,

21 move_old_file, get_version_information,

22 decode_bilby_json, docstring,

23 recursively_save_dict_contents_to_group,

24 recursively_load_dict_contents_from_group,

25 recursively_decode_bilby_json,

26 safe_file_dump,

27 random,

28)

29from .prior import Prior, PriorDict, DeltaFunction, ConditionalDeltaFunction

32EXTENSIONS = ["json", "hdf5", "h5", "pickle", "pkl"]

35def __eval_l(likelihood, params):

36 likelihood.parameters.update(params)

37 return likelihood.log_likelihood()

40def result_file_name(outdir, label, extension='json', gzip=False):

41 """ Returns the standard filename used for a result file

43 Parameters

44 ==========

45 outdir: str

46 Name of the output directory

47 label: str

48 Naming scheme of the output file

49 extension: str, optional

50 Whether to save as `hdf5`, `json`, or `pickle`

51 gzip: bool, optional

52 Set to True to append `.gz` to the extension for saving in gzipped format

54 Returns

55 =======

56 str: File name of the output file

57 """

58 if extension == 'pickle':

59 extension = 'pkl'

60 if extension in ['json', 'hdf5', 'pkl']:

61 if extension == 'json' and gzip:

62 return os.path.join(outdir, '{}_result.{}.gz'.format(label, extension))

63 else:

64 return os.path.join(outdir, '{}_result.{}'.format(label, extension))

65 else:

66 raise ValueError("Extension type {} not understood".format(extension))

69def _determine_file_name(filename, outdir, label, extension, gzip):

70 """ Helper method to determine the filename """

71 if filename is not None:

72 if isinstance(filename, os.PathLike):

73 # convert PathLike object to string

74 return str(filename)

75 else:

76 return filename

77 else:

78 if (outdir is None) and (label is None):

79 raise ValueError("No information given to load file")

80 else:

81 return result_file_name(outdir, label, extension, gzip)

84def read_in_result(filename=None, outdir=None, label=None, extension='json', gzip=False, result_class=None):

85 """ Reads in a stored bilby result object

87 Parameters

88 ==========

89 filename: str

90 Path to the file to be read (alternative to giving the outdir and label)

91 outdir, label, extension: str

92 Name of the output directory, label and extension used for the default

93 naming scheme.

94 result_class: bilby.core.result.Result, or child of

95 The result class to use. By default, `bilby.core.result.Result` is used,

96 but objects which inherit from this class can be given providing

97 additional methods.

98 """

99 filename = _determine_file_name(filename, outdir, label, extension, gzip)

100

101 if result_class is None:

102 result_class = Result

103 elif not issubclass(result_class, Result):

104 raise ValueError(f"Input result_class={result_class} not understood")

105

106 # Get the actual extension (may differ from the default extension if the filename is given)

107 extension = os.path.splitext(filename)[1].lstrip('.')

108 if extension == 'gz': # gzipped file

109 extension = os.path.splitext(os.path.splitext(filename)[0])[1].lstrip('.')

110

111 if 'json' in extension:

112 result = result_class.from_json(filename=filename)

113 elif ('hdf5' in extension) or ('h5' in extension):

114 result = result_class.from_hdf5(filename=filename)

115 elif ("pkl" in extension) or ("pickle" in extension):

116 result = result_class.from_pickle(filename=filename)

117 elif extension is None:

118 raise ValueError("No filetype extension provided")

119 else:

120 raise ValueError("Filetype {} not understood".format(extension))

121 return result

122

123

124def read_in_result_list(filename_list, invalid="warning"):

125 """ Read in a set of results

126

127 Parameters

128 ==========

129 filename_list: list

130 A list of filename paths

131 invalid: str (ignore, warning, error)

132 Behaviour if a file in filename_list is not a valid bilby result

133

134 Returns

135 -------

136 result_list: ResultList

137 A list of results

138 """

139 results_list = []

140 for filename in filename_list:

141 if (

142 not os.path.exists(filename)

143 and os.path.exists(f"{os.path.splitext(filename)[0]}.pkl")

144 ):

145 pickle_path = f"{os.path.splitext(filename)[0]}.pkl"

146 logger.warning(

147 f"Result file {filename} doesn't exist but {pickle_path} does. "

148 f"Using {pickle_path}."

149 )

150 filename = pickle_path

151 try:

152 results_list.append(read_in_result(filename=filename))

153 except Exception as e:

154 msg = f"Failed to read in file {filename} due to exception {e}"

155 if invalid == "error":

156 raise ResultListError(msg)

157 elif invalid == "warning":

158 logger.warning(msg)

159 return ResultList(results_list)

160

161

162def get_weights_for_reweighting(

163 result, new_likelihood=None, new_prior=None, old_likelihood=None,

164 old_prior=None, resume_file=None, n_checkpoint=5000, npool=1):

165 """ Calculate the weights for reweight()

166

167 See bilby.core.result.reweight() for help with the inputs

168

169 Returns

170 =======

171 ln_weights: array

172 An array of the natural-log weights

173 new_log_likelihood_array: array

174 An array of the natural-log likelihoods from the new likelihood

175 new_log_prior_array: array

176 An array of the natural-log priors

177 old_log_likelihood_array: array

178 An array of the natural-log likelihoods from the old likelihood

179 old_log_prior_array: array

180 An array of the natural-log priors

181 resume_file: string

182 filepath for the resume file which stores the weights

183 n_checkpoint: int

184 Number of samples to reweight before writing a resume file

185 """

186 from tqdm.auto import tqdm

187

188 nposterior = len(result.posterior)

189

190 if (resume_file is not None) and os.path.exists(resume_file):

191 old_log_likelihood_array, old_log_prior_array, new_log_likelihood_array, new_log_prior_array = \

192 np.genfromtxt(resume_file)

193

194 starting_index = np.argmin(np.abs(old_log_likelihood_array))

195 logger.info(f'Checkpoint resuming from {starting_index}.')

196

197 else:

198 old_log_likelihood_array = np.zeros(nposterior)

199 old_log_prior_array = np.zeros(nposterior)

200 new_log_likelihood_array = np.zeros(nposterior)

201 new_log_prior_array = np.zeros(nposterior)

202

203 starting_index = 0

204

205 dict_samples = [{key: sample[key] for key in result.posterior}

206 for _, sample in result.posterior.iterrows()]

207 n = len(dict_samples) - starting_index

208

209 # Helper function to compute likelihoods in parallel

210 def eval_pool(this_logl):

211 with multiprocessing.Pool(processes=npool) as pool:

212 chunksize = max(100, n // (2 * npool))

213 return list(tqdm(

214 pool.imap(partial(__eval_l, this_logl),

215 dict_samples[starting_index:], chunksize=chunksize),

216 desc='Computing likelihoods',

217 total=n)

218 )

219

220 if old_likelihood is None:

221 old_log_likelihood_array[starting_index:] = \

222 result.posterior["log_likelihood"][starting_index:].to_numpy()

223 else:

224 old_log_likelihood_array[starting_index:] = eval_pool(old_likelihood)

225

226 if new_likelihood is None:

227 # Don't perform likelihood reweighting (i.e. likelihood isn't updated)

228 new_log_likelihood_array[starting_index:] = old_log_likelihood_array[starting_index:]

229 else:

230 new_log_likelihood_array[starting_index:] = eval_pool(new_likelihood)

231

232 # Compute priors

233 for ii, sample in enumerate(tqdm(dict_samples[starting_index:],

234 desc='Computing priors',

235 total=n),

236 start=starting_index):

237 # prior calculation needs to not have prior or likelihood keys

238 ln_prior = sample.pop("log_prior", np.nan)

239 if "log_likelihood" in sample:

240 del sample["log_likelihood"]

241

242 if old_prior is not None:

243 old_log_prior_array[ii] = old_prior.ln_prob(sample)

244 else:

245 old_log_prior_array[ii] = ln_prior

246

247 if new_prior is not None:

248 new_log_prior_array[ii] = new_prior.ln_prob(sample)

249 else:

250 # Don't perform prior reweighting (i.e. prior isn't updated)

251 new_log_prior_array[ii] = old_log_prior_array[ii]

252

253 if (ii % (n_checkpoint) == 0) and (resume_file is not None):

254 checkpointed_index = np.argmin(np.abs(old_log_likelihood_array))

255 logger.info(f'Checkpointing with {checkpointed_index} samples')

256 np.savetxt(

257 resume_file,

258 [old_log_likelihood_array, old_log_prior_array, new_log_likelihood_array, new_log_prior_array])

259

260 ln_weights = (

261 new_log_likelihood_array + new_log_prior_array - old_log_likelihood_array - old_log_prior_array)

262

263 return ln_weights, new_log_likelihood_array, new_log_prior_array, old_log_likelihood_array, old_log_prior_array

264

265

266def rejection_sample(posterior, weights):

267 """ Perform rejection sampling on a posterior using weights

268

269 Parameters

270 ==========

271 posterior: pd.DataFrame or np.ndarray of shape (nsamples, nparameters)

272 The dataframe or array containing posterior samples

273 weights: np.ndarray

274 An array of weights

275

276 Returns

277 =======

278 reweighted_posterior: pd.DataFrame

279 The posterior resampled using rejection sampling

280

281 """

282 keep = weights > random.rng.uniform(0, max(weights), weights.shape)

283 return posterior[keep]

284

285

286def reweight(result, label=None, new_likelihood=None, new_prior=None,

287 old_likelihood=None, old_prior=None, conversion_function=None, npool=1,

288 verbose_output=False, resume_file=None, n_checkpoint=5000,

289 use_nested_samples=False):

290 """ Reweight a result to a new likelihood/prior using rejection sampling

291

292 Parameters

293 ==========

294 label: str, optional

295 An updated label to apply to the result object

296 new_likelihood: bilby.core.likelood.Likelihood, (optional)

297 If given, the new likelihood to reweight too. If not given, likelihood

298 reweighting is not applied

299 new_prior: bilby.core.prior.PriorDict, (optional)

300 If given, the new prior to reweight too. If not given, prior

301 reweighting is not applied

302 old_likelihood: bilby.core.likelihood.Likelihood, (optional)

303 If given, calculate the old likelihoods from this object. If not given,

304 the values stored in the posterior are used.

305 old_prior: bilby.core.prior.PriorDict, (optional)

306 If given, calculate the old prior from this object. If not given,

307 the values stored in the posterior are used.

308 conversion_function: function, optional

309 Function which adds in extra parameters to the data frame,

310 should take the data_frame, likelihood and prior as arguments.

311 npool: int, optional

312 Number of threads with which to execute the conversion function

313 verbose_output: bool, optional

314 Flag determining whether the weight array and associated prior and

315 likelihood evaluations are output as well as the result file

316 resume_file: string, optional

317 filepath for the resume file which stores the weights

318 n_checkpoint: int, optional

319 Number of samples to reweight before writing a resume file

320 use_nested_samples: bool, optional

321 If true reweight the nested samples instead. This can greatly improve reweighting efficiency, especially if the

322 target distribution has support beyond the proposal posterior distribution.

323

324 Returns

325 =======

326 result: bilby.core.result.Result

327 A copy of the result object with a reweighted posterior

328 new_log_likelihood_array: array, optional (if verbose_output=True)

329 An array of the natural-log likelihoods from the new likelihood

330 new_log_prior_array: array, optional (if verbose_output=True)

331 An array of the natural-log priors from the new likelihood

332 old_log_likelihood_array: array, optional (if verbose_output=True)

333 An array of the natural-log likelihoods from the old likelihood

334 old_log_prior_array: array, optional (if verbose_output=True)

335 An array of the natural-log priors from the old likelihood

336

337 """

338 from scipy.special import logsumexp

339

340 result = copy(result)

341

342 if use_nested_samples:

343 result.posterior = result.nested_samples

344

345 nposterior = len(result.posterior)

346 logger.info("Reweighting posterior with {} samples".format(nposterior))

347

348 ln_weights, new_log_likelihood_array, new_log_prior_array, old_log_likelihood_array, old_log_prior_array =\

349 get_weights_for_reweighting(

350 result, new_likelihood=new_likelihood, new_prior=new_prior,

351 old_likelihood=old_likelihood, old_prior=old_prior,

352 resume_file=resume_file, n_checkpoint=n_checkpoint, npool=npool)

353

354 if use_nested_samples:

355 ln_weights += np.log(result.posterior["weights"])

356

357 weights = np.exp(ln_weights)

358

359 # Overwrite the likelihood and prior evaluations

360 result.posterior["log_likelihood"] = new_log_likelihood_array

361 result.posterior["log_prior"] = new_log_prior_array

362

363 result.posterior = rejection_sample(result.posterior, weights=weights)

364 result.posterior = result.posterior.reset_index(drop=True)

365 logger.info("Rejection sampling resulted in {} samples".format(len(result.posterior)))

366 result.meta_data["reweighted_using_rejection_sampling"] = True

367

368 if use_nested_samples:

369 result.log_evidence += np.log(np.sum(weights))

370 else:

371 result.log_evidence += logsumexp(ln_weights) - np.log(nposterior)

372

373 if new_prior is not None:

374 for key, prior in new_prior.items():

375 result.priors[key] = prior

376

377 if conversion_function is not None:

378 data_frame = result.posterior

379 if "npool" in inspect.signature(conversion_function).parameters:

380 data_frame = conversion_function(data_frame, new_likelihood, new_prior, npool=npool)

381 else:

382 data_frame = conversion_function(data_frame, new_likelihood, new_prior)

383 result.posterior = data_frame

384

385 if label:

386 result.label = label

387 else:

388 result.label += "_reweighted"

389

390 if verbose_output:

391 return result, weights, new_log_likelihood_array, \

392 new_log_prior_array, old_log_likelihood_array, old_log_prior_array

393 else:

394 return result

395

396

397class Result(object):

398 def __init__(self, label='no_label', outdir='.', sampler=None,

399 search_parameter_keys=None, fixed_parameter_keys=None,

400 constraint_parameter_keys=None, priors=None,

401 sampler_kwargs=None, injection_parameters=None,

402 meta_data=None, posterior=None, samples=None,

403 nested_samples=None, log_evidence=np.nan,

404 log_evidence_err=np.nan, information_gain=np.nan,

405 log_noise_evidence=np.nan, log_bayes_factor=np.nan,

406 log_likelihood_evaluations=None,

407 log_prior_evaluations=None, sampling_time=None, nburn=None,

408 num_likelihood_evaluations=None, walkers=None,

409 max_autocorrelation_time=None, use_ratio=None,

410 parameter_labels=None, parameter_labels_with_unit=None,

411 version=None):

412 """ A class to store the results of the sampling run

413

414 Parameters

415 ==========

416 label, outdir, sampler: str

417 The label, output directory, and sampler used

418 search_parameter_keys, fixed_parameter_keys, constraint_parameter_keys: list

419 Lists of the search, constraint, and fixed parameter keys.

420 Elements of the list should be of type `str` and match the keys

421 of the `prior`

422 priors: dict, bilby.core.prior.PriorDict

423 A dictionary of the priors used in the run

424 sampler_kwargs: dict

425 Key word arguments passed to the sampler

426 injection_parameters: dict

427 A dictionary of the injection parameters

428 meta_data: dict

429 A dictionary of meta data to store about the run

430 posterior: pandas.DataFrame

431 A pandas data frame of the posterior

432 samples, nested_samples: array_like

433 An array of the output posterior samples and the unweighted samples

434 log_evidence, log_evidence_err, log_noise_evidence, log_bayes_factor: float

435 Natural log evidences

436 information_gain: float

437 The Kullback-Leibler divergence

438 log_likelihood_evaluations: array_like

439 The evaluations of the likelihood for each sample point

440 num_likelihood_evaluations: int

441 The number of times the likelihood function is called

442 log_prior_evaluations: array_like

443 The evaluations of the prior for each sample point

444 sampling_time: datetime.timedelta, float

445 The time taken to complete the sampling

446 nburn: int

447 The number of burn-in steps discarded for MCMC samplers

448 walkers: array_like

449 The samplers taken by a ensemble MCMC samplers

450 max_autocorrelation_time: float

451 The estimated maximum autocorrelation time for MCMC samplers

452 use_ratio: bool

453 A boolean stating whether the likelihood ratio, as opposed to the

454 likelihood was used during sampling

455 parameter_labels, parameter_labels_with_unit: list

456 Lists of the latex-formatted parameter labels

457 version: str,

458 Version information for software used to generate the result. Note,

459 this information is generated when the result object is initialized

460

461 Notes

462 =====

463 All sampling output parameters, e.g. the samples themselves are

464 typically not given at initialisation, but set at a later stage.

465

466 """

467

468 self.label = label

469 self.outdir = os.path.abspath(outdir)

470 self.sampler = sampler

471 self.search_parameter_keys = search_parameter_keys

472 self.fixed_parameter_keys = fixed_parameter_keys

473 self.constraint_parameter_keys = constraint_parameter_keys

474 self.parameter_labels = parameter_labels

475 self.parameter_labels_with_unit = parameter_labels_with_unit

476 self.priors = priors

477 self.sampler_kwargs = sampler_kwargs

478 self.meta_data = meta_data

479 self.injection_parameters = injection_parameters

480 self.posterior = posterior

481 self.samples = samples

482 if isinstance(nested_samples, dict):

483 nested_samples = pd.DataFrame(nested_samples)

484 self.nested_samples = nested_samples

485 self.walkers = walkers

486 self.nburn = nburn

487 self.use_ratio = use_ratio

488 self.log_evidence = log_evidence

489 self.log_evidence_err = log_evidence_err

490 self.information_gain = information_gain

491 self.log_noise_evidence = log_noise_evidence

492 self.log_bayes_factor = log_bayes_factor

493 self.log_likelihood_evaluations = log_likelihood_evaluations

494 self.log_prior_evaluations = log_prior_evaluations

495 self.num_likelihood_evaluations = num_likelihood_evaluations

496 if isinstance(sampling_time, float):

497 sampling_time = datetime.timedelta(seconds=sampling_time)

498 self.sampling_time = sampling_time

499 self.version = version

500 self.max_autocorrelation_time = max_autocorrelation_time

501

502 self.prior_values = None

503 self._kde = None

504

505 _load_doctstring = """ Read in a saved .{format} data file

506

507 Parameters

508 ==========

509 filename: str

510 If given, try to load from this filename

511 outdir, label: str

512 If given, use the default naming convention for saved results file

513

514 Returns

515 =======

516 result: bilby.core.result.Result

517

518 Raises

519 ======

520 ValueError: If no filename is given and either outdir or label is None

521 If no bilby.core.result.Result is found in the path

522

523 """

524

525 @staticmethod

526 @docstring(_load_doctstring.format(format="pickle"))

527 def from_pickle(filename=None, outdir=None, label=None):

528 filename = _determine_file_name(filename, outdir, label, 'hdf5', False)

529 import dill

530 with open(filename, "rb") as ff:

531 return dill.load(ff)

532

533 @classmethod

534 @docstring(_load_doctstring.format(format="hdf5"))

535 def from_hdf5(cls, filename=None, outdir=None, label=None):

536 import h5py

537 filename = _determine_file_name(filename, outdir, label, 'hdf5', False)

538 with h5py.File(filename, "r") as ff:

539 data = recursively_load_dict_contents_from_group(ff, '/')

540 data["posterior"] = pd.DataFrame(data["posterior"])

541 data["priors"] = PriorDict._get_from_json_dict(

542 json.loads(data["priors"], object_hook=decode_bilby_json)

543 )

544 try:

545 cls = getattr(import_module(data['__module__']), data['__name__'])

546 except ImportError:

547 logger.debug(

548 "Module {}.{} not found".format(data["__module__"], data["__name__"])

549 )

550 except KeyError:

551 logger.debug("No class specified, using base Result.")

552 for key in ["__module__", "__name__"]:

553 if key in data:

554 del data[key]

555 return cls(**data)

556

557 @classmethod

558 @docstring(_load_doctstring.format(format="json"))

559 def from_json(cls, filename=None, outdir=None, label=None, gzip=False):

560 from json.decoder import JSONDecodeError

561

562 filename = _determine_file_name(filename, outdir, label, 'json', gzip)

563

564 if os.path.isfile(filename):

565 try:

566 dictionary = load_json(filename, gzip)

567 except JSONDecodeError as e:

568 raise IOError(

569 "JSON failed to decode {} with message {}".format(filename, e)

570 )

571 try:

572 return cls(**dictionary)

573 except TypeError as e:

574 raise IOError("Unable to load dictionary, error={}".format(e))

575 else:

576 raise IOError("No result '{}' found".format(filename))

577

578 def __str__(self):

579 """Print a summary """

580 if getattr(self, 'posterior', None) is not None:

581 if getattr(self, 'log_noise_evidence', None) is not None:

582 return ("nsamples: {:d}\n"

583 "ln_noise_evidence: {:6.3f}\n"

584 "ln_evidence: {:6.3f} +/- {:6.3f}\n"

585 "ln_bayes_factor: {:6.3f} +/- {:6.3f}\n"

586 .format(len(self.posterior), self.log_noise_evidence, self.log_evidence,

587 self.log_evidence_err, self.log_bayes_factor,

588 self.log_evidence_err))

589 else:

590 return ("nsamples: {:d}\n"

591 "ln_evidence: {:6.3f} +/- {:6.3f}\n"

592 .format(len(self.posterior), self.log_evidence, self.log_evidence_err))

593 else:

594 return ''

595

596 @property

597 def meta_data(self):

598 return self._meta_data

599

600 @meta_data.setter

601 def meta_data(self, meta_data):

602 if meta_data is None:

603 meta_data = dict()

604 meta_data = recursively_decode_bilby_json(meta_data)

605 self._meta_data = meta_data

606

607 @property

608 def priors(self):

609 if self._priors is not None:

610 return self._priors

611 else:

612 raise ValueError('Result object has no priors')

613

614 @priors.setter

615 def priors(self, priors):

616 if isinstance(priors, dict):

617 if isinstance(priors, PriorDict):

618 self._priors = priors

619 else:

620 self._priors = PriorDict(priors)

621 if self.parameter_labels is None:

622 self.parameter_labels = [self.priors[k].latex_label for k in

623 self.search_parameter_keys]

624 if self.parameter_labels_with_unit is None:

625 self.parameter_labels_with_unit = [

626 self.priors[k].latex_label_with_unit for k in

627 self.search_parameter_keys]

628 elif priors is None:

629 self._priors = priors

630 self.parameter_labels = self.search_parameter_keys

631 self.parameter_labels_with_unit = self.search_parameter_keys

632 else:

633 raise ValueError("Input priors not understood")

634

635 @property

636 def samples(self):

637 """ An array of samples """

638 if self._samples is not None:

639 return self._samples

640 else:

641 raise ValueError("Result object has no stored samples")

642

643 @samples.setter

644 def samples(self, samples):

645 self._samples = samples

646

647 @property

648 def num_likelihood_evaluations(self):

649 """ number of likelihood evaluations """

650 if self._num_likelihood_evaluations is not None:

651 return self._num_likelihood_evaluations

652 else:

653 raise ValueError("Result object has no stored likelihood evaluations")

654

655 @num_likelihood_evaluations.setter

656 def num_likelihood_evaluations(self, num_likelihood_evaluations):

657 self._num_likelihood_evaluations = num_likelihood_evaluations

658

659 @property

660 def nested_samples(self):

661 """" An array of unweighted samples """

662 if self._nested_samples is not None:

663 return self._nested_samples

664 else:

665 raise ValueError("Result object has no stored nested samples")

666

667 @nested_samples.setter

668 def nested_samples(self, nested_samples):

669 self._nested_samples = nested_samples

670

671 @property

672 def walkers(self):

673 """" An array of the ensemble walkers """

674 if self._walkers is not None:

675 return self._walkers

676 else:

677 raise ValueError("Result object has no stored walkers")

678

679 @walkers.setter

680 def walkers(self, walkers):

681 self._walkers = walkers

682

683 @property

684 def nburn(self):

685 """" An array of the ensemble walkers """

686 if self._nburn is not None:

687 return self._nburn

688 else:

689 raise ValueError("Result object has no stored nburn")

690

691 @nburn.setter

692 def nburn(self, nburn):

693 self._nburn = nburn

694

695 @property

696 def posterior(self):

697 """ A pandas data frame of the posterior """

698 if self._posterior is not None:

699 return self._posterior

700 else:

701 raise ValueError("Result object has no stored posterior")

702

703 @posterior.setter

704 def posterior(self, posterior):

705 self._posterior = posterior

706

707 @property

708 def log_10_bayes_factor(self):

709 return self.log_bayes_factor / np.log(10)

710

711 @property

712 def log_10_evidence(self):

713 return self.log_evidence / np.log(10)

714

715 @property

716 def log_10_evidence_err(self):

717 return self.log_evidence_err / np.log(10)

718

719 @property

720 def log_10_noise_evidence(self):

721 return self.log_noise_evidence / np.log(10)

722

723 @property

724 def version(self):

725 return self._version

726

727 @version.setter

728 def version(self, version):

729 if version is None:

730 self._version = 'bilby={}'.format(utils.get_version_information())

731 else:

732 self._version = version

733

734 def _get_save_data_dictionary(self):

735 # This list defines all the parameters saved in the result object

736 save_attrs = [

737 'label', 'outdir', 'sampler', 'log_evidence', 'log_evidence_err',

738 'log_noise_evidence', 'log_bayes_factor', 'priors', 'posterior',

739 'injection_parameters', 'meta_data', 'search_parameter_keys',

740 'fixed_parameter_keys', 'constraint_parameter_keys',

741 'sampling_time', 'sampler_kwargs', 'use_ratio', 'information_gain',

742 'log_likelihood_evaluations', 'log_prior_evaluations',

743 'num_likelihood_evaluations', 'samples', 'nested_samples',

744 'walkers', 'nburn', 'parameter_labels', 'parameter_labels_with_unit',

745 'version']

746 dictionary = dict()

747 for attr in save_attrs:

748 try:

749 dictionary[attr] = getattr(self, attr)

750 except ValueError as e:

751 logger.debug("Unable to save {}, message: {}".format(attr, e))

752 pass

753 return dictionary

754

755 def save_to_file(self, filename=None, overwrite=False, outdir=None,

756 extension='json', gzip=False):

757 """

758

759 Writes the Result to a file.

760

761 Supported formats are: `json`, `hdf5`, `pickle`

762

763 Parameters

764 ==========

765 filename: optional,

766 Filename to write to (overwrites the default)

767 overwrite: bool, optional

768 Whether or not to overwrite an existing result file.

769 default=False

770 outdir: str, optional

771 Path to the outdir. Default is the one stored in the result object.

772 extension: str, optional {json, hdf5, pkl, pickle, True}

773 Determines the method to use to store the data (if True defaults

774 to json)

775 gzip: bool, optional

776 If true, and outputting to a json file, this will gzip the resulting

777 file and add '.gz' to the file extension.

778 """

779

780 if extension is True:

781 extension = "json"

782

783 outdir = self._safe_outdir_creation(outdir, self.save_to_file)

784 if filename is None:

785 filename = result_file_name(outdir, self.label, extension, gzip)

786

787 move_old_file(filename, overwrite)

788

789 # Convert the prior to a string representation for saving on disk

790 dictionary = self._get_save_data_dictionary()

791

792 # Convert callable sampler_kwargs to strings

793 if dictionary.get('sampler_kwargs', None) is not None:

794 for key in dictionary['sampler_kwargs']:

795 if hasattr(dictionary['sampler_kwargs'][key], '__call__'):

796 dictionary['sampler_kwargs'][key] = str(dictionary['sampler_kwargs'])

797

798 try:

799 # convert priors to JSON dictionary for both JSON and hdf5 files

800 if extension == 'json':

801 dictionary["priors"] = dictionary["priors"]._get_json_dict()

802 if gzip:

803 import gzip

804 # encode to a string

805 json_str = json.dumps(dictionary, cls=BilbyJsonEncoder).encode('utf-8')

806 with gzip.GzipFile(filename, 'w') as file:

807 file.write(json_str)

808 else:

809 with open(filename, 'w') as file:

810 json.dump(dictionary, file, indent=2, cls=BilbyJsonEncoder)

811 elif extension == 'hdf5':

812 import h5py

813 dictionary["__module__"] = self.__module__

814 dictionary["__name__"] = self.__class__.__name__

815 with h5py.File(filename, 'w') as h5file:

816 recursively_save_dict_contents_to_group(h5file, '/', dictionary)

817 elif extension == 'pkl':

818 safe_file_dump(self, filename, "dill")

819 else:

820 raise ValueError("Extension type {} not understood".format(extension))

821 except Exception as e:

822 filename = ".".join(filename.split(".")[:-1]) + ".pkl"

823 safe_file_dump(self, filename, "dill")

824 logger.error(

825 "\n\nSaving the data has failed with the following message:\n"

826 "{}\nData has been dumped to {}.\n\n".format(e, filename)

827 )

828

829 def save_posterior_samples(self, filename=None, outdir=None, label=None):

830 """ Saves posterior samples to a file

831

832 Generates a .dat file containing the posterior samples and auxiliary

833 data saved in the posterior. Note, strings in the posterior are

834 removed while complex numbers will be given as absolute values with

835 abs appended to the column name

836

837 Parameters

838 ==========

839 filename: str

840 Alternative filename to use. Defaults to

841 outdir/label_posterior_samples.dat

842 outdir, label: str

843 Alternative outdir and label to use

844

845 """

846 if filename is None:

847 if label is None:

848 label = self.label

849 outdir = self._safe_outdir_creation(outdir, self.save_posterior_samples)

850 filename = '{}/{}_posterior_samples.dat'.format(outdir, label)

851 else:

852 outdir = os.path.dirname(filename)

853 self._safe_outdir_creation(outdir, self.save_posterior_samples)

854

855 # Drop non-numeric columns

856 df = self.posterior.select_dtypes([np.number]).copy()

857

858 # Convert complex columns to abs

859 for key in df.keys():

860 if np.any(np.iscomplex(df[key])):

861 complex_term = df.pop(key)

862 df.loc[:, key + "_abs"] = np.abs(complex_term)

863 df.loc[:, key + "_angle"] = np.angle(complex_term)

864

865 logger.info("Writing samples file to {}".format(filename))

866 df.to_csv(filename, index=False, header=True, sep=' ')

867

868 def get_latex_labels_from_parameter_keys(self, keys):

869 """ Returns a list of latex_labels corresponding to the given keys

870

871 Parameters

872 ==========

873 keys: list

874 List of strings corresponding to the desired latex_labels

875

876 Returns

877 =======

878 list: The desired latex_labels

879

880 """

881 latex_labels = []

882 for key in keys:

883 if key in self.search_parameter_keys:

884 idx = self.search_parameter_keys.index(key)

885 label = self.parameter_labels_with_unit[idx]

886 elif key in self.parameter_labels:

887 label = key

888 else:

889 label = None

890 logger.debug(

891 'key {} not a parameter label or latex label'.format(key)

892 )

893 if label is None:

894 label = key.replace("_", " ")

895 latex_labels.append(label)

896 return latex_labels

897

898 @property

899 def covariance_matrix(self):

900 """ The covariance matrix of the samples the posterior """

901 samples = self.posterior[self.search_parameter_keys].values

902 return np.cov(samples.T)

903

904 @property

905 def posterior_volume(self):

906 """ The posterior volume """

907 if self.covariance_matrix.ndim == 0:

908 return np.sqrt(self.covariance_matrix)

909 else:

910 return 1 / np.sqrt(np.abs(np.linalg.det(

911 1 / self.covariance_matrix)))

912

913 @staticmethod

914 def prior_volume(priors):

915 """ The prior volume, given a set of priors """

916 return np.prod([priors[k].maximum - priors[k].minimum for k in priors])

917

918 def occam_factor(self, priors):

919 """ The Occam factor,

920

921 See Chapter 28, `Mackay "Information Theory, Inference, and Learning

922 Algorithms" <http://www.inference.org.uk/itprnn/book.html>`_ Cambridge

923 University Press (2003).

924

925 """

926 return self.posterior_volume / self.prior_volume(priors)

927

928 @property

929 def bayesian_model_dimensionality(self):

930 """ Characterises how many parameters are effectively constraint by the data

931

932 See <https://arxiv.org/abs/1903.06682>

933

934 Returns

935 =======

936 float: The model dimensionality

937 """

938 return 2 * (np.mean(self.posterior['log_likelihood']**2) -

939 np.mean(self.posterior['log_likelihood'])**2)

940

941 def get_one_dimensional_median_and_error_bar(self, key, fmt='.2f',

942 quantiles=(0.16, 0.84)):

943 """ Calculate the median and error bar for a given key

944

945 Parameters

946 ==========

947 key: str

948 The parameter key for which to calculate the median and error bar

949 fmt: str, ('.2f')

950 A format string

951 quantiles: list, tuple

952 A length-2 tuple of the lower and upper-quantiles to calculate

953 the errors bars for.

954

955 Returns

956 =======

957 summary: namedtuple

958 An object with attributes, median, lower, upper and string

959

960 """

961 summary = namedtuple('summary', ['median', 'lower', 'upper', 'string'])

962

963 if len(quantiles) != 2:

964 raise ValueError("quantiles must be of length 2")

965

966 quants_to_compute = np.array([quantiles[0], 0.5, quantiles[1]])

967 quants = np.percentile(self.posterior[key], quants_to_compute * 100)

968 summary.median = quants[1]

969 summary.plus = quants[2] - summary.median

970 summary.minus = summary.median - quants[0]

971

972 fmt = "{{0:{0}}}".format(fmt).format

973 string_template = r"${{{0}}}_{{-{1}}}^{{+{2}}}$"

974 summary.string = string_template.format(

975 fmt(summary.median), fmt(summary.minus), fmt(summary.plus))

976 return summary

977

978 @latex_plot_format

979 def plot_single_density(self, key, prior=None, cumulative=False,

980 title=None, truth=None, save=True,

981 file_base_name=None, bins=50, label_fontsize=16,

982 title_fontsize=16, quantiles=(0.16, 0.84), dpi=300):

983 """ Plot a 1D marginal density, either probability or cumulative.

984

985 Parameters

986 ==========

987 key: str

988 Name of the parameter to plot

989 prior: {bool (True), bilby.core.prior.Prior}

990 If true, add the stored prior probability density function to the

991 one-dimensional marginal distributions. If instead a Prior

992 is provided, this will be plotted.

993 cumulative: bool

994 If true plot the CDF

995 title: bool

996 If true, add 1D title of the median and (by default 1-sigma)

997 error bars. To change the error bars, pass in the quantiles kwarg.

998 See method `get_one_dimensional_median_and_error_bar` for further

999 details). If `quantiles=None` is passed in, no title is added.

1000 truth: {bool, float}

1001 If true, plot self.injection_parameters[parameter].

1002 If float, plot this value.

1003 save: bool:

1004 If true, save plot to disk.

1005 file_base_name: str, optional

1006 If given, the base file name to use (by default `outdir/label_` is

1007 used)

1008 bins: int

1009 The number of histogram bins

1010 label_fontsize, title_fontsize: int

1011 The fontsizes for the labels and titles

1012 quantiles: tuple

1013 A length-2 tuple of the lower and upper-quantiles to calculate

1014 the errors bars for.

1015 dpi: int

1016 Dots per inch resolution of the plot

1017

1018 Returns

1019 =======

1020 figure: matplotlib.pyplot.figure

1021 A matplotlib figure object

1022 """

1023 import matplotlib.pyplot as plt

1024 logger.info('Plotting {} marginal distribution'.format(key))

1025 label = self.get_latex_labels_from_parameter_keys([key])[0]

1026 fig, ax = plt.subplots()

1027 try:

1028 ax.hist(self.posterior[key].values, bins=bins, density=True,

1029 histtype='step', cumulative=cumulative)

1030 except ValueError as e:

1031 logger.info(

1032 'Failed to generate 1d plot for {}, error message: {}'

1033 .format(key, e))

1034 return

1035 ax.set_xlabel(label, fontsize=label_fontsize)

1036 if truth is not None:

1037 ax.axvline(truth, ls='-', color='orange')

1038

1039 summary = self.get_one_dimensional_median_and_error_bar(

1040 key, quantiles=quantiles)

1041 ax.axvline(summary.median - summary.minus, ls='--', color='C0')

1042 ax.axvline(summary.median + summary.plus, ls='--', color='C0')

1043 if title:

1044 ax.set_title(summary.string, fontsize=title_fontsize)

1045

1046 if isinstance(prior, Prior):

1047 theta = np.linspace(ax.get_xlim()[0], ax.get_xlim()[1], 300)

1048 if cumulative is False:

1049 ax.plot(theta, prior.prob(theta), color='C2')

1050 else:

1051 ax.plot(theta, prior.cdf(theta), color='C2')

1052

1053 if save:

1054 fig.tight_layout()

1055 if cumulative:

1056 file_name = file_base_name + key + '_cdf'

1057 else:

1058 file_name = file_base_name + key + '_pdf'

1059 safe_save_figure(fig=fig, filename=file_name, dpi=dpi)

1060 plt.close(fig)

1061 else:

1062 return fig

1063

1064 def plot_marginals(self, parameters=None, priors=None, titles=True,

1065 file_base_name=None, bins=50, label_fontsize=16,

1066 title_fontsize=16, quantiles=(0.16, 0.84), dpi=300,

1067 outdir=None):

1068 """ Plot 1D marginal distributions

1069

1070 Parameters

1071 ==========

1072 parameters: (list, dict), optional

1073 If given, either a list of the parameter names to include, or a

1074 dictionary of parameter names and their "true" values to plot.

1075 priors: {bool (False), bilby.core.prior.PriorDict}

1076 If true, add the stored prior probability density functions to the

1077 one-dimensional marginal distributions. If instead a PriorDict

1078 is provided, this will be plotted.

1079 titles: bool

1080 If true, add 1D titles of the median and (by default 1-sigma)

1081 error bars. To change the error bars, pass in the quantiles kwarg.

1082 See method `get_one_dimensional_median_and_error_bar` for further

1083 details). If `quantiles=None` is passed in, no title is added.

1084 file_base_name: str, optional

1085 If given, the base file name to use (by default `outdir/label_` is

1086 used)

1087 bins: int

1088 The number of histogram bins

1089 label_fontsize, title_fontsize: int

1090 The font sizes for the labels and titles

1091 quantiles: tuple

1092 A length-2 tuple of the lower and upper-quantiles to calculate

1093 the errors bars for.

1094 dpi: int

1095 Dots per inch resolution of the plot

1096 outdir: str, optional

1097 Path to the outdir. Default is the one store in the result object.

1098

1099 Returns

1100 =======

1101 """

1102 if isinstance(parameters, dict):

1103 plot_parameter_keys = list(parameters.keys())

1104 truths = parameters

1105 elif parameters is None:

1106 plot_parameter_keys = self.posterior.keys()

1107 if self.injection_parameters is None:

1108 truths = dict()

1109 else:

1110 truths = self.injection_parameters

1111 else:

1112 plot_parameter_keys = list(parameters)

1113 if self.injection_parameters is None:

1114 truths = dict()

1115 else:

1116 truths = self.injection_parameters

1117

1118 if file_base_name is None:

1119 outdir = self._safe_outdir_creation(outdir, self.plot_marginals)

1120 file_base_name = '{}/{}_1d/'.format(outdir, self.label)

1121 check_directory_exists_and_if_not_mkdir(file_base_name)

1122

1123 if priors is True:

1124 priors = getattr(self, 'priors', dict())

1125 elif isinstance(priors, dict):

1126 pass

1127 elif priors in [False, None]:

1128 priors = dict()

1129 else:

1130 raise ValueError('Input priors={} not understood'.format(priors))

1131

1132 for i, key in enumerate(plot_parameter_keys):

1133 if not isinstance(self.posterior[key].values[0], float):

1134 continue

1135 prior = priors.get(key, None)

1136 truth = truths.get(key, None)

1137 for cumulative in [False, True]:

1138 self.plot_single_density(

1139 key, prior=prior, cumulative=cumulative, title=titles,

1140 truth=truth, save=True, file_base_name=file_base_name,

1141 bins=bins, label_fontsize=label_fontsize, dpi=dpi,

1142 title_fontsize=title_fontsize, quantiles=quantiles)

1143

1144 @latex_plot_format

1145 def plot_corner(self, parameters=None, priors=None, titles=True, save=True,

1146 filename=None, dpi=300, **kwargs):

1147 """ Plot a corner-plot

1148

1149 Parameters

1150 ==========

1151 parameters: (list, dict), optional

1152 If given, either a list of the parameter names to include, or a

1153 dictionary of parameter names and their "true" values to plot.

1154 priors: {bool (False), bilby.core.prior.PriorDict}

1155 If true, add the stored prior probability density functions to the

1156 one-dimensional marginal distributions. If instead a PriorDict

1157 is provided, this will be plotted.

1158 titles: bool

1159 If true, add 1D titles of the median and (by default 1-sigma)

1160 error bars. To change the error bars, pass in the quantiles kwarg.

1161 See method `get_one_dimensional_median_and_error_bar` for further

1162 details). If `quantiles=None` is passed in, no title is added.

1163 save: bool, optional

1164 If true, save the image using the given label and outdir

1165 filename: str, optional

1166 If given, overwrite the default filename

1167 dpi: int, optional

1168 Dots per inch resolution of the plot

1169 **kwargs:

1170 Other keyword arguments are passed to `corner.corner`. We set some

1171 defaults to improve the basic look and feel, but these can all be

1172 overridden. Also optional an 'outdir' argument which can be used

1173 to override the outdir set by the absolute path of the result object.

1174

1175 Notes

1176 =====

1177 The generation of the corner plot themselves is done by the corner

1178 python module, see https://corner.readthedocs.io for more

1179 information.

1180

1181 Truth-lines can be passed in in several ways. Either as the values

1182 of the parameters dict, or a list via the `truths` kwarg. If

1183 injection_parameters where given to run_sampler, these will auto-

1184 matically be added to the plot. This behaviour can be stopped by

1185 adding truths=False.

1186

1187 Returns

1188 =======

1189 fig:

1190 A matplotlib figure instance

1191

1192 """

1193 import corner

1194 import matplotlib.pyplot as plt

1195

1196 # If in testing mode, not corner plots are generated

1197 if utils.command_line_args.bilby_test_mode:

1198 return

1199

1200 defaults_kwargs = dict(

1201 bins=50, smooth=0.9,

1202 title_kwargs=dict(fontsize=16), color='#0072C1',

1203 truth_color='tab:orange', quantiles=[0.16, 0.84],

1204 levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.)),

1205 plot_density=False, plot_datapoints=True, fill_contours=True,

1206 max_n_ticks=3)

1207

1208 if 'lionize' in kwargs and kwargs['lionize'] is True:

1209 defaults_kwargs['truth_color'] = 'tab:blue'

1210 defaults_kwargs['color'] = '#FF8C00'

1211

1212 label_kwargs_defaults = dict(fontsize=16)

1213 hist_kwargs_defaults = dict(density=True)

1214

1215 label_kwargs_input = kwargs.get("label_kwargs", dict())

1216 hist_kwargs_input = kwargs.get("hist_kwargs", dict())

1217

1218 label_kwargs_defaults.update(label_kwargs_input)

1219 hist_kwargs_defaults.update(hist_kwargs_input)

1220

1221 defaults_kwargs.update(kwargs)

1222 kwargs = defaults_kwargs

1223

1224 kwargs["label_kwargs"] = label_kwargs_defaults

1225 kwargs["hist_kwargs"] = hist_kwargs_defaults

1226

1227 # Handle if truths was passed in

1228 if 'truth' in kwargs:

1229 kwargs['truths'] = kwargs.pop('truth')

1230 if "truths" in kwargs:

1231 truths = kwargs.get('truths')

1232 if isinstance(parameters, list) and isinstance(truths, list):

1233 if len(parameters) != len(truths):

1234 raise ValueError(

1235 "Length of parameters and truths don't match")

1236 elif isinstance(truths, dict) and parameters is None:

1237 parameters = kwargs.pop('truths')

1238 elif isinstance(truths, bool):

1239 pass

1240 elif truths is None:

1241 kwargs["truths"] = False

1242 else:

1243 raise ValueError(

1244 "Combination of parameters and truths not understood")

1245

1246 # If injection parameters where stored, use these as parameter values

1247 # but do not overwrite input parameters (or truths)

1248 cond1 = getattr(self, 'injection_parameters', None) is not None

1249 cond2 = parameters is None

1250 cond3 = bool(kwargs.get("truths", True))

1251 if cond1 and cond2 and cond3:

1252 parameters = {

1253 key: self.injection_parameters.get(key, np.nan)

1254 for key in self.search_parameter_keys

1255 }

1256

1257 # If parameters is a dictionary, use the keys to determine which

1258 # parameters to plot and the values as truths.

1259 if isinstance(parameters, dict):

1260 plot_parameter_keys = list(parameters.keys())

1261 kwargs['truths'] = list(parameters.values())

1262 elif parameters is None:

1263 plot_parameter_keys = self.search_parameter_keys

1264 else:

1265 plot_parameter_keys = list(parameters)

1266

1267 # Get latex formatted strings for the plot labels

1268 kwargs['labels'] = kwargs.get(

1269 'labels', self.get_latex_labels_from_parameter_keys(

1270 plot_parameter_keys))

1271

1272 kwargs["labels"] = sanity_check_labels(kwargs["labels"])

1273

1274 # Unless already set, set the range to include all samples

1275 # This prevents ValueErrors being raised for parameters with no range

1276 kwargs['range'] = kwargs.get('range', [1] * len(plot_parameter_keys))

1277

1278 # Remove truths if it is a bool

1279 if isinstance(kwargs.get('truths'), bool):

1280 kwargs.pop('truths')

1281

1282 # Create the data array to plot and pass everything to corner

1283 xs = self.posterior[plot_parameter_keys].values

1284 if len(plot_parameter_keys) > 1:

1285 fig = corner.corner(xs, **kwargs)

1286 else:

1287 ax = kwargs.get("ax", plt.subplot())

1288 ax.hist(xs, bins=kwargs["bins"], color=kwargs["color"],

1289 histtype="step", **kwargs["hist_kwargs"])

1290 ax.set_xlabel(kwargs["labels"][0])

1291 fig = plt.gcf()

1292

1293 axes = fig.get_axes()

1294

1295 # Add the titles

1296 if titles and kwargs.get('quantiles', None) is not None:

1297 for i, par in enumerate(plot_parameter_keys):

1298 ax = axes[i + i * len(plot_parameter_keys)]

1299 if ax.title.get_text() == '':

1300 ax.set_title(self.get_one_dimensional_median_and_error_bar(

1301 par, quantiles=kwargs['quantiles']).string,

1302 **kwargs['title_kwargs'])

1303

1304 # Add priors to the 1D plots

1305 if priors is True:

1306 priors = getattr(self, 'priors', False)

1307 if isinstance(priors, dict):

1308 for i, par in enumerate(plot_parameter_keys):

1309 ax = axes[i + i * len(plot_parameter_keys)]

1310 theta = np.linspace(ax.get_xlim()[0], ax.get_xlim()[1], 300)

1311 ax.plot(theta, priors[par].prob(theta), color='C2')

1312 elif priors in [False, None]:

1313 pass

1314 else:

1315 raise ValueError('Input priors={} not understood'.format(priors))

1316

1317 if save:

1318 if filename is None:

1319 outdir = self._safe_outdir_creation(kwargs.get('outdir'), self.plot_corner)

1320 filename = '{}/{}_corner.png'.format(outdir, self.label)

1321 logger.debug('Saving corner plot to {}'.format(filename))

1322 safe_save_figure(fig=fig, filename=filename, dpi=dpi)

1323 plt.close(fig)

1324

1325 return fig

1326

1327 @latex_plot_format

1328 def plot_walkers(self, **kwargs):

1329 """ Method to plot the trace of the walkers in an ensemble MCMC plot """

1330 import matplotlib.pyplot as plt

1331 if hasattr(self, 'walkers') is False:

1332 logger.warning("Cannot plot_walkers as no walkers are saved")

1333 return

1334

1335 if utils.command_line_args.bilby_test_mode:

1336 return

1337

1338 nwalkers, nsteps, ndim = self.walkers.shape

1339 idxs = np.arange(nsteps)

1340 fig, axes = plt.subplots(nrows=ndim, figsize=(6, 3 * ndim))

1341 walkers = self.walkers[:, :, :]

1342 parameter_labels = sanity_check_labels(self.parameter_labels)

1343 for i, ax in enumerate(axes):

1344 ax.plot(idxs[:self.nburn + 1], walkers[:, :self.nburn + 1, i].T,

1345 lw=0.1, color='r')

1346 ax.set_ylabel(parameter_labels[i])

1347

1348 for i, ax in enumerate(axes):

1349 ax.plot(idxs[self.nburn:], walkers[:, self.nburn:, i].T, lw=0.1,

1350 color='k')

1351 ax.set_ylabel(parameter_labels[i])

1352

1353 fig.tight_layout()

1354 outdir = self._safe_outdir_creation(kwargs.get('outdir'), self.plot_walkers)

1355 filename = '{}/{}_walkers.png'.format(outdir, self.label)

1356 logger.debug('Saving walkers plot to {}'.format('filename'))

1357 safe_save_figure(fig=fig, filename=filename)

1358 plt.close(fig)

1359

1360 @latex_plot_format

1361 def plot_with_data(self, model, x, y, ndraws=1000, npoints=1000,

1362 xlabel=None, ylabel=None, data_label='data',

1363 data_fmt='o', draws_label=None, filename=None,

1364 maxl_label='max likelihood', dpi=300, outdir=None):

1365 """ Generate a figure showing the data and fits to the data

1366

1367 Parameters

1368 ==========

1369 model: function

1370 A python function which when called as `model(x, **kwargs)` returns

1371 the model prediction (here `kwargs` is a dictionary of key-value

1372 pairs of the model parameters.

1373 x, y: np.ndarray

1374 The independent and dependent data to plot

1375 ndraws: int

1376 Number of draws from the posterior to plot

1377 npoints: int

1378 Number of points used to plot the smoothed fit to the data

1379 xlabel, ylabel: str

1380 Labels for the axes

1381 data_label, draws_label, maxl_label: str

1382 Label for the data, draws, and max likelihood legend

1383 data_fmt: str

1384 Matpltolib fmt code, defaults to `'-o'`

1385 dpi: int

1386 Passed to `plt.savefig`

1387 filename: str

1388 If given, the filename to use. Otherwise, the filename is generated

1389 from the outdir and label attributes.

1390 outdir: str, optional

1391 Path to the outdir. Default is the one store in the result object.

1392

1393 """

1394 import matplotlib.pyplot as plt

1395

1396 # Determine model_posterior, the subset of the full posterior which

1397 # should be passed into the model

1398 model_keys = infer_parameters_from_function(model)

1399 model_posterior = self.posterior[model_keys]

1400

1401 xsmooth = np.linspace(np.min(x), np.max(x), npoints)

1402 fig, ax = plt.subplots()

1403 logger.info('Plotting {} draws'.format(ndraws))

1404 for _ in range(ndraws):

1405 s = model_posterior.sample().to_dict('records')[0]

1406 ax.plot(xsmooth, model(xsmooth, **s), alpha=0.25, lw=0.1, color='r',

1407 label=draws_label)

1408 try:

1409 if all(~np.isnan(self.posterior.log_likelihood)):

1410 logger.info('Plotting maximum likelihood')

1411 s = model_posterior.iloc[self.posterior.log_likelihood.idxmax()]

1412 ax.plot(xsmooth, model(xsmooth, **s), lw=1, color='k',

1413 label=maxl_label)

1414 except (AttributeError, TypeError):

1415 logger.debug(

1416 "No log likelihood values stored, unable to plot max")

1417

1418 ax.plot(x, y, data_fmt, markersize=2, label=data_label)

1419

1420 if xlabel is not None:

1421 ax.set_xlabel(xlabel)

1422 if ylabel is not None:

1423 ax.set_ylabel(ylabel)

1424

1425 handles, labels = plt.gca().get_legend_handles_labels()

1426 by_label = dict(zip(labels, handles))

1427 plt.legend(by_label.values(), by_label.keys())

1428 ax.legend(numpoints=3)

1429 fig.tight_layout()

1430 if filename is None:

1431 outdir = self._safe_outdir_creation(outdir, self.plot_with_data)

1432 filename = '{}/{}_plot_with_data'.format(outdir, self.label)

1433 safe_save_figure(fig=fig, filename=filename, dpi=dpi)

1434 plt.close(fig)

1435

1436 @staticmethod

1437 def _add_prior_fixed_values_to_posterior(posterior, priors):

1438 if priors is None:

1439 return posterior

1440 for key in priors:

1441 if isinstance(priors[key], DeltaFunction) and \

1442 not isinstance(priors[key], ConditionalDeltaFunction):

1443 posterior[key] = priors[key].peak

1444 elif isinstance(priors[key], float):

1445 posterior[key] = priors[key]

1446 return posterior

1447

1448 def samples_to_posterior(self, likelihood=None, priors=None,

1449 conversion_function=None, npool=1):

1450 """

1451 Convert array of samples to posterior (a Pandas data frame)

1452

1453 Also applies the conversion function to any stored posterior

1454

1455 Parameters

1456 ==========

1457 likelihood: bilby.likelihood.GravitationalWaveTransient, optional

1458 GravitationalWaveTransient likelihood used for sampling.

1459 priors: bilby.prior.PriorDict, optional

1460 Dictionary of prior object, used to fill in delta function priors.

1461 conversion_function: function, optional

1462 Function which adds in extra parameters to the data frame,

1463 should take the data_frame, likelihood and prior as arguments.

1464 """

1465

1466 data_frame = pd.DataFrame(

1467 self.samples, columns=self.search_parameter_keys)

1468 data_frame = self._add_prior_fixed_values_to_posterior(

1469 data_frame, priors)

1470 data_frame['log_likelihood'] = getattr(

1471 self, 'log_likelihood_evaluations', np.nan)

1472 if self.log_prior_evaluations is None and priors is not None:

1473 data_frame['log_prior'] = priors.ln_prob(

1474 dict(data_frame[self.search_parameter_keys]), axis=0)

1475 else:

1476 data_frame['log_prior'] = self.log_prior_evaluations

1477

1478 if conversion_function is not None:

1479 if "npool" in inspect.signature(conversion_function).parameters:

1480 data_frame = conversion_function(data_frame, likelihood, priors, npool=npool)

1481 else:

1482 data_frame = conversion_function(data_frame, likelihood, priors)

1483 self.posterior = data_frame

1484

1485 def calculate_prior_values(self, priors):

1486 """

1487 Evaluate prior probability for each parameter for each sample.

1488

1489 Parameters

1490 ==========

1491 priors: dict, PriorDict

1492 Prior distributions

1493 """

1494 self.prior_values = pd.DataFrame()

1495 for key in priors:

1496 if key in self.posterior.keys():

1497 if isinstance(priors[key], DeltaFunction):

1498 continue

1499 else:

1500 self.prior_values[key]\

1501 = priors[key].prob(self.posterior[key].values)

1502

1503 def get_all_injection_credible_levels(self, keys=None, weights=None):

1504 """

1505 Get credible levels for all parameters

1506

1507 Parameters

1508 ==========

1509 keys: list, optional

1510 A list of keys for which return the credible levels, if None,

1511 defaults to search_parameter_keys

1512 weights: array, optional

1513 A list of weights for the posterior samples to calculate a set of

1514 weighted credible intervals.

1515 If None, assumes equal weights between samples.

1516

1517 Returns

1518 =======

1519 credible_levels: dict

1520 The credible levels at which the injected parameters are found.

1521 """

1522 if keys is None:

1523 keys = self.search_parameter_keys

1524 if self.injection_parameters is None:

1525 raise (

1526 TypeError,

1527 "Result object has no 'injection_parameters'. "

1528 "Cannot compute credible levels."

1529 )

1530 credible_levels = {key: self.get_injection_credible_level(key, weights=weights)

1531 for key in keys

1532 if isinstance(self.injection_parameters.get(key, None), float)}

1533 return credible_levels

1534

1535 def get_injection_credible_level(self, parameter, weights=None):

1536 """

1537 Get the credible level of the injected parameter

1538

1539 Calculated as CDF(injection value)

1540

1541 Parameters

1542 ==========

1543 parameter: str

1544 Parameter to get credible level for

1545 weights: array, optional

1546 A list of weights for the posterior samples to calculate a

1547 weighted credible interval.

1548 If None, assumes equal weights between samples.

1549

1550 Returns

1551 =======

1552 float: credible level

1553 """

1554 if self.injection_parameters is None:

1555 raise (

1556 TypeError,

1557 "Result object has no 'injection_parameters'. "

1558 "Cannot copmute credible levels."

1559 )

1560

1561 if weights is None:

1562 weights = np.ones(len(self.posterior))

1563

1564 if parameter in self.posterior and\

1565 parameter in self.injection_parameters:

1566 credible_level =\

1567 sum(np.array(self.posterior[parameter].values <

1568 self.injection_parameters[parameter]) * weights) / (sum(weights))

1569 return credible_level

1570 else:

1571 return np.nan

1572

1573 def _check_attribute_match_to_other_object(self, name, other_object):

1574 """ Check attribute name exists in other_object and is the same

1575

1576 Parameters

1577 ==========

1578 name: str

1579 Name of the attribute in this instance

1580 other_object: object

1581 Other object with attributes to compare with

1582

1583 Returns

1584 =======

1585 bool: True if attribute name matches with an attribute of other_object, False otherwise

1586

1587 """

1588 a = getattr(self, name, False)

1589 b = getattr(other_object, name, False)

1590 logger.debug('Checking {} value: {}=={}'.format(name, a, b))

1591 if (a is not False) and (b is not False):

1592 type_a = type(a)

1593 type_b = type(b)

1594 if type_a == type_b:

1595 if type_a in [str, float, int, dict, list]:

1596 try:

1597 return a == b

1598 except ValueError:

1599 return False

1600 elif type_a in [np.ndarray]:

1601 return np.all(a == b)

1602 return False

1603

1604 @property

1605 def kde(self):

1606 """ Kernel density estimate built from the stored posterior

1607

1608 Uses `scipy.stats.gaussian_kde` to generate the kernel density

1609 """

1610 if self._kde:

1611 return self._kde

1612 else:

1613 self._kde = scipy.stats.gaussian_kde(

1614 self.posterior[self.search_parameter_keys].values.T)

1615 return self._kde

1616

1617 def posterior_probability(self, sample):

1618 """ Calculate the posterior probability for a new sample

1619

1620 This queries a Kernel Density Estimate of the posterior to calculate

1621 the posterior probability density for the new sample.

1622

1623 Parameters

1624 ==========

1625 sample: dict, or list of dictionaries

1626 A dictionary containing all the keys from

1627 self.search_parameter_keys and corresponding values at which to

1628 calculate the posterior probability

1629

1630 Returns

1631 =======

1632 p: array-like,

1633 The posterior probability of the sample

1634

1635 """

1636 if isinstance(sample, dict):

1637 sample = [sample]

1638 ordered_sample = [[s[key] for key in self.search_parameter_keys]

1639 for s in sample]

1640 return self.kde(ordered_sample)

1641

1642 def _safe_outdir_creation(self, outdir=None, caller_func=None):

1643 if outdir is None:

1644 outdir = self.outdir

1645 try:

1646 utils.check_directory_exists_and_if_not_mkdir(outdir)

1647 except PermissionError:

1648 raise FileMovedError("Can not write in the out directory.\n"

1649 "Did you move the here file from another system?\n"

1650 "Try calling " + caller_func.__name__ + " with the 'outdir' "

1651 "keyword argument, e.g. " + caller_func.__name__ + "(outdir='.')")

1652 return outdir

1653

1654 def get_weights_by_new_prior(self, old_prior, new_prior, prior_names=None):

1655 """ Calculate a list of sample weights based on the ratio of new to old priors

1656

1657 Parameters

1658 ==========

1659 old_prior: PriorDict,

1660 The prior used in the generation of the original samples.

1661

1662 new_prior: PriorDict,

1663 The prior to use to reweight the samples.

1664

1665 prior_names: list

1666 A list of the priors to include in the ratio during reweighting.

1667

1668 Returns

1669 =======

1670 weights: array-like,

1671 A list of sample weights.

1672

1673 """

1674 weights = []

1675

1676 # Shared priors - these will form a ratio

1677 if prior_names is not None:

1678 shared_parameters = {key: self.posterior[key] for key in new_prior if

1679 key in old_prior and key in prior_names}

1680 else:

1681 shared_parameters = {key: self.posterior[key] for key in new_prior if key in old_prior}

1682 parameters = [{key: self.posterior[key][i] for key in shared_parameters.keys()}

1683 for i in range(len(self.posterior))]

1684

1685 for i in range(len(self.posterior)):

1686 weight = 1

1687 for prior_key in shared_parameters.keys():

1688 val = self.posterior[prior_key][i]

1689 weight *= new_prior.evaluate_constraints(parameters[i])

1690 weight *= new_prior[prior_key].prob(val) / old_prior[prior_key].prob(val)

1691

1692 weights.append(weight)

1693

1694 return weights

1695

1696 def to_arviz(self, prior=None):

1697 """ Convert the Result object to an ArviZ InferenceData object.

1698

1699 Parameters

1700 ==========

1701 prior: int

1702 If a positive integer is given then that number of prior

1703 samples will be drawn and stored in the ArviZ InferenceData

1704 object.

1705

1706 Returns

1707 =======

1708 azdata: InferenceData

1709 The ArviZ InferenceData object.

1710 """

1711

1712 try:

1713 import arviz as az

1714 except ImportError:

1715 logger.debug(

1716 "ArviZ is not installed, so cannot convert to InferenceData"

1717 )

1718

1719 posdict = {}

1720 for key in self.posterior:

1721 posdict[key] = self.posterior[key].values

1722

1723 if "log_likelihood" in posdict:

1724 loglikedict = {

1725 "log_likelihood": posdict.pop("log_likelihood")

1726 }

1727 else:

1728 if self.log_likelihood_evaluations is not None:

1729 loglikedict = {

1730 "log_likelihood": self.log_likelihood_evaluations

1731 }

1732 else:

1733 loglikedict = None

1734

1735 priorsamples = None

1736 if prior is not None:

1737 if self.priors is None:

1738 logger.warning(

1739 "No priors are in the Result object, so prior samples "

1740 "will not be included in the output."

1741 )

1742 else:

1743 priorsamples = self.priors.sample(size=prior)

1744

1745 azdata = az.from_dict(

1746 posterior=posdict,

1747 log_likelihood=loglikedict,

1748 prior=priorsamples,

1749 )

1750

1751 # add attributes

1752 version = {

1753 "inference_library": "bilby: {}".format(self.sampler),

1754 "inference_library_version": get_version_information()

1755 }

1756

1757 azdata.posterior.attrs.update(version)

1758 if "log_likelihood" in azdata._groups:

1759 azdata.log_likelihood.attrs.update(version)

1760 if "prior" in azdata._groups:

1761 azdata.prior.attrs.update(version)

1762

1763 return azdata

1764

1765

1766class ResultList(list):

1767

1768 def __init__(self, results=None, consistency_level="warning"):

1769 """ A class to store a list of :class:`bilby.core.result.Result` objects

1770 from equivalent runs on the same data. This provides methods for

1771 outputting combined results.

1772

1773 Parameters

1774 ==========

1775 results: list

1776 A list of `:class:`bilby.core.result.Result`.

1777 consistency_level: str, [ignore, warning, error]

1778 If warning, print a warning if inconsistencies are discovered

1779 between the results. If error, raise an error if inconsistencies

1780 are discovered between the results before combining. If ignore, do

1781 nothing.

1782

1783 """

1784 super(ResultList, self).__init__()

1785 self.consistency_level = consistency_level

1786 for result in results:

1787 self.append(result)

1788

1789 def append(self, result):

1790 """

1791 Append a :class:`bilby.core.result.Result`, or set of results, to the

1792 list.

1793

1794 Parameters

1795 ==========

1796 result: :class:`bilby.core.result.Result` or filename

1797 pointing to a result object, to append to the list.

1798 """

1799

1800 if isinstance(result, Result):

1801 super(ResultList, self).append(result)

1802 elif isinstance(result, (str, os.PathLike)):

1803 super(ResultList, self).append(read_in_result(result))

1804 else:

1805 raise TypeError("Could not append a non-Result type")

1806

1807 def combine(self, shuffle=False, consistency_level="error"):

1808 """

1809 Return the combined results in a :class:bilby.core.result.Result`

1810 object.

1811

1812 Parameters

1813 ----------

1814 shuffle: bool

1815 If true, shuffle the samples when combining, otherwise they are concatenated.

1816 consistency_level: str, [ignore, warning, error]

1817 Overwrite the class level consistency_level. If warning, print a

1818 warning if inconsistencies are discovered between the results. If

1819 error, raise an error if inconsistencies are discovered between

1820 the results before combining. If ignore, do nothing.

1821

1822 Returns

1823 -------

1824 result: bilby.core.result.Result

1825 The combined result file

1826

1827 """

1828

1829 self.consistency_level = consistency_level

1830

1831 if len(self) == 0:

1832 return Result()

1833 elif len(self) == 1:

1834 return copy(self[0])

1835 else:

1836 result = copy(self[0])

1837

1838 if result.label is not None:

1839 result.label += '_combined'

1840

1841 self.check_consistent_sampler()

1842 self.check_consistent_data()

1843 self.check_consistent_parameters()

1844 self.check_consistent_priors()

1845

1846 # check which kind of sampler was used: MCMC or Nested Sampling

1847 if result._nested_samples is not None:

1848 posteriors, result = self._combine_nested_sampled_runs(result)

1849 elif result.sampler in ["bilby_mcmc", "bilbymcmc"]:

1850 posteriors, result = self._combine_mcmc_sampled_runs(result)

1851 else:

1852 posteriors = [res.posterior for res in self]

1853

1854 combined_posteriors = pd.concat(posteriors, ignore_index=True)

1855

1856 if shuffle:

1857 result.posterior = combined_posteriors.sample(len(combined_posteriors))

1858 else:

1859 result.posterior = combined_posteriors

1860

1861 logger.info(f"Combined results have {len(result.posterior)} samples")

1862

1863 return result

1864

1865 def _combine_nested_sampled_runs(self, result):

1866 """

1867 Combine multiple nested sampling runs.

1868

1869 Currently this keeps posterior samples from each run in proportion with

1870 the evidence for each individual run

1871

1872 Parameters

1873 ==========

1874 result: bilby.core.result.Result

1875 The result object to put the new samples in.

1876

1877 Returns

1878 =======

1879 posteriors: list

1880 A list of pandas DataFrames containing the reduced sample set from

1881 each run.

1882 result: bilby.core.result.Result

1883 The result object with the combined evidences.

1884 """

1885 from scipy.special import logsumexp

1886 self.check_nested_samples()

1887

1888 # Combine evidences

1889 log_evidences = np.array([res.log_evidence for res in self])

1890 result.log_evidence = logsumexp(log_evidences, b=1. / len(self))

1891 result.log_bayes_factor = result.log_evidence - result.log_noise_evidence

1892

1893 # Propagate uncertainty in combined evidence

1894 log_errs = [res.log_evidence_err for res in self if np.isfinite(res.log_evidence_err)]

1895 if len(log_errs) > 0:

1896 result.log_evidence_err = 0.5 * logsumexp(2 * np.array(log_errs), b=1. / len(self))

1897 else:

1898 result.log_evidence_err = np.nan

1899

1900 # Combined posteriors with a weighting

1901 result_weights = np.exp(log_evidences - np.max(log_evidences))

1902 posteriors = list()

1903 for res, frac in zip(self, result_weights):

1904 selected_samples = (random.rng.uniform(size=len(res.posterior)) < frac)

1905 posteriors.append(res.posterior[selected_samples])

1906

1907 # remove original nested_samples

1908 result.nested_samples = None

1909 result.sampler_kwargs = None

1910 return posteriors, result

1911

1912 def _combine_mcmc_sampled_runs(self, result):

1913 """

1914 Combine multiple MCMC sampling runs.

1915

1916 Currently this keeps all posterior samples from each run

1917

1918 Parameters

1919 ----------

1920 result: bilby.core.result.Result

1921 The result object to put the new samples in.

1922

1923 Returns

1924 -------

1925 posteriors: list

1926 A list of pandas DataFrames containing the reduced sample set from

1927 each run.

1928 result: bilby.core.result.Result

1929 The result object with the combined evidences.

1930 """

1931

1932 from scipy.special import logsumexp

1933

1934 # Combine evidences

1935 log_evidences = np.array([res.log_evidence for res in self])

1936 result.log_evidence = logsumexp(log_evidences, b=1. / len(self))

1937 result.log_bayes_factor = result.log_evidence - result.log_noise_evidence

1938

1939 # Propagate uncertainty in combined evidence

1940 log_errs = [res.log_evidence_err for res in self if np.isfinite(res.log_evidence_err)]

1941 if len(log_errs) > 0:

1942 result.log_evidence_err = 0.5 * logsumexp(2 * np.array(log_errs), b=1. / len(self))

1943 else:

1944 result.log_evidence_err = np.nan

1945

1946 # Combined posteriors with a weighting

1947 posteriors = list()

1948 for res in self:

1949 posteriors.append(res.posterior)

1950

1951 return posteriors, result

1952

1953 def check_nested_samples(self):

1954 for res in self:

1955 try:

1956 res.nested_samples

1957 except ValueError:

1958 raise ResultListError("Not all results contain nested samples")

1959

1960 def _error_or_warning_consistency(self, msg):

1961 if self.consistency_level == "error":

1962 raise ResultListError(msg)

1963 elif self.consistency_level == "warning":

1964 logger.warning(msg)

1965 elif self.consistency_level == "ignore":

1966 pass

1967 else:

1968 raise ValueError(f"Input consistency_level {self.consistency_level} not understood")

1969

1970 def check_consistent_priors(self):

1971 for res in self:

1972 for p in self[0].priors.keys():

1973 if not self[0].priors[p] == res.priors[p] or len(self[0].priors) != len(res.priors):

1974 msg = "Inconsistent priors between results"

1975 self._error_or_warning_consistency(msg)

1976

1977 def check_consistent_parameters(self):

1978 if not np.all([set(self[0].search_parameter_keys) == set(res.search_parameter_keys) for res in self]):

1979 msg = "Inconsistent parameters between results"

1980 self._error_or_warning_consistency(msg)

1981

1982 def check_consistent_data(self):

1983 if not np.allclose(

1984 [res.log_noise_evidence for res in self],

1985 self[0].log_noise_evidence,

1986 atol=1e-8,

1987 rtol=0.0,

1988 equal_nan=True,

1989 ):

1990 msg = "Inconsistent data between results"

1991 self._error_or_warning_consistency(msg)

1992

1993 def check_consistent_sampler(self):

1994 if not np.all([res.sampler == self[0].sampler for res in self]):

1995 msg = "Inconsistent samplers between results"

1996 self._error_or_warning_consistency(msg)

1997

1998

1999@latex_plot_format

2000def plot_multiple(results, filename=None, labels=None, colours=None,

2001 save=True, evidences=False, corner_labels=None, linestyles=None,

2002 **kwargs):

2003 """ Generate a corner plot overlaying two sets of results

2004

2005 Parameters

2006 ==========

2007 results: list

2008 A list of `bilby.core.result.Result` objects containing the samples to

2009 plot.

2010 filename: str

2011 File name to save the figure to. If None (default), a filename is

2012 constructed from the outdir of the first element of results and then

2013 the labels for all the result files.

2014 labels: list

2015 List of strings to use when generating a legend. If None (default), the

2016 `label` attribute of each result in `results` is used.

2017 colours: list

2018 The colours for each result. If None, default styles are applied.

2019 save: bool

2020 If true, save the figure

2021 kwargs: dict

2022 All other keyword arguments are passed to `result.plot_corner` (except

2023 for the keyword `labels` for which you should use the dedicated

2024 `corner_labels` input).

2025 However, `show_titles` and `truths` are ignored since they would be

2026 ambiguous on such a plot.

2027 evidences: bool, optional

2028 Add the log-evidence calculations to the legend. If available, the

2029 Bayes factor will be used instead.

2030 corner_labels: list, optional

2031 List of strings to be passed to the input `labels` to `result.plot_corner`.

2032

2033 Returns

2034 =======

2035 fig:

2036 A matplotlib figure instance

2037

2038 """

2039 import matplotlib.pyplot as plt

2040 import matplotlib.lines as mpllines

2041

2042 kwargs['show_titles'] = False

2043 kwargs['truths'] = None

2044 if corner_labels is not None:

2045 kwargs['labels'] = corner_labels

2046

2047 fig = results[0].plot_corner(save=False, **kwargs)

2048 default_filename = '{}/{}'.format(results[0].outdir, 'combined')

2049 lines = []

2050 default_labels = []

2051 for i, result in enumerate(results):

2052 if colours:

2053 c = colours[i]

2054 else:

2055 c = 'C{}'.format(i)

2056 if linestyles is not None:

2057 linestyle = linestyles[i]

2058 else:

2059 linestyle = 'solid'

2060 hist_kwargs = kwargs.get('hist_kwargs', dict())

2061 hist_kwargs['color'] = c

2062 hist_kwargs["linestyle"] = linestyle

2063 kwargs["hist_kwargs"] = hist_kwargs

2064 fig = result.plot_corner(fig=fig, save=False, color=c, contour_kwargs={"linestyles": linestyle}, **kwargs)

2065 default_filename += '_{}'.format(result.label)

2066 lines.append(mpllines.Line2D([0], [0], color=c, linestyle=linestyle))

2067 default_labels.append(result.label)

2068

2069 # Rescale the axes

2070 for i, ax in enumerate(fig.axes):

2071 ax.autoscale()

2072 plt.draw()

2073

2074 if labels is None:

2075 labels = default_labels

2076

2077 labels = sanity_check_labels(labels)

2078

2079 if evidences:

2080 if np.isnan(results[0].log_bayes_factor):

2081 template = r'{label} $\mathrm{{ln}}(Z)={lnz:1.3g}$'

2082 else:

2083 template = r'{label} $\mathrm{{ln}}(B)={lnbf:1.3g}$'

2084 labels = [

2085 template.format(

2086 label=label,

2087 lnz=result.log_evidence,

2088 lnbf=result.log_bayes_factor,

2089 )

2090 for label, result in zip(labels, results)

2091 ]

2092

2093 axes = fig.get_axes()

2094 ndim = int(np.sqrt(len(axes)))

2095 axes[ndim - 1].legend(lines, labels)

2096

2097 if filename is None:

2098 filename = default_filename

2099

2100 if save:

2101 safe_save_figure(fig=fig, filename=filename)

2102 return fig

2103

2104

2105@latex_plot_format

2106def make_pp_plot(results, filename=None, save=True, confidence_interval=[0.68, 0.95, 0.997],

2107 lines=None, legend_fontsize='x-small', keys=None, title=True,

2108 confidence_interval_alpha=0.1, weight_list=None,

2109 **kwargs):

2110 """

2111 Make a P-P plot for a set of runs with injected signals.

2112

2113 Parameters

2114 ==========

2115 results: list

2116 A list of Result objects, each of these should have injected_parameters

2117 filename: str, optional

2118 The name of the file to save, the default is "outdir/pp.png"

2119 save: bool, optional

2120 Whether to save the file, default=True

2121 confidence_interval: (float, list), optional

2122 The confidence interval to be plotted, defaulting to 1-2-3 sigma

2123 lines: list

2124 If given, a list of matplotlib line formats to use, must be greater

2125 than the number of parameters.

2126 legend_fontsize: float

2127 The font size for the legend

2128 keys: list

2129 A list of keys to use, if None defaults to search_parameter_keys

2130 title: bool

2131 Whether to add the number of results and total p-value as a plot title

2132 confidence_interval_alpha: float, list, optional

2133 The transparency for the background condifence interval

2134 weight_list: list, optional

2135 List of the weight arrays for each set of posterior samples.

2136 kwargs:

2137 Additional kwargs to pass to matplotlib.pyplot.plot

2138

2139 Returns

2140 =======

2141 fig, pvals:

2142 matplotlib figure and a NamedTuple with attributes `combined_pvalue`,

2143 `pvalues`, and `names`.

2144 """

2145 import matplotlib.pyplot as plt

2146

2147 if keys is None:

2148 keys = results[0].search_parameter_keys

2149

2150 if weight_list is None:

2151 weight_list = [None] * len(results)

2152

2153 credible_levels = list()

2154 for i, result in enumerate(results):

2155 credible_levels.append(

2156 result.get_all_injection_credible_levels(keys, weights=weight_list[i])

2157 )

2158 credible_levels = pd.DataFrame(credible_levels)

2159

2160 if lines is None:

2161 colors = ["C{}".format(i) for i in range(8)]

2162 linestyles = ["-", "--", ":"]

2163 lines = ["{}{}".format(a, b) for a, b in product(linestyles, colors)]

2164 if len(lines) < len(credible_levels.keys()):

2165 raise ValueError("Larger number of parameters than unique linestyles")

2166

2167 x_values = np.linspace(0, 1, 1001)

2168

2169 N = len(credible_levels)

2170 fig, ax = plt.subplots()

2171

2172 if isinstance(confidence_interval, float):

2173 confidence_interval = [confidence_interval]

2174 if isinstance(confidence_interval_alpha, float):

2175 confidence_interval_alpha = [confidence_interval_alpha] * len(confidence_interval)

2176 elif len(confidence_interval_alpha) != len(confidence_interval):

2177 raise ValueError(

2178 "confidence_interval_alpha must have the same length as confidence_interval")

2179

2180 for ci, alpha in zip(confidence_interval, confidence_interval_alpha):

2181 edge_of_bound = (1. - ci) / 2.

2182 lower = scipy.stats.binom.ppf(1 - edge_of_bound, N, x_values) / N

2183 upper = scipy.stats.binom.ppf(edge_of_bound, N, x_values) / N

2184 # The binomial point percent function doesn't always return 0 @ 0,

2185 # so set those bounds explicitly to be sure

2186 lower[0] = 0

2187 upper[0] = 0

2188 ax.fill_between(x_values, lower, upper, alpha=alpha, color='k')

2189

2190 pvalues = []

2191 logger.info("Key: KS-test p-value")

2192 for ii, key in enumerate(credible_levels):

2193 pp = np.array([sum(credible_levels[key].values < xx) /

2194 len(credible_levels) for xx in x_values])

2195 pvalue = scipy.stats.kstest(credible_levels[key], 'uniform').pvalue

2196 pvalues.append(pvalue)

2197 logger.info("{}: {}".format(key, pvalue))

2198

2199 try:

2200 name = results[0].priors[key].latex_label

2201 except (AttributeError, KeyError):

2202 name = key

2203 label = "{} ({:2.3f})".format(name, pvalue)

2204 plt.plot(x_values, pp, lines[ii], label=label, **kwargs)

2205

2206 Pvals = namedtuple('pvals', ['combined_pvalue', 'pvalues', 'names'])

2207 pvals = Pvals(combined_pvalue=scipy.stats.combine_pvalues(pvalues)[1],

2208 pvalues=pvalues,

2209 names=list(credible_levels.keys()))

2210 logger.info(

2211 "Combined p-value: {}".format(pvals.combined_pvalue))

2212

2213 if title:

2214 ax.set_title("N={}, p-value={:2.4f}".format(

2215 len(results), pvals.combined_pvalue))

2216 ax.set_xlabel("C.I.")

2217 ax.set_ylabel("Fraction of events in C.I.")

2218 ax.legend(handlelength=2, labelspacing=0.25, fontsize=legend_fontsize)

2219 ax.set_xlim(0, 1)

2220 ax.set_ylim(0, 1)

2221 fig.tight_layout()

2222 if save:

2223 if filename is None:

2224 filename = 'outdir/pp.png'

2225 safe_save_figure(fig=fig, filename=filename, dpi=500)

2226

2227 return fig, pvals

2228

2229

2230def sanity_check_labels(labels):

2231 """ Check labels for plotting to remove matplotlib errors """

2232 for ii, lab in enumerate(labels):

2233 if "_" in lab and "$" not in lab:

2234 labels[ii] = lab.replace("_", "-")

2235 return labels

2236

2237

2238class ResultError(Exception):

2239 """ Base exception for all Result related errors """

2240

2241

2242class ResultListError(ResultError):

2243 """ For Errors occurring during combining results. """

2244

2245

2246class FileMovedError(ResultError):

2247 """ Exceptions that occur when files have been moved """