Move Cython random number generation to C++ mechanism

[mstimberg]

In C++ standalone mode, we use the random generator mechanism from the C++ standard library, wrapped around in our own functions to get rand, randn, and poisson implementations. See #1559 for the PR that introduced this for C++ standalone mode. Note that the reason for using our own functions for rand and so on is that this way, we are generating exactly the same random numbers as before. In Cython, we are instead calling out to the numpy library (which under the hood does the same thing that we are doing in C++ standalone). It always asks for 20000 numbers at a time, to avoid going through numpy for every single number. It would be good to avoid this and use the same approach as C++ standalone mode. The C++ standalone code is distributed over a few files:

• the objects.h template (note that it's in the same file as objects.cpp contains the definition of the RandomGenerator class, which is wrapping the standard library function and adds the rand, randn, and seed functionality
• The generate_rand_code function in the brian2/devices/cpp_standalone/codeobject.py generates the actual rand/randn function definition – this is a bit more complicated than necessary for Cython, because in C++ mode we can run things in multiple threads and in that case we have more one RandomGenerator object per thread.
• The cpp_generator.py file adds an implementation for the poisson function (which calls out to rand).

In principle, we should be able to reuse all this from the Cython code. The main obstacle is that we don't have the objects.h file, so we should move the RandomGenerator definition out of it and into an independent file that we can then use both from Cython and from C++ standalone mode.

[Peace-png]

Interesting approach! I'm working on emotional learning SNNs in Brian2 and noticed the random number generation differences between modes.

A few questions about the unified approach:

• Would this maintain the same random sequences for reproducibility?
• Any performance benchmarks showing the numpy batch approach vs direct C++ calls?
• Does this affect the thread-safety model in any way?

This would definitely benefit research reproducibility across different Brian2 backends.

[mstimberg]

Hi @Peace-png ,a few quick replies:

• with the outlined approach, rerunning old simulations with the Cython backend would potentially give different random values. On the other hand, values should then be reproducible across (single-threaded) backends.
• we didn't implement this yet, so no benchmarks. But given that we used a buffer for the numpy call, I don't think it will make a big difference
• the Cython target is always single-threaded, so this shouldn't affect thread stuff in any way

[itu-itis21-nguyen21]

Hi! Is this task open to an external contributor?

[mstimberg]

Hi! Is this task open to an external contributor?

It sure is! I wouldn't say it is a very easy issue, since it touches quite a few different places in the code, but then we don't have many easy issues in general… If you want to have a go at it, please go ahead. The very first step (which could be in its own PR), would be to move the RandomGenerator class out of the objects.h template and into a separate file, without breaking the C++ standalone mode.

[itu-itis21-nguyen21]

Thanks for your help. I'll try doing that and let you know how it goes.

[Legend101Zz]

Hi @itu-itis21-nguyen21 , thanks for your interest to contribute on this! I wanted to give you a detailed summary of how the current random number generation works in both modes and why we'd like to change that. I'm going to show you the actual code, exact file locations, and specific functions so you can really understand what's happening under the hood. Hope this helps you:

The Current Situation: Two Different Worlds

Right now, Brian2 has two completely different approaches to generating random numbers, depending on which backend you're using. Let me walk you through both systems with the actual code so you can see exactly why we want to unify them.

How Runtime Mode (Cython) Works Today

When you run a Brian2 simulation in the default runtime mode, the random number generation uses what we call a "buffering" strategy. Let me show you exactly how this works by tracing through the actual code.

Step 1: Device Initialization

The foundation starts in brian2/devices/device.py, specifically in the RuntimeDevice class. When Brian2 starts up, it creates a RuntimeDevice instance that looks like this:

class RuntimeDevice(Device):
    """
    The default device used in Brian, state variables are stored as numpy
    arrays in memory.
    """

    def __init__(self):
        super().__init__()
        # Mapping from Variable objects to numpy arrays
        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)

