Implement Synapses

[mstimberg]

This is of course a big chunk of work but I think we should rather start it sooner than later because it is a core component and probably a lot of unforeseen problems will pop up during its implementation.

I had a quick look through the code in Brian 1 and will use this issue to jot down some observations:

• The good news: A lot of code in Synapses works around restrictions in Brian 1 and is no longer necessary in Brian 2. For example, we don't need a specific SynapticEquations class (it was only introduce to allow for the (event-driven) keyword, support for _pre and _post variables can be easily implemented using the specifiers, no need for implementing them via static equations, etc.
• We need to put some thought into connecting synapses stuff to the specifier/language/devices system, the SynapticVariable class, for example, is in fact very similar to a specifier in the current system (this also entails that Synapses can use the standard __getattr__ mechanism of Group).
• A major point is how the SpikeQueue fits into the current system: It could quite naturally be handled in a Devices system (see discussion in #23 ) because it is about storing data. On the other hand, we probably want to have a C version available in the standard runtime version for performance reasons, this would rather suggest to connect it to Language (i.e. automatically use the C version if we use C code generation).
• Another major point is the synapse creation: Should this use code generation? If yes, if we use C++ code generation but synapses are stored in a dynamical array in Python, how would that work? Should we directly access the dynamical array in C or rather return an array of "synapses to create"?

[thesamovar]

I absolutely agree, July is not so far away and we need to have Synapses working by then.

We need to put some thought into connecting synapses stuff to the specifier/language/devices system

I think this will end up being the bulk of the work (especially if we include synapse creation as below).

A major point is how the SpikeQueue fits into the current system

I think it should go in Device and should work like in Brian 1, namely if you do a pure Python install you get the pure Python version, if you do a C++ install you get the C++ version (and we try to detect this at install time). This code is static and doesn't need codegen, so it shouldn't go in Language. (Incidentally, should Language be renamed CodeGenerator or something like that?)

Another major point is the synapse creation: Should this use code generation? If yes, if we use C++ code generation but synapses are stored in a dynamical array in Python, how would that work? Should we directly access the dynamical array in C or rather return an array of "synapses to create"?

Good point, and tricky. If we want to return variable sized arrays from C++ we can do it, we simply have to do the same thing we do for threshold, namely return a numpy array. What we would do is use an std::vector (which is a dynamic array) to build the synapses, and then return a numpy array. This returned numpy array could be converted to a DynamicArray without a memory copy.

Global question: how do we want to start implementing this? Do we copy Brian 1 code across or rewrite from scratch? In principle I think rewriting from scratch is the better way, as we've done with the rest of Brian 2, but on the other hand there is probably some stuff we don't need to modify very much in Synapses (like SpikeQueue for example).

[mstimberg]

I think it should go in Device and should work like in Brian 1, namely if you do a pure Python install you get the pure Python version, if you do a C++ install you get the C++ version (and we try to detect this at install time).

But there would still be the issue of how to access this in C++ generated code (I vaguely remember we discussed this already at some point): Otherwise i would access the queue via the weave API that talks to a Python object that itself wraps around a C-object instead of directly accessing the C structure.

Global question: how do we want to start implementing this? Do we copy Brian 1 code across or rewrite from scratch? In principle I think rewriting from scratch is the better way, as we've done with the rest of Brian 2, but on the other hand there is probably some stuff we don't need to modify very much in Synapses (like SpikeQueue for example).

I'd say we start from scratch and copy over things individually that can be copied, in particular the SpikeQueue, but also some helper functions, etc.

[thesamovar]

To recap the meeting: for SpikeQueue if we implement it in C++ and connect to Python via SWIG then hopefully we would also be able to directly access it in weave directly without going through Python (although we might be able to be a bit cleverer about it and avoid some memory copies). We can then also use the same code directly in the standalone version too.

The major issue in implementing Synapses would be the construction of synapses because you have to return a dynamic array. Since stl::vector gives you a contiguous array we should be able to construct like this.

What we still need to think about is how to efficiently construct the arrays when the synapses are constructed in several steps. We shouldn't be afraid to make big changes to the way it is currently implemented (which may not be as efficient as it could be).

[mstimberg]

Ok, I started a synapses branch and you can already create a Synapses object -- it's not doing anything, though ;)

I had a look at the Synapses code, here are some things that might be worth discussing:

• For consistency with the way we are handling stuff in NeuronGroup, I think the compress method (where the dynamic arrays are replaced with normal arrays) should disappear. Instead, we should do the conversion in pre_run (together with the namespace resolution). This would have the benefit that one could add synapses after a run.
• Maybe we should reconsider some of the syntax/design decisions in Synapses (but this is something we certainly have to discuss with Romain and/or Victor): For example, I think it might be really useful in some cases to delete synapses, e.g. after learning a while one might want to prune weak synapses (this could of course be done by setting their weights to zero, but they would still be present and take time/memory). Or even more simply: During the creation phase one might want to create synapses between all neurons (or according to some criterion) except for some rare combinations. In this case it might be simpler to apply the general rule first and then remove some synapses. If deleting were possible, I'd also say that we should no longer do the "additive assignment".
• Also a syntax/implementation question: The three dimensional array access is quite cumbersome at the moment, maybe we can find a way to make this work more general (without adding a lot of special case code for all combinations of indices and slices)
• One thing we have to keep in mind: Since the number of synapses is not fixed (in contrast to the number of neurons in NeuronGroup) we have to create the "core namespace" (that provides the functions and units) in pre_run and not during object creation because we have to get functions like randn() provide the correct number of values.
• Not a core thing but also worth considering: Maybe we should not require the source group to be a spike source. I'm mentioning this because I know users who tried to connect a PoissonInput to synapses in Brian 1 -- this does not work as the PoissonInput object does not spike but directly alters variables in the target group. Maybe we should introduce a rule like "if pre (resp. post) code is provided, the source (resp. target) group has to be a spike source". Or actually we wouldn't need any source in this case at all, as the PoissonInput class "pushes" changes, Synapses doesn't actually have to know about it?

[thesamovar]

OK, I've finished going through the code for Synapses, which is something I probably should have done a long time ago. Quite a few observations so I might split it up into several posts.

I'll start with some syntax issues. I think there are a few things where we might make the syntax a bit more coherent.

1. Lumped variables. I find the syntax group.var = synapses.x to be non-obvious (it means sum the corresponding synapses each time step). What about group.var = synaptic_sum(synapses, 'x') instead? Or something like that? The reasoning being that we want everything to be explicit and readable.
2. Creating synapses. We have S[i,j]=True and S[i,j]=3 to create 1 and 3 synapses. I don't find this particularly intuitive. What about S.num_synapses[i,j]=1 or S.num_synapses[i,j]=3 instead? Or S.create_synapses(i,j,3). Do we want to keep something like the slice based syntax S[:,i]=True? One possibility is to say no, we reserve the slice notation for setting values. I don't have a particularly good idea yet, but I feel like the current one is too non-obvious.
3. Accessing variables. I think the current notation S.w[i,j,k]=... is OK, but maybe we need to make it clearer and more explicit in a document? Also, for the string based assignment, as well as symbols i and j, perhaps we should have k meaning the synapse number if there are multiple synapses?
4. Multiple pathways. Currently this is implemented within a single Synapse object using tuples of arguments. This makes the delay notation a bit confusing. There is some discussion in the document suggesting multiple Synapse objects but that it would be too complicated. That was probably true in Brian 1 but maybe not in Brian 2? i.e. we could have multiple Synapse objects that share some state?

[thesamovar]

Oh, I forgot to mention, there is this notation S[:,:]=float_val which does a probabilistic assignment. How about we get rid of that and just use string assignments or a method in this case? It feels confusing.

[thesamovar]

Implementation of SynapticVariable. I quite like the way it goes through the Synapses.synapse_index method and find it pretty well designed. Can it be refactored into specifiers more generally?

[thesamovar]

Some feature questions. Agreed with you that having deletion would be nice. We need some thought to decide if it is computationally feasible. I think with the current data structures, it would be horribly slow but maybe there is another data structure we could use?

There is also this notation that if you do S[i,j]=True and then the same again, you will get two synapses. Again, this is non-obvious. It would make sense if you had S.create_synapses(i,j) repeated twice though, so that might be an alternative.

[thesamovar]

Finally, there are implementation issues.

First, for SpikeQueue rather than having an API that is based on returning arrays, I think it would make sense if there were two APIs, one that returned arrays for use in vectorised languages, and another that is essentially an iterator for use in C++ or Java, so that you can avoid the creation of these large temporary arrays.

Second, there seem to be some HUGE inefficiencies in Synapse.setitem and Synapse.create_synapses, namely at one point there are some dictionaries created by doing zip on two arrays. I think these data structures will be very large, and created very slowly. Hopefully we can find ways to fix this.

[mstimberg]

Ok, first my comments on the creation syntax: I share your feeling that using the setitem syntax directly on the Synapses object is confusing. One major problem I see is that it introduces an assymetry between setitem and getitem:

>>> S[:, :] = True
>>> print S[:, :]
...
TypeError: can only concatenate tuple (not "int") to tuple

>>> print S[:]
Size 200 subgroup of Synapses
>>> S[:] = True
...
AttributeError: Synapses behave as 2-D objects

So one does one-dimensional getitem access to get a subgroup and a two-dimensional setitem access for creating synapses. Related to that, I'm regularly making the mistake of typing S[:, :].variable = ... instead of S.variable[:, :] = ... which leads to the nice "can only concatenate tuple to tuple" error message -- at least this message has to be more helpful. But then, I guess we don't allow subgrouping Synapses in the first place (this was a side-effect of subclassing from NeuronGroup) in Brian2 -- in your current implementation subgrouping is for spike sources and a synapse doesn't generate any spikes.

I think I'd rather opt for keeping something with a slice based syntax, i.e. rather S.num_synapses[i, j] = 1 and not S.create_synapses(i, j, 1). Otherwise, how would setting connections via strings fit in? Even if we would allow setting a string for all synapses via a special syntax (say 'None' as an index), I think that (given a presynaptic group size of 50 and a postsynaptic group size of 100)

S.num_synapses[:, :50] = 'i==j'
S.num_synapses[:, 50:] = '2*(i==j)'

is considerably clearer/more concise than

S.create_synapses(None, None, '(j<50)*(i==j) + (j>50)*(i==j-50)')

On the other hand, for standalone the latter would be more useful :-/ Currently, the syntax also allows having (sub)groups as indices, so one can do stuff like:

G = NeuronGroup(100, ...)
E = G[:50]
I = G[50:]
S = Synapses(G, G, ...)
S[E, E] = ...
S[E, I] = ...
S[I, E] = ...

If we use a S.num_synapses[..., ...] syntax, then the 'creating = adding' semantics would make even less sense than at the moment. If it's not possible to implement in a reasonable way, I'd rather say we should raise an error if the number of synapses is reduced for a connection. If we implement getitem differently, then one could still add synapses via S.num_synapses[:, :] += ... I guess.

About boolean vs. integers vs. floating point values: I'm a bit undecided here. I think that having True/False has to stay -- if we allow an expression like 'i==j' we also have to allow 'True', and if we allow that then we should also allow True (i.e. the value, not the string). I actually wrote a paragraph for our paper yesterday where I argued why boolean/integer/floating point is not a problem: In the two cases for ambiguities can occur (True/1/1.0 and False/0/0.0) all the interpretations are identical. But I wouldn't mind getting rid of the floating point shortcut, writing 'S[:, :] = 'rand() < 0.1' instead of S[:, :] = 0.1 is certainly more explicit.

[mstimberg]

Ok, on to the remaining Syntax issues:

Lumped variables

I agree that an explicit syntax would be better here. One reason is that maybe in the future we might want to support some non-linear summation. Another reason is the implementation side: To support the current syntax, synapses returns a SynapticVariable object instead of the values -- this is different from state variables in NeuronGroup. To get the values in NeuronGroup you can use print G.variable but in Synapses you have to do print S.variable[:]. Now, there is also some advantage to the Synapses way of doing this, as for standalone we could then handle state variable assignments such as G.v = G.vr -- but either way it should be consistent between NeuronGroup and Synapses.

Variable access

Accessing variables. I think the current notation S.w[i,j,k]=... is OK, but maybe we need to make it clearer and more explicit in a document? Also, for the string based assignment, as well as symbols i and j, perhaps we should have k meaning the synapse number if there are multiple synapses?

What I don't like about S.w[i, j, k] is that it pretends the weight matrix to be a three-dimensional array, which it isn't (the last "dimension" can vary across connections). If one has full connectivity with 2 synapses each, one would expect that S.w[:, :, 1] = ... would set the value of the second synapse in all cases but this doesn't work. From an implementation point of view, the only "correct" solution would be to have a two dimensional matrix, containing an array of values for each connection. So one would have to write S.w[i, j][k] = .... I think this would make the syntax clearer (albeit uglier) in the multiple synapse case. And we can use S.w[i, j] = ... for single synapses (and in general this would broadcast it to all synapses).

In general I'm in favor of having "k" as another index for string-based assignments. We should keep in mind that we are currently preventing users from using all special variables everywhere, though, so that would entail that one could not use yet another simple variable in neuronal equations.

Multiple pathways I didn't think deeply about this yet but I'd be very happy to get rid of the pre=tuple syntax -- it also makes the implementation a bit ugly because we always have to deal with lists for the pre code, even though it will have a single entry in 99% of the cases. Sharing state variables by itself is not very difficult to implement now with the specifiers system. What needs a bit more thought is how we would make the connection: Not all variables are to be shared (delays may be different) and the two Synapses classes need exactly the same connection pattern. I guess we would need a special syntax to declare shared state variables (e.g. variable : unit (shared) and would have to make one of the objects the "master" object where all the synapse creation takes place. But my feeling is that this is worth the effort and clearer then the current list of pre codes/list of pre delays approach. Another advantage is that we could implement it incrementally: We could start with a version that does not allow for multiple pathways and later add the feature without changing the core. Maybe we could even put it in a completely separate class -- I think it's possible to implement it in a separate class without the Synapses class even knowing about this feature at all.

[mstimberg]

Let's talk about implementation! I just pushed changes that actually show a first running implementation of Synapses -- and it does not even use a lot of ugly hackery (you can run examples/synapses.py, it should work in Python and C++). There are of course plenty of features missing, I avoided everything we are still discussing here (most importantly, creating synapses and setting state variable values) but it still does something: It integrates the synaptic equations and applies the pre-code with the correct delays on pre-synaptic spikes (I did not check post codes yet, but in principle they should work). The SpikeQueue i ported as it is, only adapted it a bit to our framework. I removed everything around the C version, weave and dynamically updated delays for now.

Creating synapses is currently only possible via a connect_one_to_one method which manually sets the indices (which is straightforward in that particular case). Apart from the great feeling of finally seeing the effect of a spike on a postsynaptic cell in Brian2, implementing was quite helpful in figuring out where work is necessary to make things fit nicely together:

• The synaptic code templates used a synapse_idx: This is better to read but if we bind synaptic state variables to this index, we cannot use the state updater template (which uses neuron_idx) for integrating the variables. For now, I replaced synapse_idx by neuron_idx in the templates but a nicer solution would be to use more general term, e.g. group_idx everywhere. Same for _num_neurons.
• Accessing the state variables of the pre- and postsynaptic NeuronGroup in the Synapses class didn't work as straightforward (just using the specifiers from the source/target group) as I hoped: We have to replace the index information in the specifier -- I did this by creating a new Specifier object which copies everything and replaces the index but I guess this could be done more elegantly.
• The C++ code does not currently work with dynamic arrays, I therefore replace all dynamic arrays by numpy arrays in pre_run and change the information in the specifiers. I guess we need a new DynamicArrayVariable specifier to handle this properly.
• Again related to specifiers: I'm currently injecting the spikes array into the namespace at every timestep but this would be more elegantly implemented via a specifier. The problem is that we want the same symbol _spikes to stand for either pre- or postsynaptic spikes depending on whether we are dealing with pre or post code. This does not work with the current system where GroupCodeRunner pulls all the specifiers out of the specifiers in the Synapses class. This is not a fundamental problem, it only needs refactoring of the code around GroupCodeRunner. Alternatively, we could use two templates (which would use the same basic template and just replace _spikes by _pre_spikes or something like that).

All in all, it went much more smoothly than I expected, our modular structure really starts to pay off!

[thesamovar]

But then, I guess we don't allow subgrouping Synapses in the first place

I agree, I don't see the point of subgrouping Synapses.

I think I'd rather opt for keeping something with a slice based syntax

Agreed. The only concern here is that you cannot do create_synapses but have to directly set num_synapses, which is odd given that most of the time users will probably only create one synapse. But I think it's OK and it is at least very explicit.

If we use a S.num_synapses[..., ...] syntax, then the 'creating = adding' semantics would make even less sense than at the moment. If it's not possible to implement in a reasonable way, I'd rather say we should raise an error if the number of synapses is reduced for a connection. If we implement getitem differently, then one could still add synapses via S.num_synapses[:, :] += ... I guess.

Agreed with all that!

About boolean vs. integers vs. floating point values

OK, let's see about that. For floating points, what do we do if we have e.g. 1.5? I think we actually discussed this when originally designing Synapses and I agreed with the logic of float values then, but I feel a bit more dubious about it now.

Now, there is also some advantage to the Synapses way of doing this, as for standalone we could then handle state variable assignments such as G.v = G.vr -- but either way it should be consistent between NeuronGroup and Synapses.

Well, if SynapticVariable and NeuronGroupVariable were derived from ndarray we could potentially have the best of both worlds?

What I don't like about S.w[i, j, k] is that it pretends the weight matrix to be a three-dimensional array, which it isn't (the last "dimension" can vary across connections).

In a way, the same could be said of the first two dimensions given that it is a sparse matrix. Logically speaking, you could say that if there are 3 synapses say, then S.w[i, j, 3:] are all zeros in a 3D sparse matrix.

If one has full connectivity with 2 synapses each, one would expect that S.w[:, :, 1] = ... would set the value of the second synapse in all cases but this doesn't work.

But can we make it work?

I'm thinking we make the getitem/setitem syntax work as follows: when you do e.g. S.w[i,j,k] (where here i, j and k can be arrays, ints or slices) it constructs a subset X of all synapse indices which match the criteria. The order of synapses in X is undefined as far as the user is concerned, but they can extract it via, e.g. S.i[i,j,k], S.j[i,j,k] and S.k[i,j,k] if they need to know. These would be the values substituted into a string expression if it was used, so we already have to have this functionality. If some of i, j or k is missing they are interpreted as : or slice(None). Slices are interpreted as matching any synapses that exist, whereas arrays or ints have to exist and an error is raised if they don't. (What about when there are arrays and slices? Need to think about that.) With this syntax, I think the 3D indexing is coherent and explainable. What do you think?

Maybe we could even put [multiple pathways] in a completely separate class -- I think it's possible to implement it in a separate class without the Synapses class even knowing about this feature at all.

I like this idea!

[thesamovar]

Let's talk about implementation! I just pushed changes that actually show a first running implementation of Synapses

Wow! Incredibly cool. Well done!

Small note: I might create a new subpackage synapses rather than putting it in groups. Technically Synapses is a Group (although SpikeQueue certainly shouldn't be there), but even so I think it warrants its own place.

a nicer solution would be to use more general term, e.g. group_idx everywhere

Yes that seems logical. Alternatively, maybe there is a way of redesigning codegen that makes this problem go away.

I guess we need a new DynamicArrayVariable specifier to handle this properly.

Yep. Hopefully that shouldn't be too tricky.

All in all, it went much more smoothly than I expected, our modular structure really starts to pay off!

Very, very cool indeed. Thanks for doing all this work!

[mstimberg]

Ok, syntax discussion first:

Agreed. The only concern here is that you cannot do create_synapses but have to directly set num_synapses, which is odd given that most of the time users will probably only create one synapse. But I think it's OK and it is at least very explicit.

Ah, I see -- I did not quite get that part of the argument before. I agree that most of the time users will only want one synapse between two neurons, but on the other hand they will also rarely connect single neuron pairs, except if they do excessive looping which they should avoid. We have to think about how other ways of creating synapses fit in as well, e.g. setting directly via a (sparse) matrix.

OK, let's see about that. For floating points, what do we do if we have e.g. 1.5? I think we actually discussed this when originally designing Synapses and I agreed with the logic of float values then, but I feel a bit more dubious about it now.

I share your feeling... One could argue for 1.5 being interpreted as 1 synapse + 50% chance for a second one. I think that would lead to a consistent interpretation at all points (except for adding up values, e.g. S[:, :] = 0.5;S[:, :] += 0.5 isn't the same as S[:, :] = 1.0). But then, this is isn't entirely obvious and the overhead for just writing rand() < 0.5 is not really big.

Well, if SynapticVariable and NeuronGroupVariable were derived from ndarray we could potentially have the best of both worlds?

That might be quite a bit of work to get right... We have to be very careful about the reference/value semantics here. For example, should v = G.v; v[:] = 0 overwrite the state variable values (I'd say no)? Do we want G.v += 1 to work (I'd say yes)? We will also not be able to keep the subclass information in all cases and it might therefore be better to have one system that always works (setting via strings) and not a second one that sometimes works (setting directly). Apart from the implementation issues we should decide whether it is useful to "see" directly that G.v is a state variable and not an arbitrary array or are we happy with the "group reference + string name" notation. If we decide to change it, it would be only logical to replace StateMonitor(G, 'v') by StateMonitor(G.v), actually...

What I don't like about S.w[i, j, k] is that it pretends the weight matrix to be a three-dimensional array, which it isn't (the last "dimension" can vary across connections).

In a way, the same could be said of the first two dimensions given that it is a sparse matrix. Logically speaking, you could say that if there are 3 synapses say, then S.w[i, j, 3:] are all zeros in a 3D sparse matrix.

I don't quite agree: You could speak of the first two dimensions as a sparse matrix, but that's an implementation detail that the user doesn't have to care about, all he needs to know is that it has NxM dimensions. You can't really speak of the three-dimensional structure as a 3D sparse matrix because even a sparse matrix has a fixed number of dimensions. Trying to present it in this way will always lead to some inconsistencies. But I think your proposal (treating all indices as "search parameters" and returning a flat list of results) is consistent and would work, in particular it would deal nicely with the "set values to a scalar" case. It sacrifices the 2d structure, though, and makes it difficult to determine where values provided at the RHS end up. But maybe there's no way around it and we should only guarantee that the order is as expected if one slice is used, i.e. S.w[:, j, k] = should go in order over the presynaptic cells, S.w[i, :, k] over the postsynaptic cells and S.w[i, j, :] over the synapses. Everything that involves more than one slice has to do something with S.i, S.j, S.k as you proposed (or probably just use a string instead). Well, during the course of writing I got somewhat convinced by your idea :) Maybe a little change:

Slices are interpreted as matching any synapses that exist, whereas arrays or ints have to exist and an error is raised if they don't. (What about when there are arrays and slices? Need to think about that.)

I think there's no need to treat slices vs. ints/arrays differently: Just try to match and raise an IndexError if no match is found (I think that would even ensure Slice(start, stop, step) == np.arange(start, stop, step) equivalence, where stop would be the current maximum number of synapses per connection if not provided).

What I'd suggest for all syntax questions (not only Synapses, but also refractoriness, etc.): Let's put up a wiki page with concrete examples of how it does/should look like -- this will also be a nice "executive summary" for Romain and we should then be able to decide upon these matters quickly.

[mstimberg]

Small note: I might create a new subpackage synapses rather than putting it in groups. Technically Synapses is a Group (although SpikeQueue certainly shouldn't be there), but even so I think it warrants its own place.

Agreed. I thought about putting SpikeQueue into brian2.memory at first but it's a bit too Synpases-specific for that, so a package for the two modules makes sense.

a nicer solution would be to use more general term, e.g. group_idx everywhere

Yes that seems logical. Alternatively, maybe there is a way of redesigning codegen that makes this problem go away.

Yes, we could provide arguments to the template that set the index -- on the other hand, we never deal with synapses and neural code in the same file/function so there's not much room for confusion. Let's think of this during the template redesign in #30.

I guess we need a new DynamicArrayVariable specifier to handle this properly.

Yep. Hopefully that shouldn't be too tricky.

Don't see much of a problem, maybe it needs some small change in codegen as well.

[mstimberg]

I guess we need a new DynamicArrayVariable specifier to handle this properly.

I implemented this now, it's almost trivial: The DynamicArrayVariable only overrides get_value in ArrayVariable with return self.array.data -- the only important thing is that code generation always uses the get_value method when building its namespace and not directly accesses the array attribute.

BTW: Will we ever need a dynamical array with more than one dimension in our current design where we no longer have the state matrix _S? If not, we could merge DynamicArray and DynamicArray1D and probably simplify the code quite a bit.

[mstimberg]

Oh, and one more question: SpikeQueue is now a BrianObject and gets the spikes from the presynaptic group in its update method. Currently we only have the threshold and the synapses slots in the schedule -- I currently let the SpikeQueue update in the synapses step, but I guess the correct thing would be to have a "transmission" step between thresholding and synapses? Or after synapses? The question is: If a connection has a delay of 0, should a threshold crossing of a neuron lead to an effect in the pre-synaptic cell in the same timestep or rather the next one? Possibly we could also let the synaptic state variables update in the same slot as the NeuronGroup state variables and only execute pre/post codes in the synapses step?

[thesamovar]

We have to think about how other ways of creating synapses fit in as well, e.g. setting directly via a (sparse) matrix.

That's an interesting one. I think we could maybe make a case for NOT allowing setting via a sparse matrix. Apparently the sparse module of scipy is scheduled for improvement this year, but for the moment it's still terribly inconsistent and unreliable, and it's one of the most annoying things about Brian 1 development coping with that.

That might be quite a bit of work to get right... We have to be very careful about the reference/value semantics here. For example, should v = G.v; v[:] = 0 overwrite the state variable values (I'd say no)? Do we want G.v += 1 to work (I'd say yes)?

Indeed, these are tricky issues. If we don't want the first one to work (and I can see your point!), I think it means we can't have G.v[:]=0 working either, because this is equivalent, and I think we do want that to work. Or at least, it would be surprising if it didn't do anything.

I think there's no need to treat slices vs. ints/arrays differently

Well my thinking here is that if you did, say, S.w[:,:,0] I would expect it to find all possible synapses. If you did S.w[arange(N), arange(N), 0] I would expect EITHER that all pairs (i, j) exist OR (following numpy fancy-indexing conventions) that all pairs (i, i) exist.

What I'd suggest for all syntax questions (not only Synapses, but also refractoriness, etc.): Let's put up a wiki page with concrete examples of how it does/should look like

Agreed, discussing it only in the issue (issue 42 by the way! :)) can be difficult for someone else to follow.

[thesamovar]

BTW: Will we ever need a dynamical array with more than one dimension in our current design where we no longer have the state matrix _S? If not, we could merge DynamicArray and DynamicArray1D and probably simplify the code quite a bit.

I think SpikeQueue uses a 2D dynamic array.

[thesamovar]

The question is: If a connection has a delay of 0, should a threshold crossing of a neuron lead to an effect in the pre-synaptic cell in the same timestep or rather the next one?

That probably needs some careful thought with examples. In any case, I'm happy to make a new 'transmission' slot if needed.

Making SpikeQueue a BrianObject I'm not sure about. I think of it rather as an entirely self contained utility class, while the Synapses object which manages actually getting spikes and passing them to SpikeQueue, etc. Why do you think it should be a BrianObject? (I'm willing to be convinced.)

[mstimberg]

Ok, forget about the DynamicArray question, you are of course right that SpikeQueue uses a 2D array...

Making SpikeQueue a BrianObject I'm not sure about. I think of it rather as an entirely self contained utility class, while the Synapses object which manages actually getting spikes and passing them to SpikeQueue, etc. Why do you think it should be a BrianObject? (I'm willing to be convinced.)

I didn't think much about it when I implemented it this way, it was just the easiest thing to do (I just converted SpikeQueue's propagate method into an update method). But I think you are right, modeling it as BrianObject is unnecessarily complicated, pushing/retrieving spikes should happen from outside. Now that I'm thinking about it, the code object that executes the pre/post code (currently called TargetUpdater, that's not a good name) would maybe be a more logical choice to contain it and take care of the push/retrieve. It could possibly also contain the delays array (but we'd have to still make it available in Synapses for allowing to set S.delay[:, :] = ... ). If the object contains everything that is needed for the execution of the pre/post code, including the spike queue, then actually the multiple pathway thing could be quite easily implemented, as it would just mean adding another of these objects. From an implementation point of view, it would be clearer to not write S.delay_pre[:, :] = ... and S.delay_post[:, :] = ... but rather S.pre.delay[:, :] = which could then be extended to a pre2 pathway. On the other hand, we don't want to make the syntax for the simple and common case more complex just to support the rarely used multiple pathways formulation. But we could just introduce S.delay as a shortcut for S.pre.delay. Having to write S.post.delay instead of S.delay_post would be fine in my view, this is quite rarely used, anyway (actually, none of the examples in examples/synapses uses it).

[thesamovar]

+1 - awesome idea.

[mstimberg]

Getting back to syntax questions (I already pushed the changes moving the spike queue and the delays into the updater objects):

Well my thinking here is that if you did, say, S.w[:,:,0] I would expect it to find all possible synapses. If you did S.w[arange(N), arange(N), 0] I would expect EITHER that all pairs (i, j) exist OR (following numpy fancy-indexing conventions) that all pairs (i, i) exist.

I'm not sure... Having a simple rule that is easy to explain would be great of course. The simplest rule I'm seeing is using: Interpret S.var[I, J, K] (the indices can be integer values, slices or arrays) as: "Get all synapses (i, j, k), where i in I and j in J and k in K", slices would be interpreted as arange(start, stop, step). I would then only raise an IndexError if this results in an empty list.

Any other approach means that we have to treat the third index differently from the other two (which was my whole point from the beginning....). But I have to admit, that it would be nice if we could make this 3-index notation work because in many simple cases it is quite clear.

I think the following might work: Interpret S.w[i, j, k] as

for synapses in S.w[i, j]:
   yield synapses[k]

This would mean that we use all the standard indexing for the first two dimensions (and I guess then we should also apply numpy's fancy indexing semantics) only apply the third indices if the synaptic connection exists. This is still quite clear to explain and it covers the two cases you mentioned in the way you expect them. An expression like S.w[:, :, 1] would raise an IndexError if one of the connections has only one connection, which also makes sense. What do you think of this approach?

It would be a nice advantage over the current system in brian 1, where you are only allowed to use integers for i and j if you use an index k. E.g. writing:

S.delay[:, :, 0] = 1*ms
S.delay[:, :, 1] = 2*ms

wasn't possible.

Oh, and I think the returned array (or the "array" used for setting) should be always a flat array and we should use your S.i, S.j, S.k idea. This would mean if one has a function weight(i, j, k) returning the weight for a synapse (and correctly working for index arrays). The following would be equivalent:

S.w[:, :, :] = 'weight(i, j, k)'
S.w[:, :, :] = weight(S.i[:, :, :], S.j[:, :, :], S.k[:, :, :])

That's not the way that string assignments are evaluated currently but maybe that would be a good way, actually?

[thesamovar]

I wouldn't like to raise an IndexError only if there are no matches. It might lead to unpredictable bugs in some corner cases (e.g. someone is looping through sparseness in linspace(0, 1, 100) and then they would have to special case sparseness=0 to avoid the IndexError. On the other hand, I would also want S.w[1,2,3] to raise an IndexError if that synapse doesn't exist.

I'm still not quite sure what the problem is with using : as a 'matching' operator, and arrays and ints as 'non matching' as in numpy. For example if x=array[1,2,3] then x[2:5]==array([3]) but x[arange(2,5)] raises as IndexError. So we have precedent in numpy for considering : as matching anything that exists, and for raising an IndexError when doing explicit indices that don't exist.

The interpretation of S.var[I, J, K] would then be that an existing synapse (i, j, k) matches if either (a) there are no arrays, and i=I or i in slice I; j=J or j in slice J; k=K or k in slice K; (b) there are arrays, and there is an index l say so that i=I, i in slice I, or i=I[l]; j=J, j in slice J, or j=J[l]; k=K, k in slice k, or k=K[l]. If none of I, J or K is : then we expect a result of length 1 (if they are all ints) or the length of the array if at least one of them is an array. If the length of the result doesn't match that, then we raise an IndexError.

That reads as complicated mathematically, but I think it's how numpy fancy indexing works and covers all the cases we need (including treating k in the same way as i and j). It's predictable in the sense that no error will ever be raised if one of your indices is : and an error will always be raised if none of your indices is : and one of the synapses you specify is missing. I think it's also reasonably intuitive if you know how numpy fancy indexing works, and even if you don't it's intuitive in the easy cases.

So for example S.w[:,:,1] would match all synapses with k==1 (which might be the empty set). And, we could write S.delay[:, :, 0]=1*ms.

Agreed on the last point - we should use that same mechanism for evaluating the string assignments.

[mstimberg]

I'm still not quite sure what the problem is with using : as a 'matching' operator, and arrays and ints as 'non matching' as in numpy. For example if x=array[1,2,3] then x[2:5]==array([3]) but x[arange(2,5)] raises as IndexError. So we have precedent in numpy for considering : as matching anything that exists, and for raising an IndexError when doing explicit indices that don't exist.

This is what my proposal (the last one) would do, no? You are proposing that having an explicit index for k would not raise an error, even if it doesn't exist for some i and j.

It's predictable in the sense that no error will ever be raised if one of your indices is : and an error will always be raised if none of your indices is : and one of the synapses you specify is missing.

I think my point is that this is somewhat different from numpy syntax. You can get index errors even when using : if one of your indices is out of bounds. So I think your proposal would be to ignore out of bounds for the last but not for the first two indices? What about S.w[0, :, :] if source neuron 0 does not have any connections, according to your proposal you would return an empty set here, right?

I guess it all boils down to the question, do we want S.[:, :, k] always to succeed? The difference between our proposals isn't very big in the end, I see two differences:

1. My proposal would raise an index error, where yours would return an empty set. If we deal with assignments to empty sets in the same way as numpy, this means that we would silently ignore incorrect assignments, e.g. doing S.var[:, :, 1] = ... if all connections have only one synapse would just do nothing (according to mine it would also ignore it if you do S.var[:, :, 1:], though).
2. In cases where you have different number of synapses for neurons (where different means different numbers > 0), you would allow S.var[:, :, 1] = ... to set all the 2nd synapses, whereas I would raise an error if some connections have only 1 synapse.

I'm not sure what is more intuitive from a user's point of view...

[thesamovar]

I'm not sure what I think would be more intuitive for the user either. Actually, since we're proposing major changes to Synapses we'll have to have a meeting to discuss it anyway, let's add this to the list of options to discuss?

Here's my summary of the issues so far:

New features:

• Adding/pruning synapses should be possible even if not efficient.
• Multiple data structures

Syntax changes:

• S.num_synapses[i,j]=1 instead of S[i,j]=True (more explicit)
• Add S.i, S.j, S.k and use them in string assignments
• Multiple pathways in multiple Synapse objects or multiple TargetUpdater objects (simpler implementation)
• Remove float syntax for probabilistic synapses and use strings (more explicit)
• Lumped variable syntax change (more explicit)
• Array access syntax, several options to consider

Implementation:

• What sort of object should S.w and G.v be? Should v=G.v; v[:]=0 modify the group variables? This is equivalent to G.v[:]=0. Should G.v += 1 work?

The big one is the array access syntax and we need to set out all the possible cases and how the ideas differ.

I think my point is that this is somewhat different from numpy syntax. You can get index errors even when using : if one of your indices is out of bounds. So I think your proposal would be to ignore out of bounds for the last but not for the first two indices?

Ah I see. OK, so it looks like it's impossible to get something that is consistent with numpy but also generalises to this case. Right? If so, we should focus on doing something that is (a) not too surprising, (b) easy to explain, (c) functional. In that case, I would say that my proposal is entirely symmetric in i, j and k but doesn't match numpy in that with slices you can never raise an IndexError. Some cases:

• S.w[:,:]=0. (Side note: do we assume that the missing third index is a :?) Is it more surprising if it raises an IndexError if there are no synapses or not? Personally, I feel like I would never expect this statement to fail because the intent of the user seems clear here: do it for everything, no matter what everything is.
• S.w[:, :, 1]=0, does it set all second synapses to 0 or raise an error if not all synapses have a second synapse? I can't decide which is more surprising here, but I feel like the latter is less functional in the sense that: how do you set all second synapses to 0 if not all synapses have a second synapse? I guess the point here is that raising an IndexError always reduced functionality, so we should only do it if it is justified because we think the user has made an error in their code.

The last point suggests an idea for a new notation: S.w['cond']=values. For example, S.w['i>j']=1'. IfS.w[p]means the pth synapse, then this would be equivalent toS.w[S.i>S.j]=1`. If this notation was available, then maybe it would be OK to have a more restrictive slicing syntax.

OK I'm tired maybe I should stop writing now. ;)

[mstimberg]

Thanks for the nice summary. I started listing syntax questions in the wiki (nothing for Synapses yet), I used restructured text because I thought we might want to copy some parts to the docs at a later stage. Unfortunately that has the disadvantage of no easy way to link other pages, one has to use the full URL for that.

OK, so it looks like it's impossible to get something that is consistent with numpy but also generalises to this case. Right?

Yes, I think we cannot treat it 100% like a numpy 3d array since we have the complication of a "holes" (no synapse for an (i, j) pair) and a variable third dimension. I tend to think that your approach for the S.w[:, :, 1] example is better, we will have to make this very clear in the documentation, though. The only thing I'm worried about is the case where no connection has a second synapse, here we would silently fail (setting an empty array). This might be a bit unintuitive since using an explicit integer or array in numpy will never fail silently, it will either raise an error or set something. Maybe we should have a special case for the k index here: If it is explicit (int or array), it has to match at least once (or once for every entry in the case of an array). This means that S.w[:, :, 1] would raise an index error if no connection has more than one synapse. S.w[:, :] or S.w[:, :, :] would not raise an index error, even if no synapse exists at all but S.w[:, :, 0] would. Hmm, but then we should raise an error for S.w[0, :, :] or S.w[:, 0, :] as well, I guess... I think raising error in these cases would be in line with your statement that we don't want to restrict functionality but only raise errors when we think the user made a mistake. Actually, we could also emit a warning if a slicing operation does result in an empty set, e.g. S.w[:, :] without any synapses.

I guess my preferred choice would now be something that is close to your proposal but instead of "never raise an error if one index uses a slice" I would raise errors in the numpy way with an exception for the k index where an error is only raised if the index (or every of the indices in an array) never matched. I'll try to formulate this in a more formal way. I think we now have something like 4 feasible approaches that are all of comparable complexity from the implementation point of view. We could list them all and then do a table with concrete examples, including all the corner cases we discussed here and show what the different approaches would do (e.g. raise an error, set all second synapses, etc.). Then we should either schedule a meeting with Romain or let him look over the examples. And then we can have a majority vote :) Or we ask others for feedback as well (Victor, Achilleas). We should come to a decision soon, though...

I like the idea of a string index! Actually, that would be even a great addition for synapse creation. I think writing:

S.num_synapses['i!=j'] = 'rand() < 0.1'

is considerably more readable than

S.num_synapses[:, :] = '(i !=j) * (rand() < 0.1)'

[thesamovar]

OK, sounds good. I think we should get as much feedback as possible on this one - it's pretty critical.

[mstimberg]

I added examples for the different indexing options we discussed to https://github.com/brian-team/brian2/wiki/Syntax-Synapses#state-variable-access -- your proposal from above is "Option 3", please check whether I'm not misrepresenting it.

[thesamovar]

Looks good. I added a few comments and extra bits.

[mstimberg]

Looks good. I added a few comments and extra bits.

I replied to your comments in the wiki. While writing I realized that there's yet another possiblity that is between our preferred options (3 and 4): Raise errors / return empty arrays exactly as numpy would do by assuming the size of the third dimension being the maximum number of synapses after applying the first two indices. This wouldn't raise an error anymore for S.w[0, 0, 1:] in line with your option, but it would keep raising errors for S.w[:, :, [0, 1]] if no connection has more than one synapse. Of all the proposals, I think it's the closest to numpy semantics and it's easy to explain. I think I'll switch "my option" to this new one.

[thesamovar]

I still feel like anyone writing code with variable numbers of synapses is going to have to wrap a lot of stuff in try/excepts like this. Suppose you know that for each pair (i, j) you generate between 0 and 4 synapses according to some random distribution, and then you want to do something for all these synapses. If you wrote code looping through the synapses for example for i in xrange(4): S.w[:,:,i]=... but by chance sometimes there was no synapse number 3 then it would raise an IndexError but unpredictably. OK, it's a bit contrived as an example, but I just don't like the idea of exceptions being unpredictable like that.

[mstimberg]

I still feel like anyone writing code with variable numbers of synapses is going to have to wrap a lot of stuff in try/excepts like this.

I think in general we spend 90% of our time here discussing corner cases that are very rare anyway -- using a varying number of synapses is already something that almost no one uses.

Suppose you know that for each pair (i, j) you generate between 0 and 4 synapses according to some random distribution, and then you want to do something for all these synapses. If you wrote code looping through the synapses for example for i in xrange(4): S.w[:,:,i]=... but by chance sometimes there was no synapse number 3 then it would raise an IndexError but unpredictably.

Well, that is a bit contrived indeed :) What about having code that creates 3d-arrays with varying third dimension, wouldn't the same reasoning apply there? What I like about "my solution" is that it is very close to numpy for the "0 or 1 synapses" case, you'll get exceptions whenever numpy would give you an exception for the corresponding 3d array. I think if I write S.w[:, :, 1] = ... and no connection has a second synapse then it is very likely that this is a genuine mistake. Actually, I'm still allowing for your random synapses use case (even though this is probably not obvious enough for the average user): You can write for i in xrange(4): S.w[:,:,i:i+1]=... and it will never raise an error. I'll go ahead and update the option 4, it's time for getting outsider feedback, I guess :)

[thesamovar]

I agree that we're mostly discussing corner cases, but that's OK. I guess my list of priorities is something like:

1. It's difficult to write code that silently performs wrongly (i.e. raise errors if there is any doubt).
2. Everything the user wants to do should be possible.
3. The syntax should be easy to explain.
4. The syntax should match numpy to the extent that this is possible.
5. We should optimise for the case of a single synapse.
6. Errors should be predictable and not depend on parameters, i.e. if I run it on a small example and it doesn't raise any errors, it ought to run on a big example. I think this is important because it raises frustration levels with the software if you have to debug code on a huge simulation that might take days to complete rather than on mini examples.

So I put the unpredictable errors as the lowest priority, but I still think it's important because it can really be frustrating when you get hit by unexpected errors. I know because it happens even in my Brian code quite regularly and then I need to rerun a big simulation.

Another option we haven't really discussed, but might be worth thinking about: we have a more radical change, dropping the numpy-like syntax entirely. This minimises the chance of confusion, but it would also be a shame to lose it.

[mstimberg]

I agree with this list and would of course argue that "option 4" scores best on these points :smiley: A particular advantage is that a fully-connected, same-number-of-synapses network will show exactly the same behaviour as numpy. And I don't think that number 6 counts against it, I don't see where scaling-up makes a difference.

[thesamovar]

Agreeing on the list of priorities is a good start! :) Perhaps we should move this list into the wiki too?

