brian-team / brian2

Brian is a free, open source simulator for spiking neural networks.
http://briansimulator.org
Other
918 stars 218 forks source link

Improve safety of linear integration #679

Open thesamovar opened 8 years ago

thesamovar commented 8 years ago

At the moment, we have the situation that the exact linear integration method can output wrong results (see #626). In that case, it will output all nans so you know there is a problem. Is there a possibility that it can output finite values that are wrong? If so, that's a situation we need to fairly urgently address.

One option would be to reduce the scope of exact integration rather than using sympy's solver. This will mean some equations which have exact solutions wouldn't get found, but it would remove the possibility of a catastrophic failure. I think this would be worth it. We could also have a new method choice (in addition to 'linear' which would use the more conservative/safe strategy) which forces the use of sympy's solver (e.g. method='sympy').

A related idea which may or may not be feasible, would be to get the assumptions necessary to make the exact formula work. For example, in #626 that would be that the discriminant expression is positive (because we square root it). Is there any way to get sympy to tell us this? Or to infer it safely? If not using sympy, perhaps we could do this in our own routine? With that in place, we could check at runtime if the conditions hold. That would be straightforward for shared constants, but trickier if the conditions involved variables that can change during runtime or for individual neurons (i.e. parameters).

I think the assumptions idea isn't necessary if we can be sure that our choice of integrator won't give incorrect finite values, only nans, because if it always gives nans when there is a problem then #554 will tell the user there is a problem.

mstimberg commented 8 years ago

At the moment, we have the situation that the exact linear integration method can output wrong results (see #626). In that case, it will output all nans so you know there is a problem. Is there a possibility that it can output finite values that are wrong? If so, that's a situation we need to fairly urgently address.

I don't think it can give finite values that are wrong. The equation it produces are correct, the only problem is that it involves complex algebra (e.g. square roots of negative numbers) which will not work.

One option would be to reduce the scope of exact integration rather than using sympy's solver. This will mean some equations which have exact solutions wouldn't get found, but it would remove the possibility of a catastrophic failure. I think this would be worth it. We could also have a new method choice (in addition to 'linear' which would use the more conservative/safe strategy) which forces the use of sympy's solver (e.g. method='sympy').

Even if it does not give us incorrect results, I do agree to some extent that replacing the general solver by something more restricted might be a good idea. I've come across equations where it takes literally hours to compute the solution, and given that we always try it if the user does not specify the integration method manually, this can be a problem. I don't know what a good restriction would be that is itself straightforward to check, but I could for example imagine to by default only use it for equations with a single variable. This would still cover simple integrate and fire models, on the other hand it would not work anymore for anything involving synaptic currents...

We could have completely separate state updaters (something like linear and linear_full, or sympy has you suggested -- even though I would probably rather not use the name of the library, we are using it in any case. Maybe something like analytic), or handle this via an extra argument for the integrator (#481). Somewhat related, some problematic cases will be simpler if we freeze constants, which would also be specified by such an option.

thesamovar commented 8 years ago

OK so if it doesn't return incorrect results then #554 is good enough for now.

For the time problem though:

So does it take a long time to solve for some linear equations or does the issue arise when we try to use it to solve nonlinear ones? If the latter then we could just restrict to linear ones, but I guess that's not the issue given the name linear is in there.

I think the restriction to a single variable is too strong (e.g. wouldn't work for CUBA), but I'm not sure where the line is. Do we have some examples of when it takes a long time compared to when it takes a short time? Can we work out what a good criteria would be?

For example, if we are thinking of differential equations of the form x'=Ax for some vector x and matrix A then we could think about how large A could be before the general solution is too complex. I guess it'll get quite complicated even for quite small A, but maybe we can do better than dim A = 1?

For the name, I like linear for the default, and analytic to force an attempt to find an analytic solution no matter how long it takes.

thesamovar commented 8 years ago

I guess given that to find the eigenvalues of an NxN matrix you need to solve det(A-lambda I)=0 which is a degree N polynomial that we would only be able to do it for N<5 and probably even N=3 is going to be very complicated in general. Seems a shame to lose that possibility for cases like the CUBA example though, where it is very quick. The ideal would be the timeout idea but it seems getting that to work on Windows is not at all straightforward.

thesamovar commented 8 years ago

One idea: maybe we can allow N=2 for general matrices, or a larger N for upper triangular matrices (whose eigenvalues are just the diagonal elements). This would cover CUBA for example, and probably many others.

thesamovar commented 8 years ago

Just checking this quickly on maple, indeed if the matrix defining the ODE is general, even N=3 leads to a solution with more than 10,000 terms. However, for an upper triangular matrix N=3 is quite compact, N=4 is reasonable but large, and I didn't check N=5. So we could have that as our condition maybe? We pass to sympy if we detect that the differential equation has less than or equal to 2 variables, or if it has less than or equal to 4 variables but the matrix is upper triangular. Sounds reasonable? I feel like that would cover a lot of use cases.

thesamovar commented 8 years ago

See my latest comment on #626: it looks like we can sometimes give wrong results quietly.