CHANGELOG.md

 # All notable changes will be documented in this file
 
+## [2.0.1] 2023-03-13
+Version 2.0.1 release of Bilby
+
+This is a bugfix release after the last major update.
+
+Users may notice changes in inferred binary neutron star masses after updating to match [lalsuite](https://git.ligo.org/lscsoft/lalsuite/-/merge_requests/1658).
+
+### Changes
+- Make sure quantities that need to be conserved between dynesty iterations are class-level attributes (!1225).
+- Fix massive memory usage in post-processing calculation of SNRs (!1227).
+- Update value for the solar mass (!1229).
+- Make `scikit-learn` an explicit dependence of `bilby[GW]` (!1230).
+
 ## [2.0.0] 2023-02-29
 Version 2.0.0 release of Bilby
 

bilby/core/sampler/dynesty.py

@@ -450,32 +450,32 @@ class Dynesty(NestedSampler):
                 f"An average of {2 * self.nact} steps will be accepted up to chain length "
                 f"{self.maxmcmc}."
             )
-            from .dynesty_utils import sample_rwalk_bilby
+            from .dynesty_utils import AcceptanceTrackingRWalk
 
             if self.kwargs["walks"] > self.maxmcmc:
                 raise DynestySetupError("You have maxmcmc < walks (minimum mcmc)")
             if self.nact < 1:
                 raise DynestySetupError("Unable to run with nact < 1")
-            dynesty.nestedsamplers._SAMPLING["rwalk"] = sample_rwalk_bilby
+            AcceptanceTrackingRWalk.old_act = None
+            dynesty.nestedsamplers._SAMPLING["rwalk"] = AcceptanceTrackingRWalk()
         elif sample == "acceptance-walk":
             logger.info(
                 "Using the bilby-implemented rwalk sampling with an average of "
                 f"{self.naccept} accepted steps per MCMC and maximum length {self.maxmcmc}"
             )
-            from .dynesty_utils import fixed_length_rwalk_bilby
+            from .dynesty_utils import FixedRWalk
 
-            dynesty.nestedsamplers._SAMPLING[
-                "acceptance-walk"
-            ] = fixed_length_rwalk_bilby
+            dynesty.nestedsamplers._SAMPLING["acceptance-walk"] = FixedRWalk()
         elif sample == "act-walk":
             logger.info(
                 "Using the bilby-implemented rwalk sampling tracking the "
                 f"autocorrelation function and thinning by "
                 f"{self.nact} with maximum length {self.nact * self.maxmcmc}"
             )
-            from .dynesty_utils import bilby_act_rwalk
+            from .dynesty_utils import ACTTrackingRWalk
 
-            dynesty.nestedsamplers._SAMPLING["act-walk"] = bilby_act_rwalk
+            ACTTrackingRWalk._cache = list()
+            dynesty.nestedsamplers._SAMPLING["act-walk"] = ACTTrackingRWalk()
         elif sample == "rwalk_dynesty":
             sample = sample.strip("_dynesty")
             self.kwargs["sample"] = sample
@@ -901,9 +901,11 @@ class Dynesty(NestedSampler):
         """Run the sampler very briefly as a sanity test that it works."""
         import pandas as pd
 
+        self._set_sampling_method()
+        self._setup_pool()
         self.sampler = self.sampler_init(
-            loglikelihood=self.log_likelihood,
-            prior_transform=self.prior_transform,
+            loglikelihood=_log_likelihood_wrapper,
+            prior_transform=_prior_transform_wrapper,
             ndim=self.ndim,
             **self.sampler_init_kwargs,
         )
@@ -915,6 +917,8 @@ class Dynesty(NestedSampler):
 
             self.pbar = tqdm(file=sys.stdout, initial=self.sampler.it)
         self.sampler.run_nested(**sampler_kwargs)
+        self._close_pool()
+
         if self.pbar is not None:
             self.pbar = self.pbar.close()
             print("")

bilby/core/sampler/dynesty_utils.py

@@ -181,8 +181,11 @@ class ACTTrackingRWalk:
     parallel process.
     """
 
+    # the _cache is a class level variable to avoid being forgotten at every
+    # iteration when using multiprocessing
+    _cache = list()
+
     def __init__(self):
-        self._cache = list()
         self.act = 1
         self.thin = getattr(_SamplingContainer, "nact", 2)
         self.maxmcmc = getattr(_SamplingContainer, "maxmcmc", 5000) * 50
@@ -367,10 +370,13 @@ class AcceptanceTrackingRWalk:
     corresponds to specifying :code:`sample="rwalk"`
     """
 
