Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • john-veitch/bilby
  • duncanmmacleod/bilby
  • colm.talbot/bilby
  • lscsoft/bilby
  • matthew-pitkin/bilby
  • salvatore-vitale/tupak
  • charlie.hoy/bilby
  • bfarr/bilby
  • virginia.demilio/bilby
  • vivien/bilby
  • eric-howell/bilby
  • sebastian-khan/bilby
  • rhys.green/bilby
  • moritz.huebner/bilby
  • joseph.mills/bilby
  • scott.coughlin/bilby
  • matthew.carney/bilby
  • hyungwon.lee/bilby
  • monica.rizzo/bilby
  • christopher-berry/bilby
  • lindsay.demarchi/bilby
  • kaushik.rao/bilby
  • charles.kimball/bilby
  • andrew.matas/bilby
  • juan.calderonbustillo/bilby
  • patrick-meyers/bilby
  • hannah.middleton/bilby
  • eve.chase/bilby
  • grant.meadors/bilby
  • khun.phukon/bilby
  • sumeet.kulkarni/bilby
  • daniel.reardon/bilby
  • cjhaster/bilby
  • sylvia.biscoveanu/bilby
  • james-clark/bilby
  • meg.millhouse/bilby
  • joshua.willis/bilby
  • nikhil.sarin/bilby
  • paul.easter/bilby
  • youngmin/bilby
  • daniel-williams/bilby
  • shanika.galaudage/bilby
  • bruce.edelman/bilby
  • avi.vajpeyi/bilby
  • isobel.romero-shaw/bilby
  • andrew.kim/bilby
  • dominika.zieba/bilby
  • jonathan.davies/bilby
  • marc.arene/bilby
  • srishti.tiwari/bilby-tidal-heating-eccentric
  • aditya.vijaykumar/bilby
  • michael.williams/bilby
  • cecilio.garcia-quiros/bilby
  • rory-smith/bilby
  • maite.mateu-lucena/bilby
  • wushichao/bilby
  • kaylee.desoto/bilby
  • brandon.piotrzkowski/bilby
  • rossella.gamba/bilby
  • hunter.gabbard/bilby
  • deep.chatterjee/bilby
  • tathagata.ghosh/bilby
  • arunava.mukherjee/bilby
  • philip.relton/bilby
  • reed.essick/bilby
  • pawan.gupta/bilby
  • francisco.hernandez/bilby
  • rhiannon.udall/bilby
  • leo.tsukada/bilby
  • will-farr/bilby
  • vijay.varma/bilby
  • jeremy.baier/bilby
  • joshua.brandt/bilby
  • ethan.payne/bilby
  • ka-lok.lo/bilby
  • antoni.ramos-buades/bilby
  • oliviastephany.wilk/bilby
  • jack.heinzel/bilby
  • samson.leong/bilby-psi4
  • viviana.caceres/bilby
  • nadia.qutob/bilby
  • michael-coughlin/bilby
  • hemantakumar.phurailatpam/bilby
  • boris.goncharov/bilby
  • sama.al-shammari/bilby
  • siqi.zhong/bilby
  • jocelyn-read/bilby
  • marc.penuliar/bilby
  • stephanie.letourneau/bilby
  • alexandresebastien.goettel/bilby
  • alec.gunny/bilby
  • serguei.ossokine/bilby
  • pratyusava.baral/bilby
  • sophie.hourihane/bilby
  • eunsub/bilby
  • james.hart/bilby
  • pratyusava.baral/bilby-tg
  • zhaozc/bilby
  • pratyusava.baral/bilby_SoG
  • tomasz.baka/bilby
  • nicogerardo.bers/bilby
  • soumen.roy/bilby
  • isaac.mcmahon/healpix-redundancy
  • asamakai.baker/bilby-frequency-dependent-antenna-pattern-functions
  • anna.puecher/bilby
  • pratyusava.baral/bilby-x-g
  • thibeau.wouters/bilby
  • christian.adamcewicz/bilby
  • raffi.enficiaud/bilby
