Implement "exponential euler" integration for conditionally linear equations

[mstimberg]

Integrate with the system described in #11

[mstimberg]

I started implementing this state updater in the exponential_euler branch. Similar to the current exact integrator for linear equations it is using sympy and generates abstract code for calculating the coefficients. This needs more testing, though -- it runs fine and gives "reasonably-looking" results but there are significant differences to Brian1's results so something is not quite right. Oh and we should consider including something like the vtrap function mentioned by Margarita on the mailing list a while ago: The COBAHH example contains expressions with (exp(...) - 1) in the divisor, leading to inf if the argument of exp is 0 (and eventually this leads to NaN when the value is multiplied with 0, for example).

[thesamovar]

That's worrying - I hope the Brian 1 implementation is right!

[mstimberg]

For very simple examples, both give the same results. For the complete COBAHH model the results differ a bit (and small differences can lead to a spike in one or the other) -- I'll still have to downsize the model to something simpler to get a better idea of where the difference lies. The differences are not big enough for real fundamental differences so it's possible that it's all about rounding errors that are progagated differently: In Brian 1, the A and B values are precalculated and extracted from the equations using the x=1 or x=0 trick, whereas in my current implementation for Brian2, everything is done symbollically. I'll try to test a stiff equation with a known analytical solution, next. Or do you have any other idea how to test this?

[thesamovar]

Just to recap what we said at the meeting: if it's an accumulation of small numerical differences, then we should see it by e.g. driving a HH model with a subthreshold frozen noise current - the differences should gradually increase starting from 0, not suddenly diverge right from the beginning.

[mstimberg]

Ok, I'm getting closer. A big chunk of the difference is actually not coming from rounding errors but from a more fundamental sympy issue. The HH equations contained expressions like 40*mV and in the current code, I'm not replacing mV by 0.001 but leaving it as it is (the namespace will map it to the value on execution). This is not a problem in itself, but sympy does some simplication of the equations and generates expressions like 13/4 -- which unfortunately are interpreted as integer division in Python and C... For now, the problem is circumenvented by only using floating point values in the equations (e.g. 40. * mV) but this is a more fundamental issue we should cope with in a general way. Still, some numerical differences remain, I'll check them more thoroughly now.

[thesamovar]

Just to recap our meeting: the resolution to this problem should probably be to use sympy rewriting on Equations and turn all integers into floats, and not to use rewriting on other abstract code because the user may be expecting integer arithmetic. We could go all Python 3 and insist that integer division should use the // operator.

[mstimberg]

Arg, had a major head slap moment here... After spending hours to find the remaining numerical difference between Brian 1 and Brian 2, trying to replace all constants by literal values, all literal values by constants, changing sympy simplifications, etc. I realized that the issue has nothing to do with the exponential Euler state updater at all, but can be also seen with Euler or RK integration... Integration timesteps were almost always identical, except for some rare cases (and then trajectories started to diverge, of course). After way too much time it turned out that this had nothing to do with the integration at all, but only with my ad-hoc implementation of a TimedArray for Brian 2 (which should be another point for our syntax/feature discussion, by the way) which used the dangerous i = int(t / clock.dt) expression which almost always does the correct thing... Even more embarrassing is that Victor mentioned this problem today on the mailing list and I did not think a single moment that it could apply to my problem...

But the good news: The integration is consistent with Brian 1, I tested the HH model and for a subthreshold current with frozen noise, differences are on the order of 10e-13 mV (I did not check but I think this is close to machine precision), this sounds very reasonable (same for Euler and RK integration). If I let the neuron spike, numerical differences directly at the action potential go up to ~10e-10 mV, perfectly reasonable I think (Euler and RK were not stable with the same dt). The differences are also positive and negative, i.e. there is no constant bias. In most simulations, I don't think one should see any difference. I'll clean up things and then create a pull request. Let's deal with the integer issue separately in issue #51 .

[thesamovar]

Good catch! I've run into this problem more than once I have to admit.

For the implementation/syntax of TimedArray, we should indeed probably discuss it and make an issue for it but I'm more or less happy with what we have and it's probably not a top priority.