.gitlab-ci.yml

@@ -17,17 +17,25 @@ stages:
 
 # ------------------- Initial stage -------------------------------------------
 
+.list-env: &list-env
+  - PREFIX="$(dirname $(which python))/.."
+  - if [ -d "${PREFIX}/conda-meta" ]; then
+      conda list --prefix "${PREFIX}" --show-channel-urls;
+    else
+      python -m pip list installed;
+    fi
+
 # Check author list is up to date
 authors:
   stage: initial
-  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python39
+  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python310
   script:
     - python test/check_author_list.py
 
 # Test containers scripts are up to date
 containers:
   stage: initial
-  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python39
+  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python310
   script:
     - cd containers
     - python write_dockerfiles.py #HACK
@@ -40,7 +48,7 @@ containers:
   image: python
   script:
     - python -m pip install .
-    - python -m pip list installed
+    - *list-env
     - python -c "import bilby"
     - python -c "import bilby.bilby_mcmc"
     - python -c "import bilby.core"
@@ -59,10 +67,6 @@ containers:
           ${script} --help;
       done
 
-basic-3.8:
-  <<: *test-python
-  image: python:3.8
-
 basic-3.9:
   <<: *test-python
   image: python:3.9
@@ -75,13 +79,9 @@ basic-3.10:
   stage: initial
   script:
     - python -m pip install .
-    - python -m pip list installed
+    - *list-env
     - python test/test_samplers_import.py
 
-import-samplers-3.8:
-  <<: *test-samplers-import
-  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python38
-
 import-samplers-3.9:
   <<: *test-samplers-import
   image: containers.ligo.org/lscsoft/bilby/v2-bilby-python39
@@ -102,13 +102,6 @@ import-samplers-3.10:
     # Run precommits (flake8, spellcheck, isort, no merge conflicts, etc)
     - pre-commit run --all-files --verbose --show-diff-on-failure
 
-precommits-py3.8:
-  <<: *precommits
-  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python38
-  variables:
-    CACHE_DIR: ".pip38"
-    PYVERSION: "python38"
-
 precommits-py3.9:
   <<: *precommits
   image: containers.ligo.org/lscsoft/bilby/v2-bilby-python39
@@ -128,7 +121,7 @@ install:
   parallel:
     matrix:
       - EXTRA: [gw, mcmc, all]
-  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python39
+  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python310
   script:
     - pip install .[$EXTRA]
 
@@ -138,19 +131,19 @@ install:
   stage: test
   script:
     - python -m pip install .
-    - python -m pip list installed
+    - *list-env
 
     - pytest --cov=bilby --durations 10
 
-python-3.8:
-  <<: *unit-test
-  needs: ["basic-3.8", "precommits-py3.8"]
-  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python38
-
 python-3.9:
   <<: *unit-test
   needs: ["basic-3.9", "precommits-py3.9"]
   image: containers.ligo.org/lscsoft/bilby/v2-bilby-python39
+
+python-3.10:
+  <<: *unit-test
+  needs: ["basic-3.10", "precommits-py3.10"]
+  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python310
   after_script:
     - coverage html
     - coverage xml
@@ -164,23 +157,13 @@ python-3.9:
       - htmlcov/
     expire_in: 30 days
 
-python-3.10:
-  <<: *unit-test
-  needs: ["basic-3.10", "precommits-py3.10"]
-  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python310
-
 .test-sampler: &test-sampler
   stage: test
   script:
     - python -m pip install .[all]
-    - python -m pip list installed
+    - *list-env
     - pytest test/integration/sampler_run_test.py --durations 10 -v
 
-python-3.8-samplers:
-  <<: *test-sampler
-  needs: ["basic-3.8", "precommits-py3.8"]
-  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python38
-
 python-3.9-samplers:
   <<: *test-sampler
   needs: ["basic-3.9", "precommits-py3.9"]
@@ -191,15 +174,15 @@ python-3.10-samplers:
   needs: ["basic-3.10", "precommits-py3.10"]
   image: containers.ligo.org/lscsoft/bilby/v2-bilby-python310
 
-integration-tests-python-3.9:
+integration-tests-python-3.10:
   stage: test
-  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python39
-  needs: ["basic-3.9", "precommits-py3.9"]
+  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python310
+  needs: ["basic-3.10", "precommits-py3.10"]
   only:
     - schedules
   script:
     - python -m pip install .
-    - python -m pip list installed
+    - *list-env
     # Run tests which are only done on schedule
     - pytest test/integration/example_test.py
 
@@ -209,15 +192,9 @@ integration-tests-python-3.9:
     - schedules
   script:
     - python -m pip install .
-    - python -m pip list installed
+    - *list-env
     - pytest test/gw/plot_test.py
 
-
-plotting-python-3.8:
-  <<: *plotting
-  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python38
-  needs: ["basic-3.8", "precommits-py3.8"]
-
 plotting-python-3.9:
   <<: *plotting
   image: containers.ligo.org/lscsoft/bilby/v2-bilby-python39
@@ -232,7 +209,7 @@ plotting-python-3.10:
 
 docs:
   stage: docs
-  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python39
+  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python310
   before_script:
     - conda install -c conda-forge pandoc ipython jupyter nbconvert
     - python -m pip install ipykernel
@@ -254,7 +231,7 @@ docs:
 
 pages:
   stage: deploy
-  needs: ["docs", "python-3.9"]
+  needs: ["docs", "python-3.10"]
   script:
     - mkdir public/
     - mv htmlcov/ public/
@@ -285,11 +262,6 @@ pages:
     - docker image tag v3-bilby-$PYVERSION containers.ligo.org/lscsoft/bilby/v2-bilby-$PYVERSION:latest
     - docker image push containers.ligo.org/lscsoft/bilby/v2-bilby-$PYVERSION:latest
 
-build-python38-container:
-  <<: *build-container
-  variables:
-    PYVERSION: "python38"
-
 build-python39-container:
   <<: *build-container
   variables:
@@ -302,7 +274,7 @@ build-python310-container:
 
 pypi-release:
   stage: deploy
-  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python39
+  image: containers.ligo.org/lscsoft/bilby/v2-bilby-python310
   variables:
     TWINE_USERNAME: $PYPI_USERNAME
     TWINE_PASSWORD: $PYPI_PASSWORD

AUTHORS.md

@@ -91,3 +91,4 @@ Isaac Legred
 Marc Penuliar
 Andrew Fowlie
 Martin White
+Peter Tsun-Ho Pang

CHANGELOG.md

 # All notable changes will be documented in this file
 
+## [2.2.0] 2023-07-24
+Version 2.2.0 release of Bilby
+
+This release contains one new feature and drops support for Python 3.8.
+
+### Added
+- New waveform interface to support the SEOBNRv5 family of waveforms (!1218)
+- Enable default noise + injection function for non-CBC signals (!1263)
+- Fallback to result pickle loading to match result writing (!1291)
+
+### Changes
+- Additional error catching for plotting (!1261, !1271)
+- Improve plotting options for corner plots (!1270)
+- Fix bugs in closing the pool for emcee (!1274)
+- Generalize MPI support (!1278)
+- Fix a bug with saving hdf5 results when conda isn't present (!1290)
+
+### Deprecated
+- Drop support for py38 (!1277)
+
+## [2.1.2] 2023-07-17
+Version 2.1.2 release of Bilby
+
+This is a bugfix release.
+Note that one of the changes will have a significant impact on scripts that rely on
+a seed for random data generation.
+Where users have previously used `np.random.seed` they should now call
+`bilby.core.utils.random.seed`.
+
+### Changes
+- Fix issues related to random number generation with multiprocessing (!1273)
+- Enable cosmological priors to be written/read in our plain text format (!1258)
+- Allow posterior reweighting to be performed when changing the likelihood and the prior (!1260)
+
+## [2.1.1] 2023-04-28
+Version 2.1.1 release of Bilby
+
+Bugfix release
+
+### Changes
+- Fix the matched filter SNR phase for the multiband likelihood (!1253)
+- Bugfix for Fisher matrix proposals in `bilby_mcmc` (!1251)
+- Make the changes to the spline calibration backward compatible, 2.0.2 resume files can't be read with 2.1.0 (!1250)
+
 ## [2.1.0] 2023-04-12
 Version 2.1.0 release of Bilby
 

bilby/bilby_mcmc/chain.py

@@ -136,12 +136,14 @@ class Chain(object):
 
     @property
     def _random_idx(self):
+        from ..core.utils.random import rng
+
         mindex = self._last_minimum_index[1]
         # Check if mindex exceeds current position by 10 ACT: if so use a random sample
         # otherwise we draw only from the chain past the minimum_index
         if np.isinf(self.tau_last) or self.position - mindex < 10 * self.tau_last:
             mindex = 0
-        return np.random.randint(mindex, self.position + 1)
+        return rng.integers(mindex, self.position + 1)
 
     @property
     def random_sample(self):
@@ -412,12 +414,20 @@ class Chain(object):
         for ax, key in zip(axes[:, 1], self.keys):
             if all_samples is not None:
                 yy_all = all_samples[key]
-                ax.hist(yy_all, bins=50, alpha=0.6, density=True, color="k")
-
+                if np.any(np.isinf(yy_all)):
+                    logger.warning(
+                        f"Could not plot histogram for parameter {key} due to infinite values"
+                    )
+                else:
+                    ax.hist(yy_all, bins=50, alpha=0.6, density=True, color="k")
             yy = self.get_1d_array(key)[nburn : self.position : self.thin]
-            ax.hist(yy, bins=50, alpha=0.8, density=True)
-
-            ax.set_xlabel(self._get_plot_label_by_key(key, priors))
+            if np.any(np.isinf(yy)):
+                logger.warning(
+                    f"Could not plot histogram for parameter {key} due to infinite values"
+                )
+            else:
+                ax.hist(yy, bins=50, alpha=0.8, density=True)
+                ax.set_xlabel(self._get_plot_label_by_key(key, priors))
 
         # Add x-axes labels to the traceplots
         axes[-1, 0].set_xlabel(r"Iteration $[\times 10^{3}]$")

bilby/bilby_mcmc/proposals.py

@@ -9,7 +9,7 @@ from scipy.stats import gaussian_kde
 from ..core.fisher import FisherMatrixPosteriorEstimator
 from ..core.prior import PriorDict
 from ..core.sampler.base_sampler import SamplerError
-from ..core.utils import logger, reflect
+from ..core.utils import logger, random, reflect
 from ..gw.source import PARAMETER_SETS
 
 
