# Library of functions for Brian equations

**URL:** <https://brian.discourse.group/t/library-of-functions-for-brian-equations/333>\
**Category:** Support\
**Created:** [9 March 2021 15:29 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333 "2021-03-09T15:29:53Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [9 March 2021 15:29 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/1 "2021-03-09T15:29:53Z")

</div>

Hi everyone,  
This isn’t a classic support request, sorry, but I couldn’t find a good category for my question.

With years of Brian/Brian2 usage, I accumulated a bunch of functions that I really would like to have inside equations. Moreover, I want them on all devices `cpp`, `numpy`, `cython`, etc. So is there any template for such a library?

What will be very handy is to have something like this:

```python
from brian2 import *
from brian2pls import fe

n = NeuronGroup(..., model ="dr/dt = fe(r)/tau : 1" ....)

```

It seems decorator `@implementation` is the best way to go. So a `fe` function may look like this:

```python
@implementation('cython', f'''
     cdef double fe(double x):
         return {ge}*(x-{thetae}) if x-{thetae} > 0. else 0.
     ''')
@implementation('cpp',f''''
     double fe(double x)
     {{
         return ((x-{thetae})<0.)?0.:{ge}*(x-{thetae});
     }}''')
@check_units(x=1, result=1)
def fe(x):
     return ge*(x-thetae) if x-thetae > 0. else 0.

```

This seems (I’m really not sure) works for both `cpp`, standalone, omp device, and `cython` device. It isn’t completely clear to me how to implements it for `numpy` device, and do we have `cuda`/`gpu` device or not?

Any thoughts, suggestions, and examples are highly appreciated.

P.S. As we discussed with @mstimberg before, many of my functions can be written as a combination of brian2’s equation functions, but it is too messy. On the other hand, I want to have full control over `cpp`/`cython` implementation, to be sure that it is the most optimal one.

---

<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 March 2021 16:15 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/2 "2021-03-09T16:15:13Z")

</div>

