Implement exact integration for linear equations about brian2 HOT 15 CLOSED

A first implementation of this feature is now in the branch https://github.com/brian-team/brian2/tree/linear_equations.

It uses sympy to solve the equations analytically, which only works for systems that can represented by a diagonalizable matrices (because sympy's matrix exponential is only implemented for this case). I think the current implementation has one advantage over the implementation in Brian1, it allows for a varying constant part of the equation. For example, the phase locking example (https://github.com/brian-team/brian2/blob/linear_equations/examples/phase_locking.py), uses the equations:

dv/dt=(-v+a*sin(2*pi*freq*t)+b)/tau : 1
a : 1

In Brian1, we don't use linear integration here as the equation is both time-dependent and includes a variable parameter. I think it's ok to use it, though, as the coefficient of v is constant (it's just 1/tau).

The current approach I'm using is to allow the state updater to fill in constants, instead of passing them to the code generation stage. I guess this is more efficient, I think gcc for example is only able to substitute an expression like exp(-dt/tau) by exp(-0.1) but not by 0.99 ?

With the parameters in the above mentioned phase locking example, the linear state updater now generates the following abstract code:

v = 0.00498752080731768*a*sin(62.8318530717959*t) + 0.995012479192682*v + 0.00598502496878117

There are a couple of technical difficulties/decisions to take now:

• How to handle constant parameters? I implemented that it is ok for constant parameters to be contained as coefficients. Now, the easiest way to make this work is to regenerate state updaters on every run -- I think for the typical use cases this is enough and probably much simpler than having some elaborate system to invalidate a state updater?
• How to deal with namespace changes? A minor problem is that the state updater is already filling in values, so we could get rid of some variables in the code object's namespace -- but having them there does not really hurt. More problematic is that state updaters can also add something to the namespace. In the general case, the linear state updater will introduce a use of the exponential function. I think a good solution would be to allow state updaters to not only return code, but also additions to the namespace (I think we discussed something like this a while ago). This would also make the addition of rand by the stochastic state updaters cleaner.

A major drawback of the implemented approach is the restriction to diagonalizable matrices, e.g. the following simple system can't be integrated with it:

dv/dt = (ge+gi-(v+49*mV))/(10*ms) : volt
dge/dt = -ge/(5*ms) : volt
dgi/dt = -gi/(10*ms) : volt

Probably it's the best to have to kinds of linear state updaters in the end, the sympy-based updater allows to include terms that change in time, a numpy-based one (similar to what Brian1 has) would allow for more complex equations but requires that they do not change in time.

from brian2.

Just in case you hadn't seen them already, I wrote a few notes about this issue some time ago. I don't know if they're still relevant.

from brian2.

Interesting, this suggests that rather than dividing into linear/nonlinear we should divide into analytically solvable/not solvable equations, and use the analytical solution when it's available or fallback to a numerical method. This would cover slightly more than the current linear stuff, and we could even special case a few things that sympy currently can't do maybe?

The key thing is to make sure it works straightforwardly from the user POV. To that end, I think the user shouldn't have to know if it using a sympy-based updater or numpy-based one, etc. Ideally, they should specify just whether a variable is constant or not - if it is declared as constant we can do compile time simplifications, otherwise we don't. We still allow changing constants, but when we do that it means regenerating the state updater. I think this has an important use case, namely you might want to do multiple longish runs just changing one constant each time without wanting to rebuild the whole network (e.g. change tau only). I think several people do this for optimisation problems.

There are also efficiency issues as you mentioned. This page suggests that gcc can simplify exp(0.1), but probably we should check that - exp is expensive so it should show up easily. And yeah, I think unfortunately we also need the numpy-based state updater for cases where there is no access to gcc.

And we should definitely incorporate Cyrille's work here too.

On to the technical questions:

• How difficult is the invalidation idea? I think it could be implemented just by having for each state updater a list of the constants which it depends on being constant, and then when you change a constant value it just does a quick look up in this table and recomputes it if necessary. Alternatively, your idea of keeping it fixed per run is quite neat too, but in that case we should lock constant values during a run and raise an error if the user tries to change them. In fact, maybe your solution is better because it doesn't let the user make the mistake of declaring something as constant but then frequently changing it. I think both are nice, maybe we should get an input from @romainbrette?
• I agree on allowing state updaters to modify namespaces. Actually, I thought we had this already but I guess I misunderstood. :)

from brian2.

Just in case you hadn't seen them already, I wrote a few notes about this issue some time ago. I don't know if they're still relevant.

@rossant: Thanks. I think I saw it quite a while ago but I did not remember it. How do these notes fit in with what brian 1 does -- would this be a replacement for the "inexact exact" numerical integration? I have to admit I do not yet entirely understand what I implemented (but it seems to work 😄 ) -- I basically copied the approach in brian 1 and simply replaced numerical numpy/scipy functions with sympy's symbolical ones.

from brian2.

@mstimberg Yes. There's an analytical expression for the update matrix as a function of the equations matrix, but it involves an integral with a matrix exponential over [0, dt]. This integral can be computed numerically with an integration method. This is not that different from what is done in Brian 1 where the ODE is integrated with an Euler integration method and a fixed number of time steps over [0, dt] (if I remember correctly). But with this method, we can control the precision and compute the appropriate number of steps as a function of the model parameters.

from brian2.

Interesting, this suggests that rather than dividing into linear/nonlinear we should divide into analytically solvable/not solvable equations, and use the analytical solution when it's available or fallback to a numerical method. This would cover slightly more than the current linear stuff, and we could even special case a few things that sympy currently can't do maybe?

Sympy actually has quite powerful solvers for ODEs -- but only for single equations. I'm not sure how often we have this use case in Brian. Oh, BTW: There's ongoing work (since quite a while) to have a dedicated system of (linear) ODEs solver in sympy, but it will probably not be ready in the near future: sympy/sympy#1322
Generally speaking, there are a lot of opportunities for fancy stuff here. I don't think we need a lot of support for special cases at the moment, but it's good to have the current implementation because it points to the problems we might encounter. And we should make sure that state updaters with support for special cases will be a great and easy way to contribute to Brian :)