@@ -19,7 +19,7 @@ class ProposalCycle(object):
         self.weights = [prop.weight for prop in self.proposal_list]
         self.normalized_weights = [w / sum(self.weights) for w in self.weights]
         self.weighted_proposal_list = [
-            np.random.choice(self.proposal_list, p=self.normalized_weights)
+            random.rng.choice(self.proposal_list, p=self.normalized_weights)
             for _ in range(10 * int(1 / min(self.normalized_weights)))
         ]
         self.nproposals = len(self.weighted_proposal_list)
@@ -219,7 +219,7 @@ class FixedGaussianProposal(BaseProposal):
         sample = chain.current_sample
         for key in self.parameters:
             sigma = self.prior_width_dict[key] * self.sigmas[key]
-            sample[key] += sigma * np.random.randn()
+            sample[key] += sigma * random.rng.normal(0, 1)
         log_factor = 0
         return sample, log_factor
 
@@ -256,15 +256,15 @@ class AdaptiveGaussianProposal(BaseProposal):
     def propose(self, chain):
         sample = chain.current_sample
         self.update_scale(chain)
-        if np.random.random() < 1e-3:
+        if random.rng.uniform(0, 1) < 1e-3:
             factor = 1e1
-        elif np.random.random() < 1e-4:
+        elif random.rng.uniform(0, 1) < 1e-4:
             factor = 1e2
         else:
             factor = 1
         for key in self.parameters:
             sigma = factor * self.scale * self.prior_width_dict[key] * self.sigmas[key]
-            sample[key] += sigma * np.random.randn()
+            sample[key] += sigma * random.rng.normal(0, 1)
         log_factor = 0
         return sample, log_factor
 
@@ -306,13 +306,13 @@ class DifferentialEvolutionProposal(BaseProposal):
         theta = chain.current_sample
         theta1 = chain.random_sample
         theta2 = chain.random_sample
-        if np.random.rand() > self.mode_hopping_frac:
+        if random.rng.uniform(0, 1) > self.mode_hopping_frac:
             gamma = 1
         else:
             # Base jump size
-            gamma = np.random.normal(0, 2.38 / np.sqrt(2 * self.ndim))
+            gamma = random.rng.normal(0, 2.38 / np.sqrt(2 * self.ndim))
             # Scale uniformly in log between 0.1 and 10 times
-            gamma *= np.exp(np.log(0.1) + np.log(100.0) * np.random.rand())
+            gamma *= np.exp(np.log(0.1) + np.log(100.0) * random.rng.uniform(0, 1))
 
         for key in self.parameters:
             theta[key] += gamma * (theta2[key] - theta1[key])
@@ -347,7 +347,7 @@ class UniformProposal(BaseProposal):
         for key in self.parameters:
             width = self.prior_width_dict[key]
             if np.isinf(width) is False:
-                sample[key] = np.random.uniform(
+                sample[key] = random.rng.uniform(
                     self.prior_minimum_dict[key], self.prior_maximum_dict[key]
                 )
             else:
@@ -778,14 +778,14 @@ class FixedJumpProposal(BaseProposal):
     def propose(self, chain):
         sample = chain.current_sample
         for key, jump in self.jumps.items():
-            sign = np.random.randint(2) * 2 - 1
+            sign = random.rng.integers(2) * 2 - 1
             sample[key] += sign * jump + self.epsilon * self.prior_width_dict[key]
         log_factor = 0
         return sample, log_factor
 
     @property
     def epsilon(self):
-        return self.scale * np.random.normal()
+        return self.scale * random.rng.normal()
 
 
 class FisherMatrixProposal(AdaptiveGaussianProposal):
@@ -805,7 +805,7 @@ class FisherMatrixProposal(AdaptiveGaussianProposal):
         weight=1,
         update_interval=100,
         scale_init=1e0,
-        fd_eps=1e-6,
+        fd_eps=1e-4,
         adapt=False,
     ):
         super(FisherMatrixProposal, self).__init__(
@@ -827,14 +827,14 @@ class FisherMatrixProposal(AdaptiveGaussianProposal):
             )
             try:
                 self.iFIM = fmp.calculate_iFIM(sample.dict)
-            except (RuntimeError, ValueError) as e:
+            except (RuntimeError, ValueError, np.linalg.LinAlgError) as e:
                 logger.warning(f"FisherMatrixProposal failed with {e}")
                 if hasattr(self, "iFIM") is False:
                     # No past iFIM exists, return sample
                     return sample, 0
             self.steps_since_update = 0
 
