Efficient synapse creation

[thesamovar]

OK a thought on how to efficiently implement synapse creation. The problem to solve is that if you just compute cond(i,j) for all pairs (i, j) it is O(N^2) and in many cases this is highly inefficient. We should list some of the cases and propose solutions.

The two most important that I can think of are:

• rand()<p should be replaced with binomial distribution sampling, and also compounds like rand()<p and cond(i,j) should be replaced by binomial distribution sampling followed by condition evaluation only on that subset.
• i==j should be replaced with one-to-one

So there are various possible solutions here - it would be nice to come up with something quite general that can adapt and become more efficient as we add special cases, but that doesn't require us to rewrite all the templates each time.

The most obvious implementation of i==j is just to explicitly detect it and convert it to connect_one_to_one and require that devices implement this connectivity function. Similarly for rand()<p and connect_random. I don't really like this approach because it's not very general.

Another option is to have a set of flags that are passed to codegen and the template can make use of that extra information. For example, if the condition was rand()<p and f(i,j) we could pass to codegen sparseness=p and cond=f(i,j). The code would then look something like:

for i in xrange(_num_source_neurons):
  if sparseness<1:
    j = binomial_sampling_code_here
  else:
    j = arange(_num_target_neurons)
  _cond = f(i, j)
  ...

A more general and flexible mechanism would be to handle all conditions of the form a and b and c, parse them into a set [a, b, c], decode this set a bit and then pass that to codegen. For example, a might be rand()<p, b might be abs(i-j)<=width and c might be i!=j. This would be parsed into:

a: ('rand()<p', 'binomial', p)
b: ('abs(i-j)<=width', 'abswidth', width)
c: ('i!=j', None, None)

The parsed form is (original_expr, specialisation_name, argument). The codegen template could then look through this list, check which ones it supported, implement them directly if they are supported and if not fall back on the default mechanism.

It might be the case though that it's only worth doing one specialisation, because you won't improve the order of the computation by doing multiple specialisations assuming that each one returns an array of candidate indices. In fact, you might even slow it down. So another option is that each specialisation has an estimate of the number of candidate synapses j it will generate for any given i, and we simply choose the best one. For example if there was rand()<0.5 and i==j then rand()<0.5 has an estimated 0.5*_num_target_neurons for each i, whereas i==j has precisely 1, and in most cases 1<0.5*_num_target_neurons so this will be selected.

Any thoughts on the above?

[mstimberg]

Ok, I have to put some thought into this. Might also be interesting to look at the connection-set algebra implementation for inspiration. Will be back tomorrow with some feedback.

[mstimberg]

Ok, no real feedback yet... But I'm thinking that we should maybe restrict efficient approaches to correct use of the connect function, where you already have probabilities and condition separated. We wouldn't have to deal with taking string conditions apart, those would be just executed as they are. We could still think about optimisations for expressions involving i==j or subsets like i > 10000 for example, but this would be more a long-term plan, I guess. First get something running, then worry about efficiency :)

[thesamovar]

I agree we don't need to prioritise this - I marked it for the beta release rather than the alpha/dev-preview release.

I'd rather everything worked with the string conditions, and that actually things like the p=... argument for connect would only work if you hadn't declared any of the other arguments. Otherwise the semantics are not entirely obvious. (I think there's only one consistent way to interpret it, which is to and it with the other conditions, but it's not entirely obvious to do that). I don't think analysing the conditions is too difficult to do. You simply check the ast representation, is it an and structure?, if so it gives you a list of the expressions that are being anded, you check each one if it is of the required form, etc.

[mstimberg]

I'd rather everything worked with the string conditions, and that actually things like the p=... argument for connect would only work if you hadn't declared any of the other arguments.

I think in our e-mail discussion with Romain the idea (his idea at least ;) ) was to make all the options combine, e.g. instead of writing S.connect('(i != j) and (rand() < 0.1)') one would rather write S.connect('(i != j)', p=0.1)

You simply check the ast representation, is it an and structure?, if so it gives you a list of the expressions that are being anded, you check each one if it is of the required form, etc.

Separating it into and/or components seems easy enough, I'm more concerned with analysing the form, e.g. would we recognize not only 'i == j' but also '(i-100) == j', etc.? I'm a bit worried that in the end it will not be entirely clear for a user under what circumstances expressions will be treated efficiently.

Anyway, let's discuss this in more detail post-Alpha :)

[thesamovar]

