From 0d20df25e206ffa465a64ad3758d1ae24c5e8ec2 Mon Sep 17 00:00:00 2001
From: Sam Yates <yates@cscs.ch>
Date: Thu, 26 Jul 2018 20:04:00 +0200
Subject: [PATCH] Return view from schedule, replace time_seq. (#536)
Reduce redundant functionality across event_generator, time_seq
and schedule by providing a low-heap overhead interface to schedule
and using that for time sequences in event_generator and specialized
cell groups.
* Have schedule return pair of pointers as view to generated times.
* Fix missing DEBUG/TRACE functionality.
* Use rate instead of mean_dt for Poisson schedule.
* Move merge_events() functionality to simulator.cpp.
* Migrate event_generator to event span interface.
* Migrate tourney_tree to event span interface.
* Only invoke tourney_tree merge if generators have events in the epoch.
* Use schedule for times in event_generator implementations.
* Replace seq_generator with explicit_generator that keeps a copy of events.
* Replace vector_backed_generator and poisson_generator with schedule_generator.
* Replace time_seq uses with schedule.
* Add default empty schedule.
* Move rounding error test for regular time sequence into schedule test.
* Update sampling API documentation for schedule.
---
arbor/benchmark_cell_group.cpp | 9 +-
arbor/io/trace.hpp | 6 +-
arbor/mc_cell_group.cpp | 3 +-
arbor/merge_events.cpp | 131 ++++---------
arbor/merge_events.hpp | 44 +----
arbor/schedule.cpp | 27 +--
arbor/simulation.cpp | 98 +++++++---
arbor/spike_source_cell_group.cpp | 11 +-
arbor/spike_source_cell_group.hpp | 5 +-
arbor/util/range.hpp | 5 +
doc/sampling_api.rst | 12 +-
example/bench/recipe.cpp | 11 +-
example/brunel/brunel_miniapp.cpp | 4 +-
example/generators/event_gen.cpp | 13 +-
example/generators/readme.md | 12 +-
example/miniapp/miniapp_recipes.cpp | 6 +-
include/arbor/benchmark_cell.hpp | 4 +-
include/arbor/event_generator.hpp | 285 +++++++++++-----------------
include/arbor/schedule.hpp | 78 ++++++--
include/arbor/spike_source_cell.hpp | 4 +-
include/arbor/time_sequence.hpp | 248 ------------------------
test/unit/CMakeLists.txt | 1 -
test/unit/test_event_generators.cpp | 87 +++------
test/unit/test_lif_cell_group.cpp | 3 +-
test/unit/test_merge_events.cpp | 74 ++++++--
test/unit/test_schedule.cpp | 102 +++++++---
test/unit/test_spike_source.cpp | 68 ++++---
test/unit/test_time_seq.cpp | 175 -----------------
28 files changed, 542 insertions(+), 984 deletions(-)
delete mode 100644 include/arbor/time_sequence.hpp
delete mode 100644 test/unit/test_time_seq.cpp
diff --git a/arbor/benchmark_cell_group.cpp b/arbor/benchmark_cell_group.cpp
index 2da67ef0..4c24ff42 100644
--- a/arbor/benchmark_cell_group.cpp
+++ b/arbor/benchmark_cell_group.cpp
@@ -3,7 +3,7 @@
#include <arbor/benchmark_cell.hpp>
#include <arbor/recipe.hpp>
-#include <arbor/time_sequence.hpp>
+#include <arbor/schedule.hpp>
#include "cell_group.hpp"
#include "profile/profiler_macro.hpp"
@@ -50,7 +50,6 @@ void benchmark_cell_group::advance(epoch ep,
// Micro-seconds to advance in this epoch.
auto us = 1e3*(ep.tfinal-t_);
for (auto i: util::make_span(0, gids_.size())) {
- auto& tseq = cells_[i].time_sequence;
// Expected time to complete epoch in micro seconds.
const double duration_us = cells_[i].realtime_ratio*us;
const auto gid = gids_[i];
@@ -58,9 +57,9 @@ void benchmark_cell_group::advance(epoch ep,
// Start timer.
auto start = high_resolution_clock::now();
- while (tseq.front()<ep.tfinal) {
- spikes_.push_back({{gid, 0u}, tseq.front()});
- tseq.pop();
+ auto spike_times = util::make_range(cells_[i].time_sequence.events(t_, ep.tfinal));
+ for (auto t: spike_times) {
+ spikes_.push_back({{gid, 0u}, t});
}
// Wait until the expected time to advance has elapsed. Use a busy-wait
diff --git a/arbor/io/trace.hpp b/arbor/io/trace.hpp
index 15984715..3bb71054 100644
--- a/arbor/io/trace.hpp
+++ b/arbor/io/trace.hpp
@@ -35,6 +35,8 @@
namespace arb {
namespace impl {
+ inline void debug_emit_csv(std::ostream&) {}
+
template <typename Head, typename... Tail>
void debug_emit_csv(std::ostream& out, const Head& head, const Tail&... tail) {
out << head;
@@ -44,7 +46,9 @@ namespace impl {
debug_emit_csv(out, tail...);
}
- void debug_emit_trace_leader(std::ostream&, const char* file, int line, const char* vars);
+ inline void debug_emit_trace_leader(std::ostream& out, const char* file, int line, const char* vars) {
+ out << file << ':' << line << ": " << vars << ": ";
+ }
struct emit_nl_locked: public io::locked_ostream {
emit_nl_locked(std::streambuf* buf):
diff --git a/arbor/mc_cell_group.cpp b/arbor/mc_cell_group.cpp
index 54793ddc..f0680809 100644
--- a/arbor/mc_cell_group.cpp
+++ b/arbor/mc_cell_group.cpp
@@ -19,6 +19,7 @@
#include "util/filter.hpp"
#include "util/maputil.hpp"
#include "util/partition.hpp"
+#include "util/range.hpp"
#include "util/span.hpp"
namespace arb {
@@ -127,7 +128,7 @@ void mc_cell_group::advance(epoch ep, time_type dt, const event_lane_subrange& e
sample_size_type max_samples_per_call = 0;
for (auto& sa: sampler_map_) {
- auto sample_times = sa.sched.events(tstart, ep.tfinal);
+ auto sample_times = util::make_range(sa.sched.events(tstart, ep.tfinal));
if (sample_times.empty()) {
continue;
}
diff --git a/arbor/merge_events.cpp b/arbor/merge_events.cpp
index 096dd169..04355188 100644
--- a/arbor/merge_events.cpp
+++ b/arbor/merge_events.cpp
@@ -2,21 +2,24 @@
#include <set>
#include <vector>
+#include <arbor/assert.hpp>
#include <arbor/common_types.hpp>
-#include <arbor/domain_decomposition.hpp>
-#include <arbor/recipe.hpp>
+#include <arbor/math.hpp>
-#include "cell_group.hpp"
-#include "cell_group_factory.hpp"
+#include "io/trace.hpp"
#include "merge_events.hpp"
-#include "util/filter.hpp"
-#include "util/span.hpp"
#include "profile/profiler_macro.hpp"
namespace arb {
namespace impl {
+// A postsynaptic spike event that has delivery time set to
+// terminal_time, used as a sentinel in `tourney_tree`.
+
+static constexpr spike_event terminal_pse{cell_member_type{0,0}, terminal_time, 0};
+
+
// The tournament tree data structure is used to merge k sorted lists of events.
// See online for high-level information about tournament trees.
//
@@ -30,25 +33,27 @@ namespace impl {
// 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>& input):
+
+tourney_tree::tourney_tree(std::vector<event_span>& input):
input_(input),
n_lanes_(input_.size())
{
- // Must have at least 1 queue
- arb_assert(n_lanes_);
- // Maximum value in unsigned limits how many queues we can have
- arb_assert(n_lanes_<(1u<<(sizeof(unsigned)*8u-1u)));
+ // Must have at least 1 queue.
+ arb_assert(n_lanes_>=1u);
+
+ leaves_ = math::next_pow2(n_lanes_);
- leaves_ = next_power_2(n_lanes_);
- nodes_ = 2u*(leaves_-1u)+1u; // 2*l-1 with overflow protection
+ // Must be able to fit leaves in unsigned count.
+ arb_assert(leaves_>=n_lanes_);
+ nodes_ = 2*leaves_-1;
// 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].front()):
- key_val(i, make_terminal_pse()); // null leaf node
+ key_val(i, input[i].empty()? terminal_pse: input[i].front()):
+ key_val(i, terminal_pse); // null leaf node
}
// Walk the tree to initialize the non-leaf nodes
setup(0);
@@ -70,10 +75,6 @@ bool tourney_tree::empty() const {
return event(0).time == terminal_time;
}
-bool tourney_tree::empty(time_type t) const {
- return event(0).time >= t;
-}
-
spike_event tourney_tree::head() const {
return event(0);
}
@@ -83,10 +84,15 @@ spike_event tourney_tree::head() const {
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].front();
+ auto& in = input_[lane];
+
+ if (!in.empty()) {
+ ++in.left;
+ }
+
+ event(i) = in.empty()? terminal_pse: in.front();
// re-heapify the tree with a single walk from leaf to root
while ((i=parent(i))) {
@@ -138,83 +144,16 @@ const 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>& generators,
- pse_vector& lf)
-{
- using std::distance;
- using std::lower_bound;
-
- PE(communication_enqueue_sort);
-
- // Sort events from the communicator in place.
- util::sort(events);
- // Clear lf to store merged list.
- lf.clear();
-
- PL();
-
-
- // 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â‚, ∞).
- arb_assert(generators.size()>2u);
-
- PE(communication_enqueue_setup);
- // Make an event generator with all the events in events.
- generators[0] = 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] = seq_generator<decltype(lc_range)>(lc_range);
- PL();
-
- PE(communication_enqueue_tree);
- // 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();
- }
- PL();
-
- PE(communication_enqueue_merge);
- // 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);
- PL();
- }
- else {
- PE(communication_enqueue_merge);
- // 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());
- PL();
+void tree_merge_events(std::vector<event_span>& sources, pse_vector& out) {
+ PE(communication_enqueue_tree);
+ impl::tourney_tree tree(sources);
+ while (!tree.empty()) {
+ out.push_back(tree.head());
+ tree.pop();
}
+ PL();
}
} // namespace arb
diff --git a/arbor/merge_events.hpp b/arbor/merge_events.hpp
index eb058fec..6427d753 100644
--- a/arbor/merge_events.hpp
+++ b/arbor/merge_events.hpp
@@ -8,42 +8,15 @@
#include "event_queue.hpp"
#include "profile/profiler_macro.hpp"
+#include "util/range.hpp"
+
+// Merge a collection of sorted event sequences into a sorted output sequence.
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
-// pending_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
-// [---------------] pending_events
-// [------) generators
-//
-// The output list, stored in lf, will contain all the following:
-// * all events in pending_events
-// * 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& pending_events,
- std::vector<event_generator>& generators,
- pse_vector& lf);
+using event_span = util::range<const spike_event*>;
+
+void tree_merge_events(std::vector<event_span>& sources, pse_vector& out);
namespace impl {
// The tournament tree is used internally by the merge_events method, and
@@ -53,9 +26,8 @@ namespace impl {
using key_val = std::pair<unsigned, spike_event>;
public:
- tourney_tree(std::vector<event_generator>& input);
+ tourney_tree(std::vector<event_span>& input);
bool empty() const;
- bool empty(time_type t) const;
spike_event head() const;
void pop();
friend std::ostream& operator<<(std::ostream&, const tourney_tree&);
@@ -75,7 +47,7 @@ namespace impl {
unsigned next_power_2(unsigned x) const;
std::vector<key_val> heap_;
- std::vector<event_generator>& input_;
+ std::vector<event_span>& input_;
unsigned leaves_;
unsigned nodes_;
unsigned n_lanes_;
diff --git a/arbor/schedule.cpp b/arbor/schedule.cpp
index 2bfc246e..2a0d6bef 100644
--- a/arbor/schedule.cpp
+++ b/arbor/schedule.cpp
@@ -1,5 +1,6 @@
#include <algorithm>
#include <iterator>
+#include <numeric>
#include <utility>
#include <vector>
@@ -12,12 +13,14 @@ namespace arb {
// Regular schedule implementation.
-std::vector<time_type> regular_schedule_impl::events(time_type t0, time_type t1) {
- std::vector<time_type> ts;
+time_event_span regular_schedule_impl::events(time_type t0, time_type t1) {
+ times_.clear();
+
+ t0 = std::max(t0, t0_);
+ t1 = std::min(t1, t1_);
- t0 = t0<0? 0: t0;
if (t1>t0) {
- ts.reserve(1+std::size_t((t1-t0)*oodt_));
+ times_.reserve(1+std::size_t((t1-t0)*oodt_));
long long n = t0*oodt_;
time_type t = n*dt_;
@@ -27,22 +30,24 @@ std::vector<time_type> regular_schedule_impl::events(time_type t0, time_type t1)
}
while (t<t1) {
- ts.push_back(t);
+ times_.push_back(t);
t = (++n)*dt_;
}
}
- return ts;
+ return as_time_event_span(times_);
}
// Explicit schedule implementation.
-std::vector<time_type> explicit_schedule_impl::events(time_type t0, time_type t1) {
- auto lb = std::lower_bound(times_.begin()+start_index_, times_.end(), t0);
- auto ub = std::lower_bound(times_.begin()+start_index_, times_.end(), t1);
+time_event_span explicit_schedule_impl::events(time_type t0, time_type t1) {
+ time_event_span view = as_time_event_span(times_);
+
+ const time_type* lb = std::lower_bound(view.first+start_index_, view.second, t0);
+ const time_type* ub = std::lower_bound(lb, view.second, t1);
- start_index_ = std::distance(times_.begin(), ub);
- return std::vector<time_type>(lb, ub);
+ start_index_ = ub-view.first;
+ return {lb, ub};
}
} // namespace arb
diff --git a/arbor/simulation.cpp b/arbor/simulation.cpp
index dcef80c4..24587676 100644
--- a/arbor/simulation.cpp
+++ b/arbor/simulation.cpp
@@ -2,8 +2,9 @@
#include <set>
#include <vector>
-#include <arbor/recipe.hpp>
#include <arbor/domain_decomposition.hpp>
+#include <arbor/generic_event.hpp>
+#include <arbor/recipe.hpp>
#include <arbor/schedule.hpp>
#include <arbor/simulation.hpp>
@@ -95,8 +96,6 @@ private:
// Hash table for looking up the the local index of a cell with a given gid
std::unordered_map<cell_gid_type, cell_size_type> gid_to_local_;
- util::optional<cell_size_type> local_cell_index(cell_gid_type);
-
communicator communicator_;
task_system_handle task_system_;
@@ -150,18 +149,7 @@ simulation_state::simulation_state(
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));
- }
- }
+ event_generators_[lidx] = rec.event_generators(gid);
++lidx;
}
}
@@ -196,11 +184,10 @@ void simulation_state::reset() {
}
}
- // Reset all event generators, and advance to t_.
+ // Reset all event generators.
for (auto& lane: event_generators_) {
for (auto& gen: lane) {
gen.reset();
- gen.advance(t_);
}
}
@@ -294,6 +281,15 @@ time_type simulation_state::run(time_type tfinal, time_type dt) {
return t_;
}
+template <typename Seq, typename Value, typename Less = std::less<>>
+auto split_sorted_range(Seq&& seq, const Value& v, Less cmp = Less{}) {
+ auto canon = util::canonical_view(seq);
+ auto it = std::lower_bound(canon.begin(), canon.end(), v, cmp);
+ return std::make_pair(
+ util::make_range(seq.begin(), it),
+ util::make_range(it, seq.end()));
+}
+
// Populate the event lanes for epoch+1 (i.e event_lanes_[epoch+1)]
// Update each lane in parallel, if supported by the threading backend.
// On completion event_lanes[epoch+1] will contain sorted lists of events with
@@ -302,18 +298,72 @@ time_type simulation_state::run(time_type tfinal, time_type dt) {
// event_lanes[epoch]: take all events ≥ t_from
// event_generators : take all events < t_to
// pending_events : take all events
+
+// merge_cell_events() is a separate function for unit testing purposes.
+void merge_cell_events(
+ time_type t_from,
+ time_type t_to,
+ event_span old_events,
+ event_span pending,
+ std::vector<event_generator>& generators,
+ pse_vector& new_events)
+{
+ PE(communication_enqueue_setup);
+ new_events.clear();
+ old_events = split_sorted_range(old_events, t_from, event_time_less()).second;
+ PL();
+
+ if (!generators.empty()) {
+ PE(communication_enqueue_setup);
+ // Tree-merge events in [t_from, t_to) from old, pending and generator events.
+
+ std::vector<event_span> spanbuf;
+ spanbuf.reserve(2+generators.size());
+
+ auto old_split = split_sorted_range(old_events, t_to, event_time_less());
+ auto pending_split = split_sorted_range(pending, t_to, event_time_less());
+
+ spanbuf.push_back(old_split.first);
+ spanbuf.push_back(pending_split.first);
+
+ for (auto& g: generators) {
+ event_span evs = g.events(t_from, t_to);
+ if (!evs.empty()) {
+ spanbuf.push_back(evs);
+ }
+ }
+ PL();
+
+ PE(communication_enqueue_tree);
+ tree_merge_events(spanbuf, new_events);
+ PL();
+
+ old_events = old_split.second;
+ pending = pending_split.second;
+ }
+
+ // Merge (remaining) old and pending events.
+ PE(communication_enqueue_merge);
+ auto n = new_events.size();
+ new_events.resize(n+pending.size()+old_events.size());
+ std::merge(pending.begin(), pending.end(), old_events.begin(), old_events.end(), new_events.begin()+n);
+ PL();
+}
+
void simulation_state::setup_events(time_type t_from, time_type t_to, std::size_t epoch) {
const auto n = communicator_.num_local_cells();
threading::parallel_for::apply(0, n, task_system_.get(),
[&](cell_size_type i) {
- merge_events(
- t_from, t_to,
- event_lanes(epoch)[i], // in: the current event lane
- pending_events_[i], // in: events from the communicator
- event_generators_[i], // in: event generators for this lane
- event_lanes(epoch+1)[i]); // out: the event lane for the next epoch
+ PE(communication_enqueue_sort);
+ util::sort(pending_events_[i]);
+ PL();
+
+ event_span pending = util::range_pointer_view(pending_events_[i]);
+ event_span old_events = util::range_pointer_view(event_lanes(epoch)[i]);
+
+ merge_cell_events(t_from, t_to, old_events, pending, event_generators_[i], event_lanes(epoch+1)[i]);
pending_events_[i].clear();
- });
+ });
}
sampler_association_handle simulation_state::add_sampler(
diff --git a/arbor/spike_source_cell_group.cpp b/arbor/spike_source_cell_group.cpp
index aea7deaa..5b309e40 100644
--- a/arbor/spike_source_cell_group.cpp
+++ b/arbor/spike_source_cell_group.cpp
@@ -2,7 +2,7 @@
#include <arbor/recipe.hpp>
#include <arbor/spike_source_cell.hpp>
-#include <arbor/time_sequence.hpp>
+#include <arbor/schedule.hpp>
#include "cell_group.hpp"
#include "profile/profiler_macro.hpp"
@@ -34,14 +34,14 @@ void spike_source_cell_group::advance(epoch ep, time_type dt, const event_lane_s
PE(advance_sscell);
for (auto i: util::count_along(gids_)) {
- auto& tseq = time_sequences_[i];
const auto gid = gids_[i];
- while (tseq.front()<ep.tfinal) {
- spikes_.push_back({{gid, 0u}, tseq.front()});
- tseq.pop();
+ for (auto t: util::make_range(time_sequences_[i].events(t_, ep.tfinal))) {
+ spikes_.push_back({{gid, 0u}, t});
}
}
+ t_ = ep.tfinal;
+
PL();
};
@@ -49,6 +49,7 @@ void spike_source_cell_group::reset() {
for (auto& s: time_sequences_) {
s.reset();
}
+ t_ = 0;
clear_spikes();
}
diff --git a/arbor/spike_source_cell_group.hpp b/arbor/spike_source_cell_group.hpp
index 1fcf8a88..ff2d4b49 100644
--- a/arbor/spike_source_cell_group.hpp
+++ b/arbor/spike_source_cell_group.hpp
@@ -5,8 +5,8 @@
#include <arbor/common_types.hpp>
#include <arbor/recipe.hpp>
#include <arbor/sampling.hpp>
+#include <arbor/schedule.hpp>
#include <arbor/spike.hpp>
-#include <arbor/time_sequence.hpp>
#include "cell_group.hpp"
#include "epoch.hpp"
@@ -36,9 +36,10 @@ public:
void remove_all_samplers() override {}
private:
+ time_type t_ = 0;
std::vector<spike> spikes_;
std::vector<cell_gid_type> gids_;
- std::vector<time_seq> time_sequences_;
+ std::vector<schedule> time_sequences_;
};
} // namespace arb
diff --git a/arbor/util/range.hpp b/arbor/util/range.hpp
index 3f692016..154f8c35 100644
--- a/arbor/util/range.hpp
+++ b/arbor/util/range.hpp
@@ -60,6 +60,11 @@ struct range {
left(std::forward<U1>(l)), right(std::forward<U2>(r))
{}
+ template <typename U1, typename U2>
+ range(const std::pair<U1, U2>& p):
+ left(p.first), right(p.second)
+ {}
+
template <
typename U1,
typename U2,
diff --git a/doc/sampling_api.rst b/doc/sampling_api.rst
index a5abcfde..2ec63824 100644
--- a/doc/sampling_api.rst
+++ b/doc/sampling_api.rst
@@ -225,13 +225,21 @@ A ``schedule`` object has two methods:
void schedule::reset();
- std::vector<time_type> events(time_type t0, time_type t1)
+ time_event_span events(time_type t0, time_type t1)
-The ``events(t0, t1)`` method returns a vector of monotonically
+A ``time_event_span`` is a ``std::pair`` of pointers `const time_type*`,
+representing a view into an internally maintained collection of generated
+time values.
+
+The ``events(t0, t1)`` method returns a view of monotonically
increasing time values in the half-open interval ``[t0, t1)``.
Successive calls to ``events`` — without an intervening call to ``reset()``
— must request strictly subsequent intervals.
+The data represented by the returned ``time_event_span`` view is valid
+for the lifetime of the ``schedule`` object, and is invalidated by any
+subsequent call to ``reset()`` or ``events()``.
+
The ``reset()`` method resets the state such that events can be retrieved
from again from time zero. A schedule that is reset must then produce
the same sequence of time points, that is, it must exhibit repeatable
diff --git a/example/bench/recipe.cpp b/example/bench/recipe.cpp
index b8dd5f3f..aa9c94b8 100644
--- a/example/bench/recipe.cpp
+++ b/example/bench/recipe.cpp
@@ -2,7 +2,7 @@
#include <arbor/benchmark_cell.hpp>
#include <arbor/common_types.hpp>
-#include <arbor/time_sequence.hpp>
+#include <arbor/schedule.hpp>
#include "recipe.hpp"
@@ -15,19 +15,16 @@ cell_size_type bench_recipe::num_cells() const {
}
arb::util::unique_any bench_recipe::get_cell_description(cell_gid_type gid) const {
- using rng_type = std::mt19937_64;
+ std::mt19937_64 rng(gid);
arb::benchmark_cell cell;
cell.realtime_ratio = params_.cell.realtime_ratio;
// The time_sequence of the cell produces the series of time points at
- // which it will spike. We use a poisson_time_seq with a random sequence
+ // which it will spike. We use a poisson_schedule with a random sequence
// seeded with the gid. In this way, a cell's random stream depends only
// on its gid, and will hence give reproducable results when run with
// different MPI ranks and threads.
- cell.time_sequence =
- arb::poisson_time_seq<rng_type>(
- rng_type(gid), 0, 1e-3*params_.cell.spike_freq_hz);
-
+ cell.time_sequence = arb::poisson_schedule(1e-3*params_.cell.spike_freq_hz, rng);
return std::move(cell);
}
diff --git a/example/brunel/brunel_miniapp.cpp b/example/brunel/brunel_miniapp.cpp
index 9de099c0..19f3df93 100644
--- a/example/brunel/brunel_miniapp.cpp
+++ b/example/brunel/brunel_miniapp.cpp
@@ -126,12 +126,10 @@ public:
std::mt19937_64 G;
G.seed(gid + seed_);
- using pgen = poisson_generator<std::mt19937_64>;
-
time_type t0 = 0;
cell_member_type target{gid, 0};
- gens.emplace_back(pgen(target, weight_ext_, G, t0, lambda_));
+ gens.emplace_back(poisson_generator(target, weight_ext_, t0, lambda_, G));
return gens;
}
diff --git a/example/generators/event_gen.cpp b/example/generators/event_gen.cpp
index fd7423d4..a2e3d595 100644
--- a/example/generators/event_gen.cpp
+++ b/example/generators/event_gen.cpp
@@ -76,7 +76,6 @@ public:
arb_assert(gid==0); // There is only one cell in the model
using RNG = std::mt19937_64;
- using pgen = arb::poisson_generator<RNG>;
auto hz_to_freq = [](double hz) { return hz*1e-3; };
time_type t0 = 0;
@@ -92,15 +91,15 @@ public:
// Add excitatory generator
gens.push_back(
- pgen(cell_member_type{0,0}, // Target synapse (gid, local_id).
- w_e, // Weight of events to deliver
- RNG(29562872), // Random number generator to use
- t0, // Events start being delivered from this time
- lambda_e)); // Expected frequency (events per ms)
+ poisson_generator(cell_member_type{0,0}, // Target synapse (gid, local_id).
+ w_e, // Weight of events to deliver
+ t0, // Events start being delivered from this time
+ lambda_e, // Expected frequency (kHz)
+ RNG(29562872))); // Random number generator to use
// Add inhibitory generator
gens.emplace_back(
- pgen(cell_member_type{0,0}, w_i, RNG(86543891), t0, lambda_i));
+ poisson_generator(cell_member_type{0,0}, w_i, t0, lambda_i, RNG(86543891)));
return gens;
}
diff --git a/example/generators/readme.md b/example/generators/readme.md
index c9146c1a..f3f5c870 100644
--- a/example/generators/readme.md
+++ b/example/generators/readme.md
@@ -92,9 +92,6 @@ The implementation of this with hard-coded frequencies and weights is:
// The type of random number generator to use.
using RNG = std::mt19937_64;
- // The poisson_generator is templated on the random number generator type.
- using pgen = arb::poisson_generator<RNG>;
-
auto hz_to_freq = [](double hz) { return hz*1e-3; };
time_type t0 = 0;
@@ -109,15 +106,16 @@ The implementation of this with hard-coded frequencies and weights is:
// Add excitatory generator
gens.emplace_back(
- pgen(cell_member_type{0,0}, // Target synapse (gid, local_id).
+ arb::poisson_generator(
+ cell_member_type{0,0}, // Target synapse (gid, local_id).
w_e, // Weight of events to deliver
- RNG(29562872), // Random number generator to use
t0, // Events start being delivered from this time
- lambda_e)); // Expected frequency (events per ms)
+ lambda_e, // Expected frequency (events per ms)
+ RNG(29562872))); // Random number generator to use
// Add inhibitory generator
gens.emplace_back(
- pgen(cell_member_type{0,0}, w_i, RNG(86543891), t0, lambda_i));
+ arb::poisson_generator(cell_member_type{0,0}, w_i, t0, lambda_i, RNG(86543891)));
return gens;
}
diff --git a/example/miniapp/miniapp_recipes.cpp b/example/miniapp/miniapp_recipes.cpp
index 50838d1a..62d0ac99 100644
--- a/example/miniapp/miniapp_recipes.cpp
+++ b/example/miniapp/miniapp_recipes.cpp
@@ -7,8 +7,8 @@
#include <arbor/event_generator.hpp>
#include <arbor/mc_cell.hpp>
#include <arbor/morphology.hpp>
+#include <arbor/schedule.hpp>
#include <arbor/spike_source_cell.hpp>
-#include <arbor/time_sequence.hpp>
#include "io.hpp"
@@ -87,9 +87,9 @@ public:
}
util::unique_any get_cell_description(cell_gid_type i) const override {
- // The last 'cell' is a spike source cell.
+ // The last 'cell' is a spike source cell, producing one spike at t = 0.
if (i == ncell_) {
- return util::unique_any(spike_source_cell{regular_time_seq(0.0, 0.1, 0.1)});
+ return util::unique_any(spike_source_cell{explicit_schedule({0.})});
}
auto gen = std::mt19937(i); // TODO: replace this with hashing generator...
diff --git a/include/arbor/benchmark_cell.hpp b/include/arbor/benchmark_cell.hpp
index cbb32ba1..44795027 100644
--- a/include/arbor/benchmark_cell.hpp
+++ b/include/arbor/benchmark_cell.hpp
@@ -1,6 +1,6 @@
#pragma once
-#include <arbor/time_sequence.hpp>
+#include <arbor/schedule.hpp>
namespace arb {
@@ -9,7 +9,7 @@ namespace arb {
struct benchmark_cell {
// Describes the time points at which spikes are to be generated.
- time_seq time_sequence;
+ schedule time_sequence;
// Time taken in ms to advance the cell one ms of simulation time.
// If equal to 1, then a single cell can be advanced in realtime
diff --git a/include/arbor/event_generator.hpp b/include/arbor/event_generator.hpp
index 02ad77b8..247c53a3 100644
--- a/include/arbor/event_generator.hpp
+++ b/include/arbor/event_generator.hpp
@@ -9,48 +9,67 @@
#include <arbor/common_types.hpp>
#include <arbor/generic_event.hpp>
#include <arbor/spike_event.hpp>
-#include <arbor/time_sequence.hpp>
+#include <arbor/schedule.hpp>
namespace arb {
-// Generate a postsynaptic spike event that has delivery time set to
-// terminal_time. Such events are used as sentinels, to indicate the
-// end of a sequence.
-inline constexpr
-spike_event make_terminal_pse() {
- return spike_event{cell_member_type{0,0}, terminal_time, 0};
-}
-
-inline
-bool is_terminal_pse(const spike_event& e) {
- return e.time==terminal_time;
-}
+// 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.
+//
+// An `event_generator` supports two operations:
+//
+// `void event_generator::reset()`
+//
+// Reset generator state.
+//
+// `event_seq event_generator::events(time_type to, time_type from)`
+//
+// Provide a non-owning view on to the events in the time interval
+// [to, from).
+//
+// The `event_seq` type is a pair of `spike_event` pointers that
+// provide a view onto an internally-maintained contiguous sequence
+// of generated spike event objects. This view is valid only for
+// the lifetime of the generator, and is invalidated upon a call
+// to `reset` or another call to `events`.
+//
+// Calls to the `events` method must be monotonic in time: without an
+// intervening call to `reset`, two successive calls `events(t0, t1)`
+// and `events(t2, t3)` to the same event generator must satisfy
+// 0 ≤ t0 ≤ t1 ≤ t2 ≤ t3.
+//
+// `event_generator` objects have value semantics, and use type erasure
+// to wrap implementation details. An `event_generator` can be constructed
+// from an onbject of an implementation class Impl that is copy-constructible
+// and otherwise provides `reset` and `events` methods following the
+// API described above.
+//
+// Some pre-defined event generators are included:
+// - `empty_generator`: produces no events
+// - `schedule_generator`: events to a fixed target according to a time schedule
+
+using event_seq = std::pair<const spike_event*, const spike_event*>;
// The simplest possible generator that generates no events.
// Declared ahead of event_generator so that it can be used as the default
// generator.
struct empty_generator {
- spike_event front() {
- return spike_event{cell_member_type{0,0}, terminal_time, 0};
- }
- void pop() {}
void reset() {}
- void advance(time_type t) {};
+ event_seq events(time_type, time_type) {
+ return {&no_event, &no_event};
+ }
+
+private:
+ static spike_event no_event;
};
-// 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.
class event_generator {
public:
- //
- // copy, move and constructor interface
- //
-
event_generator(): event_generator(empty_generator()) {}
template <typename Impl>
@@ -70,40 +89,18 @@ public:
return *this;
}
- //
- // event generator interface
- //
-
- // Get the current event in the stream.
- // Does not modify the state of the stream, i.e. multiple calls to
- // front() will return the same event in the absence of calls to pop(),
- // advance() or reset().
- spike_event front() {
- return impl_->front();
- }
-
- // Move the generator to the next event in the stream.
- void pop() {
- impl_->pop();
- }
-
- // Reset the generator to the same state that it had on construction.
void reset() {
impl_->reset();
}
- // Update state of the generator such that the event returned by front() is
- // the first event with delivery time >= t.
- void advance(time_type t) {
- return impl_->advance(t);
+ event_seq events(time_type t0, time_type t1) {
+ return impl_->events(t0, t1);
}
private:
struct interface {
- virtual spike_event front() = 0;
- virtual void pop() = 0;
- virtual void advance(time_type t) = 0;
virtual void reset() = 0;
+ virtual event_seq events(time_type, time_type) = 0;
virtual std::unique_ptr<interface> clone() = 0;
virtual ~interface() {}
};
@@ -115,16 +112,8 @@ private:
explicit wrap(const Impl& impl): wrapped(impl) {}
explicit wrap(Impl&& impl): wrapped(std::move(impl)) {}
- spike_event front() override {
- return wrapped.front();
- }
-
- void pop() override {
- return wrapped.pop();
- }
-
- void advance(time_type t) override {
- return wrapped.advance(t);
+ event_seq events(time_type t0, time_type t1) override {
+ return wrapped.events(t0, t1);
}
void reset() override {
@@ -139,153 +128,93 @@ private:
};
};
-// Generator that feeds events that are specified with a vector.
-// Makes a copy of the input sequence of events.
-struct vector_backed_generator {
- using pse = spike_event;
- vector_backed_generator(cell_member_type target, float weight, std::vector<time_type> samples):
- target_(target),
- weight_(weight),
- tseq_(std::move(samples))
- {}
- spike_event front() {
- return spike_event{target_, tseq_.front(), weight_};
- }
+// Generate events with a fixed target and weight according to
+// a provided time schedule.
- void pop() {
- tseq_.pop();
- }
+struct schedule_generator {
+ schedule_generator(cell_member_type target, float weight, schedule sched):
+ target_(target), weight_(weight), sched_(std::move(sched))
+ {}
void reset() {
- tseq_.reset();
+ sched_.reset();
}
- void advance(time_type t) {
- tseq_.advance(t);
- }
+ event_seq events(time_type t0, time_type t1) {
+ auto ts = sched_.events(t0, t1);
-private:
- cell_member_type target_;
- float weight_;
- vector_time_seq tseq_;
-};
+ events_.clear();
+ events_.reserve(ts.second-ts.first);
-// 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 {
- using pse = spike_event;
- seq_generator(Seq& events):
- events_(events),
- it_(std::begin(events_))
- {
- arb_assert(std::is_sorted(events_.begin(), events_.end()));
- }
-
- spike_event front() {
- return it_==events_.end()? make_terminal_pse(): *it_;
- }
-
- void pop() {
- if (it_!=events_.end()) {
- ++it_;
+ for (auto i = ts.first; i!=ts.second; ++i) {
+ events_.push_back(spike_event{target_, *i, weight_});
}
- }
- void reset() {
- it_ = events_.begin();
- }
-
- void advance(time_type t) {
- it_ = std::lower_bound(events_.begin(), events_.end(), t, event_time_less());
+ return {events_.data(), events_.data()+events_.size()};
}
private:
- const Seq& events_;
- typename Seq::const_iterator it_;
+ pse_vector events_;
+ cell_member_type target_;
+ float weight_;
+ schedule sched_;
};
-// 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 {
- using pse = spike_event;
-
- regular_generator(cell_member_type target,
- float weight,
- time_type tstart,
- time_type dt,
- time_type tstop=terminal_time):
- target_(target),
- weight_(weight),
- tseq_(tstart, dt, tstop)
- {}
+// Convenience routines for making schedule_generator:
- spike_event front() {
- return spike_event{target_, tseq_.front(), weight_};
- }
+inline event_generator regular_generator(
+ cell_member_type target, float weight, time_type tstart, time_type dt, time_type tstop=terminal_time)
+{
+ return schedule_generator(target, weight, regular_schedule(tstart, dt, tstop));
+}
- void pop() {
- tseq_.pop();
- }
+template <typename RNG>
+inline event_generator poisson_generator(
+ cell_member_type target, float weight, time_type tstart, time_type rate_kHz, const RNG& rng)
+{
+ return schedule_generator(target, weight, poisson_schedule(tstart, rate_kHz, rng));
+}
- void advance(time_type t0) {
- tseq_.advance(t0);
- }
- void reset() {
- tseq_.reset();
- }
+// Generate events from a predefined sorted event sequence.
-private:
- cell_member_type target_;
- float weight_;
- regular_time_seq tseq_;
-};
+struct explicit_generator {
+ explicit_generator() = default;
+ explicit_generator(const explicit_generator&) = default;
+ explicit_generator(explicit_generator&&) = default;
-// 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 {
- using pse = spike_event;
-
- poisson_generator(cell_member_type target,
- float weight,
- RandomNumberEngine rng,
- time_type tstart,
- time_type rate_per_ms,
- time_type tstop=terminal_time):
- target_(target),
- weight_(weight),
- tseq_(std::move(rng), tstart, rate_per_ms, tstop)
+ template <typename Seq>
+ explicit_generator(const Seq& events):
+ start_index_(0)
{
- reset();
- }
+ using std::begin;
+ using std::end;
- spike_event front() {
- return spike_event{target_, tseq_.front(), weight_};
+ events_ = pse_vector(begin(events), end(events));
+ arb_assert(std::is_sorted(events_.begin(), events_.end()));
}
- void pop() {
- tseq_.pop();
+ void reset() {
+ start_index_ = 0;
}
- void advance(time_type t0) {
- tseq_.advance(t0);
- }
+ event_seq events(time_type t0, time_type t1) {
+ const spike_event* lb = events_.data()+start_index_;
+ const spike_event* ub = events_.data()+events_.size();
- void reset() {
- tseq_.reset();
+ lb = std::lower_bound(lb, ub, t0, event_time_less{});
+ ub = std::lower_bound(lb, ub, t1, event_time_less{});
+
+ start_index_ = ub-events_.data();
+ return {lb, ub};
}
private:
- const cell_member_type target_;
- const float weight_;
- poisson_time_seq<RandomNumberEngine> tseq_;
+ pse_vector events_;
+ std::size_t start_index_ = 0;
};
+
} // namespace arb
diff --git a/include/arbor/schedule.hpp b/include/arbor/schedule.hpp
index 7ef42afc..4a9213e0 100644
--- a/include/arbor/schedule.hpp
+++ b/include/arbor/schedule.hpp
@@ -4,6 +4,7 @@
#include <iterator>
#include <memory>
#include <random>
+#include <utility>
#include <vector>
#include <arbor/assert.hpp>
@@ -14,6 +15,12 @@
namespace arb {
+using time_event_span = std::pair<const time_type*, const time_type*>;
+
+inline time_event_span as_time_event_span(const std::vector<time_type>& v) {
+ return {&v[0], &v[0]+v.size()};
+}
+
// A schedule describes a sequence of time values used for sampling. Schedules
// are queried monotonically in time: if two method calls `events(t0, t1)`
// and `events(t2, t3)` are made without an intervening call to `reset()`,
@@ -21,6 +28,8 @@ namespace arb {
class schedule {
public:
+ schedule();
+
template <typename Impl>
explicit schedule(const Impl& impl):
impl_(new wrap<Impl>(impl)) {}
@@ -40,7 +49,7 @@ public:
return *this;
}
- std::vector<time_type> events(time_type t0, time_type t1) {
+ time_event_span events(time_type t0, time_type t1) {
return impl_->events(t0, t1);
}
@@ -48,7 +57,7 @@ public:
private:
struct interface {
- virtual std::vector<time_type> events(time_type t0, time_type t1) = 0;
+ virtual time_event_span events(time_type t0, time_type t1) = 0;
virtual void reset() = 0;
virtual std::unique_ptr<interface> clone() = 0;
virtual ~interface() {}
@@ -61,7 +70,7 @@ private:
explicit wrap(const Impl& impl): wrapped(impl) {}
explicit wrap(Impl&& impl): wrapped(std::move(impl)) {}
- virtual std::vector<time_type> events(time_type t0, time_type t1) {
+ virtual time_event_span events(time_type t0, time_type t1) {
return wrapped.events(t0, t1);
}
@@ -77,27 +86,51 @@ private:
};
};
+// Default schedule is empty.
+
+class empty_schedule {
+public:
+ void reset() {}
+ time_event_span events(time_type t0, time_type t1) {
+ static time_type no_time;
+ return {&no_time, &no_time};
+ }
+};
+
+inline schedule::schedule(): schedule(empty_schedule{}) {}
// Common schedules
-// Schedule at k·dt for integral k≥0.
+// Schedule at k·dt for integral k≥0 within the interval [t0, t1).
class regular_schedule_impl {
public:
- explicit regular_schedule_impl(time_type dt):
- dt_(dt), oodt_(1./dt) {};
+ explicit regular_schedule_impl(time_type t0, time_type dt, time_type t1):
+ t0_(t0), t1_(t1), dt_(dt), oodt_(1./dt)
+ {
+ if (t0_<0) t0_ = 0;
+ };
void reset() {}
- std::vector<time_type> events(time_type t0, time_type t1);
+ time_event_span events(time_type t0, time_type t1);
private:
- time_type dt_;
+ time_type t0_, t1_, dt_;
time_type oodt_;
+
+ std::vector<time_type> times_;
};
+inline schedule regular_schedule(
+ time_type t0, time_type dt, time_type t1 = std::numeric_limits<time_type>::max())
+{
+ return schedule(regular_schedule_impl(t0, dt, t1));
+}
+
inline schedule regular_schedule(time_type dt) {
- return schedule(regular_schedule_impl(dt));
+ return regular_schedule(0, dt);
}
+
// Schedule at times given explicitly via a provided sorted sequence.
class explicit_schedule_impl {
public:
@@ -119,7 +152,7 @@ public:
start_index_ = 0;
}
- std::vector<time_type> events(time_type t0, time_type t1);
+ time_event_span events(time_type t0, time_type t1);
private:
std::ptrdiff_t start_index_;
@@ -131,13 +164,17 @@ inline schedule explicit_schedule(const Seq& seq) {
return schedule(explicit_schedule_impl(seq));
}
+inline schedule explicit_schedule(const std::initializer_list<time_type>& seq) {
+ return schedule(explicit_schedule_impl(seq));
+}
+
// Schedule at Poisson point process with rate 1/mean_dt,
// restricted to non-negative times.
template <typename RandomNumberEngine>
class poisson_schedule_impl {
public:
- poisson_schedule_impl(time_type tstart, time_type mean_dt, const RandomNumberEngine& rng):
- tstart_(tstart), exp_(1./mean_dt), rng_(rng), reset_state_(rng), next_(tstart)
+ poisson_schedule_impl(time_type tstart, time_type rate_kHz, const RandomNumberEngine& rng):
+ tstart_(tstart), exp_(rate_kHz), rng_(rng), reset_state_(rng), next_(tstart)
{
arb_assert(tstart_>=0);
step();
@@ -149,19 +186,19 @@ public:
step();
}
- std::vector<time_type> events(time_type t0, time_type t1) {
- std::vector<time_type> ts;
+ time_event_span events(time_type t0, time_type t1) {
+ times_.clear();
while (next_<t0) {
step();
}
while (next_<t1) {
- ts.push_back(next_);
+ times_.push_back(next_);
step();
}
- return ts;
+ return as_time_event_span(times_);
}
private:
@@ -174,16 +211,17 @@ private:
RandomNumberEngine rng_;
RandomNumberEngine reset_state_;
time_type next_;
+ std::vector<time_type> times_;
};
template <typename RandomNumberEngine>
-inline schedule poisson_schedule(time_type mean_dt, const RandomNumberEngine& rng) {
- return schedule(poisson_schedule_impl<RandomNumberEngine>(0., mean_dt, rng));
+inline schedule poisson_schedule(time_type rate_kHz, const RandomNumberEngine& rng) {
+ return schedule(poisson_schedule_impl<RandomNumberEngine>(0., rate_kHz, rng));
}
template <typename RandomNumberEngine>
-inline schedule poisson_schedule(time_type tstart, time_type mean_dt, const RandomNumberEngine& rng) {
- return schedule(poisson_schedule_impl<RandomNumberEngine>(tstart, mean_dt, rng));
+inline schedule poisson_schedule(time_type tstart, time_type rate_kHz, const RandomNumberEngine& rng) {
+ return schedule(poisson_schedule_impl<RandomNumberEngine>(tstart, rate_kHz, rng));
}
} // namespace arb
diff --git a/include/arbor/spike_source_cell.hpp b/include/arbor/spike_source_cell.hpp
index 74b66e5b..0e9bbf17 100644
--- a/include/arbor/spike_source_cell.hpp
+++ b/include/arbor/spike_source_cell.hpp
@@ -1,6 +1,6 @@
#pragma once
-#include <arbor/time_sequence.hpp>
+#include <arbor/schedule.hpp>
namespace arb {
@@ -8,7 +8,7 @@ namespace arb {
// recipe::cell_kind(gid) returning cell_kind::spike_source
struct spike_source_cell {
- time_seq seq;
+ schedule seq;
};
} // namespace arb
diff --git a/include/arbor/time_sequence.hpp b/include/arbor/time_sequence.hpp
deleted file mode 100644
index 9d2d0029..00000000
--- a/include/arbor/time_sequence.hpp
+++ /dev/null
@@ -1,248 +0,0 @@
-#pragma once
-
-#include <algorithm>
-#include <memory>
-#include <random>
-#include <type_traits>
-
-#include <arbor/common_types.hpp>
-
-namespace arb {
-
-struct empty_time_seq {
- time_type front() { return terminal_time; }
- void pop() {}
- void reset() {}
- void advance(time_type t) {};
-};
-
-// An time_sequence generates a sequence of time points.
-// The sequence of times is always monotonically increasing, i.e. each time is
-// not earlier than the time that proceded it.
-class time_seq {
-public:
- //
- // copy, move and constructor interface
- //
-
- time_seq(): time_seq(empty_time_seq()) {}
-
- template <
- typename Impl,
- typename = std::enable_if_t<
- !std::is_same<std::decay_t<Impl>, time_seq>::value>
- >
- time_seq(Impl&& impl):
- impl_(new wrap<Impl>(std::forward<Impl>(impl)))
- {}
-
- time_seq(time_seq&& other) = default;
- time_seq& operator=(time_seq&& other) = default;
-
- time_seq(const time_seq& other):
- impl_(other.impl_->clone())
- {}
-
- time_seq& operator=(const time_seq& other) {
- impl_ = other.impl_->clone();
- return *this;
- }
-
- //
- // time sequence interface
- //
-
- // Get the current time in the stream.
- // Does not modify the state of the stream, i.e. multiple calls to
- // front() will return the same time in the absence of calls to pop(),
- // advance() or reset().
- time_type front() {
- return impl_->front();
- }
-
- // Move the generator to the front time in the stream.
- void pop() {
- impl_->pop();
- }
-
- // Reset the generator to the same state that it had on construction.
- void reset() {
- impl_->reset();
- }
-
- // Update state of the generator such that the time returned by front() is
- // the first time with delivery time >= t.
- void advance(time_type t) {
- return impl_->advance(t);
- }
-
-private:
- struct interface {
- virtual time_type front() = 0;
- virtual void pop() = 0;
- virtual void advance(time_type t) = 0;
- virtual void reset() = 0;
- virtual std::unique_ptr<interface> clone() = 0;
- virtual ~interface() {}
- };
-
- std::unique_ptr<interface> impl_;
-
- template <typename Impl>
- struct wrap: interface {
- explicit wrap(const Impl& impl): wrapped(impl) {}
- explicit wrap(Impl&& impl): wrapped(std::move(impl)) {}
-
- time_type front() override {
- return wrapped.front();
- }
- void pop() override {
- return wrapped.pop();
- }
- void advance(time_type t) override {
- return wrapped.advance(t);
- }
- void reset() override {
- wrapped.reset();
- }
- std::unique_ptr<interface> clone() override {
- return std::unique_ptr<interface>(new wrap<Impl>(wrapped));
- }
-
- Impl wrapped;
- };
-};
-
-// Sequence of time points prescribed by a vector
-struct vector_time_seq {
- vector_time_seq(std::vector<time_type> seq):
- seq_(std::move(seq))
- {
- // Ensure that the time values are sorted.
- if (!std::is_sorted(seq_.begin(), seq_.end())) {
- std::sort(seq_.begin(), seq_.end());
- }
- reset();
- }
-
- time_type front() {
- return it_==seq_.end()? terminal_time: *it_;
- }
-
- void pop() {
- if (it_!=seq_.end()) {
- ++it_;
- }
- }
-
- void advance(time_type t0) {
- it_ = std::lower_bound(seq_.begin(), seq_.end(), t0);
- }
-
- void reset() {
- it_ = seq_.begin();
- }
-
-private:
- std::vector<time_type> seq_;
- std::vector<time_type>::const_iterator it_;
-};
-
-// Generates a set of regularly spaced time samples.
-// t=t_start+n*dt, ∀ t ∈ [t_start, t_stop)
-struct regular_time_seq {
- regular_time_seq(time_type tstart,
- time_type dt,
- time_type tstop=terminal_time):
- step_(0),
- t_start_(tstart),
- dt_(dt),
- t_stop_(tstop)
- {}
-
- time_type front() {
- const auto t = time();
- return t<t_stop_? t: terminal_time;
- }
-
- void pop() {
- ++step_;
- }
-
- void advance(time_type t0) {
- 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() {
- step_ = 0;
- }
-
-private:
- time_type time() const {
- return t_start_ + step_*dt_;
- }
-
- std::size_t step_;
- time_type t_start_;
- time_type dt_;
- time_type t_stop_;
-};
-
-// Generates a stream of time points described by a Poisson point process
-// with rate_per_ms samples per ms.
-template <typename RandomNumberEngine>
-struct poisson_time_seq {
- poisson_time_seq(RandomNumberEngine rng,
- time_type tstart,
- time_type rate_per_ms,
- time_type tstop=terminal_time):
- exp_(rate_per_ms),
- reset_state_(std::move(rng)),
- t_start_(tstart),
- t_stop_(tstop)
- {
- reset();
- }
-
- time_type front() {
- return next_<t_stop_? next_: terminal_time;
- }
-
- void pop() {
- next_ += exp_(rng_);
- }
-
- void advance(time_type t0) {
- while (next_<t0) {
- pop();
- }
- }
-
- void reset() {
- rng_ = reset_state_;
- next_ = t_start_;
- pop();
- }
-
-private:
- std::exponential_distribution<time_type> exp_;
- RandomNumberEngine rng_;
- const RandomNumberEngine reset_state_;
- const time_type t_start_;
- const time_type t_stop_;
- time_type next_;
-};
-
-} // namespace arb
-
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index b2e99d0d..42205d7d 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -71,7 +71,6 @@ set(unit_sources
test_strprintf.cpp
test_swcio.cpp
test_synapses.cpp
- test_time_seq.cpp
test_thread.cpp
test_tree.cpp
test_transform.cpp
diff --git a/test/unit/test_event_generators.cpp b/test/unit/test_event_generators.cpp
index 6004986c..d79eb423 100644
--- a/test/unit/test_event_generators.cpp
+++ b/test/unit/test_event_generators.cpp
@@ -10,15 +10,8 @@
using namespace arb;
namespace{
- pse_vector draw(event_generator& gen, time_type t0, time_type t1) {
- gen.reset();
- gen.advance(t0);
- pse_vector v;
- while (gen.front().time<t1) {
- v.push_back(gen.front());
- gen.pop();
- }
- return v;
+ pse_vector as_vector(event_seq s) {
+ return pse_vector(s.first, s.second);
}
}
@@ -30,8 +23,7 @@ TEST(event_generators, regular) {
cell_member_type target = {42, 3};
float weight = 3.14;
- //regular_generator gen(t0, dt, target, weight);
- regular_generator gen(target, weight, t0, dt);
+ event_generator gen = regular_generator(target, weight, t0, dt);
// helper for building a set of
auto expected = [&] (std::vector<time_type> times) {
@@ -42,23 +34,15 @@ TEST(event_generators, regular) {
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.front());
- gen.pop();
- }
- gen.reset();
- for (auto e: expected({2.0, 2.5, 3.0, 3.5, 4.0, 4.5})) {
- EXPECT_EQ(e, gen.front());
- gen.pop();
- }
+ EXPECT_EQ(expected({2.0, 2.5, 3.0, 3.5, 4.0, 4.5}), as_vector(gen.events(0, 5)));
+
+ // Test reset and re-generate.
gen.reset();
+ EXPECT_EQ(expected({2.0, 2.5, 3.0, 3.5, 4.0, 4.5}), as_vector(gen.events(0, 5)));
- // Test advance
- gen.advance(10.1);
- EXPECT_EQ(gen.front().time, time_type(10.5));
- gen.advance(12);
- EXPECT_EQ(gen.front().time, time_type(12));
+ // Test later intervals.
+ EXPECT_EQ(expected({10.5}), as_vector(gen.events(10.1, 10.7)));
+ EXPECT_EQ(expected({12, 12.5}), as_vector(gen.events(12, 12.7)));
}
TEST(event_generators, seq) {
@@ -76,57 +60,44 @@ TEST(event_generators, seq) {
return pse_vector(in.begin()+b, in.begin()+e);
};
- event_generator gen = seq_generator<pse_vector>(in);
-
- // Test pop, next and reset.
- for (auto e: in) {
- EXPECT_EQ(e, gen.front());
- gen.pop();
- }
+ event_generator gen = explicit_generator(in);
+ EXPECT_EQ(in, as_vector(gen.events(0, 100.)));
gen.reset();
- for (auto e: in) {
- EXPECT_EQ(e, gen.front());
- gen.pop();
- }
- // The loop above should have drained all events from gen, so we expect
- // that the front() event will be the special terminal_pse event.
- EXPECT_TRUE(is_terminal_pse(gen.front()));
-
+ EXPECT_EQ(in, as_vector(gen.events(0, 100.)));
gen.reset();
- // Update of the input sequence, and run tests again to
- // verify that results reflect the new set of input events.
+ // Check reported sub-intervals against a smaller set of 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},
};
+ gen = explicit_generator(in);
+
+ auto draw = [](event_generator& gen, time_type t0, time_type t1) {
+ gen.reset();
+ return as_vector(gen.events(t0, t1));
+ };
// 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));
@@ -134,32 +105,22 @@ TEST(event_generators, seq) {
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.front().time<t1) {
- int1.push_back(gen.front());
- gen.pop();
- }
+ event_generator gen = poisson_generator(target, weight, t0, lambda, G);
+
+ pse_vector int1 = as_vector(gen.events(0, t1));
// 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.front().time<t1) {
- int2.push_back(gen.front());
- gen.pop();
- }
-
- // Assert that the same sequence was generated
+ pse_vector int2 = as_vector(gen.events(0, t1));
EXPECT_EQ(int1, int2);
}
diff --git a/test/unit/test_lif_cell_group.cpp b/test/unit/test_lif_cell_group.cpp
index 80b1909b..25297e25 100644
--- a/test/unit/test_lif_cell_group.cpp
+++ b/test/unit/test_lif_cell_group.cpp
@@ -6,6 +6,7 @@
#include <arbor/load_balance.hpp>
#include <arbor/threadinfo.hpp>
#include <arbor/recipe.hpp>
+#include <arbor/schedule.hpp>
#include <arbor/simulation.hpp>
#include <arbor/spike_source_cell.hpp>
@@ -61,7 +62,7 @@ public:
// regularly spiking cell.
if (gid == 0) {
// Produces just a single spike at time 0ms.
- return spike_source_cell{vector_time_seq({0.f})};
+ return spike_source_cell{explicit_schedule({0.f})};
}
// LIF cell.
return lif_cell();
diff --git a/test/unit/test_merge_events.cpp b/test/unit/test_merge_events.cpp
index b906c02e..d1756f4e 100644
--- a/test/unit/test_merge_events.cpp
+++ b/test/unit/test_merge_events.cpp
@@ -1,12 +1,39 @@
#include "../gtest.h"
-#include <event_queue.hpp>
-#include <merge_events.hpp>
+#include <vector>
+#include <arbor/event_generator.hpp>
+#include <arbor/spike_event.hpp>
+
+#include "merge_events.hpp"
#include "util/rangeutil.hpp"
+namespace arb {
+void merge_cell_events(
+ time_type t_from,
+ time_type t_to,
+ event_span old_events,
+ event_span pending,
+ std::vector<event_generator>& generators,
+ pse_vector& new_events);
+} // namespace arb
+
using namespace arb;
+// Wrapper for arb::merge_cell_events.
+static void merge_events(
+ time_type t_from,
+ time_type t_to,
+ const pse_vector& old_events,
+ pse_vector& pending,
+ std::vector<event_generator>& generators,
+ pse_vector& new_events)
+{
+ util::sort(pending);
+ merge_cell_events(t_from, t_to, util::range_pointer_view(old_events), util::range_pointer_view(pending), generators, new_events);
+}
+
+
std::vector<event_generator> empty_gens;
// Test the trivial case of merging empty sets
@@ -21,6 +48,7 @@ TEST(merge_events, empty)
EXPECT_EQ(lf.size(), 0u);
}
+
// Test the case where there are no events in lc that are to be delivered
// after tfinal.
TEST(merge_events, no_overlap)
@@ -128,9 +156,9 @@ TEST(merge_events, X)
{{3, 0}, 26, 4},
};
- std::vector<event_generator> generators(2);
- generators.emplace_back(
- regular_generator(cell_member_type{4,2}, 42.f, t0, 5));
+ std::vector<event_generator> generators = {
+ regular_generator(cell_member_type{4,2}, 42.f, t0, 5)
+ };
merge_events(t0, t1, lc, events, generators, lf);
@@ -154,7 +182,7 @@ TEST(merge_events, X)
// Test the tournament tree for merging two small sequences
TEST(merge_events, tourney_seq)
{
- pse_vector g1 = {
+ pse_vector evs1 = {
{{0, 0}, 1, 1},
{{0, 0}, 2, 2},
{{0, 0}, 3, 3},
@@ -162,7 +190,7 @@ TEST(merge_events, tourney_seq)
{{0, 0}, 5, 5},
};
- pse_vector g2 = {
+ pse_vector evs2 = {
{{0, 0}, 1.5, 1},
{{0, 0}, 2.5, 2},
{{0, 0}, 3.5, 3},
@@ -170,10 +198,13 @@ TEST(merge_events, tourney_seq)
{{0, 0}, 5.5, 5},
};
- std::vector<event_generator> generators;
- generators.emplace_back(seq_generator<pse_vector>(g1));
- generators.emplace_back(seq_generator<pse_vector>(g2));
- impl::tourney_tree tree(generators);
+ event_generator g1 = explicit_generator(evs1);
+ event_generator g2 = explicit_generator(evs2);
+
+ std::vector<event_span> spans;
+ spans.emplace_back(g1.events(0, terminal_time));
+ spans.emplace_back(g2.events(0, terminal_time));
+ impl::tourney_tree tree(spans);
pse_vector lf;
while (!tree.empty()) {
@@ -182,8 +213,8 @@ TEST(merge_events, tourney_seq)
}
EXPECT_TRUE(std::is_sorted(lf.begin(), lf.end()));
- auto expected = g1;
- util::append(expected, g2);
+ auto expected = evs1;
+ util::append(expected, evs2);
util::sort(expected);
EXPECT_EQ(expected, lf);
@@ -209,17 +240,16 @@ TEST(merge_events, tourney_poisson)
// of events with the same time but different weights works properly.
rndgen G(i%(ngen-1));
generators.emplace_back(
- poisson_generator<std::mt19937_64>(tgt, weight, G, t0, lambda));
+ poisson_generator(tgt, weight, t0, lambda, G));
}
// 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.front().time<tfinal) {
- expected.push_back(gen.front());
- gen.pop();
- }
+ event_span evs = gen.events(t0, tfinal);
+ util::append(expected, evs);
+
// Reset the generator so that it is ready to generate the same
// events again for the tournament tree test.
gen.reset();
@@ -228,9 +258,13 @@ TEST(merge_events, tourney_poisson)
util::sort(expected);
// Generate output using tournament tree in lf.
- impl::tourney_tree tree(generators);
+ std::vector<event_span> spans;
+ for (auto& gen: generators) {
+ spans.emplace_back(gen.events(t0, tfinal));
+ }
+ impl::tourney_tree tree(spans);
pse_vector lf;
- while (!tree.empty(tfinal)) {
+ while (!tree.empty()) {
lf.push_back(tree.head());
tree.pop();
}
diff --git a/test/unit/test_schedule.cpp b/test/unit/test_schedule.cpp
index 0796352b..c49a1a0e 100644
--- a/test/unit/test_schedule.cpp
+++ b/test/unit/test_schedule.cpp
@@ -6,8 +6,8 @@
#include <arbor/common_types.hpp>
#include <arbor/schedule.hpp>
-#include <util/partition.hpp>
-#include <util/rangeutil.hpp>
+#include "util/partition.hpp"
+#include "util/rangeutil.hpp"
#include "common.hpp"
#include "stats.hpp"
@@ -15,6 +15,8 @@
using namespace arb;
using namespace testing;
+using time_range = util::range<const time_type*>;
+
// Pull events from n non-contiguous subintervals of [t0, t1)
// and check for monotonicity and boundedness.
@@ -31,7 +33,7 @@ void run_invariant_checks(schedule S, time_type t0, time_type t1, unsigned n, in
bool skip = false;
for (auto ival: util::partition_view(divisions)) {
if (!skip) {
- auto ts = S.events(ival.first, ival.second);
+ time_range ts = S.events(ival.first, ival.second);
EXPECT_TRUE(std::is_sorted(ts.begin(), ts.end()));
if (!ts.empty()) {
@@ -63,31 +65,37 @@ void run_reset_check(schedule S, time_type t0, time_type t1, unsigned n, int see
std::vector<time_type> first;
for (auto ival: util::partition_view(first_div)) {
- util::append(first, S.events(ival.first, ival.second));
+ time_range ts = S.events(ival.first, ival.second);
+ util::append(first, ts);
}
S.reset();
std::vector<time_type> second;
for (auto ival: util::partition_view(second_div)) {
- util::append(second, S.events(ival.first, ival.second));
+ time_range ts = S.events(ival.first, ival.second);
+ util::append(second, ts);
}
EXPECT_EQ(first, second);
}
+static std::vector<time_type> as_vector(std::pair<const time_type*, const time_type*> ts) {
+ return std::vector<time_type>(ts.first, ts.second);
+}
+
TEST(schedule, regular) {
// Use exact fp representations for strict equality testing.
std::vector<time_type> expected = {0, 0.25, 0.5, 0.75, 1.0};
schedule S = regular_schedule(0.25);
- EXPECT_EQ(expected, S.events(0, 1.25));
+ EXPECT_EQ(expected, as_vector(S.events(0, 1.25)));
S.reset();
- EXPECT_EQ(expected, S.events(0, 1.25));
+ EXPECT_EQ(expected, as_vector(S.events(0, 1.25)));
S.reset();
expected = {0.25, 0.5, 0.75, 1.0};
- EXPECT_EQ(expected, S.events(0.1, 1.01));
+ EXPECT_EQ(expected, as_vector(S.events(0.1, 1.01)));
}
TEST(schedule, regular_invariants) {
@@ -100,19 +108,54 @@ TEST(schedule, regular_reset) {
run_reset_check(regular_schedule(0.3), 3, 12, 7);
}
+TEST(schedule, regular_rounding) {
+ // Test for consistent behaviour in the face of rounding at large time values.
+ // Example: with t1, dt below, and int n = floor(t0/dt),
+ // then n*dt is not the smallest multiple of dt greater than or equal to t0.
+ // In fact, (n-4)*dt is still greater than t0.
+
+ time_type t1 = 1802667.f;
+ time_type dt = 0.024999f;
+
+ time_type t0 = t1-10*dt;
+ time_type t2 = t1+10*dt;
+
+ schedule S = regular_schedule(t0, dt);
+ auto int_l = as_vector(S.events(t0, t1));
+ auto int_r = as_vector(S.events(t1, t2));
+
+ S.reset();
+ auto int_a = as_vector(S.events(t0, t2));
+
+ EXPECT_GE(int_l.front(), t0);
+ EXPECT_LT(int_l.back(), t1);
+
+ EXPECT_GE(int_r.front(), t1);
+ EXPECT_LT(int_r.back(), t2);
+
+ EXPECT_GE(int_a.front(), t0);
+ EXPECT_LT(int_a.back(), t2);
+
+ std::vector<time_type> int_merged = int_l;
+ util::append(int_merged, int_r);
+
+ EXPECT_EQ(int_merged, int_a);
+ EXPECT_TRUE(util::is_sorted(int_a));
+}
+
TEST(schedule, explicit_schedule) {
time_type times[] = {0.1, 0.3, 1.0, 1.25, 1.7, 2.2};
std::vector<time_type> expected = {0.1, 0.3, 1.0};
schedule S = explicit_schedule(times);
- EXPECT_EQ(expected, S.events(0, 1.25));
+ EXPECT_EQ(expected, as_vector(S.events(0, 1.25)));
S.reset();
- EXPECT_EQ(expected, S.events(0, 1.25));
+ EXPECT_EQ(expected, as_vector(S.events(0, 1.25)));
S.reset();
expected = {0.3, 1.0, 1.25, 1.7};
- EXPECT_EQ(expected, S.events(0.3, 1.71));
+ EXPECT_EQ(expected, as_vector(S.events(0.3, 1.71)));
}
TEST(schedule, explicit_invariants) {
@@ -157,11 +200,11 @@ private:
};
template <typename RNG>
-double poisson_schedule_dispersion(int nbin, double mean_dt, RNG& G) {
- schedule S = poisson_schedule(mean_dt, G);
+double poisson_schedule_dispersion(int nbin, double rate_kHz, RNG& G) {
+ schedule S = poisson_schedule(rate_kHz, G);
std::vector<int> bin(nbin);
- for (auto t: S.events(0, nbin)) {
+ for (auto t: time_range(S.events(0, nbin))) {
int j = (int)t;
if (j<0 || j>=nbin) {
throw std::logic_error("poisson schedule result out of bounds");
@@ -196,7 +239,7 @@ TEST(schedule, poisson_uniformity) {
constexpr double chi2_ub = 1118.9480663231843;
std::mt19937_64 G;
- double dispersion = poisson_schedule_dispersion(N, 1.23, G);
+ double dispersion = poisson_schedule_dispersion(N, .813, G);
double test_value = N*dispersion;
EXPECT_GT(test_value, chi2_lb);
EXPECT_LT(test_value, chi2_ub);
@@ -204,8 +247,8 @@ TEST(schedule, poisson_uniformity) {
// Run one sample K-S test for uniformity, with critical
// value for the finite K-S statistic Dn of α=0.01.
- schedule S = poisson_schedule(1./100, G);
- auto events = S.events(0,1);
+ schedule S = poisson_schedule(100., G);
+ auto events = as_vector(S.events(0,1));
int n = (int)events.size();
double dn = ks::dn_statistic(events);
@@ -215,13 +258,13 @@ TEST(schedule, poisson_uniformity) {
// source.
skew_adaptor<std::mt19937_64> W(1.5);
- dispersion = poisson_schedule_dispersion(N, 1.23, W);
+ dispersion = poisson_schedule_dispersion(N, .813, W);
test_value = N*dispersion;
EXPECT_FALSE(test_value>=chi2_lb && test_value<=chi2_ub);
- S = poisson_schedule(1./100, W);
- events = S.events(0,1);
+ S = poisson_schedule(100., W);
+ events = as_vector(S.events(0,1));
n = (int)events.size();
dn = ks::dn_statistic(events);
@@ -242,8 +285,8 @@ TEST(schedule, poisson_rate) {
constexpr double lambda = 123.4;
std::mt19937_64 G;
- schedule S = poisson_schedule(1./lambda, G);
- int n = (int)S.events(0,1).size();
+ schedule S = poisson_schedule(lambda, G);
+ int n = (int)time_range(S.events(0, 1)).size();
double cdf = poisson::poisson_cdf_approx(n, lambda);
EXPECT_GT(cdf, alpha/2);
@@ -253,8 +296,8 @@ TEST(schedule, poisson_rate) {
// source.
skew_adaptor<std::mt19937_64> W(1.5);
- S = poisson_schedule(1./lambda, W);
- n = (int)S.events(0,1).size();
+ S = poisson_schedule(lambda, W);
+ n = (int)time_range(S.events(0, 1)).size();
cdf = poisson::poisson_cdf_approx(n, lambda);
EXPECT_FALSE(cdf>=alpha/2 && cdf<=1-alpha/2);
@@ -264,14 +307,14 @@ TEST(schedule, poisson_invariants) {
SCOPED_TRACE("poisson_invariants");
std::mt19937_64 G;
G.discard(100);
- run_invariant_checks(poisson_schedule(12.3, G), 5.1, 15.3, 7);
+ run_invariant_checks(poisson_schedule(0.81, G), 5.1, 15.3, 7);
}
TEST(schedule, poisson_reset) {
SCOPED_TRACE("poisson_reset");
std::mt19937_64 G;
G.discard(200);
- run_reset_check(poisson_schedule(9.1, G), 1, 10, 7);
+ run_reset_check(poisson_schedule(.11, G), 1, 10, 7);
}
TEST(schedule, poisson_offset) {
@@ -284,7 +327,7 @@ TEST(schedule, poisson_offset) {
G1.discard(300);
std::vector<time_type> expected;
- for (auto t: poisson_schedule(3.4, G1).events(0., 100.)) {
+ for (auto t: as_vector(poisson_schedule(.234, G1).events(0., 100.))) {
t += offset;
if (t<100.) {
expected.push_back(t);
@@ -294,13 +337,14 @@ TEST(schedule, poisson_offset) {
std::mt19937_64 G2;
G2.discard(300);
- EXPECT_TRUE(seq_almost_eq<time_type>(expected, poisson_schedule(offset, 3.4, G2).events(0., 100.)));
+ EXPECT_TRUE(seq_almost_eq<time_type>(expected,
+ as_vector(poisson_schedule(offset, .234, G2).events(0., 100.))));
}
TEST(schedule, poisson_offset_reset) {
SCOPED_TRACE("poisson_reset");
std::mt19937_64 G;
G.discard(400);
- run_reset_check(poisson_schedule(0.3, 9.1, G), 1, 10, 7);
+ run_reset_check(poisson_schedule(3.3, 9.1, G), 1, 10, 7);
}
diff --git a/test/unit/test_spike_source.cpp b/test/unit/test_spike_source.cpp
index b877a64a..8892b774 100644
--- a/test/unit/test_spike_source.cpp
+++ b/test/unit/test_spike_source.cpp
@@ -1,7 +1,8 @@
#include "../gtest.h"
+#include <arbor/schedule.hpp>
+#include <arbor/spike.hpp>
#include <arbor/spike_source_cell.hpp>
-#include <arbor/time_sequence.hpp>
#include <arbor/util/unique_any.hpp>
#include "spike_source_cell_group.hpp"
@@ -10,57 +11,61 @@
using namespace arb;
using ss_recipe = homogeneous_recipe<cell_kind::spike_source, spike_source_cell>;
-using pseq = arb::poisson_time_seq<std::mt19937_64>;
// Test that a spike_source_cell_group identifies itself with the correct
// cell_kind enum value.
TEST(spike_source, cell_kind)
{
- ss_recipe rec(1u, spike_source_cell{vector_time_seq({})});
+ ss_recipe rec(1u, spike_source_cell{explicit_schedule({})});
spike_source_cell_group group({0}, rec);
EXPECT_EQ(cell_kind::spike_source, group.get_cell_kind());
}
-// Test that a spike_source_cell_group produces a sequence spikes with spike
+static std::vector<time_type> as_vector(std::pair<const time_type*, const time_type*> ts) {
+ return std::vector<time_type>(ts.first, ts.second);
+}
+
+static std::vector<time_type> spike_times(const std::vector<spike>& evs) {
+ std::vector<time_type> ts;
+ for (auto& s: evs) {
+ ts.push_back(s.time);
+ }
+ return ts;
+}
+
+// Test that a spike_source_cell_group produces a sequence of spikes with spike
// times corresponding to the underlying time_seq.
TEST(spike_source, matches_time_seq)
{
- auto test_seq = [](arb::time_seq seq) {
+ auto test_seq = [](schedule seq) {
ss_recipe rec(1u, spike_source_cell{seq});
spike_source_cell_group group({0}, rec);
// epoch ending at 10ms
epoch ep(0, 10);
group.advance(ep, 1, {});
- for (auto s: group.spikes()) {
- EXPECT_EQ(s.time, seq.front());
- seq.pop();
- }
- EXPECT_TRUE(seq.front()>=ep.tfinal);
+ EXPECT_EQ(spike_times(group.spikes()), as_vector(seq.events(0, 10)));
+
group.clear_spikes();
// advance to 20 ms and repeat
ep.advance(20);
group.advance(ep, 1, {});
- for (auto s: group.spikes()) {
- EXPECT_EQ(s.time, seq.front());
- seq.pop();
- }
- EXPECT_TRUE(seq.front()>=ep.tfinal);
+ EXPECT_EQ(spike_times(group.spikes()), as_vector(seq.events(10, 20)));
};
std::mt19937_64 G;
- test_seq(arb::regular_time_seq(0,1));
- test_seq(pseq(G, 0., 10)); // produce many spikes in each interval
- test_seq(pseq(G, 0., 1e-6)); // very unlikely to produce any spikes in either interval
+ test_seq(regular_schedule(0, 1));
+ test_seq(poisson_schedule(10, G)); // produce many spikes in each interval
+ test_seq(poisson_schedule(1e-6, G)); // very unlikely to produce any spikes in either interval
}
// Test that a spike_source_cell_group will produce the same sequence of spikes
// after being reset.
TEST(spike_source, reset)
{
- auto test_seq = [](arb::time_seq seq) {
+ auto test_seq = [](schedule seq) {
ss_recipe rec(1u, spike_source_cell{seq});
spike_source_cell_group group({0}, rec);
@@ -80,9 +85,9 @@ TEST(spike_source, reset)
};
std::mt19937_64 G;
- test_seq(arb::regular_time_seq(0,1));
- test_seq(pseq(G, 0., 10)); // produce many spikes in each interval
- test_seq(pseq(G, 0., 1e-6)); // very unlikely to produce any spikes in either interval
+ test_seq(regular_schedule(0, 1));
+ test_seq(poisson_schedule(10, G)); // produce many spikes in each interval
+ test_seq(poisson_schedule(1e-6, G)); // very unlikely to produce any spikes in either interval
}
// Test that a spike_source_cell_group will produce the expected
@@ -90,26 +95,19 @@ TEST(spike_source, reset)
TEST(spike_source, exhaust)
{
// This test assumes that seq will exhaust itself before t=10 ms.
- auto test_seq = [](arb::time_seq seq) {
+ auto test_seq = [](schedule seq) {
ss_recipe rec(1u, spike_source_cell{seq});
spike_source_cell_group group({0}, rec);
// epoch ending at 10ms
epoch ep(0, 10);
group.advance(ep, 1, {});
- auto spikes = group.spikes();
- for (auto s: group.spikes()) {
- EXPECT_EQ(s.time, seq.front());
- seq.pop();
- }
- // The sequence shoule be exhausted, in which case the next value in the
- // sequence should be marked as time_max.
- EXPECT_EQ(seq.front(), arb::terminal_time);
+ EXPECT_EQ(spike_times(group.spikes()), as_vector(seq.events(0, 10)));
+
// Check that the last spike was before the end of the epoch.
- EXPECT_LT(spikes.back().time, time_type(10));
+ EXPECT_LT(group.spikes().back().time, time_type(10));
};
- std::mt19937_64 G;
- test_seq(arb::regular_time_seq(0,1,5));
- test_seq(pseq(G, 0., 10, 5));
+ test_seq(regular_schedule(0, 1, 5));
+ test_seq(explicit_schedule({0.3, 2.3, 4.7}));
}
diff --git a/test/unit/test_time_seq.cpp b/test/unit/test_time_seq.cpp
deleted file mode 100644
index b4166279..00000000
--- a/test/unit/test_time_seq.cpp
+++ /dev/null
@@ -1,175 +0,0 @@
-#include "../gtest.h"
-
-#include <vector>
-
-#include <arbor/time_sequence.hpp>
-
-#include "util/rangeutil.hpp"
-
-#include "common.hpp"
-
-using namespace arb;
-
-namespace{
- // Helper function that draws all samples in the half open interval
- // t ∈ [t0, t1) from a time_seq.
- std::vector<time_type> draw(time_seq& gen, time_type t0, time_type t1) {
- gen.reset();
- gen.advance(t0);
- std::vector<time_type> v;
- while (gen.front()<t1) {
- v.push_back(gen.front());
- gen.pop();
- }
- return v;
- }
-}
-
-TEST(time_seq, vector) {
- std::vector<time_type> times = {0.1, 1.0, 1.0, 1.5, 2.3, 3.0, 3.5, };
-
- vector_time_seq seq(times);
-
- // Test pop, next and reset.
- for (auto t: times) {
- EXPECT_EQ(t, seq.front());
- seq.pop();
- }
- seq.reset();
- for (auto t: times) {
- EXPECT_EQ(t, seq.front());
- seq.pop();
- }
-
- // The loop above should have drained all samples from seq, so we expect
- // that the front() time sample will be terminal_time.
- EXPECT_EQ(seq.front(), terminal_time);
-}
-
-TEST(time_seq, regular) {
- // make a regular generator that generates its first event at t=2ms and subsequent
- // events regularly spaced 0.5 ms apart.
- regular_time_seq seq(2, 0.5);
-
- // Test pop, next and reset.
- for (auto e: {2.0, 2.5, 3.0, 3.5, 4.0, 4.5}) {
- EXPECT_EQ(e, seq.front());
- seq.pop();
- }
- seq.reset();
- for (auto e: {2.0, 2.5, 3.0, 3.5, 4.0, 4.5}) {
- EXPECT_EQ(e, seq.front());
- seq.pop();
- }
- seq.reset();
-
- // Test advance()
-
- seq.advance(10.1);
- // next event greater ≥ 10.1 should be 10.5
- EXPECT_EQ(seq.front(), time_type(10.5));
-
- seq.advance(12);
- // next event greater ≥ 12 should be 12
- EXPECT_EQ(seq.front(), time_type(12));
-}
-
-// Test for rounding problems with large time values and the regular sequence
-TEST(time_seq, regular_rounding) {
- // make a regular generator that generates its first time point at t=2ms
- // and subsequent times regularly spaced 0.5 ms apart.
- time_type t0 = 2.0;
- time_type dt = 0.5;
-
- // 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;
- time_seq seq = regular_time_seq(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 points from each interval then merge them, we expect same set
- // of points as when we draw from that large interval.
- std::vector<time_type> int_l = draw(seq, t0, t1);
- std::vector<time_type> int_r = draw(seq, t1, t2);
- std::vector<time_type> int_a = draw(seq, t0, t2);
-
- std::vector<time_type> int_merged = int_l;
- util::append(int_merged, int_r);
-
- EXPECT_TRUE(int_l.front() >= t0);
- EXPECT_TRUE(int_l.back() < t1);
- EXPECT_TRUE(int_r.front() >= t1);
- EXPECT_TRUE(int_r.back() < t2);
- EXPECT_EQ(int_a, int_merged);
- EXPECT_TRUE(std::is_sorted(int_a.begin(), int_a.end()));
-}
-
-TEST(time_seq, poisson) {
- std::mt19937_64 G;
- using pseq = poisson_time_seq<std::mt19937_64>;
-
- time_type t0 = 0;
- time_type t1 = 10;
- time_type lambda = 10; // expect 10 samples per ms
-
- pseq seq(G, t0, lambda);
-
- std::vector<time_type> int1;
- while (seq.front()<t1) {
- int1.push_back(seq.front());
- seq.pop();
- }
- // Test that the output is sorted
- EXPECT_TRUE(std::is_sorted(int1.begin(), int1.end()));
-
- // Reset and generate the same sequence of time points
- seq.reset();
- std::vector<time_type> int2;
- while (seq.front()<t1) {
- int2.push_back(seq.front());
- seq.pop();
- }
-
- // Assert that the same sequence was generated
- EXPECT_EQ(int1, int2);
-}
-
-// Test a poisson generator that has a tstop past which no samples
-// should be generated.
-TEST(time_seq, poisson_terminates) {
- std::mt19937_64 G;
- using pseq = poisson_time_seq<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 samples per ms
-
- // construct sequence with explicit end time t1
- pseq seq(G, t0, lambda, t1);
-
- std::vector<time_type> sequence;
- // pull samples off the sequence well past the end of the end time t1
- while (seq.front()<t2) {
- sequence.push_back(seq.front());
- seq.pop();
- }
-
- // the sequence should be exhausted
- EXPECT_EQ(seq.front(), terminal_time);
-
- // the last sample should be less than the end time
- EXPECT_TRUE(sequence.back()<t1);
-}
--
GitLab