+    # to retain state between calls to pool.Map, this needs to be a class
+    # level attribute
+    old_act = None
+
     def __init__(self, old_act=None):
         self.maxmcmc = getattr(_SamplingContainer, "maxmcmc", 5000)
         self.nact = getattr(_SamplingContainer, "nact", 40)
-        self.old_act = old_act
 
     def __call__(self, args):
         rstate = get_random_generator(args.rseed)
@@ -437,7 +443,7 @@ class AcceptanceTrackingRWalk:
             logl = args.loglikelihood(v)
 
         blob = {"accept": accept, "reject": reject + nfail, "scale": args.scale}
-        self.old_act = act
+        AcceptanceTrackingRWalk.old_act = act
 
         ncall = accept + reject
         return u, v, logl, ncall, blob
@@ -708,7 +714,3 @@ def apply_boundaries_(u_prop, periodic, reflective):
 
 
 proposal_funcs = dict(diff=propose_differetial_evolution, volumetric=propose_volumetric)
-
-fixed_length_rwalk_bilby = FixedRWalk()
-bilby_act_rwalk = ACTTrackingRWalk()
-sample_rwalk_bilby = AcceptanceTrackingRWalk()

bilby/core/utils/constants.py

@@ -2,6 +2,6 @@
 
 speed_of_light = 299792458.0  # m/s
 parsec = 3.085677581491367e+16  # m
-solar_mass = 1.9884099021470415e+30  # Kg
+solar_mass = 1.988409870698050731911960804878414216e30  # Kg
 radius_of_earth = 6378136.6  # m
 gravitational_constant = 6.6743e-11  # m^3 kg^-1 s^-2

bilby/gw/conversion.py

@@ -1465,7 +1465,7 @@ def _compute_snrs(args):
     )
     snrs = list()
     for ifo in likelihood.interferometers:
-        snrs.append(likelihood.calculate_snrs(signal_polarizations, ifo))
+        snrs.append(likelihood.calculate_snrs(signal_polarizations, ifo, return_array=False))
     return snrs
 
 

bilby/gw/likelihood/base.py

@@ -236,16 +236,26 @@ class GravitationalWaveTransient(Likelihood):
                     "waveform_generator.".format(attribute))
             setattr(self.waveform_generator, attribute, ifo_attr)
 
