EXMP: update RNG

examples/corner_plot_example.py

@@ -24,7 +24,8 @@
 truth = {"mu": 1.7, "sigma": 0.7}
 bounds = {"mu": [-3, 3], "sigma": [0.01, 3]}
 n_points = 1000
-data = np.random.normal(truth["mu"], truth["sigma"], size=n_points)
+rng = np.random.default_rng(1234)
+data = rng.normal(truth["mu"], truth["sigma"], size=n_points)
 
 
 class GaussianLikelihood(Model):

examples/discrete_parameter.py

@@ -13,6 +13,9 @@
 output = "./outdir/discrete_example/"
 logger = setup_logger(output=output, log_level="INFO")
 
+# Set the random number generator
+rng = np.random.default_rng(1234)
+
 
 # Define the model class, this model has two signal models and a discrete
 # variable that determines which is used for the computing the log-likelihood
@@ -32,9 +35,9 @@ def new_point(self, N=1):
         # Redefine the new point method so that it correctly samples the
         # discrete parameter
         x = empty_structured_array(N, self.names)
-        x["m"] = np.random.uniform(*self.bounds["m"], size=N)
-        x["c"] = np.random.uniform(*self.bounds["c"], size=N)
-        x["model"] = np.random.choice([0, 1], size=N)
+        x["m"] = rng.uniform(*self.bounds["m"], size=N)
+        x["c"] = rng.uniform(*self.bounds["c"], size=N)
+        x["model"] = rng.choice([0, 1], size=N)
         return x
 
     def new_point_log_prob(self, x):
@@ -77,7 +80,7 @@ def log_likelihood(self, x):
 # Generate some fake data.
 # We use the straight line model for this example.
 x_data = np.linspace(0, 10, 100)
-y_data = 2.5 * x_data + 1.4 + np.random.randn(len(x_data))
+y_data = 2.5 * x_data + 1.4 + rng.standard_normal(len(x_data))
 
 # Define the sampler and specify the 'dequantise' reparameterisation for the
 # discrete parameter.

examples/gw/basic_gw_example.py

@@ -8,7 +8,6 @@
 """
 
 import bilby
-import numpy as np
 
 outdir = "./outdir/"
 label = "basic_gw_example"
@@ -18,7 +17,7 @@
 duration = 4.0
 sampling_frequency = 2048.0
 
-np.random.seed(170817)
+bilby.core.utils.random.seed(170817)
 
 # Use an injection that is similar to GW150914
 injection_parameters = dict(

examples/gw/calibration_example.py

@@ -9,7 +9,6 @@
 """
 
 import bilby
-import numpy as np
 
 # Standard configuration using bilby
 duration = 4.0
@@ -19,7 +18,7 @@
 label = "calibration_example"
 bilby.core.utils.setup_logger(outdir=outdir, label=label)
 