The key thing is to make sure it works straightforwardly from the user POV. To that end, I think the user shouldn't have to know if it using a sympy-based updater or numpy-based one, etc. Ideally, they should specify just whether a variable is constant or not - if it is declared as constant we can do compile time simplifications, otherwise we don't.

Sure, but that's very much in the line with how the new state updater mechanism works. Every state updater says whether it can integrate the equations, so the sympy updater could reject equations where the matrix is not diagonalizable and the numpy updater would reject ones where non-constant terms appear.

There are also efficiency issues as you mentioned. This page suggests that gcc can simplify exp(0.1), but probably we should check that - exp is expensive so it should show up easily. And yeah, I think unfortunately we also need the numpy-based state updater for cases where there is no access to gcc.

You mean from the performance point of view, right? To have the "withdot" performance instead of the "withpython" performance in your little benchmark?

About invalidation vs. fixed-per-run:
The invalidation principle has the advantage that it could be used more broadly, e.g. it could be combined with the refractory mechanism and we would recalculate/store the matrices when a neuron changes from active to refractory. The disadvantage is that it introduces some complexity to the state updater system in NeuronGroup. Currently, the state update method is basically a function converting an equation into abstract code, it does not have any state. The advantage of the fixed-per-run strategy is that it keeps the state updater very simple, we simply regenerate the state updater code object on every run. For a use case where you have a run for something like a couple of 100ms and then change a parameter, this would not introduce too much penalty, I guess.

I agree on allowing state updaters to modify namespaces. Actually, I thought we had this already but I guess I misunderstood. :)

No, currently the state updater only returns abstract code. I guess the cleanest (but this would have to be very well documented) would be if the state updater directly changes the namespace it's given as an argument?

from brian2.

This is not that different from what is done in Brian 1 where the ODE is integrated with an Euler integration method and a fixed number of time steps over [0, dt] (if I remember correctly). But with this method, we can control the precision and compute the appropriate number of steps as a function of the model parameters.

@rossant Alright, that's what I thought. There's indeed the get_linear_equations_solution_numerically method in Brian 1 which does 100 Euler steps. I always thought that this looked a bit hacky, having a principled way to chose the number of iterations sounds certainly like a better approach!

from brian2.

Yes, the 100 steps solution was never anything than hacky and supposed to be temporary! But, it's much better than using one of the nonlinear methods at least, so I don't feel too ashamed. ;)

The single equation thing is quite a sad limitation of sympy. Very much agreed that we want to make a very easy route for people to add new exact integrators as an easy way to contribute to Brian! Let's make sure the API is simple and well explained.