For (6) my point is just that what seems like a very unlikely corner case can actually become more likely when you're doing things like big parameter searches where unlikely cases will very likely turn up occasionally.

[mstimberg]

I moved the list to the wiki.

[thesamovar]

An implementation issue that I just thought of. If we want to allow people to do stuff like S.num_synapses['cond(i,j,k)']=1 we probably don't want to construct all possible i, j and k in one vast triple of arrays and then evaluate that because there wouldn't be enough memory. For setting weights there ought to be enough memory in most cases, but it's still wasteful to allocate these huge arrays for i, j, and k and then throw them away. I suggest instead that we will want to buffer this in some way. In Brian 1 we loop through all i and vectorise over j and k essentially. I think this is reasonable but we might be able to improve on it. Not sure if we can improve on it much for synapse creation, but for setting values for existing synapses we will have access already to the i, j, k arrays so we can just choose a fixed buffer size and go through it in chunks of that size. This would be more efficient when there are not many j and k for each i.

[mstimberg]

Agreed. And of course it's the same issue for S.num_synapses['cond(i,j,k)']=1 and S.num_synapses[:, :, :] = 'cond(i, j, k)

[mstimberg]

Two more things I just realized:

1. We don't have three indices for synapse creation, k doesn't make sense here
2. The problem isn't actually generating all possible combinations but only the evaluation -- the rest numpy can do by broadcasting. E.g. while I, J = numpy.meshgrid(np.arange(10000), np.arange(10000)) would return two 10000x10000 arrays, I, J = numpy.meshgrid(np.arange(10000), np.arange(10000), sparse=True) would only return a 10000x1 and a 1x10000 array -- in both cases cond(I, J) would evaluate the same.

