Allow the definition of scalar parameters/subexpressions in model equations

[mstimberg]

[Maybe I shouldn't add more and more things to the alpha milestone, but this one might still lead to some quite fundamental changes in codegen (though I hope not), so we should get it sorted out rather sooner than later]

Internally, Variable objects have a scalar attribute and I think we handle it correctly in various places. We do not expose this functionality in the equations yet, though. One use case for scalar variables is a stimulus parameter (e.g. the orientation of a bar) that is the same for a NeuronGroup but changes dynamically during the simulation. I think we agreed a while ago on the syntax to have a (scalar) flag.

One remaining question: do we want to allow this flag for all equations? I did not think this through yet, but my feeling is that allowing it for differential equations might cause some problems, i.e. I'd only allow it for parameters and subexpressions (where subexpressions are only allowed to refer to other scalar parameters/subexpressions).

Related to that and #174: if the state updaters were allowed to return additional variables (in particular constants and scalar subexpression) in addition to the abstract code, they'd be able to assist in generating more efficient code.

[thesamovar]

One remaining question: do we want to allow this flag for all equations?

I think not for differential equations as you suggest. It's something we can add later if we see a need for it since as far as user syntax is concerned nothing would change if we allow it. We could even leave out support for subexpressions to start with.

Related to that and #174: if the state updaters were allowed to return additional variables (in particular constants and scalar subexpression) in addition to the abstract code, they'd be able to assist in generating more efficient code.

This sounds potentially interesting, but I'm not quite clear about what exactly you have in mind. For the loop invariant optimisations, I was imagining more along the lines of having a generic optimisation function that applies to any block of abstract code, so not specific to state updaters.

I'd be happy to include this in the next alpha release but I wonder if it's essential or not? Let's think about the implementation and if it involves some fundamental changes maybe we should include it, otherwise happy to let it slide to beta.

For the initial implementation, we don't need to do any optimisations with the scalars, just allow for their use. I think this would be very little work indeed. The trickier thing is pulling scalars out of loops, etc. Even that shouldn't be too tricky, right?

[mstimberg]

I started working on this in the scalar_parameters branch (only for parameters, let's postpone everything else) and I got 95% of the way quite easily. One trick I used that is slightly hackish but means that codegen doesn't have to change is to use '0' as an index for scalar variables. There's only one remaining problem (exemplified by the one failing test): conditional updates, e.g. the reset or pre/post statements. I guess one valid use case for a scalar parameter is the following:

G = NeuronGroup(N, '''dv/dt = 100*Hz : 1
                      n_spikes : 1 (scalar)''',
                threshold='v>1',
                reset='v=0; n_spikes +=1')

where in the end, G.n_spikes should hold the total number of spikes. This works in generated C++ code, the code is

for(int _index_spikes=0; _index_spikes<_num_spikes; _index_spikes++)
        {
            const int _idx = _ptr_array_neurongroup__spikespace[_index_spikes];
            double n_spikes = _ptr_array_neurongroup_n_spikes[0];
            double v;
            v = 0;
            n_spikes += 1;
            _ptr_array_neurongroup_1_n_spikes[0] = n_spikes;
            _ptr_array_neurongroup_1_v[_idx] = v;
        }

But in Python, this does not work -- we are always increasing n_spikes by 1 at every time step:

_idx = _array_neurongroup__spikespace[:_array_neurongroup__spikespace[-1]]
n_spikes = _array_neurongroup_n_spikes[0]
v = 0
n_spikes += 1
_array_neurongroup_n_spikes[0] = n_spikes
_array_neurongroup_v[_idx] = v

Any idea how we could deal with it without doing excessive special-casing of scalar variables? We could of course not allow to write to scalar variables in abstract code at all, but that would also mean no string expressions (which are of not much use since we can only refer to other scalar variables).

[thesamovar]

I don't think we can consistently allow writing to scalar variables in the same way as vectors. For example, what should scalar_variable = vector_variable do? It's undefined. What would be well defined would be scalar_variable += vector_variable or scalar_variable += constant (*= works too), but I'm not sure these are worth spending a lot of effort on. The latter is just counting neurons or spikes, and while the former is something we talked about, I think it would be better handled in the context of reductions (#180).

My feeling is that we just don't allow writing to scalar variables. You raise a very good point though that we do want to be able to use string expressions to write to them, probably. Not sure what the best design is here. It sort of depends on what we are going to use them for. The first thing will be for state updaters and loop invariant optimisations, and stuff like that, right? What else?

[mstimberg]

My feeling is that we just don't allow writing to scalar variables. You raise a very good point though that we do want to be able to use string expressions to write to them, probably. Not sure what the best design is here. It sort of depends on what we are going to use them for. The first thing will be for state updaters and loop invariant optimisations, and stuff like that, right? What else?

I guess you are right that for the moment we should simply not allow writing to them. And I didn't intend to talk about statements mixing scalar and vector variables for now, as you say this is something we should get into when dealing with #180.

To illustrate what I'd like them to use for, you can have a look at my (not yet working) Brian2 version of the barrelcortex example: https://github.com/brian-team/brian2/blob/barrelcortex_example/examples/synapses_barrelcortex.py

In the Brian1 code, there is a network operation run every time step and I want to avoid this in the Brian2 code by having only one bit of abstract code run in a CodeRunner that draws a new random stimulus. Everything else should be expressed related to a couple of scalar state variables that encode stimulus information. The runner code then looks like:

direction = rand()*2*pi
stim_start_x = 2.5 - cos(direction)*stimradius
stim_start_y = 2.5 - sin(direction)*stimradius
stim_start_time = t

Now, apart from wasting memory if direction, stim_start_x, etc. are vectors, the direction = rand()*2*pi will draw a different direction for each neuron, which is not what we want. But making this work would mean some fundamental change to codegen. I guess templates would need two code sections, one for scalar variables and one (the currently existing one) for vectors. Still, what that would mean in the context of conditional statements (resets, pre/post, ...) is a bit unclear...

Coming back to my use case from above: even without being able to write to scalar variables in code, I can of course implement this in a network operation, but then I will no longer be able to run the example in standalone. Actually, in the current situation there is a hackish workaround to get everything working, using linked variables (which one had to do manually until the implementation of #189): in a separate NeuronGroup of size 1, direction etc. can be declared as scalar parameters and updated using a runner as shown above -- since the group has size 1, this will do the right thing. Linking the variables into the original group would then solve the problem.

Back to practical matters. My proposal:

• Short-term solution: Do not allow to write to scalar variables in abstract code, not even with string initialization (which is of limited use in this case anyway)
• Long-term solution: Templates have an additional "scalar code" block that is outside of the loop in C++, sets the right vectorisation_idx, etc. If a template doesn't have such a block, this is meant to indicate that writing to scalar variables in this context is not well defined (e.g. in a reset). Non-scalar variables appearing on the right-hand-side of an assignment to a scalar variable raise an error.

[thesamovar]

From the programming side, I'm happy with your proposal. I'm wondering what the user syntax should be though, because it sounds like it could easily get ugly.

[mstimberg]

I'm wondering what the user syntax should be though, because it sounds like it could easily get ugly.

I don't think we necessarily need any new user syntax. It's true that it's not entirely obvious at the first glance whether a statement like a = b + c is a vectorized or scalar operation but at least it is not ambiguous: if a is a scalar, then b and c have to be scalar as well.

[thesamovar]

Might it not lead to potentially unexpected results? e.g. if you do count += 1 and expect it to tell you how many neurons/spikes there are.

[mstimberg]

Might it not lead to potentially unexpected results? e.g. if you do count += 1 and expect it to tell you how many neurons/spikes there are.

Well, would you expect it to do that :) ? I think if you write scalar += 1 in a CodeRunner that is executed every once in a while, you'd expect scalar to increase by 1 every time it is executed and nothing else. I agree that it would be unclear in a reset statement or a synaptic statement but according to my proposal, you are not allowed to assign to scalar variables there.

[thesamovar]

OK I see what you mean. Two possibilities to further reduce the possibility of unexpected behaviour:

• More restrictive: You can't mix scalar and vector statements in a single abstract code block.
• Less restrictive: Scalar and vector statements must be separate in a single code block, i.e. you can't do scalar_statement_1; vector_statement_1; scalar_statement_2; vector_statement_2 you'd have to do scalar_statement_1; scalar_statement_2; vector_statement_1; vector_statement_2. Maybe also make it so that all scalar statements have to come before the vector statements? The first time the user does it wrong, they get an error which tells them that they have to think about the difference between scalar and vector statements.

I think my slight preference is for the more restrictive option. What do you think?

[mstimberg]

OK I see what you mean. Two possibilities to further reduce the possibility of unexpected behaviour:

I have to say I still don't see the problem: what kind of unexpected behaviour are we talking about? If we do all the restrictions I mentioned (no conditional update semantics, no vector variables on the RHS of a scalar assignment) then I don't see any room for confusion. In fact, except for implicitly vectorized function calls like rand(), code should work exactly the same, regardless of whether the user actually thinks of declaring variables identical across neurons as scalar. In code where c is a vector variable:

a += 1
b += 1
c = a + b

it doesn't matter whether a and b are scalars or vectors, it only means they store the same value once or multiple times. I'd therefore not introduce any restrictions.

[thesamovar]

Good point! OK, let's go without restrictions. So what's the status of this then? Almost finished?

[mstimberg]

So what's the status of this then? Almost finished?

The "short-term solution" (no writing in code allowed) is implemented, I didn't start with the rest yet. Probably we should merge this already but keep this issue open until we have the final implementation. I don't think it's an incredible amount of work -- do you think we should try to get this in for the alpha or rather postpone it to beta?

[thesamovar]

Depends how much work it is. If it's not too much, let's include it because it's a really nice thing to have.