Threshold conditions that depend on subexpressions fail

[brian-team]

Threshold conditions that depend on subexpressions fail

[tnowotny]

What type of content is allowed in a subexpression? Will it refer to variables? Then option 3 would necessitate creating and storing additional arrays for the evaluated subexpression?

I wonder whether not
4. Push the sub-expression into a support code function would be easy and efficient enough?

But I may not fully appreciate the idea and range of the subexpression concept.

[mstimberg]

Subexpressions can refer to other subexpressions or to other variables, the code generation system normally deals with all this automatically.

Then option 3 would necessitate creating and storing additional arrays for the evaluated subexpression?

Yes, the constant over dt flag is basically a shorthand, if you write:

group = NeuronGroup(..., '''subexpr = ... : volt (constant over dt)''', ...)

then this is equivalent to

group = NeuronGroup(..., '''subexpr: volt''')
group.run_regularly('subexpr = ...', when='before_start')

So for brian2genn, subexpr would be handled just like any other array.

4. Push the sub-expression into a support code function would be easy and efficient enough?

Converting a sub-expression into a function would mean identifying its dependencies and rewriting the expressions where it is called (e.g. in my example above, rate would turn into rate(t) when used). This would be nice, but is far from trivial (I mentioned this approach in a Brian issue: brian-team/brian2#412, but the fact that subexpressions can refer to other subexpressions complicates things quite a bit...).

[tnowotny]

Ok, I can see that the function idea is quite ambitious and we might not want to spend too much time on it.
Short of this, option 2 would probably be the most performant in the GPU context but carries some of the same problems. If sub-expressions include each other, substituting everything into everything could become a bit of a job.
Short of doing that, I am thinking, again with my GPU hat on, option 1 is probably more efficient in most scenarios than 3 because most likely local variables would be in registers and no global memory access would be necessary. Whereas in 3, the subexpressions would be written to arrays in global memory and read from them - two additional costly memory accesses.
To do version 1, the cheapest way, even if it's looking a bit hacky, would be to add the early lines of threshold code to the neuron update code. This, however, only works if GENN_PREFERENCES::autoRefractory is false (which is our default position in brian2genn for compatibility to Brian 2 conventions I believe?)
(if GENN_PREFERENCES::autoRefractory is on, the spike condition is also tested before the neuron update code is executed)

[mstimberg]

Hmm, I'm not sure of the performance penalty. In most cases, constant over dt should be better for performance since it will be only evaluated once even if the subexpression is used in several places. But this question is a bit separate, in general the ideal thing would be to allow both types of subexpression as for non-brian2genn Brian code.
If I understand your proposed approach for (1), then this would be purely on the brian2genn level, right? This might indeed work, I think we do not use the autoRefractory feature. However, I think we might run into problem if the neuron update code also refers to the subexpression, but I'll check it, maybe our codegen code would actually be smart enough to figure this out.
Thinking about this, we could probably have a slightly cleaner and somewhat less hacky solution: add the full Brian threshold code to the neuron update code (which will assign to a boolean variable _cond) and then just use a fixed threshold condition code _cond for GeNN. This should work, doesn't it?

[tnowotny]

Yes, that should work indeed (autoRefractory off, obviously).
I am not sure what you meant with there would be a problem if the neuron update code also refers to the subexpression. If it's the same subexpression (not a name clash with a differently defined one) then adding the subexpressions to the beginning of the neuron update code in GeNN should work?
Or are you cencerned it would be added twice?

[mstimberg]

Adding it to the beginning would be incorrect, since the values of the subexpression used in the neuron update code are the values before the update while the values used in the threshold should be the values after the update. If we just combine Brian's generated C++ code for the updater and the threshold, then we'd end up with duplicate local variable declarations. So I think the correct solution would be to combine the abstract code of the neuron update and the thresholder, then Brian's codegen mechanism would handle everything automatically (re-evaluate the subexpressions if necessary, but not re-declare them). This is a bit different from the approach we generally use in brian2genn (where we already deal with the CodeObjects generated by Brian), but I think it shouldn't be too complicated. I'll look into this.

[mstimberg]

I had a go at this and managed to implement it from the brian2genn side, but then I realized this won't work: GeNN uses separate "namespace brackets" for the neuron update step and the thresholder so we can't use variables from the former in the latter...
Maybe adding support for the "constant over dt" mechanism would be a simpler solution for now (in simple cases like my example, users might be able to simply insert the expression directly into the threshold, anyway, so this is more a matter of clear error messages).

[tnowotny]

hum ... I don't think the namespaces are a problem. Both the neuron code and the thresholding happen in the same namespace. I know for a fact that it works for the oldspike variable that is defined early on in the neuron kernel and then used in the spike detection much later ... Are you sure you are having trouble with this?

[mstimberg]

The namespace itself is not a problem, but it introduces a new local variable scope.
Here's what my version of brian2genn produced for the merged stateupdater+threshold (in this example, the state update code is empty:

// calculate membrane potential
             {
                 using namespace neurongroup_neuron;
            bool _cond = lv > 1;
                }
             // namespace bracket closed
// test for and register a true spike
             {
                 using namespace neurongroup_neuron;
            if (_cond)  {
                    glbSpkneurongroup[glbSpkCntneurongroup[0]++] = n;
                    }
                }
             // namespace bracket closed

The variable _cond is not accessible in the second scope.

I know for a fact that it works for the oldspike variable that is defined early on in the neuron kernel and then used in the spike detection much later ...

I had a look at GeNN's generateKernels.cc and it seems that the "namespace brackets" are only added if there's support code. I did not find any GeNN example in the userprojects folder that uses support code, so it might actually be that the oldspike variable leads to an error in GeNN as well as soon as neuronModel.supportCode is not an empty string...?

[tnowotny]

argh ... we need to look into this. It is slightly annoying that one cannot do using namespace xxx; in a limited region without creating the unnamed local variable scope ... or can one? Sorry if I am clumsy with this.

[tnowotny]

Of course, one solution would be for GeNN to stop doing the
using namespace neurongroup_neuron;
and instead for all calls to support code, substitute
neurongroup_neuron::call_to_support_code
which is unfortunately not entirely trivial.

Yet, alternatively, I am actually now wondering whether there is any potential for conflict within the kernel that would necessitate the (essentially 3) local variable scopes or whether it may be safe to just combine it into one. I am now thinking why would there be as we are using the same support code namespace for all elements of the neuron model... this may be a very easy fix in GeNN.

[mstimberg]

It is slightly annoying that one cannot do using namespace xxx; in a limited region without creating the unnamed local variable scope ... or can one?

I don't think you can do this, I think the standard recommendation would be to use fully qualified names, or define your local variables outside the new scope.

But looking at the GeNN code again, couldn't we merge all the "using namespace" directives that refer to the same namespace (neurongroup_neuron in this case)? I.e. in the simple example from my initial comment, we currently (with my changes to merge stateupdater and thresholder code) have:

extern "C" __global__ void calcNeurons(double t)
{
    //[...]

    // neuron group neurongroup
    if (id < 128) {

        // only do this for existing neurons
        if (id < 100) {
            // pull neuron variables in a coalesced access
            int32_t li = dd_ineurongroup[id];
            uint64_t l_seed = dd__seedneurongroup[id];

            // test whether spike condition was fulfilled previously
            // calculate membrane potential
            {
                using namespace neurongroup_neuron;
                double rate = ((_brian_pow(sin((((2 * (3.14159)) * t) * 1) * (1.00000)), 2)) * 100) * (1.00000);
                bool _cond = _rand(l_seed) < (rate * DT);
            }
            // namespace bracket closed
            // test for and register a true spike
            {
                using namespace neurongroup_neuron;
                if (_cond) {
                    spkIdx = atomicAdd((unsigned int*)&spkCount, 1);
                    shSpk[spkIdx] = id;
                }
            }
            // namespace bracket closed
            dd_ineurongroup[id] = li;
            dd__seedneurongroup[id] = l_seed;
        }
        __syncthreads();
        //[...]
    }
}

Instead we could use:

extern "C" __global__ void calcNeurons(double t)
{
    //[...]

    // neuron group neurongroup
    if (id < 128) {
        // only do this for existing neurons
        if (id < 100) {
            // pull neuron variables in a coalesced access
            int32_t li = dd_ineurongroup[id];
            uint64_t l_seed = dd__seedneurongroup[id];

            using namespace neurongroup_neuron;
            // test whether spike condition was fulfilled previously
            // calculate membrane potential
            double rate = ((_brian_pow(sin((((2 * (3.14159)) * t) * 1) * (1.00000)), 2)) * 100) * (1.00000);
            bool _cond = _rand(l_seed) < (rate * DT);
            // namespace bracket closed
            // test for and register a true spike
            if (_cond) {
                spkIdx = atomicAdd((unsigned int*)&spkCount, 1);
                shSpk[spkIdx] = id;
            }
            // namespace bracket closed
            dd_ineurongroup[id] = li;
            dd__seedneurongroup[id] = l_seed;
        }
        __syncthreads();
        //[...]
    }
}

The using namespace doesn't necessarily need brackets here, it is simply valid until the end of the scope (given by the if (id < 100) here). As a bonus, this would make the generated code shorter and more readable :)

[tnowotny]

Yes, agreed. I am going to fix this in GeNN. I was thinking one bracket instead of three but you are right, all of these are already encapsulated by an if statement on the neuron ID range, so none of the brackets are necessary. It will all become a lot simpler. Will do this Monday first thing.