> [@rth](#):
>
> This isn’t a classic support request, sorry, but I couldn’t find a good category for my question.

I personally don’t mind having it here in this category, otherwise #feature-requests might be a good fit?

> [@rth](#):
>
> With years of Brian/Brian2 usage, I accumulated a bunch of functions that I really would like to have inside equations. Moreover, I want them on all devices `cpp` , `numpy` , `cython` , etc. So is there any template for such a library?

With your specific example, I’d probably rather write it in the equations as `clip(ge*(x - thetae), 0, inf)`, and then ideally we’d have “macros” in the equations ([Macro definitions in equations · Issue #378 · brian-team/brian2 · GitHub](https://github.com/brian-team/brian2/issues/378)) so you could define it once and then refer to it under some name in multiple places. But I get the general idea, and we made the function system extensible to support use cases where you want to have detailed control about the implementation. If you want to generate functions in a generic way, then it would probably be more straightforward to not use the decorators but to create a `Function` object yourself. E.g. your example above would look like this:

```Python
def make_fe(ge, thetae):
    # Python implementation
    def fe(x):
        # operates on an array of values
        return np.where(x - thetae > 0, ge * (x - thetae), 0.)

    # Cython implementation
    cython_code = f'''
     cdef double fe(double x):
         return {ge}*(x-{thetae}) if x-{thetae} > 0. else 0.
     '''

    # C++ implementation
    cpp_code = f'''
     double fe(double x)
     {{
         return ((x-{thetae})<0.)?0.:{ge}*(x-{thetae});
     }}'''
    func = Function(fe, arg_units=[1], arg_names=['x'], return_unit=1)
    func.implementations.add_implementation('cython', cython_code)
    func.implementations.add_implementation('cpp', cpp_code)
    return func

fe = make_fe(5., 0.1)

```

> [@rth](#):
>
> This seems (I’m really not sure) works for both `cpp` , standalone, omp device, and `cython` device. It isn’t completely clear to me how to implements it for `numpy` device, and do we have `cuda` / `gpu` device or not?

Yes, your example will work for C++ standalone with or without OpenMP and for Cython. For numpy you have to be aware that functions operate on arrays and not on single values as the other targets (see [Functions — Brian 2 2.5.4 documentation](https://brian2.readthedocs.io/en/stable/advanced/functions.html#arrays-vs-scalar-values-in-user-provided-functions)).

> [@rth](#):
>
> do we have `cuda` / `gpu` device or not?

The `brian2cuda` project is still in alpha, but in principle could be used here I guess. You can also run code on the GPU via the GeNN simulator with `brian2genn` and you can define functions for it. You only have to be aware that you prefix your function header with `SUPPORT_CODE_FUNC` ([GeNN: Neuron models](http://genn-team.github.io/genn/documentation/4/html/de/ded/sectNeuronModels.html#neuron_support_code)), so that GeNN can generate both code for CPU and GPU. See [https://github.com/brian-team/brian2genn/blob/master/brian2genn/genn\_generator.py](https://github.com/brian-team/brian2genn/blob/master/brian2genn/genn_generator.py) for examples of function definitions in Brian2GeNN.

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [9 March 2021 16:37 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/3 "2021-03-09T16:37:30Z")

</div>

@mstimberg thank you for explanations, very useful.

A side question, I have 99% of my functions implemented as C/C++ macros in a header file. Very useful, when you code in C/C++/Cython. How can I add `#include` statement in `cpp` implementation? Well now I realized this isn’t a great idea at all. I will have to include h file into the module 🙄

---

<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 March 2021 16:50 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/4 "2021-03-09T16:50:36Z")

</div>

You can ask Brian to include additional header files: [Functions — Brian 2 2.5.4 documentation](https://brian2.readthedocs.io/en/stable/advanced/functions.html#additional-compiler-arguments) If this header file implements a function as a macro, you might not have to implement the function at all, you can simply state that the implementation exists (with a different name potentially).

> [@rth](#):
>
> I will have to include h file into the module

Not quite sure that I understand. You mean that you’d have to package it with the Python module?

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [9 March 2021 17:00 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/5 "2021-03-09T17:00:57Z")

</div>

> Blockquote Not quite sure that I understand. You mean that you’d have to package it with the Python module?

Yes, but I thought to have a simple, ‘single file’ module, not a _whole python’s module package_.

BTW, Can one Function() object calls another Function() object inside? Will it work for all devices?

```python
fe = Function(....)
fi = Function(which_uses_fe)

```

It seems, it should work, but I’ve tested that yet.

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [9 March 2021 17:07 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/6 "2021-03-09T17:07:04Z")

</div>

Ah! @mstimberg I found everything in documentation. Thank you!

---

<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 March 2021 17:08 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/7 "2021-03-09T17:08:06Z")

</div>

> [@rth](#):
>
> BTW, Can one Function() object calls another Function() object inside? Will it work for all devices?

You will have to declare it in the implementation’s `dependencies`, see e.g. how the `poisson` function depends on the `rand` function: [https://github.com/brian-team/brian2/blob/master/brian2/codegen/generators/cpp\_generator.py](https://github.com/brian-team/brian2/blob/master/brian2/codegen/generators/cpp_generator.py)

---

<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 March 2021 17:08 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/8 "2021-03-09T17:08:30Z")

</div>

> [@rth](#):
>
> Ah! @mstimberg I found everything in documentation. Thank you!

Ah, saw this too late 😀

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [9 March 2021 18:39 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/9 "2021-03-09T18:39:36Z")

</div>

> [@mstimberg](#):
>
> You will have to declare it in the implementation’s `dependencies` , see e.g. how the `poisson` function depends on the `rand` function: [brian2/cpp\_generator.py at master · brian-team/brian2 · GitHub](https://github.com/brian-team/brian2/blob/master/brian2/codegen/generators/cpp_generator.py)

Well now I hit a wall 😕

Here is a simple example, and none of `cython` or `cpp` works.

```python
# brian2pls.py
from brian2 import *

def _P1(X,X0):
    return X0-X
P1 = Function(_P1,arg_units=[1,1], return_unit=1)
P1.implementations.add_implementation('cpp','#define P1(X,X0) ((X0)-(X))')
P1.implementations.add_implementation('cython','''
cdef double P1(double X, double X0):
    return (X0-X)
''')
def _P2(X,X0,X1):
    return P1(X,X0)*P1(X,X1)
P2 = Function(_P2,arg_units=[1,1,1], return_unit=1)
P2.implementations.add_implementation('cpp','#define P2(X,X0,X1) (P1((X),(X0))*P1((X),(X1)))',dependencies={'P1':P1})
P2.implementations.add_implementation('cython','''
cdef double P2(double X, double X0, double X1):
    return P1(X,X0)*P1(X,X1)
''',dependencies={'P1':P1})

```

and usage"

```python
from brian2 import *
from brian2pls import *

set_device('cpp_standalone', directory='CODE')
prefs.devices.cpp_standalone.openmp_threads = 1

equ="""
dv/dt = (P2(v,-50,-30) + I )/ms : 1
I : 1
"""

n = NeuronGroup(1,equ,threshold="v>0",reset="v=-80",refractory=1*ms)
r = StateMonitor(n,'v',record=True)
n.v=-70.
n.I=0

run(200*ms)

plot(r.t/ms,r.v[0])
show()

```

Cython generator (default) returns error:

```auto
INFO No numerical integration method specified for group 'neurongroup', using method 'euler' (took 0.01s, trying other methods took 0.01s). [brian2.stateupdaters.base.method_choice]

Error compiling Cython file:
------------------------------------------------------------
...

# support code

cdef double P2(double X, double X0, double X1):
    return _P1(X,X0)*_P1(X,X1)
          ^
------------------------------------------------------------

```

C++ generator returns

```auto
code_objects/neurongroup_stateupdater_codeobject.cpp: In function ‘void _run_neurongroup_stateupdater_codeobject()’:
code_objects/neurongroup_stateupdater_codeobject.cpp:15:29: error: ‘_P1’ was not declared in this scope; did you mean ‘P1’?
   15 | #define P2(X,X0,X1) (_P1((X),(X0))*_P1((X),(X1)))
      | ^~~
code_objects/neurongroup_stateupdater_codeobject.cpp:125:42: note: in expansion of macro ‘P2’
  125 | const double _v = (_lio_2 * (I + P2(v, - 50, - 30))) + v;
      | ^~
make: *** [makefile:31: code_objects/neurongroup_stateupdater_codeobject.o] Error 1

```

What did I mess up?

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [9 March 2021 19:33 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/10 "2021-03-09T19:33:02Z")

</div>

but decorators work pretty well!

```python
@implementation('cpp', '''
#define P1(X,X0) ((X0)-(X))
''')
@implementation('cython', '''
cdef double P1(double X, double X0):
    return (X0-X)
''')
@check_units(arg=[1,1], result=1)
def P1(X,X0):
    return X0-X

@implementation('cpp', '''
#define P2(X,X0,X1) (P1((X),(X0))*P1((X),(X1)))
''',dependencies={'P1':P1})
@implementation('cython', '''
cdef double P2(double X, double X0, double X1):
    return P1(X,X0)*P1(X,X1)
''',dependencies={'P1':P1})
@check_units(arg=[1,1,1], result=1)
def P2(X,X0,X1):
    return P1(X,X0)*P1(X,X1)

```

```python
from brian2 import *
from brian2pls import *

set_device('cpp_standalone', directory='CODE')
prefs.devices.cpp_standalone.openmp_threads = 1

equ="""
dv/dt = (P2(v,-50,-30)*1e-5 + I )/ms : 1
I : 1
"""

n = NeuronGroup(1,equ,threshold="v>0",reset="v=-80",refractory=1*ms)
r = StateMonitor(n,'v',record=True)
n.v=-70.
n.I=0

run(200*ms)

plot(r.t/ms,r.v[0])
show()

```

---

<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 March 2021 09:34 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/11 "2021-03-10T09:34:04Z")

</div>

In your first code, you named the Python function `_P1` with an underscore, and Brian takes this as the name of the function implementations (it does not parse the Cython/C++ code to figure out the actual names, there). Then, when you use `{'P1': P1}` as a dependency, it will understand that `P1` refers to the function implementation `_P1`, and renames the function calls in your code (in the error messages, it therefore complains about not finding `_P1`).

To fix this, you have two options. The first option is to rename your Python functions so that the name matches the name of the functions in the Cython/C++ code. When you then create a `Function` object of the same name this would overwrite the Python function name, but that’s not a problem. This is actually what happens in your second example where you annotate the `P1` function with decorators – `P1` is now a `Function` object. The second option is to add `name='P1'` and `name='P2'` to the `add_implementation` calls for Cython/C++ to tell Brian the name of the functions in the code.

But of course, if the decorator approach is sufficient for you and you prefer it, no need to go the `Function.add_implementation` route.

---

<div class="post-metadata">

**Author:** ![rth](https://yyz2.discourse-cdn.com/free1/user_avatar/brian.discourse.group/rth/32/16_2.png) [@rth](https://brian.discourse.group/u/rth)\
**Post date:** [10 March 2021 13:27 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/12 "2021-03-10T13:27:12Z")

</div>

Thank you, @mstimberg! Yes I realized letter that difference in python function name and in Brian’s function name creates a problem.

I went with decorators, it is much simpler, and I still have the python function with the same name. That is really handy for phase-plan analysis, because it is so simple to reconstruct nullclines in the python code and still have a very fast simulation with the same functions in Brian. I’m not sure that I can call Brian’s function as a regular python function outside of equations statement. Correct me if I’m wrong.

Thank you for your help!

---

<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 March 2021 13:34 UTC](https://brian.discourse.group/t/library-of-functions-for-brian-equations/333/13 "2021-03-10T13:34:10Z")

</div>

> [@rth](#):
>
> I’m not sure that I can call Brian’s function as a regular python function outside of equations statement. Correct me if I’m wrong.

If they have have a Python/numpy implementation, you can, but not if there is only a Cython/C++ implementation. We could actually write a simple Python wrapper in such cases, might be useful for testing etc.  
In general, the decorator and the `Function` / `add_implementation` approach are mostly the same under the hood (e.g. the `implementation` decorator simply calls `add_implementation` on the `Function` object).