Modifying presynaptic variables doesn't work using numpy

[thesamovar]

Bug reported by Owen Mackwood on the Brian development list:

If you increment a presynaptic variable in the Synapse pre/post code, it is only incremented once per presynaptic neuron, rather than once per synapse (as it presumably should). The following code exhibits this behaviour:

import os
from brian2 import *

on_pre = True
standalone = True
if standalone: set_device('cpp_standalone')
else: prefs.codegen.target = 'numpy'

tau_A = 20*ms
gl = 10.0*nsiemens
el = -60*mV
vt = -50.*mV
memc = 200.0*pfarad
bgcurrent = 200*pA

eqs_neurons='''
dA/dt = -A / tau_A : 1
dv/dt=(-gl*(v-el)+bgcurrent)/memc : volt (unless refractory)'''
neurons = NeuronGroup(10, model=eqs_neurons, threshold='v > vt',
                      reset='v=el', refractory=5*ms)

kwds = {'connect':True, ('pre' if on_pre else 'post'): 'A_pre += 1'}
con = Synapses(neurons[:1], neurons[1:], **kwds)
sm = StateMonitor(neurons[:1], "A", record=True)

run(1*second, report='text')
if standalone:
    device.build(compile=True, run=True, debug=False)

plot(sm.t, sm.A.T)
show()

The plot should show increments of 9, whether on_pre is True or False. The numpy version exhibits increments of 1 in both cases.

[mstimberg]

This is obviously a serious bug. The "good" news: this is not a regression, the same problem occurs with the first alpha version.

[mstimberg]

This is actually quite tricky. We do a complicated algorithm for repeated indices on the post-synaptic side but nothing equivalent for the pre-synaptic indices. To do this properly we'd need to check whether pre- or post-synaptic variables are used in the code and accordingly partition synapses according to pre- or post-synaptic indices or both. This might get quite complicated... Probably it is better to directly go for a solution using ufunc.at (#161), we now require numpy 1.8, so that's not a problem. In examples as the one above, this is straightforward but there are tricky corner cases if the changed variable is used later on. In the issue #161, the idea was to use ufunc.at for all the cases where it is relatively straightforward and fall back to the current solution otherwise, but this is no longer an option since the current solution does not always give correct results. So, either we have to fix the current solution (which leaves the use of ufunc.at as a later option to increase performance) or we try to come up with a general solution based on ufunc.at.

I think I'd prefer the latter, I think we could come up with something reasonable that covers the corner cases but is not necessarily the most performant solution possible, e.g. we could do new copies whenever a variable is used on the RHS, e.g. something like:

A_pre += A_pre + w
w += A_pre

would become

A_pre = _array_neurongroup_A[_presynaptic_idx]
w = _array_synapses_w[_idx]
# This operation marks A_pre as "dirty", so it has to be reloaded when it is used later
add.at(_array_neurongroup_A, _presynaptic_idx, A_pre + w)
A_pre = _array_neurongroup_A[_presynaptic_idx]
add.at(_array_synapses_w, _idx, A_pre)

What do you think?

[thesamovar]

Does this problem also occur in Brian 1? If so, do we need to also fix the problem there and/or issue an advisory warning?

[thesamovar]

I agree it would be good to get a general solution using ufunc.at that works in all cases. There is one case that cannot work with ufunc.at though, where the user does something like v += f(v+w) for some arbitrary nonlinear function w. As I noted in #161 these cases are actually ill-defined (they depend on the order in which we consider the synapses), but users may write them. One option would be that if the user writes something like this, we issue a warning that their operation is ill-defined, and we fall back on a super-inefficient loop. I wouldn't find this solution too terrible because in the vast majority of cases the user won't do something like this, and in the small number of cases where they actually do want to, it will still be efficient for every target except numpy.

[thesamovar]

Who should work on this by the way? We need to decide soon because I'm on holiday next week and may not have internet access, but I have some long train journeys at the beginning and end where I could work on it.

[mstimberg]

There is one case that cannot work with ufunc.at though, where the user does something like v += f(v+w) for some arbitrary nonlinear function w. As I noted in #161 these cases are actually ill-defined (they depend on the order in which we consider the synapses), but users may write them.

Is this actually ill-defined? The approach I outlined above would interpret this as:

v_new = v_old + f(v_old + w)

by doing

v = _array_neurongroup_v[some_indices]  # this is a copy
add.at(_array_neurongroup_v, some_indices, f(v + w))

Isn't this the the correct interpretation, or am I missing something?

[mstimberg]

Who should work on this by the way? We need to decide soon because I'm on holiday next week and may not have internet access, but I have some long train journeys at the beginning and end where I could work on it.

I'd be happy for you to work on this, since you already did quite a bit of work in this direction in the ufunc_at branch.

[mstimberg]