Well we didn't really discuss that idea much. What do you think about it? I feel it's less explicit and doesn't save much typing.

For analysing the form, I think it's ok if we incrementally improve the recognised forms over time. The only critical ones are rand()<p and i==j I would say, but we could try to be more general over time.

But yes, let's leave it for the moment and just go with the O(N^2) technique for now.

[mstimberg]

Well we didn't really discuss that idea much. What do you think about it? I feel it's less explicit and doesn't save much typing.

I kind of like the idea, it separates strict conditions (e.g. "no self-connections") from probabilities (e.g. "Gaussian connection probability depending on distance"), joining all together with and makes this less obvious, I think.

[thesamovar]

OK, if you both prefer it I don't mind being overruled on this!

[thesamovar]

An even simpler suggestion. We allow the user to specify a function range(i) which returns a pair (low, high) and only values of low<=j<high are checked. I think I mildly prefer my general suggestion above (which would include this as a specific optimisation), but it's an alternative to consider.

[mstimberg]

I'm not sure I 100% understand this proposal -- how does this connect to your general suggestion of parsing and analyzing the string?

Just to copy here my remark I made in an email a while ago so that it doesn't get lost:
I guess I would be more in favour of some explicit method, e.g. something like a connect_each_source_to method which would allow you to write something like:

# One-to-one-connectivity, loops over all i
S.connect_each_source_to('i')
# Connect each cell to its two neighbours, "counter" runs automatically
# from 0 to n-1
S.connect_each_source_to('i + (-1)**counter', n=2)

[thesamovar]

It doesn't connect to it - it's an alternative. The advantage of specifying a range of values to consider is that you can do combined operations like connect i to j if i-W<=j<=i+W AND rand()<0.1 and get the speed advantage of only needing to evaluate rand() 2W+1 times.

[mstimberg]

It doesn't connect to it - it's an alternative.

Ok, I just wondered because you said "I mildly prefer my general suggestion above (which would include this as a specific optimisation)" above. Anyway, this would be a big addition that really needs some careful thought.

[thesamovar]

Indeed - let's think about it another time!

[thesamovar]

I just had a few thoughts on this.

Firstly, we could make a potentially enormous optimisation that is quite general. If the set of indices used in an expression are only i or only j (and there are no stateful functions like rand and randn) then we can evaluate in two passes by computing the set of possible values of i and j first and then only looping over those. So for example if someone had the condition i<1000 and rand()<p then we note that the subexpression i<1000 has to be true for the whole condition to be true (it is a sequence of and statements) and so in phase 1 we loop over all candidates i and compute the set where i<1000 before evaluating rand()<p on all pairs i, j (where i is drawn from the subset we computed). This works for arbitrary expressions involving only i or j (or variables indexed by them) so we don't need to pattern match specific cases like i<K.

As a followup to that, we can also look for expressions of the form i==... (or j=...) where the RHS only depends on indices j (or i respectively). Can sympy be used to do that? I think so probably. This would allow us to automatically handle many cases like i==j or i==j-10. We'd just compute the RHS for all the j values, and then check which of those values were valid i values (or vice versa).

Can we do something that is similarly general that copes with expressions like abs(i-j)<10 and i<=j? I don't immediately see how but maybe some more things are possible? Obviously we can't deal with the fully general case without computing all pairs i and j, but maybe we can deal with enough of them?

[thesamovar]

I would love to have this in 2.0 but maybe it's too much to do on short notice. I think maybe the suggestions in my previous post are workable and could be done reasonably quickly. Thoughts?

To add to the suggestion in my last paragraph, for dealing with things like abs(i-j)<10 we could imagine that this could be handled relatively simply in the following way:

1. We note that the expression i-j only has a limited range of values because i and j are integers, so it can only vary between min(i)-max(j) and max(i)-min(j).
2. Now we can iterate through that range and check if the conditions holds (this only works if the RHS of the expression doesn't mention i and j of course.
3. Now we have all the values of i-j where the expression is true, and we can reconstruct all the pairs (i,j) that lead to this value of i-j quite easily.

Steps (1) and (2) hold for arbitrary integer subexpressions, like i+j or i-2*j, but i haven't thought about whether or not step (3) holds in this way or not. I think maybe it does. In any case, I think if we implemented some of this, a lot of typical use cases would be very efficient compared to how they are at the moment.

[thesamovar]

The linear case a*i+b*j can at least be handled, as it's a single linear diophantine equation in two variables: http://mathworld.wolfram.com/DiophantineEquation.html

[thesamovar]

I just thought of another potentially easier method. It's not as automatic or as general, but it would be easy to implement. As well as the user giving a condition cond(i,j) they can optionally give extra functions jmin(i), jmax(i) and we only evaluate cond(i,j) for jmin(i)<=j<=jmax(i). Or maybe even more flexible, they can give a jrange(i)=(start, stop, step). This puts the burden on the user to compute the boundaries, but then covers most typical use cases.

[thesamovar]

Update from CNS meeting. We didn't agree on a syntax or a structure but we did agree that it would be good to have a way of providing code (in python/c++/etc. that would generate the synaptic connectivity). In Python this could be a generator function (with a yield statement), but with C++/Cython it would be a loop.

[thesamovar]

I've been thinking about an implementation for some of these ideas and I think I've got something that will be doable in time for 2.0, isn't too complicated and allows it to be efficient (with a bit of user cooperation).

First though, a quick review of some of the ideas:

The general idea is to break down the condition into a disjunction cond_1 and cond_2 and cond_3 and analyse the form of each condition to produce a smaller set of i, j pairs to iterate over. The simplest of these include:

• j<=f(i) and similar means we have a reduced set j to iterate over for each value of i, particularly if we have j<=f(i) and j>=g(i).
• rand()<p(i) means we can generate candidate pairs by generating a binomial random number and shuffle.

So this leads to the strategy I think would be doable for 2.0.

Analyse the conditions. For each i, compute jmin(i) and jmax(i) where possible. If there is no rand()<p(i) term then iterate from jmin to jmax checking the condition. If there is a rand()<p(i) term generate the set of j between jmin and jmax and iterate over those, checking the full condition minus the rand()<p(i) part.

This covers a lot of cases, like i==j would be efficient, as would rand()<p. One thing that is missing is abs(i-j)<K would still be inefficient but if the user would rewrite it as -K<=i-j and i-j<=K then it would be efficient. We could actually think of handling just this very specific case of abs(i-j) explicitly since it's probably quite common. The point about it is that it allows the user to do things efficiently, whereas the current Synapses mechanism makes it impossible to be efficient. That's why I think it's important that something be in 2.0.

A big thing that is missed by this but I don't think there's any way to solve the problem is rand()<f(i,j) (e.g. f(i,j) is some Gaussian distribution). The only way we can make this efficient is for the user to additionally specify limits for j.

Other things that would be nice for a more general system (post 2.0):

• cond(f(a*i+b*j)) where cond is a condition and a and b are integers. This can be solved using the diophantine equation trick as described above at cost O(N log N) + however long it takes to iterate over all the solutions (which could be O(N^2) but only if there are O(N^2) solutions).
• f(i)<g(j) can be computed quite efficiently by computing f(i) once for all i, and g(j) once for all j, then getting a sorted map from g(j)->j. Now for each i we can find the set of j satisfying the condition by doing a bisection on this sorted map and then iterating over the corresponding j values.

The way I anticipate implementing this is that we will just do some computations and pass additional variables to the template. This will be backwards compatible because templates are free to ignore this information and just use the condition directly. In addition, when we add the more detailed information in the future we can just further add additional variables and keep backwards compatibility again.

[thesamovar]

A user could somewhat efficiently implement rand()<f(i,j) if they computed analytically g(j)=max(f(i,j) for all j) and then changed the condition to rand()<g(j) and rand()<f(i,j)/g(j). The first of these two conditions would efficiently generate a set of candidates i and then you'd evaluate the second condition for each of these candidates, and the total probability would be the same (independent random variables). This would allow efficiency in the case that you had a good bound on f(i,j). If you have something like a Gaussian probability distribution around a centre value, you could do it in three connect calls, one that has the condition say i>j-sigma and i<j+sigma and rand()<exp(-(i-j)**2/(2*sigma**2)) and the second has the condition i<j-sigma and rand()<f(sigma) and rand()<f(i-j)/f(sigma) where f(x) is the Gaussian pdf (and similarly for the third with i>j+sigma.

Again, this isn't quite as nice as it could be, but the point is that it allows users to make it efficient if we write a little guide to how it works.

[mstimberg]

We did not quite follow the approach described here, but I think for all practical purposes everything mentioned here has been dealt with in #600. If something remains to do, let's open a new more specific issue.