diff --git a/brian2/codegen/generators/cython_generator.py b/brian2/codegen/generators/cython_generator.py index baf6000b3..6fd6ce1a0 100644 --- a/brian2/codegen/generators/cython_generator.py +++ b/brian2/codegen/generators/cython_generator.py @@ -555,30 +555,23 @@ def determine_keywords(self): name="_exprel", availability_check=C99Check("exprel"), ) -_BUFFER_SIZE = 20000 + +# ============================================================================== +# Random Number Generation +# ============================================================================== +# We use a pre-compiled global RandomGenerator via cythonrng module. +# This ensures: +# 1. All code objects share the same RNG (like C++ standalone) +# 2. Proper state save/restore is possible +# 3. Identical sequences between Cython and C++ standalone backends rand_code = """ -cdef double _rand(int _idx): - cdef double **buffer_pointer = _namespace_rand_buffer - cdef double *buffer = buffer_pointer[0] - cdef _numpy.ndarray _new_rand - - if(_namespace_rand_buffer_index[0] == 0): - if buffer != NULL: - free(buffer) - _new_rand = _numpy.random.rand(_BUFFER_SIZE) - buffer = _numpy.PyArray_DATA(_new_rand) - PyArray_CLEARFLAGS(<_numpy.PyArrayObject*>_new_rand, _numpy.NPY_ARRAY_OWNDATA) - buffer_pointer[0] = buffer - - cdef double val = buffer[_namespace_rand_buffer_index[0]] - _namespace_rand_buffer_index[0] += 1 - if _namespace_rand_buffer_index[0] == _BUFFER_SIZE: - _namespace_rand_buffer_index[0] = 0 - return val -""".replace("_BUFFER_SIZE", str(_BUFFER_SIZE)) - -randn_code = rand_code.replace("rand", "randn").replace("randnom", "random") +# Import the global random number generator functions +from brian2.random.cythonrng cimport _rand, _randn +""" + +# randn is included in the same import +randn_code = rand_code poisson_code = """ cdef double _loggam(double x): @@ -665,20 +658,14 @@ def determine_keywords(self): CythonCodeGenerator, code=rand_code, name="_rand", - namespace={ - "_rand_buffer": device.rand_buffer, - "_rand_buffer_index": device.rand_buffer_index, - }, + namespace={}, ) DEFAULT_FUNCTIONS["randn"].implementations.add_implementation( CythonCodeGenerator, code=randn_code, name="_randn", - namespace={ - "_randn_buffer": device.randn_buffer, - "_randn_buffer_index": device.randn_buffer_index, - }, + namespace={}, ) DEFAULT_FUNCTIONS["poisson"].implementations.add_implementation( CythonCodeGenerator, diff --git a/brian2/codegen/runtime/cython_rt/templates/common_group.pyx b/brian2/codegen/runtime/cython_rt/templates/common_group.pyx index b9ea2572e..6db51f9c0 100644 --- a/brian2/codegen/runtime/cython_rt/templates/common_group.pyx +++ b/brian2/codegen/runtime/cython_rt/templates/common_group.pyx @@ -69,6 +69,7 @@ cdef extern from "dynamic_array.h": size_t rows() size_t cols() size_t stride() + {% endmacro %} {% macro before_run() %} diff --git a/brian2/devices/cpp_standalone/brianlib/randomgenerator.h b/brian2/devices/cpp_standalone/brianlib/randomgenerator.h new file mode 100644 index 000000000..54d928492 --- /dev/null +++ b/brian2/devices/cpp_standalone/brianlib/randomgenerator.h @@ -0,0 +1,157 @@ +#ifndef _BRIAN_RANDOMGENERATOR_H +#define _BRIAN_RANDOMGENERATOR_H + +#include +#include +#include +#include +#include +#include + +/** + * @brief Random number generator class that provides reproducible + * random sequences across different Brian2 backends. + * + * Uses std::mt19937 (Mersenne Twister) with the Jean-Sebastien Roy + * algorithm for converting to uniform doubles, ensuring identical + * sequences when seeded with the same value. + * + * This class is used by both C++ standalone mode and Cython runtime mode + * to ensure cross-backend reproducibility. + */ +class RandomGenerator +{ +private: + std::mt19937 gen; + double stored_gauss; + bool has_stored_gauss; + +public: + RandomGenerator() : has_stored_gauss(false) + { + seed(); + } + + /** + * @brief Seed with a random value from the system. + */ + void seed() + { + std::random_device rd; + gen.seed(rd()); + has_stored_gauss = false; + } + + /** + * @brief Seed with a specific value for reproducibility. + * @param seed The seed value. + */ + void seed(unsigned long seed_value) + { + gen.seed(seed_value); + has_stored_gauss = false; + } + + /** + * @brief Generate a uniform random double in [0, 1) + * + * Uses the Jean-Sebastien Roy algorithm to extract 53 bits + * of randomness from two consecutive MT19937 outputs. + * This ensures reproducibility with older Brian2 versions + * and across different backends. + * + * The algorithm: + * - Takes two 32-bit outputs from MT19937 + * - Extracts 27 bits from the first (shift right 5) + * - Extracts 26 bits from the second (shift right 6) + * - Combines them into a 53-bit integer (full double precision mantissa) + * - Divides by 2^53 to get a value in [0, 1) + */ + double rand() + { + /* shifts : 67108864 = 0x4000000 = 2^26 + * 9007199254740992 = 0x20000000000000 = 2^53 */ + const long a = gen() >> 5; // Upper 27 bits + const long b = gen() >> 6; // Upper 26 bits + return (a * 67108864.0 + b) / 9007199254740992.0; + } + + /** + * @brief Generate a standard normal random double (mean=0, std=1) + * + * Uses the Box-Muller transform with rejection sampling. + * Generates two values at once and caches one for the next call, + * making every other call essentially free. + */ + double randn() + { + if (has_stored_gauss) + { + const double tmp = stored_gauss; + has_stored_gauss = false; + return tmp; + } + else + { + double f, x1, x2, r2; + + do + { + x1 = 2.0 * rand() - 1.0; + x2 = 2.0 * rand() - 1.0; + r2 = x1 * x1 + x2 * x2; + } while (r2 >= 1.0 || r2 == 0.0); + + /* Box-Muller transform */ + f = sqrt(-2.0 * log(r2) / r2); + /* Keep for next call */ + stored_gauss = f * x1; + has_stored_gauss = true; + return f * x2; + } + } + + /** + * @brief Get the complete internal state as a string + * + * This captures the MT19937 state plus the Box-Muller cache, + * allowing exact restoration of the generator state. + */ + std::string get_state() const { + std::ostringstream oss; + oss << gen << " " << stored_gauss << " " << has_stored_gauss; + return oss.str(); + } + + /** + * @brief Set the complete internal state from a string. + * + * Restores the exact state previously captured by get_state(). + */ + void set_state(const std::string& state){ + std::istringstream iss(state); + iss >> gen >> stored_gauss >> has_stored_gauss; + } + + // Allow exporting/setting the internal state of the random generator + friend std::ostream &operator<<(std::ostream &out, const RandomGenerator &rng); + friend std::istream &operator>>(std::istream &in, RandomGenerator &rng); +}; + +/** + * @brief Stream output operator for serializing generator state. + */ +inline std::ostream &operator<<(std::ostream &out, const RandomGenerator &rng) +{ + return out << rng.gen << " " << rng.stored_gauss << " " << rng.has_stored_gauss; +} + +/** + * @brief Stream input operator for deserializing generator state. + */ +inline std::istream &operator>>(std::istream &in, RandomGenerator &rng) +{ + return in >> rng.gen >> rng.stored_gauss >> rng.has_stored_gauss; +} + +#endif // _BRIAN_RANDOMGENERATOR_H diff --git a/brian2/devices/cpp_standalone/templates/objects.cpp b/brian2/devices/cpp_standalone/templates/objects.cpp index 1762089e9..1a78230f9 100644 --- a/brian2/devices/cpp_standalone/templates/objects.cpp +++ b/brian2/devices/cpp_standalone/templates/objects.cpp @@ -19,6 +19,7 @@ set_variable_from_value(name, {{array_name}}, var_size, (char)atoi(s_value.c_str #include "brianlib/clocks.h" #include "brianlib/dynamic_array.h" #include "brianlib/stdint_compat.h" +#include "brianlib/randomgenerator.h" #include "network.h" #include #include @@ -37,16 +38,6 @@ std::string results_dir = "results/"; // can be overwritten by --results_dir co // For multhreading, we need one generator for each thread. std::vector< RandomGenerator > _random_generators; -std::ostream& operator<<(std::ostream& out, const RandomGenerator& rng) -{ - return out << rng.gen; -} - -std::istream& operator>>(std::istream& in, RandomGenerator& rng) -{ - return in >> rng.gen; -} - //////////////// networks ///////////////// {% for net in networks | sort(attribute='name') %} Network {{net.name}}; @@ -413,6 +404,7 @@ void _dealloc_arrays() #include "brianlib/clocks.h" #include "brianlib/dynamic_array.h" #include "brianlib/stdint_compat.h" +#include "brianlib/randomgenerator.h" #include "network.h" #include #include @@ -423,64 +415,6 @@ namespace brian { extern std::string results_dir; -class RandomGenerator { - private: - std::mt19937 gen; - double stored_gauss; - bool has_stored_gauss = false; - public: - RandomGenerator() { - seed(); - } - void seed() { - std::random_device rd; - gen.seed(rd()); - has_stored_gauss = false; - } - void seed(unsigned long seed) { - gen.seed(seed); - has_stored_gauss = false; - } - // Allow exporting/setting the internal state of the random generator - friend std::ostream& operator<<(std::ostream& out, const RandomGenerator& rng); - friend std::istream& operator>>(std::istream& in, RandomGenerator& rng); - - double rand() { - /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */ - const long a = gen() >> 5; - const long b = gen() >> 6; - return (a * 67108864.0 + b) / 9007199254740992.0; - } - - double randn() { - if (has_stored_gauss) { - const double tmp = stored_gauss; - has_stored_gauss = false; - return tmp; - } - else { - double f, x1, x2, r2; - - do { - x1 = 2.0*rand() - 1.0; - x2 = 2.0*rand() - 1.0; - r2 = x1*x1 + x2*x2; - } - while (r2 >= 1.0 || r2 == 0.0); - - /* Box-Muller transform */ - f = sqrt(-2.0*log(r2)/r2); - /* Keep for next call */ - stored_gauss = f*x1; - has_stored_gauss = true; - return f*x2; - } - } -}; - -extern std::ostream& operator<<(std::ostream& out, const RandomGenerator& rng); -extern std::istream& operator>>(std::istream& in, RandomGenerator& rng); - // In OpenMP we need one state per thread extern std::vector< RandomGenerator > _random_generators; diff --git a/brian2/devices/device.py b/brian2/devices/device.py index ea6298757..d115d9ce1 100644 --- a/brian2/devices/device.py +++ b/brian2/devices/device.py @@ -479,12 +479,6 @@ def __init__(self): #: objects). Arrays in this dictionary will disappear as soon as the #: last reference to the `Variable` object used as a key is gone self.arrays = WeakKeyDictionary() - # Note that the buffers only store a pointer to the actual random - # numbers -- the buffer will be filled in Cython code - self.randn_buffer = np.zeros(1, dtype=np.intp) - self.rand_buffer = np.zeros(1, dtype=np.intp) - self.randn_buffer_index = np.zeros(1, dtype=np.int32) - self.rand_buffer_index = np.zeros(1, dtype=np.int32) def __getstate__(self): state = dict(self.__dict__) @@ -594,25 +588,52 @@ def seed(self, seed=None): The seed value for the random number generator, or ``None`` (the default) to set a random seed. """ + # Seed the global Cython RandomGenerator + from brian2.random.cythonrng import seed as rng_seed + + rng_seed(seed) + + # Also seed numpy for any code that might use it directly np.random.seed(seed) - self.rand_buffer_index[:] = 0 - self.randn_buffer_index[:] = 0 def get_random_state(self): + """ + Return a representation of the current random number generator state. + + This captures the exact internal state of the RandomGenerator, + allowing precise restoration to this point in the random sequence. + + Returns + ------- + dict + A dictionary containing: + - 'rng_state': The internal state of the RandomGenerator + - 'numpy_state': The state of NumPy's random generator + """ + from brian2.random.cythonrng import ( + get_state as rng_get_state, + ) + return { + "rng_state": rng_get_state(), "numpy_state": np.random.get_state(), - "rand_buffer_index": np.array(self.rand_buffer_index), - "rand_buffer": np.array(self.rand_buffer), - "randn_buffer_index": np.array(self.randn_buffer_index), - "randn_buffer": np.array(self.randn_buffer), } def set_random_state(self, state): + """ + Restore the random number generator to a previously saved state. + + Parameters + ---------- + state : dict + A state dictionary previously returned by get_random_state(). + """ + from brian2.random.cythonrng import ( + set_state as rng_set_state, + ) + + rng_set_state(state["rng_state"]) np.random.set_state(state["numpy_state"]) - self.rand_buffer_index[:] = state["rand_buffer_index"] - self.rand_buffer[:] = state["rand_buffer"] - self.randn_buffer_index[:] = state["randn_buffer_index"] - self.randn_buffer[:] = state["randn_buffer"] class Dummy: diff --git a/brian2/random/__init__.py b/brian2/random/__init__.py index e69de29bb..c914ddf39 100644 --- a/brian2/random/__init__.py +++ b/brian2/random/__init__.py @@ -0,0 +1,11 @@ +# brian2/random/__init__.py +""" +Random number generation for Brian2. + +This module provides a unified random number generator that produces +identical sequences in both Cython runtime and C++ standalone modes. +""" + +from .cythonrng import seed, get_state, set_state, rand, randn + +__all__ = ["seed", "get_state", "set_state", "rand", "randn"] diff --git a/brian2/random/cythonrng.pxd b/brian2/random/cythonrng.pxd new file mode 100644 index 000000000..c6ae1bfdd --- /dev/null +++ b/brian2/random/cythonrng.pxd @@ -0,0 +1,6 @@ +# Cython declaration file for the global random number generator +# Other Cython modules can use: from brian2.random.cythonrng cimport _rand, _randn + + # Note we accept (and ignore) the _idx parameter for backwards compatibility +cdef double _rand(int _idx) noexcept nogil +cdef double _randn(int _idx) noexcept nogil diff --git a/brian2/random/cythonrng.pyx b/brian2/random/cythonrng.pyx new file mode 100644 index 000000000..4682c5c37 --- /dev/null +++ b/brian2/random/cythonrng.pyx @@ -0,0 +1,94 @@ +# cython: language_level=3 +# distutils: language = c++ +# distutils: include_dirs = brian2/devices/cpp_standalone/brianlib + +""" +Global random number generator for Brian2's Cython runtime. + +This module provides a single global RandomGenerator instance that is shared +by all generated Cython code. This ensures: +1. Consistent random number sequences across all code objects +2. Identical behavior to C++ standalone mode (which also uses a global generator) +3. Proper state save/restore functionality +""" + +from libcpp.string cimport string + +cdef extern from "randomgenerator.h": + cdef cppclass RandomGenerator: + Randomgenerator() except + + void seed() except + + void seed(unsigned long) except + + double rand() noexcept nogil + double randn() noexcept nogil + string get_state() + void set_state(string) except + + +# The ONE global random generator instance +# This is created when the module is first imported and lives for the entire process +cdef RandomGenerator _global_rng + + +##### C-level functions (for use by generated Cython code via cimport) ##### +cdef double _rand(int _idx) noexcept nogil: + return _global_rng.rand() # Note we accept (and ignore) the _idx parameter for backwards compatibility + +cdef double _randn(int _idx) noexcept nogil: + return _global_rng.randn() # Note we accept (and ignore) the _idx parameter for backwards compatibility + + + +##### Python-level functions (for use by Brian2's Python code) ##### +def seed(seed_value=None): + """ + Seed the global random number generator. + + Parameters + ---------- + seed_value : int or None + If None, seed with system entropy (non-deterministic). + If an integer, seed with that value (deterministic/reproducible). + """ + if seed_value is None: + _global_rng.seed() + else: + _global_rng.seed(seed_value) + +def get_state(): + """ + Get the complete internal state of the random generator. + + Returns a string that captures: + - The full MT19937 internal state (624 x 32-bit words) + - The Box-Muller cached value (for randn) + - Whether there's a cached value + + This allows exact restoration to this point in the sequence. + + Returns + ------- + str + The serialized state that can be passed to set_state(). + """ + cdef string state = _global_rng.get_state() + return state.decode('utf-8') + +def set_state(state_str): + """ + Restore the random generator to a previously saved state. + + Parameters + ---------- + state_str : str + A state string previously returned by get_state(). + """ + cdef bytes state_bytes = state_str.encode('utf-8') + cdef string state = state_bytes + _global_rng.set_state(state) + + +def rand(): + return _global_rng.rand() + +def randn(): + return _global_rng.randn() diff --git a/brian2/tests/test_neurongroup.py b/brian2/tests/test_neurongroup.py index d9de59192..f80ce1c9c 100644 --- a/brian2/tests/test_neurongroup.py +++ b/brian2/tests/test_neurongroup.py @@ -2017,9 +2017,7 @@ def test_random_values_fixed_seed(): ), ("RuntimeDevice", "cython", None): ( [0.1636023, 0.76229608, 0.74945305, 0.82121212, 0.82669968], - # Cython uses a buffer for the random values that it gets from numpy, the - # values for the second call are therefore different - [-0.24349748, 1.1164414, -1.97421849, 1.58092889, -0.06444478], + [-0.7758696, 0.13295831, 0.87360834, -1.21879122, 0.62980314], ), ("CPPStandaloneDevice", None, 1): ( [0.1636023, 0.76229608, 0.74945305, 0.82121212, 0.82669968], diff --git a/setup.py b/setup.py index 816e1b92e..6fb7d240c 100644 --- a/setup.py +++ b/setup.py @@ -41,6 +41,7 @@ def require_cython_extension(module_path, module_name,extra_include_dirs=None): extensions.append(spike_queue_ext) +# Dynamic Array extension dynamic_array_ext = require_cython_extension( module_path=["brian2", "memory"], module_name="cythondynamicarray", @@ -49,5 +50,21 @@ def require_cython_extension(module_path, module_name,extra_include_dirs=None): extensions.append(dynamic_array_ext) +# Random number generator extension +rng_ext = require_cython_extension( + module_path=["brian2", "random"], + module_name="cythonrng", + extra_include_dirs=["brian2/devices/cpp_standalone/brianlib"] +) +extensions.append(rng_ext) -setup(ext_modules=extensions) + +setup( + ext_modules=extensions, + # Include .pxd files so they get installed + package_data={ + 'brian2.random': ['*.pxd'], + }, + # Make sure package data is included + include_package_data=True, +)