Per-synapse RNGs

[kernfel]

Working on brian-team/brian2genn#142, I noticed that not all instances of $(gennrand_**) are fully translated. I understand that translation is omitted entirely in support code snippets, which I guess may have to do with ambiguity about whether it runs on host or device. However, translation is awkwardly partial in the following case, where it seems like it should be fully supported.

Brian2 code, from a failing test:

source = NeuronGroup(10, '', threshold='True')
target = NeuronGroup(10, 'dv/dt = 0/second : 1 (unless refractory)',
                     threshold='i>=5', refractory=defaultclock.dt)
S = Synapses(source, target, on_pre='v += rand()')
S.connect(j='i')

Brian2GeNN translation into GeNN code, cutting to the relevant chunk:

class synapsesWEIGHTUPDATE : public WeightUpdateModels::Base
{
public:
    DECLARE_MODEL(synapsesWEIGHTUPDATE, 0, 0);

    SET_SIM_CODE("$(addToInSyn,(int($(not_refractory_post)) * $(gennrand_uniform)));");
// snip
}

GeNN fails during compilation as follows:

/tmp/genn.sboZYgvY/magicnetwork_model_CODE/synapseUpdateCUDAOptim.cc(62): error: identifier "rng" is undefined

/tmp/genn.sboZYgvY/magicnetwork_model_CODE/synapseUpdateCUDAOptim.cc(62): error: identifier "$" is undefined

2 errors detected in the compilation of "/tmp/genn.sboZYgvY/magicnetwork_model_CODE/synapseUpdateCUDAOptim.cc".
terminate called after throwing an instance of 'std::runtime_error'
  what():  optimizeBlockSize: NVCC failed
/home/felix/projects/lib/genn/bin/genn-buildmodel.sh: line 101: 2614840 Aborted                 (core dumped) "$GENERATOR" "$BASEDIR/../" "$OUT_PATH" "$FORCE_REBUILD"

subprocess.CalledProcessError: Command '['/home/felix/projects/lib/genn/bin/genn-buildmodel.sh', '-i', '/home/felix/projects/brian2genn_DEV:/home/felix/projects/brian2genn_DEV/GeNNworkspace:/home/felix/projects/brian2genn_DEV/GeNNworkspace/brianlib/randomkit', 'magicnetwork_model.cpp']' returned non-zero exit status 50.

... and finally, the offending line in context:

                // loop through all incoming spikes
                for (unsigned int j = 0; j < numSpikesInBlock; j++) {
                    // only work on existing neurons
                    if (lid < group->rowStride) {
                        using namespace PresynapticUpdateSupportCode0;
                        const unsigned int synAddress = (shSpk[j] * group->rowStride) + lid;
                        const unsigned int npost = shRowLength[j];
                        if (lid < npost) {
                            const unsigned int ipost = group->ind[synAddress];
                            shLg[ipost] += (int(group->not_refractoryPost[ipost]) * curand_uniform_double($(rng)));}
                    }
                }

Note, this is not an urgent issue on my end, my immediate use case is covered.

[neworderofjamie]

Ahh, I had forgotten about this....GeNN doesn't currently support random number generation in any type of per-synapse code (proper error handling to report this in the SynapseGroup constructor would be very nice). The simplest and most efficient approach would be to have a XORWOW RNG per thread (this is the RNG we use in neuron kernels) but, the code you pasted illustrates the issue with this approach which is that, with postsynaptic parallelism, each thread processes an unordered list of spikes so you would not get repeatable results across multiple simulation runs. If this ok (@tnowotny thoughts please!) I can totally implement it this way.

The less performant and easy alternatives would be:

• XORWOW RNG per-synapse - this would be something of a memory/memory bandwidth hog as the cuRAND implementation has 44 bytes of state (including caching of box muller coefficients etc)
• Single Phillox counting RNG - we currently use these for initialisation and procedural connectivity as you only need one per model and can 'skip' through streams in constant time. However, here, we would need to further divide the streams between timesteps and I think it would get pretty complex pretty quickly.

What does Brian2GeNN currently do by the way?

[tnowotny]

I think non-reproducible runs are acceptable but if possible we should attach some warnings about it in several places. Even without the explicit breakage the RNG issue implemented with one RNG per thread will introduce, in more subtle cases the reproducibility is already broken simply by spikes being processed in unpredictable order.

[kernfel]

Brian2GeNN currently uses a dedicated RNG per synapse:

                for (unsigned int j = 0; j < numSpikesInBlock; j++) {
                    // only work on existing neurons
                    if (lid < group->rowStride) {
                        using namespace PresynapticUpdateSupportCode0;
                        const unsigned int synAddress = (shSpk[j] * group->rowStride) + lid;
                        const unsigned int npost = shRowLength[j];
                        if (lid < npost) {
                            const unsigned int ipost = group->ind[synAddress];
                            shLg[ipost] += (int(group->not_refractoryPost[ipost]) * _rand(group->_seed[synAddress]));}
                    }
                }

However, the RNG state here is a single uint64, so considerably lighter than cuRAND.

I guess Brian users who absolutely need reproducibility always have the deterministic, single-threaded backends to fall back to, but on the other hand, I can just keep synaptic random number generation on the old RNG to avoid any arguments for now. We can always re-address this if/when you decide to explicitly support synaptic RNG in GeNN.

[neworderofjamie]

It would definitely be a good start to switch over the neuron RPGs. I will try and implement the synaptic RNGs relatively soon though. I vaguely remember there being another issue with Brian requiring an implementation of the Binomial distribution...Is that still true?

[kernfel]

Yeah, that's right. PoissonInput objects need this; it's currently implemented in support code, so I can't easily port it to the GeNN RNG either.

[neworderofjamie]

Give me a second - I'll just try and rebase my last attempt two years forward 😨

[kernfel]

No need to hurry. The model that gave the initial impetus for this effort is covered with just the per-neuron RNGs. Everything else is nice-to-have, takes work away from any potential GSoC students ;)

[neworderofjamie]

No worries - reducing the number of excess branches is always a good thing. Have a go with #498

[tnowotny]

Ok, folks, have run into this issue independently doing my own naive thing in pyGeNN. Maybe, when time allows, doing the per-thread RNG would still be a good thing?