These four arrays are crucial to understanding the system. The rand_buffer and randn_buffer are single-element arrays that store memory addresses (pointers) to where the actual random numbers will live. The rand_buffer_index and randn_buffer_index are counters that track which position in the buffer we're currently at. They start at zero, meaning the buffer is empty and needs to be filled.

When you call the seed() method on the device, it does this:

def seed(self, seed=None):
    """
    Set the seed for the random number generator.

    Parameters
    ----------
    seed : int, optional
        The seed value for the random number generator, or ``None`` (the
        default) to set a random seed.
    """
    np.random.seed(seed)
    self.rand_buffer_index[:] = 0
    self.randn_buffer_index[:] = 0

Notice that it seeds NumPy's random number generator and resets both buffer indices to zero. This forces a refill of the buffers on the next random number request. The actual buffers themselves aren't reset because they'll be overwritten anyway when they get refilled.

Step 2: Registering the Cython Implementation

Now let's look at where the actual random number generation code is defined. This happens in brian2/codegen/generators/cython_generator.py. First, there's a constant that defines the buffer size:

_BUFFER_SIZE = 20000

This means every time we refill the buffer, we're asking NumPy for twenty thousand random numbers in one shot. Now here's the actual Cython code that implements the rand() function:

rand_code = """
cdef double _rand(int _idx):
    cdef double **buffer_pointer = <double**>_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 = <double *>_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))

The _namespace_rand_buffer will be replaced at code generation time with the actual device.rand_buffer array. We're casting it to a double pointer pointer because device.rand_buffer is a NumPy array that contains one element which is itself a pointer. So we need a pointer to that pointer to be able to modify it.

Then cdef double *buffer = buffer_pointer[0] dereferences once to get the actual pointer to the buffer of random numbers. If this is the first call or the buffer has been refilled, this will point to the actual array of twenty thousand doubles.

The conditional if(_namespace_rand_buffer_index[0] == 0) checks whether we need to refill the buffer. When the index is zero, it means either this is the first call ever, or we just wrapped around after using all twenty thousand numbers from the previous batch.

Inside that conditional, if buffer != NULL: free(buffer) checks if there was a previous buffer and frees it to avoid memory leaks. Then comes the key line: _new_rand = _numpy.random.rand(_BUFFER_SIZE). This is where we cross into Python and call NumPy's random number generator to give us twenty thousand uniform random numbers. This returns a NumPy array.

The next line buffer = <double *>_numpy.PyArray_DATA(_new_rand) gets the raw C pointer to the underlying data in that NumPy array. NumPy arrays are just wrappers around C arrays, and this function gives us direct access to that C array.

The following line PyArray_CLEARFLAGS(<_numpy.PyArrayObject*>_new_rand, _numpy.NPY_ARRAY_OWNDATA) tells NumPy that it doesn't own this data anymore. This is necessary because we're keeping the pointer around after the NumPy array object goes out of scope. Without this flag change, NumPy would free the memory when _new_rand gets garbage collected, which would cause our pointer to become invalid.

Finally buffer_pointer[0] = buffer stores this new pointer back into the device's buffer holder so future calls can access it.

After the conditional, cdef double val = buffer[_namespace_rand_buffer_index[0]] retrieves the random number at the current index position. Then _namespace_rand_buffer_index[0] += 1 moves to the next position. The check if _namespace_rand_buffer_index[0] == _BUFFER_SIZE sees if we've used all twenty thousand numbers, and if so, wraps back to zero with _namespace_rand_buffer_index[0] = 0. This will trigger a refill on the next call.

After defining this code string, it gets registered with Brian2's function system :

device = all_devices["runtime"]
DEFAULT_FUNCTIONS["rand"].implementations.add_implementation(
    CythonCodeGenerator,
    code=rand_code,
    name="_rand",
    namespace={
        "_rand_buffer": device.rand_buffer,
        "_rand_buffer_index": device.rand_buffer_index,
    },
)

This registration tells Brian2 that whenever it's generating Cython code and encounters a call to the rand() function in user code, it should use this implementation. The namespace dictionary provides the actual device buffer and index arrays, which will be injected into the generated code. The exact same pattern is used for randn(), just with different buffer variables.

Step 3: Code Generation and Execution

Now let's see what happens when you actually write Brian2 code that uses random numbers. Suppose you write:

from brian2 import *

G = NeuronGroup(100, 'v : 1')
G.v = 'rand()'
run(1*ms)

During code generation, Brian2 analyzes the statement G.v = 'rand()' and determines it needs to generate code to set each neuron's voltage to a random value. It looks up the implementation for rand() in the CythonCodeGenerator and finds the code we just examined. It then generates a Cython file (typically in a temporary directory like ~/.cython/brian_extensions/) that might look something like this:

# This is a simplified version of what gets generated
# Actual generated code has more boilerplate

import numpy as _numpy
cimport numpy as _numpy
from libc.stdlib cimport free
from numpy cimport PyArray_CLEARFLAGS, PyArray_DATA, NPY_ARRAY_OWNDATA, PyArrayObject

# The namespace variables get injected here with actual memory addresses
# These are pointers to the device's buffer holders
cdef long long _namespace_rand_buffer_ptr = 140523847392848  # example address
cdef long long _namespace_rand_buffer_index_ptr = 140523847392912  # example address

# Cast them to the appropriate types
cdef double **_namespace_rand_buffer = <double**><void*>_namespace_rand_buffer_ptr
cdef int32_t *_namespace_rand_buffer_index = <int32_t*><void*>_namespace_rand_buffer_index_ptr

# The _rand function from rand_code is inserted here
cdef double _rand(int _idx):
    cdef double **buffer_pointer = <double**>_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(20000)
        buffer = <double *>_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] == 20000:
        _namespace_rand_buffer_index[0] = 0
    return val

# The main execution function
def _run_stateupdater():
    # Get the array for neuron voltages
    cdef double[:] _v = _namespace_v
    
    # Loop over all neurons
    cdef int _idx
    for _idx in range(100):
        # Call _rand and assign to v
        _v[_idx] = _rand(_idx)

This Cython code then gets compiled by the Cython compiler into C code, which then gets compiled by your C compiler into a shared library (a .so file on Linux, .pyd on Windows). When you call run(), Brian2 imports this compiled module and calls the _run_stateupdater() function.

During execution, here's the flow for the first few neurons:

First neuron (index zero): The function _rand(0) is called. The index check _namespace_rand_buffer_index[0] == 0 is true because this is the first call. It enters the refill block, calls _numpy.random.rand(20000) which generates twenty thousand random numbers using NumPy's random number generator. Let's say the first few numbers are [0.374540, 0.950714, 0.731994, ...]. It stores the pointer to this array and returns 0.374540, then increments the index to one.

Second neuron (index one): The function _rand(1) is called. Now the index is one (not zero), so it skips the refill block. It directly accesses buffer[1] which gives 0.950714, returns it, and increments the index to two.

Third neuron (index two): Similar flow, returns 0.731994 and increments to three.

This continues for all one hundred neurons, using numbers from indices zero through ninety-nine of the twenty thousand number buffer. The remaining 19,900 numbers stay in the buffer for future calls.

How C++ Standalone Mode Works currently

The C++ standalone mode takes a fundamentally different architectural approach. Instead of executing code interactively through Python, it generates complete C++ source files, compiles them into a standalone executable, and then runs that executable. Let me show you exactly how the random number generation fits into this.

Step 1: The RandomGenerator Class

The heart of the system is the RandomGenerator class, which is currently embedded in the objects.h template. You can find this in brian2/devices/cpp_standalone/templates/objects.cpp. Here's the exact code:

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;
        }
        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;
            }
        }
};

Let me explain what's happening in each part of this class, because the implementation details matter for understanding why we want to reuse this code:

The private member std::mt19937 gen is the actual random number generator from the C++ standard library. MT19937 refers to the Mersenne Twister algorithm with a period of 2^19937-1, which is an extremely well-tested and widely-used pseudorandom number generator. It's the same basic algorithm that many implementations use, though the details of how you extract and transform its output matter for reproducibility.

The stored_gauss and has_stored_gauss members are for the Box-Muller optimization in the normal distribution generator, which I'll explain when we get to the randn() method.

The default constructor just calls seed() with no arguments, which uses std::random_device to get a truly random seed from the operating system. The std::random_device typically reads from /dev/urandom on Linux or uses the Windows cryptographic API, giving you a non-deterministic starting point.

Now let's look at the rand() method in detail, because this is where the Jean-Sebastien Roy algorithm comes in:

double rand() {
    /* shifts : 67108864 = 0x4000000, 9007199254740992 = 0x20000000000000 */
    const long a = gen() >> 5;
    const long b = gen() >> 6;
    return (a * 67108864.0 + b) / 9007199254740992.0;
}

The call to gen() returns an unsigned 32-bit integer in the range from zero to 2^32-1. We make two consecutive calls to get two independent random 32-bit values. The first one gets shifted right by five bits with >> 5, which discards the five least significant bits and leaves us with the upper 27 bits. This value is stored in a. The second call gets shifted right by six bits with >> 6, keeping the upper 26 bits, stored in b.

Now we have 27 bits in a and 26 bits in b, for a total of 53 bits of randomness. This is important because a double-precision floating point number has a 53-bit mantissa, so we're using the full precision available. The formula (a * 67108864.0 + b) / 9007199254740992.0 combines these bits. Let's break down the constants:

The constant 67108864 is 2^26, which is exactly the right value to shift a up by 26 bit positions. When we multiply a by this, we're effectively placing the 27 bits of a in bit positions 26 through 52, and then adding the 26 bits of b in positions 0 through 25. This gives us a 53-bit integer.

The constant 9007199254740992 is 2^53, which is the number of possible values in our 53-bit space. Dividing by this converts our integer to a floating point number in the range from zero to one. Specifically, the range is [0, 1) because the numerator can be at most 2^53-1, giving us (2^53-1)/2^53 which is just slightly less than one.

This specific algorithm is what gives us reproducibility with older versions of Brian2 that used this implementation. If we just called the standard library's distribution functions directly, we'd get different sequences even with the same seed, because different implementations might extract bits differently or use different transformation algorithms.

Now for the randn() method, which implements the Box-Muller transform:

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;
    }
}

The Box-Muller transform is a method for converting uniform random numbers into normally distributed random numbers. The clever part is that it generates two independent normal random values at once from two uniform values. The first thing the method does is check whether we have a stored value from the previous call. If we do, we return it immediately and mark that we've used it. This makes every other call essentially free because we're just returning a pre-computed value.

If we don't have a stored value, we need to generate a new pair. We generate two uniform random numbers x1 and x2 in the range (-1, 1) by calling our rand() method and transforming the result. We compute r2 = x1*x1 + x2*x2, which is the squared distance from the origin. The do-while loop keeps generating pairs until we get one that falls inside the unit circle (r2 < 1.0) and isn't exactly at the origin (r2 ≠ 0). This rejection sampling ensures we have a uniform distribution over the disk.

Once we have a valid pair, we apply the Box-Muller formula. The factor f = sqrt(-2.0*log(r2)/r2) is derived from the transformation equations. We multiply this factor by both x1 and x2 to get two independent normally distributed values. We store one in stored_gauss for next time and return the other immediately. This means if you call randn() six times, you're actually only going through the full generation process three times, with the other three calls just returning stored values.

Step 2: Template Integration and Initialization

This RandomGenerator class sits inside a larger template structure ...

namespace brian {

extern std::string results_dir;

class RandomGenerator {
    // ... the class definition we just examined ...
};

// In OpenMP we need one state per thread
extern std::vector< RandomGenerator > _random_generators;

//////////////// clocks ///////////////////
{% for clock in clocks | sort(attribute='name') %}
extern Clock {{clock.name}};
{% endfor %}

// ... more declarations ...

} // namespace brian

The line extern std::vector< RandomGenerator > _random_generators; is declaring (but not defining) a vector that will hold one RandomGenerator for each thread. The extern keyword means the actual definition lives elsewhere, specifically in the generated objects.cpp file.

Looking at the corresponding cpp_file() ...

namespace brian {

std::string results_dir = "results/";

// For multithreading, we need one generator for each thread.
std::vector< RandomGenerator > _random_generators;

// ... more definitions ...

} // namespace brian

This creates the actual vector object. At this point it's empty; the actual RandomGenerator objects get added during initialization. That happens in the generated run.cpp file, in the _init_arrays() function. Looking at the template in brian2/devices/cpp_standalone/templates/run.cpp, around line 40 in the cpp_file() macro:

void _init_arrays()
{
    using namespace brian;

    // ... other initialization code ...

    // Random number generator states
    std::random_device rd;
    for (int i=0; i<{{openmp_pragma('get_num_threads')}}; i++)
        _random_generators.push_back(RandomGenerator());
}

The template expression {{openmp_pragma('get_num_threads')}} gets replaced at generation time with either a literal number (like 1 for single-threaded) or with a call to omp_get_max_threads() if OpenMP is enabled. This creates exactly as many generators as we have threads. Each RandomGenerator() constructor call will seed itself with std::random_device, giving each thread an independent random sequence.

Step 3: Dynamic Wrapper Generation

Now here's where it gets interesting. Brian2 doesn't know at import time whether it will be running with OpenMP or not, because that's determined by user preferences that might not be set until later. So the wrapper functions that call the RandomGenerator methods need to be generated dynamically during code generation.

This happens in brian2/devices/cpp_standalone/codeobject.py, in the generate_rand_code() function ....

def generate_rand_code(rand_func, owner):
    """
    Generate wrapper code for rand/randn that works with threading.
    This is called during code generation for each code object that uses rand/randn.
    """
    # Check the user's preference for OpenMP threads
    nb_threads = prefs.devices.cpp_standalone.openmp_threads
    
    if nb_threads == 0:  # No OpenMP - single threaded
        thread_number = "0"  # Always use generator 0
    else:  # OpenMP enabled - multiple threads
        thread_number = "omp_get_thread_num()"  # Get current thread ID
    
    # Validate the function name
    if rand_func not in ["rand", "randn"]:
        raise AssertionError(rand_func)
    
    # Generate the wrapper function code
    code = """
           inline double _%RAND_FUNC%(const int _vectorisation_idx) {
               return brian::_random_generators[%THREAD_NUMBER%].%RAND_FUNC%();
           }
           """
    
    code = replace(
        code,
        {
            "%THREAD_NUMBER%": thread_number,
            "%RAND_FUNC%": rand_func,
        },
    )
    
    return {"support_code": code}

Let's trace through this with an example. Suppose the user has set prefs.devices.cpp_standalone.openmp_threads = 0 (single-threaded mode). When Brian2 generates code for a code object that uses rand(), it calls generate_rand_code("rand", owner). The function checks nb_threads and finds it's zero, so it sets thread_number = "0". It then does the string substitution to create:

inline double _rand(const int _vectorisation_idx) {
    return brian::_random_generators[0].rand();
}

This is a simple wrapper that always calls the rand() method on the first (and only) generator in the vector. The _vectorisation_idx parameter is accepted but not used, because in single-threaded mode we don't need it.

Now suppose instead the user has enabled OpenMP with four threads. The same function would set thread_number = "omp_get_thread_num()" and generate:

inline double _rand(const int _vectorisation_idx) {
    return brian::_random_generators[omp_get_thread_num()].rand();
}

Now when this code runs in parallel, each thread calls omp_get_thread_num() which returns zero, one, two, or three depending on which thread is executing. This ensures each thread uses its own generator from the vector, avoiding any race conditions or need for synchronization.

This generated wrapper code gets inserted into each code object's support code section. Then during the final code generation phase , these implementations are registered:

rand_impls = DEFAULT_FUNCTIONS["rand"].implementations
rand_impls.add_dynamic_implementation(
    CPPStandaloneCodeObject,
    code=lambda owner: generate_rand_code("rand", owner),
    namespace=lambda owner: {},
    name="_rand",
)

randn_impls = DEFAULT_FUNCTIONS["randn"].implementations
randn_impls.add_dynamic_implementation(
    CPPStandaloneCodeObject,
    code=lambda owner: generate_rand_code("randn", owner),
    namespace=lambda owner: {},
    name="_randn",
)

The add_dynamic_implementation method means the code gets generated fresh for each code object, calling the lambda function each time to get the current preferences.

Step 4: Actual Code Generation and Compilation

Let's see what happens with the same user code we used before:

from brian2 import *
set_device('cpp_standalone')

G = NeuronGroup(100, 'v : 1')
G.v = 'rand()'
run(1*ms)

Brian2 analyzes this and generates several C++ files in the output directory (or wherever you've specified). For the state updater that handles G.v = 'rand()', it generates a file like stateupdater_codeobject.cpp that looks approximately like this:

#include "objects.h"
#include "brianlib/common_math.h"
#include<cmath>

namespace {
    // The wrapper from generate_rand_code() is inserted here
    inline double _rand(const int _vectorisation_idx) {
        return brian::_random_generators[0].rand();
    }
}

void _run_stateupdater_codeobject()
{
    using namespace brian;
    
    // Get the array pointer
    double* const _array_neurongroup_v = _ptr_array_neurongroup_v;
    const int _N = 100;
    
    // Loop over all neurons
    for(int _idx=0; _idx<_N; _idx++)
    {
        // Call the wrapper which calls the RandomGenerator
        const double _v_value = _rand(_idx);
        
        // Assign to the array
        _array_neurongroup_v[_idx] = _v_value;
    }
}

All these generated files then get compiled together. The build system (typically using make with a generated Makefile) compiles each .cpp file to an object file, then links them all together into an executable, typically called main or main.exe.

Step 5: Execution

When you run this executable, the main() function (from the generated main.cpp) does this sequence:

int main(int argc, char **argv)
{
    // ... command line argument parsing ...
    
    brian_start();  // Initialization phase
    
    {
        using namespace brian;
        // ... various run functions get called here ...
        _run_stateupdater_codeobject();  // This is where our code runs
    }
    
    brian_end();  // Cleanup phase
    return 0;
}

The brian_start() function calls _init_arrays(), which as we saw earlier creates the RandomGenerator objects. Let's trace through what happens when the actual random numbers get generated:

The first neuron (iteration zero of the loop) calls _rand(0). This calls brian::_random_generators[0].rand(), which executes the rand() method on the first (and only) generator. Inside rand(), it calls gen() twice. The first call to gen() might return something like 2357136044 (some random 32-bit value). Shifting right by five gives 73660501. The second call might return 2546248239, and shifting right by six gives 39784503. Then it computes (73660501 * 67108864.0 + 39784503) / 9007199254740992.0, which equals approximately 0.550055. This gets stored in _array_neurongroup_v[0].

The second neuron calls _rand(1), which again calls gen() twice with the next state of the Mersenne Twister, performs the same mathematical operations, and returns a different random number, maybe 0.715189, which gets stored in _array_neurongroup_v[1].

This continues for all one hundred neurons, with each call advancing the state of the single Mersenne Twister generator.

The Problem This Creates

Now you can see the fundamental issue. Let's trace through what happens if you set a seed and run the same simulation with both backends.

Suppose you set seed(12345) and then create a neuron group with G.v = 'rand()' to set one hundred voltages.

With the runtime (Cython) backend: The seed(12345) call executes np.random.seed(12345) and resets the buffer indices to zero. When the code object runs and the first neuron calls _rand(0), the index is zero so it refills the buffer by calling numpy.random.rand(20000). With NumPy's random number generator seeded at 12345, this produces twenty thousand numbers, let's say starting with [0.929616, 0.318763, 0.980764, ...]. The first neuron gets 0.929616, the second gets 0.318763, the third gets 0.980764, and so on.

With the C++ standalone backend: When the generated executable runs, _init_arrays() creates a RandomGenerator and you'd need to explicitly seed it by calling the generated seed function with 12345. This seeds the std::mt19937 generator with 12345. Even though this is also a Mersenne Twister implementation, the initial state will be different because the seeding algorithm differs between implementations. When the first neuron calls _rand(0), it calls the RandomGenerator::rand() method, which calls gen() twice and applies the Jean-Sebastien Roy formula. The first call to gen() with the C++ implementation seeded at 12345 might return something like 3992670690. Shifted right by five, that's 124770959. The second call might return 3823185381, shifted right by six gives 59737271. The formula gives (124770959 * 67108864.0 + 59737271) / 9007199254740992.0, which is approximately 0.933024.

Even though both backends used seed 12345, the first random number is different: 0.929616 versus 0.933024. The sequences diverge immediately and remain different for all subsequent numbers.

This happens for multiple reasons. First, NumPy and the C++ standard library use different seeding algorithms to initialize their internal states from the seed value. Second, even if the internal states were identical, NumPy's rand() function and our Jean-Sebastien Roy implementation extract and combine the bits differently. Third, the order of operations differs: NumPy generates twenty thousand numbers at once from consecutive states, while C++ standalone generates each number individually.

This creates real practical problems. Suppose you're developing a complex network model and notice that one particular neuron exhibits interesting dynamics. You want to debug this to understand why. If you develop and test in runtime mode (because it's faster for iteration), then switch to standalone mode for the final production run (because it's faster overall), the random numbers will all be different, and that interesting neuron might not exhibit the same behavior anymore. Your debugging session becomes useless because you're essentially running a different simulation.

What We Want to Achieve

The goal of this issue is to make both backends use exactly the same random number generation code so that identical seeds produce identical sequences regardless of which backend you use. The path forward is clear: make the Cython runtime backend use the same RandomGenerator class that C++ standalone mode already has.

Since Cython can compile C++ code directly, once we extract the RandomGenerator class into a standalone header file, we can include that header in the generated Cython code. The Cython code will instantiate RandomGenerator objects (one per simulation, since runtime mode doesn't use threading) and call their rand() and randn() methods directly.

This eliminates all the sources of divergence. Both backends will use std::mt19937 with identical seeding, both will use the identical Jean-Sebastien Roy bit manipulation, and both will generate numbers in exactly the same order. A simulation with seed 12345 will produce identical results whether you run it in runtime mode or standalone mode.

There's a potential performance question here. The current Cython implementation with NumPy buffering was designed to minimize overhead from repeatedly calling Python code. However, in practice this might not matter much. The Cython implementation we'll create will be calling C++ methods directly from compiled Cython code, so there's essentially no Python overhead. The RandomGenerator::rand() method is quite fast (just two generator calls and some arithmetic), and modern compilers will likely inline it. We might even see comparable or better performance because we're eliminating the complexity of buffer management and memory allocation that the current system does.

Additionally, the current buffering approach has a subtle inefficiency: when you ask for twenty thousand random numbers from NumPy, NumPy actually has to generate those numbers sequentially anyway. It's not like it can magically create them all in parallel. So the buffering really just shifts when the work happens (all at once at buffer refill time) rather than reducing the total work. By eliminating the buffer, we're just spreading that work evenly across all the calls instead of bunching it up.

---

Let me know if you have any questions about the architecture, the specific files you'll be modifying, the testing process, or anything else! Feel free to ask about any part of the code that's unclear. And cc: @mstimberg in case I made any mistakes with the explanation 😅 or if there are additional details that would be helpful.

Also, if you don't get the time to work on this, please don't worry at all—just ping me and I'd be happy to take this on. We really appreciate your interest in contributing!

[Legend101Zz]

Hi @itu-itis21-nguyen21 just wanted to check in to see if you had any chance to look into the issue :)

[Legend101Zz]

@mstimberg @itu-itis21-nguyen21 as there's been no response yet , I'll take a dig at the issue and start working on it :)