.gitlab-ci.yml

@@ -59,9 +59,9 @@ containers:
           ${script} --help;
       done
 
-basic-3.8:
-  <<: *test-python
-  image: python:3.8
+# basic-3.8:
+#   <<: *test-python
+#   image: python:3.8
 
 basic-3.9:
   <<: *test-python
@@ -78,9 +78,9 @@ basic-3.10:
     - python -m pip list installed
     - 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.8:
+#   <<: *test-samplers-import
+#   image: containers.ligo.org/lscsoft/bilby/v2-bilby-python38
 
 import-samplers-3.9:
   <<: *test-samplers-import
@@ -102,12 +102,12 @@ 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.8:
+#   <<: *precommits
+#   image: containers.ligo.org/lscsoft/bilby/v2-bilby-python38
+#   variables:
+#     CACHE_DIR: ".pip38"
+#     PYVERSION: "python38"
 
 precommits-py3.9:
   <<: *precommits
@@ -142,10 +142,10 @@ install:
 
     - 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.8:
+#   <<: *unit-test
+#   needs: ["basic-3.8", "precommits-py3.8"]
+#   image: containers.ligo.org/lscsoft/bilby/v2-bilby-python38
 
 python-3.9:
   <<: *unit-test
@@ -176,10 +176,10 @@ python-3.10:
     - python -m pip list installed
     - 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.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
@@ -213,10 +213,10 @@ integration-tests-python-3.9:
     - 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.8:
+#   <<: *plotting
+#   image: containers.ligo.org/lscsoft/bilby/v2-bilby-python38
+#   needs: ["basic-3.8", "precommits-py3.8"]
 
 plotting-python-3.9:
   <<: *plotting
@@ -285,10 +285,10 @@ 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-python38-container:
+#   <<: *build-container
+#   variables:
+#     PYVERSION: "python38"
 
 build-python39-container:
   <<: *build-container

CHANGELOG.md

 # All notable changes will be documented in this file
 
+## [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):

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
@@ -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
 
@@ -219,10 +220,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 +265,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]
 
 
@@ -1870,7 +1876,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

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)

bilby/core/sampler/dnest4.py

@@ -42,9 +42,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]
 

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/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/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

bilby/core/utils/__init__.py

+from . import random
 from .calculus import *
 from .cmd import *
 from .colors import *

bilby/core/utils/random.py

+from numpy.random import default_rng, SeedSequence
+
+
+def __getattr__(name):
+    if name == "rng":
+        return Generator.rng
+
+
+class Generator:
+    rng = default_rng()
+
+
+def seed(seed):
+    Generator.rng = default_rng(seed)
+
+
+def generate_seeds(nseeds):
+    return SeedSequence(Generator.rng.integers(0, 2**63 - 1, size=4)).spawn(nseeds)

bilby/core/utils/series.py

 import numpy as np
 
-
 _TOL = 14
 
 
@@ -165,29 +164,23 @@ def create_white_noise(sampling_frequency, duration):
     array_like: white noise
     array_like: frequency array
     """
+    from .random import rng
 
     number_of_samples = duration * sampling_frequency
     number_of_samples = int(np.round(number_of_samples))
 
-    delta_freq = 1. / duration
-
     frequencies = create_frequency_series(sampling_frequency, duration)
 
-    norm1 = 0.5 * (1. / delta_freq)**0.5
-    re1 = np.random.normal(0, norm1, len(frequencies))
-    im1 = np.random.normal(0, norm1, len(frequencies))
-    htilde1 = re1 + 1j * im1
+    norm1 = 0.5 * duration**0.5
+    re1, im1 = rng.normal(0, norm1, (2, len(frequencies)))
+    white_noise = re1 + 1j * im1
 
-    # convolve data with instrument transfer function
-    otilde1 = htilde1 * 1.
     # set DC and Nyquist = 0
-    otilde1[0] = 0
+    white_noise[0] = 0
     # no Nyquist frequency when N=odd
     if np.mod(number_of_samples, 2) == 0:
-        otilde1[-1] = 0
+        white_noise[-1] = 0
 
-    # normalise for positive frequencies and units of strain/rHz
-    white_noise = otilde1
     # python: transpose for use with infft
     white_noise = np.transpose(white_noise)
     frequencies = np.transpose(frequencies)

bilby/gw/conversion.py

@@ -1619,8 +1619,9 @@ def binary_love_lambda_symmetric_to_lambda_1_lambda_2_automatic_marginalisation(
     lambda_2: float
         Tidal parameter of less massive neutron star.
     """
+    from ..core.utils.random import rng
 
-    binary_love_uniform = np.random.uniform(0, 1, len(lambda_symmetric))
+    binary_love_uniform = rng.uniform(0, 1, len(lambda_symmetric))
 
     lambda_1, lambda_2 = binary_love_lambda_symmetric_to_lambda_1_lambda_2_manual_marginalisation(
         binary_love_uniform, lambda_symmetric, mass_ratio)
@@ -2422,6 +2423,8 @@ def generate_posterior_samples_from_marginalized_likelihood(
     sample: DataFrame
         Returns the posterior with new samples.
     """
+    from ..core.utils.random import generate_seeds
+
     marginalized_parameters = getattr(likelihood, "_marginalized_parameters", list())
     if len(marginalized_parameters) == 0:
         return samples
@@ -2485,7 +2488,8 @@ def generate_posterior_samples_from_marginalized_likelihood(
         _sampling_convenience_dump.likelihood = likelihood
         pool = None
 
-    fill_args = [(ii, row) for ii, row in samples.iterrows()]
+    seeds = generate_seeds(len(samples))
+    fill_args = [(ii, row, seed) for (ii, row), seed in zip(samples.iterrows(), seeds)]
     ii = 0
     pbar = tqdm(total=len(samples), file=sys.stdout)
     while ii < len(samples):
@@ -2544,8 +2548,11 @@ def generate_sky_frame_parameters(samples, likelihood):
 
 def fill_sample(args):
     from ..core.sampler.base_sampler import _sampling_convenience_dump
+    from ..core.utils.random import seed
+
+    ii, sample, rseed = args
+    seed(rseed)
     likelihood = _sampling_convenience_dump.likelihood
-    ii, sample = args
     marginalized_parameters = getattr(likelihood, "_marginalized_parameters", list())
     sample = dict(sample).copy()
     likelihood.parameters.update(dict(sample).copy())