From 7331befcc3052b4f54c8792096c4e0cf887e1cae Mon Sep 17 00:00:00 2001 From: "eitan_perlin@brown.edu" Date: Sun, 30 Aug 2026 01:34:49 -0400 Subject: [PATCH] Pin the trial-to-trial variability draws to random_state - set_seed seeds C's srand and a module-level NumPy Generator, but the *_dist simulator param mappings (e.g. ddm_st's t_dist = partial(sps.uniform.rvs, ...)) are scipy calls with no explicit random_state, so scipy drew them from NumPy's legacy global RNG, which nothing seeded. Two simulator() calls with the same random_state returned different data for every model carrying sv, sz or st, while plain ddm - which has no *_dist - stayed reproducible, hiding the gap. Any "same data across arms" comparison relying on random_state alone was not actually holding the data fixed. - Fix: bind those partials to a Generator derived from random_state, after adapt_parameters, so the draws are pinned without touching process-global state. Verified that np.random's stream position and its next draws are unchanged across seeded simulator() calls, and that reproducibility holds while the global stream is being consumed between calls. - Only integer seeds participate; any other random_state object passes through and fails where it already failed (the Cython layer requires an integer). The modulo keeps the validated negative-seed range inside default_rng's domain without changing that pre-existing failure point. - Regression test folded into the existing random_state block in tests/test_simulator.py, over ddm, ddm_st, full_ddm_rv and ddm_sdv. It fails for the three variability models without the fix and passes with it, and it also asserts that a call leaves NumPy's global RNG state byte-identical, which a global reseed or a global draw would not - so it pins where the draws come from, not only that they repeat. full_ddm is deliberately not used: it has no *_dist mappings and draws sz/sv/st internally from the generator set_seed already seeds, so it cannot exercise this path. - Alternatives considered: reseeding np.random globally inside simulator() (rejected - mutates process-global state for callers who did not ask for it), and adding the legacy seed inside set_seed itself (one line, and broader since sequential_models.pyx has the same bare np.random.* pattern, but it needs a Cython rebuild and widens the contract for every direct cssm caller). The sequential_models.pyx gap is left as a separate issue. Co-Authored-By: Claude Fable 5 - Binding is restricted to callables that accept random_state: a model config may map a parameter to any callable, and one that does not take the keyword would otherwise raise on the simulator call. scipy's rvs qualifies through **kwds. - The regression test snapshots and restores NumPy's global RNG state, so its deliberate consumption of that stream cannot make later tests order-dependent. - A positional-only random_state does not count as accepting the keyword: the value would be diverted into **kwargs while the parameter kept its default, or raise. A test pins both the exclusion and that divert-and-keep-default behaviour. - A distribution that already binds its own random_state is left alone rather than rebound, so a caller's deliberate per-distribution seeding survives; the docstring says so, since such draws follow that generator, not the seed. --- ssms/basic_simulators/simulator.py | 51 ++++++++++++++++++-- tests/test_simulator.py | 75 +++++++++++++++++++++++++++++- 2 files changed, 122 insertions(+), 4 deletions(-) diff --git a/ssms/basic_simulators/simulator.py b/ssms/basic_simulators/simulator.py index 5f59089f..738f6a22 100755 --- a/ssms/basic_simulators/simulator.py +++ b/ssms/basic_simulators/simulator.py @@ -6,6 +6,8 @@ """ from copy import deepcopy +import functools +import inspect import numbers from threading import Lock @@ -60,8 +62,8 @@ def _validate_random_state_for_c_rng(random_state: object) -> None: The Cython layer casts seeds to ``long``. On Windows that is 32-bit signed; larger values raise ``OverflowError`` inside Cython. We validate here with a clear - Python error. Non-integer seeds (e.g. ``numpy.random.Generator``) are skipped; - those paths do not use this cast in the same way. + Python error. Non-integer seeds (e.g. ``numpy.random.Generator``) skip this + range check and are rejected by the C-level seeding itself. """ if random_state is None: return @@ -664,7 +666,12 @@ def simulator( Integer passed to the C-level RNG seeding. Must lie in ``[-2**31, 2**31 - 1]`` so it fits in a 32-bit signed C ``long`` (required on Windows). ``None`` draws a seed from that range automatically. - Non-integer RNG objects may be supported on specific code paths. + Non-integer RNG objects (e.g. ``numpy.random.Generator``) are rejected + by the C-level seeding. An integer seed also pins the trial-to-trial + variability draws (``sv``/``sz``/``st``), which are taken from a + generator derived from it rather than from NumPy's global RNG. A + variability distribution that already binds its own ``random_state`` + keeps that generator, so its draws follow it rather than the seed. return_option: str Determines what the function returns. Can be either 'full' or 'minimal'. If 'full' the function returns @@ -774,6 +781,25 @@ def simulator( theta, model_config_local, n_trials ) + # Bind the trial-to-trial variability distributions to an explicit generator: + # theta's ``*_dist`` entries are scipy ``rvs`` partials, whose default draw + # source is NumPy's global RNG. Binding order is irrelevant - the simulator + # fixes the call order. The modulo is load-bearing: the validated seed range + # admits negative integers while ``default_rng`` takes only ``[0, 2**32)``, + # and the map is injective over that range. + if isinstance(random_state, numbers.Integral): + dist_rng = default_rng(int(random_state) % (2**32)) + for key in set(model_config_local.get("simulator_param_mappings", {})) | set( + model_config_local.get("simulator_fixed_params", {}) + ): + entry = theta.get(key) + if ( + callable(entry) + and _accepts_random_state(entry) + and "random_state" not in getattr(entry, "keywords", {}) + ): + theta[key] = functools.partial(entry, random_state=dist_rng) + # Make boundary dictionary boundary_dict = make_boundary_dict(model_config_local, theta) # Make drift dictionary @@ -862,3 +888,22 @@ def simulator( bin_simulator_output(x, nbins=256, max_t=-1, freq_cnt=True), axis=0 ) return x + + +def _accepts_random_state(func) -> bool: + """Report whether ``func`` takes a ``random_state`` keyword. + + A model configuration may map a parameter to any callable, and only those + that accept ``random_state`` can be bound to an explicit generator. scipy's + ``rvs`` takes it through ``**kwds``, so a variadic keyword also qualifies. + """ + try: + params = inspect.signature(func).parameters + except (TypeError, ValueError): + return False + declared = params.get("random_state") + if declared is not None: + # A positional-only parameter cannot be filled by keyword: the value + # would land in **kwargs while the parameter kept its default, or raise. + return declared.kind is not inspect.Parameter.POSITIONAL_ONLY + return any(p.kind is inspect.Parameter.VAR_KEYWORD for p in params.values()) diff --git a/tests/test_simulator.py b/tests/test_simulator.py index 25533cc2..5e276af1 100755 --- a/tests/test_simulator.py +++ b/tests/test_simulator.py @@ -2,11 +2,14 @@ import logging from unittest.mock import patch +import functools + import numpy as np +import scipy.stats as sps import pandas as pd import pytest -from ssms.basic_simulators.simulator import simulator +from ssms.basic_simulators.simulator import _accepts_random_state, simulator from ssms.config import model_config logger = logging.getLogger(__name__) @@ -330,3 +333,73 @@ def test_random_state_boundary_max_ok(): random_state=2**31 - 1, ) assert "rts" in out + + +@pytest.mark.rng_validation +@pytest.mark.parametrize("model", ["ddm", "ddm_st", "full_ddm_rv", "ddm_sdv"]) +def test_random_state_pins_variability_draws(sim_input_data, model): + """A repeated integer random_state reproduces the trial-to-trial variability + draws (sv/sz/st), not just the diffusion path, and is unaffected by unrelated + consumption of NumPy's global RNG between calls.""" + theta = dict(sim_input_data[model]["theta_dict_all_scalars"]) + # The registry defaults put sv/sz/st at 1e-3, below what a float32 RT + # resolves; raise them so a differing draw shows up in the output. + theta.update( + {k: v for k, v in {"sv": 0.5, "sz": 0.1, "st": 0.13}.items() if k in theta} + ) + + state = np.random.get_state() + try: + first = simulator(model=model, theta=theta, n_samples=500, random_state=7) + # The variability draws come from a generator derived from the seed, so + # the call leaves the process-global RNG exactly where it found it. A + # global reseed or a global draw would move it. + after_first = np.random.get_state() + np.testing.assert_array_equal(after_first[1], state[1]) + assert after_first[2] == state[2] + # unrelated consumption of the global RNG must not affect the draws + np.random.uniform(size=1234) + second = simulator(model=model, theta=theta, n_samples=500, random_state=7) + finally: + np.random.set_state(state) + + np.testing.assert_array_equal(first["rts"], second["rts"]) + np.testing.assert_array_equal(first["choices"], second["choices"]) + + +@pytest.mark.rng_validation +def test_accepts_random_state_rejects_positional_only(): + """A positional-only ``random_state`` cannot be filled by keyword. + + Binding one would either divert the value into ``**kwargs`` while the + parameter kept its default, or raise, so such callables are left alone. + """ + + def positional_only(size=1, random_state=None, /, **kwargs): + """Take random_state positionally only; keyword use lands in kwargs.""" + return random_state + + def keyword_ok(size=1, random_state=None): + """Accept random_state by keyword.""" + return random_state + + def variadic_only(**kwargs): + """Accept random_state only through **kwargs, as scipy rvs does.""" + return kwargs.get("random_state") + + assert not _accepts_random_state(positional_only) + assert _accepts_random_state(keyword_ok) + assert _accepts_random_state(variadic_only) + # the failure mode the exclusion prevents: value diverted, default kept + assert functools.partial(positional_only, random_state="RNG")() is None + + +@pytest.mark.rng_validation +def test_bound_random_state_is_not_overridden(): + """A distribution that already carries its own random_state keeps it.""" + own = np.random.default_rng(123) + dist = functools.partial(sps.uniform.rvs, loc=0.0, scale=1.0, random_state=own) + # accepted by the predicate, but already bound, so the simulator leaves it + assert _accepts_random_state(dist) + assert "random_state" in dist.keywords + assert dist.keywords["random_state"] is own