Issue/241 regular spiking source cell (#287)

Replace the 20 artificial spikes we inserted until now with spikes generated with a frequency spiking cell.
This is a cell that spikes regularly for a set period.
This connections outgoing from this cell mirror all the outgoing connections of cells with gid % 20 == 0.
This means that 5% extra connections have been added to the model to implement the equivalent behaviour with an 'implemented' neuron.

To implement this behavoir a new cell_group type has been introduced.
The fs_cell is extremely small and is not implemented on the backends.

Includes:
Unit test for the rss_cell functionality

miniapp/miniapp.cpp

@@ -105,11 +105,6 @@ int main(int argc, char** argv) {
 
         m.set_binning_policy(binning_policy, options.bin_dt);
 
-        // Inject some artificial spikes, 1 per 20 neurons.
-        for (cell_gid_type c=0; c<recipe->num_cells(); c+=20) {
-            m.add_artificial_spike({c, 0});
-        }
-
         // Attach samplers to all probes
         std::vector<std::unique_ptr<sample_trace_type>> traces;
         const time_type sample_dt = options.sample_dt;
@@ -195,6 +190,8 @@ std::unique_ptr<recipe> make_recipe(const io::cl_options& options, const probe_d
     p.morphology_round_robin = options.morph_rr;
 
     p.num_compartments = options.compartments_per_segment;
+
+    // TODO: Put all recipe parameters in the recipes file
     p.num_synapses = options.all_to_all? options.cells-1: options.synapses_per_cell;
     p.synapse_type = options.syn_type;
 
@@ -272,7 +269,7 @@ void report_compartment_stats(const recipe& rec) {
 
     for (std::size_t i = 0; i<ncell; ++i) {
         std::size_t ncomp = 0;
-        auto c = rec.get_cell(i);
+        auto c = rec.get_cell_description(i);
         if (auto ptr = util::any_cast<cell>(&c)) {
             ncomp = ptr->num_compartments();
         }

miniapp/miniapp_recipes.cpp

@@ -4,9 +4,9 @@
 #include <utility>
 
 #include <cell.hpp>
+#include <rss_cell.hpp>
 #include <morphology.hpp>
 #include <util/debug.hpp>
-#include <util/unique_any.hpp>
 
 #include "miniapp_recipes.hpp"
 #include "morphology_pool.hpp"
@@ -79,9 +79,17 @@ public:
                             - param_.min_connection_delay_ms};
     }
 
-    cell_size_type num_cells() const override { return ncell_; }
+    cell_size_type num_cells() const override {
+        return ncell_ + 1;  // We automatically add a fake cell to each recipe!
+    }
+
+    util::unique_any get_cell_description(cell_gid_type i) const override {
+        // The last 'cell' is a rss_cell with one spike at t=0
+        if (i == ncell_) {
+            return util::unique_any(std::move(
+                rss_cell::rss_cell_description(0.0, 0.1, 0.1) ));
+        }
 
-    util::unique_any get_cell(cell_gid_type i) const override {
         auto gen = std::mt19937(i); // TODO: replace this with hashing generator...
 
         auto cc = get_cell_count_info(i);
@@ -113,8 +121,11 @@ public:
         return util::unique_any(std::move(cell));
     }
 
-    cell_kind get_cell_kind(cell_gid_type) const override {
-        // The basic_cell_recipe only produces mc cells, so return cable1d_neuron for now
+    cell_kind get_cell_kind(cell_gid_type i ) const override {
+        // The last 'cell' is a rss_cell with one spike at t=0
+        if (i == ncell_) {
+            return cell_kind::regular_spike_source;
+        }
         return cell_kind::cable1d_neuron;
     }
 
@@ -175,6 +186,12 @@ public:
 
     std::vector<cell_connection> connections_on(cell_gid_type i) const override {
         std::vector<cell_connection> conns;
+
+        // The rss_cell does not have inputs
+        if (i == ncell_) {
+            return conns;
+        }
+
         auto gen = std::mt19937(i); // TODO: replace this with hashing generator...
 
         cell_gid_type prev = i==0? ncell_-1: i-1;
@@ -183,6 +200,13 @@ public:
             cc.source = {prev, 0};
             cc.dest = {i, t};
             conns.push_back(cc);
+
+            // The rss_cell spikes at t=0, with these connections it looks like
+            // (source % 20) == 0 spikes at that moment.
+            if (prev % 20 == 0) {
+                cc.source = {ncell_, 0}; // also add connection from reg spiker!
+                conns.push_back(cc);
+            }
         }
 
         return conns;
@@ -207,6 +231,11 @@ public:
 
     std::vector<cell_connection> connections_on(cell_gid_type i) const override {
         std::vector<cell_connection> conns;
+
+        // The rss_cell does not have inputs
+        if (i == ncell_) {
+            return conns;
+        }
         auto conn_param_gen = std::mt19937(i); // TODO: replace this with hashing generator...
         auto source_gen = std::mt19937(i*123+457); // ditto
 
@@ -220,6 +249,13 @@ public:
             cc.source = {source, 0};
             cc.dest = {i, t};
             conns.push_back(cc);
+
+            // The rss_cell spikes at t=0, with these connections it looks like
+            // (source % 20) == 0 spikes at that moment.
+            if (source % 20 == 0) {
+                cc.source = {ncell_, 0};
+                conns.push_back(cc);
+            }
         }
 
         return conns;
@@ -249,6 +285,10 @@ public:
 
     std::vector<cell_connection> connections_on(cell_gid_type i) const override {
         std::vector<cell_connection> conns;
+        // The rss_cell does not have inputs
+        if (i == ncell_) {
+            return conns;
+        }
         auto conn_param_gen = std::mt19937(i); // TODO: replace this with hashing generator...
 
         for (unsigned t=0; t<param_.num_synapses; ++t) {
@@ -259,6 +299,13 @@ public:
             cc.source = {source, 0};
             cc.dest = {i, t};
             conns.push_back(cc);
+
+            // The rss_cell spikes at t=0, with these connections it looks like
+            // (source % 20) == 0 spikes at that moment.
+            if (source % 20 == 0) {
+                cc.source = {ncell_, 0};
+                conns.push_back(cc);
+            }
         }
 
         return conns;

src/cell_group_factory.cpp

@@ -2,6 +2,7 @@
 
 #include <backends.hpp>
 #include <cell_group.hpp>
+#include <rss_cell_group.hpp>
 #include <fvm_multicell.hpp>
 #include <mc_cell_group.hpp>
 #include <util/unique_any.hpp>
@@ -15,21 +16,21 @@ using mc_fvm_cell = mc_cell_group<fvm::fvm_multicell<multicore::backend>>;
 cell_group_ptr cell_group_factory(
         cell_kind kind,
         cell_gid_type first_gid,
-        const std::vector<util::unique_any>& cells,
+        const std::vector<util::unique_any>& cell_descriptions,
         backend_policy backend)
 {
-    if (backend==backend_policy::prefer_gpu) {
-        switch (kind) {
-        case cell_kind::cable1d_neuron:
-            return make_cell_group<gpu_fvm_cell>(first_gid, cells);
-        default:
-            throw std::runtime_error("unknown cell kind");
-        }
-    }
-
     switch (kind) {
     case cell_kind::cable1d_neuron:
-        return make_cell_group<mc_fvm_cell>(first_gid, cells);
+        if (backend == backend_policy::prefer_gpu) {
+            return make_cell_group<gpu_fvm_cell>(first_gid, cell_descriptions);
+        }
+        else {
+            return make_cell_group<mc_fvm_cell>(first_gid, cell_descriptions);
+        }
+
+    case cell_kind::regular_spike_source:
+        return make_cell_group<rss_cell_group>(first_gid, cell_descriptions);
+
     default:
         throw std::runtime_error("unknown cell kind");
     }

src/common_types.hpp

@@ -57,7 +57,8 @@ using time_type = float;
 // group equal kinds in the same cell group.
 
 enum cell_kind {
-    cable1d_neuron           // Our own special mc neuron
+    cable1d_neuron,           // Our own special mc neuron
+    regular_spike_source,     // Regular spiking source
 };
 
 } // namespace mc

src/mc_cell_group.hpp

@@ -35,25 +35,25 @@ public:
     mc_cell_group() = default;
 
     template <typename Cells>
-    mc_cell_group(cell_gid_type first_gid, const Cells& cells):
+    mc_cell_group(cell_gid_type first_gid, const Cells& cell_descriptions):
         gid_base_{first_gid}
     {
         // Create lookup structure for probe and target ids.
-        build_handle_partitions(cells);
+        build_handle_partitions(cell_descriptions);
         std::size_t n_probes = probe_handle_divisions_.back();
         std::size_t n_targets = target_handle_divisions_.back();
         std::size_t n_detectors = algorithms::sum(util::transform_view(
-            cells, [](const cell& c) { return c.detectors().size(); }));
+            cell_descriptions, [](const cell& c) { return c.detectors().size(); }));
 
         // Allocate space to store handles.
         target_handles_.resize(n_targets);
         probe_handles_.resize(n_probes);
 
-        lowered_.initialize(cells, target_handles_, probe_handles_);
+        lowered_.initialize(cell_descriptions, target_handles_, probe_handles_);
 
         // Create a list of the global identifiers for the spike sources
         auto source_gid = cell_gid_type{gid_base_};
-        for (const auto& cell: cells) {
+        for (const auto& cell: cell_descriptions) {
             for (cell_lid_type lid=0u; lid<cell.detectors().size(); ++lid) {
                 spike_sources_.push_back(source_id_type{source_gid, lid});
             }
@@ -63,9 +63,9 @@ public:
 
         // Create the enumeration of probes attached to cells in this cell group
         probes_.reserve(n_probes);
-        for (auto i: util::make_span(0, cells.size())){
+        for (auto i: util::make_span(0, cell_descriptions.size())){
             const cell_gid_type probe_gid = gid_base_ + i;
-            const auto probes_on_cell = cells[i].probes();
+            const auto probes_on_cell = cell_descriptions[i].probes();
             for (cell_lid_type lid: util::make_span(0, probes_on_cell.size())) {
                 // get the unique global identifier of this probe
                 cell_member_type id{probe_gid, lid};
@@ -79,11 +79,11 @@ public:
         }
     }
 
-    mc_cell_group(cell_gid_type first_gid, const std::vector<util::unique_any>& cells):
+    mc_cell_group(cell_gid_type first_gid, const std::vector<util::unique_any>& cell_descriptions):
         mc_cell_group(
             first_gid,
             util::transform_view(
-                cells,
+                cell_descriptions,
                 [](const util::unique_any& c) -> const cell& {return util::any_cast<const cell&>(c);})
         )
     {}

src/model.cpp

@@ -31,15 +31,15 @@ model::model(const recipe& rec, const domain_decomposition& decomp):
             PE("setup", "cells");
 
             auto group = domain_.get_group(i);
-            std::vector<util::unique_any> cells(group.end-group.begin);
+            std::vector<util::unique_any> cell_descriptions(group.end-group.begin);
 
             for (auto gid: util::make_span(group.begin, group.end)) {
                 auto i = gid-group.begin;
-                cells[i] = rec.get_cell(gid);
+                cell_descriptions[i] = rec.get_cell_description(gid);
             }
 
             cell_groups_[i] = cell_group_factory(
-                    group.kind, group.begin, cells, domain_.backend());
+                    group.kind, group.begin, cell_descriptions, domain_.backend());
             PL(2);
         });
 

src/recipe.hpp

@@ -48,7 +48,7 @@ class recipe {
 public:
     virtual cell_size_type num_cells() const =0;
 
-    virtual util::unique_any get_cell(cell_gid_type) const =0;
+    virtual util::unique_any get_cell_description(cell_gid_type) const =0;
     virtual cell_kind get_cell_kind(cell_gid_type) const = 0;
 
     virtual cell_count_info get_cell_count_info(cell_gid_type) const =0;
@@ -70,7 +70,7 @@ public:
         return 1;
     }
 
-    util::unique_any get_cell(cell_gid_type) const override {
+    util::unique_any get_cell_description(cell_gid_type) const override {
         return util::unique_any(cell(clone_cell, cell_));
     }
 

src/rss_cell.hpp

+#pragma once
+#pragma once
+
+#include <vector>
+
+#include <common_types.hpp>
+#include <util/debug.hpp>
+
+namespace nest {
+namespace mc {
+
+/// Regular spike source: A cell that generated spikes with a certain
+/// period for a set time range.
+class rss_cell {
+public:
+    using index_type = cell_lid_type;
+    using size_type = cell_local_size_type;
+    using value_type = double;
+
+    struct rss_cell_description {
+        time_type start_time;
+        time_type period;
+        time_type stop_time;
+
+        rss_cell_description(time_type st, time_type per, time_type stop):
+            start_time(st),
+            period(per),
+            stop_time(stop)
+        {}
+    };
+
+    /// Construct a rss cell from its description
+    rss_cell(rss_cell_description descr) :
+        start_time_(descr.start_time),
+        period_(descr.period),
+        stop_time_(descr.stop_time),
+        time_(0.0)
+    {}
+
+    /// Construct a rss cell from individual arguments
+    rss_cell(time_type start_time=0.0, time_type period=1.0, time_type stop_time=0.0):
+        start_time_(start_time),
+        period_(period),
+        stop_time_(stop_time),
+        time_(0.0)
+    {}
+
+    /// Return the kind of cell, used for grouping into cell_groups
+    cell_kind  get_cell_kind() const  {
+        return cell_kind::regular_spike_source;
+    }
+
+    /// Collect all spikes until tfinal.
+    // updates the internal time state to tfinal as a side effect
+    std::vector<time_type> spikes_until(time_type tfinal) {
+        std::vector<time_type> spike_times;
+
+        // If we should be spiking in this 'period'
+        if (tfinal > start_time_) {
+            // We have to spike till tfinal or the stop_time_of the neuron
+            auto end_time = stop_time_ < tfinal ? stop_time_ : tfinal;
+            // Generate all possible spikes in this time frame typically only
+            for (time_type time = start_time_ > time_ ? start_time_ : time_;
+                time < end_time;
+                time += period_) {
+                spike_times.push_back(time);
+            }
+        }
+        // Save our current time we generate exclusive a possible spike at tfinal
+        time_ = tfinal;
+
+        return spike_times;
+    }
+
+    /// reset internal time to 0.0
+    void reset() {
+        time_ = 0.0;
+    }
+
+private:
+    // When to start spiking
+    time_type start_time_;
+
+    // with what period
+    time_type period_;
+
+    // untill when
+    time_type stop_time_;
+
+    // internal time, storage
+    time_type time_;
+};
+
+} // namespace mc
+} // namespace nest

src/rss_cell_group.hpp

+#pragma once
+
+#include <cell_group.hpp>
+#include <rss_cell.hpp>
+#include <util/span.hpp>
+#include <util/unique_any.hpp>
+
+
+#include <iostream>
+namespace nest {
+namespace mc {
+
+/// Cell_group to collect cells that spike at a set frequency
+/// Cell are lightweight and are not executed in anybackend implementation
+class rss_cell_group : public cell_group {
+public:
+    using source_id_type = cell_member_type;
+
+    rss_cell_group(cell_gid_type first_gid, const std::vector<util::unique_any>& cell_descriptions):
+        gid_base_(first_gid)
+    {
+        using util::make_span;
+
+        for (cell_gid_type i: make_span(0, cell_descriptions.size())) {
+            // Copy all the rss_cells
+            cells_.push_back(rss_cell(
+                util::any_cast<rss_cell::rss_cell_description>(cell_descriptions[i])
+            ));
+
+            // create a lid to gid map
+            spike_sources_.push_back({gid_base_+i, 0});
+        }
+    }
+
+    virtual ~rss_cell_group() = default;
+
+    cell_kind get_cell_kind() const override {
+        return cell_kind::regular_spike_source;
+    }
+
+    void reset() override {
+        for (auto cell: cells_) {
+            cell.reset();
+        }
+    }
+
+    void set_binning_policy(binning_kind policy, time_type bin_interval) override
+    {} // Nothing to do?
+
+    void advance(time_type tfinal, time_type dt) override {
+        // TODO: Move source information to rss_cell implementation
+        for (auto i: util::make_span(0, cells_.size())) {
+            for (auto spike_time: cells_[i].spikes_until(tfinal)) {
+                spikes_.push_back({spike_sources_[i], spike_time});
+            }
+        }
+    };
+
+    void enqueue_events(const std::vector<postsynaptic_spike_event>& events) override {
+        std::logic_error("The rss_cells do not support incoming events!");
+    }
+
+    const std::vector<spike>& spikes() const override {
+        return spikes_;
+    }
+
+    void clear_spikes() override {
+        spikes_.clear();
+    }
+
+    std::vector<probe_record> probes() const override {
+        return probes_;
+    }
+
+    void add_sampler(cell_member_type probe_id, sampler_function s, time_type start_time = 0) override {
+        std::logic_error("The rss_cells do not support sampling of internal state!");
+    }
+
+private:
+    // gid of first cell in group.
+    cell_gid_type gid_base_;
+
+    // Spikes that are generated.
+    std::vector<spike> spikes_;
+
+    // Spike generators attached to the cell
+    std::vector<source_id_type> spike_sources_;
+
+    // Store a reference to the cell actually implementing the spiking
+    std::vector<rss_cell> cells_;
+
+    std::vector<probe_record> probes_;
+};
+
+} // namespace mc
+} // namespace nest
+

tests/CMakeLists.txt

@@ -16,6 +16,8 @@ add_subdirectory(performance)
 # Microbenchmarks.
 add_subdirectory(ubench)
 
+# regression / delta tests
+# Employing the full simulator. validated using deltas on output data
 
 # modcc tests
 if(NOT use_external_modcc)
@@ -25,9 +27,6 @@ endif()
 
 # Proposed additional test types:
 
-# Large test, employing the full simulator. validated using deltas on output data
-
 # Test to check integration between components
 
-
 # Numbered tests based on bugs in the tracker

tests/unit/test_domain_decomposition.cpp

@@ -17,7 +17,7 @@ public:
         return size_;
     }
 
-    util::unique_any get_cell(cell_gid_type) const override {
+    util::unique_any get_cell_description(cell_gid_type) const override {
         return {};
     }
     cell_kind get_cell_kind(cell_gid_type) const override {

tests/unit/test_rss_cell.cpp

+#include "../gtest.h"
+
+#include "rss_cell.hpp"
+
+using namespace  nest::mc;
+
+TEST(rss_cell, constructor)
+{
+    rss_cell test(0.0, 0.01, 1.0);
+}
+
+TEST(rss_cell, basic_usage)
+{
+    rss_cell sut(0.1, 0.01, 0.2);
+
+    // no spikes in this time frame
+    auto spikes = sut.spikes_until(0.09);
+    EXPECT_EQ(size_t(0), spikes.size());
+
+    //only on in this time frame
+    spikes = sut.spikes_until(0.11);
+    EXPECT_EQ(size_t(1), spikes.size());
+
+    // Reset the internal state to null
+    sut.reset();
+
+    // Expect 10 excluding the 0.2
+    spikes = sut.spikes_until(0.2);
+    EXPECT_EQ(size_t(10), spikes.size());
+}
+
+TEST(rss_cell, poll_time_after_end_time)
+{
+    rss_cell sut(0.1, 0.01, 0.2);
+
+    // no spikes in this time frame
+    auto spikes = sut.spikes_until(0.3);
+    EXPECT_EQ(size_t(10), spikes.size());
+
+    // now ask for spikes for a time slot already passed.
+    spikes = sut.spikes_until(0.2);
+    // It should result in zero spikes because of the internal state!
+    EXPECT_EQ(size_t(0), spikes.size());
+
+    sut.reset();
+
+    // Expect 10 excluding the 0.2
+    spikes = sut.spikes_until(0.2);
+    EXPECT_EQ(size_t(10), spikes.size());
+}
+
+TEST(rss_cell, cell_kind_correct)
+{
+    rss_cell sut(0.1, 0.01, 0.2);
+
+    EXPECT_EQ(cell_kind::regular_spike_source, sut.get_cell_kind());
+}