-    def calculate_snrs(self, waveform_polarizations, interferometer):
+    def calculate_snrs(self, waveform_polarizations, interferometer, return_array=True):
         """
         Compute the snrs
 
         Parameters
-        ==========
+        ----------
         waveform_polarizations: dict
             A dictionary of waveform polarizations and the corresponding array
         interferometer: bilby.gw.detector.Interferometer
             The bilby interferometer object
+        return_array: bool
+            If true, calculate and return internal array objects
+            (d_inner_h_array and optimal_snr_squared_array), otherwise
+            these are returned as None.
+
+        Returns
+        -------
+        calculated_snrs: _CalculatedSNRs
+            An object containing the SNR quantities and (if return_array=True)
+            the internal array objects.
 
         """
         signal = self._compute_full_waveform(
@@ -266,7 +276,10 @@ class GravitationalWaveTransient(Likelihood):
 
         normalization = 4 / self.waveform_generator.duration
 
-        if self.time_marginalization and self.calibration_marginalization:
+        if return_array is False:
+            d_inner_h_array = None
+            optimal_snr_squared_array = None
+        elif self.time_marginalization and self.calibration_marginalization:
 
             d_inner_h_integrand = np.tile(
                 interferometer.frequency_domain_strain.conjugate() * signal /

bilby/gw/likelihood/multiband.py

@@ -693,18 +693,26 @@ class MBGravitationalWaveTransient(GravitationalWaveTransient):
             else:
                 setattr(self, key, value)
 
-    def calculate_snrs(self, waveform_polarizations, interferometer):
+    def calculate_snrs(self, waveform_polarizations, interferometer, return_array=False):
         """
-        Compute the snrs for multi-banding
+        Compute the snrs
 
         Parameters
         ----------
-        waveform_polarizations: waveform
+        waveform_polarizations: dict
+            A dictionary of waveform polarizations and the corresponding array
         interferometer: bilby.gw.detector.Interferometer
+            The bilby interferometer object
+        return_array: bool
+            If true, calculate and return internal array objects
+            (d_inner_h_array and optimal_snr_squared_array), otherwise
+            these are returned as None. This parameter is ignored for the multiband
+            model as these arrays are never calculated.
 
         Returns
         -------
-        snrs: named tuple of snrs
+        calculated_snrs: _CalculatedSNRs
+            An object containing the SNR quantities.
 
         """
         strain = np.zeros(len(self.banded_frequency_points), dtype=complex)

bilby/gw/likelihood/relative.py

@@ -382,7 +382,7 @@ class RelativeBinningGravitationalWaveTransient(GravitationalWaveTransient):
         full_waveform_ratio = duplicated_r0 + duplicated_r1 * (f - duplicated_fm)
         return fiducial_waveform * full_waveform_ratio
 
-    def calculate_snrs(self, waveform_polarizations, interferometer):
+    def calculate_snrs(self, waveform_polarizations, interferometer, return_array=True):
         r0, r1 = self.compute_waveform_ratio_per_interferometer(
             waveform_polarizations=waveform_polarizations,
             interferometer=interferometer,
@@ -393,7 +393,7 @@ class RelativeBinningGravitationalWaveTransient(GravitationalWaveTransient):
         optimal_snr_squared = h_inner_h
         complex_matched_filter_snr = d_inner_h / (optimal_snr_squared ** 0.5)
 
-        if self.time_marginalization:
+        if return_array and self.time_marginalization:
             full_waveform = self._compute_full_waveform(
                 signal_polarizations=waveform_polarizations,
                 interferometer=interferometer,

bilby/gw/likelihood/roq.py

@@ -394,7 +394,7 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
     def waveform_generator(self, waveform_generator):
         self._waveform_generator = waveform_generator
 
-    def calculate_snrs(self, waveform_polarizations, interferometer):
+    def calculate_snrs(self, waveform_polarizations, interferometer, return_array=True):
         """
         Compute the snrs for ROQ
 
@@ -458,7 +458,7 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
             with np.errstate(invalid="ignore"):
                 complex_matched_filter_snr = d_inner_h / (optimal_snr_squared**0.5)
 
-        if self.time_marginalization:
+        if return_array and self.time_marginalization:
             ifo_times = self._times - interferometer.strain_data.start_time
             ifo_times += dt
             if self.jitter_time:

gw_requirements.txt

@@ -3,3 +3,4 @@ lalsuite
 gwpy
 tables
 pyfftw
+scikit-learn

test/core/sampler/dynesty_test.py

 import unittest
 from copy import deepcopy
-from unittest.mock import MagicMock
 
 from attr import define
 import bilby
@@ -22,9 +21,22 @@ class Dummy:
     loglstar: float = -1
 
 
+class DummyLikelihood(bilby.core.likelihood.Likelihood):
+    """
+    A trivial likelihood used for testing. Add some randomness so the likelihood
+    isn't flat everywhere as that can cause issues for nested samplers.
+    """
+
+    def __init__(self):
+        super().__init__(dict())
+
+    def log_likelihood(self):
+        return np.random.uniform(0, 0.01)
+
+
 class TestDynesty(unittest.TestCase):
     def setUp(self):
-        self.likelihood = MagicMock()
+        self.likelihood = DummyLikelihood()
         self.priors = bilby.core.prior.PriorDict(
             dict(a=bilby.core.prior.Uniform(0, 1), b=bilby.core.prior.Uniform(0, 1))
         )
@@ -83,6 +95,9 @@ class TestDynesty(unittest.TestCase):
         self.assertEqual([0, 4], self.sampler.kwargs["periodic"])
         self.assertEqual([1, 3], self.sampler.kwargs["reflective"])
 
+    def test_run_test_runs(self):
+        self.sampler._run_test()
+
 
 class ProposalsTest(unittest.TestCase):
 
@@ -242,6 +257,7 @@ class TestEstimateNMCMC(unittest.TestCase):
         safety * (2 / accept_ratio - 1)
         """
         sampler = dynesty_utils.AcceptanceTrackingRWalk()
+        dynesty_utils.AcceptanceTrackingRWalk.old_act = None
         for _ in range(10):
             accept_ratio = np.random.uniform()
             safety = np.random.randint(2, 8)

test/gw/utils_test.py

@@ -8,6 +8,7 @@ import lalsimulation as lalsim
 from gwpy.timeseries import TimeSeries
 from gwpy.detector import Channel
 from scipy.stats import ks_2samp
+import pytest
 
 import bilby
 from bilby.gw import utils as gwutils
@@ -72,6 +73,7 @@ class TestGWUtils(unittest.TestCase):
         )
         self.assertEqual(mfsnr, 25.510869054168282)
 
+    @pytest.mark.skip(reason="GWOSC unstable: avoiding this test")
     def test_get_event_time(self):
         from urllib3.exceptions import NewConnectionError
         events = [