Oh, and about Brian 1: it suffers from the same problem... For that, I'd probably just raise an error if pre-synaptic variables are changed in synaptic pre/post statements (and we should probably do a bug-fix release at some point...).

[thesamovar]

The new/old interpretation is one possible interpretation, and it would work for numpy, but it wouldn't work for weave (without having to make a copy of the arrays). I'm also not sure that it's consistent with our general semantics for abstract code, which is that a reference to a variable gives its current value.

OK for Brian 1: the error it raises could say "Please use Brian 2" once we've solved it for Brian 2. We should do an announcement though (even though it's painful to do so), maybe today actually.

[mstimberg]

The new/old interpretation is one possible interpretation, and it would work for numpy, but it wouldn't work for weave (without having to make a copy of the arrays). I'm also not sure that it's consistent with our general semantics for abstract code, which is that a reference to a variable gives its current value.

Ah ok, now I'm getting it -- I didn't really have the situation of repeated indices in mind, but this is of course what this is all about. I'm happy with any solution that correctly addresses the most common use case and in some way deals with any other case. From my side, I'd be even happy with a solution that just raises an error (instead of doing an inefficient loop). Can we find some general criterion for statements that are unproblematic? If we are very strict and only allow statements that do not depend on the order of evaluation (e.g. order of synapses) then something like v = w would be already forbidden. Do we also want to check the indices? Everything works without problems if we have at most one synapse per source and target (e.g. 1-to-1 connectivity)...

OK for Brian 1: the error it raises could say "Please use Brian 2" once we've solved it for Brian 2. We should do an announcement though (even though it's painful to do so), maybe today actually.

Given that the error has been around for quite a while and that changing pre-synaptic variables in a synaptic operation is quite rare, I'm not sure we really need an emergency announcement...

[thesamovar]

We have two issues, one is whether or not the statements can be efficiently implemented. So for this case, v=w would be allowed because it can be efficiently implemented in numpy. The other case is whether or not it's meaningful. For this case, v=w would be disallowed. I think all the cases that cannot be efficiently implemented are also not meaningful, so the efficient implementation check could be prior to the meaningfulness check. (Although we need to think carefully about this.)

One of the complexities is that people might write v = f(w)+v+g(w), v+=f(w)+g(w), v=v+f(w)+g(w), etc. These are all meaningful, and should be translated into v+=f(w)+g(w), but checking that they are is not entirely straightforward. I guess it's not too much of a problem if we only efficiently deal with v+=... or v*=... expressions, and raise a warning for the other cases (that mentions that you should rewrite as v+=...).

I do think we should make an announcement about the bug since it silently gives wrong results. We can say that it is unlikely to affect anyone. We've done this twice before: https://groups.google.com/forum/#!searchin/briansupport/bug$20dan$20goodman$20-%22re%22/briansupport/R6Bs4EDxuNM/reMIFFLkLNQJ and https://groups.google.com/forum/#!searchin/briansupport/bug$20dan$20goodman$20-%22re%22/briansupport/vOGwh3YHrbY/Hey9qXpNA2wJ

[thesamovar]

I'll do the announcement later unless you have strong feelings against it.

[mstimberg]

I'll do the announcement later unless you have strong feelings against it.

No, of course not, please go ahead.

[thesamovar]

Done!

I meant to add: I don't think we should check the indices.

[thesamovar]

OK I'm trying to work out what the conditions should be for meaningfulness and for vectorisability with ufunc.at, maybe you could take a look at the reasoning below and see if you agree or if I've missed something.

Let I_syn be the synaptic indices, and I_pre, I_post be the corresponding pre- and post-synaptic neuron indices.

The definition for meaningfulness of an abstract code block is that permuting I_syn doesn't change the results.

If we look at each line of abstract code that has the form var op expr we can analyse the following cases.

If var is a synaptic variable, then expr must be made up of constants, synaptic variables, and non-synaptic variables whose values are permutation-dependent. So for example, if it reads a non-synaptic variable that is never written to, that value is not permutation-dependent. So u_syn = v_pre is allowed if v_pre is never written to. However, if we had v_pre = 0; u_syn = v_pre this would also be OK, because even though v_pre is written to, it takes a permutation-independent value. Something like v_pre += w_syn; u_syn = v_pre is no good though, because the value of u_syn will depend on the order in which the synaptic indices were looped through. This reasoning holds both for meaningfulness and for vectorisability with ufunc.at.

If var is a non-synaptic variable, then there are two cases depending on the operation.

First case: the operation is commutative (i.e. +=, *=, ...), or equivalent to a commutative operation (e.g. v = v+w is equivalent to v += w). In this case, if the expression only contains constants, synaptic variables, or any non-synaptic variable whose values are permutation-independent, then this is OK. It is meaningful because the value of the expression is permutation-independent and the commutativity of the operation means the modified value of var will also be permutation-independent. In this case, the code is meaningful (because the final value is permutation-independent) and vectorisable (because we can replace it with a ufunc.at expression). The existence of this case makes the value of var be non-permutation independent with respect to the other operations being considered here (e.g. it could not appear on the right hand side of u_syn = ...).