-np.random.seed(150914)
+bilby.core.utils.random.seed(150914)
 
 # Injection parameters for a GW150914-like BBH.
 injection_parameters = dict(

examples/gw/full_gw_example.py

@@ -9,7 +9,6 @@
 """
 
 import bilby
-import numpy as np
 
 outdir = "./outdir/"
 label = "full_gw_example"
@@ -19,7 +18,7 @@
 duration = 4.0
 sampling_frequency = 2048.0
 
-np.random.seed(151226)
+bilby.core.utils.random.seed(151226)
 
 # Use an injection that is similar to GW150914
 injection_parameters = dict(

examples/gw/ins_gw_example.py

@@ -0,0 +1,141 @@
+#!/usr/bin/env python
+
+"""
+Example of running nessai with bilby on a gravitational wave likelihood. This
+examples includes all 15 parameters for CBC and should take around 2 hours to
+run.
+
+Based on the Bilby example: https://git.ligo.org/lscsoft/bilby
+"""
+
+import bilby
+
+outdir = "./outdir/"
+label = "ins_gw_example_full"
+
+bilby.core.utils.setup_logger(outdir=outdir, label=label)
+
+duration = 4.0
+sampling_frequency = 2048.0
+
+bilby.core.utils.random.seed(151226)
+
+# Use an injection that is similar to GW150914
+injection_parameters = dict(
+    total_mass=66.0,
+    mass_ratio=0.9,
+    a_1=0.4,
+    a_2=0.3,
+    tilt_1=0.5,
+    tilt_2=1.0,
+    phi_12=1.7,
+    phi_jl=0.3,
+    luminosity_distance=2000,
+    theta_jn=0.4,
+    psi=2.659,
+    phase=1.3,
+    geocent_time=1126259642.413,
+    ra=1.375,
+    dec=-1.2108,
+)
+
+waveform_arguments = dict(
+    waveform_approximant="IMRPhenomPv2", reference_frequency=50.0
+)
+
+# Create the waveform_generator
+waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
+    sampling_frequency=sampling_frequency,
+    duration=duration,
+    frequency_domain_source_model=bilby.gw.source.lal_binary_black_hole,
+    parameter_conversion=(
+        bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters
+    ),
+    waveform_arguments=waveform_arguments,
+)
+
+# Set up interferometers
+ifos = bilby.gw.detector.InterferometerList(["H1", "L1", "V1"])
+ifos.set_strain_data_from_power_spectral_densities(
+    sampling_frequency=sampling_frequency,
+    duration=duration,
+    start_time=injection_parameters["geocent_time"] - 3,
+)
+ifos.inject_signal(
+    waveform_generator=waveform_generator, parameters=injection_parameters
+)
+
+# Set up prior
+priors = bilby.gw.prior.BBHPriorDict()
+priors["geocent_time"] = bilby.core.prior.Uniform(
+    minimum=injection_parameters["geocent_time"] - 0.1,
+    maximum=injection_parameters["geocent_time"] + 0.1,
+    name="geocent_time",
+    latex_label="$t_c$",
+    unit="$s$",
+)
+
+# Initialise the likelihood
+# nessai supports the marginalisation included in bilby
+likelihood = bilby.gw.likelihood.GravitationalWaveTransient(
+    interferometers=ifos,
+    waveform_generator=waveform_generator,
+    priors=priors,
+    phase_marginalization=True,
+    distance_marginalization=False,
+)
+priors["chirp_mass"].maximum = 40
+
+# for key in [
+#     "a_1",
+#     "a_2",
+#     "tilt_1",
+#     "tilt_2",
+#     "phi_12",
+#     "phi_jl",
+#     # "luminosity_distance",
+#     "psi",
+#     "geocent_time",
+#     "ra",
+#     "dec",
+# ]:
+#     priors[key] = injection_parameters[key]
+
+
+# Run sampler
+
+# The `flow_class` should be set to `GWFlowProposal` for GW PE. This includes
+# specific default reparameterisations for certain parameters. For example,
+# it knows that theta_jn is angle with a sine prior.
+
+result = bilby.core.sampler.run_sampler(
+    likelihood=likelihood,
+    priors=priors,
+    outdir=outdir,
+    injection_parameters=injection_parameters,
+    label=label,
+    conversion_function=bilby.gw.conversion.generate_all_bbh_parameters,
+    sampler="inessai",
+    resume=False,
+    plot=True,
+    nlive=8000,
+    min_samples=1000,
+    seed=150914,
+    flow_config=dict(
+        n_blocks=6,
+        n_neurons=32,
+        batch_norm_between_layers=True,
+    ),
+    reset_flow=4,
+    n_pool=8,
+    threshold_kwargs=dict(q=0.66),
+    draw_iid_live=True,
+    rebalance_interval=1,
+    stopping_criterion=["ratio", "fractional_error"],
+    check_criteria="all",
+    tolerance=[-1, 0.1],
+    min_iteration=5,
+)
+
+# Produce corner plots
+result.plot_corner()

examples/parallelisation_example.py

@@ -24,7 +24,8 @@
 truth = {"mu": 1.7, "sigma": 0.7}
 bounds = {"mu": [-3, 3], "sigma": [0.01, 3]}
 n_points = 1000
-data = np.random.normal(truth["mu"], truth["sigma"], size=n_points)
+rng = np.random.default_rng(1234)
+data = rng.normal(truth["mu"], truth["sigma"], size=n_points)
 
 
 class GaussianLikelihood(Model):

examples/unbounded_prior.py

@@ -12,6 +12,7 @@
 
 output = "./outdir/unbounded_prior/"
 logger = setup_logger(output=output)
+rng = np.random.default_rng(1234)
 
 # We define the model as usual but also need to redefine `new_point` since by
 # default this tries to draw from within the prior bounds which will fail if
@@ -53,10 +54,8 @@ def new_point(self, N=1):
         # There are various ways to create live points in nessai, such as
         # from dictionaries and numpy arrays. See nessai.livepoint for options
         d = {
-            "x": np.random.uniform(
-                self.bounds["x"][0], self.bounds["x"][1], N
-            ),
-            "y": norm(scale=5).rvs(size=N),
+            "x": rng.uniform(self.bounds["x"][0], self.bounds["x"][1], N),
+            "y": norm(scale=5).rvs(size=N, random_state=rng),
         }
         return dict_to_live_points(d)