[mstimberg]

Coming back to fundamental implementation issues: In the Synapses class, the connection information is stored in a duplicate fashion: We have two arrays stating the pre/postsynaptic target for every synapse and we have two lists of synapse indices for every pre/postsynaptic cell. Both can be calculated from the other but the update code uses the first structure while the SpikeQueue uses the second. I'm wondering whether we can't use the same structure for both, probably rather the former (flat dynamic array) instead of the latter (list of dynamic arrays).

I don't see any principle reason why this shouldn't work but maybe it would be less efficient (or just horribly to code)?

[thesamovar]

Good points. For the implementation, let's give it a go and see how efficient it is.

[mstimberg]

I had a quick look at the redundant information storage in synapses and, well, it's there for a good reason. If a pre-synaptic neuron i spikes it allows to find the relevant synapses via synapses[i], i.e. with constant cost. From the other datastructure (pre-synaptic number for every synapse) one would have to do something like (synapses == i).nonzero()[0], i.e. scaling with the number of synapses.

We could make this more efficient by having the synapses sorted and then using binary search (logarithmic in number of synapses), or making it even more efficient (constant lookup) by storing an array with the start indices for every neuron.

The only problem is that this only works for the presynaptic indices, in the case of a postcode we would have to either resort to the slow synapses == i expression or store a datastructure in the same way it is currently done. But at least for now I think that having some wasted memory is preferrable over having two different code paths for pre- and postsynaptic spikes. We can always optimise this later.

