Current sources

[andrewgait]

Resolves #1060

Allows the user to specify a CurrentSource and inject it into a Population or PopulationView (see e.g. http://neuralensemble.org/docs/PyNN/injecting_current.html).

Edit: This should now work for multiple CurrentSources into the same core and indeed into the same neuron.

Another unfinished piece of this is to deal with the case where the dt parameter in the NoisyCurrentSource is not the same as the (simulation) timestep value. For now the code throws an error if these two parameters are different.

An example is shown in SpiNNakerManchester/PyNNExamples#66; more detailed tests could possibly be written to check numerics etc. but not necessarily in this PR.

Merge with:
SpiNNakerManchester/sPyNNakerNewModelTemplate#76
SpiNNakerManchester/PyNNExamples#66 (which shows an example)

Tested by:
SpiNNakerManchester/IntegrationTests#51

[rowleya]

This looks good in general. The only issue I have is the repetition of code. A thought would be to have something like the full-neuron vs. components impl. This might have current_source_inst.h which defines a single current source instance (then with current_source_inst_noisy.h etc.) and then a current_source_impl.h to bring those things together. It may need a bit of thought though to get it right e.g. if there is a way to detect what sources are allowed etc. in some more automatic way.

[rowleya]

Thinking about this a bit more, I think it is likely that the main thing that brings the current sources together needs to know about all possible sources anyway, as you have done with the IDs. This could then maybe do something with #define and #ifdef in the right places, then do different things depending on which sources end up being included.

An example of this might be something like the following (note lots of copy and paste and totally untested!):

current_source_dc_impl.h

#ifndef CURRENT_SOURCE_DC_IMPL_H
#define CURRENT_SOURCE_DC_IMPL_H
typedef struct dc_source_t {
    REAL amplitude;
    uint32_t start;
    uint32_t stop;
} dc_source_t;

static bool current_source_noisy_init(void **data, void **dc_source) {
    dc_source_t *sdram_data = *data;
    *data = &sdram_data[1];
    dc_source_t *dtcm_data = spin1_malloc(sizeof(dc_source_t));
    *dtcm_data = *sdram_data;
    *dc_source = dtcm_data;
    return true;
}

static REAL current_source_get_dc_offset(void* current_source_data, uint32_t time, uint32_t neuron_index) {
    dc_source_t *dc_source = current_source_data;
    if ((time >= dc_source->start) && (time < dc_source->stop)) {
        return dc_source->amplitude;
    }
    return ZERO;
}

#endif

current_sources.h

typedef struct current_source_item {
    uint32_t source_id;
    uint32_t n_neuron_ids;
    uint32_t neuron_ids[];
    // Followed by data for the item
} current_source_item_t;

typedef struct current_source_data {
    uint32_t n_current_sources;
    current_source_item_t item[];
} current_source_data_t;

typedef struct current_source {
    uint32_t id;
    void *data;
    uint32_t n_neuron_ids;
    uint32_t *neuron_ids;
} current_source_t;

static uint32_t n_current_sources;

static current_source_t *current_sources;

static bool current_source_impl_initialise(current_source_data_t *cs_data) {

    // Read from cs_address; the first value is the number of current sources
    n_current_sources = cs_data->n_current_sources;
    current_sources = spin1_malloc(n_current_sources * sizeof(current_source_t));

    // Loop over number of current sources and get ID list for each
    current_source_item_t *item = &cs_data->item[0];
    for (uint32_t ncs=0; ncs < n_current_sources; ncs++)  {
        // Get neuron ids
        uint32_t n_neuron_ids =item->n_neuron_ids;
        current_sources[ncs]->n_neuron_ids = n_neuron_ids;
        current_sources[ncs]->neuron_ids = spin1_malloc(n_neuron_ids * sizeof(uint32_t));
        spin1_memcpy(current_sources[ncs]->neuron_ids, item->neuron_ids, n_neuron_ids * sizeof(uint32_t));

        // Data comes after neuron ids
        void *data = &item->neuron_ids[n_neuron_ids];

        // Setup source
        if (cs_data->item[ncs].source_id == 1) {
            #ifdef CURRENT_SOURCE_DC_IMPL_H
            if (!current_source_dc_impl_init(&data, &current_sources[ncs].data)) {
                return false;
            }
            #else
            log_error("DC current source is not supported at source %u", ncs);
            return false;
            #endif
        } else if (current_source[ncs]->current_source_id == 2) {
            // TODO
        } else if (current_source[ncs]->current_source_id == 3) {
            // TODO
        } else if (current_source[ncs]->current_source_id == 4) {
            // TODO
        }
    
        // Next item comes after the data, updated when the source was initialised
        item = (current_source_item_t *) data;
    }
}

static inline REAL current_source_get_offset(uint32_t time, uint32_t neuron_index) {
    // Avoid the loops if no current sources
    // TODO: Add the rest
    #if !defined(CURRENT_SOURCE_DC_IMPL_H) && !defined(CURRENT_SOURCE_AC_IMPL_H) && ...
    return ZERO;
    #else
    
    // Could simply just have the different cases in here
    // TODO: use an enum or something like that instead for the IDs
    REAL current_offset = ZERO;

    // Loop over current sources
    for (uint32_t n_cs=0; n_cs < n_current_sources; n_cs++) {
        current_source_t *source = &current_sources[n_cs];
        // Loop over neuron IDs associated with this current source
        uint32_t n_neuron_ids = source->n_neuron_ids;
        for (uint32_t n=0; n < n_neuron_ids; n++) {
            // If neuron ID matches the index value we are currently at
            if (neuron_index == source->neuron_ids[n]) {
                #ifdef CURRENT_SOURCE_DC_IMPL_H
                // Now do the appropriate calculation based on the ID value
                if (cs_id == 1) {  // DCSource
                    current_offset += current_source_get_dc_offset(source->data, time, n);
                }
                #endif
                // TODO: Fill in the rest
                #ifdef CURRENT_SOURCE_AC_IMPL_H
                #endif
            }
        }
    }

    return current_offset;
    #endif
}

[rowleya]

I was also just thinking how to avoid looping every neuron ID when trying the source. Another idea is below (again untested):

current_sources.h

typedef struct current_source_item {
    uint32_t source_id;
    // Followed by data for the item
} current_source_item_t;

typedef struct neuron_to_sources {
    uint32_t n_sources;
    uint32_t source_index[];
} neuron_to_sources_t;

typedef struct current_source_data {
    uint32_t n_current_sources;
    current_source_item_t item[];
    // Followed by
    // neuron_to_sources_t[n_neurons];
} current_source_data_t;

typedef struct current_source {
    uint32_t id;
    void *data;
} current_source_t;

typedef struct neuron_current_sources {
    uint32_t n_sources;
    current_source_t *source_indices;
} neuron_current_sources_t;

static uint32_t n_current_sources;

static current_source_t *current_sources;

static neuron_current_sources_t *neuron_current_sources;

static bool current_source_impl_initialise(current_source_data_t *cs_data, uint32_t n_neurons) {
    // Avoid the loops if no current sources
    // TODO: Add the rest
    #if !defined(CURRENT_SOURCE_DC_IMPL_H) && !defined(CURRENT_SOURCE_AC_IMPL_H) && ...
    return ZERO;
    #else

    // Read from cs_address; the first value is the number of current sources
    n_current_sources = cs_data->n_current_sources;
    current_sources = spin1_malloc(n_current_sources * sizeof(current_source_t));
    neuron_current_sources = spin1_malloc(n_neurons * sizeof(neuron_current_sources_t));

    // Loop over number of current sources and get ID list for each
    current_source_item_t *item = &cs_data->item[0];
    for (uint32_t ncs=0; ncs < n_current_sources; ncs++)  {

        // Data comes after neuron ids
        void *data = &item[1];

        // Setup source
        if (cs_data->item[ncs].source_id == 1) {
            #ifdef CURRENT_SOURCE_DC_IMPL_H
            if (!current_source_dc_impl_init(&data, &current_sources[ncs].data)) {
                return false;
            }
            #else
            log_error("DC current source is not supported at source %u", ncs);
            return false;
            #endif
        } else if (current_source[ncs]->current_source_id == 2) {
            // TODO
        } else if (current_source[ncs]->current_source_id == 3) {
            // TODO
        } else if (current_source[ncs]->current_source_id == 4) {
            // TODO
        }
    
        // Next item comes after the data, updated when the source was initialised
        item = (current_source_item_t *) data;
    }

    // Neuron mapping data follows
    neuron_to_sources_t *neuron_mapping = (neuron_to_sources_t *) item;
    for (uint32_t n = 0; n < n_neurons; n++) {
        uint32_t n_sources = neuron_mapping->n_sources;
        neuron_current_sources[n].n_sources = n_sources;
        neuron_current_sources[n].source_indices = spin1_malloc(n_sources * sizeof(uint32_t));
        spin1_memcpy(neuron_current_sources[n].source_indices, neuron_mapping->source_index, n_sources * sizeof(uint32_t));
        neuron_mapping = (neuron_to_sources_t *) &neuron_mapping->source_index[n_sources];
    }
    
    #endif
}

static inline REAL current_source_get_offset(uint32_t time, uint32_t neuron_index) {
    // Avoid the loops if no current sources
    // TODO: Add the rest
    #if !defined(CURRENT_SOURCE_DC_IMPL_H) && !defined(CURRENT_SOURCE_AC_IMPL_H) && ...
    return ZERO;
    #else
    
    // Could simply just have the different cases in here
    // TODO: use an enum or something like that instead for the IDs
    REAL current_offset = ZERO;

    // Loop over current sources
    for (uint32_t n_cs=0; n_cs < neuron_current_sources[neuron_id].n_sources; n_cs++) {
        current_source_t *source = &current_sources[neuron_current_sources[neuron_id].source_indices[n_cs]];
       
        #ifdef CURRENT_SOURCE_DC_IMPL_H
        // Now do the appropriate calculation based on the ID value
        if (cs_id == 1) {  // DCSource
            current_offset += current_source_get_dc_offset(source->data, time, n);
        }
        #endif
        // TODO: Fill in the rest
        #ifdef CURRENT_SOURCE_AC_IMPL_H
        #endif
    }

    return current_offset;
    #endif
}

[andrewgait]

Both suggestions seem plausible - will investigate and return with either updated code or comments on the suggestions.

[andrewgait]

This now has much less repetition in than previously - it's now an amalgamation between your first suggestion and what was originally on the branch. I couldn't get the various structs you'd suggested to work properly but I think it must be possible, and the second idea to avoid looping every neuron ID must also be possible but again I haven't yet got the memory structure working correctly in that instance.

[andrewgait]

I think I now have this running correctly to avoid the loop to check the neuron ID each time. I did this slightly differently to the suggestion in that I changed the array being written in from the Python side to have the current source IDs per neuron rather than the neuron IDs per current source as it previously was.

[rowleya]

I think this is fine now. I made one comment about a bit of code that looks a bit odd, but it isn't critical.