-        jump = self.scale * np.random.multivariate_normal(
+        jump = self.scale * random.rng.multivariate_normal(
             self.mean, self.iFIM, check_valid="ignore"
         )
 
@@ -897,11 +897,11 @@ class CorrelatedPolarisationPhaseJump(BaseGravitationalWaveTransientProposal):
         alpha = sample["psi"] + phase
         beta = sample["psi"] - phase
 
-        draw = np.random.random()
+        draw = random.rng.random()
         if draw < 0.5:
-            alpha = 3.0 * np.pi * np.random.random()
+            alpha = 3.0 * np.pi * random.rng.random()
         else:
-            beta = 3.0 * np.pi * np.random.random() - 2 * np.pi
+            beta = 3.0 * np.pi * random.rng.random() - 2 * np.pi
 
         # Update
         sample["psi"] = (alpha + beta) * 0.5
@@ -936,7 +936,7 @@ class PhaseReversalProposal(BaseGravitationalWaveTransientProposal):
     @property
     def epsilon(self):
         if self.fuzz:
-            return np.random.normal(0, self.fuzz_sigma)
+            return random.rng.normal(0, self.fuzz_sigma)
         else:
             return 0
 
@@ -1001,7 +1001,7 @@ class StretchProposal(BaseProposal):
 
 def _stretch_move(sample, complement, scale, ndim, parameters):
     # Draw z
-    u = np.random.rand()
+    u = random.rng.uniform(0, 1)
     z = (u * (scale - 1) + 1) ** 2 / scale
 
     log_factor = (ndim - 1) * np.log(z)
@@ -1045,7 +1045,7 @@ class EnsembleStretch(EnsembleProposal):
     def propose(self, chain, chain_complement):
         sample = chain.current_sample
         completement = chain_complement[
-            np.random.randint(len(chain_complement))
+            random.rng.integers(len(chain_complement))
         ].current_sample
         return _stretch_move(
             sample, completement, self.scale, self.ndim, self.parameters

bilby/bilby_mcmc/sampler.py

@@ -16,7 +16,12 @@ from ..core.sampler.base_sampler import (
     _sampling_convenience_dump,
     signal_wrapper,
 )
-from ..core.utils import check_directory_exists_and_if_not_mkdir, logger, safe_file_dump
+from ..core.utils import (
+    check_directory_exists_and_if_not_mkdir,
+    logger,
+    random,
+    safe_file_dump,
+)
 from . import proposals
 from .chain import Chain, Sample
 from .utils import LOGLKEY, LOGPKEY, ConvergenceInputs, ParallelTemperingInputs
@@ -383,7 +388,7 @@ class Bilby_MCMC(MCMCSampler):
 
         with open(self.resume_file, "rb") as file:
             ptsampler = dill.load(file)
-            if type(ptsampler) != BilbyPTMCMCSampler:
+            if not isinstance(ptsampler, BilbyPTMCMCSampler):
                 logger.debug("Malformed resume file, ignoring")
                 return False
             self.ptsampler = ptsampler
@@ -818,7 +823,7 @@ class BilbyPTMCMCSampler(object):
                 with np.errstate(over="ignore"):
                     alpha_swap = np.exp(dbeta * (logli - loglj))
 
-                if np.random.uniform(0, 1) <= alpha_swap:
+                if random.rng.uniform(0, 1) <= alpha_swap:
                     sampleri.chain[-1] = vj
                     samplerj.chain[-1] = vi
                     self.sampler_dictionary[Tindex][Eindex] = sampleri
@@ -852,7 +857,7 @@ class BilbyPTMCMCSampler(object):
                         - curr[LOGPKEY]
                     )
 
-                    if np.random.uniform(0, 1) <= alpha:
+                    if random.rng.uniform(0, 1) <= alpha:
                         sampler.accept_proposal(prop, proposal)
                     else:
                         sampler.reject_proposal(curr, proposal)
@@ -1021,7 +1026,7 @@ class BilbyPTMCMCSampler(object):
         ln_z_realisations = []
         try:
             for _ in range(repeats):
-                idxs = [np.random.randint(i, i + ll) for i in range(steps - ll)]
+                idxs = [random.rng.integers(i, i + ll) for i in range(steps - ll)]
                 ln_z_realisations.append(
                     self._calculate_stepping_stone(betas, ln_likes[idxs, :])[0]
                 )
@@ -1256,7 +1261,7 @@ class BilbyMCMCSampler(object):
                     - curr[LOGPKEY]
                 )
 
-            if np.random.uniform(0, 1) <= alpha:
+            if random.rng.uniform(0, 1) <= alpha:
                 internal_accepted += 1
                 proposal.accepted += 1
                 curr = prop

bilby/core/fisher.py

 import numpy as np
-import scipy.linalg
-
 import pandas as pd
+import scipy.linalg
 from scipy.optimize import minimize
 
 
@@ -61,12 +60,14 @@ class FisherMatrixPosteriorEstimator(object):
         return iFIM
 
     def sample_array(self, sample, n=1):
+        from .utils.random import rng
+
         if sample == "maxL":
             sample = self.get_maximum_likelihood_sample()
 
         self.mean = np.array(list(sample.values()))
         self.iFIM = self.calculate_iFIM(sample)
-        return np.random.multivariate_normal(self.mean, self.iFIM, n)
+        return rng.multivariate_normal(self.mean, self.iFIM, n)
 
     def sample_dataframe(self, sample, n=1):
         samples = self.sample_array(sample, n)

bilby/core/prior/base.py

@@ -7,8 +7,13 @@ import numpy as np
 import scipy.stats
 from scipy.interpolate import interp1d
 
-from bilby.core.utils import infer_args_from_method, BilbyJsonEncoder, decode_bilby_json, logger, \
-    get_dict_with_properties
+from ..utils import (
+    infer_args_from_method,
+    BilbyJsonEncoder,
+    decode_bilby_json,
+    logger,
+    get_dict_with_properties,
+)
 
 
 class Prior(object):
@@ -124,7 +129,9 @@ class Prior(object):
         float: A random number between 0 and 1, rescaled to match the distribution of this Prior
 
         """
-        self.least_recently_sampled = self.rescale(np.random.uniform(0, 1, size))
+        from ..utils.random import rng
+
+        self.least_recently_sampled = self.rescale(rng.uniform(0, 1, size))
         return self.least_recently_sampled
 
     def rescale(self, val):
@@ -222,18 +229,6 @@ class Prior(object):
         else:
             return f"{prior_module}.{prior_name}({args})"
 
-    @property
-    def _repr_dict(self):
-        """
-        Get a dictionary containing the arguments needed to reproduce this object.
-        """
-        property_names = {p for p in dir(self.__class__) if isinstance(getattr(self.__class__, p), property)}
-        subclass_args = infer_args_from_method(self.__init__)
-        dict_with_properties = self.__dict__.copy()
-        for key in property_names.intersection(subclass_args):
-            dict_with_properties[key] = getattr(self, key)
-        return {key: dict_with_properties[key] for key in subclass_args}
-
     @property
     def is_fixed(self):
         """