The only thing I'm considering: The datastructure is only used within SpikeQueue and relatively simple to construct from the other structure. Maybe we should construct it from scratch every time in the SpikeQueue.compress method? That would make things in Synapses a bit simpler, because we have not to bother with two structures during synapse creation.

[thesamovar]

For constructing from scratch it's not a bad idea but as long as we don't reconstruct when not necessary. A fairly quick way to test if its necessary is to check the sizes of the array and do a hash of the index arrays - if these are unchanged we don't need to reconstruct. (Or we can keep a modified flag but that's potentially more bug-prone.)

So let's proceed with the basic data structure as it is, but keep in mind that there are possible optimisations for the future. Seems good?

[mstimberg]

Agreed on all points. Currently, the checking is simple since we do not allow to delete synapses, so all that could happen is that some were added.

[mstimberg]

I worked on Synapses and made quite a bit progress (of course it's missing a lot of documentation, tests, etc...): The state variable indexing works (I didn't bother with raising any indexing errors for now) for most of 1-d, 2-d and 3-d indexing and also for string indices. It's all quite straightforward and I did not bother with optimizing it yet (I didn't take exactly the same approach as in the old Synapses class), it seems to be fast enough for now at least. I refactored all the indexing into its own class and also refactored SpikeQueue a bit: It is now a "dumb object", in brian1 there was a bit of inconsistency between the C version and the Python version, the C version had to be wrapped in a Python class that implemented the compress and push method -- this is no longer necessary. Accessing state variables and delays seems to work as expected, I also implemented the new multiple pathway syntax.

The big chunk that is missing is the synapse creation (there's a simple connect method to which you can provide an array of pre- and an array of postsynaptic neuron indices, that's what I use for testing right now) and making everything work with the new Devices framework.

[thesamovar]

Great stuff! I'll try to take a look at it later today. So can we already run Synapses code in both Python and C++ for example? If so, it's already an improvement on Brian 1 in that respect (we never did get propagation code running in C++).

[mstimberg]

So can we already run Synapses code in both Python and C++ for example?

Yes we can! And with the named pathways we even have a second feature that Brian 1 didn't have. Oh, and I also allowed for multiple "post" pathways (3rd feature!) -- I'm not sure there is any real use for it but someone might come up with a fancy STDP rule that needs it. And to be honest, it meant less code to support it because we now treat pre and post the same.

[thesamovar]

OK I just had a look. Nice work indeed! Nothing to add for the moment.