I have a very complicated model with over 60 equations. I’m fitting some parameters of this model with Brian2modelfitting. ZeroDivisionError can be encountered for some samples with unsuitable parameters. For example, the design space is \{(a,b)|1\leq a\leq 2,1\leq b\leq 2\}. (a=1 and b=1), (a=2 and b=2), or (a=1, b=2) will yield reasonable results, but (a=2, b=1) will not converge and evoke ZeroDivision. Ideally, the optimizer will learn that (a=2, b=1) is not a good set of parameters due to divergence, and try other parameters far from (a=2, b=1). Unfortunately, the ZeroDivisionError will halt the optimization.
One plausible remedy may be constraining the parameters to be optimized within proper bounds. But it is not feasible in my case because (a=2, b=1) is within the respective bounds.
For brian2 only, if I encounter this kind of error, I may monitor some variables to debug. But for model fitting, I cannot even know what parameters lead to such an error, adding to the difficulty in debugging.
I wonder how to catch the ZeroDivisionError and change it into a very large error, so that the optimization can continue while keeping away from the bad parameter choices.
Users can customize the Metric to bypass the unexpected behavior owing to NaN value. When it comes to ZeroDivisionError, however, I cannot come up with any solutions.
Minimal code for reproduction
I think this is a very general problem and only need some new thoughts and advice (no code is OK). The code is very complicated that it is hard for me to simplify it as minimal code which can yield ZeroDivisionError. If you really need the code, see chambers-fit.py (6.3 KB)
This seems to be a numerical issue – when I used a smaller timestep, i.e., 0.001 * ms, the ZeroDivisionError could be avoided, which can be computationally expensive. I want to find out which equation caused the error so that I can continue to use 0.01 * ms as the timestep. But this is not possible for the current brian2modelfitting simply because I could not get the parameter set during optimization to reproduce the error.
Hi @Jasmine9691. I wonder whether it would be good to have an (optional) pre-simulation stage where one could reject/accept parameter combinations. Currently, we can only give a range for the parameters, but sometimes you might need to avoid certain parameters or parameter combination, and we don’t have a way to express that.
Regarding your specific problem: the ZeroDivisionError is actually only an issue in Cython (the default runtime backend), since it emulates Python’s handling of division by zero. If you run things with the numpy target (prefs.codegen.target = 'numpy') or with the C++ standalone mode (set_device("cpp_standalone")), things should actually work (with numpy you’d get a warning about the division by zero but not an error). If you can, I’d definitely recommend using the C++ standalone mode, since it should be quite a bit faster than the Cython runtime.
Here you can see the difference between numpy and Python (and Cython, which emulates this behaviour):
>>> import numpy as np
>>> 1/np.arange(3)
<stdin>:1: RuntimeWarning: divide by zero encountered in divide
array([inf, 1. , 0.5])
>>> 1/0
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
ZeroDivisionError: division by zero
On a more general note: if you want to debug equation issues, it is often the best approach to use the numpy target. In your case, it only emits a warning, but you can have it raise an error with np.seterr(divide='raise'). This will give you at least the actual line that raised the error, instead of something not useful like compiled_code.main(). You also have access to the variable values if you run it in a debugger, but it is not super-straightforward since you are dealing with a “Code object” that uses exec with Python code, and not a straightforward Python file.
Hope that helps!
PS: I don’t actually remember in detail why we use Python and not C division semantics in Cython, but I think it had to do with the fact that is also linked to modulo semantics and not-obvious corner cases (-1 % 3 equals 2 in Python but -1 in C…). We have a very old issue open about (optionally) enabling C semantics: Add option for cdivision:True in Cython · Issue #993 · brian-team/brian2 · GitHub
I tried the numpy codegen, it is much slower than cython (30 min vs 30 sec for a 10-second simulation) so I have to wait for a long time before the error is raised.
Anyway, my problem has been resolved before, and you taught me a new way for debugging. I may try it next time!
This model is a part of a more comprehensive model which has to be communicated with another Python library to get the stimulus (see the link above). It is not trivial to do the coupling in the cpp-standalone mode because many Python features are not supported. I am not an expert in cpp, I do not have enough time and ability to fix the potential bugs. So I prefer a time-consuming but lazy way, that is cython.
Also, note that the most time-consuming part is the repeated restart of computation due to coupling, instead of choice of cython/cpp.
Ok, I see, thanks for the explanation. If you don’t mind editing the Brian source manually, you could change the cdividion=False line in brian2/codegen/runtime/cython_rt/templates/common_group.pyx to cdivision=True which should make Cython handle divisions by 0 in the same way as the numpy target.