Second case: the operation is non-commutative (i.e. op is =). In this case, the expr needs to consist only of constants or variables that have the same index as var and are permutation-independent. In this case, the code is meaningful (it's permutation-independent because it only depends on indices belonging to the variable), and it can be vectorised in a separate block from the synapses. The existence of this case leaves var being permutation-independent (e.g. we could write v_pre = 0; u_syn = v_pre).

I think that covers all cases, but the logic is quite difficult to work through so I'd appreciate it if you could go through it carefully and check if you agree. If I'm right then I think there is a reasonably straightforward algorithm for checking if its meaningful and vectorising it based on the reasoning above.

Do we want to allow for optimisations in the case that something can be efficiently vectorised but isn't meaningful (e.g. v=w works because in numpy v[[0,0,0]]=[1,2,3] will give v[0]==3 which is what you'd get if you looped)? My feeling is that this isn't important to do.

[thesamovar]

There are some cases that are missed out by the analysis above. For example, if we wrote v_pre = 2*(v_pre+w_syn) then this is meaningful (permutation-independent) but vectorising it is non-trivial. We'd need to write add.at(v_pre, indices, w_syn); v_pre *= 2**len(indices). More complicated cases would include v_pre=v_pre*w_syn*(v_pre+w_syn). Note that we could write these two cases also as v_pre+=v_pre+2*w_syn and v_pre*=w_syn*(v_pre+w_syn) so restricting to statements with the += and *= operators doesn't help us here. In general, any symmetric function f(a,b) where f(a,b)==f(b,a) will let us write code of the form a=f(a,b) and it will be meaningful (but not necessarily vectorisable). This includes all functions of the form f(a,b)=p(ab,a+b) for any polynomial p. Note that this covers all f that are polynomials, but there will be other cases, e.g. f(a,b)=g(a+b) for an arbitrary function g.

So the question is: do we want to worry about any of this or not? My feeling is no, that all examples other than x += expr or x *= expr are very likely not meaningful (and even x *= expr is probably not very important).

So I think the reasoning in my previous post gives the condition for vectorisability of meaningful code using ufunc.at, but it misses out (a) non-meaningful code that can be easily vectorised (e.g. v=w), and (b) meaningful code that cannot be easily vectorised (as in the examples in this post). We could probably check for both of these conditions and either handle/warn about them, but perhaps it's not necessary? Thoughts? (Checking for meaningfulness is actually non-trivial, although sympy might be able to help here. It means checking for a user given function f(a,b) whether f(a,b)==f(b,a) for all a,b.)

[mstimberg]

I need a bit of time to think about this in detail. Actually maybe something like a wiki page might be better to put the information together? For now, I have the feeling that it might be better to come up with something less general and more restrictive and rather enhance the scope when we find new use cases (but then again one of the ideas behind Brian is that we support things that we didn't think of, so...)

A more general thought: maybe we should allow statements that depend on the order of synapses? I don't quite remember where, but I do remember that we once had a discussion about the order of synapses and we decided to follow some simple rules (something like: connect statements with concrete source/target arrays will create synapses in the specified order, statements with string expressions will create synapses in some specified order (ordered by pre- or by post), several connect calls will simply result in concatenated synapses). This means that abstract code is slightly less abstract than we want it to be in general and that synaptic propagation cannot be easily paralleled for OpenMP/GPU, but maybe this is more robust and easy to explain than a complex set of rules such as the ones you outline above? The reference implementation would be a single loop over the indices, i.e. what weave does. We'd make sure that numpy matches this.

I'm not sure about this, but I have brian2hears in mind where linked variables are used in a way that depends on the order of execution which is not well-defined according to a strict interpretation of abstract code. The issue is basically the same as for synapses (you can re-create synaptic pathway statements -- except for delays -- with linked variables and custom operations, without using any Synapses object). Then again, I think the need for this in brian2hears only arises because we don't have any support for "buffer" arrays, using such a data structure instead would also be more efficient (the current solution unnecessarily copies over a lot of data every time step).

[thesamovar]

I think it would be good to have buffering be more explicit for Brian Hears 2, like in BH1. Let's discuss that another time though.

I quite like having the freedom to reorder synapses if necessary (and this could be very necessary on GPU as you say, and on SpiNNaker we can't make any guarantees about the order in which synapses are considered). So I think I'd prefer to stick to the current approach of warning if the ordering makes a difference, but not insisting on it and not guaranteeing that they'll be evaluated in a particular order.

I've started work on this in a new branch. I'll make a pull request for discussing the implementation.