GSL integration with CPP/standalone

[mdepitta]

Hi folks,
What's the implementation status of GSL integrators with CPP/standalone?
I have a tricky HH-type model that, if I integrate by Euler in Brian2, it shows firing. In contrast, in a separate implementation outside Brian2, it fires -- and further analysis proves that it is supposed to fire (periodically).

I tried to run my code using GSL integrators, but they don't work:

from brian2 import *
from brian2.devices.device import reinit_devices
from brian2.utils.logger import clean_up_logging
code_dir = './codegen'
prefs.GSL.directory = '/usr/include/'
set_device('cpp_standalone',directory=code_dir,build_on_run=False)
prefs.devices.cpp_standalone.openmp_threads = 4
prefs.logging.file_log = False
prefs.logging.delete_log_on_exit = True`
##
device.delete(force=True)
## Build the neuron model and monitors
params = lpc5_parameters(model='hh-neuron',T_ramp=duration*second,**kwargs)
cell = lpc5_neuron(1,params,name='HH',dt=1*us) #--> This return a NeuronGroup where internally, method='gsl'
cell.v = -70*mV
v_mon = StateMonitor(cell,variables='v',record=True,dt=0.1*ms,name='vmon')
network = Network([cell,v_mon])
## Run the simulator
network.run(duration=duration*second,report='text')
device.build(directory=code_dir, clean=True)
```

The above code generates:

```
Traceback (most recent call last):
  File "/home/xxx/Python3.10-VENV/lib/python3.10/site-packages/brian2/core/network.py", line 1003, in before_run
    obj.before_run(run_namespace)
  File "/home/xxx/Python3.10-VENV/lib/python3.10/site-packages/brian2/groups/group.py", line 1266, in before_run
    self.create_code_objects(run_namespace)
  File "/home/xxx/Python3.10-VENV/lib/python3.10/site-packages/brian2/groups/group.py", line 1259, in create_code_objects
    code_object = self.create_default_code_object(run_namespace)
  File "/home/xxx/Python3.10-VENV/lib/python3.10/site-packages/brian2/groups/group.py", line 1240, in create_default_code_object
    self.codeobj = create_runner_codeobj(
  File "/home/xxx/Python3.10-VENV/lib/python3.10/site-packages/brian2/codegen/codeobject.py", line 484, in create_runner_codeobj
    return device.code_object(
  File "/home/xxx/Python3.10-VENV/lib/python3.10/site-packages/brian2/devices/cpp_standalone/device.py", line 704, in code_object
    codeobj = super().code_object(
  File "/home/xxx/Python3.10-VENV/lib/python3.10/site-packages/brian2/devices/device.py", line 368, in code_object
    codeobj.compile()
  File "/home/xxx/Python3.10-VENV/lib/python3.10/site-packages/brian2/codegen/codeobject.py", line 125, in compile
    self.compiled_code[block] = self.compile_block(block)
  File "/home/xxx/Python3.10-VENV/lib/python3.10/site-packages/brian2/codegen/codeobject.py", line 121, in compile_block
    raise NotImplementedError("Implement compile_block method")
NotImplementedError: Implement compile_block method
```

Any workaround?
Thank you.

M

[mstimberg]

Hi Maurizio,
thank you for letting us know! It turned out, that changes in C++ standalone mode broke the GSL support, but due to some silliness in our test suite, these errors became skipped tests instead of failures, so we didn't notice them. Both the test suite and the actual errors should be fixed in PR #1687. I'll wait for the test suite to run through before merging.

All that said, are you sure you need the GSL integration for your model? That forward Euler does not work correctly for an HH model is not surprising, but did you try exponential_euler or rk_4?

[mdepitta]

Thanks, Marcel.
It definitely works either by 'rk4' (you have a typo in your reply) or 'exponential_euler' -- still, I wanted to use GSL for some gain in adaptive steps. On the same page: any chance you will implement any backward integration methods in Brian2? Having a simple backward Euler or the more powerful RADAU would really assure that we can integrate any problem -- in fact, eventually you would use Brian2's interface to incorporate anything quickly -- or maybe just port Brian2 to Julia :)

[mstimberg]

No plans to implement backward Euler, and also we'd have to decide whether we'd implement it by solving the resulting equation system or by using some iterative method (or both) – in either case a non-trivial change. I'd like to first see a real advantage over the existing methods in a relevant neuron model before undertaking such an effort. Note that multicompartmental neurons already use a backward Euler integration for the propagation of currents between the compartments.