bilby/core/prior/conditional.py

-import numpy as np
-
 from .base import Prior, PriorException
 from .interpolated import Interped
 from .analytical import DeltaFunction, PowerLaw, Uniform, LogUniform, \
@@ -76,7 +74,9 @@ def conditional_prior_factory(prior_class):
             float: See superclass
 
             """
-            self.least_recently_sampled = self.rescale(np.random.uniform(0, 1, size), **required_variables)
+            from ..utils.random import rng
+
+            self.least_recently_sampled = self.rescale(rng.uniform(0, 1, size), **required_variables)
             return self.least_recently_sampled
 
         def rescale(self, val, **required_variables):

bilby/core/prior/joint.py

@@ -6,6 +6,7 @@ from scipy.special import erfinv
 
 from .base import Prior, PriorException
 from ..utils import logger, infer_args_from_method, get_dict_with_properties
+from ..utils import random
 
 
 class BaseJointPriorDist(object):
@@ -583,7 +584,7 @@ class MultivariateGaussianDist(BaseJointPriorDist):
             if self.nmodes == 1:
                 mode = 0
             else:
-                mode = np.argwhere(self.cumweights - np.random.rand() > 0)[0][0]
+                mode = np.argwhere(self.cumweights - random.rng.uniform(0, 1) > 0)[0][0]
 
         samp = erfinv(2.0 * samp - 1) * 2.0 ** 0.5
 
@@ -604,12 +605,12 @@ class MultivariateGaussianDist(BaseJointPriorDist):
                 mode = 0
             else:
                 if size == 1:
-                    mode = np.argwhere(self.cumweights - np.random.rand() > 0)[0][0]
+                    mode = np.argwhere(self.cumweights - random.rng.uniform(0, 1) > 0)[0][0]
                 else:
                     # pick modes
                     mode = [
                         np.argwhere(self.cumweights - r > 0)[0][0]
-                        for r in np.random.rand(size)
+                        for r in random.rng.uniform(0, 1, size)
                     ]
 
         samps = np.zeros((size, len(self)))
@@ -617,7 +618,7 @@ class MultivariateGaussianDist(BaseJointPriorDist):
             inbound = False
             while not inbound:
                 # sample the multivariate Gaussian keys
-                vals = np.random.uniform(0, 1, len(self))
+                vals = random.rng.uniform(0, 1, len(self))
 
                 if isinstance(mode, list):
                     samp = np.atleast_1d(self.rescale(vals, mode=mode[i]))

bilby/core/result.py

@@ -24,6 +24,7 @@ from .utils import (
     recursively_load_dict_contents_from_group,
     recursively_decode_bilby_json,
     safe_file_dump,
+    random,
 )
 from .prior import Prior, PriorDict, DeltaFunction, ConditionalDeltaFunction
 
@@ -133,6 +134,16 @@ def read_in_result_list(filename_list, invalid="warning"):
     """
     results_list = []
     for filename in filename_list:
+        if (
+            not os.path.exists(filename)
+            and os.path.exists(f"{os.path.splitext(filename)[0]}.pkl")
+        ):
+            pickle_path = f"{os.path.splitext(filename)[0]}.pkl"
+            logger.warning(
+                f"Result file {filename} doesn't exist but {pickle_path} does. "
+                f"Using {pickle_path}."
+            )
+            filename = pickle_path
         try:
             results_list.append(read_in_result(filename=filename))
         except Exception as e:
@@ -219,10 +230,15 @@ def get_weights_for_reweighting(
                                      desc='Computing priors',
                                      total=n),
                                 start=starting_index):
+        # prior calculation needs to not have prior or likelihood keys
+        ln_prior = sample.pop("log_prior", np.nan)
+        if "log_likelihood" in sample:
+            del sample["log_likelihood"]
+
         if old_prior is not None:
             old_log_prior_array[ii] = old_prior.ln_prob(sample)
         else:
-            old_log_prior_array[ii] = sample["log_prior"]
+            old_log_prior_array[ii] = ln_prior
 
         if new_prior is not None:
             new_log_prior_array[ii] = new_prior.ln_prob(sample)
