diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index be2b07e476fac16a241e2ddd0615658db1605826..f1aee7ee0c7a8ce7b5d31543feaedf33fdb60cf4 100644 --- a/.gitlab-ci.yml +++ b/.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 diff --git a/bilby/bilby_mcmc/chain.py b/bilby/bilby_mcmc/chain.py index 9aec43421e1ab4c14c315505fb432d369ca4b3e6..9a6e36717708ec1d7ff235a23f5875a9d23689d2 100644 --- a/bilby/bilby_mcmc/chain.py +++ b/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): diff --git a/bilby/bilby_mcmc/proposals.py b/bilby/bilby_mcmc/proposals.py index edba8fda2eee3a22f55857e46104217af714c56d..17892e050b33b89d8b391203b72140056ea29254 100644 --- a/bilby/bilby_mcmc/proposals.py +++ b/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): @@ -834,7 +834,7 @@ class FisherMatrixProposal(AdaptiveGaussianProposal): 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 diff --git a/bilby/bilby_mcmc/sampler.py b/bilby/bilby_mcmc/sampler.py index efa2fcb85039cc494a1c1be381b2f9bcf0fadc9a..1b528ec78cad6bd929feda0658b4b7f1d5c222d7 100644 --- a/bilby/bilby_mcmc/sampler.py +++ b/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 diff --git a/bilby/core/fisher.py b/bilby/core/fisher.py index fd452f887f36741a140ed82d53dbcadd5bc69cd5..11b9af386461d95ee7824d7c276c538d76907911 100644 --- a/bilby/core/fisher.py +++ b/bilby/core/fisher.py @@ -1,7 +1,6 @@ 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) diff --git a/bilby/core/prior/base.py b/bilby/core/prior/base.py index eef710a603321cab16f8dba196754e53b0658a9e..4ac924f74639837a2c323e5769e53c3e0a6a2871 100644 --- a/bilby/core/prior/base.py +++ b/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): diff --git a/bilby/core/prior/conditional.py b/bilby/core/prior/conditional.py index c4dcc36827fd844853095dcbbcafd69f2585c40a..cbc99f269f6dcb7acbaa75ef9ea377dc2510cbcf 100644 --- a/bilby/core/prior/conditional.py +++ b/bilby/core/prior/conditional.py @@ -1,5 +1,3 @@ -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): diff --git a/bilby/core/prior/joint.py b/bilby/core/prior/joint.py index 7f35b8716526e8811571409c2aba6d476a9958cc..cec6fefa288d48d0b71cd73ed36692e239cd44df 100644 --- a/bilby/core/prior/joint.py +++ b/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])) diff --git a/bilby/core/result.py b/bilby/core/result.py index 563ef85815fc7fbd924a2469e2b6e8dc9e1bb51d..a5e45be5cea81f5038afc8bc53451b36bcd4c763 100644 --- a/bilby/core/result.py +++ b/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 @@ -264,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] @@ -1875,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 diff --git a/bilby/core/sampler/base_sampler.py b/bilby/core/sampler/base_sampler.py index 32bde7c783a8478c578c075c2522d392e595c7b4..f4a717787219fafdf3d72edf29fb9d46fa1a5263 100644 --- a/bilby/core/sampler/base_sampler.py +++ b/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) diff --git a/bilby/core/sampler/dnest4.py b/bilby/core/sampler/dnest4.py index 5c3d7566e729fb92e5427e9c8b5968c3df6c6abe..dda360ca80211cfe6f87bd08bf8c93d885ee3e36 100644 --- a/bilby/core/sampler/dnest4.py +++ b/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] diff --git a/bilby/core/sampler/dynesty.py b/bilby/core/sampler/dynesty.py index d2e3a85b1f8cf0cd3a9b8091f980c38ce10a6ec9..2a4d73d738effde9265e0e8a450b8e5a17d4f69a 100644 --- a/bilby/core/sampler/dynesty.py +++ b/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( diff --git a/bilby/core/sampler/proposal.py b/bilby/core/sampler/proposal.py index 023caac5744de968d338cea30125574248b91f95..d23d19b4c27c5793f4cb6a074d70c805e6d5dce5 100644 --- a/bilby/core/sampler/proposal.py +++ b/bilby/core/sampler/proposal.py @@ -1,9 +1,9 @@ -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) diff --git a/bilby/core/sampler/ptemcee.py b/bilby/core/sampler/ptemcee.py index 8fce83996c37868c95708a8d227e69a2e0d1a4bc..531fb102a0614015f4053e7326b87c5ec923dbfe 100644 --- a/bilby/core/sampler/ptemcee.py +++ b/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"]), diff --git a/bilby/core/sampler/ultranest.py b/bilby/core/sampler/ultranest.py index 4cc14a9fa7ff9ba1ce4f19979571f8dde2c71e7b..542f862468c290d007f556ca6033bb30e78f2729 100644 --- a/bilby/core/sampler/ultranest.py +++ b/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 diff --git a/bilby/core/utils/__init__.py b/bilby/core/utils/__init__.py index 10509641a3e71b10620e51ba3ec98ca788b6e9d7..25d1eda934151fdcb85806dd336406a9f5ea7a90 100644 --- a/bilby/core/utils/__init__.py +++ b/bilby/core/utils/__init__.py @@ -1,3 +1,4 @@ +from . import random from .calculus import * from .cmd import * from .colors import * diff --git a/bilby/core/utils/random.py b/bilby/core/utils/random.py new file mode 100644 index 0000000000000000000000000000000000000000..dbee5b091181555e991d5d0229b147b8f31a03ac --- /dev/null +++ b/bilby/core/utils/random.py @@ -0,0 +1,18 @@ +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) diff --git a/bilby/core/utils/series.py b/bilby/core/utils/series.py index 7593d0ff8d45fde273bd56f8b53d6a553d8ea516..ccd59372245b56087d77a9b45ca0da645363f3b7 100644 --- a/bilby/core/utils/series.py +++ b/bilby/core/utils/series.py @@ -1,6 +1,5 @@ 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) diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py index 091727c5d1f6fcda39449788ce31c86265371f7a..00d957f98fea3b0eeba1459bf25f16cc47c39ab7 100644 --- a/bilby/gw/conversion.py +++ b/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()) diff --git a/bilby/gw/likelihood/base.py b/bilby/gw/likelihood/base.py index ab8eadc1d96e4ee7c7c30f47faf34a3cb5302705..256237bfb9f8395a67a7152adb66fbd18892b9d1 100644 --- a/bilby/gw/likelihood/base.py +++ b/bilby/gw/likelihood/base.py @@ -521,6 +521,8 @@ class GravitationalWaveTransient(Likelihood): new_calibration: dict Sample set from the calibration posterior """ + from ...core.utils.random import rng + if 'recalib_index' in self.parameters: self.parameters.pop('recalib_index') self.parameters.update(self.get_sky_frame_parameters()) @@ -533,7 +535,7 @@ class GravitationalWaveTransient(Likelihood): calibration_post = np.exp(log_like - max(log_like)) calibration_post /= np.sum(calibration_post) - new_calibration = np.random.choice(self.number_of_response_curves, p=calibration_post) + new_calibration = rng.choice(self.number_of_response_curves, p=calibration_post) return new_calibration diff --git a/bilby/gw/likelihood/roq.py b/bilby/gw/likelihood/roq.py index b3f93a8472200e76f460ba3623f2436ab2323183..12d5212f2230202df4410d2374c185665194c72b 100644 --- a/bilby/gw/likelihood/roq.py +++ b/bilby/gw/likelihood/roq.py @@ -1151,6 +1151,8 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient): signal[kind][mode] *= self._ref_dist / new_distance def generate_time_sample_from_marginalized_likelihood(self, signal_polarizations=None): + from ...core.utils.random import rng + self.parameters.update(self.get_sky_frame_parameters()) if signal_polarizations is None: signal_polarizations = \ @@ -1180,7 +1182,7 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient): time_prior_array = self.priors['geocent_time'].prob(times) time_post = np.exp(time_log_like - max(time_log_like)) * time_prior_array time_post /= np.sum(time_post) - return np.random.choice(times, p=time_post) + return rng.choice(times, p=time_post) class BilbyROQParamsRangeError(Exception): diff --git a/bilby/gw/prior.py b/bilby/gw/prior.py index f8b919c45f23204f04e048d486bd9d81c3128b0a..63b51790eaa5f05be5fe11f3c7b2c592f1abc32f 100644 --- a/bilby/gw/prior.py +++ b/bilby/gw/prior.py @@ -12,7 +12,7 @@ from ..core.prior import ( ConditionalPriorDict, ConditionalBasePrior, BaseJointPriorDist, JointPrior, JointPriorDistError, ) -from ..core.utils import infer_args_from_method, logger +from ..core.utils import infer_args_from_method, logger, random from .conversion import ( convert_to_lal_binary_black_hole_parameters, convert_to_lal_binary_neutron_star_parameters, generate_mass_parameters, @@ -1503,7 +1503,7 @@ class HealPixMapPriorDist(BaseJointPriorDist): """ pixel_choices = np.arange(self.npix) pixel_probs = self._check_norm(self.prob) - sample_pix = np.random.choice(pixel_choices, size=size, p=pixel_probs, replace=True) + sample_pix = random.rng.choice(pixel_choices, size=size, p=pixel_probs, replace=True) sample = np.empty((size, self.num_vars)) for samp in range(size): theta, ra = self.hp.pix2ang(self.nside, sample_pix[samp]) @@ -1535,7 +1535,7 @@ class HealPixMapPriorDist(BaseJointPriorDist): """ if self.distmu[pix] == np.inf or self.distmu[pix] <= 0: return 0 - dist = self.distance_icdf(np.random.uniform(0, 1)) + dist = self.distance_icdf(random.rng.uniform(0, 1)) name = self.names[-1] if (dist > self.bounds[name][1]) | (dist < self.bounds[name][0]): self.draw_distance(pix) @@ -1564,8 +1564,8 @@ class HealPixMapPriorDist(BaseJointPriorDist): self.draw_from_pixel(ra, dec, pix) return np.array( [ - np.random.uniform(ra - self.pixel_length, ra + self.pixel_length), - np.random.uniform(dec - self.pixel_length, dec + self.pixel_length), + random.rng.uniform(ra - self.pixel_length, ra + self.pixel_length), + random.rng.uniform(dec - self.pixel_length, dec + self.pixel_length), ] ) diff --git a/examples/gw_examples/injection_examples/fast_tutorial.py b/examples/gw_examples/injection_examples/fast_tutorial.py index f462008a70994982d47c9cbd69179161c96aae3c..9503fb04d9250f7a1f1d2f42af2ecc10846fab3d 100644 --- a/examples/gw_examples/injection_examples/fast_tutorial.py +++ b/examples/gw_examples/injection_examples/fast_tutorial.py @@ -9,7 +9,6 @@ between luminosity distances of 100Mpc and 5Gpc, the cosmology is Planck15. """ import bilby -import numpy as np # Set the duration and sampling frequency of the data segment that we're # going to inject the signal into @@ -23,7 +22,7 @@ label = "fast_tutorial" bilby.core.utils.setup_logger(outdir=outdir, label=label) # Set up a random seed for result reproducibility. This is optional! -np.random.seed(88170235) +bilby.core.utils.random.seed(88170235) # We are going to inject a binary black hole waveform. We first establish a # dictionary of parameters that includes all of the different waveform diff --git a/test/core/grid_test.py b/test/core/grid_test.py index 45cd0e1afb6a9e5172995b7cb673e2e87cbbd235..82e44c5cc656bf9e96e3cbc2ea71d0ce441e49ad 100644 --- a/test/core/grid_test.py +++ b/test/core/grid_test.py @@ -27,7 +27,7 @@ class MultiGaussian(bilby.Likelihood): class TestGrid(unittest.TestCase): def setUp(self): - np.random.seed(7) + bilby.core.utils.random.seed(7) # set 2D multivariate Gaussian (zero mean, unit variance) self.mus = [0.0, 0.0] diff --git a/test/core/prior/conditional_test.py b/test/core/prior/conditional_test.py index 94a869936fc6cd71cc3647e51a5a372e08ed4544..5ee4efd60618f89547ec4731407ca97aebf062a7 100644 --- a/test/core/prior/conditional_test.py +++ b/test/core/prior/conditional_test.py @@ -275,20 +275,22 @@ class TestConditionalPriorDict(unittest.TestCase): self.conditional_priors.ln_prob(sample=self.test_sample) def test_sample_subset_all_keys(self): - with mock.patch("numpy.random.uniform") as m: - m.return_value = 0.5 - self.assertDictEqual( - dict(var_0=0.5, var_1=0.5 ** 2, var_2=0.5 ** 3, var_3=0.5 ** 4), - self.conditional_priors.sample_subset( - keys=["var_0", "var_1", "var_2", "var_3"] - ), - ) + bilby.core.utils.random.seed(5) + self.assertDictEqual( + dict( + var_0=0.8050029237453802, + var_1=0.6503946979510289, + var_2=0.33516501262044845, + var_3=0.09579062316418356, + ), + self.conditional_priors.sample_subset( + keys=["var_0", "var_1", "var_2", "var_3"] + ), + ) def test_sample_illegal_subset(self): - with mock.patch("numpy.random.uniform") as m: - m.return_value = 0.5 - with self.assertRaises(bilby.core.prior.IllegalConditionsException): - self.conditional_priors.sample_subset(keys=["var_1"]) + with self.assertRaises(bilby.core.prior.IllegalConditionsException): + self.conditional_priors.sample_subset(keys=["var_1"]) def test_sample_multiple(self): def condition_func(reference_params, a): diff --git a/test/core/prior/dict_test.py b/test/core/prior/dict_test.py index f2045ed227eddbc63d2038bc806a4ea442b63c08..7945b5fba5682a52d1ba1b54769fd60714e38e14 100644 --- a/test/core/prior/dict_test.py +++ b/test/core/prior/dict_test.py @@ -258,11 +258,11 @@ class TestPriorDict(unittest.TestCase): def test_sample(self): size = 7 - np.random.seed(42) + bilby.core.utils.random.seed(42) samples1 = self.prior_set_from_dict.sample_subset( keys=self.prior_set_from_dict.keys(), size=size ) - np.random.seed(42) + bilby.core.utils.random.seed(42) samples2 = self.prior_set_from_dict.sample(size=size) self.assertEqual(set(samples1.keys()), set(samples2.keys())) for key in samples1: @@ -305,12 +305,12 @@ class TestPriorDict(unittest.TestCase): Note that the format of inputs/outputs is different between the two methods. """ sample = self.prior_set_from_dict.sample() - self.assertEqual( - self.prior_set_from_dict.rescale( - sample.keys(), - self.prior_set_from_dict.cdf(sample=sample).values() - ), list(sample.values()) - ) + original = np.array(list(sample.values())) + new = np.array(self.prior_set_from_dict.rescale( + sample.keys(), + self.prior_set_from_dict.cdf(sample=sample).values() + )) + self.assertLess(max(abs(original - new)), 1e-10) def test_redundancy(self): for key in self.prior_set_from_dict.keys(): diff --git a/test/core/sampler/proposal_test.py b/test/core/sampler/proposal_test.py index ce678d3fe1d1fdaa5a2750c664a3e6d7a284377f..1b2331af0febf8cebe2e730b1e22eb3a66062201 100644 --- a/test/core/sampler/proposal_test.py +++ b/test/core/sampler/proposal_test.py @@ -1,9 +1,9 @@ import unittest from unittest import mock -import random import numpy as np +import bilby from bilby.core import prior from bilby.core.sampler import proposal @@ -132,6 +132,7 @@ class TestNormJump(unittest.TestCase): ) ) self.jump_proposal = proposal.NormJump(step_size=3.0, priors=self.priors) + bilby.core.utils.random.seed(5) def tearDown(self): del self.priors @@ -145,12 +146,14 @@ class TestNormJump(unittest.TestCase): self.assertEqual(1.0, self.jump_proposal.step_size) def test_jump_proposal_call(self): - with mock.patch("numpy.random.normal") as m: - m.return_value = 0.5 - sample = proposal.Sample(dict(reflective=0.0, periodic=0.0, default=0.0)) - new_sample = self.jump_proposal(sample) - expected = proposal.Sample(dict(reflective=0.5, periodic=0.5, default=0.5)) - self.assertDictEqual(expected, new_sample) + sample = proposal.Sample(dict(reflective=0.0, periodic=0.0, default=0.0)) + new_sample = self.jump_proposal(sample) + expected = proposal.Sample(dict( + reflective=0.5942057242396577, + periodic=-0.02692301311556511, + default=-0.7450848662857457, + )) + self.assertDictEqual(expected, new_sample) class TestEnsembleWalk(unittest.TestCase): @@ -164,9 +167,11 @@ class TestEnsembleWalk(unittest.TestCase): default=prior.Uniform(minimum=-0.5, maximum=1), ) ) + bilby.core.utils.random.seed(5) self.jump_proposal = proposal.EnsembleWalk( - random_number_generator=random.random, n_points=4, priors=self.priors + random_number_generator=bilby.core.utils.random.rng.uniform, n_points=4, priors=self.priors ) + self.coordinates = [proposal.Sample(self.priors.sample()) for _ in range(10)] def tearDown(self): del self.priors @@ -180,7 +185,7 @@ class TestEnsembleWalk(unittest.TestCase): self.assertEqual(3, self.jump_proposal.n_points) def test_random_number_generator_init(self): - self.assertEqual(random.random, self.jump_proposal.random_number_generator) + self.assertEqual(bilby.core.utils.random.rng.uniform, self.jump_proposal.random_number_generator) def test_get_center_of_mass(self): samples = [ @@ -193,17 +198,15 @@ class TestEnsembleWalk(unittest.TestCase): self.assertAlmostEqual(expected[key], actual[key]) def test_jump_proposal_call(self): - with mock.patch("random.sample") as m: - self.jump_proposal.random_number_generator = lambda: 2 - m.return_value = [ - proposal.Sample(dict(periodic=0.3, reflective=0.3, default=0.3)), - proposal.Sample(dict(periodic=0.1, reflective=0.1, default=0.1)), - ] - sample = proposal.Sample(dict(periodic=0.1, reflective=0.1, default=0.1)) - new_sample = self.jump_proposal(sample, coordinates=None) - expected = proposal.Sample(dict(periodic=0.1, reflective=0.1, default=0.1)) - for key, value in new_sample.items(): - self.assertAlmostEqual(expected[key], value) + sample = proposal.Sample(dict(periodic=0.1, reflective=0.1, default=0.1)) + new_sample = self.jump_proposal(sample, coordinates=self.coordinates) + expected = proposal.Sample(dict( + periodic=0.437075089594473, + reflective=-0.18027731528487945, + default=-0.17570046901727415, + )) + for key, value in new_sample.items(): + self.assertAlmostEqual(expected[key], value) class TestEnsembleEnsembleStretch(unittest.TestCase): @@ -217,7 +220,9 @@ class TestEnsembleEnsembleStretch(unittest.TestCase): default=prior.Uniform(minimum=-0.5, maximum=1), ) ) + bilby.core.utils.random.seed(5) self.jump_proposal = proposal.EnsembleStretch(scale=3.0, priors=self.priors) + self.coordinates = [proposal.Sample(self.priors.sample()) for _ in range(10)] def tearDown(self): del self.priors @@ -231,47 +236,24 @@ class TestEnsembleEnsembleStretch(unittest.TestCase): self.assertEqual(5.0, self.jump_proposal.scale) def test_jump_proposal_call(self): - with mock.patch("random.choice") as m: - with mock.patch("random.uniform") as n: - second_sample = proposal.Sample( - dict(periodic=0.3, reflective=0.3, default=0.3) - ) - random_number = 0.5 - m.return_value = second_sample - n.return_value = random_number - sample = proposal.Sample( - dict(periodic=0.1, reflective=0.1, default=0.1) - ) - new_sample = self.jump_proposal(sample, coordinates=None) - coords = 0.3 - 0.2 * np.exp( - random_number * np.log(self.jump_proposal.scale) - ) - expected = proposal.Sample( - dict(periodic=coords, reflective=coords, default=coords) - ) - for key, value in new_sample.items(): - self.assertAlmostEqual(expected[key], value) + sample = proposal.Sample( + dict(periodic=0.1, reflective=0.1, default=0.1) + ) + new_sample = self.jump_proposal(sample, coordinates=self.coordinates) + expected = proposal.Sample(dict( + periodic=0.5790181653312239, + reflective=-0.028378746842481914, + default=-0.23534241783479043, + )) + for key, value in new_sample.items(): + self.assertAlmostEqual(expected[key], value) def test_log_j_after_call(self): - with mock.patch("random.uniform") as m1: - with mock.patch("numpy.log") as m2: - with mock.patch("numpy.exp") as m3: - m1.return_value = 1 - m2.return_value = 1 - m3.return_value = 1 - coordinates = [ - proposal.Sample( - dict(periodic=0.3, reflective=0.3, default=0.3) - ), - proposal.Sample( - dict(periodic=0.3, reflective=0.3, default=0.3) - ), - ] - sample = proposal.Sample( - dict(periodic=0.2, reflective=0.2, default=0.2) - ) - self.jump_proposal(sample=sample, coordinates=coordinates) - self.assertEqual(3, self.jump_proposal.log_j) + sample = proposal.Sample( + dict(periodic=0.2, reflective=0.2, default=0.2) + ) + self.jump_proposal(sample=sample, coordinates=self.coordinates) + self.assertAlmostEqual(-3.2879289432183088, self.jump_proposal.log_j, 10) class TestDifferentialEvolution(unittest.TestCase): @@ -285,9 +267,11 @@ class TestDifferentialEvolution(unittest.TestCase): default=prior.Uniform(minimum=-0.5, maximum=1), ) ) + bilby.core.utils.random.seed(5) self.jump_proposal = proposal.DifferentialEvolution( sigma=1e-3, mu=0.5, priors=self.priors ) + self.coordinates = [proposal.Sample(self.priors.sample()) for _ in range(10)] def tearDown(self): del self.priors @@ -305,22 +289,17 @@ class TestDifferentialEvolution(unittest.TestCase): self.assertEqual(2, self.jump_proposal.sigma) def test_jump_proposal_call(self): - with mock.patch("random.sample") as m: - with mock.patch("random.gauss") as n: - m.return_value = ( - proposal.Sample(dict(periodic=0.2, reflective=0.2, default=0.2)), - proposal.Sample(dict(periodic=0.3, reflective=0.3, default=0.3)), - ) - n.return_value = 1 - sample = proposal.Sample( - dict(periodic=0.1, reflective=0.1, default=0.1) - ) - expected = proposal.Sample( - dict(periodic=0.2, reflective=0.2, default=0.2) - ) - new_sample = self.jump_proposal(sample, coordinates=None) - for key, value in new_sample.items(): - self.assertAlmostEqual(expected[key], value) + sample = proposal.Sample( + dict(periodic=0.1, reflective=0.1, default=0.1) + ) + expected = proposal.Sample(dict( + periodic=0.09440864471444077, + reflective=0.567962015300636, + default=0.0657296821780595, + )) + new_sample = self.jump_proposal(sample, coordinates=self.coordinates) + for key, value in new_sample.items(): + self.assertAlmostEqual(expected[key], value) class TestEnsembleEigenVector(unittest.TestCase): @@ -334,7 +313,9 @@ class TestEnsembleEigenVector(unittest.TestCase): default=prior.Uniform(minimum=-0.5, maximum=1), ) ) + bilby.core.utils.random.seed(5) self.jump_proposal = proposal.EnsembleEigenVector(priors=self.priors) + self.coordinates = [proposal.Sample(self.priors.sample()) for _ in range(10)] def tearDown(self): del self.priors @@ -386,26 +367,17 @@ class TestEnsembleEigenVector(unittest.TestCase): self.assertEqual(2, self.jump_proposal.eigen_vectors) def test_jump_proposal_call(self): - self.jump_proposal.update_eigenvectors = lambda x: None - self.jump_proposal.eigen_values = np.array([1, np.nan, np.nan]) - self.jump_proposal.eigen_vectors = np.array( - [[0.1, np.nan, np.nan], [0.4, np.nan, np.nan], [0.7, np.nan, np.nan]] - ) - with mock.patch("random.randrange") as m: - with mock.patch("random.gauss") as n: - m.return_value = 0 - n.return_value = 1 - expected = proposal.Sample() - expected["periodic"] = 0.2 - expected["reflective"] = 0.5 - expected["default"] = 0.8 - sample = proposal.Sample() - sample["periodic"] = 0.1 - sample["reflective"] = 0.1 - sample["default"] = 0.1 - new_sample = self.jump_proposal(sample, coordinates=None) - for key, value in new_sample.items(): - self.assertAlmostEqual(expected[key], value) + expected = proposal.Sample() + expected["periodic"] = 0.10318172002873117 + expected["reflective"] = 0.11177972036165257 + expected["default"] = 0.10053457100669783 + sample = proposal.Sample() + sample["periodic"] = 0.1 + sample["reflective"] = 0.1 + sample["default"] = 0.1 + new_sample = self.jump_proposal(sample, coordinates=self.coordinates) + for key, value in new_sample.items(): + self.assertAlmostEqual(expected[key], value) if __name__ == "__main__": diff --git a/test/gw/likelihood/marginalization_test.py b/test/gw/likelihood/marginalization_test.py index 42fa0332fce6c3b1d0f79baa089d5013a16a2729..d7b4f8a59785e695b5b6d58cb5f036c60b4e5101 100644 --- a/test/gw/likelihood/marginalization_test.py +++ b/test/gw/likelihood/marginalization_test.py @@ -14,7 +14,7 @@ from scipy.special import logsumexp class TestMarginalizedLikelihood(unittest.TestCase): def setUp(self): - np.random.seed(500) + bilby.core.utils.random.seed(500) self.duration = 4 self.sampling_frequency = 2048 self.parameters = dict( @@ -177,7 +177,7 @@ class TestMarginalizations(unittest.TestCase): path_to_roq_weights = "weights.npz" def setUp(self): - np.random.seed(500) + bilby.core.utils.random.seed(200) self.duration = 4 self.sampling_frequency = 2048 self.parameters = dict( @@ -215,7 +215,7 @@ class TestMarginalizations(unittest.TestCase): waveform_arguments=dict( reference_frequency=20.0, minimum_frequency=20.0, - approximant="IMRPhenomPv2", + waveform_approximant="IMRPhenomPv2", ) ) self.interferometers.inject_signal( @@ -249,7 +249,7 @@ class TestMarginalizations(unittest.TestCase): waveform_arguments=dict( reference_frequency=20.0, minimum_frequency=20.0, - approximant="IMRPhenomPv2", + waveform_approximant="IMRPhenomPv2", frequency_nodes_linear=np.load(f"{roq_dir}/fnodes_linear.npy"), frequency_nodes_quadratic=np.load(f"{roq_dir}/fnodes_quadratic.npy"), ) @@ -265,7 +265,7 @@ class TestMarginalizations(unittest.TestCase): waveform_arguments=dict( reference_frequency=20.0, minimum_frequency=20.0, - approximant="IMRPhenomPv2", + waveform_approximant="IMRPhenomPv2", ) ) @@ -333,7 +333,7 @@ class TestMarginalizations(unittest.TestCase): cls_ = bilby.gw.likelihood.ROQGravitationalWaveTransient elif kind == "relbin": cls_ = bilby.gw.likelihood.RelativeBinningGravitationalWaveTransient - kwargs["epsilon"] = 0.3 + kwargs["epsilon"] = 0.1 self.parameters["fiducial"] = 0 else: raise ValueError(f"kind {kind} not understood") @@ -366,7 +366,15 @@ class TestMarginalizations(unittest.TestCase): marg_like, marginalized.log_likelihood_ratio(), delta=0.5 ) - @parameterized.expand(_parameters) + @parameterized.expand( + _parameters, + name_func=lambda func, num, param: ( + f"{func.__name__}_{num}__{param.args[0]}_{param.args[1]}_" + "_".join([ + ["D", "T", "P"][ii] for ii, val + in enumerate(param.args[-3:]) if val + ]) + ) + ) def test_marginalisation(self, kind, key, distance, time, phase): if all([distance, time, phase]): pytest.skip() diff --git a/test/gw/likelihood_test.py b/test/gw/likelihood_test.py index 16b9efab9d035611a2d4d3caacd4038d4662f0e5..87a007072f7dff2f0ee00976356eef43a3d5d15b 100644 --- a/test/gw/likelihood_test.py +++ b/test/gw/likelihood_test.py @@ -12,7 +12,7 @@ from bilby.gw.likelihood import BilbyROQParamsRangeError class TestBasicGWTransient(unittest.TestCase): def setUp(self): - np.random.seed(500) + bilby.core.utils.random.seed(500) self.parameters = dict( mass_1=31.0, mass_2=29.0, @@ -56,13 +56,13 @@ class TestBasicGWTransient(unittest.TestCase): """Test noise log likelihood matches precomputed value""" self.likelihood.noise_log_likelihood() self.assertAlmostEqual( - -4037.0994372143414, self.likelihood.noise_log_likelihood(), 3 + -4014.1787704539474, self.likelihood.noise_log_likelihood(), 3 ) def test_log_likelihood(self): """Test log likelihood matches precomputed value""" self.likelihood.log_likelihood() - self.assertAlmostEqual(self.likelihood.log_likelihood(), -4054.047229508672, 3) + self.assertAlmostEqual(self.likelihood.log_likelihood(), -4032.4397343470005, 3) def test_log_likelihood_ratio(self): """Test log likelihood ratio returns the correct value""" @@ -87,7 +87,7 @@ class TestBasicGWTransient(unittest.TestCase): class TestGWTransient(unittest.TestCase): def setUp(self): - np.random.seed(500) + bilby.core.utils.random.seed(500) self.duration = 4 self.sampling_frequency = 2048 self.parameters = dict( @@ -141,14 +141,14 @@ class TestGWTransient(unittest.TestCase): """Test noise log likelihood matches precomputed value""" self.likelihood.noise_log_likelihood() self.assertAlmostEqual( - -4037.0994372143414, self.likelihood.noise_log_likelihood(), 3 + -4014.1787704539474, self.likelihood.noise_log_likelihood(), 3 ) def test_log_likelihood(self): """Test log likelihood matches precomputed value""" self.likelihood.log_likelihood() self.assertAlmostEqual(self.likelihood.log_likelihood(), - -4054.047229508673, 3) + -4032.4397343470005, 3) def test_log_likelihood_ratio(self): """Test log likelihood ratio returns the correct value""" @@ -760,7 +760,8 @@ class TestROQLikelihoodHDF5(unittest.TestCase): if add_cal_errors: spline_calibration_nodes = 10 - np.random.seed(170817) + bilby.core.utils.random.seed(170817) + rng = bilby.core.utils.random.rng for ifo in interferometers: prefix = f"recalib_{ifo.name}_" ifo.calibration_model = bilby.gw.calibration.CubicSpline( @@ -772,9 +773,9 @@ class TestROQLikelihoodHDF5(unittest.TestCase): for i in range(spline_calibration_nodes): # 5% in amplitude, 5deg in phase self.injection_parameters[f"{prefix}amplitude_{i}"] = \ - np.random.normal(loc=0, scale=0.05) + rng.normal(loc=0, scale=0.05) self.injection_parameters[f"{prefix}phase_{i}"] = \ - np.random.normal(loc=0, scale=5 * np.pi / 180) + rng.normal(loc=0, scale=5 * np.pi / 180) waveform_generator = bilby.gw.WaveformGenerator( duration=self.duration, @@ -1171,7 +1172,8 @@ class TestMBLikelihood(unittest.TestCase): ) # Network SNR is ~50 self.ifos = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"]) - np.random.seed(170817) + bilby.core.utils.random.seed(70817) + rng = bilby.core.utils.random.rng self.ifos.set_strain_data_from_power_spectral_densities( sampling_frequency=self.sampling_frequency, duration=self.duration, start_time=self.test_parameters['geocent_time'] - self.duration + 2. @@ -1193,9 +1195,9 @@ class TestMBLikelihood(unittest.TestCase): self.test_parameters[f"recalib_{ifo.name}_phase_{i}"] = 0 # Calibration errors of 5% in amplitude and 5 degrees in phase self.calibration_parameters[f"recalib_{ifo.name}_amplitude_{i}"] = \ - np.random.normal(loc=0, scale=0.05) + rng.normal(loc=0, scale=0.05) self.calibration_parameters[f"recalib_{ifo.name}_phase_{i}"] = \ - np.random.normal(loc=0, scale=5 * np.pi / 180) + rng.normal(loc=0, scale=5 * np.pi / 180) self.priors = bilby.gw.prior.BBHPriorDict() self.priors.pop("mass_1") @@ -1230,7 +1232,7 @@ class TestMBLikelihood(unittest.TestCase): duration=self.duration, sampling_frequency=self.sampling_frequency, frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole, waveform_arguments=dict( - reference_frequency=self.fmin, approximant=approximant + reference_frequency=self.fmin, waveform_approximant=approximant ) ) self.ifos.inject_signal(parameters=self.test_parameters, waveform_generator=wfg) @@ -1239,7 +1241,7 @@ class TestMBLikelihood(unittest.TestCase): duration=self.duration, sampling_frequency=self.sampling_frequency, frequency_domain_source_model=bilby.gw.source.binary_black_hole_frequency_sequence, waveform_arguments=dict( - reference_frequency=self.fmin, approximant=approximant + reference_frequency=self.fmin, waveform_approximant=approximant ) ) likelihood = bilby.gw.likelihood.GravitationalWaveTransient( @@ -1270,7 +1272,7 @@ class TestMBLikelihood(unittest.TestCase): duration=self.duration, sampling_frequency=self.sampling_frequency, frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole, waveform_arguments=dict( - reference_frequency=self.fmin, approximant=approximant + reference_frequency=self.fmin, waveform_approximant=approximant ) ) self.ifos.inject_signal(parameters=self.test_parameters, waveform_generator=wfg) @@ -1279,7 +1281,7 @@ class TestMBLikelihood(unittest.TestCase): duration=self.duration, sampling_frequency=self.sampling_frequency, frequency_domain_source_model=bilby.gw.source.binary_black_hole_frequency_sequence, waveform_arguments=dict( - reference_frequency=self.fmin, approximant=approximant + reference_frequency=self.fmin, waveform_approximant=approximant ) ) likelihood = bilby.gw.likelihood.GravitationalWaveTransient( diff --git a/test/integration/make_standard_data.py b/test/integration/make_standard_data.py index 85266173496d44e875ea576013afc109db70db84..4f5d4073d44f6ce049477d56beda59524bca739d 100644 --- a/test/integration/make_standard_data.py +++ b/test/integration/make_standard_data.py @@ -5,7 +5,7 @@ import numpy as np import bilby from bilby.gw.waveform_generator import WaveformGenerator -np.random.seed(10) +bilby.core.utils.random.seed(10) time_duration = 4.0 sampling_frequency = 4096.0 diff --git a/test/integration/sample_from_the_prior_test.py b/test/integration/sample_from_the_prior_test.py index 64bb0bed105f4b983632eef689903830ccec6411..7fce8b1621109253b12165885f7bbb7d4952264a 100644 --- a/test/integration/sample_from_the_prior_test.py +++ b/test/integration/sample_from_the_prior_test.py @@ -4,7 +4,6 @@ import logging from packaging import version import unittest -import numpy as np import bilby import scipy from scipy.stats import ks_2samp, kstest @@ -40,7 +39,7 @@ class Test(unittest.TestCase): duration = 4.0 sampling_frequency = 2048.0 label = "full_15_parameters" - np.random.seed(8817021) + bilby.core.utils.random.seed(8817021) waveform_arguments = dict( waveform_approximant="IMRPhenomPv2", diff --git a/test/integration/sampler_run_test.py b/test/integration/sampler_run_test.py index adfd9f3e1149097ef684acd7f66e2a3054aef035..fa682565453c112784bf78a2edb044de17a978a7 100644 --- a/test/integration/sampler_run_test.py +++ b/test/integration/sampler_run_test.py @@ -79,12 +79,13 @@ def model(x, m, c): class TestRunningSamplers(unittest.TestCase): def setUp(self): - np.random.seed(42) + bilby.core.utils.random.seed(42) bilby.core.utils.command_line_args.bilby_test_mode = False + rng = bilby.core.utils.random.rng self.x = np.linspace(0, 1, 11) self.injection_parameters = dict(m=0.5, c=0.2) self.sigma = 0.1 - self.y = model(self.x, **self.injection_parameters) + np.random.normal( + self.y = model(self.x, **self.injection_parameters) + rng.normal( 0, self.sigma, len(self.x) ) self.likelihood = bilby.likelihood.GaussianLikelihood(