-
Notifications
You must be signed in to change notification settings - Fork 216
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Implement exact integration for linear equations #13
Comments
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:
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 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 With the parameters in the above mentioned phase locking example, the linear state updater now generates the following abstract code:
There are a couple of technical difficulties/decisions to take now:
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:
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. |
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. |
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:
|
@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. |
@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 |
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
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.
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:
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? |
@rossant Alright, that's what I thought. There's indeed the |
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. |
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
And it should also take care of the possibility that 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:
Now, some issues would disappear when steps 2 and 3 would be exchanged. The state updaters wouldn't have to bother to add 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 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 Phew, long text, I'm sorry... hopefully I made my points more or less clear :-/ |
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. |
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.
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
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 :) |
Just wanted to add something I wrote on the brian2 list here regarding doing refractoriness with linear equations:
I wonder how inefficient this is in practice? For example, if you had: dV/dt=(ge+gi-V)/tau : 1 (active) 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). |
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. |
Cool, thanks -- that should be all I need. This will again simplify the code in NeuronGroup a bit! |
And integrate with the system described in #11
See the notes about several issues surrounding linear and "multi-linear" equations in the document refractoriness and linearity
The text was updated successfully, but these errors were encountered: