Code generation syntax for reductions?

[thesamovar]

One thing we would like to do in monitors and potentially in other situations is allow for reductions, e.g. mean(V) for monitoring. We already do this for summed synaptic variables, but we only allow the sum. In that case, it seemed sufficient, but for monitors and possibly more widely there might be a case for a more general syntax.

One option is to look at existing packages that include reductions, including Theano, numexpr and pycuda:

• http://deeplearning.net/software/theano/library/scan.html
• http://documen.tician.de/pycuda/array.html#module-pycuda.reduction
• http://code.google.com/p/numexpr/wiki/UsersGuide#Supported_reduction_operations

Theano goes for a complicated general system, pycuda just lets you write your own code basically, and numexpr has a small list of just two reductions. My initial feeling is that the pycuda option is the best with maybe a few defaults like sum and prod built in.

[mstimberg]

It's not exactly the same thing, but for parallel loops numba identifies reductions implicitly whenever inplace operations are used: http://docs.continuum.io/numbapro/prange.html

[thesamovar]

Interesting, but maybe difficult for us to implement? I wouldn't want to introduce a dependency on a non-free product (numbapro) even if they are giving away academic licenses.

[mstimberg]

Sorry, I did not mean that we should use numbapro for this task, I thought you were only listing syntax options that other packages use. Or do you actually want to use e.g. the pycuda package for the task?

[thesamovar]

Yep, understood. I was just thinking that if we were to use that syntax it might be rather hard for us to parse, and the first thought that came into my head was that they are already doing the parsing maybe we could use that. And then I thought, no! So basically you just got the last part of a stream of consciousness without any explanation. ;)

[thesamovar]

Christmas morning thought about code generation and reductions...

We are probably also going to need reductions and more for Brian Hears when we update that to work with Brian 2. In fact, for BH you have a wide variety of possibilities for each filterbank, i.e. N to M for any N and M.

So proposal: we make the idea of filterbanks more generic and core to Brian 2. Now, if you wanted to record the sum or mean of a variable, you would have the variable G.v, a filterbank F:v->sum(v)/N and a monitor recording the output of F.

There are other things you could do with this too. For example, the problem that Wieland just posted where he wanted only the neuron with the highest voltage to fire could be solved by having a filterbank F:v->argmax(v) and adding i==F(v) to the threshold condition.

What do you think?

To support Brian Hears, we would need some fairly detailed possible filterbanks. For example, the BH library makes regular use of Repeat and Tile. There's also N:1 mappings like sum obviously. I've also regularly used (for binaural models) something like f(x)=x[:len(x)/2]+x[len(x)/2:] which is a 2N:N mapping which is quite hard to specify with a nicer syntax than just giving the Python code.

[thesamovar]

On further thought, I think the way to do it is to have nice syntax for the simple, highly generic cases like repeat, tile, sum, and for complex cases like the last one above, you have to write code in the target language directly. So, it would be a bit like Functions.

[mstimberg]

Just a small remark: f(x)=x[:len(x)/2]+x[len(x)/2:] is equivalent to f(x) = x.reshape(2, -1).sum(axis=0), maybe we can actually find something highly general based on a two-step description with a first step of reshaping (including repeat/tile/...) and a second, optional step of reduction (e.g. sum).

[thesamovar]

Interesting. I think I'd like to stay away from reshaping because it means potentially reimplementing ND arrays in target languages. I wonder if it can all be done with index expressions though? For example, tile would be j % N, repeat would be j // N and so forth (here j is the target index, and the expression gives the source index). Not immediately sure how you fit reductions into that framework though. Maybe, inspired by your reshape idea, you could specify by an index expression which target it should be reduced into, so the binaural example would be i % N where i is the source index. A simple sum of everything would be the expression 0. I think this works logically, but not sure how easy it would be to implement efficiently.

[thesamovar]

And you could combine both types, as you suggest.

[thesamovar]

I realise in retrospect that essentially what I'm suggesting here is equivalent to Synapses but with expressions instead of arrays. I wonder if maybe we ought to support that in Synapses and then we get reductions and expansions for free. Some examples:

# sum all
G = NeuronGroup(N, ...)
H = NeuronGroup(1, ...)
S = Synapses(G, H, init='V_post=0', pre='V_post += V_pre', connect='True')
# binaural sum
G = NeuronGroup(2*N, ...)
H = NeuronGroup(N, ...)
S = Synapses(G, H, init='V_post=0', pre='V_post += V_pre', connect='i%N==j')
# tile, note that G and H are swapped in Synapses and we set the pre value not the post value
G = NeuronGroup(N, ...)
H = NeuronGroup(N*M, ...)
S = Synapses(H, G, 'V_pre = V_post', connect='i%M==j')
# repeat
S = Synapses(H, G, 'V_pre = V_post', connect='i/M==j')