The more I think about it, the more I like the generate-state-updater-per-run technique. It's very simple to explain, easy to test for errors, not too restrictive and fairly efficient - seems like an excellent balance to me. I would like to suggest it to @romainbrette and others to see if anyone can think of any potential problems, but as far as I can see it sounds great.

from brian2.

I looked into the namespace issue again, and I think I'd now argue for the opposite of what I said before ;)

The reason is the following: Currently, the state updaters are nice and simple, they turn equations+namespace+specifiers (the simpler ones ignore namespace and specifiers entirely) into abstract code. If we now let them take care of modifying the namespace, they become a lot more linked with other parts of brian. For example, stochastic state updaters need the randn function which is not present in the equations. Adding randn to the namespace is non-trivial, though, as the function depends on the number of neurons in the group. So the state updater would have to add the following to the namespace as randn:

brian2.codegen.functions.numpyfunctions.RandnFunction(specifiers['_num_neurons'].get_value())

And it should also take care of the possibility that randn is possibly already included in the namespace, maybe raising a warning if it points to something unexpected.

Another issue is that a sympy-based state updater might introduce new functions -- we would have to go through the generated code to find them.

All this is about the function namespace. In some situations (e.g. in a numpy-based linear state updater) we might want to add variables (pre-calculated coefficient matrix) to the namespace, though.

So, how can get rid of the issue with the functions? We could revise my long-discussed implementation of the namespace interface between NeuronGroup and codegen (again)...

Currently:

1. We have a potential namespace, containing all units etc. in a CompoundNamespace object [NeuronGroup]
2. We resolve all the identifiers that are used for a code object, resulting in a namespace dictionary [NeuronGroup]
3. This resolved dictionary is given to the state updater (so it can replace constant values, for example)
4. The namespace dictionary is used for running the code [codegen]

Now, some issues would disappear when steps 2 and 3 would be exchanged. The state updaters wouldn't have to bother to add randn to the namespace: They'd just generate code containing randn, and the namespace is resolved for all identifiers in the code, including randn (and excluding any constants that are already replaced).

I think this would make things clearer and simpler (it would also allow for one additional small simplification in NeuronGroup), what do you think? One thing we would need for this is a reliable "identifier detection" function for abstract code. The current get_identifiers method in stringtools is not exactly what we want, as it would return also all the temporary variables, e.g. for _cond = v > 1 it would return {'_cond', 'v'}, but _cond does not exist in the namespace. So we only want identifiers that are not introduced in the code itself. I don't know the codegen code well enough, but this is something that is already solved for the translation to intermediate code, isn't it? Maybe this could be refactored into a general function and made part of stringtools that can then be used both by code generation and NeuronGroup (or ObjectWithNamespace)?

A final question: How do we do the additions to the namespace for variables (e.g. pre-calculated coefficient matrices)? We could have the state updater add things to the namespace object it receives (according to the above procedure that would be a CompoundNamespace object). I think I'd rather avoid this, as I would like to keep the NeuronGroup's namespace unchanged after its construction. Also, a stateupdater-specific variable shouldn't be in the namespace of the whole group, but only in the namespace related to the stateupdater's codeobject. The alternative would be to give state updaters a choice: Either they return abstract code (a string), or they return abstract code and an additional namespace (a tuple of string and dictionary). I think this is still quite simple and reasonable elegant.

Phew, long text, I'm sorry... hopefully I made my points more or less clear :-/

from brian2.

I have to admit I didn't 100% follow all the issues. I do agree that if we can do it, it's nicer to have the StateUpdater construction have no side-effects on the namespace to keep the code cleaner and more followable. If all that's needed is better identifier detection then we're OK, this is indeed already done in codegen and I could refactor this function. My only question is, is there any forseeable effect on other things we want to be able to do like: user functions, noise generation, etc.? I don't see any obvious problem but I might be missing something.

How to introduce additional variables for the pre-calculated coefficients is an important question. Previously I would have said that we had to add new state variables because we might want to change those pre-calculated coefficient matrices later. However, if we go with your idea of having StateUpdater computed once per run, then this isn't so important so we have a bit more flexibility on how to solve the problem. In this case, I'm happy with state updaters returning a pair of abstract code and additional namespace, I agree it's fairly simple and elegant.

from brian2.