@@ -259,7 +275,7 @@ def rejection_sample(posterior, weights):
         The posterior resampled using rejection sampling
 
     """
-    keep = weights > np.random.uniform(0, max(weights), weights.shape)
+    keep = weights > random.rng.uniform(0, max(weights), weights.shape)
     return posterior[keep]
 
 
@@ -1177,22 +1193,33 @@ class Result(object):
         if utils.command_line_args.bilby_test_mode:
             return
 
-        # bilby default corner kwargs. Overwritten by anything passed to kwargs
         defaults_kwargs = dict(
-            bins=50, smooth=0.9, label_kwargs=dict(fontsize=16),
+            bins=50, smooth=0.9,
             title_kwargs=dict(fontsize=16), color='#0072C1',
             truth_color='tab:orange', quantiles=[0.16, 0.84],
             levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.)),
             plot_density=False, plot_datapoints=True, fill_contours=True,
-            max_n_ticks=3, hist_kwargs=dict(density=True))
+            max_n_ticks=3)
 
         if 'lionize' in kwargs and kwargs['lionize'] is True:
             defaults_kwargs['truth_color'] = 'tab:blue'
             defaults_kwargs['color'] = '#FF8C00'
 
+        label_kwargs_defaults = dict(fontsize=16)
+        hist_kwargs_defaults = dict(density=True)
+
+        label_kwargs_input = kwargs.get("label_kwargs", dict())
+        hist_kwargs_input = kwargs.get("hist_kwargs", dict())
+
+        label_kwargs_defaults.update(label_kwargs_input)
+        hist_kwargs_defaults.update(hist_kwargs_input)
+
         defaults_kwargs.update(kwargs)
         kwargs = defaults_kwargs
 
+        kwargs["label_kwargs"] = label_kwargs_defaults
+        kwargs["hist_kwargs"] = hist_kwargs_defaults
+
         # Handle if truths was passed in
         if 'truth' in kwargs:
             kwargs['truths'] = kwargs.pop('truth')
@@ -1870,7 +1897,7 @@ class ResultList(list):
         result_weights = np.exp(log_evidences - np.max(log_evidences))
         posteriors = list()
         for res, frac in zip(self, result_weights):
-            selected_samples = (np.random.uniform(size=len(res.posterior)) < frac)
+            selected_samples = (random.rng.uniform(size=len(res.posterior)) < frac)
             posteriors.append(res.posterior[selected_samples])
 
         # remove original nested_samples
@@ -1967,7 +1994,8 @@ class ResultList(list):
 
 @latex_plot_format
 def plot_multiple(results, filename=None, labels=None, colours=None,
-                  save=True, evidences=False, corner_labels=None, **kwargs):
+                  save=True, evidences=False, corner_labels=None, linestyles=None,
+                  **kwargs):
     """ Generate a corner plot overlaying two sets of results
 
     Parameters