The issue is still efficiency, but this does have the nice property that we don't need to create a whole new type of object. Well, it needs a lot of details working out but what do you think overall?

[thesamovar]

OK, thinking a bit more about reductions. I think several possible things come together: Synapses, continuous state monitors, filterbanks for Brian Hears, and possibly even more. What I'm thinking is that we have the concepts of expanding and reducing vector variables.

An expansion would include things like numpy repeat and tile, but in general it could be specified by any index expression source_index = f(target_index), so that repeat would be specified by source_index=target_index/M and tile would be source_index=target_index%M.

Summarising what we said about reductions I see three things:

1. Things like sum, max, etc. If we use @mstimberg's idea of x.reshape(...).sum(axis=...) then this is nicely general. Could be made efficient in parallel.
2. Synapses-like reductions as I was talking about in the pevious post, with index expressions or directly using arrays of indices. Not super-efficient in general.
3. Most generally, provide code for reductions, like in PyCUDA.

All of these can be easily implemented in both numpy/Python and C++.

Already, there's a lot of overlap with Synapses.

Using multiple stages as above, you can do the filterbanks you need for Brian Hears. By applying operations and then recording the output, you can do the generalised continuous state monitors we want in Brian 2.

It would also be nice to make this efficient by having a way of chaining operations into a single code block. I think this should be doable, although I haven't thought it through all the way. Suppose you had an expansion following by an operation followed by a reduction. I think this can be combined into a single loop in C++ where you iterate over the largest values that the indices take. I think we'll really want to get this right to make Brian Hears as fast as possible. It might be more complicated if you multiple stages of expansion and reduction - not sure.

OK so that's an outline of a way we could go. What do you think? It's quite an effort, but on the other hand it makes for a super-general system AND it would mean that we would get Brian Hears almost for free as a bonus.

[thesamovar]

In a meeting we came up with this:

Reduction syntax and continous monitoring:

• StateMonitor(G, record='sum(v**2)/N')
• StateMonitor(G, record='v**2', reduction='sum')
• StateMonitor(G, record='v2', reduction='sum')

Possibly, this is enough. We can already do, e.g. 2N:N using linked variables, and this takes care of N:1. Maybe not as beautiful, but does the job.

[mstimberg]

I know this is a very old topic, but some of this came up again in #1445, where @schmitts used a network operation to take the dot product between two neuronal variables (as a linear readout to get a scalar variable). I realized that a straightforward extension of our existing (summed) syntax could allow for (summed, shared) in neuronal equations, i.e.

eqs = '''
x : 1
y : 1
z = x*y : 1 (summed, shared)
'''

would have z store the dot product of x and y (not sure whether we need the extra shared – it is redundant but more explicit). I also thought that by making summed a special-case of more general reductions, we could get away by only changing the syntax for equations in NeuronGroup and Synapses, and leave things like StateMonitor untouched (it could simply record z in the above example). Here's what I could imagine: (summed) would be a special case of (reduction: sum) or (reduced: sum). I'm not sure "reduction"/"reduced" is the most intuitive term for less technical people not familiar with map/reduce, etc. It might also be something like "merge", "combine", "aggregate", etc. We could then (both in NeuronGroup and Synapses) have (reduce: max), etc., with a common basic template that initializes with the neutral value (0 for sum, -∞ for max, etc.) and then applies the respective function. It would be nice to support e.g. argmin, argmax as well, but it would be slightly more complex (the final value is different from the values that are compared). This would cover all N:1 mapping as part of NeuronGroup, and any(?) general N:M mapping could be done via Synapses. As a second step, we could then add convenience functions that hide a bit of the machinery (in particular for Synapses). The example that @schmitts contributed (https://brian2.readthedocs.io/en/latest/examples/frompapers.Nicola_Clopath_2017.html) would be a nice use case.

[mstimberg]

I came across this again via this post on the discussion forum. I think one additional advantage of the approach outlined in my comment above over an addition to StateMonitor is that we don't need to come up with a naming convention. For a StateMonitor that records v, we access the values with mon.v, but how would we access the recorded values for StateMonitor(G, record='v**2', reduction='sum')? If we only allowed reductions on simple variables, then we could probably access them as something like mon.v_sum, but it would still be less explicit than just defining a name in the NeuronGroup.