If all that's needed is better identifier detection then we're OK, this is indeed already done in codegen and I could refactor this function.

This would be great (feel free to push changes to this branch) -- I think this would be a good approach, the abstract code implicitly defines what identifiers have to be present in the namespace or specifiers. This would save a couple of lines in NeuronGroup, generating resets, thresholds or the state updater would be even simpler than it is now.

My only question is, is there any forseeable effect on other things we want to be able to do like: user functions, noise generation, etc.? I don't see any obvious problem but I might be missing something.

I don't see any problems here, on the contrary: For example, replacing the random number generator with a different function would just mean providing a randn function in the namespace. If the state updater manipulates the namespace, it would have to check for an already existing randn function and deal with it specifically.

How to introduce additional variables for the pre-calculated coefficients is an important question. Previously I would have said that we had to add new state variables because we might want to change those pre-calculated coefficient matrices later.

Yes, for this use case, dealing with it via specifiers would be more flexible. On the other hand, we only have one specifiers dictionary for the whole NeuronGroup and conceptually, I don't think that a state updater coefficient matrix is something that should be part of the model (and has to be accessible from everywhere). The "resolved namespaces" are something that are bound to one code object, so that seems like the best place for them to live.

So it all boils down to the question how much flexibility we want to have. I'd be fine with the "state updaters are fixed per run" approach but if we want to have something more complicated, like two coefficient matrices per neuron (that are switched between active and refractory) then we need a different state updater system in general, one where the state updater is called at every timestep to have the opportunity to change matrices in the namespace.

I think I'll write to the mailing list for feedback from @romainbrette, I think he's not following his github notifications :)

from brian2.

Just wanted to add something I wrote on the brian2 list here regarding doing refractoriness with linear equations:

• The active vs. non-active coefficient matrices can also be implemented
  (albeit not very efficiently), by formulating the update using the
  is_active variable, so an update step would look on the lines of:
  S = is_active * dot(A_act, S) + (1 - is_active) * dot(A_ref, S)

I wonder how inefficient this is in practice? For example, if you had:

dV/dt=(ge+gi-V)/tau : 1 (active)
dge/dt=-ge/taue : 1
dgi/dt=-gi/taui : 1

In this case, the update step lines for ge and gi wouldn't have to depend on is_active because they have no dependence on V. The step for V would be just to modify the last line V += _V_dt to V += _is_active__V*dt.

So at least in this case, it wouldn't be too inefficient, in fact I think it would make no perceptible difference to the total computation time. The question is how efficient would it be in more complicated examples, and how easy would it be to formulate these little tricks?

I guess in general you would compute two matrices one for active, one for refractory, and if the coefficients were c_active and c_refractory then define c_diff=c_active-c_refractory then the general coefficient would be c_general=c_refractory+is_active*c_diff. When c_diff=0 you could remove this line. So the cost would be, for each coefficient where c_diff is nonzero, one additional read (of c_diff), one multiply and one add.

The optimisation where you can just replace V+=_V_dt to V+=_is_active__V_dt could also mean you could ignore some coefficients. For example, for _V above you ought to have two sets of coefficients, one for active and one for inactive, meaning you would get the read/multiply/add cost three times. However, since you know that the whole line _V=... is irrelevant (since it is multiplied by 0 later when doing V+=_is_active__V*dt) you can just use the active coefficients and ignore the refractory coefficients (the value of _V generated would be wrong for the refractory indices, but it doesn't matter).

from brian2.

One thing we would need for this is a reliable "identifier detection" function for abstract code. The current get_identifiers method in stringtools is not exactly what we want, as it would return also all the temporary variables, e.g. for _cond = v > 1 it would return {'_cond', 'v'}, but _cond does not exist in the namespace. So we only want identifiers that are not introduced in the code itself.

OK, I added an analyse_identifiers function to codegen.translation which does this - let me know if it works for you. It takes two arguments, the code and the already known identifiers, and returns three things, the newly defined identifiers, the known identifiers that were actually used, and any identifiers that were used but not known or defined in the code. The union of these three sets will be the output of get_identifiers.

from brian2.

OK, I added an analyse_identifiers function to codegen.translation which does this - let me know if it works for you.

Cool, thanks -- that should be all I need. This will again simplify the code in NeuronGroup a bit!

from brian2.

Closed by ef4a663, merging #32.

from brian2.