109 results
Show changes
Commits on Source (18)
Showing
with 375 additions and 142 deletions
...@@ -13,6 +13,6 @@ MANIFEST ...@@ -13,6 +13,6 @@ MANIFEST
*.dat *.dat
*.version *.version
*.ipynb_checkpoints *.ipynb_checkpoints
outdir/* **/outdir
.idea/* .idea/*
bilby/_version.py bilby/_version.py
# All notable changes will be documented in this file # All notable changes will be documented in this file
## [1.4.1] 2022-12-07
Version 1.4.1 release of Bilby
This is a bugfix release to address some minor issues identified after v1.4.0.
### Changes
- Documentation updates (!1181, !1183)
- Fix some of the examples in the repository (!1182)
- Make sure conversion to symmetric mass ratio always returns a valid value (!1184)
- Provide a default nlive for dynamic dynesty (!1185)
- Enable the relative binning likelihood to be initialized with ra/dec when sampling in a different sky parameterization (!1186)
- Make sure that all dumping pickle files is done safely (!1189)
- Make error catching for `dynesty` checkpointing more robust (!1190)
## [1.4.0] 2022-11-18 ## [1.4.0] 2022-11-18
Version 1.4.0 release of Bilby Version 1.4.0 release of Bilby
......
...@@ -22,6 +22,7 @@ from .utils import ( ...@@ -22,6 +22,7 @@ from .utils import (
recursively_save_dict_contents_to_group, recursively_save_dict_contents_to_group,
recursively_load_dict_contents_from_group, recursively_load_dict_contents_from_group,
recursively_decode_bilby_json, recursively_decode_bilby_json,
safe_file_dump,
) )
from .prior import Prior, PriorDict, DeltaFunction, ConditionalDeltaFunction from .prior import Prior, PriorDict, DeltaFunction, ConditionalDeltaFunction
...@@ -735,16 +736,12 @@ class Result(object): ...@@ -735,16 +736,12 @@ class Result(object):
with h5py.File(filename, 'w') as h5file: with h5py.File(filename, 'w') as h5file:
recursively_save_dict_contents_to_group(h5file, '/', dictionary) recursively_save_dict_contents_to_group(h5file, '/', dictionary)
elif extension == 'pkl': elif extension == 'pkl':
import dill safe_file_dump(self, filename, "dill")
with open(filename, "wb") as ff:
dill.dump(self, ff)
else: else:
raise ValueError("Extension type {} not understood".format(extension)) raise ValueError("Extension type {} not understood".format(extension))
except Exception as e: except Exception as e:
import dill
filename = ".".join(filename.split(".")[:-1]) + ".pkl" filename = ".".join(filename.split(".")[:-1]) + ".pkl"
with open(filename, "wb") as ff: safe_file_dump(self, filename, "dill")
dill.dump(self, ff)
logger.error( logger.error(
"\n\nSaving the data has failed with the following message:\n" "\n\nSaving the data has failed with the following message:\n"
"{}\nData has been dumped to {}.\n\n".format(e, filename) "{}\nData has been dumped to {}.\n\n".format(e, filename)
......
...@@ -17,7 +17,16 @@ class DynamicDynesty(Dynesty): ...@@ -17,7 +17,16 @@ class DynamicDynesty(Dynesty):
@property @property
def nlive(self): def nlive(self):
return self.kwargs["nlive_init"] """
Users can either specify :code:`nlive_init` or :code:`nlive` (with
that precedence) or specify no value, in which case 500 is used.
"""
if self.kwargs["nlive_init"] is not None:
return self.kwargs["nlive_init"]
elif self.kwargs["nlive"] is not None:
return self.kwargs["nlive"]
else:
return 500
@property @property
def sampler_init(self): def sampler_init(self):
......
...@@ -343,14 +343,11 @@ class Dynesty(NestedSampler): ...@@ -343,14 +343,11 @@ class Dynesty(NestedSampler):
self.kwargs[key] = selected self.kwargs[key] = selected
def nestcheck_data(self, out_file): def nestcheck_data(self, out_file):
import pickle
import nestcheck.data_processing import nestcheck.data_processing
ns_run = nestcheck.data_processing.process_dynesty_run(out_file) ns_run = nestcheck.data_processing.process_dynesty_run(out_file)
nestcheck_result = f"{self.outdir}/{self.label}_nestcheck.pickle" nestcheck_result = f"{self.outdir}/{self.label}_nestcheck.pickle"
with open(nestcheck_result, "wb") as file_nest: safe_file_dump(ns_run, nestcheck_result, "pickle")
pickle.dump(ns_run, file_nest)
@property @property
def nlive(self): def nlive(self):
...@@ -370,7 +367,6 @@ class Dynesty(NestedSampler): ...@@ -370,7 +367,6 @@ class Dynesty(NestedSampler):
@signal_wrapper @signal_wrapper
def run_sampler(self): def run_sampler(self):
import dill
import dynesty import dynesty
logger.info(f"Using dynesty version {dynesty.__version__}") logger.info(f"Using dynesty version {dynesty.__version__}")
...@@ -436,8 +432,7 @@ class Dynesty(NestedSampler): ...@@ -436,8 +432,7 @@ class Dynesty(NestedSampler):
self.nestcheck_data(out) self.nestcheck_data(out)
dynesty_result = f"{self.outdir}/{self.label}_dynesty.pickle" dynesty_result = f"{self.outdir}/{self.label}_dynesty.pickle"
with open(dynesty_result, "wb") as file: safe_file_dump(out, dynesty_result, "dill")
dill.dump(out, file)
self._generate_result(out) self._generate_result(out)
self.result.sampling_time = self.sampling_time self.result.sampling_time = self.sampling_time
...@@ -669,10 +664,14 @@ class Dynesty(NestedSampler): ...@@ -669,10 +664,14 @@ class Dynesty(NestedSampler):
np.linalg.linalg.LinAlgError, np.linalg.linalg.LinAlgError,
ValueError, ValueError,
OverflowError, OverflowError,
Exception,
) as e: ) as e:
logger.warning(e) logger.warning(e)
logger.warning("Failed to create dynesty state plot at checkpoint") logger.warning("Failed to create dynesty state plot at checkpoint")
except Exception as e:
logger.warning(
f"Unexpected error {e} in dynesty plotting. "
"Please report at git.ligo.org/lscsoft/bilby/-/issues"
)
finally: finally:
plt.close("all") plt.close("all")
try: try:
...@@ -691,41 +690,60 @@ class Dynesty(NestedSampler): ...@@ -691,41 +690,60 @@ class Dynesty(NestedSampler):
np.linalg.linalg.LinAlgError, np.linalg.linalg.LinAlgError,
ValueError, ValueError,
OverflowError, OverflowError,
Exception,
) as e: ) as e:
logger.warning(e) logger.warning(e)
logger.warning("Failed to create dynesty unit state plot at checkpoint") logger.warning("Failed to create dynesty unit state plot at checkpoint")
except Exception as e:
logger.warning(
f"Unexpected error {e} in dynesty plotting. "
"Please report at git.ligo.org/lscsoft/bilby/-/issues"
)
finally: finally:
plt.close("all") plt.close("all")
try: try:
filename = f"{self.outdir}/{self.label}_checkpoint_run.png" filename = f"{self.outdir}/{self.label}_checkpoint_run.png"
fig, axs = dyplot.runplot( fig, _ = dyplot.runplot(
self.sampler.results, logplot=False, use_math_text=False self.sampler.results, logplot=False, use_math_text=False
) )
fig.tight_layout() fig.tight_layout()
plt.savefig(filename) plt.savefig(filename)
except (RuntimeError, np.linalg.linalg.LinAlgError, ValueError) as e: except (
RuntimeError,
np.linalg.linalg.LinAlgError,
ValueError,
OverflowError,
) as e:
logger.warning(e) logger.warning(e)
logger.warning("Failed to create dynesty run plot at checkpoint") logger.warning("Failed to create dynesty run plot at checkpoint")
except Exception as e:
logger.warning(
f"Unexpected error {e} in dynesty plotting. "
"Please report at git.ligo.org/lscsoft/bilby/-/issues"
)
finally: finally:
plt.close("all") plt.close("all")
try: try:
filename = f"{self.outdir}/{self.label}_checkpoint_stats.png" filename = f"{self.outdir}/{self.label}_checkpoint_stats.png"
fig, axs = dynesty_stats_plot(self.sampler) fig, _ = dynesty_stats_plot(self.sampler)
fig.tight_layout() fig.tight_layout()
plt.savefig(filename) plt.savefig(filename)
except (RuntimeError, ValueError) as e: except (RuntimeError, ValueError, OverflowError) as e:
logger.warning(e) logger.warning(e)
logger.warning("Failed to create dynesty stats plot at checkpoint") logger.warning("Failed to create dynesty stats plot at checkpoint")
except DynestySetupError: except DynestySetupError:
logger.debug("Cannot create Dynesty stats plot with dynamic sampler.") logger.debug("Cannot create Dynesty stats plot with dynamic sampler.")
except Exception as e:
logger.warning(
f"Unexpected error {e} in dynesty plotting. "
"Please report at git.ligo.org/lscsoft/bilby/-/issues"
)
finally: finally:
plt.close("all") plt.close("all")
def _run_test(self): def _run_test(self):
import pandas as pd import pandas as pd
self.sampler = self.sampler_class( self.sampler = self.sampler_init(
loglikelihood=self.log_likelihood, loglikelihood=self.log_likelihood,
prior_transform=self.prior_transform, prior_transform=self.prior_transform,
ndim=self.ndim, ndim=self.ndim,
...@@ -734,7 +752,14 @@ class Dynesty(NestedSampler): ...@@ -734,7 +752,14 @@ class Dynesty(NestedSampler):
sampler_kwargs = self.sampler_function_kwargs.copy() sampler_kwargs = self.sampler_function_kwargs.copy()
sampler_kwargs["maxiter"] = 2 sampler_kwargs["maxiter"] = 2
if self.print_method == "tqdm" and self.kwargs["print_progress"]:
from tqdm.auto import tqdm
self.pbar = tqdm(file=sys.stdout, initial=self.sampler.it)
self.sampler.run_nested(**sampler_kwargs) self.sampler.run_nested(**sampler_kwargs)
if self.pbar is not None:
self.pbar = self.pbar.close()
print("")
N = 100 N = 100
self.result.samples = pd.DataFrame(self.priors.sample(N))[ self.result.samples = pd.DataFrame(self.priors.sample(N))[
self.search_parameter_keys self.search_parameter_keys
......
...@@ -7,7 +7,7 @@ from shutil import copyfile ...@@ -7,7 +7,7 @@ from shutil import copyfile
import numpy as np import numpy as np
from pandas import DataFrame from pandas import DataFrame
from ..utils import check_directory_exists_and_if_not_mkdir, logger from ..utils import check_directory_exists_and_if_not_mkdir, logger, safe_file_dump
from .base_sampler import MCMCSampler, SamplerError, signal_wrapper from .base_sampler import MCMCSampler, SamplerError, signal_wrapper
from .ptemcee import LikePriorEvaluator from .ptemcee import LikePriorEvaluator
...@@ -285,20 +285,20 @@ class Emcee(MCMCSampler): ...@@ -285,20 +285,20 @@ class Emcee(MCMCSampler):
return self.sampler.chain[:, :nsteps, :] return self.sampler.chain[:, :nsteps, :]
def write_current_state(self): def write_current_state(self):
"""Writes a pickle file of the sampler to disk using dill""" """
import dill Writes a pickle file of the sampler to disk using dill
Overwrites the stored sampler chain with one that is truncated
to only the completed steps
"""
logger.info( logger.info(
f"Checkpointing sampler to file {self.checkpoint_info.sampler_file}" f"Checkpointing sampler to file {self.checkpoint_info.sampler_file}"
) )
with open(self.checkpoint_info.sampler_file, "wb") as f: self.sampler._chain = self.sampler_chain
# Overwrites the stored sampler chain with one that is truncated _pool = self.sampler.pool
# to only the completed steps self.sampler.pool = None
self.sampler._chain = self.sampler_chain safe_file_dump(self._sampler, self.checkpoint_info.sampler_file, "dill")
_pool = self.sampler.pool self.sampler.pool = _pool
self.sampler.pool = None
dill.dump(self._sampler, f)
self.sampler.pool = _pool
def _initialise_sampler(self): def _initialise_sampler(self):
from emcee import EnsembleSampler from emcee import EnsembleSampler
......
...@@ -8,7 +8,7 @@ from collections import namedtuple ...@@ -8,7 +8,7 @@ from collections import namedtuple
import numpy as np import numpy as np
import pandas as pd import pandas as pd
from ..utils import check_directory_exists_and_if_not_mkdir, logger from ..utils import check_directory_exists_and_if_not_mkdir, logger, safe_file_dump
from .base_sampler import ( from .base_sampler import (
MCMCSampler, MCMCSampler,
SamplerError, SamplerError,
...@@ -1189,8 +1189,6 @@ def checkpoint( ...@@ -1189,8 +1189,6 @@ def checkpoint(
Q_list, Q_list,
time_per_check, time_per_check,
): ):
import dill
logger.info("Writing checkpoint and diagnostics") logger.info("Writing checkpoint and diagnostics")
ndim = sampler.dim ndim = sampler.dim
...@@ -1223,8 +1221,7 @@ def checkpoint( ...@@ -1223,8 +1221,7 @@ def checkpoint(
pos0=pos0, pos0=pos0,
) )
with open(resume_file, "wb") as file: safe_file_dump(data, resume_file, "dill")
dill.dump(data, file, protocol=4)
del data, sampler_copy del data, sampler_copy
logger.info("Finished writing checkpoint") logger.info("Finished writing checkpoint")
......
...@@ -370,10 +370,11 @@ def safe_file_dump(data, filename, module): ...@@ -370,10 +370,11 @@ def safe_file_dump(data, filename, module):
data to dump data to dump
filename: str filename: str
The file to dump to The file to dump to
module: pickle, dill module: pickle, dill, str
The python module to use The python module to use. If a string, the module will be imported
""" """
if isinstance(module, str):
module = import_module(module)
temp_filename = filename + ".temp" temp_filename = filename + ".temp"
with open(temp_filename, "wb") as file: with open(temp_filename, "wb") as file:
module.dump(data, file) module.dump(data, file)
......
...@@ -12,7 +12,7 @@ import numpy as np ...@@ -12,7 +12,7 @@ import numpy as np
from pandas import DataFrame, Series from pandas import DataFrame, Series
from ..core.likelihood import MarginalizedLikelihoodReconstructionError from ..core.likelihood import MarginalizedLikelihoodReconstructionError
from ..core.utils import logger, solar_mass, command_line_args from ..core.utils import logger, solar_mass, command_line_args, safe_file_dump
from ..core.prior import DeltaFunction from ..core.prior import DeltaFunction
from .utils import lalsim_SimInspiralTransformPrecessingNewInitialConditions from .utils import lalsim_SimInspiralTransformPrecessingNewInitialConditions
from .eos.eos import SpectralDecompositionEOS, EOSFamily, IntegrateTOV from .eos.eos import SpectralDecompositionEOS, EOSFamily, IntegrateTOV
...@@ -597,7 +597,7 @@ def component_masses_to_symmetric_mass_ratio(mass_1, mass_2): ...@@ -597,7 +597,7 @@ def component_masses_to_symmetric_mass_ratio(mass_1, mass_2):
Symmetric mass ratio of the binary Symmetric mass ratio of the binary
""" """
return (mass_1 * mass_2) / (mass_1 + mass_2) ** 2 return np.minimum((mass_1 * mass_2) / (mass_1 + mass_2) ** 2, 1 / 4)
def component_masses_to_mass_ratio(mass_1, mass_2): def component_masses_to_mass_ratio(mass_1, mass_2):
...@@ -1687,8 +1687,7 @@ def generate_posterior_samples_from_marginalized_likelihood( ...@@ -1687,8 +1687,7 @@ def generate_posterior_samples_from_marginalized_likelihood(
cached_samples_dict[ii] = subset_samples cached_samples_dict[ii] = subset_samples
if use_cache: if use_cache:
with open(cache_filename, "wb") as f: safe_file_dump(cached_samples_dict, cache_filename, "pickle")
pickle.dump(cached_samples_dict, f)
ii += block ii += block
pbar.update(len(subset_samples)) pbar.update(len(subset_samples))
......
...@@ -8,7 +8,7 @@ from bilby_cython.geometry import ( ...@@ -8,7 +8,7 @@ from bilby_cython.geometry import (
) )
from ...core import utils from ...core import utils
from ...core.utils import docstring, logger, PropertyAccessor from ...core.utils import docstring, logger, PropertyAccessor, safe_file_dump
from .. import utils as gwutils from .. import utils as gwutils
from .calibration import Recalibrate from .calibration import Recalibrate
from .geometry import InterferometerGeometry from .geometry import InterferometerGeometry
...@@ -798,11 +798,9 @@ class Interferometer(object): ...@@ -798,11 +798,9 @@ class Interferometer(object):
format="pickle", extra=".. versionadded:: 1.1.0" format="pickle", extra=".. versionadded:: 1.1.0"
)) ))
def to_pickle(self, outdir="outdir", label=None): def to_pickle(self, outdir="outdir", label=None):
import dill
utils.check_directory_exists_and_if_not_mkdir('outdir') utils.check_directory_exists_and_if_not_mkdir('outdir')
filename = self._filename_from_outdir_label_extension(outdir, label, extension="pkl") filename = self._filename_from_outdir_label_extension(outdir, label, extension="pkl")
with open(filename, "wb") as ff: safe_file_dump(self, filename, "dill")
dill.dump(self, ff)
@classmethod @classmethod
@docstring(_load_docstring.format(format="pickle")) @docstring(_load_docstring.format(format="pickle"))
......
...@@ -4,7 +4,7 @@ import numpy as np ...@@ -4,7 +4,7 @@ import numpy as np
import math import math
from ...core import utils from ...core import utils
from ...core.utils import logger from ...core.utils import logger, safe_file_dump
from .interferometer import Interferometer from .interferometer import Interferometer
from .psd import PowerSpectralDensity from .psd import PowerSpectralDensity
...@@ -273,15 +273,12 @@ class InterferometerList(list): ...@@ -273,15 +273,12 @@ class InterferometerList(list):
""" """
def to_pickle(self, outdir="outdir", label="ifo_list"): def to_pickle(self, outdir="outdir", label="ifo_list"):
import dill
utils.check_directory_exists_and_if_not_mkdir("outdir") utils.check_directory_exists_and_if_not_mkdir("outdir")
label = label + "_" + "".join(ifo.name for ifo in self) label = label + "_" + "".join(ifo.name for ifo in self)
filename = self._filename_from_outdir_label_extension( filename = self._filename_from_outdir_label_extension(
outdir, label, extension="pkl" outdir, label, extension="pkl"
) )
with open(filename, "wb") as ff: safe_file_dump(self, filename, "dill")
dill.dump(self, ff)
@classmethod @classmethod
def from_pickle(cls, filename=None): def from_pickle(cls, filename=None):
......
...@@ -1112,21 +1112,55 @@ class GravitationalWaveTransient(Likelihood): ...@@ -1112,21 +1112,55 @@ class GravitationalWaveTransient(Likelihood):
else: else:
raise ValueError("Unable to parse reference frame {}".format(frame)) raise ValueError("Unable to parse reference frame {}".format(frame))
def get_sky_frame_parameters(self): def get_sky_frame_parameters(self, parameters=None):
time = self.parameters['{}_time'.format(self.time_reference)] """
Generate ra, dec, and geocenter time for :code:`parameters`
This method will attempt to convert from the reference time and sky
parameters, but if they are not present it will fall back to ra and dec.
Parameters
==========
parameters: dict, optional
The parameters to be converted.
If not specified :code:`self.parameters` will be used.
Returns
=======
dict: dictionary containing ra, dec, and geocent_time
"""
if parameters is None:
parameters = self.parameters
time = parameters.get(f'{self.time_reference}_time', None)
if time is None and "geocent_time" in parameters:
logger.warning(
f"Cannot find {self.time_reference}_time in parameters. "
"Falling back to geocent time"
)
if not self.reference_frame == "sky": if not self.reference_frame == "sky":
ra, dec = zenith_azimuth_to_ra_dec( try:
self.parameters['zenith'], self.parameters['azimuth'], ra, dec = zenith_azimuth_to_ra_dec(
time, self.reference_frame) parameters['zenith'], parameters['azimuth'],
time, self.reference_frame)
except KeyError:
if "ra" in parameters and "dec" in parameters:
ra = parameters["ra"]
dec = parameters["dec"]
logger.warning(
"Cannot convert from zenith/azimuth to ra/dec falling "
"back to provided ra/dec"
)
else:
raise
else: else:
ra = self.parameters["ra"] ra = parameters["ra"]
dec = self.parameters["dec"] dec = parameters["dec"]
if "geocent" not in self.time_reference: if "geocent" not in self.time_reference and f"{self.time_reference}_time" in parameters:
geocent_time = time - self.reference_ifo.time_delay_from_geocenter( geocent_time = time - self.reference_ifo.time_delay_from_geocenter(
ra=ra, dec=dec, time=time ra=ra, dec=dec, time=time
) )
else: else:
geocent_time = self.parameters["geocent_time"] geocent_time = parameters["geocent_time"]
return dict(ra=ra, dec=dec, geocent_time=geocent_time) return dict(ra=ra, dec=dec, geocent_time=geocent_time)
@property @property
......
...@@ -23,7 +23,10 @@ class RelativeBinningGravitationalWaveTransient(GravitationalWaveTransient): ...@@ -23,7 +23,10 @@ class RelativeBinningGravitationalWaveTransient(GravitationalWaveTransient):
given some set of parameters given some set of parameters
fiducial_parameters: dict, optional fiducial_parameters: dict, optional
A starting guess for initial parameters of the event for finding the A starting guess for initial parameters of the event for finding the
maximum likelihood (fiducial) waveform. maximum likelihood (fiducial) waveform. These should be specified in
the same parameter basis as the one that sampling is carried out in.
For example, if sampling in `mass_1` and `mass_2`, the fiducial
parameters should also be provided in `mass_1` and `mass_2.`
parameter_bounds: dict, optional parameter_bounds: dict, optional
Dictionary of bounds (lists) for the initial parameters when finding Dictionary of bounds (lists) for the initial parameters when finding
the initial maximum likelihood (fiducial) waveform. the initial maximum likelihood (fiducial) waveform.
...@@ -216,8 +219,9 @@ class RelativeBinningGravitationalWaveTransient(GravitationalWaveTransient): ...@@ -216,8 +219,9 @@ class RelativeBinningGravitationalWaveTransient(GravitationalWaveTransient):
) )
def set_fiducial_waveforms(self, parameters): def set_fiducial_waveforms(self, parameters):
_fiducial = parameters["fiducial"] parameters = parameters.copy()
parameters["fiducial"] = 1 parameters["fiducial"] = 1
parameters.update(self.get_sky_frame_parameters(parameters=parameters))
self.fiducial_polarizations = self.waveform_generator.frequency_domain_strain( self.fiducial_polarizations = self.waveform_generator.frequency_domain_strain(
parameters) parameters)
...@@ -236,8 +240,6 @@ class RelativeBinningGravitationalWaveTransient(GravitationalWaveTransient): ...@@ -236,8 +240,6 @@ class RelativeBinningGravitationalWaveTransient(GravitationalWaveTransient):
wf[interferometer.frequency_array > self.maximum_frequency] = 0 wf[interferometer.frequency_array > self.maximum_frequency] = 0
self.per_detector_fiducial_waveforms[interferometer.name] = wf self.per_detector_fiducial_waveforms[interferometer.name] = wf
parameters["fiducial"] = _fiducial
def find_maximum_likelihood_parameters(self, parameter_bounds, def find_maximum_likelihood_parameters(self, parameter_bounds,
iterations=5, maximization_kwargs=None): iterations=5, maximization_kwargs=None):
if maximization_kwargs is None: if maximization_kwargs is None:
......
...@@ -7,7 +7,7 @@ import numpy as np ...@@ -7,7 +7,7 @@ import numpy as np
from ..core.result import Result as CoreResult from ..core.result import Result as CoreResult
from ..core.utils import ( from ..core.utils import (
infft, logger, check_directory_exists_and_if_not_mkdir, infft, logger, check_directory_exists_and_if_not_mkdir,
latex_plot_format, safe_save_figure latex_plot_format, safe_file_dump, safe_save_figure,
) )
from .utils import plot_spline_pos, spline_angle_xform, asd_from_freq_series from .utils import plot_spline_pos, spline_angle_xform, asd_from_freq_series
from .detector import get_empty_interferometer, Interferometer from .detector import get_empty_interferometer, Interferometer
...@@ -781,8 +781,7 @@ class CompactBinaryCoalescenceResult(CoreResult): ...@@ -781,8 +781,7 @@ class CompactBinaryCoalescenceResult(CoreResult):
logger.info('Initialising skymap class') logger.info('Initialising skymap class')
skypost = confidence_levels(pts, trials=trials, jobs=jobs) skypost = confidence_levels(pts, trials=trials, jobs=jobs)
logger.info('Pickling skymap to {}'.format(default_obj_filename)) logger.info('Pickling skymap to {}'.format(default_obj_filename))
with open(default_obj_filename, 'wb') as out: safe_file_dump(skypost, default_obj_filename, "pickle")
pickle.dump(skypost, out)
else: else:
if isinstance(load_pickle, str): if isinstance(load_pickle, str):
......
...@@ -42,7 +42,7 @@ def lal_binary_black_hole( ...@@ -42,7 +42,7 @@ def lal_binary_black_hole(
theta_jn: float theta_jn: float
Angle between the total binary angular momentum and the line of sight Angle between the total binary angular momentum and the line of sight
phase: float phase: float
The phase at coalescence The phase at reference frequency or peak amplitude (depends on waveform)
kwargs: dict kwargs: dict
Optional keyword arguments Optional keyword arguments
Supported arguments: Supported arguments:
...@@ -125,7 +125,7 @@ def lal_binary_neutron_star( ...@@ -125,7 +125,7 @@ def lal_binary_neutron_star(
theta_jn: float theta_jn: float
Orbital inclination Orbital inclination
phase: float phase: float
The phase at coalescence The phase at reference frequency or peak amplitude (depends on waveform)
lambda_1: float lambda_1: float
Dimensionless tidal deformability of mass_1 Dimensionless tidal deformability of mass_1
lambda_2: float lambda_2: float
...@@ -197,7 +197,7 @@ def lal_eccentric_binary_black_hole_no_spins( ...@@ -197,7 +197,7 @@ def lal_eccentric_binary_black_hole_no_spins(
theta_jn: float theta_jn: float
Orbital inclination Orbital inclination
phase: float phase: float
The phase at coalescence The phase at reference frequency or peak amplitude (depends on waveform)
kwargs: dict kwargs: dict
Optional keyword arguments Optional keyword arguments
Supported arguments: Supported arguments:
...@@ -275,7 +275,7 @@ def _base_lal_cbc_fd_waveform( ...@@ -275,7 +275,7 @@ def _base_lal_cbc_fd_waveform(
theta_jn: float theta_jn: float
Orbital inclination Orbital inclination
phase: float phase: float
The phase at coalescence The phase at reference frequency or peak amplitude (depends on waveform)
eccentricity: float eccentricity: float
Binary eccentricity Binary eccentricity
lambda_1: float lambda_1: float
...@@ -548,7 +548,7 @@ def _base_roq_waveform( ...@@ -548,7 +548,7 @@ def _base_roq_waveform(
theta_jn: float theta_jn: float
Orbital inclination Orbital inclination
phase: float phase: float
The phase at coalescence The phase at reference frequency or peak amplitude (depends on waveform)
Waveform arguments Waveform arguments
=================== ===================
...@@ -649,7 +649,7 @@ def binary_black_hole_frequency_sequence( ...@@ -649,7 +649,7 @@ def binary_black_hole_frequency_sequence(
theta_jn: float theta_jn: float
Angle between the total binary angular momentum and the line of sight Angle between the total binary angular momentum and the line of sight
phase: float phase: float
The phase at coalescence The phase at reference frequency or peak amplitude (depends on waveform)
kwargs: dict kwargs: dict
Required keyword arguments Required keyword arguments
- frequencies: - frequencies:
...@@ -735,7 +735,7 @@ def binary_neutron_star_frequency_sequence( ...@@ -735,7 +735,7 @@ def binary_neutron_star_frequency_sequence(
theta_jn: float theta_jn: float
Angle between the total binary angular momentum and the line of sight Angle between the total binary angular momentum and the line of sight
phase: float phase: float
The phase at coalescence The phase at reference frequency or peak amplitude (depends on waveform)
kwargs: dict kwargs: dict
Required keyword arguments Required keyword arguments
- frequencies: - frequencies:
...@@ -811,7 +811,7 @@ def _base_waveform_frequency_sequence( ...@@ -811,7 +811,7 @@ def _base_waveform_frequency_sequence(
theta_jn: float theta_jn: float
Orbital inclination Orbital inclination
phase: float phase: float
The phase at coalescence The phase at reference frequency or peak amplitude (depends on waveform)
waveform_kwargs: dict waveform_kwargs: dict
Optional keyword arguments Optional keyword arguments
......
"""
This tutorial demonstrates how we can sample a prior in the shape of a ball.
Note that this will not end up sampling uniformly in that space, constraint
priors are more suitable for that.
This implementation will draw a value for the x-coordinate from p(x), and
given that draw a value for the
y-coordinate from p(y|x), and given that draw a value for the z-coordinate
from p(z|x,y).
Only the x-coordinate will end up being uniform for this
"""
import bilby
import numpy as np
class ZeroLikelihood(bilby.core.likelihood.Likelihood):
"""Flat likelihood. This always returns 0.
This way our posterior distribution is exactly the prior distribution."""
def log_likelihood(self):
return 0
def condition_func_y(reference_params, x):
"""Condition function for our p(y|x) prior."""
radius = 0.5 * (reference_params["maximum"] - reference_params["minimum"])
y_max = np.sqrt(radius**2 - x**2)
return dict(minimum=-y_max, maximum=y_max)
def condition_func_z(reference_params, x, y):
"""Condition function for our p(z|x, y) prior."""
radius = 0.5 * (reference_params["maximum"] - reference_params["minimum"])
z_max = np.sqrt(radius**2 - x**2 - y**2)
return dict(minimum=-z_max, maximum=z_max)
# Set up the conditional priors and the flat likelihood
priors = bilby.core.prior.ConditionalPriorDict()
priors["x"] = bilby.core.prior.Uniform(minimum=-1, maximum=1, latex_label="$x$")
priors["y"] = bilby.core.prior.ConditionalUniform(
condition_func=condition_func_y, minimum=-1, maximum=1, latex_label="$y$"
)
priors["z"] = bilby.core.prior.ConditionalUniform(
condition_func=condition_func_z, minimum=-1, maximum=1, latex_label="$z$"
)
likelihood = ZeroLikelihood(parameters=dict(x=0, y=0, z=0))
# Sample the prior distribution
res = bilby.run_sampler(
likelihood=likelihood,
priors=priors,
sampler="dynesty",
nlive=5000,
label="conditional_prior",
outdir="outdir",
resume=False,
clean=True,
)
res.plot_corner()
...@@ -32,6 +32,10 @@ start_time = end_time - duration ...@@ -32,6 +32,10 @@ start_time = end_time - duration
# The fiducial parameters are taken to me the max likelihood sample from the # The fiducial parameters are taken to me the max likelihood sample from the
# posterior sample release of LIGO-Virgo # posterior sample release of LIGO-Virgo
# https://www.gw-openscience.org/eventapi/html/O3_Discovery_Papers/GW190425/ # https://www.gw-openscience.org/eventapi/html/O3_Discovery_Papers/GW190425/
# The fiducial parameters should always be in provided in the same basis as
# the sampling basis. For example, if sampling in `mass_1` and `mass_2` instead of
# `chirp_mass` and `mass_ratio`, the fiducial parameters should also be provided in
# `mass_1` and `mass_2` below.
fiducial_parameters = { fiducial_parameters = {
"a_1": 0.018, "a_1": 0.018,
......
...@@ -7,6 +7,7 @@ This example estimates the masses using a uniform prior in both component masses ...@@ -7,6 +7,7 @@ This example estimates the masses using a uniform prior in both component masses
and distance using a uniform in comoving volume prior on luminosity distance and distance using a uniform in comoving volume prior on luminosity distance
between luminosity distances of 100Mpc and 5Gpc, the cosmology is Planck15. between luminosity distances of 100Mpc and 5Gpc, the cosmology is Planck15.
""" """
from copy import deepcopy
import bilby import bilby
import numpy as np import numpy as np
...@@ -109,6 +110,17 @@ priors["fiducial"] = 0 ...@@ -109,6 +110,17 @@ priors["fiducial"] = 0
# Perform a check that the prior does not extend to a parameter space longer than the data # Perform a check that the prior does not extend to a parameter space longer than the data
priors.validate_prior(duration, minimum_frequency) priors.validate_prior(duration, minimum_frequency)
# Set up the fiducial parameters for the relative binning likelihood to be the
# injected parameters. Note that because we sample in chirp mass and mass ratio
# but injected with mass_1 and mass_2, we need to convert the mass parameters
fiducial_parameters = injection_parameters.copy()
m1 = fiducial_parameters.pop("mass_1")
m2 = fiducial_parameters.pop("mass_2")
fiducial_parameters["chirp_mass"] = bilby.gw.conversion.component_masses_to_chirp_mass(
m1, m2
)
fiducial_parameters["mass_ratio"] = m2 / m1
# Initialise the likelihood by passing in the interferometer data (ifos) and # Initialise the likelihood by passing in the interferometer data (ifos) and
# the waveform generator # the waveform generator
likelihood = bilby.gw.likelihood.RelativeBinningGravitationalWaveTransient( likelihood = bilby.gw.likelihood.RelativeBinningGravitationalWaveTransient(
...@@ -116,15 +128,15 @@ likelihood = bilby.gw.likelihood.RelativeBinningGravitationalWaveTransient( ...@@ -116,15 +128,15 @@ likelihood = bilby.gw.likelihood.RelativeBinningGravitationalWaveTransient(
waveform_generator=waveform_generator, waveform_generator=waveform_generator,
priors=priors, priors=priors,
distance_marginalization=True, distance_marginalization=True,
fiducial_parameters=injection_parameters, fiducial_parameters=fiducial_parameters,
) )
# Run sampler. In this case we're going to use the `dynesty` sampler # Run sampler. In this case, we're going to use the `nestle` sampler
result = bilby.run_sampler( result = bilby.run_sampler(
likelihood=likelihood, likelihood=likelihood,
priors=priors, priors=priors,
sampler="nestle", sampler="nestle",
npoints=100, npoints=1000,
injection_parameters=injection_parameters, injection_parameters=injection_parameters,
outdir=outdir, outdir=outdir,
label=label, label=label,
...@@ -158,5 +170,15 @@ print( ...@@ -158,5 +170,15 @@ print(
) )
print(f"Binned vs unbinned log Bayes factor {np.log(np.mean(weights)):.2f}") print(f"Binned vs unbinned log Bayes factor {np.log(np.mean(weights)):.2f}")
# Make a corner plot. # Generate result object with the posterior for the regular likelihood using
# result.plot_corner() # rejection sampling
alt_result = deepcopy(result)
keep = weights > np.random.uniform(0, max(weights), len(weights))
alt_result.posterior = result.posterior.iloc[keep]
# Make a comparison corner plot.
bilby.core.result.plot_multiple(
[result, alt_result],
labels=["Binned", "Reweighted"],
filename=f"{outdir}/{label}_corner.png",
)
%% Cell type:markdown id: tags:
# Conditional prior demonstration
Conditional priors enable inference to be performed with priors that correlate different parameters.
In this notebook, we demonstrate two uses of this: maintaining a two-dimensional distribution while changing the parameterization to be more efficient for sampling, and enforcing an ordering between parameters.
Many cases where `Conditional` priors are useful can also be expressed with `Constraint` priors, however the conditional approach can improve sampling efficiency.
%% Cell type:code id: tags:
```
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from bilby.core.prior import (
Prior, PriorDict, ConditionalPriorDict,
Uniform, ConditionalUniform, Constraint,
)
from corner import corner
from scipy.stats import semicircular
%matplotlib inline
```
%% Cell type:markdown id: tags:
### Sampling from a disc
Our first example is sampling uniformly from a disc
$$p(x, y) = \frac{1}{\pi}; x^2 + y+2 \leq 1.$$
A naive implementation of this would define a uniform prior over a square over `[-1, 1]` and then reject points that don't satisfy the radius constraint.
Naive sampling from this parameterization would have an efficiency of $\pi / 4$.
If we instead consider the marginal distribution $p(x)$ and conditional distribution $p(y | x)$, we can achieve a sampling efficiency of 100%.
$$
p(x) = \int_{-1 + \sqrt{x^2 + y^2}}^{1 - \sqrt{x^2 + y^2}} dy p(x, y) = \frac{2 (1 - \sqrt{x^2 + y^2})}{\pi} \\
p(y | x) = \frac{1}{2 (1 - \sqrt{x^2})}
$$
The marginal distribution for $x$ is the [Wigner semicircle distribution](https://en.wikipedia.org/wiki/Wigner_semicircle_distribution), this distribution is not currently defined in `Bilby`, but we can wrap the `scipy` implementation.
The conditional distribution for $y$ is implemented in `Bilby` as the `ConditionUnifrom`, we just need to define the condition function.
%% Cell type:code id: tags:
```
class SemiCircular(Prior):
def __init__(self, radius=1, center=0, name=None, latex_label=None, unit=None, boundary=None):
super(SemiCircular, self).__init__(
minimum=center - radius,
maximum=center + radius,
name=name,
latex_label=latex_label,
unit=unit,
boundary=boundary,
)
self.radius = radius
self.center = center
self._dist = semicircular(loc=center, scale=radius)
def prob(self, val):
return self._dist.pdf(val)
def ln_prob(self, val):
return self._dist.logpdf(val)
def cdf(self, val):
return self._dist.cdf(val)
def rescale(self, val):
return self._dist.ppf(val)
def conditional_func_y(reference_parameters, x):
condition = np.sqrt(reference_parameters["maximum"]-x**2)
return dict(minimum=-condition, maximum=condition)
```
%% Cell type:markdown id: tags:
#### Sample from the distribution
To demonstrate the equivalence of the two methods, we will draw samples from the distribution using the two methods and verify that they agree.
%% Cell type:code id: tags:
```
N = int(2e4)
CORNER_KWARGS = dict(
plot_contours=False,
plot_density=False,
fill_contours=False,
max_n_ticks=3,
verbose=False,
use_math_text=True,
)
def convert_to_radial(parameters):
p = parameters.copy()
p['r'] = p['x']**2 + p['y']**2
return p
def sample_circle_with_constraint():
d = PriorDict(
dictionary=dict(
x=Uniform(-1, 1),
y=Uniform(-1, 1),
r=Constraint(0, 1),
),
conversion_function=convert_to_radial
)
return pd.DataFrame(d.sample(N))
def sample_circle_with_conditional():
d = ConditionalPriorDict(
dictionary=dict(
x=SemiCircular(),
y=ConditionalUniform(
condition_func=conditional_func_y,
minimum=-1, maximum=1
)
)
)
return pd.DataFrame(d.sample(N))
s1 = sample_circle_with_constraint()
s2 = sample_circle_with_conditional()
fig = corner(s1.values, **CORNER_KWARGS, color="tab:blue", labels=["$x$", "$y$"])
corner(s2.values, **CORNER_KWARGS, color="tab:green", fig=fig)
plt.show()
plt.close()
```
%% Cell type:markdown id: tags:
### Sampling from ordered distributions
As our second example, we demonstrate defining a prior distribution over a set of strictly ordered parameters.
We note that in this case, we do not require that the marginal distributions for each of the parameters are independently and identically disributed, although this can be fairly simply remedied.
%% Cell type:code id: tags:
```
class BoundedUniform(ConditionalUniform):
"""Conditional Uniform prior where prior sample < previous prior sample
This is ensured by fixing the maximum bound to be the previous prior sample value.
"""
def __init__(self, idx: int, minimum, maximum, name=None, latex_label=None,
unit=None, boundary=None):
super(BoundedUniform, self).__init__(
minimum=minimum, maximum=maximum, name=name,
latex_label=latex_label, unit=unit,
boundary=boundary, condition_func=self.bounds_condition
)
self.idx = idx
self.previous_name = f"{name[:-1]}{self.idx - 1}"
self._required_variables = [self.previous_name]
# this is used in prior.sample(... **required_variables)
def bounds_condition(self, reference_params, **required_variables):
previous_sample = required_variables[self.previous_name]
return dict(maximum=previous_sample)
def make_uniform_conditonal_priordict(n_priors=3):
priors = ConditionalPriorDict()
for i in range(n_priors):
if i==0:
priors[f"uni{i}"] = Uniform(minimum=0, maximum=1, name=f"uni{i}")
else:
priors[f"uni{i}"] = BoundedUniform(idx=i, minimum=0, maximum=1, name=f"uni{i}")
return priors
```
%% Cell type:code id: tags:
```
samples = pd.DataFrame(make_uniform_conditonal_priordict(3).sample(10000))
fig = corner(samples.values, **CORNER_KWARGS, color="tab:blue", labels=[f"$A_{ii}$" for ii in range(3)])
plt.show()
plt.close()
```