User-provided function in run_regularly

Support
,
1

Hi Marcel

I have problem compiling the following example code in standalone mode. Any hints would be appreciated.

  from brian2 import *

  set_device('cpp_standalone', directory='standalone')

  @implementation('cpp', '''
       double spike_rate(t, i) {
          return ta2d(t, i);
       }
       ''')
  @check_units(t=second, i=1, result=1)
  def spike_rate(t, i):
      raise NotImplementedError('use standalone mode')

  #

  N = 6
  input_neurons = NeuronGroup(N, 'rate : Hz', threshold='rand()<rate*dt')
  ta2d = TimedArray([[1, 2, 3, 4, 5, 6], [100, 200, 300, 400, 500, 600]], dt=0.2*ms)
  #input_neurons.run_regularly('rate=ta2d(t, i)*Hz', dt=0.2*ms)
  input_neurons.run_regularly('rate=spike_rate(t, i)*Hz', dt=0.2*ms)
  run(0.3*ms, report='text')
  print(input_neurons.rate[:])

The error is:

    run(0.3*ms, report='text')
    File "/home/wxie/local_pkgs/brian2/brian2/units/fundamentalunits.py", line 2428, in new_f
      result = f(*args, **kwds)
    File "/home/wxie/local_pkgs/brian2/brian2/core/magic.py", line 377, in run
      return magic_network.run(duration, report=report, report_period=report_period,
    File "/home/wxie/local_pkgs/brian2/brian2/core/magic.py", line 230, in run
      Network.run(self, duration, report=report, report_period=report_period,
    File "/home/wxie/local_pkgs/brian2/brian2/core/base.py", line 291, in device_override_decorated_function
      return getattr(curdev, name)(*args, **kwds)
    File "/home/wxie/local_pkgs/brian2/brian2/devices/cpp_standalone/device.py", line 1557, in network_run
      self.build(direct_call=False, **self.build_options)
    File "/home/wxie/local_pkgs/brian2/brian2/devices/cpp_standalone/device.py", line 1290, in build
      self.compile_source(directory, compiler, debug, clean)
    File "/home/wxie/local_pkgs/brian2/brian2/devices/cpp_standalone/device.py", line 1031, in compile_source
      raise RuntimeError(error_message)
    RuntimeError: Project compilation failed (error code: 512). Consider running with "clean=True" to force a complete rebuild.
2

Hi @DavidKing2020 . If you are getting compilation errors, it can be helpful to add debug=True to the set_device call. It will then print out the compiler error message. In your case, the function is not using correct C++ syntax, since the arguments of the function are missing the types. It should read double spike_rate(double t, int i). But after fixing this, it will still complain about ta2d, which is not defined in the C++ code. To make the ta2d function available here, you need to declare it as a dependency of the function: Functions ā€” Brian 2 2.5.1 documentation
I.e., in your example:

ta2d = TimedArray([[1, 2, 3, 4, 5, 6], [100, 200, 300, 400, 500, 600]], dt=0.2*ms)

@implementation('cpp', '''
   double spike_rate(double t, int i) {
      return ta2d(t, i);
   }
   ''', dependencies={'ta2d': ta2d})

Of course, you could directly use ta2d instead of spike_rate in the equations here, but I guess this was just an example for testing.

3

Thanks Marcel. I tried to update the code to using external source according to:

https://brian2.readthedocs.io/en/stable/advanced/functions.html#external-source-files

  @implementation('cpp',
          sources=[os.path.join(current_dir,'spike_rate.cpp')],
          headers=['"spike_rate.h"'], include_dirs=[current_dir],
          dependencies={'ta2d': ta2d})
  @check_units(t=second, i=1, result=1)
  def spike_rate(t, i):
      raise NotImplementedError('use standalone mode')

and the cpp file is:

double spike_rate(double t, int i) {
   return ta2d(t, i);
}

The issue here is that ā€œā€˜ta2dā€™ was not declared in this scopeā€. What is the best way to declare the python function ta2d(t, i) to cpp code?

4

Hmm, Iā€™m afraid it is impossible to use the ta2d function from an external source file. The idea for external source files is to use them for code that already exists (e.g. algorithms), i.e. for code that is independent of the Brian simulation. Here, the file would need an additional #include and also change the way it calls the function (which is not called ta2d internally). If you want to keep the C++ code in a separate file just to organize things more nicely, Iā€™d suggest to simply use Python to load the file, and then plug it into the usual function implementation definition:

with open(os.path.join(current_dir,'spike_rate.cpp')) as f:
    @implementation('cpp', f.read(),
                    dependencies={'ta2d': ta2d})
    @check_units(t=second, i=1, result=1)
    def spike_rate(t, i):
      raise NotImplementedError('use standalone mode')
5

Thanks again. Iā€™m trying to add the spikemon into example code as follows:

 spikemon = SpikeMonitor(input_neurons)

 ta2d = TimedArray([[1, 2, 3, 4, 5, 6], [100, 200, 300, 400, 500, 600]], dt=DT)

 with open(os.path.join(current_dir,'spike_rate.cpp')) as f:
     @implementation('cpp',
             f.read(),
             dependencies={'ta2d': ta2d, 'spikemon':spikemon})
     @check_units(t=second, i=1, result=1)
     def spike_rate(t, i):
         raise NotImplementedError('use standalone mode')

but more the compilation shows the following errors:

      Traceback (most recent call last):
        File "/home/wxie/local_pkgs/brian2/brian2/groups/group.py", line 351, in state
          var = self.variables[name]
        File "/home/wxie/local_pkgs/brian2/brian2/core/variables.py", line 1453, in __getitem__
          return self._variables[item]
      KeyError: 'implementations'

The full code and the error message are attached. Further help is appreciated.
err.txt (4.2 KB)
example.py (802 Bytes)
spike_rate.cpp.txt (127 Bytes)

6

Hi @DavidKing2020 . A SpikeMonitor is not a function, so it cannot be added as a dependency. This mechanism is not meant as a general interface to arbitrary Brian objects. Since standalone mode uses global variables for all group variables (including variables of a SpikeMonitor), you can access them from functions using their name if necessary. See this discussino for more information on how to do that: User-defined functions - #4 by mdepitta
That said, maybe you could explain what you want to achieve? Iā€™d try to implement as much as possible with Brianā€™s standard mechanisms, and only use hand-written C++ code if absolutely necessary.

7

Hi Marcel

Iā€™m trying to implement the suggestion you made in the following thread:

Re-run the model in a previous time segment during a long run

 Quote: " If you want to do it in C++ standalone mode, then the only way I can think of is via C++ code that 
  implements the logic. Instead of `rate = ta_img(t, i)*Hz` youā€™d use something like `rate = my_func(t, 
  i)*Hz` and implement `my_func` in C++ so that it looks at the recent firing rate and decides whether 
  to present the next stimulus or the same stimulus with a higher firing rate."

In that thread you pointed to the same thread of discussion with @mdepitta. I can try to digest further by going through the discussion but it would be very helpful if could provide inputs to my specific issues.

8

Hi again. Apologies for not having remembered the problem you were trying to solve, I got lost in all the threadsā€¦ If you want to/need to go the manual route, hereā€™s the general approach: in your C++ code, you add placeholders for the variables that you need to access. Then, you replace these placeholders by the name as provided by device.get_array_name. E.g., to make your previous example print out the spike count for each neuron, you could use:

// spike_rate.cpp
double spike_rate(double t, int i) {
   std::cout<<"t = " << t << "s" <<std::endl;
   std::cout<<"spike count for neuron " << i << ": " <<brian::%%SPIKECOUNT%%[i]<<std::endl;
   return ta2d(t, i);
}

and then in your Python code:

import os

from brian2 import *

current_dir = os.path.dirname(__file__)
set_device('cpp_standalone', directory='standalone', debug=True)

ta2d = TimedArray([[1, 2, 3, 4, 5, 6], [100, 200, 300, 400, 500, 600]], dt=0.2*ms)

N = 6
input_neurons = NeuronGroup(N, 'rate : Hz', threshold='rand()<rate*dt')
spikemon = SpikeMonitor(input_neurons)

with open(os.path.join(current_dir,'spike_rate.cpp')) as f:
    spike_rate_code = f.read()
spike_rate_code = spike_rate_code.replace('%%SPIKECOUNT%%',
                                          device.get_array_name(spikemon.variables['count']))

@implementation('cpp', spike_rate_code,
                dependencies={'ta2d': ta2d})
@check_units(t=second, i=1, result=1)
def spike_rate(t, i):
    raise NotImplementedError('use standalone mode')


input_neurons.run_regularly('rate=spike_rate(t, i)*Hz', dt=0.2*ms)
run(0.3*ms, report='text')
print(input_neurons.rate[:])

Hope that makes things clearer.

9

Thanks Marcel. Here are two follow up questions:

  1. if I replace

     device.get_array_name(spikemon.variables['count'])

by

  device.get_array_name(spikemon.variables['num_spikes'])

I got the following errors:

device.get_array_name(spikemon.variables['num_spikes']))
File "/home/wxie/local_pkgs/anaconda3/envs/brian2/lib/python3.8/site-packages/brian2/core/variables.py", line 1440, in __getitem__
return self._variables[item]
KeyError: 'num_spikes'

Could you provide further clue on this?

There are also methods like all_values(), is it possible to access to those method using ā€œreplaceā€?

  1. The function spike_rate(double t, int i) process information for a single neuron at a certain instant in time. I did the following to check the number of spikes from all neurons at a certain time, e.g. a 0.4*ms

     std::cout<<"count for "<<i<<": "<<brian::%%SPIKECOUNT%%[i]<<std::endl;
 int tot = 0;
 static bool done = false;
 if(t==0.0004 && done == false) {
   for(int m=0; m<6; m++) {
     tot += brian::%%SPIKECOUNT%%[i];
   }
   done = true;
 }
 std::cout<<" tot = "<<tot<<endl;

This seems rather inefficient. Is there a better way of doing it.

10

Hi. You can have a look at the variables dictionary to see the internal variables. Functions like all_values are written in Python, you cannot access them from C++. The num_spikes attribute is just a more explicit name for the internal variable N. I donā€™t quite understand what you try to achieve with the check of the count at a specific time point ā€“ why do you call that function at every dt, then? Note that checks like if t == 0.0004 are also very fragile due to the way floating point numbers work. You might want to call the function with the integer value t_in_timesteps instead.

A few more general remarks:
If you want to do complex stuff to calculate the rates in PoissonGroup, then it might make sense to substitute it with an equivalent NeuronGroup. Instead of PoissonGroup(N, rates='something'), you can instead use NeuronGroup(N, 'rates = something : Hz', threshold='rand() < rates*dt'). The advantage is that you can then easily add new variables to the model, which can make things a lot more readable. You might even get around the use of C++ functions in the first place if you lay out things nicely. Iā€™ll try to come up with a quick example of what I mean.

11

Sorry, just realized that you already replaced your PoissonGroup with a NeuronGroupā€¦

12

I just posted an example showing what I had in mind: Changing stimuli programmatically in C++ standalone
I hope this gives you an idea what you can do with Brian code alone ā€“ if you need something more complex, you can still augment it by a C++ function of course.

1 Like
13

This is amazing and just shows how flexible brian2 is in this case. Thanks for all the help!

2 Likes