User-defined functions

Support
1

My code target is cpp_standalone, and I am trying to implement some function with if/else and for/while statements that modify variables inside my NeuronGroup. Two issues come up. First, it seems that for some reasons, Brian is not happy with the mere declaration of

@implementation('cpp', '''
 void myfunction(){ ...})

But also asks for a specification of that function in Python and the need to then pre-pend the @check_units() decorator, which in my case is not necessary. Second, in CPP, the only way I can easily pass modifiable variables is to pass them by pointer and declare myfunction of void type. But this is not possible if I need to match my CPP function with a Python declaration. How can I work it around?

This is the basic code:

pars = {'q_PD': 100, 'q_DP': 10.}

code = '''
       tau        : 1 
       rval       : 1
       fraction   : 1
       X_O        : integer
       p_O        : 1
       transition = tau<Q*dt/second : boolean
       Q = X_O*q_DP+(N_r-X_O)*q_PD : 1
       '''

cell = NeuronGroup(1,
             code,namespace=pars,
             events={'transitioned': 'transition', 'no_transition': 'not transition'})
cell.tau = '-log(rand())'

# "Transition" code (q_DP is a parameter in the namespace)
code_upon_transition = """
                       rval = rand()
                       p_O = q_DP/Q
                       X_O = X_O + (rval<=X_O*p_O)*1 - (rval>X_O*p_O)*1
                       fraction = 1 - tau/(Q*dt/second)
                       nMGA(dt/second,Q,p_O,rval,tau,fraction,X_O) <----- This is the function I want to implement
                       """
code_upon_notransition = """
                         tau -= Q*dt/second
                         """
# Run on events
astro_compartment.run_on_event('transitioned',code_upon_transition)
astro_compartment.run_on_event('notransition',code_upon_notransition)

As for my user-defined function:

@implementation('cpp', '''
 void nMGA(double dt_,double Q,double p_O,double *rval,double *tau,double *fraction,double *X_O) {
    while(true):
        tau = -log((double) (rand() / (RAND_MAX)));
        if(tau>=fraction*Q*dt_){
            tau -= fraction*Q*dt_;
            break;
        }
        else {
            rval = (double) (rand() / (RAND_MAX));
            X_O += int(rval<=X_O*p_O) - int(rval>X_O*p_O);
            fraction -= tau/(Q*dt_);
        {        
 }
 ''')

If I run the code above, what I get is an error of the type:

ValueError: The Python function "AstroCell" does not specify how it deals with units, need to specify "arg_units" or use the "@check_units" decorator.

Where AstroCell is the function that contains the definition of my above cell NeuronGroup. I believe that this is because Brian’s parser cannot identify nMGA in the scope. But then, from what I know, the only way to declare a Python function that does something similar to my nMGA in C++ would be to have something like

def nMGA(dt_,Q,p_O,rval,tau,fraction,X_O):
    ...
    return rval,tau,fraction,X_O
2

I am elaborating a bit more actually on the issue. I experimented a bit with Python functions. There is this ambiguity in the language that when I define a function, if the input arguments are mutable objects, then it is like they are passed by reference. If they are not, they are automatically passed by value. Hence, there is another potential question mark here: in my cell right now, I am dealing only with NeuronGroup(1,...) but what happens if I deal with NeuronGroup(N,...) with N>1?

The docs on the user-function definition, towards the end, address this issue somehow cryptically. But I cannot figure it out well. Concomitantly, unless I use C+±specific data containers like vectors, my C++ function right now is not dealing well with arguments meant to be arrays: I should include some check inside to see if I am dealing with scalars passed by reference or actual arrays of size >1, and in that case run for cycles for along the whole array.

Finally, there is the problem of rand() inside my C++ code. Again, in principle, it will work fine if my arguments are scalars. Still, suppose I am getting an array. In that case, I should generate a scalar value passing the index from whatever element shows a transition in my cell, which is not so obvious how to be obtained from the internal Brian code. It does not seem there is a simple solution, or Brian-friendly solution here…

3

If you look at the documentation you will see that functions that are used in the code can never be void, because we do not support statements with side effects. The only statements in Brian code you can execute are assignments, therefore you need at least something like dummy = function(...). You cannot pass values by pointer/reference either. And finally, you always need check_units (or the equivalent arguments of the Function object) to tell Brian the units of your input and output arguments so that it can check that things are used correctly.

As a side note, there was a discussion about a similar model to yours here: Implementation of Markov channel models · Issue #1103 · brian-team/brian2 · GitHub
We’d need to support this kind of thing in Brian itself, it goes beyond our current infrastructure.

So at the current state, you have two options:

  1. Work on a more general support for loop structures in Brian (syntax, code generation, and all that) and earn eternal gratitude from developers and users alike.
  2. Ignore most of Brian’s code generation (doable it you only want it to work with C++ standalone), and do your own thing.

For 2, the easiest approach would be to have your function only get the “read-only” arguments as actual arguments + the index of the neuron, and then you access the global arrays storing the other variables yourself. You can either look up the names of these variables in the generated code or you use the Device.get_array_name function to get their name with something like device.get_array_name(cell.variables['tau']). You should also declare Brian’s rand function as a dependency and use it instead of the stdlib integer version which gives very bad random numbers.

So roughly things could look like this:

@implementation('cpp', '''
double nMGA(double dt, double Q, double p_0, int index) {
   while (True) {
     %tau_NAME%[index] = -log(_rand(index));
     ...
   }
}
'''.replace('%tau_NAME%', device.get_array_name(cell.variables['tau'])),
dependencies={'_rand': DEFAULT_FUNCTIONS['rand']})
@check_units(dt=second, Q=1, p0=1, index='integer', result=1)
def nmGA(dt, Q, p0, index):
   raise NotImplementedError('use standalone mode')

and used as

dummy = nMDG(dt, Q, p_0, i)

in the Brian code.

For debugging and such, you will have to look at the generated code most likely, but then it should be pretty straightforward. Our generated standalone code is not that complex.

As stated in the docs, in standalone mode functions always get a single value.

The variable i is the index of the neuron, so you should simply pass this as an argument.

4

@mstimberg Thanks. This is really an interesting workaround. Before I elaborate concerning “option 1,” – thanks for trusting my programming skills – but for now, I am still on the end-user side :slight_smile:

I am experimenting with your suggestion. However, two issues seem to emerge, plus a third one that actually falls in a side question.

First: device.get_array_name
Looking into the device.get_array_name docs, the method requires a global variable to retrieve the actual array name from. What if I am actually nesting my cell definition inside a method that returns cell as an argument of a class? What is the actual global variable? Global wherein? Python or CPP? And global concerning what scope? Global the whole script scope? Or global to a parent method? The following is the hierarchy I am dealing with:

def cell_compartment(N,pars):
    ...
    compartment = NeuronGroup(N,...,events={'transitioned': 'transition', 'no_transition': 'not transition'}))
    code_upon_transition = """
                           rval = rand()
                           p_O = q_DP/Q
                           X_O = X_O + (rval<=X_O*p_O)*1 - (rval>X_O*p_O)*1
                           fraction = 1 - tau/(Q*dt/second)
                           dummy = nMGA(dt,Q,p_O,i)  <----- Referring to this line also in the **Second point**
                           """
    code_upon_notransition = """
                             tau -= Q*dt/second
                             """
    compartment.run_on_event('transitioned',code_upon_transition)
    compartment.run_on_event('notransition',code_upon_notransition)
    return compartment

def compartment_geometry(N,...):
    comps = cell_compartment(N,...)
    ... 
    return Network([comps,...])

class ComplexCell(object):
    def __init__(self, N,...):
        # Pre-defs of standalone
        set_device('cpp_standalone', directory=code_dir, build_on_run=False)
        device.delete(force=True)  # Clean codegen directory
        self._code_dir = code_dir
        self.__build = False

        # Generate model
        self.celltree = ComplexCell(N,...)

    def simulate(self,duration=1*second,build=True):
        prefs.devices.cpp_standalone.openmp_threads = 10
        prefs.logging.file_log = False
        prefs.logging.delete_log_on_exit = True
        defaultclock.dt = 0.01*msecond
        ## Run simulation
        self.celltree.run(duration, report='text')

        if build:
            # Build at the end
            device.build(directory=self._code_dir,clean=clean)
            self.__build = True

Hence if my CPP code is now:

@implementation('cpp', '''
double nMGA(double dt,double Q,double p_O,int index) {
    while(true) {
        %tau_NAME%[index] = -log(_rand(index));
        if (%tau_NAME%[index]>=%fraction_NAME%[index]*Q*dt_){
            %tau_NAME%[index] -= %fraction_NAME%[index]*Q*dt_;
        break;
        }
        else {
            %rval_NAME%[index] = _rand(index);
            %X_0_NAME%[index] += int(%rval_NAME%[index]<=%X_0_NAME%[index]*p_O) - int(%rval_NAME%[index]>%X_0_NAME%[index]*p_O);
            %fraction_NAME%[index] -= %tau_NAME%[index]/(Q*dt_);
        }
    }
    return 1.0;
}
'''.replace('%tau_NAME%', device.get_array_name(<???>.variables['tau'])).replace('%fraction_NAME%',\
            device.get_array_name(<???>.variables['fraction'])).replace('%rval_NAME%',\
            device.get_array_name(<???>.variables['rval'])).replace('%X_0_NAME%',\
            device.get_array_name(<???>.variables['X_0'])),
dependencies={'_rand': DEFAULT_FUNCTIONS['rand']})
@check_units(dt=second, Q=1, p_O=1, index='integer', result=1)
def nMGA(dt, Q, p_O, index):
   raise NotImplementedError('use standalone mode')

if __name__=='__main__':
    astro = ComplexCell(10,...)
    astro.simulate(duration=5*second)

What is my cell global variable (replaced by <???>)?

Second: index error

I am not sure. But my understanding from your explanation is that ‘i’ is a reserved symbol in the Brian scope or within the code and always refers to the index of elements in the NeuronGroup.

Then when I attempt to compile the above, I get and error when the parser gets to the the line where I am supposed to call the nMGA function in the code_upon_transition. Here is the full error:

Traceback (most recent call last):
File "/home/maurizio/local/brian2/brian2/core/network.py", line 897, in before_run
obj.before_run(run_namespace)
File "/home/maurizio/local/brian2/brian2/groups/group.py", line 1143, in before_run
self.create_code_objects(run_namespace)
File "/home/maurizio/local/brian2/brian2/groups/group.py", line 1136, in create_code_objects
code_object = self.create_default_code_object(run_namespace)
File "/home/maurizio/local/brian2/brian2/groups/group.py", line 1118, in create_default_code_object
self.codeobj = create_runner_codeobj(group=self.group,
File "/home/maurizio/local/brian2/brian2/codegen/codeobject.py", line 396, in create_runner_codeobj
check_units_statements(c, variables)
File "/home/maurizio/local/brian2/brian2/equations/unitcheck.py", line 99, in check_units_statements
expr_unit = parse_expression_dimensions(expr, variables)
File "/home/maurizio/local/brian2/brian2/parsing/expressions.py", line 315, in parse_expression_dimensions
arg_idx = func._arg_names.index(expected_unit)
AttributeError: 'NoneType' object has no attribute 'index'

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
File "/home/maurizio/Ongoing.Projects/FSU.Project/code/astrocyte_model_class.py", line 574, in <module>
cell.simulate(duration=5*second)
File "/home/maurizio/Ongoing.Projects/FSU.Project/code/astrocyte_model_class.py", line 554, in simulate
self.celltree.run(duration, report='text')
File "/home/maurizio/local/brian2/brian2/core/base.py", line 276, in device_override_decorated_function
return getattr(curdev, name)(*args, **kwds)
File "/home/maurizio/local/brian2/brian2/devices/cpp_standalone/device.py", line 1375, in network_run
net.before_run(namespace)
File "/home/maurizio/local/brian2/brian2/core/base.py", line 278, in device_override_decorated_function
return func(*args, **kwds)
File "/home/maurizio/local/brian2/brian2/core/network.py", line 899, in before_run
raise BrianObjectException("An error occurred when preparing an object.", obj) from ex
brian2.core.base.BrianObjectException: Error encountered with object named "comp_resetter".
Object was created here (most recent call only, full details in debug log):
File "/home/maurizio/Ongoing.Projects/FSU.Project/code/astrocyte_model_class.py", line 406, in AstroCompartment
astro_compartment.run_on_event('transitioned',code_upon_transition)

So for some reason, it seems that the index is NoneType – could that be because right not my N in terms of comps instantiation is 1 – which would imply treatment of scalar NumPy arrays differently from actual arrays (i.e., not atleast_1d for scalar arrays).

Third: Replacing internal variables by unitless versions

Looking at your suggestion, I recall that Brian offers this option of retrieving the unitless variables by appending the underscore _ character to the variable name. However, if this allegedly works in the CPP code, when I replaced dt/second by dt_, for example, in my above code_upon_transition I get the unidentified identifier error. So it seems that Brian does not like the underscored version of its internal variables. Is that right? I am a bit confused.

5

Sometimes I wonder whether you have access to the same documentation that I have :slight_smile: Where does the documentation talk about a “global variable”? You have to hand the Variable object over to the get_array_name function, and this object is stored in the variables dictionary of the respective object, e.g. the NeuronGroup or Synapses . That’s why I used cell.variables in my example, since cell was the name of the NeuronGroup that you used in the code you provided. It does not matter how you organize your code in classes/methods, you’ll just need access to the NeuronGroup object (could also be via network['name_of_object']).

Yes, see the documentation on state variables or the list of special symbols.

You are reading this error message incorrectly, it is complaining about func._arg_names being None. It is trying to call Python’s index method on this list, i.e. this has nothing to do with the index from your code. I am not 100% sure why this happens, but this seems to be a bug in Brian (and potentially a regression that we introduced recently when we improved the handling of units in functions, in particular the error messages).
For now, you can use the alternative syntax that does not use decorators (also, sorry about the index='integer', this should not be in the unit definition and just be index=1). The alternative syntax works like this:

nMGA = Function(None,  # for "no Python implementation"
                arg_units=[second, 1, 1, 1],
                arg_names=['dt', 'Q', 'p_O', 'index'],
                return_unit=1)
nMGA.implementations.add_implementation('cpp', '''
double nMGA(double dt, double Q, double p_O, int index) {
   ...
}
'''.replace(...), dependencies={'_rand': DEFAULT_FUNCTIONS['rand']})

The unit system is only used in Python, outside of all generated code (so that it does not slow down computations). The C++ code does not know anything about units. The unit checks happen before any code is generated. A function that receives e.g. dt as an argument always receives a simple double, regardless of whether it gets passed dt or dt/second (the latter will be translated into dt/1.0 in the generated code).

6

Thanks @mstimberg. So the solution works. The only minor correction to your code is that in the replace instance you need to explicitly resolve Brian’s scope, i.e. by brain::. That is, for example:

.replace('%X_O_NAME%','brian::' + device.get_array_name(<brian_obejct>.variables['X_O']))

I can foresee a way to actually implement the non-Markov Gillespie algorithm along the lines of this post. And I hope to post an update of a sample implementation here in the next few days. The only restriction set by Brian is the time-dependent rates for the different transitions among species need to be the same for all molecules in the population. Regretfully, after some trial and error, I figure out that this is not my case, in as much as some transitions are dependent on the molecule state. In this fashion, I resorted to using a NeuronGroup for N molecules and model stochastic transitions in the “usual way” meaning by single molecules random drawing and rejections. It is not going to be exact as the GA, since the accuracy of procedure depends on dt, but for sufficiently small dt I can hopefully still get acceptable results.