Simplify exact integration

[mstimberg]

The "exact" integration algorithm currently has some limitations and is not able to integrate some linear equations with constant coefficients. For example, it fails to integrate this (Tsodyks, Pawelzik, Markram, Neural Computation, 1998):

dx/dt = z / tau_rec          : 1
dy/dt = -y / tau_in          : 1
dz/dt = y/tau_in - z/tau_rec : 1

This is because it hardcodes the approach for non-homogeneous equations, and homogenous ones like the one above need to be handled slightly differently. Even simpler, equations like dv/dt = rate : 1 are not supported, which can be especially annoying if you add such an equation to an existing set of equations that can be integrated with exact. I first tried to tackle this by making exceptions for certain types of equations, but I finally found it simpler to convert all non-homogeneous equations into homogenous equations by transforming an equation like

dv/dt = -v/tau + inp : 1

into something along the lines of

dv/dt = -v/tau + v_inp : 1
dv_inp/dt = 0/second : 1

(which happens to be the approach used in Rotter & Diesmann (1999) for piecewise-constant input).

This simplifies the code quite a bit and also reduces our dependency on sympy details (we only call the matrix exponential, but there's no need any more to call the solve_linear_system method), which sometimes break with updates.
There's one disadvantage: the process of generating the update code takes ~30% longer, e.g. for the CUBA model 2.7s instead of 1.9s. There's no change in execution speed and the results are of course unchanged. I think the change is worth it despite the inconvenience of slightly longer code generation, and there are also ways to make this faster in the future. In particular, the equations are complicated by the int(not_refractory) term, and we might change this in the future to separate code generation for the True/False cases.

[thesamovar]

These sorts of changes seem important to get right. In principle, this all sounds fine and the code looks good, but I'd like to go through the maths before we commit it if that's ok? Could you explain the reasoning?

[thesamovar]

Also, what does the generated code look like for a couple of equations before and after? Just trying to get a feel for this change.

[mstimberg]

So the maths for the new version is quite easy to explain: we use the fact that the solution to dx/dt = Ax for x(t_0)=x_0 is x_0·e^A(t-t_0). In our case, t-t_0 is usually dt, except for the event-driven integration in synapses where it is the time since the last update. Now, this assumes a homogenous system, but we can easily transform a non-homogeneous one into a homogenous one by introducing a new variable. E.g. instead of dv/dt = -v/tau + c and v(t_0) = v_0, we reformulate things as dv/dt = -v/tau + v_c with dv_c/dt = 0 and v_c(t_0) = c.

The currently implemented solution is slightly different (and basically the same as the one used in Brian 1) and solves a linear system to get the matrix B to formulate things as dx/dt = (M - B)x. For some trivial cases this approach fails and in general it has two places (matrix exponential and solve_linear_system) where sympy can fail, which IMO makes things a bit more complicated.

I am off to my holidays unfortunately (maybe not the right word 😉 ), so more details on this have to wait until I'm back!

[thesamovar]

It would be good to think carefully about this as there are some hiccups that do occasionally crop up. One of the most common being that if you have two time constants for membrane and synapse, and those time constants have different symbols but the same value, you get nans (because there is a divide by the difference of the two time constants in the solution). I guess this corresponds to a singular matrix. That can be dealt with, but it's not straightforward. What if the time constants are different for every neuron? What if they change each time step? etc.