@@ -2021,11 +2049,17 @@ def plot_multiple(results, filename=None, labels=None, colours=None,
             c = colours[i]
         else:
             c = 'C{}'.format(i)
+        if linestyles is not None:
+            linestyle = linestyles[i]
+        else:
+            linestyle = 'solid'
         hist_kwargs = kwargs.get('hist_kwargs', dict())
         hist_kwargs['color'] = c
-        fig = result.plot_corner(fig=fig, save=False, color=c, **kwargs)
+        hist_kwargs["linestyle"] = linestyle
+        kwargs["hist_kwargs"] = hist_kwargs
+        fig = result.plot_corner(fig=fig, save=False, color=c, contour_kwargs={"linestyle": linestyle}, **kwargs)
         default_filename += '_{}'.format(result.label)
-        lines.append(mpllines.Line2D([0], [0], color=c))
+        lines.append(mpllines.Line2D([0], [0], color=c, linestyle=linestyle))
         default_labels.append(result.label)
 
     # Rescale the axes
@@ -2160,7 +2194,7 @@ def make_pp_plot(results, filename=None, save=True, confidence_interval=[0.68, 0
 
         try:
             name = results[0].priors[key].latex_label
-        except AttributeError:
+        except (AttributeError, KeyError):
             name = key
         label = "{} ({:2.3f})".format(name, pvalue)
         plt.plot(x_values, pp, lines[ii], label=label, **kwargs)

bilby/core/sampler/base_sampler.py

@@ -18,6 +18,7 @@ from ..utils import (
     command_line_args,
     logger,
 )
+from ..utils.random import seed as set_seed
 
 
 @attr.s
@@ -305,6 +306,7 @@ class Sampler(object):
             for equiv in self.sampling_seed_equiv_kwargs:
                 if equiv in kwargs:
                     kwargs[self.sampling_seed_key] = kwargs.pop(equiv)
+                    set_seed(kwargs[self.sampling_seed_key])
         return kwargs
 
     @property
@@ -568,12 +570,14 @@ class Sampler(object):
             likelihood evaluations.
 
         """
+        from ..utils.random import rng
+
         logger.info("Generating initial points from the prior")
         unit_cube = []
         parameters = []
         likelihood = []
         while len(unit_cube) < npoints:
-            unit = np.random.rand(self.ndim)
+            unit = rng.uniform(0, 1, self.ndim)
             theta = self.prior_transform(unit)
             if self.check_draw(theta, warning=False):
                 unit_cube.append(unit)
@@ -678,11 +682,11 @@ class Sampler(object):
         if self.cached_result is None:
             kwargs_print = self.kwargs.copy()
             for k in kwargs_print:
-                if type(kwargs_print[k]) in (list, np.ndarray):
+                if isinstance(kwargs_print[k], (list, np.ndarray)):
                     array_repr = np.array(kwargs_print[k])
                     if array_repr.size > 10:
                         kwargs_print[k] = f"array_like, shape={array_repr.shape}"
-                elif type(kwargs_print[k]) == DataFrame:
+                elif isinstance(kwargs_print[k], DataFrame):
                     kwargs_print[k] = f"DataFrame, shape={kwargs_print[k].shape}"
             logger.info(
                 f"Using sampler {self.__class__.__name__} with kwargs {kwargs_print}"
@@ -898,7 +902,19 @@ class _TemporaryFileSamplerMixin:
 
     def __init__(self, temporary_directory, **kwargs):
         super(_TemporaryFileSamplerMixin, self).__init__(**kwargs)
-        self.use_temporary_directory = temporary_directory
+        try:
+            from mpi4py import MPI
+
+            using_mpi = MPI.COMM_WORLD.Get_size() > 1
+        except ImportError:
+            using_mpi = False
+
+        if using_mpi and temporary_directory:
+            logger.info(
+                "Temporary directory incompatible with MPI, "
+                "will run in original directory"
+            )
+        self.use_temporary_directory = temporary_directory and not using_mpi
         self._outputfiles_basename = None
         self._temporary_outputfiles_basename = None
 

bilby/core/sampler/cpnest.py

@@ -182,7 +182,7 @@ class Cpnest(NestedSampler):
         if "proposals" in self.kwargs:
             if self.kwargs["proposals"] is None:
                 return
-            if type(self.kwargs["proposals"]) == JumpProposalCycle:
+            if isinstance(self.kwargs["proposals"], JumpProposalCycle):
                 self.kwargs["proposals"] = dict(
                     mhs=self.kwargs["proposals"], hmc=self.kwargs["proposals"]
                 )

bilby/core/sampler/dnest4.py

@@ -2,7 +2,6 @@ import datetime
 import time
 
 import numpy as np
-import pandas as pd
 
 from ..utils import logger
 from .base_sampler import NestedSampler, _TemporaryFileSamplerMixin, signal_wrapper
@@ -42,9 +41,11 @@ class _DNest4Model(object):
 
     def perturb(self, coords):
         """The perturb function to perform Monte Carlo trial moves."""
-        idx = np.random.randint(self._n_dim)
+        from ..utils.random import rng
 
-        coords[idx] += self._widths[idx] * (np.random.uniform(size=1) - 0.5)
+        idx = rng.integers(self._n_dim)
+
+        coords[idx] += self._widths[idx] * (rng.uniform(size=1) - 0.5)
         cw = self._widths[idx]
         cc = self._centers[idx]
 
@@ -150,7 +151,6 @@ class DNest4(_TemporaryFileSamplerMixin, NestedSampler):
         self.start_time = np.nan
         self.sampler = None
         self._information = np.nan
-        self._last_live_sample_info = np.nan
 
         # Get the estimates of the prior distributions' widths and centers.
         widths = []
@@ -229,22 +229,7 @@ class DNest4(_TemporaryFileSamplerMixin, NestedSampler):
         self.result.log_evidence = stats["log_Z"]
         self._information = stats["H"]
         self.result.log_evidence_err = np.sqrt(self._information / self.num_particles)
-
-        if self._backend == "memory":
-            self._last_live_sample_info = pd.DataFrame(
-                self.sampler.backend.sample_info[-1]
-            )
-            self.result.log_likelihood_evaluations = self._last_live_sample_info[
-                "log_likelihood"
-            ]
-            self.result.samples = np.array(self.sampler.backend.posterior_samples)
-        else:
-            sample_info_path = (
-                "./" + self.kwargs["outputfiles_basename"] + "/sample_info.txt"
-            )
-            sample_info = np.genfromtxt(sample_info_path, comments="#", names=True)
-            self.result.log_likelihood_evaluations = sample_info["log_likelihood"]
-            self.result.samples = np.array(self.sampler.backend.posterior_samples)
+        self.result.samples = np.array(self.sampler.backend.posterior_samples)
 
         self.result.sampler_output = out
         self.result.outputfiles_basename = self.outputfiles_basename

bilby/core/sampler/dynesty.py

@@ -568,6 +568,8 @@ class Dynesty(NestedSampler):
         import dynesty
         from scipy.special import logsumexp
 
+        from ..utils.random import rng
+
         logwts = out["logwt"]
         weights = np.exp(logwts - out["logz"][-1])
         nested_samples = DataFrame(out.samples, columns=self.search_parameter_keys)
@@ -575,7 +577,7 @@ class Dynesty(NestedSampler):
         nested_samples["log_likelihood"] = out.logl
         self.result.nested_samples = nested_samples
         if self.rejection_sample_posterior:
-            keep = weights > np.random.uniform(0, max(weights), len(weights))
+            keep = weights > rng.uniform(0, max(weights), len(weights))
             self.result.samples = out.samples[keep]
             self.result.log_likelihood_evaluations = out.logl[keep]
             logger.info(

bilby/core/sampler/emcee.py

@@ -403,6 +403,7 @@ class Emcee(MCMCSampler):
         if self.verbose:
             iterator.close()
         self.write_current_state()
+        self._close_pool()
 
         self.result.sampler_output = np.nan
         self.calculate_autocorrelation(self.sampler.chain.reshape((-1, self.ndim)))

bilby/core/sampler/proposal.py

-import random
 from inspect import isclass
 
 import numpy as np
 
 from ..prior import Uniform
+from ..utils import random
 
 
 class Sample(dict):
@@ -135,7 +135,7 @@ class JumpProposalCycle(object):
         return len(self.proposal_functions)
 
     def update_cycle(self):
-        self._cycle = np.random.choice(
+        self._cycle = random.rng.choice(
             self.proposal_functions,
             size=self.cycle_length,
             p=self.weights,
@@ -195,14 +195,14 @@ class NormJump(JumpProposal):
 
     def __call__(self, sample, **kwargs):
         for key in sample.keys():
-            sample[key] = np.random.normal(sample[key], self.step_size)
+            sample[key] = random.rng.normal(sample[key], self.step_size)
         return super(NormJump, self).__call__(sample)
 
 
 class EnsembleWalk(JumpProposal):
     def __init__(
         self,
-        random_number_generator=random.random,
+        random_number_generator=random.rng.uniform,
         n_points=3,
         priors=None,
         **random_number_generator_args
@@ -227,7 +227,7 @@ class EnsembleWalk(JumpProposal):
         self.random_number_generator_args = random_number_generator_args
 
     def __call__(self, sample, **kwargs):
-        subset = random.sample(kwargs["coordinates"], self.n_points)
+        subset = random.rng.choice(kwargs["coordinates"], self.n_points, replace=False)
         for i in range(len(subset)):
             subset[i] = Sample.from_external_type(
                 subset[i], kwargs.get("sampler_name", None)
@@ -258,11 +258,11 @@ class EnsembleStretch(JumpProposal):
         self.scale = scale
 
     def __call__(self, sample, **kwargs):
-        second_sample = random.choice(kwargs["coordinates"])
+        second_sample = random.rng.choice(kwargs["coordinates"])
         second_sample = Sample.from_external_type(
             second_sample, kwargs.get("sampler_name", None)
         )
-        step = random.uniform(-1, 1) * np.log(self.scale)
+        step = random.rng.uniform(-1, 1) * np.log(self.scale)
         sample = second_sample + (sample - second_sample) * np.exp(step)
         self.log_j = len(sample) * step
         return super(EnsembleStretch, self).__call__(sample)
@@ -286,8 +286,8 @@ class DifferentialEvolution(JumpProposal):
         self.mu = mu
 
     def __call__(self, sample, **kwargs):
-        a, b = random.sample(kwargs["coordinates"], 2)
-        sample = sample + (b - a) * random.gauss(self.mu, self.sigma)
+        a, b = random.rng.choice(kwargs["coordinates"], 2, replace=False)
+        sample = sample + (b - a) * random.rng.normal(self.mu, self.sigma)
         return super(DifferentialEvolution, self).__call__(sample)
 
 
@@ -334,8 +334,8 @@ class EnsembleEigenVector(JumpProposal):
 
     def __call__(self, sample, **kwargs):
         self.update_eigenvectors(kwargs["coordinates"])
-        i = random.randrange(len(sample))
-        jump_size = np.sqrt(np.fabs(self.eigen_values[i])) * random.gauss(0, 1)
+        i = random.rng.integers(len(sample))
+        jump_size = np.sqrt(np.fabs(self.eigen_values[i])) * random.rng.normal(0, 1)
         for j, key in enumerate(sample.keys()):
             sample[key] += jump_size * self.eigen_vectors[j, i]
         return super(EnsembleEigenVector, self).__call__(sample)

bilby/core/sampler/ptemcee.py

@@ -323,6 +323,8 @@ class Ptemcee(MCMCSampler):
 
         from scipy.optimize import minimize
 
+        from ..utils.random import rng
+
         # Set up the minimize list: keys not in this list will have initial
         # positions drawn from the prior
         if minimize_list is None:
@@ -375,7 +377,7 @@ class Ptemcee(MCMCSampler):
             pos0_max = np.max(success[:, i])
             logger.info(f"Initialize {key} walkers from {pos0_min}->{pos0_max}")
             j = self.search_parameter_keys.index(key)
-            pos0[:, :, j] = np.random.uniform(
+            pos0[:, :, j] = rng.uniform(
                 pos0_min,
                 pos0_max,
                 size=(self.kwargs["ntemps"], self.kwargs["nwalkers"]),

bilby/core/sampler/pymultinest.py

@@ -5,7 +5,6 @@ import time
 
 import numpy as np
 
-from ..utils import logger
 from .base_sampler import NestedSampler, _TemporaryFileSamplerMixin, signal_wrapper
 
 
@@ -77,12 +76,6 @@ class Pymultinest(_TemporaryFileSamplerMixin, NestedSampler):
         temporary_directory=True,
         **kwargs
     ):
-        try:
-            from mpi4py import MPI
-
-            using_mpi = MPI.COMM_WORLD.Get_size() > 1
-        except ImportError:
-            using_mpi = False
         super(Pymultinest, self).__init__(
             likelihood=likelihood,
             priors=priors,
@@ -96,13 +89,6 @@ class Pymultinest(_TemporaryFileSamplerMixin, NestedSampler):
             **kwargs
         )
         self._apply_multinest_boundaries()
-        self.exit_code = exit_code
-        if using_mpi and temporary_directory:
-            logger.info(
-                "Temporary directory incompatible with MPI, "
-                "will run in original directory"
-            )
-        self.use_temporary_directory = temporary_directory and not using_mpi
 
     def _translate_kwargs(self, kwargs):
         kwargs = super()._translate_kwargs(kwargs)

bilby/core/sampler/ultranest.py

@@ -262,11 +262,13 @@ class Ultranest(_TemporaryFileSamplerMixin, NestedSampler):
 
     def _generate_result(self, out):
         # extract results
+        from ..utils.random import rng
+
         data = np.array(out["weighted_samples"]["points"])
         weights = np.array(out["weighted_samples"]["weights"])
 
         scaledweights = weights / weights.max()
-        mask = np.random.rand(len(scaledweights)) < scaledweights
+        mask = rng.uniform(0, 1, len(scaledweights)) < scaledweights
 
         nested_samples = DataFrame(data, columns=self.search_parameter_keys)
         nested_samples["weights"] = weights