# User-defined functions

**URL:** <https://brian.discourse.group/t/user-defined-functions/271>\
**Category:** Support\
**Created:** [29 January 2021 18:06 UTC](https://brian.discourse.group/t/user-defined-functions/271 "2021-01-29T18:06:18Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![mdepitta](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mdepitta/32/21_2.png) [@mdepitta](https://brian.discourse.group/u/mdepitta)\
**Post date:** [29 January 2021 18:06 UTC](https://brian.discourse.group/t/user-defined-functions/271/1 "2021-01-29T18:06:18Z")

</div>

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
```

---

<div class="post-metadata">

**Author:** ![mdepitta](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mdepitta/32/21_2.png) [@mdepitta](https://brian.discourse.group/u/mdepitta)\
**Post date:** [30 January 2021 07:37 UTC](https://brian.discourse.group/t/user-defined-functions/271/2 "2021-01-30T07:37:23Z")

</div>

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…

---

<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:** [1 February 2021 11:11 UTC](https://brian.discourse.group/t/user-defined-functions/271/3 "2021-02-01T11:11:46Z")

</div>

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](https://github.com/brian-team/brian2/issues/1103)  
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:

```Python
@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

```auto
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.

> [@mdepitta](#):
>
> 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` ?

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

> [@mdepitta](#):
>
> 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.

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

---

<div class="post-metadata">

**Author:** ![mdepitta](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mdepitta/32/21_2.png) [@mdepitta](https://brian.discourse.group/u/mdepitta)\
**Post date:** [2 February 2021 09:47 UTC](https://brian.discourse.group/t/user-defined-functions/271/4 "2021-02-02T09:47:04Z")

</div>

@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 🙂

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.

---

<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:** [2 February 2021 14:57 UTC](https://brian.discourse.group/t/user-defined-functions/271/5 "2021-02-02T14:57:02Z")

</div>

> [@mdepitta](#):
>
> 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?

Sometimes I wonder whether you have access to the same documentation that I have 🙂 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']`).

> [@mdepitta](#):
>
> 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` .

Yes, see the documentation on [state variables](https://brian2.readthedocs.io/en/stable/user/models.html#state-variables) or the [list of special symbols](https://brian2.readthedocs.io/en/stable/user/equations.html#list-of-special-symbols).

> [@mdepitta](#):
>
> So for some reason, it seems that the index is `NoneType`

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:

```Python
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']})

```

> [@mdepitta](#):
>
> 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.

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).

---

<div class="post-metadata">

**Author:** ![mdepitta](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/mdepitta/32/21_2.png) [@mdepitta](https://brian.discourse.group/u/mdepitta)\
**Post date:** [9 February 2021 08:55 UTC](https://brian.discourse.group/t/user-defined-functions/271/6 "2021-02-09T08:55:31Z")

</div>

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.