diff --git a/miniapps/miniapp/miniapp_recipes.cpp b/miniapps/miniapp/miniapp_recipes.cpp
index cae791fcb4761485409015399d6e417b288c3100..9ab4cf3ae0a46f36fb769a34f3de7bdf0ae7fec1 100644
--- a/miniapps/miniapp/miniapp_recipes.cpp
+++ b/miniapps/miniapp/miniapp_recipes.cpp
@@ -5,6 +5,7 @@
#include <cell.hpp>
#include <dss_cell_description.hpp>
+#include <event_generator.hpp>
#include <rss_cell.hpp>
#include <morphology.hpp>
#include <util/debug.hpp>
@@ -175,6 +176,10 @@ public:
}
}
+ std::vector<event_generator_ptr> event_generators(cell_gid_type) const override {
+ return {};
+ }
+
protected:
template <typename RNG>
cell_connection draw_connection_params(RNG& rng) const {
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 4e1b4db31a6c17b18f59f812d625f0c4d3ead903..fd29c43e070fbac8d35fe9e1bbc45934a1e72838 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -9,6 +9,7 @@ set(BASE_SOURCES
hardware/memory.cpp
hardware/node_info.cpp
hardware/power.cpp
+ merge_events.cpp
model.cpp
morphology.cpp
partition_load_balance.cpp
diff --git a/src/common_types.hpp b/src/common_types.hpp
index f817ea6e9cd92ead859f38d5a4ddfb24100375ad..2cc3203454943f20691caef3ffa109eaa4948c3d 100644
--- a/src/common_types.hpp
+++ b/src/common_types.hpp
@@ -7,6 +7,7 @@
#include <cstddef>
#include <functional>
+#include <limits>
#include <iosfwd>
#include <type_traits>
@@ -53,6 +54,7 @@ DEFINE_LEXICOGRAPHIC_ORDERING(cell_member_type,(a.gid,a.index),(b.gid,b.index))
// For storing time values [ms]
using time_type = float;
+constexpr time_type max_time = std::numeric_limits<time_type>::max();
// Extra contextual information associated with a probe.
diff --git a/src/event_generator.hpp b/src/event_generator.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2a461198c9f599e297ac4d5449659250de62fbbb
--- /dev/null
+++ b/src/event_generator.hpp
@@ -0,0 +1,243 @@
+#pragma once
+
+#include <cstdint>
+#include <memory>
+#include <random>
+
+#include <common_types.hpp>
+#include <event_queue.hpp>
+#include <util/range.hpp>
+#include <util/rangeutil.hpp>
+
+namespace arb {
+
+// An event_generator generates a sequence of events to be delivered to a cell.
+// The sequence of events is always in ascending order, i.e. each event will be
+// greater than the event that proceded it, where events are ordered by:
+// - delivery time;
+// - then target id for events with the same delivery time;
+// - then weight for events with the same delivery time and target.
+struct event_generator {
+ // Return the next event in the stream.
+ // Returns the same event if called multiple times without calling pop().
+ virtual postsynaptic_spike_event next() = 0;
+
+ // Move the generator to the next event in the stream.
+ virtual void pop() = 0;
+
+ // Reset the generator to the same state that it had on construction.
+ virtual void reset() = 0;
+
+ // Update state of the generator such that the event returned by next() is
+ // the first event with delivery time >= t.
+ virtual void advance(time_type t) = 0;
+
+ virtual ~event_generator() {};
+};
+
+inline
+postsynaptic_spike_event terminal_pse() {
+ return postsynaptic_spike_event{cell_member_type{0,0}, max_time, 0};
+}
+
+// Generator that feeds events that are specified with a vector.
+// Makes a copy of the input sequence of events.
+struct vector_backed_generator: public event_generator {
+ using pse = postsynaptic_spike_event;
+ vector_backed_generator(pse_vector events):
+ events_(std::move(events)),
+ it_(events_.begin())
+ {
+ if (!std::is_sorted(events_.begin(), events_.end())) {
+ util::sort(events_);
+ }
+ }
+
+ postsynaptic_spike_event next() override {
+ return it_==events_.end()? terminal_pse(): *it_;
+ }
+
+ void pop() override {
+ if (it_!=events_.end()) {
+ ++it_;
+ }
+ }
+
+ void reset() override {
+ it_ = events_.begin();
+ }
+
+ void advance(time_type t) override {
+ it_ = std::lower_bound(events_.begin(), events_.end(), t, event_time_less());
+ }
+
+private:
+ std::vector<postsynaptic_spike_event> events_;
+ std::vector<postsynaptic_spike_event>::const_iterator it_;
+};
+
+// Generator for events in a generic sequence.
+// The generator keeps a reference to a Seq, i.e. it does not own the sequence.
+// Care must be taken to avoid lifetime issues, to ensure that the generator
+// does not outlive the sequence.
+template <typename Seq>
+struct seq_generator: public event_generator {
+ using pse = postsynaptic_spike_event;
+ seq_generator(Seq& events):
+ events_(events),
+ it_(std::begin(events_))
+ {
+ EXPECTS(std::is_sorted(events_.begin(), events_.end()));
+ }
+
+ postsynaptic_spike_event next() override {
+ return it_==events_.end()? terminal_pse(): *it_;
+ }
+
+ void pop() override {
+ if (it_!=events_.end()) {
+ ++it_;
+ }
+ }
+
+ void reset() override {
+ it_ = events_.begin();
+ }
+
+ void advance(time_type t) override {
+ it_ = std::lower_bound(events_.begin(), events_.end(), t, event_time_less());
+ }
+
+private:
+
+ const Seq& events_;
+ typename Seq::const_iterator it_;
+};
+
+// Generates a set of regularly spaced events:
+// * with delivery times t=t_start+n*dt, ∀ t ∈ [t_start, t_stop)
+// * with a set target and weight
+struct regular_generator: public event_generator {
+ using pse = postsynaptic_spike_event;
+
+ regular_generator(cell_member_type target,
+ float weight,
+ time_type tstart,
+ time_type dt,
+ time_type tstop=max_time):
+ target_(target),
+ weight_(weight),
+ step_(0),
+ t_start_(tstart),
+ dt_(dt),
+ t_stop_(tstop)
+ {}
+
+ postsynaptic_spike_event next() override {
+ const auto t = time();
+ return t<t_stop_?
+ postsynaptic_spike_event{target_, t, weight_}:
+ terminal_pse();
+ }
+
+ void pop() override {
+ ++step_;
+ }
+
+ void advance(time_type t0) override {
+ t0 = std::max(t0, t_start_);
+ step_ = (t0-t_start_)/dt_;
+
+ // Finding the smallest value for step_ that satisfies the condition
+ // that time() >= t0 is unfortunately a horror show because floating
+ // point precission.
+ while (step_ && time()>=t0) {
+ --step_;
+ }
+ while (time()<t0) {
+ ++step_;
+ }
+ }
+
+ void reset() override {
+ step_ = 0;
+ }
+
+private:
+ time_type time() const {
+ return t_start_ + step_*dt_;
+ }
+
+ cell_member_type target_;
+ float weight_;
+ std::size_t step_;
+ time_type t_start_;
+ time_type dt_;
+ time_type t_stop_;
+};
+
+// Generates a stream of events at times described by a Poisson point process
+// with rate_per_ms spikes per ms.
+template <typename RandomNumberEngine>
+struct poisson_generator: public event_generator {
+ using pse = postsynaptic_spike_event;
+
+ poisson_generator(cell_member_type target,
+ float weight,
+ RandomNumberEngine rng,
+ time_type tstart,
+ time_type rate_per_ms,
+ time_type tstop=max_time):
+ exp_(rate_per_ms),
+ reset_state_(std::move(rng)),
+ target_(target),
+ weight_(weight),
+ t_start_(tstart),
+ t_stop_(tstop)
+ {
+ reset();
+ }
+
+ postsynaptic_spike_event next() override {
+ return next_<t_stop_?
+ postsynaptic_spike_event{target_, next_, weight_}:
+ terminal_pse();
+ }
+
+ void pop() override {
+ next_ += exp_(rng_);
+ }
+
+ void advance(time_type t0) override {
+ while (next_<t0) {
+ pop();
+ }
+ }
+
+ void reset() override {
+ rng_ = reset_state_;
+ next_ = t_start_;
+ pop();
+ }
+
+private:
+ std::exponential_distribution<time_type> exp_;
+ RandomNumberEngine rng_;
+ const RandomNumberEngine reset_state_;
+ const cell_member_type target_;
+ const float weight_;
+ const time_type t_start_;
+ const time_type t_stop_;
+ time_type next_;
+
+};
+
+using event_generator_ptr = std::unique_ptr<event_generator>;
+
+template <typename T, typename... Args>
+event_generator_ptr make_event_generator(Args&&... args) {
+ return event_generator_ptr(new T(std::forward<Args>(args)...));
+}
+
+} // namespace arb
+
diff --git a/src/merge_events.cpp b/src/merge_events.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e5cda151b3e4be3bbfc64e1cc7c1891a626fdfb7
--- /dev/null
+++ b/src/merge_events.cpp
@@ -0,0 +1,205 @@
+#include <set>
+#include <vector>
+
+#include <backends.hpp>
+#include <cell_group.hpp>
+#include <cell_group_factory.hpp>
+#include <domain_decomposition.hpp>
+#include <merge_events.hpp>
+#include <model.hpp>
+#include <recipe.hpp>
+#include <util/filter.hpp>
+#include <util/span.hpp>
+#include <util/unique_any.hpp>
+#include <profiling/profiler.hpp>
+
+namespace arb {
+
+namespace impl {
+
+// The tournament tree data structure is used to merge k sorted lists of events.
+// See online for high-level information about tournament trees.
+//
+// This implementation maintains a heap-like data structure, with entries of type:
+// std::pair<unsigned, post_synaptic_event>
+// where the unsigned ∈ [0, k-1] is the id of the list from which the event was
+// drawn. The id is stored so that the operation of removing the most recent event
+// knows which leaf node needs to be updated (i.e. the leaf node of the list from
+// which the most recent event was drawn).
+//
+// unsigned is used for storing the index, because if drawing events from more
+// event generators than can be counted using an unsigned a complete redesign
+// will be needed.
+tourney_tree::tourney_tree(std::vector<event_generator_ptr>& input):
+ input_(input),
+ n_lanes_(input_.size())
+{
+ // Must have at least 1 queue
+ EXPECTS(n_lanes_);
+ // Maximum value in unsigned limits how many queues we can have
+ EXPECTS(n_lanes_<(1u<<(sizeof(unsigned)*8u-1u)));
+
+ leaves_ = next_power_2(n_lanes_);
+ nodes_ = 2u*(leaves_-1u)+1u; // 2*l-1 with overflow protection
+
+ // Allocate space for the tree nodes
+ heap_.resize(nodes_);
+ // Set the leaf nodes
+ for (auto i=0u; i<leaves_; ++i) {
+ heap_[leaf(i)] = i<n_lanes_?
+ key_val(i, input[i]->next()):
+ key_val(i, terminal_pse()); // null leaf node
+ }
+ // Walk the tree to initialize the non-leaf nodes
+ setup(0);
+}
+
+void tourney_tree::print() const {
+ auto nxt=1u;
+ for (auto i=0u; i<nodes_; ++i) {
+ if (i==nxt-1) { nxt*=2; std::cout << "\n";}
+ std::cout << "{" << heap_[i].first << "," << heap_[i].second << "}\n";
+ }
+}
+
+bool tourney_tree::empty() const {
+ return event(0).time == max_time;
+}
+
+bool tourney_tree::empty(time_type t) const {
+ return event(0).time >= t;
+}
+
+postsynaptic_spike_event tourney_tree::head() const {
+ return event(0);
+}
+
+// Remove the smallest (most recent) event from the tree, then update the
+// tree so that head() returns the next event.
+void tourney_tree::pop() {
+ unsigned lane = id(0);
+ unsigned i = leaf(lane);
+ // draw the next event from the input lane
+ input_[lane]->pop();
+ // place event the leaf node for this lane
+ event(i) = input_[lane]->next();
+
+ // re-heapify the tree with a single walk from leaf to root
+ while ((i=parent(i))) {
+ merge_up(i);
+ }
+ merge_up(0); // handle the root
+}
+
+void tourney_tree::setup(unsigned i) {
+ if (is_leaf(i)) return;
+ setup(left(i));
+ setup(right(i));
+ merge_up(i);
+};
+
+// Update the value at node i of the tree to be the smallest
+// of its left and right children.
+// The result is undefined for leaf nodes.
+void tourney_tree::merge_up(unsigned i) {
+ const auto l = left(i);
+ const auto r = right(i);
+ heap_[i] = event(l)<event(r)? heap_[l]: heap_[r];
+}
+
+// The tree is stored using the standard heap indexing scheme.
+
+unsigned tourney_tree::parent(unsigned i) const {
+ return (i-1)>>1;
+}
+unsigned tourney_tree::left(unsigned i) const {
+ return (i<<1) + 1;
+}
+unsigned tourney_tree::right(unsigned i) const {
+ return left(i)+1;
+}
+unsigned tourney_tree::leaf(unsigned i) const {
+ return i+leaves_-1;
+}
+bool tourney_tree::is_leaf(unsigned i) const {
+ return i>=leaves_-1;
+}
+const unsigned& tourney_tree::id(unsigned i) const {
+ return heap_[i].first;
+}
+postsynaptic_spike_event& tourney_tree::event(unsigned i) {
+ return heap_[i].second;
+}
+const postsynaptic_spike_event& tourney_tree::event(unsigned i) const {
+ return heap_[i].second;
+}
+
+unsigned tourney_tree::next_power_2(unsigned x) const {
+ unsigned n = 1;
+ while (n<x) n<<=1;
+ return n;
+}
+
+} // namespace impl
+
+void merge_events(time_type t0, time_type t1,
+ const pse_vector& lc, pse_vector& events,
+ std::vector<event_generator_ptr>& generators,
+ pse_vector& lf)
+{
+ using std::distance;
+ using std::lower_bound;
+
+ // Sort events from the communicator in place.
+ util::sort(events);
+
+ // Clear lf to store merged list.
+ lf.clear();
+
+ // Merge the incoming event sequences into a single vector in lf
+ if (generators.size()) {
+ // Handle the case where the cell has event generators.
+ // This is performed in two steps:
+ // Step 1 : Use tournament tree to merge all events in lc, events and
+ // generators to be delivered in the time interval [tâ‚€, tâ‚).
+ // Step 2 : Use std::merge to append events in lc and events with
+ // delivery times in the interval [tâ‚, ∞).
+ EXPECTS(generators.size()>2u);
+
+ // Make an event generator with all the events in events.
+ generators[0] = make_event_generator<seq_generator<pse_vector>>(events);
+
+ // Make an event generator with all the events in lc with time >= t0
+ auto lc_it = lower_bound(lc.begin(), lc.end(), t0, event_time_less());
+ auto lc_range = util::make_range(lc_it, lc.end());
+ generators[1] = make_event_generator<seq_generator<decltype(lc_range)>>(lc_range);
+
+ // Perform k-way merge of all events in events, lc and the generators
+ // that are due to be delivered in the interval [tâ‚€, tâ‚)
+ impl::tourney_tree tree(generators);
+ while (!tree.empty(t1)) {
+ lf.push_back(tree.head());
+ tree.pop();
+ }
+
+ // Find first event in lc with delivery time >= t1
+ lc_it = lower_bound(lc.begin(), lc.end(), t1, event_time_less());
+ // Find first event in events with delivery time >= t1
+ auto ev_it = lower_bound(events.begin(), events.end(), t1, event_time_less());
+ const auto m = lf.size();
+ const auto n = m + distance(lc_it, lc.end()) + distance(ev_it, events.end());
+ lf.resize(n);
+ std::merge(ev_it, events.end(), lc_it, lc.end(), lf.begin()+m);
+ }
+ else {
+ // Handle the case where the cell has no event generators: only events
+ // in lc and lf with delivery times >= t0 must be merged, which can be
+ // handles with a single call to std::merge.
+ auto pos = std::lower_bound(lc.begin(), lc.end(), t0, event_time_less());
+ lf.resize(events.size()+distance(pos, lc.end()));
+ std::merge(events.begin(), events.end(), pos, lc.end(), lf.begin());
+ }
+}
+
+} // namespace arb
+
diff --git a/src/merge_events.hpp b/src/merge_events.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..bba08cbd54f651d26043080c3f27fb8a1e2d1c8e
--- /dev/null
+++ b/src/merge_events.hpp
@@ -0,0 +1,83 @@
+#pragma once
+
+#include <algorithm>
+#include <vector>
+
+#include <event_generator.hpp>
+#include <event_queue.hpp>
+#include <profiling/profiler.hpp>
+
+namespace arb {
+
+// merge_events generates a sorted list of postsynaptic events that are to be
+// delivered after the current epoch ends. It merges events from multiple
+// sources:
+// lc : the list of currently enqueued events
+// events : an unsorted list of events from the communicator
+// generators : a set of event_generators
+//
+// The time intervales are illustrated below, along with the range of times
+// for events in each of lc, events and generators
+// * t the start of the current epoch (not required to perform the merge).
+// * tâ‚€ the start of the next epoch
+// * tâ‚ the end of the next epoch
+//
+// t tâ‚€ tâ‚
+// |------|------|
+//
+// [----------------------] lc
+// [---------------] events
+// [------) generators
+//
+// The output list, stored in lf, will contain all the following:
+// * all events in events list
+// * events in lc with time >= tâ‚€
+// * events from each generator with time < tâ‚
+// All events in lc that are to be delivered before tâ‚€ are discared, along with
+// events from generators after tâ‚. The generators are left in a state where
+// the next event in the generator is the first event with deliver time >= tâ‚.
+void merge_events(time_type t0,
+ time_type t1,
+ const pse_vector& lc,
+ pse_vector& events,
+ std::vector<event_generator_ptr>& generators,
+ pse_vector& lf);
+
+namespace impl {
+ // The tournament tree is used internally by the merge_events method, and
+ // it is not intended for use elsewhere. It is exposed here for unit testing
+ // of its functionality.
+ class tourney_tree {
+ using key_val = std::pair<unsigned, postsynaptic_spike_event>;
+
+ public:
+ tourney_tree(std::vector<event_generator_ptr>& input);
+ bool empty() const;
+ bool empty(time_type t) const;
+ postsynaptic_spike_event head() const;
+ void pop();
+ void print() const;
+
+ private:
+ void setup(unsigned i);
+ void merge_up(unsigned i);
+ void update_lane(unsigned lane);
+ unsigned parent(unsigned i) const;
+ unsigned left(unsigned i) const;
+ unsigned right(unsigned i) const;
+ unsigned leaf(unsigned i) const;
+ bool is_leaf(unsigned i) const;
+ const unsigned& id(unsigned i) const;
+ postsynaptic_spike_event& event(unsigned i);
+ const postsynaptic_spike_event& event(unsigned i) const;
+ unsigned next_power_2(unsigned x) const;
+
+ std::vector<key_val> heap_;
+ const std::vector<event_generator_ptr>& input_;
+ unsigned leaves_;
+ unsigned nodes_;
+ unsigned n_lanes_;
+ };
+}
+
+} // namespace arb
diff --git a/src/model.cpp b/src/model.cpp
index 0495cbdab42ff536c94acad09b5420d796429df6..87144a3461e63f9076ed4ec7aa1ef4c0326fa555 100644
--- a/src/model.cpp
+++ b/src/model.cpp
@@ -5,6 +5,7 @@
#include <cell_group.hpp>
#include <cell_group_factory.hpp>
#include <domain_decomposition.hpp>
+#include <merge_events.hpp>
#include <model.hpp>
#include <recipe.hpp>
#include <util/filter.hpp>
@@ -14,15 +15,30 @@
namespace arb {
-void merge_events(time_type tfinal, const pse_vector& lc, pse_vector& events, pse_vector& lf);
-
model::model(const recipe& rec, const domain_decomposition& decomp):
communicator_(rec, decomp)
{
+ event_generators_.resize(communicator_.num_local_cells());
cell_local_size_type lidx = 0;
- for (auto i: util::make_span(0, decomp.groups.size())) {
- for (auto gid: decomp.groups[i].gids) {
+ const auto& grps = decomp.groups;
+ for (auto i: util::make_span(0, grps.size())) {
+ for (auto gid: grps[i].gids) {
+ // Store mapping of gid to local cell index.
gid_to_local_[gid] = lidx++;
+
+ // Set up the event generators for cell gid.
+ auto rec_gens = rec.event_generators(gid);
+ auto& gens = event_generators_[lidx];
+ if (rec_gens.size()) {
+ // Allocate two empty event generators that will be used to
+ // merge events from the communicator and those already queued
+ // for delivery in future epochs.
+ gens.reserve(2+rec_gens.size());
+ gens.resize(2);
+ for (auto& g: rec_gens) {
+ gens.push_back(std::move(g));
+ }
+ }
}
}
@@ -56,6 +72,14 @@ void model::reset() {
}
}
+ for (auto& lane: event_generators_) {
+ for (auto& gen: lane) {
+ if (gen) {
+ gen->reset();
+ }
+ }
+ }
+
communicator_.reset();
current_spikes().clear();
@@ -120,8 +144,10 @@ time_type model::run(time_type tfinal, time_type dt) {
const auto epid = epoch_.id;
merge_events(
epoch_.tfinal,
+ epoch_.tfinal+std::min(t_+t_interval, tfinal),
event_lanes(epid)[i],
events[i],
+ event_generators_[i],
event_lanes(epid+1)[i]);
});
PL(2);
@@ -246,35 +272,4 @@ void model::inject_events(const pse_vector& events) {
}
}
-// Merge events that are to be delivered from two lists into a sorted list.
-// Events are sorted by delivery time, then target, then weight.
-//
-// tfinal: The time at which the current epoch finishes. The output list, `lc`,
-// will contain all events with delivery times equal to or greater than
-// tfinal.
-// lc: Sorted set of events to be delivered before and after `tfinal`.
-// events: Unsorted list of events with delivery time greater than or equal to
-// tfinal. May be modified by the call.
-// lf: Will hold a list of all postsynaptic events in `events` and `lc` that
-// have delivery times greater than or equal to `tfinal`.
-void merge_events(time_type tfinal, const pse_vector& lc, pse_vector& events, pse_vector& lf) {
- // Merge the incoming events with events in lc that are not to be delivered
- // in this epoch, and store the result in lf.
-
- // STEP 1: sort events in place in events[l]
- PE("sort");
- util::sort(events);
- PL();
-
- // STEP 2: clear lf to store merged list
- lf.clear();
-
- // STEP 3: merge new events and future events from lc into lf
- PE("merge");
- auto pos = std::lower_bound(lc.begin(), lc.end(), tfinal, event_time_less());
- lf.resize(events.size()+std::distance(pos, lc.end()));
- std::merge(events.begin(), events.end(), pos, lc.end(), lf.begin());
- PL();
-}
-
} // namespace arb
diff --git a/src/model.hpp b/src/model.hpp
index b4abcdabe9576c17f316c797269ad5ba2143c8c4..268e31e99c594adacdb77209a26ec7168ec08b6a 100644
--- a/src/model.hpp
+++ b/src/model.hpp
@@ -70,6 +70,9 @@ private:
time_type t_ = 0.;
std::vector<cell_group_ptr> cell_groups_;
+ // one set of event_generators for each local cell
+ std::vector<std::vector<event_generator_ptr>> event_generators_;
+
using local_spike_store_type = thread_private_spike_store;
util::double_buffer<local_spike_store_type> local_spikes_;
diff --git a/src/recipe.hpp b/src/recipe.hpp
index 106fac4be60252528c0ec9f7da24ab5bc127ec8c..54ba235d3521a34ce46a7b840883c408384bd202 100644
--- a/src/recipe.hpp
+++ b/src/recipe.hpp
@@ -7,6 +7,7 @@
#include <cell.hpp>
#include <common_types.hpp>
+#include <event_generator.hpp>
#include <util/unique_any.hpp>
namespace arb {
@@ -63,6 +64,8 @@ public:
virtual cell_size_type num_targets(cell_gid_type) const = 0;
virtual cell_size_type num_probes(cell_gid_type) const = 0;
+ virtual std::vector<event_generator_ptr> event_generators(cell_gid_type) const = 0;
+
virtual std::vector<cell_connection> connections_on(cell_gid_type) const = 0;
virtual probe_info get_probe(cell_member_type probe_id) const = 0;
diff --git a/tests/global_communication/test_communicator.cpp b/tests/global_communication/test_communicator.cpp
index 147daead5124b2c7acc34022505f724be40efb38..21d7be0ed1429c7ea2c5716adcf9627f9580ee76 100644
--- a/tests/global_communication/test_communicator.cpp
+++ b/tests/global_communication/test_communicator.cpp
@@ -213,6 +213,10 @@ namespace {
1.0f)}; // delay
}
+ std::vector<event_generator_ptr> event_generators(cell_gid_type) const override {
+ return {};
+ }
+
probe_info get_probe(cell_member_type) const override {
throw std::logic_error("no probes");
}
@@ -282,6 +286,10 @@ namespace {
return cons;
}
+ std::vector<event_generator_ptr> event_generators(cell_gid_type) const override {
+ return {};
+ }
+
probe_info get_probe(cell_member_type) const override {
throw std::logic_error("no probes");
}
diff --git a/tests/global_communication/test_domain_decomposition.cpp b/tests/global_communication/test_domain_decomposition.cpp
index 745d1641fbde5be54e30c632e698067f5f2960be..119dded9f51c4209faaf3f5df44af00bb91278c7 100644
--- a/tests/global_communication/test_domain_decomposition.cpp
+++ b/tests/global_communication/test_domain_decomposition.cpp
@@ -54,6 +54,10 @@ namespace {
return {};
}
+ std::vector<event_generator_ptr> event_generators(cell_gid_type) const override {
+ return {};
+ }
+
probe_info get_probe(cell_member_type) const override {
throw std::logic_error("no probes");
}
diff --git a/tests/simple_recipes.hpp b/tests/simple_recipes.hpp
index 5417c0fbe9c12a04fb93a3a3e9d38dc83a980a1d..509f6dab90a7e4171b34a70fa47ff7a93e23238f 100644
--- a/tests/simple_recipes.hpp
+++ b/tests/simple_recipes.hpp
@@ -6,6 +6,7 @@
#include <vector>
#include <cell.hpp>
+#include <event_generator.hpp>
#include <recipe.hpp>
namespace arb {
@@ -19,6 +20,10 @@ public:
return probes_.count(i)? probes_.at(i).size(): 0;
}
+ std::vector<event_generator_ptr> event_generators(cell_gid_type) const override {
+ return {};
+ }
+
std::vector<cell_connection> connections_on(cell_gid_type) const override {
return {};
}
diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt
index 7410c7249a933bd798bd79e508f4086fc655fd6c..11ad7927747e4489907c11347f73e1d9de5e929a 100644
--- a/tests/unit/CMakeLists.txt
+++ b/tests/unit/CMakeLists.txt
@@ -44,8 +44,9 @@ set(TEST_SOURCES
test_domain_decomposition.cpp
test_dss_cell_group.cpp
test_either.cpp
- test_event_queue.cpp
test_event_binner.cpp
+ test_event_generators.cpp
+ test_event_queue.cpp
test_filter.cpp
test_fvm_multi.cpp
test_mc_cell_group.cpp
diff --git a/tests/unit/test_domain_decomposition.cpp b/tests/unit/test_domain_decomposition.cpp
index 0b988fe37622f3ec8e78fa6efa5d9a54f721c7c2..c6cf69fa02b1f831a3fb33a047fb91a7c0e77979 100644
--- a/tests/unit/test_domain_decomposition.cpp
+++ b/tests/unit/test_domain_decomposition.cpp
@@ -46,6 +46,10 @@ namespace {
return {};
}
+ std::vector<event_generator_ptr> event_generators(cell_gid_type) const override {
+ return {};
+ }
+
probe_info get_probe(cell_member_type) const override {
throw std::logic_error("no probes");
}
diff --git a/tests/unit/test_event_generators.cpp b/tests/unit/test_event_generators.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..76b6c6c30976e6ee0dd4bac2747bf5b0c123b546
--- /dev/null
+++ b/tests/unit/test_event_generators.cpp
@@ -0,0 +1,264 @@
+#include "../gtest.h"
+#include "common.hpp"
+
+#include <event_generator.hpp>
+#include <util/rangeutil.hpp>
+
+using namespace arb;
+using pse = postsynaptic_spike_event;
+
+namespace{
+ pse_vector draw(event_generator& gen, time_type t0, time_type t1) {
+ gen.reset();
+ gen.advance(t0);
+ pse_vector v;
+ while (gen.next().time<t1) {
+ v.push_back(gen.next());
+ gen.pop();
+ }
+ return v;
+ }
+}
+
+TEST(event_generators, vector_backed) {
+ std::vector<pse> in = {
+ {{0, 0}, 0.1, 1.0},
+ {{0, 0}, 1.0, 2.0},
+ {{0, 0}, 1.0, 3.0},
+ {{0, 0}, 1.5, 4.0},
+ {{0, 0}, 2.3, 5.0},
+ {{0, 0}, 3.0, 6.0},
+ {{0, 0}, 3.5, 7.0},
+ };
+
+ vector_backed_generator gen(in);
+
+ // Test pop, next and reset.
+ for (auto e: in) {
+ EXPECT_EQ(e, gen.next());
+ gen.pop();
+ }
+ gen.reset();
+ for (auto e: in) {
+ EXPECT_EQ(e, gen.next());
+ gen.pop();
+ }
+
+ // The loop above should have drained all events from gen, so we expect
+ // that the next() event will be the special terminal_pse event.
+ EXPECT_EQ(gen.next(), terminal_pse());
+}
+
+TEST(event_generators, regular) {
+ // make a regular generator that generates its first event at t=2ms and subsequent
+ // events regularly spaced 0.5 ms apart.
+ time_type t0 = 2.0;
+ time_type dt = 0.5;
+ cell_member_type target = {42, 3};
+ float weight = 3.14;
+
+ //regular_generator gen(t0, dt, target, weight);
+ regular_generator gen(target, weight, t0, dt);
+
+ // helper for building a set of
+ auto expected = [&] (std::vector<time_type> times) {
+ pse_vector events;
+ for (auto t: times) {
+ events.push_back({target, t, weight});
+ }
+ return events;
+ };
+
+ // Test pop, next and reset.
+ for (auto e: expected({2.0, 2.5, 3.0, 3.5, 4.0, 4.5})) {
+ EXPECT_EQ(e, gen.next());
+ gen.pop();
+ }
+ gen.reset();
+ for (auto e: expected({2.0, 2.5, 3.0, 3.5, 4.0, 4.5})) {
+ EXPECT_EQ(e, gen.next());
+ gen.pop();
+ }
+ gen.reset();
+
+ // Test advance
+ gen.advance(10.1);
+ EXPECT_EQ(gen.next().time, time_type(10.5));
+ gen.advance(12);
+ EXPECT_EQ(gen.next().time, time_type(12));
+}
+
+// Test for rounding problems with large time values and the regular generator
+TEST(event_generator, regular_rounding) {
+ // make a regular generator that generates its first event at t=2ms and subsequent
+ // events regularly spaced 0.5 ms apart.
+ time_type t0 = 2.0;
+ time_type dt = 0.5;
+ cell_member_type target = {42, 3};
+ float weight = 3.14;
+
+ // Test for rounding problems with large time values.
+ // To better understand why this is an issue, uncomment the following:
+ // float T = 1802667.0f, DT = 0.024999f;
+ // std::size_t N = std::floor(T/DT);
+ // std::cout << "T " << T << " DT " << DT << " N " << N
+ // << " T-N*DT " << T - (N*DT) << " P " << (T - (N*DT))/DT << "\n";
+ t0 = 1802667.0f;
+ dt = 0.024999f;
+ time_type int_len = 5*dt;
+ time_type t1 = t0 + int_len;
+ time_type t2 = t1 + int_len;
+ auto gen = regular_generator(target, weight, t0, dt);
+
+ // Take the interval I_a: t ∈ [t0, t2)
+ // And the two sub-interavls
+ // I_l: t ∈ [t0, t1)
+ // I_r: t ∈ [t1, t2)
+ // Such that I_a = I_l ∪ I_r.
+ // If we draw events from each interval then merge them, we expect same set
+ // of events as when we draw from that large interval.
+ pse_vector int_l = draw(gen, t0, t1);
+ pse_vector int_r = draw(gen, t1, t2);
+ pse_vector int_a = draw(gen, t0, t2);
+
+ pse_vector int_merged = int_l;
+ util::append(int_merged, int_r);
+
+ EXPECT_TRUE(int_l.front().time >= t0);
+ EXPECT_TRUE(int_l.back().time < t1);
+ EXPECT_TRUE(int_r.front().time >= t1);
+ EXPECT_TRUE(int_r.back().time < t2);
+ EXPECT_EQ(int_a, int_merged);
+ EXPECT_TRUE(std::is_sorted(int_a.begin(), int_a.end()));
+}
+
+TEST(event_generators, seq) {
+ std::vector<pse> in = {
+ {{0, 0}, 0.1, 1.0},
+ {{0, 0}, 1.0, 2.0},
+ {{0, 0}, 1.0, 3.0},
+ {{0, 0}, 1.5, 4.0},
+ {{0, 0}, 2.3, 5.0},
+ {{0, 0}, 3.0, 6.0},
+ {{0, 0}, 3.5, 7.0},
+ };
+
+ auto events = [&in] (int b, int e) {
+ return pse_vector(in.begin()+b, in.begin()+e);
+ };
+
+ seq_generator<pse_vector> gen(in);
+
+ // Test pop, next and reset.
+ for (auto e: in) {
+ EXPECT_EQ(e, gen.next());
+ gen.pop();
+ }
+ gen.reset();
+ for (auto e: in) {
+ EXPECT_EQ(e, gen.next());
+ gen.pop();
+ }
+ // The loop above should have drained all events from gen, so we expect
+ // that the next() event will be the special terminal_pse event.
+ EXPECT_EQ(gen.next(), terminal_pse());
+
+ gen.reset();
+
+ // Update of the input sequence, and run tests again to
+ // verify that results reflect the new set of input events.
+ in = {
+ {{0, 0}, 1.5, 4.0},
+ {{0, 0}, 2.3, 5.0},
+ {{0, 0}, 3.0, 6.0},
+ {{0, 0}, 3.5, 7.0},
+ };
+
+ // a range that includes all the events
+ EXPECT_EQ(in, draw(gen, 0, 4));
+
+
+ // a strict subset including the first event
+ EXPECT_EQ(events(0, 2), draw(gen, 0, 3));
+
+
+ // a strict subset including the last event
+ EXPECT_EQ(events(2, 4), draw(gen, 3.0, 5));
+
+
+ // subset that excludes first and last entries
+ EXPECT_EQ(events(1, 3), draw(gen, 2, 3.2));
+
+
+ // empty subset in the middle of range
+ EXPECT_EQ(pse_vector{}, draw(gen, 2, 2));
+
+
+ // empty subset before first event
+ EXPECT_EQ(pse_vector{}, draw(gen, 0, 0.05));
+
+
+ // empty subset after last event
+ EXPECT_EQ(pse_vector{}, draw(gen, 10, 11));
+
+}
+
+TEST(event_generators, poisson) {
+ std::mt19937_64 G;
+ using pgen = poisson_generator<std::mt19937_64>;
+
+ time_type t0 = 0;
+ time_type t1 = 10;
+ time_type lambda = 10; // expect 10 events per ms
+ cell_member_type target{4, 2};
+ float weight = 42;
+ pgen gen(target, weight, G, t0, lambda);
+
+ pse_vector int1;
+ while (gen.next().time<t1) {
+ int1.push_back(gen.next());
+ gen.pop();
+ }
+ // Test that the output is sorted
+ EXPECT_TRUE(std::is_sorted(int1.begin(), int1.end()));
+
+ // Reset and generate the same sequence of events
+ gen.reset();
+ pse_vector int2;
+ while (gen.next().time<t1) {
+ int2.push_back(gen.next());
+ gen.pop();
+ }
+
+ // Assert that the same sequence was generated
+ EXPECT_EQ(int1, int2);
+}
+
+// Test a poisson generator that has a tstop past which no events
+// should be generated.
+TEST(event_generators, poisson_terminates) {
+ std::mt19937_64 G;
+ using pgen = poisson_generator<std::mt19937_64>;
+
+ time_type t0 = 0;
+ time_type t1 = 10;
+ time_type t2 = 1e7; // pick a time far past the end of the interval [t0, t1)
+ time_type lambda = 10; // expect 10 events per ms
+ cell_member_type target{4, 2};
+ float weight = 42;
+ // construct generator with explicit end time t1
+ pgen gen(target, weight, G, t0, lambda, t1);
+
+ pse_vector events;
+ // pull events off the generator well past the end of the end time t1
+ while (gen.next().time<t2) {
+ events.push_back(gen.next());
+ gen.pop();
+ }
+
+ // the generator should be exhausted
+ EXPECT_EQ(gen.next(), terminal_pse());
+
+ // the last event should be less than the end time
+ EXPECT_TRUE(events.back().time<t1);
+}
diff --git a/tests/unit/test_merge_events.cpp b/tests/unit/test_merge_events.cpp
index f60fca0944128a20fe44bf5a48de623bc8a78fc0..481f8cf9379813192b812415b3e1e0297937f9a6 100644
--- a/tests/unit/test_merge_events.cpp
+++ b/tests/unit/test_merge_events.cpp
@@ -1,25 +1,12 @@
#include "../gtest.h"
#include <event_queue.hpp>
+#include <merge_events.hpp>
#include <model.hpp>
-namespace arb {
- // Declare prototype of the merge_events function, because it is only
- // defined in the TU of model.cpp
- void merge_events(time_type tfinal, const pse_vector& lc, pse_vector& events, pse_vector& lf);
-} // namespace arb
-
using namespace arb;
-std::ostream& operator<<(std::ostream& o, const pse_vector& events) {
- o << "{{"; for (auto e: events) o << " " << e;
- o << "}}";
- return o;
-}
-
-using pse = postsynaptic_spike_event;
-auto ev_bind = [] (const pse& e){ return std::tie(e.time, e.target, e.weight); };
-auto ev_less = [] (const pse& l, const pse& r){ return ev_bind(l)<ev_bind(r); };
+std::vector<event_generator_ptr> empty_gens;
// Test the trivial case of merging empty sets
TEST(merge_events, empty)
@@ -28,7 +15,7 @@ TEST(merge_events, empty)
pse_vector lc;
pse_vector lf;
- merge_events(0, lc, events, lf);
+ merge_events(0, max_time, lc, events, empty_gens, lf);
EXPECT_EQ(lf.size(), 0u);
}
@@ -43,7 +30,7 @@ TEST(merge_events, no_overlap)
{{0, 0}, 3, 3},
};
// Check that the inputs satisfy the precondition that lc is sorted.
- EXPECT_TRUE(std::is_sorted(lc.begin(), lc.end(), ev_less));
+ EXPECT_TRUE(std::is_sorted(lc.begin(), lc.end()));
// These events should be removed from lf by merge_events, and replaced
// with events to be delivered after t=10
@@ -60,7 +47,7 @@ TEST(merge_events, no_overlap)
{{0, 0}, 11, 1},
};
- merge_events(10, lc, events, lf);
+ merge_events(10, max_time, lc, events, empty_gens, lf);
pse_vector expected = {
{{8, 0}, 10, 4},
@@ -69,7 +56,7 @@ TEST(merge_events, no_overlap)
{{0, 0}, 12, 1},
};
- EXPECT_TRUE(std::is_sorted(lf.begin(), lf.end(), ev_less));
+ EXPECT_TRUE(std::is_sorted(lf.begin(), lf.end()));
EXPECT_EQ(expected, lf);
}
@@ -85,7 +72,7 @@ TEST(merge_events, overlap)
{{8, 0}, 10, 2},
{{0, 0}, 11, 3},
};
- EXPECT_TRUE(std::is_sorted(lc.begin(), lc.end(), ev_less));
+ EXPECT_TRUE(std::is_sorted(lc.begin(), lc.end()));
pse_vector lf;
@@ -98,7 +85,7 @@ TEST(merge_events, overlap)
{{7, 0}, 10, 8},
};
- merge_events(10, lc, events, lf);
+ merge_events(10, max_time, lc, events, empty_gens, lf);
pse_vector expected = {
{{7, 0}, 10, 8}, // from events
@@ -110,6 +97,147 @@ TEST(merge_events, overlap)
{{0, 0}, 12, 1}, // from events
};
- EXPECT_TRUE(std::is_sorted(lf.begin(), lf.end(), ev_less));
+ EXPECT_TRUE(std::is_sorted(lf.begin(), lf.end()));
+ EXPECT_EQ(expected, lf);
+}
+
+// Test the merge_events method with event generators.
+TEST(merge_events, X)
+{
+ const time_type t0 = 10;
+ const time_type t1 = 20;
+
+ pse_vector lc = {
+ {{0, 0}, 1, 1},
+ {{0, 0}, 5, 1},
+ // The current epoch ends at t=10, so all events from here down are expected in lf.
+ {{8, 0}, 10, 2},
+ {{0, 0}, 11, 3},
+ {{8, 0}, 20, 2},
+ {{0, 0}, 21, 3},
+ };
+ EXPECT_TRUE(std::is_sorted(lc.begin(), lc.end()));
+
+ pse_vector lf;
+
+ pse_vector events = {
+ {{0, 0}, 12, 1},
+ {{1, 0}, 15, 2},
+ {{2, 0}, 22, 3},
+ {{3, 0}, 26, 4},
+ };
+
+ std::vector<event_generator_ptr> generators(2);
+ generators.push_back(
+ make_event_generator<regular_generator>
+ (cell_member_type{4,2}, 42.f, t0, 5));
+
+ merge_events(t0, t1, lc, events, generators, lf);
+
+ pse_vector expected = {
+ {{4, 2}, 10, 42}, // from generator
+ {{8, 0}, 10, 2}, // from lc
+ {{0, 0}, 11, 3}, // from lc
+ {{0, 0}, 12, 1}, // from events
+ {{1, 0}, 15, 2}, // from events
+ {{4, 2}, 15, 42}, // from generator
+ {{8, 0}, 20, 2}, // from lc
+ {{0, 0}, 21, 3}, // from lc
+ {{2, 0}, 22, 3}, // from events
+ {{3, 0}, 26, 4}, // from events
+ };
+
+ EXPECT_TRUE(std::is_sorted(lf.begin(), lf.end()));
+ EXPECT_EQ(expected, lf);
+}
+
+// Test the tournament tree for merging two small sequences
+TEST(merge_events, tourney_seq)
+{
+ pse_vector g1 = {
+ {{0, 0}, 1, 1},
+ {{0, 0}, 2, 2},
+ {{0, 0}, 3, 3},
+ {{0, 0}, 4, 4},
+ {{0, 0}, 5, 5},
+ };
+
+ pse_vector g2 = {
+ {{0, 0}, 1.5, 1},
+ {{0, 0}, 2.5, 2},
+ {{0, 0}, 3.5, 3},
+ {{0, 0}, 4.5, 4},
+ {{0, 0}, 5.5, 5},
+ };
+
+ std::vector<event_generator_ptr> generators;
+ generators.push_back(make_event_generator<seq_generator<pse_vector>>(g1));
+ generators.push_back(make_event_generator<seq_generator<pse_vector>>(g2));
+ impl::tourney_tree tree(generators);
+
+ pse_vector lf;
+ while (!tree.empty()) {
+ lf.push_back(tree.head());
+ tree.pop();
+ }
+
+ EXPECT_TRUE(std::is_sorted(lf.begin(), lf.end()));
+ auto expected = g1;
+ util::append(expected, g2);
+ util::sort(expected);
+
EXPECT_EQ(expected, lf);
}
+
+// Test the tournament tree on a large set of Poisson generators.
+TEST(merge_events, tourney_poisson)
+{
+ using rndgen = std::mt19937_64;
+ // Number of poisson generators.
+ // Not a power of 2, so that there will be "null" leaf nodes in the
+ // tournament tree.
+ auto ngen = 100u;
+ time_type tfinal = 10;
+ time_type t0 = 0;
+ time_type lambda = 10; // expected: tfinal*lambda=1000 events per generator
+
+ std::vector<event_generator_ptr> generators;
+ for (auto i=0u; i<ngen; ++i) {
+ cell_member_type tgt{0, i};
+ float weight = i;
+ // the first and last generators have the same seed to test that sorting
+ // of events with the same time but different weights works properly.
+ rndgen G(i%(ngen-1));
+ generators.push_back(
+ make_event_generator<
+ poisson_generator<std::mt19937_64>>
+ (tgt, weight, G, t0, lambda));
+ }
+
+ // manually generate the expected output
+ pse_vector expected;
+ for (auto& gen: generators) {
+ // Push all events before tfinal in gen to the expected values.
+ while (gen->next().time<tfinal) {
+ expected.push_back(gen->next());
+ gen->pop();
+ }
+ // Reset the generator so that it is ready to generate the same
+ // events again for the tournament tree test.
+ gen->reset();
+ }
+ // Manually sort the expected events.
+ util::sort(expected);
+
+ // Generate output using tournament tree in lf.
+ impl::tourney_tree tree(generators);
+ pse_vector lf;
+ while (!tree.empty(tfinal)) {
+ lf.push_back(tree.head());
+ tree.pop();
+ }
+
+ // Test output of tournament tree.
+ EXPECT_TRUE(std::is_sorted(lf.begin(), lf.end()));
+ EXPECT_EQ(lf, expected);
+}