Accessing custom array during runtime

[michelem...@gmail.com]

Dear Developers and Community,

I searched a bit in the older threads but I couldn't find something that solves my problem. Unfortunately, I am not really practical with code so please forgive any mistakes you will see here.

I have an array of currents (analog values) registered in an event-driven fashion, therefore with nonuniform time distribution.

Something like: I_inj = [[1pA , 2pA ,  1pA],[10ms, 50ms, 100ms]]

I would like to input said currents into a neuron in an efficient way. My first guess was to use a custom event that would compare the current time t with the values in the array, at every event the index increase, moving the element to which I would compare.

Something like:

I_neuron = I_inj[0,i] (clock driven)

on event t = I_inj[1,i] => i += 1 (event-driven)

The problem I am facing is that I don't know how to handle the array.

The first approach I followed was this:

allParam = {
    '# lot of parameters} 

#defining the model

L1leaky_iaf = Equations('''
                         dv/dt = (v_rest-v)/(c_m*R_m) + (i_syn + i_offset + i_inj)/c_m : volt 
                         die/dt = -ie/L1tau_syn_e : amp
                         dii/dt = -ii/L1tau_syn_i : amp
                         i_inj : amp
                         i_syn = ie + ii : amp
                         tau_syn_e : second
                         tau_syn_i : second
                         ''')
#instancing the group with the custom event (now it's just a dummy condition)
layer1 = NeuronGroup(nNeurons_layer1, model=L1leaky_iaf,
                     threshold='v > v_thresh', refractory=cell_params['tau_refrac'],
                     reset='v = v_reset', name='Layer1', method='exact', events={'aaa': 'v>0*mV'})
#setting an initial current

layer1.i_inj = 1*nA
#adding an array
layer1.variables.add_array('time_array', dimensions=nA.dim,
                                 size=3, values = [10, 50,100])

#create the custom event
layer1.run_on_event('aaa', 'i_inj = time_array[0]') #note that I put i=0 just for debugging

run(1*second, namespace=allParam)

The runtime replies with an error: brian2.core.base.BrianObjectException: Unknown syntax: Subscript

So I thought that maybe the problem is related to time_array[0]. I, therefore, created a function:

@implementation('cython', '''
    cdef array_access(int i, float array):
        return array[i]
    ''')
@check_units(i=1, array = amp, result=amp)
def array_access(i,array):
    return array[i]

allParam = {
    '# lot of parameters, 'array_access' : array_access} # I added the function to the dict

#defining the model
L1leaky_iaf = Equations('''
                         dv/dt = (v_rest-v)/(c_m*R_m) + (i_syn + i_offset + i_inj)/c_m : volt 
                         die/dt = -ie/L1tau_syn_e : amp
                         dii/dt = -ii/L1tau_syn_i : amp
                         i_inj : amp
                         i_syn = ie + ii : amp
                         tau_syn_e : second
                         tau_syn_i : second
                         ''')
#instancing the group with the custom event (now it's just a dummy condition)
layer1 = NeuronGroup(nNeurons_layer1, model=L1leaky_iaf,
                     threshold='v > v_thresh', refractory=cell_params['tau_refrac'],
                     reset='v = v_reset', name='Layer1', method='exact', events={'dummy': 'v>0*mV'})
#setting an initial current
layer1.i_inj = 1*nA
#adding an array
layer1.variables.add_array('time_array', dimensions=nA.dim,
                                 size=3, values = [10, 50,100])
#create the custom event
layer1.run_on_event('dummy', 'i_inj = array_access(0,time_array)')

run(1*second, namespace=allParam)

But it replies with AttributeError: 'NoneType' object has no attribute 'main' whenever I try the operator something = array[I]   in the cython function.

If I now change the function to:

@implementation('cython', '''
    cdef array_access(int I,float array):
        print(array)
        return array
    ''')

It prints me the first value of the time_array whenever the event is triggered with no errors.

How can I access the values in the said array? 

Thank you very much, 

Regards,

Michele

[Marcel Stimberg]

Hi Michele,

Brian currently does not allow you to directly access arrays with index notation from Brian equations/statements. This is mostly for technical reasons, but partly also because it avoids users trying to think in programming language terms instead of describing things mathematically for a single neuron, etc.

That all being said, there are certainly ways to get it to work.

The easiest solution is to consider whether you cannot sacrifice some efficiency for simplicity. If your actual example is of comparable complexity as your example current:

I_inj = [[1pA , 2pA ,  1pA],[10ms, 50ms, 100ms]]

then you could also simply use a TimedArray like this:

I_inj = TimedArray([0*pA, 1*pA, 1*pA, 1*pA, 1*pA, 2*pA, 2*pA, 2*pA, 2*pA, 2*pA, 2*pA], dt=10*ms)

and then refer to I_inj(t) in your equations.

If that does not work for your use case, you can indeed use a user-defined function, but the way to hand over the array is a bit different. The "add_array" function is meant for internal use, in principle you could use it in a way that you need, but you'd have to do some complicated stuff with indexing variables and things would look quite odd in your equations. Instead, you want to hand over the array in the functions "namespace" (https://brian2.readthedocs.io/en/stable/advanced/functions.html#additional-namespace) – this is what TimedArray uses as well.

So it would look rather like this (I did not run the code, so there might be an error in there):

currents = [1, 2, 1]*pA
@implementation('cython', '''
    cdef get_current(int event_idx):
        return _namespace_currents[event_idx]
    ''',
    namespace={'_currents': currents})
@check_units(event_idx=1, result=amp)
@declare_types(event_idx='integer')
def get_current(event_idx):
    pass  # no Python implementation needed

You'd then also need to add some sort of counter variable to your equations, e.g.

event_idx : integer (shared)  # leave away (shared) if you want one counter per neuron

and then you can set the current and increase the counter on every event:

layer1.run_on_event('dummy', '''i_inj = get_current(event_idx)
                                event_idx += 1''')

Hope that works for you, best

  Marcel