# Possibility of speeding up Poisson spike generation

**URL:** <https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342>\
**Category:** Support\
**Tags:** performance, poisson\
**Created:** [5 December 2024 16:52 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342 "2024-12-05T16:52:05Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![DavidKing2020](https://avatars.discourse-cdn.com/v4/letter/d/e0b2c6/32.png) [@DavidKing2020](https://brian.discourse.group/u/DavidKing2020)\
**Post date:** [5 December 2024 16:52 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/1 "2024-12-05T16:52:05Z")

</div>

I have a layer of input neurons where each neuron generates Poisson spikes based on the intensity of a specific pixel of an image. The profiling with running in cpp\_standalone mode seems to show this procedure takes most of the running time, i.e.

```
  Profiling summary
  =================
  input_neurons_spike_thresholder_codeobject_2 15.25 ms 66.87 %

```

Is it possible to speed this up in some way?

The python code is as follows:

```
    self.input_neurons = brian2.NeuronGroup(input_param.Num_input_neuron, 
                                     '''x: 1
                                        y: 1
                                        switch_stimulus: boolean
                                        correct_identification: boolean
                                        higher_than_min_spike_rate: boolean
                                        rate: Hz
                                        label: integer
                                        start_t : second # Start time of the current trial
                                        stimulus_idx : integer # Index of the stimulus to show
                                        stimulus_strength : 1 # Factor to scale up stimulus
                                        repetitions : integer # Number of times the stimulus has been presented
                                        ''', threshold='rand()<rate*dt', name='input_neurons',  
                                        namespace = {'rows':input_param.input_rows, 'cols':input_param.input_cols})
    
    self.input_neurons.x = 'i%cols'  
    self.input_neurons.y = 'rows - i//rows' # // is a floored division

```

The corresponding c++ code is below:

```
void _run_input_neurons_spike_thresholder_codeobject_2()
{
    using namespace brian;

    const std::clock_t _start_time = std::clock();

    ///// CONSTANTS ///////////
    const int64_t N = 784;
    const size_t _num_spikespace = 785;
    const size_t _numdt = 1;
    const size_t _numrate = 784;
    ///// POINTERS ////////////

    int32_t* __restrict _ptr_array_input_neurons__ spikespace = _array_input_neurons__spikespace;
    double* _ptr_array_defaultclock_dt = _array_defaultclock_dt;
    double* __restrict _ptr_array_input_neurons_rate = _array_input_neurons_rate;

    //// MAIN CODE ////////////
    // scalar code
    const int _vectorisation_idx = -1;

    const double dt = _ptr_array_defaultclock_dt[0];

    long _count = 0;
    for(size_t _idx=0; _idx<N; _idx++)
    {
        const size_t _vectorisation_idx = _idx;

        const double rate = _ptr_array_input_neurons_rate[_idx];
        const char _cond = _rand(_vectorisation_idx) < (dt * rate);

        if(_cond) {
            _ptr_array_input_neurons__spikespace[_count++] = _idx;
        }
    }
    _ptr_array_input_neurons__spikespace[N] = _count;

    const double _run_time = (double)(std::clock() -_start_time)/CLOCKS_PER_SEC;
    input_neurons_spike_thresholder_codeobject_2_profiling_info += _run_time;
}

```

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [5 December 2024 17:11 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/2 "2024-12-05T17:11:03Z")

</div>

Hi @DavidKing2020. I’ve recently pushed a big improvement for the C++ random number generation into the development version of Brian2 (see [Use Mersenne Twister from C++ standard library by mstimberg · Pull Request #1559 · brian-team/brian2 · GitHub](https://github.com/brian-team/brian2/pull/1559)), this should accelerate `PoissonGroup` quite a bit. Maybe you could try testing with this version (see [Installation — Brian 2 2.7.1 documentation](https://brian2.readthedocs.io/en/stable/introduction/install.html#development-install))? I hope to do a new release before the end of the year in either case.

---

<div class="post-metadata">

**Author:** ![DavidKing2020](https://avatars.discourse-cdn.com/v4/letter/d/e0b2c6/32.png) [@DavidKing2020](https://brian.discourse.group/u/DavidKing2020)\
**Post date:** [5 December 2024 20:21 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/3 "2024-12-05T20:21:42Z")

</div>

Thanks @mstimberg. It speeds it up by almost a factor of 2.

```
 input_neurons_spike_thresholder_codeobject 9.04 ms 51.79 %

```

It still takes most of the running time in my case. Could there be other waysto further speed it up?

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [6 December 2024 21:01 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/4 "2024-12-06T21:01:21Z")

</div>

I don’t see any easy way to speed things up more than that for now. But if your whole simulation already takes 20ms, how much more performance do you need 🙂 ?

---

<div class="post-metadata">

**Author:** ![DavidKing2020](https://avatars.discourse-cdn.com/v4/letter/d/e0b2c6/32.png) [@DavidKing2020](https://brian.discourse.group/u/DavidKing2020)\
**Post date:** [6 December 2024 21:45 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/5 "2024-12-06T21:45:51Z")

</div>

Hi @mstimberg, it could be significant when running through large sample of images with many epochs, but a factor two reduction is already a huge help. Looking forward to the next release 👍

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [7 December 2024 11:34 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/6 "2024-12-07T11:34:47Z")

</div>

I understand. Actually, if the rates are constant during the run, there is potentially a way to speed up the Poisson spike generation further. I’ll look into this on Monday and report back here.

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [9 December 2024 14:04 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/7 "2024-12-09T14:04:35Z")

</div>

Hi again @DavidKing2020. So, here’s what I had in mind: our standard `rand() < rates*dt` formulation is very flexible, since it applies even when the rates are changing over time, but it needs to generate a random number for each neuron at every timestep. If your rates are constant (within a single run), then you can use the following formulation instead, which makes use of the fact that the times between the spikes of a Poisson process follow an exponential distribution:

```python
g = NeuronGroup(1000,
                '''
                rates : Hz (constant)
                next_spike_timestep : integer
                ''',
                threshold='t_in_timesteps >= next_spike_timestep',
                reset='''
                next_spike_timestep = t_in_timesteps + timestep(-log(rand())/rates, dt)
                ''')
g.rates = 10*Hz
g.next_spike_timestep = "timestep(-log(rand())/rates, dt)"
M = SpikeMonitor(g)

```

For low rates, this should be much more efficient since it only needs to generate a single random number per spike.  
Hope that works for you!

---

<div class="post-metadata">

**Author:** ![DavidKing2020](https://avatars.discourse-cdn.com/v4/letter/d/e0b2c6/32.png) [@DavidKing2020](https://brian.discourse.group/u/DavidKing2020)\
**Post date:** [9 December 2024 14:31 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/8 "2024-12-09T14:31:51Z")

</div>

Thanks @mstimberg. Just to make sure I got it right. In my case, the spike rate for each neuron in the input layer is proportional to the brightness of the corresponding pixel in an image. That means that for each neuron, the spike rate is unique but is constant in each run. Will the new method apply to this scenario?

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [9 December 2024 14:39 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/9 "2024-12-09T14:39:09Z")

</div>

Hi @DavidKing2020. Yes, this is the situation it was meant for: the rate can be different for each neuron, but it needs to stay constant (for the run).

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [9 December 2024 14:40 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/10 "2024-12-09T14:40:08Z")

</div>

In your use case, `g.rates = 10*Hz` would instead set a different rate for each neuron, but nothing else would need to change.

---

<div class="post-metadata">

**Author:** ![DavidKing2020](https://avatars.discourse-cdn.com/v4/letter/d/e0b2c6/32.png) [@DavidKing2020](https://brian.discourse.group/u/DavidKing2020)\
**Post date:** [9 December 2024 21:01 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/11 "2024-12-09T21:01:41Z")

</div>

Hi @mstimberg, I tried to implement the new method, The “input\_neurons\_spike\_thresholder\_codeobject” now cost much less time, but it somehow leads to different processes that consume more CPUs:

```
  s_input2e0_pre_plastic_codeobject 244.53 ms 59.08 %
  s_input2e0_pre_plastic_push_spikes 115.95 ms 28.01 %
  input_neurons_spike_resetter_codeobject 29.09 ms 7.03 %
  s_input2e0_post_codeobject 15.54 ms 3.75 %
  input_neurons_spike_thresholder_codeobject 2.43 ms 0.59 %

```

The following is the python code. I also attached the c++ code for the top 4 codes that consume the most CPUs. These codes did not rank high when using the original ‘rand()\<rate\*dt’ method.

```
    self.input_neurons = brian2.NeuronGroup(input_param.Num_input_neuron, 
                                            '''
                                            x: 1
                                            y: 1
                                            switch_stimulus: boolean
                                            correct_identification: boolean
                                            higher_than_min_spike_rate: boolean
                                            rate: Hz (constant)
                                            next_spike_timestep : integer
                                            label: integer
                                            start_t : second # Start time of the current trial
                                            stimulus_idx : integer # Index of the stimulus to show
                                            stimulus_strength : 1 # Factor to scale up stimulus
                                            repetitions : integer # Number of times the stimulus has been presented
                                            ''', 
                                            threshold='t_in_timesteps >= next_spike_timestep',
                                            reset='''
                                                  next_spike_timestep = t_in_timesteps + timestep(-log(rand())/rate, dt)
                                                  ''', name='input_neurons',  
                                            namespace = {'rows':input_param.input_rows, 'cols':input_param.input_cols})
    
    self.input_neurons.x = 'i%cols'  
    self.input_neurons.y = 'rows - i//rows' # // is a floored division
    self.input_neurons.next_spike_timestep = "timestep(-log(rand())/rate, dt)"
   arguments = {self.input_neurons.rate:scaled_img*brian2.Hz} 

```

This argument was feeded into ‘run\_args=arguments’

[input\_neurons\_spike\_resetter\_codeobject.cpp](https://brian.discourse.group/uploads/short-url/iysUTyxjTHXzrkgpt5cY0NDh2gI.cpp) (6.7 KB)  
[s\_input2e0\_post\_codeobject.cpp](https://brian.discourse.group/uploads/short-url/u9RFwiXr9GLJaeZG51GcxdHohXI.cpp) (12.3 KB)  
[s\_input2e0\_pre\_plastic\_codeobject.cpp](https://brian.discourse.group/uploads/short-url/lElpu040Hos6aHEH8JQtZD7O4T3.cpp) (10.3 KB)  
[s\_input2e0\_pre\_plastic\_push\_spikes.cpp](https://brian.discourse.group/uploads/short-url/962W3MM0UZ96GMxu8UIcsIcIaGF.cpp) (5.9 KB)

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [10 December 2024 12:27 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/12 "2024-12-10T12:27:53Z")

</div>

Hi @DavidKing2020. The slow times for the synapse-related code objects indicate that you have many more spikes in your network. I think the problem is the initialization of `next_spike_timestep`, since it refers to the value of `rate` that will only be set later with the `run_args` mechanism. Since the value if 0 by default, it does something wrong for the first spike – I think all your neurons spike in the first time step.  
We discuss this situation in the [documentation](https://brian2.readthedocs.io/en/stable/user/computation.html#running-multiple-simulations-with-different-parameters):

> a more general solution is to move the point where the `run_args` are applied. You can do this by calling the device’s [`apply_run_args`](https://brian2.readthedocs.io/en/stable/reference/brian2.devices.cpp_standalone.device.CPPStandaloneDevice.html#brian2.devices.cpp_standalone.device.CPPStandaloneDevice.apply_run_args) function

Could you try adding `device.apply_run_args()` before the initialization of `next_spike_timestep`? It should look like similar to this after the change:

```python
    self.input_neurons.x = 'i%cols'  
    self.input_neurons.y = 'rows - i//rows' # // is a floored division
    brian2.device.apply_run_args()
    self.input_neurons.next_spike_timestep = "timestep(-log(rand())/rate, dt)"

```

---

<div class="post-metadata">

**Author:** ![DavidKing2020](https://avatars.discourse-cdn.com/v4/letter/d/e0b2c6/32.png) [@DavidKing2020](https://brian.discourse.group/u/DavidKing2020)\
**Post date:** [10 December 2024 14:16 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/13 "2024-12-10T14:16:12Z")

</div>

Hi @mstimberg , the source of the slow timing is indeed because all input neuron spikes. I tried to implement “device.apply\_run\_args()” but it didn’t help. All input neurons somehow still spiks. I’ve attached the input\_neuron code , there are two of them and essentially the same:

```
 input_neurons_spike_resetter_codeobject_1.cpp
 input_neurons_spike_resetter_codeobject.cpp

```

[input\_neurons\_spike\_resetter\_codeobject.cpp](https://brian.discourse.group/uploads/short-url/iysUTyxjTHXzrkgpt5cY0NDh2gI.cpp) (6.7 KB)  
[input\_neurons\_spike\_resetter\_codeobject\_1.cpp](https://brian.discourse.group/uploads/short-url/AtgpRHNnEygyelyUSdZsBuQCp9t.cpp) (6.7 KB)

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [10 December 2024 14:19 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/14 "2024-12-10T14:19:03Z")

</div>

Hi @DavidKing2020, could you please share the generated `main.cpp` file as well? This is where the `apply_run_args` call would make a difference.

---

<div class="post-metadata">

**Author:** ![DavidKing2020](https://avatars.discourse-cdn.com/v4/letter/d/e0b2c6/32.png) [@DavidKing2020](https://brian.discourse.group/u/DavidKing2020)\
**Post date:** [10 December 2024 14:46 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/15 "2024-12-10T14:46:32Z")

</div>

Hi @mstimberg , here it is attached. Thanks  
[main.cpp](https://brian.discourse.group/uploads/short-url/ivULiwEoDT6cYoPcy3d4LPKLDIU.cpp) (15.1 KB)

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [10 December 2024 15:02 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/16 "2024-12-10T15:02:57Z")

</div>

Hmm, that looks correct, assuming that `code_objects/input_neurons_group_variable_set_conditional_codeobject_2.h` sets the `next_spike_timestep` variable as expected.

Could you maybe present a single sample for **`0*ms`** (i.e., it will only do the initializations and nothing else), and then print out the values of `input_neurons.rate` and `input_neurons.next_spike_timestep` to see whether they are as expected?

---

<div class="post-metadata">

**Author:** ![DavidKing2020](https://avatars.discourse-cdn.com/v4/letter/d/e0b2c6/32.png) [@DavidKing2020](https://brian.discourse.group/u/DavidKing2020)\
**Post date:** [10 December 2024 15:28 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/17 "2024-12-10T15:28:10Z")

</div>

I tried to run 0_ms and 10_ms. The outputs are attached. The “input\_neurons\_group\_variable\_set\_conditional\_codeobject\_2…cpp” are also attached. When running time is 0_ms, the next\_spike\_time look correct but when runtime = 10_ms, the neurons which should not spike has lots of 99 values. These 99 value change to larger value with longer running time. Looks like a initialization issue. Thanks  
[input\_neurons\_group\_variable\_set\_conditional\_codeobject\_2.cpp](https://brian.discourse.group/uploads/short-url/388zq55CrhnGWTOEO3AwsqKlaXg.cpp) (6.0 KB)  
[out\_runTime\_10ms.txt](https://brian.discourse.group/uploads/short-url/fuyScZbKEhOwLZWYDyz1BLgeU6z.txt) (12.5 KB)  
[out\_runTime\_0ms.txt](https://brian.discourse.group/uploads/short-url/pOT0lOLkt3tMm8UGBnTggIlfjmC.txt) (12.0 KB)  
[input\_neurons\_group\_variable\_set\_conditional\_codeobject\_2.cpp](https://brian.discourse.group/uploads/short-url/388zq55CrhnGWTOEO3AwsqKlaXg.cpp) (6.0 KB)

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [10 December 2024 15:36 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/18 "2024-12-10T15:36:41Z")

</div>

Oh sorry, I just realize something: in your case you have neurons that have rate 0, right? The current formulation does not work for these neurons…  
Could you try replacing

```python
    brian2.device.apply_run_args()
    self.input_neurons.next_spike_timestep = "timestep(-log(rand())/rate, dt)"

```

by

```python
    self.input_neurons.next_spike_timestep = 1_000_000
    brian2.device.apply_run_args()
    self.input_neurons.next_spike_timestep['rate>0*Hz'] = "timestep(-log(rand())/rate, dt)"

```

The `1_000_000` value for the neurons that have rate 0 is arbitrary – it should just be a value bigger than the number of timesteps in a single run, so it never gets triggered.

---

<div class="post-metadata">

**Author:** ![mstimberg](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mstimberg/32/11_2.png) [@mstimberg](https://brian.discourse.group/u/mstimberg)\
**Post date:** [10 December 2024 15:37 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/19 "2024-12-10T15:37:27Z")

</div>

I think with the previous formulation, neurons that should never spike (rate=0Hz), spiked at every timestep…

---

<div class="post-metadata">

**Author:** ![DavidKing2020](https://avatars.discourse-cdn.com/v4/letter/d/e0b2c6/32.png) [@DavidKing2020](https://brian.discourse.group/u/DavidKing2020)\
**Post date:** [10 December 2024 15:54 UTC](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342/20 "2024-12-10T15:54:30Z")

</div>

The issue for “spike at every timestep” is gone :smiley, but the overall running time is similar as the the orignal ‘rand()\<rate\*dt’ method. The reason is from “s\_input2e0\_post\_codeobject” as shown in the following profile:

```
 s_input2e0_post_codeobject 10.75 ms 45.14 %
 input_neurons_spike_thresholder_codeobject 2.64 ms 11.10 %
 s_input2e0_pre_plastic_codeobject 1.61 ms 6.77 %
 exci_neurons0_stateupdater_codeobject 1.37 ms 5.77 %
 s_input2e0_pre_plastic_push_spikes 1.28 ms 5.37 %

```

The “s\_input2e0\_post\_codeobject .cpp” is attached. With the ‘rand()\<rate\*dt’ method, “s\_input2e0\_post\_codeobject .cpp” didn’t take as much time. The code is attached. Thanks  
[s\_input2e0\_post\_codeobject.cpp](https://brian.discourse.group/uploads/short-url/u9RFwiXr9GLJaeZG51GcxdHohXI.cpp) (12.3 KB)

[Next page](https://brian.discourse.group/t/possibility-of-speeding-up-poisson-spike-generation/1342.md?page=2)