From 67c7247cdf059b09a1fd5c17330331c1ae1c6e26 Mon Sep 17 00:00:00 2001 From: Seth Gossage Date: Mon, 16 Mar 2026 11:11:19 -0500 Subject: [PATCH 1/8] add rng arguments to rejection_sampler, inverse_sampler, and histogram_sampler plus pass RNG instances to those samplers. This allows e.g., RNG instances with fixed seed to be used by these samplers. --- posydon/popsyn/independent_sample.py | 3 +- posydon/popsyn/star_formation_history.py | 4 +-- posydon/utils/common_functions.py | 44 +++++++++++++++--------- 3 files changed, 32 insertions(+), 19 deletions(-) diff --git a/posydon/popsyn/independent_sample.py b/posydon/popsyn/independent_sample.py index 13e99063c1..7bb0ee0f11 100644 --- a/posydon/popsyn/independent_sample.py +++ b/posydon/popsyn/independent_sample.py @@ -164,7 +164,8 @@ def pdf(logp): orbital_periods_M_gt_15 = 10**(rejection_sampler( size=number_of_binaries, x_lim=[np.log10(orbital_period_min), np.log10(orbital_period_max)], - pdf=pdf)) + pdf=pdf, + rng=RNG)) orbital_periods = np.where(primary_masses <= 15.0, orbital_periods_M_lt_15, diff --git a/posydon/popsyn/star_formation_history.py b/posydon/popsyn/star_formation_history.py index 29ba21be9a..b962ff61df 100644 --- a/posydon/popsyn/star_formation_history.py +++ b/posydon/popsyn/star_formation_history.py @@ -1102,7 +1102,7 @@ def get_formation_times(N_binaries, star_formation="constant", **kwargs): # prag if scenario in ["custom_linear", "custom_log10"]: custom_ages_file = kwargs.get("custom_ages_file") x, y = np.loadtxt(custom_ages_file, unpack=True) - current_binary_ages = rejection_sampler(x, y, N_binaries) + current_binary_ages = rejection_sampler(x, y, N_binaries, rng=RNG) if "log10" in scenario: current_binary_ages = 10.0**current_binary_ages return max_time - current_binary_ages @@ -1110,7 +1110,7 @@ def get_formation_times(N_binaries, star_formation="constant", **kwargs): # prag if scenario in ["custom_linear_histogram", "custom_log10_histogram"]: custom_ages_file = kwargs.get("custom_ages_file") x, y = read_histogram_from_file(custom_ages_file) - current_binary_ages = histogram_sampler(x, y) + current_binary_ages = histogram_sampler(x, y, rng=RNG) if "log10" in scenario: current_binary_ages = 10.0**current_binary_ages return max_time - current_binary_ages diff --git a/posydon/utils/common_functions.py b/posydon/utils/common_functions.py index 9716a20051..97c7e0cfc1 100644 --- a/posydon/utils/common_functions.py +++ b/posydon/utils/common_functions.py @@ -669,7 +669,7 @@ def bondi_hoyle(binary, accretor, donor, idx=-1, wind_disk_criteria=True, return np.squeeze(mdot_acc) -def rejection_sampler(x=None, y=None, size=1, x_lim=None, pdf=None): +def rejection_sampler(x=None, y=None, size=1, x_lim=None, pdf=None, rng=None): """Generate a sample from a 1d PDF using the acceptance-rejection method. Parameters @@ -684,6 +684,8 @@ def rejection_sampler(x=None, y=None, size=1, x_lim=None, pdf=None): The boundary where the pdf function is defined if passed as pdf. pdf : func The pdf function + rng : numpy.random.Generator, optional + Random number generator. If None, uses np.random.default_rng(). Returns ------- @@ -691,6 +693,8 @@ def rejection_sampler(x=None, y=None, size=1, x_lim=None, pdf=None): An array with `size` random numbers generated from the PDF. """ + if rng is None: + rng = np.random.default_rng() if pdf is None: if ((x is None) or (y is None)): raise ValueError("x and y PDF values must be specified if no PDF" @@ -702,32 +706,32 @@ def rejection_sampler(x=None, y=None, size=1, x_lim=None, pdf=None): idxs = np.argsort(x) pdf = interp1d(x.take(idxs), y.take(idxs)) - x_rand = np.random.uniform(x.min(), x.max(), size) - y_rand = np.random.uniform(0, y.max(), size) + x_rand = rng.uniform(x.min(), x.max(), size) + y_rand = rng.uniform(0, y.max(), size) values = x_rand[y_rand <= pdf(x_rand)] while values.shape[0] < size: n = size - values.shape[0] - x_rand = np.random.uniform(x.min(), x.max(), n) - y_rand = np.random.uniform(0, y.max(), n) + x_rand = rng.uniform(x.min(), x.max(), n) + y_rand = rng.uniform(0, y.max(), n) values = np.hstack([values, x_rand[y_rand <= pdf(x_rand)]]) else: if x_lim is None: raise ValueError("x_lim must be specified for passed PDF function" " in rejection sampling") - x_rand = np.random.uniform(x_lim[0], x_lim[1], size) - pdf_max = max(pdf(np.random.uniform(x_lim[0], x_lim[1], 50000))) - y_rand = np.random.uniform(0, pdf_max, size) + x_rand = rng.uniform(x_lim[0], x_lim[1], size) + pdf_max = max(pdf(rng.uniform(x_lim[0], x_lim[1], 50000))) + y_rand = rng.uniform(0, pdf_max, size) values = x_rand[y_rand <= pdf(x_rand)] while values.shape[0] < size: n = size - values.shape[0] - x_rand = np.random.uniform(x_lim[0], x_lim[1], n) - y_rand = np.random.uniform(0, pdf_max, n) + x_rand = rng.uniform(x_lim[0], x_lim[1], n) + y_rand = rng.uniform(0, pdf_max, n) values = np.hstack([values, x_rand[y_rand <= pdf(x_rand)]]) return values -def inverse_sampler(x, y, size=1): +def inverse_sampler(x, y, size=1, rng=None): """Sample from a PDF using the inverse-sampling method. Parameters @@ -738,6 +742,8 @@ def inverse_sampler(x, y, size=1): The probablity density at `x` (or a scaled version of it). size : int Number of samples to generate. + rng : numpy.random.Generator, optional + Random number generator. If None, uses np.random.default_rng(). Returns ------- @@ -745,6 +751,8 @@ def inverse_sampler(x, y, size=1): The sample drawn from the PDF. """ + if rng is None: + rng = np.random.default_rng() # x should be sorted assert np.all(np.diff(x) > 0) # y should be above 0 @@ -757,7 +765,7 @@ def inverse_sampler(x, y, size=1): total_area = cumsum_areas[-1] # start the inverse sampling - u = np.random.uniform(size=size) + u = rng.uniform(size=size) # index of "bin" where each sampled value corresponds too u_indices = np.searchsorted(cumsum_areas, u * total_area) # the area that should be covered from the end of the previous bin @@ -779,14 +787,14 @@ def inverse_sampler(x, y, size=1): if n_where_nan: assert np.all(dy_bins[where_nan] == 0) sample[where_nan] = x[u_indices][where_nan] + ( - dx_bins[where_nan] * np.random.uniform(size=n_where_nan)) + dx_bins[where_nan] * rng.uniform(size=n_where_nan)) # make sure that everything worked as expected assert np.all(np.isfinite(sample)) return sample -def histogram_sampler(x_edges, y, size=1): +def histogram_sampler(x_edges, y, size=1, rng=None): """Sample from an empirical PDF represented by a histogram. Parameters @@ -797,6 +805,8 @@ def histogram_sampler(x_edges, y, size=1): The counts (or a scaled version of them) of the histogram. size : int Number of random values to produce. + rng : numpy.random.Generator, optional + Random number generator. If None, uses np.random.default_rng(). Returns ------- @@ -804,13 +814,15 @@ def histogram_sampler(x_edges, y, size=1): The random sample. """ + if rng is None: + rng = np.random.default_rng() assert np.all(y >= 0.0) # make sure that the lengths of the input arrays are correct n_bins = len(y) assert n_bins > 0 and len(x_edges) == n_bins + 1 # first decide which will be the bin of each element in the sample - bin_sample = np.random.choice(n_bins, replace=True, p=y/sum(y), size=size) + bin_sample = rng.choice(n_bins, replace=True, p=y/sum(y), size=size) sample = np.ones(size) * np.nan @@ -818,7 +830,7 @@ def histogram_sampler(x_edges, y, size=1): bins_found = set(bin_sample) for bin_index in bins_found: in_this_bin = np.argwhere(bin_sample == bin_index)[:, 0] - sample[in_this_bin] = np.random.uniform( + sample[in_this_bin] = rng.uniform( x_edges[bin_index], x_edges[bin_index+1], size=len(in_this_bin)) assert(np.all(np.isfinite(sample))) From 8b00955a0d6811d3a850b12fd171d744198e38f6 Mon Sep 17 00:00:00 2001 From: Seth Gossage Date: Mon, 16 Mar 2026 11:30:59 -0500 Subject: [PATCH 2/8] fix unit tests --- .../unit_tests/utils/test_common_functions.py | 24 ++++++++++--------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/posydon/unit_tests/utils/test_common_functions.py b/posydon/unit_tests/utils/test_common_functions.py index ef405107b6..08239683ce 100644 --- a/posydon/unit_tests/utils/test_common_functions.py +++ b/posydon/unit_tests/utils/test_common_functions.py @@ -754,8 +754,10 @@ def mock_rand2(shape): approx(5.13970075150e-8, abs=6e-20) def test_rejection_sampler(self, monkeypatch): - def mock_uniform(low=0.0, high=1.0, size=1): - return np.linspace(low, high, num=size) + class MockRNG: + def uniform(self, low=0.0, high=1.0, size=1): + return np.linspace(low, high, num=size) + mock_rng = MockRNG() def mock_interp1d(x, y): if x[0] > x[-1]: raise ValueError @@ -786,7 +788,6 @@ def mock_pdf(x): +" function in rejection sampling"): totest.rejection_sampler(pdf=mock_pdf) # examples: - monkeypatch.setattr(np.random, "uniform", mock_uniform) monkeypatch.setattr(totest, "interp1d", mock_interp1d) tests = [(np.array([0.0, 1.0]), np.array([0.4, 0.6]), 5,\ np.array([0.0, 0.25, 0.5, 0.75, 1.0])),\ @@ -795,7 +796,7 @@ def mock_pdf(x): (np.array([1.0, 0.0]), np.array([0.2, 0.8]), 6,\ np.array([0.0, 0.2, 0.4, 0.0, 0.5, 0.0]))] for (x, y, s, r) in tests: - assert np.allclose(totest.rejection_sampler(x=x, y=y, size=s),\ + assert np.allclose(totest.rejection_sampler(x=x, y=y, size=s, rng=mock_rng),\ r) assert np.allclose(totest.rejection_sampler(x_lim=np.array([0.0,\ 1.0]),\ @@ -805,12 +806,14 @@ def mock_pdf(x): y=np.array([0.2, 0.8]),\ x_lim=np.array([0.0,\ 1.0]),\ - pdf=mock_pdf, size=5),\ + pdf=mock_pdf, size=5, rng=mock_rng),\ np.array([0.0, 0.25, 0.0, 0.0, 0.0])) def test_inverse_sampler(self, monkeypatch): - def mock_uniform(low=0.0, high=1.0, size=1): - return np.linspace(low, high, num=size) + class MockRNG: + def uniform(self, low=0.0, high=1.0, size=1): + return np.linspace(low, high, num=size) + mock_rng = MockRNG() # missing argument with raises(TypeError, match="missing 2 required positional "\ +"arguments: 'x' and 'y'"): @@ -826,21 +829,20 @@ def mock_uniform(low=0.0, high=1.0, size=1): totest.inverse_sampler(x=np.array([0.0, 1.0]),\ y=np.array([-0.4, 0.6])) # examples: - monkeypatch.setattr(np.random, "uniform", mock_uniform) assert np.allclose(totest.inverse_sampler(x=np.array([0.0, 1.0]),\ y=np.array([0.4, 0.6]),\ - size=5),\ + size=5, rng=mock_rng),\ np.array([0.0, 0.29128785, 0.54950976, 0.78388218,\ 1.0])) assert np.allclose(totest.inverse_sampler(x=np.array([0.0, 1.0]),\ y=np.array([0.6, 0.4]),\ - size=4),\ + size=4, rng=mock_rng),\ np.array([0.0, 0.2919872, 0.61952386, 1.0])) with warns(RuntimeWarning, match="invalid value encountered in "\ +"divide"): assert np.allclose(totest.inverse_sampler(x=np.array([0.0, 1.0]),\ y=np.array([0.5, 0.5]),\ - size=5),\ + size=5, rng=mock_rng),\ np.array([0.0, 0.25, 0.5, 0.75, 1.0])) def test_histogram_sampler(self, monkeypatch): From fbef3bbd9a32cc08dd3cb491b374d4a676ab2708 Mon Sep 17 00:00:00 2001 From: Seth Gossage Date: Mon, 16 Mar 2026 11:32:15 -0500 Subject: [PATCH 3/8] fix histogram_sampler unit test --- posydon/unit_tests/utils/test_common_functions.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/posydon/unit_tests/utils/test_common_functions.py b/posydon/unit_tests/utils/test_common_functions.py index 08239683ce..6c51ba812b 100644 --- a/posydon/unit_tests/utils/test_common_functions.py +++ b/posydon/unit_tests/utils/test_common_functions.py @@ -846,8 +846,10 @@ def uniform(self, low=0.0, high=1.0, size=1): np.array([0.0, 0.25, 0.5, 0.75, 1.0])) def test_histogram_sampler(self, monkeypatch): - def mock_uniform(low=0.0, high=1.0, size=1): - return np.linspace(low, high, num=size) + class MockRNG: + def uniform(self, low=0.0, high=1.0, size=1): + return np.linspace(low, high, num=size) + mock_rng = MockRNG() def mock_choice(a, size=None, replace=True, p=None): if isinstance(a, int): a=np.arange(a) @@ -872,20 +874,19 @@ def mock_choice(a, size=None, replace=True, p=None): totest.histogram_sampler(x_edges=np.array([0.0, 1.0]),\ y=np.array([0.4, 0.6])) # examples: - monkeypatch.setattr(np.random, "uniform", mock_uniform) monkeypatch.setattr(np.random, "choice", mock_choice) assert np.allclose(totest.histogram_sampler(x_edges=np.array([0.0,\ 0.5,\ 1.0]),\ y=np.array([0.2, 0.8]),\ - size=5),\ + size=5, rng=mock_rng),\ np.array([0.0, 0.5, 0.66666667, 0.83333333, 1.0])) assert np.allclose(totest.histogram_sampler(x_edges=np.array([0.0,\ 0.5,\ 1.0]\ ), y=np.array([0.2,\ 0.8]),\ - size=4),\ + size=4, rng=mock_rng),\ np.array([0.0, 0.5, 0.75, 1.0])) def test_read_histogram_from_file(self, csv_path_failing_3_data_lines,\ From db3ef654eb3d935cc95bb61fc053141fc4a5df03 Mon Sep 17 00:00:00 2001 From: Seth Gossage Date: Mon, 16 Mar 2026 11:35:47 -0500 Subject: [PATCH 4/8] adding choice method to MockRNG for test_histogram_sampler --- posydon/unit_tests/utils/test_common_functions.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/posydon/unit_tests/utils/test_common_functions.py b/posydon/unit_tests/utils/test_common_functions.py index 6c51ba812b..ce045199f4 100644 --- a/posydon/unit_tests/utils/test_common_functions.py +++ b/posydon/unit_tests/utils/test_common_functions.py @@ -849,6 +849,16 @@ def test_histogram_sampler(self, monkeypatch): class MockRNG: def uniform(self, low=0.0, high=1.0, size=1): return np.linspace(low, high, num=size) + def choice(self, a, size=None, replace=True, p=None): + if isinstance(a, int): + a = np.arange(a) + sample = [] + for (v, q) in zip(a, p): + sample += round(size*q) * [v] + if len(sample) < size: + sample += (size-len(sample)) * [a[-1]] + return np.array(sample[:size]) + mock_rng = MockRNG() def mock_choice(a, size=None, replace=True, p=None): if isinstance(a, int): From b3a35daf74e904b3278e01e3a5b6905fbc9a2900 Mon Sep 17 00:00:00 2001 From: Seth Gossage Date: Mon, 16 Mar 2026 11:40:47 -0500 Subject: [PATCH 5/8] fix missing mock_rng arg --- posydon/unit_tests/utils/test_common_functions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/posydon/unit_tests/utils/test_common_functions.py b/posydon/unit_tests/utils/test_common_functions.py index ce045199f4..8438c7d858 100644 --- a/posydon/unit_tests/utils/test_common_functions.py +++ b/posydon/unit_tests/utils/test_common_functions.py @@ -800,7 +800,7 @@ def mock_pdf(x): r) assert np.allclose(totest.rejection_sampler(x_lim=np.array([0.0,\ 1.0]),\ - pdf=mock_pdf, size=5),\ + pdf=mock_pdf, size=5, rng=mock_rng),\ np.array([0.0, 0.25, 0.0, 0.0, 0.0])) assert np.allclose(totest.rejection_sampler(x=np.array([1.0, 0.0]),\ y=np.array([0.2, 0.8]),\ From 72665436f52f4b2ac9ce7dd94b0bce2d923da317 Mon Sep 17 00:00:00 2001 From: Seth Gossage Date: Mon, 16 Mar 2026 11:45:10 -0500 Subject: [PATCH 6/8] bump version to 2.2.7 --- conda/meta.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/conda/meta.yaml b/conda/meta.yaml index eed55444e2..5a1290666a 100644 --- a/conda/meta.yaml +++ b/conda/meta.yaml @@ -1,5 +1,5 @@ {% set name = "posydon" %} -{% set version = "2.2.6" %} +{% set version = "2.2.7" %} package: name: "{{ name|lower }}" From 5f89b6a4f1b2aac54f7a48d8b54b2d2ac79660f1 Mon Sep 17 00:00:00 2001 From: Seth Gossage Date: Mon, 16 Mar 2026 13:56:38 -0500 Subject: [PATCH 7/8] update Moe+DiStefano initial sampler to use RNG if it exists --- posydon/popsyn/Moes_distributions.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/posydon/popsyn/Moes_distributions.py b/posydon/popsyn/Moes_distributions.py index b55135931f..43be3a9dd7 100644 --- a/posydon/popsyn/Moes_distributions.py +++ b/posydon/popsyn/Moes_distributions.py @@ -114,6 +114,8 @@ def __init__(self, n_M1=101, n_logP=158, n_q=91, n_e=200, self.numq = n_q self.nume = n_e + self.RNG = kwargs.get('RNG', np.random.default_rng()) + # ranges where M+D17 has statistics corrected for selection effects: # 0.8 < M1/Msun < 40 with self.numM1 steps in log space self.M1v = 10**(np.linspace(np.log10(0.8), np.log10(40.0), self.numM1)) @@ -429,7 +431,7 @@ def __call__(self, M1, M_min=0.08, M_max=150.0, all_binaries=True): else: mybinfrac = np.max(mycumPbindist()) # Generate random number myrand between 0 and 1 - myrand = np.random.rand() + myrand = self.RNG.random() # If random number < binary star fraction, generate a binary if(myrand < mybinfrac): # Given myrand, select P and corresponding index in logPv @@ -438,7 +440,7 @@ def __call__(self, M1, M_min=0.08, M_max=150.0, all_binaries=True): == min(abs(mylogP-self.logPv))) indlogP = indlogP[0] # Given M1 & P, select e from eccentricity distribution - mye = np.interp(np.random.rand(), + mye = np.interp(self.RNG.random(), self.cumedist[:,indlogP,indM1].flatten(), self.ev) # Given M1 & P, determine mass ratio distribution. # If M1 < 0.8 Msun, truncate q distribution and consider only mass @@ -454,7 +456,7 @@ def __call__(self, M1, M_min=0.08, M_max=150.0, all_binaries=True): # Set probability = 0 where q < q_min mycumqdist[self.qv <= q_min] = 0.0 # Given M1 & P, select q from cumulative mass ratio distribution - myq = np.interp(np.random.rand(), mycumqdist, self.qv) + myq = np.interp(self.RNG.random(), mycumqdist, self.qv) else: # If instead random number > binary star fraction, generate single # star @@ -463,7 +465,7 @@ def __call__(self, M1, M_min=0.08, M_max=150.0, all_binaries=True): mye = np.nan # Get metallicity Zsun = 0.02 - logZ = -2.3 + (0.176-(-2.3)) * np.random.rand() + logZ = -2.3 + (0.176-(-2.3)) * self.RNG.random() Z = Zsun*10**logZ M2s.append(M1*myq) logPs.append(mylogP) From 76108a55e370fc33222e318454a65d2e1fe3564b Mon Sep 17 00:00:00 2001 From: Seth Gossage Date: Tue, 17 Mar 2026 10:06:18 -0500 Subject: [PATCH 8/8] Update posydon/unit_tests/utils/test_common_functions.py Co-authored-by: Max <14039563+maxbriel@users.noreply.github.com> --- posydon/unit_tests/utils/test_common_functions.py | 1 - 1 file changed, 1 deletion(-) diff --git a/posydon/unit_tests/utils/test_common_functions.py b/posydon/unit_tests/utils/test_common_functions.py index 8438c7d858..c66229fdd6 100644 --- a/posydon/unit_tests/utils/test_common_functions.py +++ b/posydon/unit_tests/utils/test_common_functions.py @@ -884,7 +884,6 @@ def mock_choice(a, size=None, replace=True, p=None): totest.histogram_sampler(x_edges=np.array([0.0, 1.0]),\ y=np.array([0.4, 0.6])) # examples: - monkeypatch.setattr(np.random, "choice", mock_choice) assert np.allclose(totest.histogram_sampler(x_edges=np.array([0.0,\ 0.5,\ 1.0]),\