diff --git a/doc/.gitignore b/doc/.gitignore
deleted file mode 100644
index 208d1dced07ec6fe895b65fbb03c282390b031cf..0000000000000000000000000000000000000000
--- a/doc/.gitignore
+++ /dev/null
@@ -1,4 +0,0 @@
-# exclude working directories
-
-build
-static
diff --git a/doc/conf.py b/doc/conf.py
index 1de522769fc6818c5779d259ceee00873a609628..e81cb603bed70a431153222676f2c871f60822c0 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -21,6 +21,8 @@
# import sys
# sys.path.insert(0, os.path.abspath('.'))
+def setup(app):
+ app.add_stylesheet('custom.css')
# -- General configuration ------------------------------------------------
diff --git a/doc/index.rst b/doc/index.rst
index 269ff4b65fc57e8b556dc57a3605be64d0a03da0..32b94e403c57de44338c32dc599400f2ce335f72 100644
--- a/doc/index.rst
+++ b/doc/index.rst
@@ -25,4 +25,5 @@ Arbor is a high-performance library for computationa neurscience simulations.
:maxdepth: 1
library
+ sampling_api
diff --git a/doc/sampling_api.rst b/doc/sampling_api.rst
new file mode 100644
index 0000000000000000000000000000000000000000..79fa620dc27b7e428a48c64f0c1abeb874c3f5e4
--- /dev/null
+++ b/doc/sampling_api.rst
@@ -0,0 +1,275 @@
+Sampling API
+============
+
+The new API replaces the flexible but irreducibly inefficient scheme
+where the next sample time for a sampling was determined by the
+return value of the sampler callback.
+
+
+Definitions
+-----------
+
+probe
+ A location or component of a cell that is available for monitoring.
+
+sample
+ A record of data corresponding to the value at a specific *probe* at a specific time.
+
+sampler
+ A function or function object that receives a sequence of *sample* records.
+
+schedule
+ A function or function object that, given a time interval, returns a list of sample times within that interval.
+
+
+
+Probes
+------
+
+Probes are specified in the recipe objects that are used to initialize a
+model; the specification of the item or value that is subjected to a
+probe will be specific to a particular cell type.
+
+.. container:: api-code
+
+ .. code-block:: cpp
+
+ using probe_tag = int;
+
+ struct probe_info {
+ cell_member_type id; // cell gid, index of probe
+ probe_tag tag; // opaque key, returned in sample record
+ any address; // cell-type specific location info
+ };
+
+ probe_info recipe::get_probe(cell_member_type probe_id);
+
+
+The ``id`` field in the ``probe_info`` struct will be the same value as
+the ``probe_id`` used in the ``get_probe`` call.
+
+The ``get_probe()`` method is intended for use by cell group
+implementations to set up sampling data structures ahead of time and for
+putting in place any structures or information in the concrete cell
+implementations to allow monitoring.
+
+The ``tag`` field has no semantics for the engine. It is provided merely
+as a way of passing additional metadata about a probe to any sampler
+that polls it, with a view to samplers that handle multiple probes,
+possibly with different value types.
+
+Probe addresses are now decoupled from the cell descriptions themselves —
+this allows a recipe implementation to construct probes independently
+of the cells themselves. It is the responsibility of a cell group implementation
+to parse the probe address objects wrapped in the ``any address`` field.
+
+
+Samplers and sample records
+---------------------------
+
+Data collected from probes (according to a schedule described below)
+will be passed to a sampler function or function object:
+
+.. container:: api-code
+
+ .. code-block:: cpp
+
+ using sampler_function =
+ std::function<void (cell_member_type, probe_tag, size_t, const sample_record*)>;
+
+where the parameters are respectively the probe id, the tag, the number
+of samples and a pointer to the sequence of sample records.
+
+The ``probe_tag`` is the key given in the ``probe_info`` returned by
+the recipe.
+
+One ``sample_record`` struct contains one sample of the probe data at a
+given simulation time point:
+
+.. container:: api-code
+
+ .. code-block:: cpp
+
+ struct sample_record {
+ time_type time; // simulation time of sample
+ any_ptr data; // sample data
+ };
+
+The ``data`` field points to the sample data, wrapped in ``any_ptr`` for
+type-checked access. The exact representation will depend on the nature of
+the object that is being probed, but it should depend only on the cell type and
+probe address.
+
+The data pointed to by ``data``, and the sample records themselves, are
+only guaranteed to be valid for the duration of the call to the sampler
+function. A simple sampler implementation for ``double`` data might be:
+
+.. container:: example-code
+
+ .. code-block:: cpp
+
+ using sample_data = std::map<cell_member_type, std::vector<std::pair<double, double>>>;
+
+ struct scalar_sampler {
+ sample_data& samples;
+
+ explicit scalar_sample(sample_data& samples): samples(samples) {}
+
+ void operator()(cell_member_type id, probe_tag, size_t n, const sample_record* records) {
+ for (size_t i=0; i<n; ++i) {
+ const auto& rec = records[i];
+
+ const double* data = any_cast<const double*>(rec.data);
+ assert(data);
+ samples[id].emplace_back(rec.time, *data);
+ }
+ }
+ };
+
+The use of ``any_ptr`` allows type-checked access to the sample data, which
+may differ in type from probe to probe.
+
+
+Model and cell group interface
+------------------------------
+
+Polling rates, policies and sampler functions are set through the
+``model`` interface, after construction from a recipe.
+
+.. container:: api-code
+
+ .. code-block:: cpp
+
+ using sampler_association_handle = std::size_t;
+ using cell_member_predicate = std::function<bool (cell_member_type)>;
+
+ sampler_association_handle model::add_sampler(
+ cell_member_predicate probe_ids,
+ schedule sched,
+ sampler_function fn,
+ sampling_policy policy = sampling_policy::lax);
+
+ void model::remove_sampler(sampler_association_handle);
+
+ void model::remove_all_samplers();
+
+Multiple samplers can then be associated with the same probe locations.
+The handle returned is only used for managing the lifetime of the
+association. The ``cell_member_predicate`` parameter defines the
+set of probe ids in terms of a membership test.
+
+Two helper functions are provided for making ``cell_member_predicate`` objects:
+
+.. container:: api-code
+
+ .. code-block:: cpp
+
+ // Match all probe ids.
+ cell_member_predicate all_probes = [](cell_member_type pid) { return true; };
+
+ // Match just one probe id.
+ cell_member_predicate one_probe(cell_member_type pid) {
+ return [pid](cell_member_type x) { return pid==x; };
+ }
+
+
+The ``sampling_policy`` policy is used to modify sampling behaviour: by
+default, the ``lax`` policy is to perform a best-effort sampling that
+minimizes sampling overhead and which will not change the numerical
+behaviour of the simulation. Other policies may be implemented in the
+future, e.g. ``interpolated`` or ``exact``.
+
+The model object will pass on the sampler setting request to the cell
+group that owns the given probe id. The ``cell_group`` interface will be
+correspondingly extended:
+
+.. container:: api-code
+
+ .. code-block:: cpp
+
+ void cell_group::add_sampler(sampler_association_handle h, cell_member_predicate probe_ids, sample_schedule sched, sampler_function fn, sampling_policy policy);
+
+ void cell_group::remove_sampler(sampler_association_handle);
+
+ void cell_group::remove_all_samplers();
+
+Cell groups will invoke the corresponding sampler function directly, and
+may aggregate multiple samples with the same probe id in one call to the
+sampler. Calls to the sampler are synchronous, in the sense that
+processing of the cell group state does not proceed while the sampler
+function is being executed, but the times of the samples given to the
+sampler will typically precede the time corresponding to the current
+state of the cell group. It should be expected that this difference in
+time should be no greater the the duration of the integration period
+(i.e. ``mindelay/2``).
+
+If a cell group does not support a given ``sampling_policy``, it should
+raise an exception. All cell groups should support the ``lax`` policy,
+if they support probes at all.
+
+
+Schedules
+---------
+
+Schedules represent a non-negative, monotonically increasing sequence
+of time points, and are used to specify the sampling schedule in any
+given association of a sampler function to a set of probes.
+
+A ``schedule`` object has two methods:
+
+.. container:: api-code
+
+ .. code-block:: cpp
+
+ void schedule::reset();
+
+ std::vector<time_type> events(time_type t0, time_type t1)
+
+The ``events(t0, t1)`` method returns a vector 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 ``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
+and deterministic behaviour.
+
+The ``schedule`` object itself uses type-erasure to wrap any schedule
+implementation class, which can be any copy--constructible class that
+provides the methods ``reset()`` and ``events(t0, t1)`` above. Three
+schedule implementations are provided by the engine:
+
+.. container:: api-code
+
+ .. code-block:: cpp
+
+
+ // Schedule at integer multiples of dt:
+ schedule regular_schedule(time_type dt);
+
+ // Schedule at a predetermined (sorted) sequence of times:
+ template <typename Seq>
+ schedule explicit_schedule(const Seq& seq);
+
+ // Schedule according to Poisson process with lambda = 1/mean_dt
+ template <typename RandomNumberEngine>
+ schedule poisson_schedule(time_type mean_dt, const RandomNumberEngine& rng);
+
+The ``schedule`` class and its implementations are found in ``schedule.hpp``.
+
+
+Helper classes for probe/sampler management
+-------------------------------------------
+
+The ``model`` and ``mc_cell_group`` classes use classes defined in ``scheduler_map.hpp`` to simplify
+the management of sampler--probe associations and probe metdata.
+
+``sampler_association_map`` wraps an ``unordered_map`` between sampler assocation
+handles and tuples (*schedule*, *sampler*, *probe set*), with thread-safe
+accessors.
+
+``probe_assocation_map<Handle>`` is a type alias for an unordered map between
+probe ids and tuples (*probe handle*, *probe tag*), where the *probe handle*
+is a cell-group specific accessor that allows efficient polling.
+
diff --git a/doc/static/custom.css b/doc/static/custom.css
new file mode 100644
index 0000000000000000000000000000000000000000..005eac316b8d6274be5c4749129a7b8780a692fe
--- /dev/null
+++ b/doc/static/custom.css
@@ -0,0 +1,22 @@
+.example-code>div[class^='highlight']:before {
+ font-weight: bold;
+ display: block;
+ padding-bottom: 2pt;
+ padding-left: 3pt;
+ content: "Example";
+ background: #6ab0de;
+ color: #ffffff
+}
+
+.api-code>div[class^='highlight']:before {
+ font-weight: bold;
+ display: block;
+ padding-bottom: 2pt;
+ padding-left: 3pt;
+ content: "API";
+ background: #e0e0e0
+}
+
+.api-code>div[class^='highlight'] {
+ background: #f4f4f4
+}
diff --git a/miniapp/CMakeLists.txt b/miniapp/CMakeLists.txt
index c41f7c5137147ccc7724f05bdcfed29a75fccaef..95c42b1f15536202fb3ef6b2ffcc049d296e2d69 100644
--- a/miniapp/CMakeLists.txt
+++ b/miniapp/CMakeLists.txt
@@ -5,12 +5,14 @@ set(MINIAPP_SOURCES
io.cpp
miniapp_recipes.cpp
morphology_pool.cpp
+ trace.cpp
)
set(MINIAPP_SOURCES_CUDA
miniapp.cu
io.cpp
miniapp_recipes.cpp
morphology_pool.cpp
+ trace.cpp
)
add_executable(miniapp.exe ${MINIAPP_SOURCES} ${HEADERS})
diff --git a/miniapp/miniapp.cpp b/miniapp/miniapp.cpp
index 05b1682b431208f8472911da0f255eb772f86e45..66b0ff2eaf3f8259aa43968bcada4fa5373ccea2 100644
--- a/miniapp/miniapp.cpp
+++ b/miniapp/miniapp.cpp
@@ -19,7 +19,10 @@
#include <model.hpp>
#include <profiling/profiler.hpp>
#include <profiling/meter_manager.hpp>
+#include <sampling.hpp>
+#include <schedule.hpp>
#include <threading/threading.hpp>
+#include <util/any.hpp>
#include <util/config.hpp>
#include <util/debug.hpp>
#include <util/ioutil.hpp>
@@ -27,27 +30,28 @@
#include "io.hpp"
#include "miniapp_recipes.hpp"
-#include "trace_sampler.hpp"
+#include "trace.hpp"
using namespace nest::mc;
+using util::any_cast;
+using util::make_span;
+
using global_policy = communication::global_policy;
-using sample_trace_type = sample_trace<time_type, double>;
using file_export_type = io::exporter_spike_file<global_policy>;
+using communicator_type = communication::communicator<communication::global_policy>;
+
void banner(hw::node_info);
std::unique_ptr<recipe> make_recipe(const io::cl_options&, const probe_distribution&);
-std::unique_ptr<sample_trace_type> make_trace(probe_record probe);
-using communicator_type = communication::communicator<communication::global_policy>;
+sample_trace make_trace(const probe_info& probe);
-void write_trace_json(const sample_trace_type& trace, const std::string& prefix = "trace_");
-void write_trace_csv(const sample_trace_type& trace, const std::string& prefix = "trace_");
void report_compartment_stats(const recipe&);
int main(int argc, char** argv) {
- nest::mc::communication::global_policy_guard global_guard(argc, argv);
+ communication::global_policy_guard global_guard(argc, argv);
try {
- nest::mc::util::meter_manager meters;
+ util::meter_manager meters;
meters.start();
std::cout << util::mask_stream(global_policy::id()==0);
@@ -85,6 +89,9 @@ int main(int argc, char** argv) {
pdist.all_segments = !options.probe_soma_only;
auto recipe = make_recipe(options, pdist);
+ if (options.report_compartments) {
+ report_compartment_stats(*recipe);
+ }
auto register_exporter = [] (const io::cl_options& options) {
return
@@ -94,11 +101,28 @@ int main(int argc, char** argv) {
};
auto decomp = partition_load_balance(*recipe, nd);
-
model m(*recipe, decomp);
- if (options.report_compartments) {
- report_compartment_stats(*recipe);
+ // Set up samplers for probes on local cable cells, as requested
+ // by command line options.
+ std::vector<sample_trace> sample_traces;
+ for (const auto& g: decomp.groups) {
+ if (g.kind==cable1d_neuron) {
+ for (auto gid: g.gids) {
+ if (options.trace_max_gid && gid>*options.trace_max_gid) {
+ continue;
+ }
+
+ for (cell_lid_type j: make_span(0, recipe->num_probes(gid))) {
+ sample_traces.push_back(make_trace(recipe->get_probe({gid, j})));
+ }
+ }
+ }
+ }
+
+ auto ssched = regular_schedule(options.sample_dt);
+ for (auto& trace: sample_traces) {
+ m.add_sampler(one_probe(trace.probe_id), ssched, make_simple_sampler(trace.samples));
}
// Specify event binning/coalescing.
@@ -109,18 +133,6 @@ int main(int argc, char** argv) {
m.set_binning_policy(binning_policy, options.bin_dt);
- // Attach samplers to all probes
- std::vector<std::unique_ptr<sample_trace_type>> traces;
- const time_type sample_dt = options.sample_dt;
- for (auto probe: m.probes()) {
- if (options.trace_max_gid && probe.id.gid>*options.trace_max_gid) {
- continue;
- }
-
- traces.push_back(make_trace(probe));
- m.attach_sampler(probe.id, make_trace_sampler(traces.back().get(), sample_dt));
- }
-
// Initialize the spike exporting interface
std::unique_ptr<file_export_type> file_exporter;
if (options.spike_file_output) {
@@ -153,8 +165,8 @@ int main(int argc, char** argv) {
// save traces
auto write_trace = options.trace_format=="json"? write_trace_json: write_trace_csv;
- for (const auto& trace: traces) {
- write_trace(*trace.get(), options.trace_prefix);
+ for (const auto& trace: sample_traces) {
+ write_trace(trace, options.trace_prefix);
}
auto report = util::make_meter_report(meters);
@@ -223,59 +235,25 @@ std::unique_ptr<recipe> make_recipe(const io::cl_options& options, const probe_d
}
}
-std::unique_ptr<sample_trace_type> make_trace(probe_record probe) {
+sample_trace make_trace(const probe_info& probe) {
std::string name = "";
std::string units = "";
- switch (probe.kind) {
- case probeKind::membrane_voltage:
+ auto addr = any_cast<cell_probe_address>(probe.address);
+ switch (addr.kind) {
+ case cell_probe_address::membrane_voltage:
name = "v";
units = "mV";
break;
- case probeKind::membrane_current:
+ case cell_probe_address::membrane_current:
name = "i";
units = "mA/cm²";
break;
default: ;
}
- name += probe.location.segment? "dend" : "soma";
-
- return util::make_unique<sample_trace_type>(probe.id, name, units);
-}
-
-void write_trace_csv(const sample_trace_type& trace, const std::string& prefix) {
- auto path = prefix + std::to_string(trace.probe_id.gid) +
- "." + std::to_string(trace.probe_id.index) + "_" + trace.name + ".csv";
-
- std::ofstream file(path);
- file << "# cell: " << trace.probe_id.gid << "\n";
- file << "# probe: " << trace.probe_id.index << "\n";
- file << "time_ms, " << trace.name << "_" << trace.units << "\n";
-
- for (const auto& sample: trace.samples) {
- file << util::strprintf("% 20.15f, % 20.15f\n", sample.time, sample.value);
- }
-}
-
-void write_trace_json(const sample_trace_type& trace, const std::string& prefix) {
- auto path = prefix + std::to_string(trace.probe_id.gid) +
- "." + std::to_string(trace.probe_id.index) + "_" + trace.name + ".json";
+ name += addr.location.segment? "dend" : "soma";
- nlohmann::json jrep;
- jrep["name"] = trace.name;
- jrep["units"] = trace.units;
- jrep["cell"] = trace.probe_id.gid;
- jrep["probe"] = trace.probe_id.index;
-
- auto& jt = jrep["data"]["time"];
- auto& jy = jrep["data"][trace.name];
-
- for (const auto& sample: trace.samples) {
- jt.push_back(sample.time);
- jy.push_back(sample.value);
- }
- std::ofstream file(path);
- file << std::setw(1) << jrep << std::endl;
+ return sample_trace{probe.id, name, units};
}
void report_compartment_stats(const recipe& rec) {
@@ -287,7 +265,7 @@ void report_compartment_stats(const recipe& rec) {
for (std::size_t i = 0; i<ncell; ++i) {
std::size_t ncomp = 0;
auto c = rec.get_cell_description(i);
- if (auto ptr = util::any_cast<cell>(&c)) {
+ if (auto ptr = any_cast<cell>(&c)) {
ncomp = ptr->num_compartments();
}
ncomp_total += ncomp;
diff --git a/miniapp/miniapp_recipes.cpp b/miniapp/miniapp_recipes.cpp
index 6fe6896d97a05d2ea9d4f3e39e150605e2471c53..a2039a7e19e65308e23a54a7167c1d1634cb4adc 100644
--- a/miniapp/miniapp_recipes.cpp
+++ b/miniapp/miniapp_recipes.cpp
@@ -17,7 +17,7 @@
namespace nest {
namespace mc {
-// TODO: split cell description into separate morphology, stimulus, probes, mechanisms etc.
+// TODO: split cell description into separate morphology, stimulus, mechanisms etc.
// description for greater data reuse.
template <typename RNG>
@@ -95,41 +95,56 @@ public:
return util::unique_any(dss_cell_description(spike_times));
}
- return util::unique_any(rss_cell::rss_cell_description(0.0, 0.1, 0.1));
+ return util::unique_any(rss_cell{0.0, 0.1, 0.1});
}
auto gen = std::mt19937(i); // TODO: replace this with hashing generator...
- auto cc = get_cell_count_info(i);
const auto& morph = get_morphology(i);
unsigned cell_segments = morph.components();
- auto cell = make_basic_cell(morph, param_.num_compartments, cc.num_targets,
+ auto cell = make_basic_cell(morph, param_.num_compartments, param_.num_synapses,
param_.synapse_type, gen);
EXPECTS(cell.num_segments()==cell_segments);
- EXPECTS(cell.probes().size()==0);
- EXPECTS(cell.synapses().size()==cc.num_targets);
- EXPECTS(cell.detectors().size()==cc.num_sources);
-
- // add probes
- if (cc.num_probes) {
- unsigned n_probe_segs = pdist_.all_segments? cell_segments: 1u;
- for (unsigned i = 0; i<n_probe_segs; ++i) {
- if (pdist_.membrane_voltage) {
- cell.add_probe({{i, i? 0.5: 0.0}, mc::probeKind::membrane_voltage});
- }
- if (pdist_.membrane_current) {
- cell.add_probe({{i, i? 0.5: 0.0}, mc::probeKind::membrane_current});
- }
- }
- }
- EXPECTS(cell.probes().size()==cc.num_probes);
+ EXPECTS(cell.synapses().size()==num_targets(i));
+ EXPECTS(cell.detectors().size()==num_sources(i));
return util::unique_any(std::move(cell));
}
- cell_kind get_cell_kind(cell_gid_type i ) const override {
+ probe_info get_probe(cell_member_type probe_id) const override {
+ if (probe_id.index>=num_probes(probe_id.gid)) {
+ throw invalid_recipe_error("invalid probe id");
+ }
+
+ // if we have both voltage and current probes, then order them
+ // voltage compartment 0, current compartment 0, voltage compartment 1, ...
+
+ cell_probe_address::probe_kind kind;
+
+ int stride = pdist_.membrane_voltage+pdist_.membrane_current;
+
+ if (stride==1) {
+ // Just one kind of probe.
+ kind = pdist_.membrane_voltage?
+ cell_probe_address::membrane_voltage: cell_probe_address::membrane_current;
+ }
+ else {
+ EXPECTS(stride==2);
+ // Both kinds available.
+ kind = (probe_id.index%stride==0)?
+ cell_probe_address::membrane_voltage: cell_probe_address::membrane_current;
+ }
+
+ cell_lid_type compartment = probe_id.index/stride;
+ segment_location loc{compartment, compartment? 0.5: 0.0};
+
+ // Use probe kind as the token to be passed to a sampler.
+ return {probe_id, (int)kind, cell_probe_address{loc, kind}};
+ }
+
+ cell_kind get_cell_kind(cell_gid_type i) const override {
// The last 'cell' is a rss_cell with one spike at t=0
if (i == ncell_) {
if (param_.input_spike_path) {
@@ -141,21 +156,24 @@ public:
return cell_kind::cable1d_neuron;
}
- cell_count_info get_cell_count_info(cell_gid_type i) const override {
- cell_count_info cc = {1, param_.num_synapses, 0 };
- unsigned cell_segments = get_morphology(i).components();
+ cell_size_type num_sources(cell_gid_type i) const override {
+ return 1;
+ }
- // probe this cell?
- if (std::floor(i*pdist_.proportion)!=std::floor((i-1.0)*pdist_.proportion)) {
- std::size_t np = pdist_.membrane_voltage + pdist_.membrane_current;
- if (pdist_.all_segments) {
- np *= cell_segments;
- }
+ cell_size_type num_targets(cell_gid_type i) const override {
+ return param_.num_synapses;
+ }
- cc.num_probes = np;
- }
+ cell_size_type num_probes(cell_gid_type i) const override {
+ bool has_probe = (std::floor(i*pdist_.proportion)!=std::floor((i-1.0)*pdist_.proportion));
- return cc;
+ if (!has_probe) {
+ return 0;
+ }
+ else {
+ cell_size_type np = pdist_.all_segments? get_morphology(i).components(): 1;
+ return np*(pdist_.membrane_voltage+pdist_.membrane_current);
+ }
}
protected:
diff --git a/miniapp/trace.cpp b/miniapp/trace.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..32903c757e92904cdeb00da6f2dba6c8fe95f4a7
--- /dev/null
+++ b/miniapp/trace.cpp
@@ -0,0 +1,50 @@
+#include <fstream>
+#include <string>
+
+#include <json/json.hpp>
+
+#include <common_types.hpp>
+#include <util/strprintf.hpp>
+
+#include "trace.hpp"
+
+using namespace nest::mc;
+
+static std::string to_string(cell_member_type m) {
+ return std::to_string(m.gid)+"."+std::to_string(m.index);
+}
+
+void write_trace_csv(const sample_trace& trace, const std::string& prefix) {
+ auto path = prefix + to_string(trace.probe_id) + "_" + trace.name + ".csv";
+
+ std::ofstream file(path);
+ file << "# cell: " << trace.probe_id.gid << "\n";
+ file << "# probe: " << trace.probe_id.index << "\n";
+ file << "time_ms, " << trace.name << "_" << trace.units << "\n";
+
+ for (const auto& sample: trace.samples) {
+ file << util::strprintf("% 20.15f, % 20.15f\n", sample.t, sample.v);
+ }
+}
+
+void write_trace_json(const sample_trace& trace, const std::string& prefix) {
+ auto path = prefix + to_string(trace.probe_id) + "_" + trace.name + ".json";
+
+ nlohmann::json jrep;
+ jrep["name"] = trace.name;
+ jrep["units"] = trace.units;
+ jrep["cell"] = trace.probe_id.gid;
+ jrep["probe"] = trace.probe_id.index;
+
+ auto& jt = jrep["data"]["time"];
+ auto& jy = jrep["data"][trace.name];
+
+ for (const auto& sample: trace.samples) {
+ jt.push_back(sample.t);
+ jy.push_back(sample.v);
+ }
+
+ std::ofstream file(path);
+ file << std::setw(1) << jrep << "\n";
+}
+
diff --git a/miniapp/trace.hpp b/miniapp/trace.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8bb5bf664001997c3eb9fb0c800a5d875f08e5ab
--- /dev/null
+++ b/miniapp/trace.hpp
@@ -0,0 +1,21 @@
+#pragma once
+
+/*
+ * Store trace data from samplers with metadata.
+ */
+
+#include <string>
+#include <vector>
+
+#include <common_types.hpp>
+#include <simple_sampler.hpp>
+
+struct sample_trace {
+ nest::mc::cell_member_type probe_id;
+ std::string name;
+ std::string units;
+ nest::mc::trace_data<double> samples;
+};
+
+void write_trace_csv(const sample_trace& trace, const std::string& prefix);
+void write_trace_json(const sample_trace& trace, const std::string& prefix);
diff --git a/miniapp/trace_sampler.hpp b/miniapp/trace_sampler.hpp
deleted file mode 100644
index 6bdc51a5bbf6f1c10bac2b8b5c5748faca888c18..0000000000000000000000000000000000000000
--- a/miniapp/trace_sampler.hpp
+++ /dev/null
@@ -1,75 +0,0 @@
-#pragma once
-
-/*
- * Simple(st?) implementation of a recorder of scalar
- * trace data from a cell probe, with some metadata.
- */
-
-#include <cstdlib>
-#include <vector>
-
-#include <common_types.hpp>
-#include <cell.hpp>
-#include <util/optional.hpp>
-
-#include <iostream>
-
-namespace nest {
-namespace mc {
-
-template <typename Time=float, typename Value=double>
-struct sample_trace {
- using time_type = Time;
- using value_type = Value;
-
- struct sample_type {
- time_type time;
- value_type value;
- };
-
- std::string name;
- std::string units;
- cell_member_type probe_id;
- std::vector<sample_type> samples;
-
- sample_trace() =default;
- sample_trace(cell_member_type probe_id, const std::string& name, const std::string& units):
- name(name), units(units), probe_id(probe_id)
- {}
-};
-
-template <typename Time=float, typename Value=double>
-struct trace_sampler {
- using time_type = Time;
- using value_type = Value;
-
- time_type next_sample_t() const { return t_next_sample_; }
-
- util::optional<time_type> operator()(time_type t, value_type v) {
- if (t<t_next_sample_) {
- return t_next_sample_;
- }
-
- trace_->samples.push_back({t,v});
- return t_next_sample_+=sample_dt_;
- }
-
- trace_sampler(sample_trace<time_type, value_type> *trace, time_type sample_dt, time_type tfrom=0):
- trace_(trace), sample_dt_(sample_dt), t_next_sample_(tfrom)
- {}
-
-private:
- sample_trace<time_type, value_type> *trace_;
-
- time_type sample_dt_;
- time_type t_next_sample_;
-};
-
-// with type deduction ...
-template <typename Time, typename Value>
-trace_sampler<Time, Value> make_trace_sampler(sample_trace<Time, Value> *trace, Time sample_dt, Time tfrom=0) {
- return trace_sampler<Time, Value>(trace, sample_dt, tfrom);
-}
-
-} // namespace mc
-} // namespace nest
diff --git a/scripts/tsplot b/scripts/tsplot
index 64186c90f0c529786b136743ff7216498146a2f6..390b5bfd266a35afa657a164875c6445a6b66055 100755
--- a/scripts/tsplot
+++ b/scripts/tsplot
@@ -57,6 +57,8 @@ def parse_clargs():
help='use colour COLOUR a base for series matching EXPR')
P.add_argument('-o', '--output', metavar='FILE', dest='outfile',
help='save plot to file FILE')
+ P.add_argument('-l', '--list', action='store_true',
+ help='list selected time-series')
P.add_argument('--dpi', metavar='NUM', dest='dpi',
type=int,
help='set dpi for output image')
@@ -446,19 +448,28 @@ for filename in args.inputs:
j = json.load(f)
tss.extend(read_json_timeseries(j, axis, select))
-if args.trange:
+if args.list:
for ts in tss:
- ts.trestrict(args.trange)
+ print 'name:', ts.meta['name']
+ print 'label:', ts.meta['label']
+ for k in [x for x in sorted(ts.meta.keys()) if x not in ['name', 'label']]:
+ print k+':', ts.meta[k]
+ print
-if args.exclude:
- for ts in tss:
- ts.exclude_outliers(args.exclude)
+else:
+ if args.trange:
+ for ts in tss:
+ ts.trestrict(args.trange)
+
+ if args.exclude:
+ for ts in tss:
+ ts.exclude_outliers(args.exclude)
-groupby = args.groupby if args.groupby else []
-plots = gather_ts_plots(tss, groupby)
+ groupby = args.groupby if args.groupby else []
+ plots = gather_ts_plots(tss, groupby)
-if not args.outfile:
- M.interactive(False)
+ if not args.outfile:
+ M.interactive(False)
-colours = args.colours if args.colours else []
-plot_plots(plots, axis=axis, colour_overrides=colours, save=args.outfile, dpi=args.dpi, scale=args.scale)
+ colours = args.colours if args.colours else []
+ plot_plots(plots, axis=axis, colour_overrides=colours, save=args.outfile, dpi=args.dpi, scale=args.scale)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index d189142b6d485c33290f65b478348eed127668a4..8429bac998cf01963b596fb50bda7a0efd70d101 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -16,6 +16,7 @@ set(BASE_SOURCES
profiling/meter_manager.cpp
profiling/power_meter.cpp
profiling/profiler.cpp
+ schedule.cpp
swcio.cpp
threading/threading.cpp
util/debug.cpp
diff --git a/src/cell.hpp b/src/cell.hpp
index fbcbe5047b4c636f16105ed5c71dc0490482884c..40139722691bce880508dcd59bcdb730532dc58a 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -8,7 +8,6 @@
#include <common_types.hpp>
#include <cell_tree.hpp>
#include <morphology.hpp>
-#include <probes.hpp>
#include <segment.hpp>
#include <stimulus.hpp>
#include <util/debug.hpp>
@@ -31,9 +30,13 @@ int find_compartment_index(
compartment_model const& graph
);
-struct probe_spec {
+struct cell_probe_address {
+ enum probe_kind {
+ membrane_voltage, membrane_current
+ };
+
segment_location location;
- probeKind kind;
+ probe_kind kind;
};
// used in constructor below
@@ -66,14 +69,12 @@ public:
// constructor
cell();
- // sometimes we really do want a copy (pending big morphology
- // refactor)
+ // Sometimes we really do want a copy (pending big morphology refactor).
cell(clone_cell_t, const cell& other):
parents_(other.parents_),
stimuli_(other.stimuli_),
synapses_(other.synapses_),
- spike_detectors_(other.spike_detectors_),
- probes_(other.probes_)
+ spike_detectors_(other.spike_detectors_)
{
// unique_ptr's cannot be copy constructed, do a manual assignment
segments_.reserve(other.segments_.size());
@@ -181,19 +182,6 @@ public:
return spike_detectors_;
}
- //////////////////
- // probes
- //////////////////
- index_type add_probe(probe_spec p) {
- probes_.push_back(p);
- return probes_.size()-1;
- }
-
- const std::vector<probe_spec>&
- probes() const {
- return probes_;
- }
-
private:
// storage for connections
std::vector<index_type> parents_;
@@ -209,9 +197,6 @@ private:
// the sensors
std::vector<detector_instance> spike_detectors_;
-
- // the probes
- std::vector<probe_spec> probes_;
};
// Checks that two cells have the same
diff --git a/src/cell_group.hpp b/src/cell_group.hpp
index 1473411bc2563a37c157a07adcf0f812a24fd2c7..7c9bdd3c449cc97a90f7bc917fbbe775859787b2 100644
--- a/src/cell_group.hpp
+++ b/src/cell_group.hpp
@@ -7,8 +7,8 @@
#include <common_types.hpp>
#include <event_binner.hpp>
#include <event_queue.hpp>
-#include <probes.hpp>
-#include <sampler_function.hpp>
+#include <sampling.hpp>
+#include <schedule.hpp>
#include <spike.hpp>
namespace nest {
@@ -26,8 +26,13 @@ public:
virtual void enqueue_events(const std::vector<postsynaptic_spike_event>& events) = 0;
virtual const std::vector<spike>& spikes() const = 0;
virtual void clear_spikes() = 0;
- virtual void add_sampler(cell_member_type probe_id, sampler_function s, time_type start_time = 0) = 0;
- virtual std::vector<probe_record> probes() const = 0;
+
+ // Sampler association methods below should be thread-safe, as they might be invoked
+ // from a sampler call back called from a different cell group running on a different thread.
+
+ virtual void add_sampler(sampler_association_handle, cell_member_predicate, schedule, sampler_function, sampling_policy) = 0;
+ virtual void remove_sampler(sampler_association_handle) = 0;
+ virtual void remove_all_samplers() = 0;
};
using cell_group_ptr = std::unique_ptr<cell_group>;
diff --git a/src/cell_group_factory.cpp b/src/cell_group_factory.cpp
index 857aa26035dc1f5a9ebc477c2f27968a62e95eae..e3cf5b7632d1906ae2729844a6e36ea52f47e5c4 100644
--- a/src/cell_group_factory.cpp
+++ b/src/cell_group_factory.cpp
@@ -17,28 +17,20 @@ using gpu_fvm_cell = mc_cell_group<fvm::fvm_multicell<gpu::backend>>;
using mc_fvm_cell = mc_cell_group<fvm::fvm_multicell<multicore::backend>>;
cell_group_ptr cell_group_factory(const recipe& rec, const group_description& group) {
- // Make a list of all the cell descriptions to be forwarded
- // to the appropriate cell_group constructor.
- std::vector<util::unique_any> descriptions;
- descriptions.reserve(group.gids.size());
- for (auto gid: group.gids) {
- descriptions.push_back(rec.get_cell_description(gid));
- }
-
switch (group.kind) {
case cell_kind::cable1d_neuron:
if (group.backend == backend_kind::gpu) {
- return make_cell_group<gpu_fvm_cell>(group.gids, descriptions);
+ return make_cell_group<gpu_fvm_cell>(group.gids, rec);
}
else {
- return make_cell_group<mc_fvm_cell>(group.gids, descriptions);
+ return make_cell_group<mc_fvm_cell>(group.gids, rec);
}
case cell_kind::regular_spike_source:
- return make_cell_group<rss_cell_group>(group.gids, descriptions);
+ return make_cell_group<rss_cell_group>(group.gids, rec);
case cell_kind::data_spike_source:
- return make_cell_group<dss_cell_group>(group.gids, descriptions);
+ return make_cell_group<dss_cell_group>(group.gids, rec);
default:
throw std::runtime_error("unknown cell kind");
diff --git a/src/common_types.hpp b/src/common_types.hpp
index 74bb33f6c6bf55c13de90e813c3cb7fe9bccc7c0..9a4530a82835873599791c2c80bc3f12c1ae2ca3 100644
--- a/src/common_types.hpp
+++ b/src/common_types.hpp
@@ -5,6 +5,8 @@
* library. (Expect analogues in future versions to be template parameters?)
*/
+#include <cstddef>
+#include <functional>
#include <iosfwd>
#include <type_traits>
@@ -53,6 +55,10 @@ DEFINE_LEXICOGRAPHIC_ORDERING(cell_member_type,(a.gid,a.index),(b.gid,b.index))
using time_type = float;
+// Extra contextual information associated with a probe.
+
+using probe_tag = int;
+
// Enumeration used to indentify the cell type/kind, used by the model to
// group equal kinds in the same cell group.
@@ -66,3 +72,16 @@ enum cell_kind {
} // namespace nest
std::ostream& operator<<(std::ostream& O, nest::mc::cell_member_type m);
+
+namespace std {
+ template <> struct hash<nest::mc::cell_member_type> {
+ std::size_t operator()(const nest::mc::cell_member_type& m) const {
+ using namespace nest::mc;
+ static_assert(sizeof(std::size_t)>sizeof(cell_gid_type), "invalid size assumptions for hash of cell_member_type");
+
+ std::size_t k = ((std::size_t)m.gid << (8*sizeof(cell_gid_type))) + m.index;
+ return std::hash<std::size_t>{}(k);
+ }
+ };
+}
+
diff --git a/src/dss_cell_description.hpp b/src/dss_cell_description.hpp
index 1fddebb75867071e916ed9b8141ea07c94091f41..bfb9fd68e2668b7aaf842dc88b68f3ea3c81f2ed 100644
--- a/src/dss_cell_description.hpp
+++ b/src/dss_cell_description.hpp
@@ -1,25 +1,23 @@
#pragma once
-#pragma once
#include <vector>
#include <common_types.hpp>
-#include <util/debug.hpp>
namespace nest {
namespace mc {
-/// Description for a data spike source: A cell that generates spikes provided as a vector of
-/// floating point valued spiketimes at the start of a run
+/// Description for a data spike source: a cell that generates spikes provided as a vector of
+/// spike times at the start of a run.
+
struct dss_cell_description {
std::vector<time_type> spike_times;
/// The description needs a vector of doubles for the description
- dss_cell_description(std::vector<time_type> & spike_times) :
- spike_times(spike_times)
+ dss_cell_description(std::vector<time_type> spike_times):
+ spike_times(std::move(spike_times))
{}
};
-
} // namespace mc
} // namespace nest
diff --git a/src/dss_cell_group.hpp b/src/dss_cell_group.hpp
index d51a4e1f7bd847239a2682b8e7de46dd250f56d3..5f3e32d68548c453d65fae9ce2372826797c480d 100644
--- a/src/dss_cell_group.hpp
+++ b/src/dss_cell_group.hpp
@@ -2,6 +2,7 @@
#include <cell_group.hpp>
#include <dss_cell_description.hpp>
+#include <recipe.hpp>
#include <util/span.hpp>
#include <util/unique_any.hpp>
@@ -11,24 +12,21 @@ namespace mc {
/// Cell_group to collect spike sources
class dss_cell_group: public cell_group {
public:
- dss_cell_group(std::vector<cell_gid_type> gids,
- const std::vector<util::unique_any>& cell_descriptions):
+ dss_cell_group(std::vector<cell_gid_type> gids, const recipe& rec):
gids_(std::move(gids))
{
- using util::make_span;
- for (cell_gid_type i: make_span(0, cell_descriptions.size())) {
+ for (auto gid: gids_) {
+ auto desc = util::any_cast<dss_cell_description>(rec.get_cell_description(gid));
// store spike times from description
- auto times = util::any_cast<dss_cell_description>(cell_descriptions[i]).spike_times;
+ auto times = desc.spike_times;
util::sort(times);
spike_times_.push_back(std::move(times));
// Take a reference to the first spike time
- not_emit_it_.push_back(spike_times_[i].begin());
+ not_emit_it_.push_back(spike_times_.back().begin());
}
}
- virtual ~dss_cell_group() = default;
-
cell_kind get_cell_kind() const override {
return cell_kind::data_spike_source;
}
@@ -45,8 +43,7 @@ public:
clear_spikes();
}
- void set_binning_policy(binning_kind policy, time_type bin_interval) override
- {}
+ void set_binning_policy(binning_kind policy, time_type bin_interval) override {}
void advance(time_type tfinal, time_type dt) override {
for (auto i: util::make_span(0, not_emit_it_.size())) {
@@ -78,14 +75,14 @@ public:
spikes_.clear();
}
- std::vector<probe_record> probes() const override {
- return probes_;
- }
-
- void add_sampler(cell_member_type probe_id, sampler_function s, time_type start_time = 0) override {
+ void add_sampler(sampler_association_handle h, cell_member_predicate probe_ids, schedule sched, sampler_function fn, sampling_policy policy) override {
std::logic_error("The dss_cells do not support sampling of internal state!");
}
+ void remove_sampler(sampler_association_handle h) override {}
+
+ void remove_all_samplers() override {}
+
private:
// Spikes that are generated.
std::vector<spike> spikes_;
@@ -93,8 +90,6 @@ private:
// Map of local index to gid
std::vector<cell_gid_type> gids_;
- std::vector<probe_record> probes_;
-
// The dss_cell is simple: Put all logic in the cellgroup cause accelerator support
// is not expected. We need storage for the cell state
diff --git a/src/fvm_multicell.hpp b/src/fvm_multicell.hpp
index 571b600d409b75617e272fff1ca83e8cb7c58aef..37aba8de4ff4805bddb2f3417088f7a2834a3cad 100644
--- a/src/fvm_multicell.hpp
+++ b/src/fvm_multicell.hpp
@@ -16,6 +16,8 @@
#include <matrix.hpp>
#include <memory/memory.hpp>
#include <profiling/profiler.hpp>
+#include <recipe.hpp>
+#include <sampler_map.hpp>
#include <segment.hpp>
#include <stimulus.hpp>
#include <util/meta.hpp>
@@ -69,11 +71,16 @@ public:
resting_potential_ = potential_mV;
}
- template <typename Cells, typename Targets, typename Probes>
+ // Set up data structures for a fixed collection of cells identified by `gids`
+ // with descriptions taken from the recipe `rec`.
+ //
+ // Lowered-cell specific handles for targets and probes are stored in the
+ // caller-provided vector `target_handles` and map `probe_map`.
void initialize(
- const Cells& cells, // collection of nest::mc::cell descriptions
- Targets& target_handles, // (write) where to store target handles
- Probes& probe_handles); // (write) where to store probe handles
+ const std::vector<cell_gid_type>& gids,
+ const recipe& rec,
+ std::vector<target_handle>& target_handles,
+ probe_association_map<probe_handle>& probe_map);
void reset();
@@ -226,8 +233,6 @@ public:
return it==mechanisms_.end() ? util::nothing: util::just(*it);
}
- std::size_t num_probes() const { return probes_.size(); }
-
//
// Threshold crossing interface.
// Used by calling code to perform spike detection
@@ -342,8 +347,6 @@ private:
/// the ion species
std::map<mechanisms::ionKind, ion> ions_;
- std::vector<std::pair<const array fvm_multicell::*, size_type>> probes_;
-
/// Compact representation of the control volumes into which a segment is
/// decomposed. Used to reconstruct the weights used to convert current
/// densities to currents for density channels.
@@ -520,13 +523,14 @@ fvm_multicell<Backend>::compute_cv_area_capacitance(
}
template <typename Backend>
-template <typename Cells, typename Targets, typename Probes>
void fvm_multicell<Backend>::initialize(
- const Cells& cells,
- Targets& target_handles,
- Probes& probe_handles)
+ const std::vector<cell_gid_type>& gids,
+ const recipe& rec,
+ std::vector<target_handle>& target_handles,
+ probe_association_map<probe_handle>& probe_map)
{
using memory::make_const_view;
+ using util::any_cast;
using util::assign_by;
using util::make_partition;
using util::make_span;
@@ -535,13 +539,19 @@ void fvm_multicell<Backend>::initialize(
using util::transform_view;
using util::subrange_view;
- // count total targets and probes for validation of handle container sizes
+ ncell_ = size(gids);
std::size_t targets_count = 0u;
- std::size_t probes_count = 0u;
- auto targets_size = size(target_handles);
- auto probes_size = size(probe_handles);
- ncell_ = size(cells);
+ // Take cell descriptions from recipe. These are used initially
+ // to count compartments for allocation of data structures, and
+ // then interrogated again for the details for each cell in turn.
+
+ std::vector<cell> cells;
+ cells.reserve(gids.size());
+ for (auto gid: gids) {
+ cells.push_back(any_cast<cell>(rec.get_cell_description(gid)));
+ }
+
auto cell_num_compartments =
transform_view(cells, [](const cell& c) { return c.num_compartments(); });
@@ -583,8 +593,7 @@ void fvm_multicell<Backend>::initialize(
// setup per-cell event stores.
events_ = util::make_unique<multi_event_stream>(ncell_);
- // create each cell:
- auto probe_hi = probe_handles.begin();
+ // Create each cell:
// Allocate scratch storage for calculating quantities used to build the
// linear system: these will later be copied into target-specific storage
@@ -611,6 +620,7 @@ void fvm_multicell<Backend>::initialize(
// - each probe, stimulus and detector is attached to its compartment.
for (auto i: make_span(0, ncell_)) {
const auto& c = cells[i];
+ auto gid = gids[i];
auto comp_ival = cell_comp_part[i];
auto graph = c.model();
@@ -643,8 +653,6 @@ void fvm_multicell<Backend>::initialize(
}
for (const auto& syn: c.synapses()) {
- EXPECTS(targets_count < targets_size);
-
const auto& name = syn.mechanism.name();
auto& map_entry = syn_mech_map[name];
@@ -693,22 +701,26 @@ void fvm_multicell<Backend>::initialize(
thresholds.push_back(detector.threshold);
}
- // record probe locations by index into corresponding state vector
- for (const auto& probe: c.probes()) {
- EXPECTS(probes_count < probes_size);
+ // Retrieve probe addresses and tags from recipe for this cell.
+ for (cell_lid_type j: make_span(0, rec.num_probes(gid))) {
+ probe_info pi = rec.get_probe({gid, j});
+ auto where = any_cast<cell_probe_address>(pi.address);
+
+ auto comp = comp_ival.first+find_cv_index(where.location, graph);
+ probe_handle handle;
- auto comp = comp_ival.first+find_cv_index(probe.location, graph);
- switch (probe.kind) {
- case probeKind::membrane_voltage:
- *probe_hi++ = {&fvm_multicell::voltage_, comp};
+ switch (where.kind) {
+ case cell_probe_address::membrane_voltage:
+ handle = {&fvm_multicell::voltage_, comp};
break;
- case probeKind::membrane_current:
- *probe_hi++ = {&fvm_multicell::current_, comp};
+ case cell_probe_address::membrane_current:
+ handle = {&fvm_multicell::current_, comp};
break;
default:
throw std::logic_error("unrecognized probeKind");
}
- ++probes_count;
+
+ probe_map.insert({pi.id, {handle, pi.tag}});
}
}
@@ -716,9 +728,6 @@ void fvm_multicell<Backend>::initialize(
threshold_watcher_ =
threshold_watcher(cv_to_cell_, time_, time_to_, voltage_, spike_detector_index, thresholds);
- // confirm user-supplied container probes were appropriately sized.
- EXPECTS(probes_size==probes_count);
-
// store the geometric information in target-specific containers
cv_areas_ = make_const_view(tmp_cv_areas);
@@ -778,6 +787,8 @@ void fvm_multicell<Backend>::initialize(
mech_to_cv_index[mech_name] = mech_cv_index;
}
+ target_handles.resize(targets_count);
+
// Create point (synapse) mechanisms.
for (auto& syni: syn_mech_map) {
size_type mech_id = mechanisms_.size();
@@ -806,9 +817,6 @@ void fvm_multicell<Backend>::initialize(
}
}
- // confirm user-supplied containers for targets are appropriately sized
- EXPECTS(targets_size==targets_count);
-
// build the ion species
for (auto ion : mechanisms::ion_kinds()) {
// find the compartment indexes of all compartments that have a
diff --git a/src/mc_cell_group.hpp b/src/mc_cell_group.hpp
index 4e25d15c88a3768caa62d1efa2d9bdae3241bb0a..5fa97cc2e6632ab1069ea0c0a3154e0c1d7a688d 100644
--- a/src/mc_cell_group.hpp
+++ b/src/mc_cell_group.hpp
@@ -13,9 +13,11 @@
#include <event_binner.hpp>
#include <event_queue.hpp>
#include <recipe.hpp>
-#include <sampler_function.hpp>
+#include <sampler_map.hpp>
+#include <sampling.hpp>
#include <spike.hpp>
#include <util/debug.hpp>
+#include <util/filter.hpp>
#include <util/partition.hpp>
#include <util/range.hpp>
#include <util/unique_any.hpp>
@@ -31,69 +33,38 @@ public:
using lowered_cell_type = LoweredCell;
using value_type = typename lowered_cell_type::value_type;
using size_type = typename lowered_cell_type::value_type;
- using source_id_type = cell_member_type;
mc_cell_group() = default;
- template <typename Cells>
- mc_cell_group(std::vector<cell_gid_type> gids, const Cells& cell_descriptions):
+ mc_cell_group(std::vector<cell_gid_type> gids, const recipe& rec):
gids_(std::move(gids))
{
// Build lookup table for gid to local index.
for (auto i: util::make_span(0, gids_.size())) {
- gid2lid_[gids_[i]] = i;
+ gid_index_map_[gids_[i]] = i;
}
- // Create lookup structure for probe and target ids.
- build_handle_partitions(cell_descriptions);
- std::size_t n_probes = probe_handle_divisions_.back();
+ // Create lookup structure for target ids.
+ build_target_handle_partition(rec);
std::size_t n_targets = target_handle_divisions_.back();
- std::size_t n_detectors = algorithms::sum(util::transform_view(
- cell_descriptions, [](const cell& c) { return c.detectors().size(); }));
- // Allocate space to store handles.
- target_handles_.resize(n_targets);
- probe_handles_.resize(n_probes);
+ // Pre-allocate space to store handles, probe map.
+ auto n_probes = util::sum_by(gids_, [&rec](cell_gid_type i) { return rec.num_probes(i); });
+ probe_map_.reserve(n_probes);
+ target_handles_.reserve(n_targets);
- lowered_.initialize(cell_descriptions, target_handles_, probe_handles_);
+ // Construct cell implementation, retrieving handles and maps.
+ lowered_.initialize(gids_, rec, target_handles_, probe_map_);
- // Create a list of the global identifiers for the spike sources.
- auto source_gid = gids_.begin();
- for (const auto& cell: cell_descriptions) {
- for (cell_lid_type lid=0u; lid<cell.detectors().size(); ++lid) {
- spike_sources_.push_back(source_id_type{*source_gid, lid});
+ // Create a list of the global identifiers for the spike sources
+ for (auto source_gid: gids_) {
+ for (cell_lid_type lid = 0; lid<rec.num_sources(source_gid); ++lid) {
+ spike_sources_.push_back({source_gid, lid});
}
- ++source_gid;
}
- EXPECTS(spike_sources_.size()==n_detectors);
-
- // Create the enumeration of probes attached to cells in this cell group.
- probes_.reserve(n_probes);
- for (auto i: util::make_span(0, cell_descriptions.size())){
- const cell_gid_type probe_gid = gids_[i];
- const auto probes_on_cell = cell_descriptions[i].probes();
- for (cell_lid_type lid: util::make_span(0, probes_on_cell.size())) {
- // Get the unique global identifier of this probe.
- cell_member_type id{probe_gid, lid};
-
- // Get the location and kind information of the probe.
- const auto p = probes_on_cell[lid];
-
- // Record the combined identifier and probe details.
- probes_.push_back(probe_record{id, p.location, p.kind});
- }
- }
- }
- mc_cell_group(std::vector<cell_gid_type> gids,
- const std::vector<util::unique_any>& cell_descriptions):
- mc_cell_group(
- std::move(gids),
- util::transform_view(cell_descriptions,
- [](const util::unique_any& c) -> const cell& {
- return util::any_cast<const cell&>(c);
- }))
- {}
+ spike_sources_.shrink_to_fit();
+ }
cell_kind get_cell_kind() const override {
return cell_kind::cable1d_neuron;
@@ -113,6 +84,7 @@ public:
void advance(time_type tfinal, time_type dt) override {
EXPECTS(lowered_.state_synchronized());
+ time_type tstart = lowered_.min_time();
// Bin pending events and enqueue on lowered state.
time_type ev_min_time = lowered_.max_time(); // (but we're synchronized here)
@@ -124,6 +96,34 @@ public:
lowered_.setup_integration(tfinal, dt);
+ // Set up sample event queue.
+ struct sampler_entry {
+ typename lowered_cell_type::probe_handle handle;
+ probe_tag tag;
+ cell_member_type probe_id;
+ sampler_function sampler;
+ };
+ std::vector<sampler_entry> samplers;
+
+ for (auto &assoc: sampler_map_) {
+ auto ts = assoc.sched.events(tstart, tfinal);
+ if (ts.empty()) {
+ continue;
+ }
+
+ for (auto p: assoc.probe_ids) {
+ EXPECTS(probe_map_.count(p)>0);
+ sample_event::size_type idx = samplers.size();
+ auto pinfo = probe_map_[p];
+
+ samplers.push_back({pinfo.handle, pinfo.tag, p, assoc.sampler});
+
+ for (auto t: ts) {
+ sample_events_.push({idx, t});
+ }
+ }
+ }
+
util::optional<time_type> first_sample_time = sample_events_.time_if_before(tfinal);
std::vector<sample_event> requeue_sample_events;
while (!lowered_.integration_complete()) {
@@ -136,20 +136,18 @@ public:
requeue_sample_events.clear();
while (auto m = sample_events_.pop_if_before(cell_max_time)) {
- auto& s = samplers_[m->sampler_index];
+ auto& s = samplers[m->sampler_index];
EXPECTS((bool)s.sampler);
- time_type cell_time = lowered_.time(gid2lid(s.cell_gid));
+ time_type cell_time = lowered_.time(gid_to_index(s.probe_id.gid));
if (cell_time<m->time) {
// This cell hasn't reached this sample time yet.
requeue_sample_events.push_back(*m);
}
else {
- auto next = s.sampler(cell_time, lowered_.probe(s.handle));
- if (next) {
- m->time = std::max(*next, cell_time);
- requeue_sample_events.push_back(*m);
- }
+ const double value = lowered_.probe(s.handle);
+ sample_record record = {cell_time, &value};
+ s.sampler(s.probe_id, s.tag, 1u, &record);
}
}
for (auto& ev: requeue_sample_events) {
@@ -201,36 +199,41 @@ public:
spikes_.clear();
}
- const std::vector<source_id_type>& spike_sources() const {
+ const std::vector<cell_member_type>& spike_sources() const {
return spike_sources_;
}
- void add_sampler(cell_member_type probe_id, sampler_function s, time_type start_time) override {
- auto handle = get_probe_handle(probe_id);
+ void add_sampler(sampler_association_handle h, cell_member_predicate probe_ids,
+ schedule sched, sampler_function fn, sampling_policy policy) override
+ {
+ std::vector<cell_member_type> probeset =
+ util::assign_from(util::filter(util::keys(probe_map_), probe_ids));
+
+ if (!probeset.empty()) {
+ sampler_map_.add(h, sampler_association{std::move(sched), std::move(fn), std::move(probeset)});
+ }
+ }
- using size_type = sample_event::size_type;
- auto sampler_index = size_type(samplers_.size());
- samplers_.push_back({handle, probe_id.gid, s});
- sampler_start_times_.push_back(start_time);
- sample_events_.push({sampler_index, start_time});
+ void remove_sampler(sampler_association_handle h) override {
+ sampler_map_.remove(h);
}
- std::vector<probe_record> probes() const override {
- return probes_;
+ void remove_all_samplers() override {
+ sampler_map_.clear();
}
private:
- // List of the gids of the cells in the group
+ // List of the gids of the cells in the group.
std::vector<cell_gid_type> gids_;
// Hash table for converting gid to local index
- std::unordered_map<cell_gid_type, cell_gid_type> gid2lid_;
+ std::unordered_map<cell_gid_type, cell_gid_type> gid_index_map_;
// The lowered cell state (e.g. FVM) of the cell.
lowered_cell_type lowered_;
// Spike detectors attached to the cell.
- std::vector<source_id_type> spike_sources_;
+ std::vector<cell_member_type> spike_sources_;
// Spikes that are generated.
std::vector<spike> spikes_;
@@ -243,73 +246,43 @@ private:
// Pending samples to be taken.
event_queue<sample_event> sample_events_;
- std::vector<time_type> sampler_start_times_;
// Handles for accessing lowered cell.
using target_handle = typename lowered_cell_type::target_handle;
std::vector<target_handle> target_handles_;
+ // Maps probe ids to probe handles (from lowered cell) and tags (from probe descriptions).
using probe_handle = typename lowered_cell_type::probe_handle;
- std::vector<probe_handle> probe_handles_;
-
- struct sampler_entry {
- typename lowered_cell_type::probe_handle handle;
- cell_gid_type cell_gid;
- sampler_function sampler;
- };
+ probe_association_map<probe_handle> probe_map_;
// Collection of samplers to be run against probes in this group.
- std::vector<sampler_entry> samplers_;
-
- // Lookup table for probe ids -> local probe handle indices.
- std::vector<std::size_t> probe_handle_divisions_;
+ sampler_association_map sampler_map_;
// Lookup table for target ids -> local target handle indices.
std::vector<std::size_t> target_handle_divisions_;
- // Enumeration of the probes that are attached to the cells in the cell group
- std::vector<probe_record> probes_;
-
// Build handle index lookup tables.
- template <typename Cells>
- void build_handle_partitions(const Cells& cells) {
- auto probe_counts =
- util::transform_view(cells, [](const cell& c) { return c.probes().size(); });
- auto target_counts =
- util::transform_view(cells, [](const cell& c) { return c.synapses().size(); });
-
- make_partition(probe_handle_divisions_, probe_counts);
- make_partition(target_handle_divisions_, target_counts);
- }
-
- // Use handle partition to get index from id.
- template <typename Divisions>
- std::size_t handle_partition_lookup(const Divisions& divisions, cell_member_type id) const {
- return divisions[gid2lid(id.gid)]+id.index;
- }
-
- // Get probe handle from probe id.
- probe_handle get_probe_handle(cell_member_type probe_id) const {
- return probe_handles_[handle_partition_lookup(probe_handle_divisions_, probe_id)];
+ void build_target_handle_partition(const recipe& rec) {
+ util::make_partition(target_handle_divisions_,
+ util::transform_view(gids_, [&rec](cell_gid_type i) { return rec.num_targets(i); }));
}
// Get target handle from target id.
- target_handle get_target_handle(cell_member_type target_id) const {
- return target_handles_[handle_partition_lookup(target_handle_divisions_, target_id)];
+ target_handle get_target_handle(cell_member_type id) const {
+ return target_handles_[target_handle_divisions_[gid_to_index(id.gid)]+id.index];
}
void reset_samplers() {
// clear all pending sample events and reset to start at time 0
sample_events_.clear();
- using size_type = sample_event::size_type;
- for(size_type i=0; i<samplers_.size(); ++i) {
- sample_events_.push({i, sampler_start_times_[i]});
+ for (auto &assoc: sampler_map_) {
+ assoc.sched.reset();
}
}
- cell_gid_type gid2lid(cell_gid_type gid) const {
- auto it = gid2lid_.find(gid);
- EXPECTS(it!=gid2lid_.end());
+ cell_gid_type gid_to_index(cell_gid_type gid) const {
+ auto it = gid_index_map_.find(gid);
+ EXPECTS(it!=gid_index_map_.end());
return it->second;
}
};
diff --git a/src/model.cpp b/src/model.cpp
index 0e84a1b8a8137f012fdf6995b5b39cca99daafe3..b74bdf6a6d57db3a4dbc9a5234e86b2b672a565a 100644
--- a/src/model.cpp
+++ b/src/model.cpp
@@ -1,11 +1,10 @@
-#include <model.hpp>
-
#include <vector>
#include <backends.hpp>
#include <cell_group.hpp>
#include <cell_group_factory.hpp>
#include <domain_decomposition.hpp>
+#include <model.hpp>
#include <recipe.hpp>
#include <util/span.hpp>
#include <util/unique_any.hpp>
@@ -32,11 +31,6 @@ model::model(const recipe& rec, const domain_decomposition& decomp):
PL(2);
});
- // Store probes.
- for (const auto& c: cell_groups_) {
- util::append(probes_, c->probes());
- }
-
// Allocate an empty queue buffer for each cell group
// These must be set initially to ensure that a queue is available for each
// cell group for the first time step.
@@ -148,17 +142,33 @@ time_type model::run(time_type tfinal, time_type dt) {
return t_;
}
-void model::attach_sampler(cell_member_type probe_id, sampler_function f, time_type tfrom) {
- // TODO: remove the gid_groups data structure completely when/if no longer needed
- // for the probe interface.
- auto it = gid_groups_.find(probe_id.gid);
- if (it!=gid_groups_.end()) {
- cell_groups_[it->second]->add_sampler(probe_id, f, tfrom);
- }
+sampler_association_handle model::add_sampler(cell_member_predicate probe_ids, schedule sched, sampler_function f, sampling_policy policy) {
+ sampler_association_handle h = sassoc_handles_.acquire();
+
+ threading::parallel_for::apply(0, cell_groups_.size(),
+ [&](std::size_t i) {
+ cell_groups_[i]->add_sampler(h, probe_ids, sched, f, policy);
+ });
+
+ return h;
}
-const std::vector<probe_record>& model::probes() const {
- return probes_;
+void model::remove_sampler(sampler_association_handle h) {
+ threading::parallel_for::apply(0, cell_groups_.size(),
+ [&](std::size_t i) {
+ cell_groups_[i]->remove_sampler(h);
+ });
+
+ sassoc_handles_.release(h);
+}
+
+void model::remove_all_samplers() {
+ threading::parallel_for::apply(0, cell_groups_.size(),
+ [&](std::size_t i) {
+ cell_groups_[i]->remove_all_samplers();
+ });
+
+ sassoc_handles_.clear();
}
std::size_t model::num_spikes() const {
diff --git a/src/model.hpp b/src/model.hpp
index 6fe5a106dc07e6c3c8f7b76e9b13ce06e4940233..9c0c1ed65f6931b263f1ded3e2793b3186b7c99d 100644
--- a/src/model.hpp
+++ b/src/model.hpp
@@ -1,5 +1,6 @@
#pragma once
+#include <mutex>
#include <vector>
#include <backends.hpp>
@@ -9,10 +10,10 @@
#include <communication/communicator.hpp>
#include <communication/global_policy.hpp>
#include <recipe.hpp>
-#include <sampler_function.hpp>
+#include <sampling.hpp>
#include <thread_private_spike_store.hpp>
#include <util/nop.hpp>
-#include <util/rangeutil.hpp>
+#include <util/handle_set.hpp>
#include <util/unique_any.hpp>
namespace nest {
@@ -29,9 +30,15 @@ public:
time_type run(time_type tfinal, time_type dt);
- void attach_sampler(cell_member_type probe_id, sampler_function f, time_type tfrom = 0);
+ // Note: sampler functions may be invoked from a different thread than that
+ // which called the `run` method.
- const std::vector<probe_record>& probes() const;
+ sampler_association_handle add_sampler(cell_member_predicate probe_ids,
+ schedule sched, sampler_function f, sampling_policy policy = sampling_policy::lax);
+
+ void remove_sampler(sampler_association_handle);
+
+ void remove_all_samplers();
std::size_t num_spikes() const;
@@ -39,7 +46,7 @@ public:
void set_binning_policy(binning_kind policy, time_type bin_interval);
// access cell_group directly
- // TODO: depricate. Currently used in some validation tests to inject
+ // TODO: deprecate. Currently used in some validation tests to inject
// events directly into a cell group. This should be done with a spiking
// neuron.
cell_group& group(int i);
@@ -57,7 +64,6 @@ private:
time_type t_ = 0.;
std::vector<cell_group_ptr> cell_groups_;
- std::vector<probe_record> probes_;
using event_queue_type = typename communicator_type::event_queue;
util::double_buffer<std::vector<event_queue_type>> event_queues_;
@@ -95,6 +101,9 @@ private:
std::vector<event_queue_type>& current_events() { return event_queues_.get(); }
std::vector<event_queue_type>& future_events() { return event_queues_.other(); }
+
+ // Sampler associations handles are managed by a helper class.
+ util::handle_set<sampler_association_handle> sassoc_handles_;
};
} // namespace mc
diff --git a/src/probes.hpp b/src/probes.hpp
deleted file mode 100644
index 36b27c3bf19227924fca9cb2a814de6911953837..0000000000000000000000000000000000000000
--- a/src/probes.hpp
+++ /dev/null
@@ -1,22 +0,0 @@
-#pragma once
-
-#include <cell.hpp>
-#include <morphology.hpp>
-#include <segment.hpp>
-
-namespace nest {
-namespace mc {
-
-enum class probeKind {
- membrane_voltage,
- membrane_current
-};
-
-struct probe_record {
- cell_member_type id;
- segment_location location;
- probeKind kind;
-};
-
-} // namespace mc
-} // namespace nest
diff --git a/src/recipe.hpp b/src/recipe.hpp
index b2b4f6b04bf45e49914a8c04eb5a15880253cec4..1eda5252aa935800396a6fc46cc201a0b612842e 100644
--- a/src/recipe.hpp
+++ b/src/recipe.hpp
@@ -2,18 +2,20 @@
#include <cstddef>
#include <memory>
+#include <unordered_map>
#include <stdexcept>
#include <cell.hpp>
+#include <common_types.hpp>
#include <util/unique_any.hpp>
namespace nest {
namespace mc {
-struct cell_count_info {
- cell_size_type num_sources;
- cell_size_type num_targets;
- cell_size_type num_probes;
+struct probe_info {
+ cell_member_type id;
+ probe_tag tag;
+ util::any address;
};
class invalid_recipe_error: public std::runtime_error {
@@ -50,53 +52,17 @@ struct cell_connection {
class recipe {
public:
- virtual cell_size_type num_cells() const =0;
+ virtual cell_size_type num_cells() const = 0;
- virtual util::unique_any get_cell_description(cell_gid_type) const =0;
+ virtual util::unique_any get_cell_description(cell_gid_type) const = 0;
virtual cell_kind get_cell_kind(cell_gid_type) const = 0;
- virtual cell_count_info get_cell_count_info(cell_gid_type) const =0;
- virtual std::vector<cell_connection> connections_on(cell_gid_type) const =0;
-};
-
-
-/*
- * Recipe consisting of a single, unconnected cell
- * is particularly simple. Note keeps a reference to
- * the provided cell, so be aware of life time issues.
- */
-
-class singleton_recipe: public recipe {
-public:
- singleton_recipe(const cell& the_cell): cell_(the_cell) {}
-
- cell_size_type num_cells() const override {
- return 1;
- }
-
- util::unique_any get_cell_description(cell_gid_type) const override {
- return util::unique_any(cell(clone_cell, cell_));
- }
-
- cell_kind get_cell_kind(cell_gid_type) const override {
- return cell_.get_cell_kind();
- }
-
- cell_count_info get_cell_count_info(cell_gid_type) const override {
- cell_count_info k;
- k.num_sources = cell_.detectors().size();
- k.num_targets = cell_.synapses().size();
- k.num_probes = cell_.probes().size();
-
- return k;
- }
-
- std::vector<cell_connection> connections_on(cell_gid_type) const override {
- return std::vector<cell_connection>{};
- }
+ virtual cell_size_type num_sources(cell_gid_type) const = 0;
+ virtual cell_size_type num_targets(cell_gid_type) const = 0;
+ virtual cell_size_type num_probes(cell_gid_type) const = 0;
-private:
- const cell& cell_;
+ virtual std::vector<cell_connection> connections_on(cell_gid_type) const = 0;
+ virtual probe_info get_probe(cell_member_type probe_id) const = 0;
};
} // namespace mc
diff --git a/src/rss_cell.hpp b/src/rss_cell.hpp
index 7f627ffc965d61733692ec665fa7bbb0146b580a..6ee5d07c4082a492aa13f5113b6bcb69b90665dd 100644
--- a/src/rss_cell.hpp
+++ b/src/rss_cell.hpp
@@ -1,94 +1,17 @@
#pragma once
-#pragma once
-
-#include <vector>
#include <common_types.hpp>
-#include <util/debug.hpp>
namespace nest {
namespace mc {
-/// Regular spike source: A cell that generated spikes with a certain
-/// period for a set time range.
-class rss_cell {
-public:
- using index_type = cell_lid_type;
- using size_type = cell_local_size_type;
- using value_type = double;
-
- struct rss_cell_description {
- time_type start_time;
- time_type period;
- time_type stop_time;
-
- rss_cell_description(time_type st, time_type per, time_type stop):
- start_time(st),
- period(per),
- stop_time(stop)
- {}
- };
-
- /// Construct a rss cell from its description
- rss_cell(rss_cell_description descr) :
- start_time_(descr.start_time),
- period_(descr.period),
- stop_time_(descr.stop_time),
- time_(0.0)
- {}
-
- /// Construct a rss cell from individual arguments
- rss_cell(time_type start_time=0.0, time_type period=1.0, time_type stop_time=0.0):
- start_time_(start_time),
- period_(period),
- stop_time_(stop_time),
- time_(0.0)
- {}
-
- /// Return the kind of cell, used for grouping into cell_groups
- cell_kind get_cell_kind() const {
- return cell_kind::regular_spike_source;
- }
-
- /// Collect all spikes until tfinal.
- // updates the internal time state to tfinal as a side effect
- std::vector<time_type> spikes_until(time_type tfinal) {
- std::vector<time_type> spike_times;
-
- // If we should be spiking in this 'period'
- if (tfinal > start_time_) {
- // We have to spike till tfinal or the stop_time_of the neuron
- auto end_time = stop_time_ < tfinal ? stop_time_ : tfinal;
- // Generate all possible spikes in this time frame typically only
- for (time_type time = start_time_ > time_ ? start_time_ : time_;
- time < end_time;
- time += period_) {
- spike_times.push_back(time);
- }
- }
- // Save our current time we generate exclusive a possible spike at tfinal
- time_ = tfinal;
-
- return spike_times;
- }
-
- /// reset internal time to 0.0
- void reset() {
- time_ = 0.0;
- }
-
-private:
- // When to start spiking
- time_type start_time_;
-
- // with what period
- time_type period_;
-
- // untill when
- time_type stop_time_;
+/// Description class for a regular spike source: a cell that generates
+/// spikes with a fixed period over a given time interval.
- // internal time, storage
- time_type time_;
+struct rss_cell {
+ time_type start_time;
+ time_type period;
+ time_type stop_time;
};
} // namespace mc
diff --git a/src/rss_cell_group.hpp b/src/rss_cell_group.hpp
index 6fc50d3fc1dcd4f0c60395a7de855f6ac0371242..201b7b9f8b58e5131503532ad049991fefec5292 100644
--- a/src/rss_cell_group.hpp
+++ b/src/rss_cell_group.hpp
@@ -1,61 +1,55 @@
#pragma once
+#include <utility>
+
#include <cell_group.hpp>
+#include <recipe.hpp>
#include <rss_cell.hpp>
-#include <util/optional.hpp>
-#include <util/span.hpp>
#include <util/unique_any.hpp>
-
-#include <iostream>
namespace nest {
namespace mc {
-/// Cell_group to collect cells that spike at a set frequency
-/// Cell are lightweight and are not executed in anybackend implementation
+/// Cell group implementing RSS cells.
+
class rss_cell_group: public cell_group {
public:
- using source_id_type = cell_member_type;
-
- rss_cell_group(std::vector<cell_gid_type> gids,
- const std::vector<util::unique_any>& cell_descriptions):
- gids_(gids)
- {
- using util::make_span;
-
- for (cell_gid_type i: make_span(0, cell_descriptions.size())) {
- cells_.push_back(rss_cell(
- util::any_cast<rss_cell::rss_cell_description>(cell_descriptions[i])
- ));
+ rss_cell_group(std::vector<cell_gid_type> gids, const recipe& rec) {
+ cells_.reserve(gids.size());
+ for (auto gid: gids) {
+ cells_.emplace_back(
+ util::any_cast<rss_cell>(rec.get_cell_description(gid)),
+ gid);
}
+ reset();
}
- virtual ~rss_cell_group() = default;
-
cell_kind get_cell_kind() const override {
return cell_kind::regular_spike_source;
}
void reset() override {
- for (auto cell: cells_) {
- cell.reset();
- }
+ clear_spikes();
+ time_ = 0;
}
- void set_binning_policy(binning_kind policy, time_type bin_interval) override
- {}
+ void set_binning_policy(binning_kind policy, time_type bin_interval) override {}
void advance(time_type tfinal, time_type dt) override {
- // TODO: Move rss_cell implementation into the rss_cell_group
- for (auto i: util::make_span(0, cells_.size())) {
- for (auto spike_time: cells_[i].spikes_until(tfinal)) {
- spikes_.push_back({{gids_[i], 0}, spike_time});
+ for (const auto& cell: cells_) {
+ auto t = std::max(cell.start_time, time_);
+ auto t_end = std::min(cell.stop_time, tfinal);
+
+ while (t < t_end) {
+ spikes_.push_back({{cell.gid, 0}, t});
+ t += cell.period;
}
}
- };
+ time_ = tfinal;
+ }
void enqueue_events(const std::vector<postsynaptic_spike_event>& events) override {
- std::logic_error("The rss_cells do not support incoming events!");
+ std::logic_error("rss_cell cannot deliver events");
}
const std::vector<spike>& spikes() const override {
@@ -66,23 +60,32 @@ public:
spikes_.clear();
}
- std::vector<probe_record> probes() const override {
- return {};
+ virtual void add_sampler(sampler_association_handle, cell_member_predicate, schedule, sampler_function, sampling_policy) {
+ std::logic_error("rss_cell does not support sampling");
}
- void add_sampler(cell_member_type, sampler_function, time_type ts=0) override {
- std::logic_error("The rss_cells do not support sampling of internal state!");
- }
+ void remove_sampler(sampler_association_handle) override {}
+
+ void remove_all_samplers() override {}
private:
- // List of the gids of the cells in the group
- std::vector<cell_gid_type> gids_;
+ // RSS description plus gid for each RSS cell.
+ struct rss_info: public rss_cell {
+ rss_info(rss_cell desc, cell_gid_type gid):
+ rss_cell(std::move(desc)), gid(gid)
+ {}
+
+ cell_gid_type gid;
+ };
+
+ // RSS cell descriptions.
+ std::vector<rss_info> cells_;
+
+ // Simulation time for all RSS cells in the group.
+ time_type time_;
// Spikes that are generated.
std::vector<spike> spikes_;
-
- // Store a reference to the cell actually implementing the spiking
- std::vector<rss_cell> cells_;
};
} // namespace mc
diff --git a/src/sampler_function.hpp b/src/sampler_function.hpp
deleted file mode 100644
index 6a042b585deb264a0378217b4ff10f0f9d1164df..0000000000000000000000000000000000000000
--- a/src/sampler_function.hpp
+++ /dev/null
@@ -1,14 +0,0 @@
-#pragma once
-
-#include <functional>
-
-#include <common_types.hpp>
-#include <util/optional.hpp>
-
-namespace nest {
-namespace mc {
-
-using sampler_function = std::function<util::optional<time_type>(time_type, double)>;
-
-} // namespace mc
-} // namespace nest
diff --git a/src/sampler_map.hpp b/src/sampler_map.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..34caadc7af34fdcd24fda895e484b5360846b750
--- /dev/null
+++ b/src/sampler_map.hpp
@@ -0,0 +1,77 @@
+#pragma once
+
+/*
+ * Helper classes for managing sampler/schedule associations in
+ * cell group classes (see sampling_api doc).
+ */
+
+#include <functional>
+#include <mutex>
+#include <unordered_map>
+
+#include <common_types.hpp>
+#include <sampling.hpp>
+#include <schedule.hpp>
+#include <util/deduce_return.hpp>
+#include <util/transform.hpp>
+
+namespace nest {
+namespace mc {
+
+// An association between a samplers, schedule, and set of probe ids, as provided
+// to e.g. `model::add_sampler()`.
+
+struct sampler_association {
+ schedule sched;
+ sampler_function sampler;
+ std::vector<cell_member_type> probe_ids;
+};
+
+// Maintain a set of associations paired with handles used for deletion.
+
+class sampler_association_map {
+public:
+ void add(sampler_association_handle h, sampler_association assoc) {
+ std::lock_guard<std::mutex> lock(m_);
+ map_.insert({h, std::move(assoc)});
+ }
+
+ void remove(sampler_association_handle h) {
+ std::lock_guard<std::mutex> lock(m_);
+ map_.erase(h);
+ }
+
+ void clear() {
+ std::lock_guard<std::mutex> lock(m_);
+ map_.clear();
+ }
+
+private:
+ using assoc_map = std::unordered_map<sampler_association_handle, sampler_association>;
+ assoc_map map_;
+ std::mutex m_;
+
+ static sampler_association& second(assoc_map::value_type& p) { return p.second; }
+ auto assoc_view() DEDUCED_RETURN_TYPE((util::transform_view(map_, &sampler_association_map::second)))
+
+public:
+ // Range-like view presents just the associations, omitting the handles.
+
+ auto begin() DEDUCED_RETURN_TYPE(assoc_view().begin());
+ auto end() DEDUCED_RETURN_TYPE(assoc_view().end());
+};
+
+// Manage associations between probe ids, probe tags, and (lowered cell) probe handles.
+
+template <typename Handle>
+struct probe_association {
+ using probe_handle_type = Handle;
+ probe_handle_type handle;
+ probe_tag tag;
+};
+
+template <typename Handle>
+using probe_association_map = std::unordered_map<cell_member_type, probe_association<Handle>>;
+
+} // namespace mc
+} // namespace nest
diff --git a/src/sampling.hpp b/src/sampling.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f1ca642e67a68ec6b655d61a8e9e4a3e0ddc3c69
--- /dev/null
+++ b/src/sampling.hpp
@@ -0,0 +1,36 @@
+#pragma once
+
+#include <cstddef>
+#include <functional>
+
+#include <common_types.hpp>
+#include <util/any_ptr.hpp>
+
+namespace nest {
+namespace mc {
+
+using cell_member_predicate = std::function<bool (cell_member_type)>;
+
+static cell_member_predicate all_probes = [](cell_member_type pid) { return true; };
+
+static inline cell_member_predicate one_probe(cell_member_type pid) {
+ return [pid](cell_member_type x) { return pid==x; };
+}
+
+struct sample_record {
+ time_type time;
+ util::any_ptr data;
+};
+
+using sampler_function = std::function<void (cell_member_type, probe_tag, std::size_t, const sample_record*)>;
+
+using sampler_association_handle = std::size_t;
+
+enum class sampling_policy {
+ lax,
+ // interpolated, // placeholder: unsupported
+ // exact // placeholder: unsupported
+};
+
+} // namespace mc
+} // namespace nest
diff --git a/src/schedule.cpp b/src/schedule.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c1f569a685c26ad040a2f4b5c1ce3fe16341d7f0
--- /dev/null
+++ b/src/schedule.cpp
@@ -0,0 +1,50 @@
+#include <algorithm>
+#include <iterator>
+#include <utility>
+#include <vector>
+
+#include <common_types.hpp>
+#include <schedule.hpp>
+
+// Implementations for specific schedules.
+
+namespace nest {
+namespace mc {
+
+// Regular schedule implementation.
+
+std::vector<time_type> regular_schedule_impl::events(time_type t0, time_type t1) {
+ std::vector<time_type> ts;
+
+ t0 = t0<0? 0: t0;
+ if (t1>t0) {
+ ts.reserve(1+std::size_t((t1-t0)*oodt_));
+
+ long long n = t0*oodt_;
+ time_type t = n*dt_;
+
+ while (t<t0) {
+ t = (++n)*dt_;
+ }
+
+ while (t<t1) {
+ ts.push_back(t);
+ t = (++n)*dt_;
+ }
+ }
+
+ return ts;
+}
+
+// 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);
+
+ start_index_ = std::distance(times_.begin(), ub);
+ return std::vector<time_type>(lb, ub);
+}
+
+} // namespace mc
+} // namespace nest
diff --git a/src/schedule.hpp b/src/schedule.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4a3212e0c975be82ad468f5713e6ad0a67d7b5e9
--- /dev/null
+++ b/src/schedule.hpp
@@ -0,0 +1,185 @@
+#pragma once
+
+#include <algorithm>
+#include <memory>
+#include <random>
+#include <vector>
+
+#include <common_types.hpp>
+#include <util/compat.hpp>
+#include <util/debug.hpp>
+#include <util/meta.hpp>
+
+// Time schedules for probe–sampler associations.
+
+namespace nest {
+namespace mc {
+
+// 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()`,
+// then 0 ≤ _t0_ ≤ _t1_ ≤ _t2_ ≤ _t3_.
+
+class schedule {
+public:
+ template <typename Impl>
+ explicit schedule(const Impl& impl):
+ impl_(new wrap<Impl>(impl)) {}
+
+ template <typename Impl>
+ explicit schedule(Impl&& impl):
+ impl_(new wrap<Impl>(std::move(impl))) {}
+
+ schedule(schedule&& other) = default;
+ schedule& operator=(schedule&& other) = default;
+
+ schedule(const schedule& other):
+ impl_(other.impl_->clone()) {}
+
+ schedule& operator=(const schedule& other) {
+ impl_ = other.impl_->clone();
+ return *this;
+ }
+
+ std::vector<time_type> events(time_type t0, time_type t1) {
+ return impl_->events(t0, t1);
+ }
+
+ void reset() { impl_->reset(); }
+
+private:
+ struct interface {
+ virtual std::vector<time_type> events(time_type t0, time_type t1) = 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)) {}
+
+ virtual std::vector<time_type> events(time_type t0, time_type t1) {
+ return wrapped.events(t0, t1);
+ }
+
+ virtual void reset() {
+ wrapped.reset();
+ }
+
+ virtual std::unique_ptr<interface> clone() {
+ return std::unique_ptr<interface>(new wrap<Impl>(wrapped));
+ }
+
+ Impl wrapped;
+ };
+};
+
+
+// Common schedules
+
+// Schedule at k·dt for integral k≥0.
+class regular_schedule_impl {
+public:
+ explicit regular_schedule_impl(time_type dt):
+ dt_(dt), oodt_(1./dt) {};
+
+ void reset() {}
+ std::vector<time_type> events(time_type t0, time_type t1);
+
+private:
+ time_type dt_;
+ time_type oodt_;
+};
+
+inline schedule regular_schedule(time_type dt) {
+ return schedule(regular_schedule_impl(dt));
+}
+
+// Schedule at times given explicitly via a provided sorted sequence.
+class explicit_schedule_impl {
+public:
+ template <typename Seq, typename = util::enable_if_sequence_t<Seq>>
+ explicit explicit_schedule_impl(const Seq& seq):
+ start_index_(0),
+ times_(std::begin(seq), compat::end(seq))
+ {
+ EXPECTS(std::is_sorted(times_.begin(), times_.end()));
+ }
+
+ void reset() {
+ start_index_ = 0;
+ }
+
+ std::vector<time_type> events(time_type t0, time_type t1);
+
+private:
+ std::ptrdiff_t start_index_;
+ std::vector<time_type> times_;
+};
+
+template <typename Seq>
+inline schedule explicit_schedule(const Seq& 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)
+ {
+ EXPECTS(tstart_>=0);
+ step();
+ }
+
+ void reset() {
+ rng_ = reset_state_;
+ next_ = tstart_;
+ step();
+ }
+
+ std::vector<time_type> events(time_type t0, time_type t1) {
+ std::vector<time_type> ts;
+
+ while (next_<t0) {
+ step();
+ }
+
+ while (next_<t1) {
+ ts.push_back(next_);
+ step();
+ }
+
+ return ts;
+ }
+
+private:
+ void step() {
+ next_ += exp_(rng_);
+ }
+
+ time_type tstart_;
+ std::exponential_distribution<time_type> exp_;
+ RandomNumberEngine rng_;
+ RandomNumberEngine reset_state_;
+ time_type next_;
+};
+
+template <typename RandomNumberEngine>
+inline schedule poisson_schedule(time_type mean_dt, const RandomNumberEngine& rng) {
+ return schedule(poisson_schedule_impl<RandomNumberEngine>(0., mean_dt, 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));
+}
+
+} // namespace mc
+} // namespace nest
diff --git a/src/simple_sampler.hpp b/src/simple_sampler.hpp
index 5ffc41d2470dbd55f0b1af732c1f0aa43db11530..e096756414d9813b8e0ed0e9caa79f343f40fd7c 100644
--- a/src/simple_sampler.hpp
+++ b/src/simple_sampler.hpp
@@ -5,73 +5,74 @@
* trace data from a cell probe, with some metadata.
*/
-#include <functional>
#include <vector>
#include <common_types.hpp>
-#include <sampler_function.hpp>
-#include <util/optional.hpp>
+#include <sampling.hpp>
+#include <util/any_ptr.hpp>
#include <util/deduce_return.hpp>
+#include <util/span.hpp>
#include <util/transform.hpp>
+#include <iostream>
+
namespace nest {
namespace mc {
+template <typename V>
struct trace_entry {
- float t;
- double v;
+ time_type t;
+ V v;
};
-using trace_data = std::vector<trace_entry>;
+template <typename V>
+using trace_data = std::vector<trace_entry<V>>;
// NB: work-around for lack of function return type deduction
// in C++11; can't use lambda within DEDUCED_RETURN_TYPE.
namespace impl {
- inline float time(const trace_entry& x) { return x.t; }
- inline float value(const trace_entry& x) { return x.v; }
+ template <typename V>
+ inline float time(const trace_entry<V>& x) { return x.t; }
+
+ template <typename V>
+ inline const V& value(const trace_entry<V>& x) { return x.v; }
}
-inline auto times(const trace_data& trace) DEDUCED_RETURN_TYPE(
- util::transform_view(trace, impl::time)
+template <typename V>
+inline auto times(const trace_data<V>& trace) DEDUCED_RETURN_TYPE(
+ util::transform_view(trace, impl::time<V>)
)
-inline auto values(const trace_data& trace) DEDUCED_RETURN_TYPE(
- util::transform_view(trace, impl::value)
+template <typename V>
+inline auto values(const trace_data<V>& trace) DEDUCED_RETURN_TYPE(
+ util::transform_view(trace, impl::value<V>)
)
+template <typename V, typename = util::enable_if_trivially_copyable_t<V>>
class simple_sampler {
public:
- trace_data trace;
+ explicit simple_sampler(trace_data<V>& trace): trace_(trace) {}
- simple_sampler(time_type dt, time_type t0=0):
- t0_(t0),
- sample_dt_(dt),
- t_next_sample_(t0)
- {}
-
- void reset() {
- trace.clear();
- t_next_sample_ = t0_;
- }
-
- sampler_function sampler() {
- return [&](time_type t, double v) -> util::optional<time_type> {
- if (t<t_next_sample_) {
- return t_next_sample_;
+ void operator()(cell_member_type probe_id, probe_tag tag, std::size_t n, const sample_record* recs) {
+ for (auto i: util::make_span(0, n)) {
+ if (auto p = util::any_cast<const V*>(recs[i].data)) {
+ trace_.push_back({recs[i].time, *p});
}
else {
- trace.push_back({t, v});
- return t_next_sample_+=sample_dt_;
+ throw std::runtime_error("unexpected sample type in simple_sampler");
}
- };
+ }
}
private:
- time_type t0_ = 0;
- time_type sample_dt_ = 0;
- time_type t_next_sample_ = 0;
+ trace_data<V>& trace_;
};
+template <typename V>
+inline simple_sampler<V> make_simple_sampler(trace_data<V>& trace) {
+ return simple_sampler<V>(trace);
+}
+
} // namespace mc
} // namespace nest
diff --git a/src/util/any_ptr.hpp b/src/util/any_ptr.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cbd7433cf691e7519d8a174cf9db4713bdc6c592
--- /dev/null
+++ b/src/util/any_ptr.hpp
@@ -0,0 +1,98 @@
+#pragma once
+
+/* Specialied type erasure for pointer types.
+ *
+ * `any_ptr` represents a non-owning pointer to an arbitrary type
+ * that can be confirmed at run-time.
+ *
+ * Semantics:
+ *
+ * 1. An `any_ptr` value p represents either a null pointer, or
+ * a non-null pointer of a specific but arbitrary type T.
+ *
+ * 2. The value of the pointer as a `void*` value can be retrieved
+ * with the member function `as<void*>()`.
+ *
+ * 3. The value of the pointer as a type T which is not `void*` is
+ * retrieved with the member function `as<T>()`. If the represented
+ * pointer is the null pointer or a pointer to a different type,
+ * `as<T>()` will return the null pointer.
+ */
+
+#include <cstddef>
+#include <type_traits>
+
+#include <util/lexcmp_def.hpp>
+#include <util/meta.hpp>
+
+namespace nest {
+namespace mc {
+namespace util {
+
+struct any_ptr {
+ any_ptr() {}
+
+ any_ptr(std::nullptr_t) {}
+
+ template <typename T>
+ any_ptr(T* ptr):
+ ptr_((void *)ptr), type_ptr_(&typeid(T*)) {}
+
+ const std::type_info& type() const noexcept { return *type_ptr_; }
+
+ bool has_value() const noexcept { return ptr_; }
+
+ operator bool() const noexcept { return has_value(); }
+
+ void reset() noexcept { ptr_ = nullptr; }
+
+ void reset(std::nullptr_t) noexcept { ptr_ = nullptr; }
+
+ template <typename T>
+ void reset(T* ptr) noexcept {
+ ptr_ = (void*)ptr;
+ type_ptr_ = &typeid(T*);
+ }
+
+ template <typename T, typename = enable_if_t<std::is_pointer<T>::value>>
+ T as() const noexcept {
+ if (std::is_same<T, void*>::value) {
+ return (T)ptr_;
+ }
+ else {
+ return typeid(T)==type()? (T)ptr_: nullptr;
+ }
+ }
+
+ any_ptr& operator=(const any_ptr& other) noexcept {
+ type_ptr_ = other.type_ptr_;
+ ptr_ = other.ptr_;
+ return *this;
+ }
+
+ any_ptr& operator=(std::nullptr_t) noexcept {
+ reset();
+ return *this;
+ }
+
+ template <typename T>
+ any_ptr& operator=(T* ptr) noexcept {
+ reset(ptr);
+ return *this;
+ }
+
+private:
+ void* ptr_ = nullptr;
+ const std::type_info* type_ptr_ = &typeid(void);
+};
+
+// Order, compare by pointer value:
+DEFINE_LEXICOGRAPHIC_ORDERING_BY_VALUE(any_ptr, (a.as<void*>()), (b.as<void*>()))
+
+// Overload `util::any_cast` for these pointers.
+template <typename T>
+T any_cast(any_ptr p) noexcept { return p.as<T>(); }
+
+} // namespace util
+} // namespace mc
+} // namespace nest
diff --git a/src/util/handle_set.hpp b/src/util/handle_set.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..04086c1ea07ec71a2fa79c7e6a8e75133848470e
--- /dev/null
+++ b/src/util/handle_set.hpp
@@ -0,0 +1,63 @@
+#pragma once
+
+#include <mutex>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+/*
+ * Manage a set of integer-valued handles.
+ *
+ * An error is thrown if no unused handle value is available on acquire.
+ *
+ * Simple implementation below does not try hard to reuse released
+ * handles; smarter versions to be implemented as required.
+ */
+
+namespace nest {
+namespace mc {
+namespace util {
+
+template <typename Handle>
+class handle_set {
+public:
+ using value_type = Handle;
+
+ value_type acquire() {
+ lock_guard lock(mex_);
+
+ if (top_==std::numeric_limits<Handle>::max()) {
+ throw std::out_of_range("no more handles");
+ }
+ return top_++;
+ }
+
+ // Pre-requisite: h is a handle returned by
+ // `acquire`, which has not been subject
+ // to a subsequent `release`.
+ void release(value_type h) {
+ lock_guard lock(mex_);
+
+ if (h+1==top_) {
+ --top_;
+ }
+ }
+
+ // Release all handles.
+ void clear() {
+ lock_guard lock(mex_);
+
+ top_ = 0;
+ }
+
+private:
+ value_type top_ = 0;
+
+ using lock_guard = std::lock_guard<std::mutex>;
+ std::mutex mex_;
+};
+
+} // namespace util
+} // namespace mc
+} // namespace nest
+
diff --git a/src/util/meta.hpp b/src/util/meta.hpp
index 35b1b3e3363c14b6fcb9f0162c0845b748667ea3..8984bc6b9df6c8a877674cd9bfb8f2f4b365fee1 100644
--- a/src/util/meta.hpp
+++ b/src/util/meta.hpp
@@ -123,6 +123,10 @@ template <typename T>
using enable_if_move_assignable_t =
enable_if_t<std::is_move_assignable<T>::value>;
+template <typename T>
+using enable_if_trivially_copyable_t =
+ enable_if_t<std::is_trivially_copyable<T>::value>;
+
// Iterator class test
// (might not be portable before C++17)
@@ -271,8 +275,6 @@ struct second_t {
};
constexpr second_t second{};
-
-
} // namespace util
} // namespace mc
} // namespace nest
diff --git a/src/util/rangeutil.hpp b/src/util/rangeutil.hpp
index 8d03be492d52ade3b90c5007e1bd83919cc2636c..b5d4d040ebbc848e87f1b44d0926ca85602a5594 100644
--- a/src/util/rangeutil.hpp
+++ b/src/util/rangeutil.hpp
@@ -10,6 +10,7 @@
#include <ostream>
#include <numeric>
+#include <util/deduce_return.hpp>
#include <util/meta.hpp>
#include <util/range.hpp>
#include <util/transform.hpp>
@@ -314,6 +315,12 @@ std::pair<Value, Value> minmax_value(const Seq& seq, Compare cmp = Compare{}) {
return {lower, upper};
}
+// View over the keys in an associative container.
+
+template <typename Map>
+auto keys(Map& m) DEDUCED_RETURN_TYPE(util::transform_view(m, util::first));
+
+
template <typename C, typename Seq>
C make_copy(Seq const& seq) {
return C{std::begin(seq), std::end(seq)};
diff --git a/src/util/transform.hpp b/src/util/transform.hpp
index da9165cf085cde1766f651ab3a3f8f8a5eedb6c5..ed9a31d6bc2aab70ae61a898fb2e65e0b4d05905 100644
--- a/src/util/transform.hpp
+++ b/src/util/transform.hpp
@@ -35,7 +35,7 @@ class transform_iterator: public iterator_adaptor<transform_iterator<I, F>, I> {
const I& inner() const { return inner_; }
I& inner() { return inner_; }
- using inner_value_type = util::decay_t<decltype(*inner_)>;
+ using inner_value_type = decltype(*inner_);
using raw_value_type = typename std::result_of<F (inner_value_type)>::type;
static constexpr bool present_lvalue = std::is_reference<raw_value_type>::value;
@@ -129,25 +129,25 @@ transform_iterator<I, util::decay_t<F>> make_transform_iterator(const I& i, cons
template <
typename Seq,
typename F,
- typename seq_citer = typename sequence_traits<Seq>::const_iterator,
- typename seq_csent = typename sequence_traits<Seq>::const_sentinel,
- typename = enable_if_t<std::is_same<seq_citer, seq_csent>::value>
+ typename seq_iter = typename sequence_traits<Seq>::iterator,
+ typename seq_sent = typename sequence_traits<Seq>::sentinel,
+ typename = enable_if_t<std::is_same<seq_iter, seq_sent>::value>
>
-range<transform_iterator<seq_citer, util::decay_t<F>>>
-transform_view(const Seq& s, const F& f) {
- return {make_transform_iterator(util::cbegin(s), f), make_transform_iterator(util::cend(s), f)};
+range<transform_iterator<seq_iter, util::decay_t<F>>>
+transform_view(Seq&& s, const F& f) {
+ return {make_transform_iterator(std::begin(s), f), make_transform_iterator(std::end(s), f)};
}
template <
typename Seq,
typename F,
- typename seq_citer = typename sequence_traits<Seq>::const_iterator,
- typename seq_csent = typename sequence_traits<Seq>::const_sentinel,
- typename = enable_if_t<!std::is_same<seq_citer, seq_csent>::value>
+ typename seq_iter = typename sequence_traits<Seq>::iterator,
+ typename seq_sent = typename sequence_traits<Seq>::sentinel,
+ typename = enable_if_t<!std::is_same<seq_iter, seq_sent>::value>
>
-range<transform_iterator<seq_citer, util::decay_t<F>>, seq_csent>
-transform_view(const Seq& s, const F& f) {
- return {make_transform_iterator(util::cbegin(s), f), util::cend(s)};
+range<transform_iterator<seq_iter, util::decay_t<F>>, seq_sent>
+transform_view(Seq&& s, const F& f) {
+ return {make_transform_iterator(std::begin(s), f), std::end(s)};
}
} // namespace util
diff --git a/tests/test_common_cells.hpp b/tests/common_cells.hpp
similarity index 94%
rename from tests/test_common_cells.hpp
rename to tests/common_cells.hpp
index e5f1a49bc40147822ad9e22e18c26bb87d63ecbe..feae376db2b26534fea3f97413b961ce30cb0501 100644
--- a/tests/test_common_cells.hpp
+++ b/tests/common_cells.hpp
@@ -1,6 +1,7 @@
#include <cmath>
#include <cell.hpp>
+#include <recipe.hpp>
#include <segment.hpp>
#include <math.hpp>
#include <parameter_list.hpp>
@@ -297,22 +298,5 @@ inline cell make_cell_simple_cable(bool with_stim = true) {
return c;
}
-
-/*
- * Attach voltage probes at each cable mid-point and end-point,
- * and at soma mid-point.
- */
-
-inline cell& add_common_voltage_probes(cell& c) {
- auto ns = c.num_segments();
- for (auto i=0u; i<ns; ++i) {
- c.add_probe({{i, 0.5}, probeKind::membrane_voltage});
- if (i>0) {
- c.add_probe({{i, 1.0}, probeKind::membrane_voltage});
- }
- }
- return c;
-}
-
} // namespace mc
} // namespace nest
diff --git a/tests/simple_recipes.hpp b/tests/simple_recipes.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0dd687362b09958c4caf8fc53540ec25665242c0
--- /dev/null
+++ b/tests/simple_recipes.hpp
@@ -0,0 +1,114 @@
+#pragma once
+
+// Simple recipe classes for use in unit and validation tests.
+
+#include <unordered_map>
+#include <vector>
+
+#include <cell.hpp>
+#include <recipe.hpp>
+
+namespace nest {
+namespace mc {
+
+// Common functionality: maintain an unordered map of probe data
+// per gid, built with `add_probe()`.
+
+class simple_recipe_base: public recipe {
+public:
+ cell_size_type num_probes(cell_gid_type i) const override {
+ return probes_.count(i)? probes_.at(i).size(): 0;
+ }
+
+ std::vector<cell_connection> connections_on(cell_gid_type) const override {
+ return {};
+ }
+
+ virtual probe_info get_probe(cell_member_type probe_id) const override {
+ return probes_.at(probe_id.gid).at(probe_id.index);
+ }
+
+ virtual void add_probe(cell_gid_type gid, probe_tag tag, util::any address) {
+ auto& pvec_ = probes_[gid];
+
+ cell_member_type probe_id{gid, cell_lid_type(pvec_.size())};
+ pvec_.push_back({probe_id, tag, std::move(address)});
+ }
+
+protected:
+ std::unordered_map<cell_gid_type, std::vector<probe_info>> probes_;
+};
+
+// Convenience derived recipe class for wrapping n copies of a single
+// cell description, with no sources or targets. (Derive the class to
+// add sources, targets, connections etc.)
+//
+// Probes are simply stored in a multimap, keyed by gid and can be manually
+// added with 'add_probe()'.
+//
+// Wrapped description class must be both move- and copy-constructable.
+
+template <cell_kind Kind, typename Description>
+class homogeneous_recipe: public simple_recipe_base {
+public:
+ homogeneous_recipe(cell_size_type n, Description desc):
+ n_(n), desc_(std::move(desc))
+ {}
+
+ cell_size_type num_cells() const override { return n_; }
+ cell_kind get_cell_kind(cell_gid_type) const override { return Kind; }
+
+ cell_size_type num_sources(cell_gid_type) const override { return 0; }
+ cell_size_type num_targets(cell_gid_type) const override { return 0; }
+
+ util::unique_any get_cell_description(cell_gid_type) const override {
+ return util::make_unique_any<Description>(desc_);
+ }
+
+protected:
+ cell_size_type n_;
+ Description desc_;
+};
+
+// Recipe for a set of `cable1d_neuron` neurons without connections,
+// and probes which can be added by `add_probe()` (similar to above).
+//
+// Cell descriptions passed to the constructor are cloned.
+
+class cable1d_recipe: public simple_recipe_base {
+public:
+ template <typename Seq>
+ explicit cable1d_recipe(const Seq& cells) {
+ for (const auto& c: cells) {
+ cells_.emplace_back(clone_cell, c);
+ }
+ }
+
+ explicit cable1d_recipe(const cell& c) {
+ cells_.reserve(1);
+ cells_.emplace_back(clone_cell, c);
+ }
+
+ cell_size_type num_cells() const override { return cells_.size(); }
+ cell_kind get_cell_kind(cell_gid_type) const override { return cell_kind::cable1d_neuron; }
+
+ cell_size_type num_sources(cell_gid_type i) const override {
+ return cells_.at(i).detectors().size();
+ }
+
+ cell_size_type num_targets(cell_gid_type i) const override {
+ return cells_.at(i).synapses().size();
+ }
+
+ util::unique_any get_cell_description(cell_gid_type i) const override {
+ return util::make_unique_any<cell>(clone_cell, cells_[i]);
+ }
+
+protected:
+ std::vector<cell> cells_;
+};
+
+
+} // namespace mc
+} // namespace nest
+
diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt
index 40427c1ab001a231e1001cafa096b24be472e3f8..61b36aa2a2b79acc7a7ec36e44d206873807dafc 100644
--- a/tests/unit/CMakeLists.txt
+++ b/tests/unit/CMakeLists.txt
@@ -61,12 +61,14 @@ set(TEST_SOURCES
test_path.cpp
test_point.cpp
test_probe.cpp
- test_segment.cpp
test_range.cpp
+ test_segment.cpp
+ test_schedule.cpp
test_rss_cell.cpp
test_span.cpp
test_spikes.cpp
test_spike_store.cpp
+ test_stats.cpp
test_stimulus.cpp
test_strprintf.cpp
test_swcio.cpp
@@ -79,6 +81,9 @@ set(TEST_SOURCES
# unit test driver
test.cpp
+
+ # common routines
+ stats.cpp
)
diff --git a/tests/unit/stats.cpp b/tests/unit/stats.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0e4a114df63e378b7257594355942aeb54ec1889
--- /dev/null
+++ b/tests/unit/stats.cpp
@@ -0,0 +1,157 @@
+#include <cmath>
+#include <numeric>
+#include <vector>
+
+#include "stats.hpp"
+
+namespace testing {
+namespace ks {
+
+/* Compute the exact cdf F(d) of the one-sample, two-sided Kolmogorov–Smirnov
+ * statisitc Dn for n observations, using the Dubin 1973 algorithm with the
+ * faster implementation of Carvalho 2015.
+ *
+ * Bounds on numerical applicability, and the lower bound of 1-10^-15 (giving a
+ * result of 1 in double precision) when n*d^2 ≥ 18 is from Simard and L'Ecuyer
+ * 2011.
+ *
+ * This function should implement the Pelz-Good asymptopic approximation, but
+ * doesn't. Instead it will throw an error if the parameters fall into a bad
+ * region.
+ *
+ * References:
+ *
+ * Durbin, J. 1973. "Distribution Theory for Tests Based on the Sample
+ * Distribution Function". Society for Industrial Mathematics.
+ * ISBN: 9780898710076
+ * DOI: 10.1137/1.9781611970568.ch1
+ *
+ * Carvalho, L. 2015. "An improved evaluation of Kolmogorov's distribution."
+ * Journal of Statistical Software, Code Samples 65 (3), pp. 1-8.
+ * DOI: 10.18637/jss.v065.c03
+ *
+ * Simard, R. and P. L'Ecuyer, 2011. "Computing the two-sided Kolmogorov-Smirnov distributiobnn."
+ * Journal of Statistical Softwarem, Articles 39 (11), pp. 1-18.
+ * DOI: 10.18637/jss.v039.i11
+ */
+
+double dn_cdf(double d, int n) {
+ // Tail cases:
+
+ double nd = n*d;
+
+ if (d>=1 || nd*d>=18.37) {
+ return 1;
+ }
+
+ if (nd<=0.5) {
+ return 0;
+ }
+ else if (nd<=1) {
+ double x = 2*d-1./n;
+ double s = x;
+ for (int i = 2; i<=n; ++i) {
+ s *= x*i;
+ }
+ return s;
+ }
+ else if (nd>=n-1) {
+ return 1-2*std::pow(1-d, n);
+ }
+
+ // If n>=10^5, or n>140 and n·d²>0.3, throw an
+ // exception rather than risk bad numeric behaviour.
+ // Fix this if required by implementing Pelz-Good.
+
+ if (n>=1e5 || (n>140 && nd*d>0.3)) {
+ throw std::range_error("approximation for parameters not implemented");
+ }
+
+ int k = std::ceil(nd);
+ int m = 2*k-1;
+ double h = k-nd;
+
+ // Representation of matrix H [Dubin eq 2.4.3, Carvalho eq 6].
+ std::vector<double> data(3*m-2);
+ double* v = data.data(); // length m
+ double* q = v+m; // length m
+ double* w = q+m; // length m-2
+
+ // Initialize H representation.
+ double r=1;
+ double p=h;
+ for (int i=0; i<m-1; ++i) {
+ if (i<m-2) {
+ w[i] = r;
+ }
+ v[i] = (1-p)*r;
+ p *= h;
+ r /= i+2;
+ }
+ double pp = h<=0.5? 0: std::pow(2*h-1, m);
+ v[m-1] = (1-2*p+pp)*r;
+ q[k-1] = 1;
+
+ // Compute n!/n^n H(k,k).
+ auto dot = [](double* a, double* b, int n) {
+ return std::inner_product(a, a+n, b, 0.);
+ };
+
+ constexpr double scale = 1e140;
+ constexpr double ooscale = 1/scale;
+ int scale_pow = 0;
+
+ double s = 1;
+ double oon = 1.0/n;
+ for (int i = 1; i<=n; ++i) {
+ s = i*oon;
+
+ double qprev = q[0];
+ q[0] = s*dot(v, q, m);
+ for (int j = 1; j<m-1; ++j) {
+ double a = qprev;
+ qprev = q[j];
+ q[j] = s*(dot(w, q+j, m-j-1) + v[m-j-1]*q[m-1] + a);
+ }
+ q[m-1] = s*(v[0]*q[m-1] + qprev);
+
+ if (q[k-1]>scale) {
+ for (int i = 0; i<m; ++i) q[i] *= ooscale;
+ ++scale_pow;
+ }
+ else if (q[k-1]<ooscale) {
+ for (int i = 0; i<m; ++i) q[i] *= scale;
+ --scale_pow;
+ }
+ }
+
+ return scale_pow? q[k-1]*std::pow(scale, scale_pow): q[k-1];
+}
+
+} // namespace ks
+
+
+namespace poisson {
+
+/* Approximate cdf of Poisson(μ) by using Wilson-Hilferty transform
+ * and then the normal distribution cdf.
+ *
+ * See e.g.:
+ *
+ * Terrell, G. 2003. "The Wilson–Hilferty transformation is locally saddlepoint,"
+ * Biometrika. 90 (2), pp. 445-453.
+ * DOI: 10.1093/biomet/90.2.445
+ */
+
+double poisson_cdf_approx(int n, double mu) {
+ double x = std::pow(0.5+n, 2./3.);
+ double norm_mu = std::pow(mu, 2./3.)*(1-1./(9*mu));
+ double norm_sd = 2./3.*std::pow(mu, 1./6.)*(1+1./(24*mu));
+
+ return 0.5+0.5*std::erf((x-norm_mu)/(norm_sd*std::sqrt(2.)));
+}
+
+
+} // namespace poisson
+
+} // namespace testing
diff --git a/tests/unit/stats.hpp b/tests/unit/stats.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b2430738055e16fa91cebed5cc1b203393a18458
--- /dev/null
+++ b/tests/unit/stats.hpp
@@ -0,0 +1,101 @@
+#pragma once
+
+#include <limits>
+
+#include <util/meta.hpp>
+
+/* For randomness, distribution testing we need some
+ * statistical tests... */
+
+namespace testing {
+
+// Simple summary statistics (min, max, mean, variance)
+// computed with on-line algorithm.
+//
+// Note that variance is population variance (i.e. with
+// denominator n, not n-1).
+
+struct summary_stats {
+ int n = 0;
+ double max = -std::numeric_limits<double>::infinity();
+ double min = std::numeric_limits<double>::infinity();
+ double mean = 0;
+ double variance = 0;
+};
+
+template <typename Seq>
+summary_stats summarize(const Seq& xs) {
+ summary_stats s;
+
+ for (auto x: xs) {
+ ++s.n;
+ s.min = x<s.min? x: s.min;
+ s.max = x>s.max? x: s.max;
+ double d1 = x-s.mean;
+ s.mean += d1/s.n;
+ double d2 = x-s.mean;
+ s.variance += d1*d2;
+ }
+
+ if (s.n) {
+ s.variance /= s.n;
+ }
+ return s;
+}
+
+// One-sample Kolmogorov-Smirnov test: test probability sample
+// comes from a continuous, 1-dimensional distribution.
+
+namespace ks {
+
+// Compute K-S test statistic:
+//
+// Input: n quantiles of the sample data F(x_1) ... F(x_n) in
+// increasing order, where F is the cdf of the specified distribution.
+//
+// Output: K-S test statistic Dn of the x_i.
+
+template <typename Seq>
+double dn_statistic(const Seq& qs) {
+ double n = static_cast<double>(nest::mc::util::size(qs));
+ double d = 0;
+ int j = 0;
+ for (auto q: qs) {
+ auto nq = n*q;
+ double d1 = nq-j;
+ ++j;
+ double d2 = j-nq;
+ d = std::max(d, std::max(d1, d2));
+ }
+ return d/n;
+}
+
+// One sample K-S test.
+//
+// Input: K-S test statistic d, number of observations n.
+// Output: P(Dn < d) where Dn ~ one-sample K–S statistic for n observations.
+//
+// Note: the implementation does not cover the full parameter space;
+// if N>140 and n·u² in [0.3, 18], the algorithm used may give a poor
+// result. The Pelz-Good approximation should be used here, but it has
+// not yet been implemented.
+
+double dn_cdf(double d, int n);
+
+} // namespace ks
+
+// Functions related to the Poisson distribution.
+
+namespace poisson {
+
+// Approximate cdf of Poisson distribution using Wilson-Hilferty transform.
+//
+// Input: integer sample n.
+// Output: P(X < n) where X ~ Poisson(mu).
+
+double poisson_cdf_approx(int n, double mu);
+
+} // namespace poisson
+
+} // namespace testing
+
diff --git a/tests/unit/test_any.cpp b/tests/unit/test_any.cpp
index 1bfb0627d0e8ce50696c812fb9cb5b63554cc6f6..594cb63a61404f43b23d2657498ceb550221eefc 100644
--- a/tests/unit/test_any.cpp
+++ b/tests/unit/test_any.cpp
@@ -1,9 +1,13 @@
+//#include <iostream>
+#include <type_traits>
+
#include "../gtest.h"
#include "common.hpp"
-#include <iostream>
-
#include <util/any.hpp>
+#include <util/any_ptr.hpp>
+
+#include <typeinfo>
using namespace nest::mc;
@@ -322,5 +326,120 @@ TEST(any, make_any) {
std::vector<int> ref{1, 2, 3};
EXPECT_EQ(ref, *vec);
}
+}
+
+// Tests below are for `any_ptr` class.
+
+TEST(any_ptr, ctor_and_assign) {
+ using util::any_ptr;
+
+ any_ptr p;
+
+ EXPECT_FALSE(p);
+ EXPECT_FALSE(p.has_value());
+
+ int x;
+ any_ptr q(&x);
+
+ EXPECT_TRUE(q);
+ EXPECT_TRUE(q.has_value());
+
+ p = q;
+
+ EXPECT_TRUE(p);
+ EXPECT_TRUE(p.has_value());
+
+ p = nullptr;
+
+ EXPECT_FALSE(p);
+ EXPECT_FALSE(p.has_value());
+
+ p = &x;
+
+ EXPECT_TRUE(p);
+ EXPECT_TRUE(p.has_value());
+
+ p.reset();
+
+ EXPECT_FALSE(p);
+ EXPECT_FALSE(p.has_value());
+
+ p.reset(&x);
+
+ EXPECT_TRUE(p);
+ EXPECT_TRUE(p.has_value());
+
+ p.reset(nullptr);
+
+ EXPECT_FALSE(p);
+ EXPECT_FALSE(p.has_value());
+
+ p = nullptr;
+ EXPECT_FALSE(p);
+ EXPECT_FALSE(p.has_value());
}
+
+
+TEST(any_ptr, ordering) {
+ using util::any_ptr;
+
+ int x[2];
+ double y;
+
+ any_ptr a(&x[0]);
+ any_ptr b(&x[1]);
+
+ EXPECT_LT(a, b);
+ EXPECT_LE(a, b);
+ EXPECT_NE(a, b);
+ EXPECT_GE(b, a);
+ EXPECT_GT(b, a);
+ EXPECT_FALSE(a==b);
+
+ any_ptr c(&y);
+
+ EXPECT_NE(c, a);
+ EXPECT_TRUE(a<c || a>c);
+ EXPECT_FALSE(a==c);
+}
+
+TEST(any_ptr, as_and_type) {
+ using util::any_ptr;
+
+ int x = 0;
+ const int y = 0;
+ any_ptr p;
+
+ EXPECT_FALSE(p.as<int*>());
+
+ p = &y;
+ EXPECT_FALSE(p.as<int*>());
+ EXPECT_TRUE(p.as<const int*>());
+ EXPECT_EQ(typeid(const int*), p.type());
+
+ p = &x;
+ EXPECT_TRUE(p.as<int*>());
+ EXPECT_FALSE(p.as<const int*>());
+ EXPECT_EQ(typeid(int*), p.type());
+
+ *p.as<int*>() = 3;
+ EXPECT_EQ(3, x);
+}
+
+TEST(any_ptr, any_cast) {
+ using util::any_ptr;
+ using util::any_cast;
+
+ int x = 0;
+ any_ptr p;
+
+ auto c1 = any_cast<int*>(p);
+ EXPECT_FALSE(c1);
+ EXPECT_TRUE((std::is_same<int*, decltype(c1)>::value));
+
+ p = &x;
+ auto c2 = any_cast<int*>(p);
+ EXPECT_TRUE(c2);
+}
+
diff --git a/tests/unit/test_domain_decomposition.cpp b/tests/unit/test_domain_decomposition.cpp
index 90c00f5e6e4757369117f38f4f6bbf03fd63f165..f9724af498840e37779176d170f136411a922f78 100644
--- a/tests/unit/test_domain_decomposition.cpp
+++ b/tests/unit/test_domain_decomposition.cpp
@@ -1,52 +1,28 @@
#include "../gtest.h"
+#include <stdexcept>
+
#include <backends.hpp>
#include <domain_decomposition.hpp>
#include <hardware/node_info.hpp>
#include <load_balance.hpp>
+#include "../simple_recipes.hpp"
+
using namespace nest::mc;
namespace {
- //
- // Dummy recipes type for testing.
- //
-
- // Homogenous cell population of cable cells.
- class homo_recipe: public recipe {
- public:
- homo_recipe(cell_size_type s): size_(s)
- {}
-
- cell_size_type num_cells() const override {
- return size_;
- }
-
- util::unique_any get_cell_description(cell_gid_type) const override {
- return {};
- }
- cell_kind get_cell_kind(cell_gid_type) const override {
- return cell_kind::cable1d_neuron;
- }
+ // Dummy recipes types for testing.
- cell_count_info get_cell_count_info(cell_gid_type) const override {
- return {0, 0, 0};
- }
- std::vector<cell_connection> connections_on(cell_gid_type) const override {
- return {};
- }
-
- private:
- cell_size_type size_;
- };
+ struct dummy_cell {};
+ using homo_recipe = homogeneous_recipe<cell_kind::cable1d_neuron, dummy_cell>;
// Heterogenous cell population of cable and rss cells.
- // Interleaved so that cells with even gid are cable cells, and even gid are
+ // Interleaved so that cells with even gid are cable cells, and odd gid are
// rss cells.
class hetero_recipe: public recipe {
public:
- hetero_recipe(cell_size_type s): size_(s)
- {}
+ hetero_recipe(cell_size_type s): size_(s) {}
cell_size_type num_cells() const override {
return size_;
@@ -55,19 +31,25 @@ namespace {
util::unique_any get_cell_description(cell_gid_type) const override {
return {};
}
+
cell_kind get_cell_kind(cell_gid_type gid) const override {
return gid%2?
cell_kind::regular_spike_source:
cell_kind::cable1d_neuron;
}
- cell_count_info get_cell_count_info(cell_gid_type) const override {
- return {0, 0, 0};
- }
+ cell_size_type num_sources(cell_gid_type) const override { return 0; }
+ cell_size_type num_targets(cell_gid_type) const override { return 0; }
+ cell_size_type num_probes(cell_gid_type) const override { return 0; }
+
std::vector<cell_connection> connections_on(cell_gid_type) const override {
return {};
}
+ probe_info get_probe(cell_member_type) const override {
+ throw std::logic_error("no probes");
+ }
+
private:
cell_size_type size_;
};
@@ -82,7 +64,7 @@ TEST(domain_decomposition, homogenous_population)
hw::node_info nd(1, 0);
unsigned num_cells = 10;
- const auto D = partition_load_balance(homo_recipe(num_cells), nd);
+ const auto D = partition_load_balance(homo_recipe(num_cells, dummy_cell{}), nd);
EXPECT_EQ(D.num_global_cells, num_cells);
EXPECT_EQ(D.num_local_cells, num_cells);
@@ -110,7 +92,7 @@ TEST(domain_decomposition, homogenous_population)
hw::node_info nd(1, 1);
unsigned num_cells = 10;
- const auto D = partition_load_balance(homo_recipe(num_cells), nd);
+ const auto D = partition_load_balance(homo_recipe(num_cells, dummy_cell{}), nd);
EXPECT_EQ(D.num_global_cells, num_cells);
EXPECT_EQ(D.num_local_cells, num_cells);
diff --git a/tests/unit/test_dss_cell_group.cpp b/tests/unit/test_dss_cell_group.cpp
index 5d853b10cb4f1e6d5b3ded429187805d46574321..8c841bc2f8ffe2eef3586a967c12b1358b8fed4c 100644
--- a/tests/unit/test_dss_cell_group.cpp
+++ b/tests/unit/test_dss_cell_group.cpp
@@ -1,75 +1,60 @@
#include "../gtest.h"
-#include "dss_cell_description.hpp"
-#include "dss_cell_group.hpp"
+#include <dss_cell_description.hpp>
+#include <dss_cell_group.hpp>
#include <util/unique_any.hpp>
+#include "../simple_recipes.hpp"
-using namespace nest::mc;
+using namespace nest::mc;
-TEST(dss_cell, constructor)
-{
- std::vector<time_type> spikes;
-
- std::vector<util::unique_any> cell_descriptions(1);
- cell_descriptions[0] = util::unique_any(dss_cell_description(spikes));
-
- dss_cell_group sut({0}, cell_descriptions);
-}
+using dss_recipe = homogeneous_recipe<cell_kind::data_spike_source, dss_cell_description>;
TEST(dss_cell, basic_usage)
{
- std::vector<time_type> spikes_to_emit;
-
- time_type spike_time = 0.1;
- spikes_to_emit.push_back(spike_time);
+ const time_type spike_time = 0.1;
+ dss_recipe rec(1u, dss_cell_description({spike_time}));
+ dss_cell_group sut({0}, rec);
- std::vector<util::unique_any> cell_descriptions(1);
- cell_descriptions[0] = util::unique_any(dss_cell_description(spikes_to_emit));
-
- dss_cell_group sut({0}, cell_descriptions);
-
- // no spikes in this time frame
- sut.advance(0.09, 0.01); // The dt (0,01) is not used
+ // No spikes in this time frame.
+ time_type dt = 0.01; // (note that dt is ignored in dss_cell_group).
+ sut.advance(0.09, dt);
auto spikes = sut.spikes();
- EXPECT_EQ(size_t(0), spikes.size());
+ EXPECT_EQ(0u, spikes.size());
- // only one in this time frame
+ // Only one in this time frame.
sut.advance(0.11, 0.01);
spikes = sut.spikes();
- EXPECT_EQ(size_t(1), spikes.size());
+ EXPECT_EQ(1u, spikes.size());
ASSERT_FLOAT_EQ(spike_time, spikes[0].time);
- // Clear the spikes after 'processing' them
+ // Clear the spikes after 'processing' them.
sut.clear_spikes();
spikes = sut.spikes();
- EXPECT_EQ(size_t(0), spikes.size());
+ EXPECT_EQ(0u, spikes.size());
- // No spike to be emitted
- sut.advance(0.12, 0.01);
+ // No spike to be emitted.
+ sut.advance(0.12, dt);
spikes = sut.spikes();
- EXPECT_EQ(size_t(0), spikes.size());
+ EXPECT_EQ(0u, spikes.size());
- // Reset the internal state to null
+ // Reset the internal state.
sut.reset();
- // Expect 10 excluding the 0.2
- sut.advance(0.2, 0.01);
+ // Expect to have the one spike again after reset.
+ sut.advance(0.2, dt);
spikes = sut.spikes();
- EXPECT_EQ(size_t(1), spikes.size());
+ EXPECT_EQ(1u, spikes.size());
ASSERT_FLOAT_EQ(spike_time, spikes[0].time);
}
TEST(dss_cell, cell_kind_correct)
{
- std::vector<time_type> spikes_to_emit;
-
- std::vector<util::unique_any> cell_descriptions(1);
- cell_descriptions[0] = util::unique_any(dss_cell_description(spikes_to_emit));
-
- dss_cell_group sut({0}, cell_descriptions);
+ const time_type spike_time = 0.1;
+ dss_recipe rec(1u, dss_cell_description({spike_time}));
+ dss_cell_group sut({0}, rec);
EXPECT_EQ(cell_kind::data_spike_source, sut.get_cell_kind());
}
diff --git a/tests/unit/test_fvm_multi.cpp b/tests/unit/test_fvm_multi.cpp
index 6a9863191e2f71d4bc1733caabb15fc4311cda55..78057b856024884a3576f70a2c369a4c9dc2a3c6 100644
--- a/tests/unit/test_fvm_multi.cpp
+++ b/tests/unit/test_fvm_multi.cpp
@@ -7,11 +7,14 @@
#include <cell.hpp>
#include <common_types.hpp>
#include <fvm_multicell.hpp>
+#include <recipe.hpp>
+#include <sampler_map.hpp>
#include <util/meta.hpp>
#include <util/rangeutil.hpp>
#include "../test_util.hpp"
-#include "../test_common_cells.hpp"
+#include "../common_cells.hpp"
+#include "../simple_recipes.hpp"
using fvm_cell =
nest::mc::fvm::fvm_multicell<nest::mc::multicore::backend>;
@@ -23,10 +26,10 @@ TEST(fvm_multi, cable)
nest::mc::cell cell=make_cell_ball_and_3stick();
std::vector<fvm_cell::target_handle> targets;
- std::vector<fvm_cell::probe_handle> probes;
+ probe_association_map<fvm_cell::probe_handle> probe_map;
fvm_cell fvcell;
- fvcell.initialize(util::singleton_view(cell), targets, probes);
+ fvcell.initialize({0}, cable1d_recipe(cell), targets, probe_map);
auto& J = fvcell.jacobian();
@@ -66,10 +69,10 @@ TEST(fvm_multi, init)
cell.segment(1)->set_compartments(10);
std::vector<fvm_cell::target_handle> targets;
- std::vector<fvm_cell::probe_handle> probes;
+ probe_association_map<fvm_cell::probe_handle> probe_map;
fvm_cell fvcell;
- fvcell.initialize(util::singleton_view(cell), targets, probes);
+ fvcell.initialize({0}, cable1d_recipe(cell), targets, probe_map);
// This is naughty: removing const from the matrix reference, but is needed
// to test the build_matrix() method below (which is only accessable
@@ -103,7 +106,6 @@ TEST(fvm_multi, init)
};
EXPECT_TRUE(is_neg(mat.u(1, J.size())));
EXPECT_TRUE(is_pos(mat.d));
-
}
TEST(fvm_multi, multi_init)
@@ -129,11 +131,13 @@ TEST(fvm_multi, multi_init)
cells[1].add_detector({0, 0}, 3.3);
- std::vector<fvm_cell::target_handle> targets(4);
- std::vector<fvm_cell::probe_handle> probes;
+ std::vector<fvm_cell::target_handle> targets;
+ probe_association_map<fvm_cell::probe_handle> probe_map;
fvm_cell fvcell;
- fvcell.initialize(cells, targets, probes);
+ fvcell.initialize({0, 1}, cable1d_recipe(cells), targets, probe_map);
+
+ EXPECT_EQ(4u, targets.size());
auto& J = fvcell.jacobian();
EXPECT_EQ(J.size(), 5u+13u);
@@ -169,7 +173,6 @@ TEST(fvm_multi, multi_init)
TEST(fvm_multi, stimulus)
{
using namespace nest::mc;
- using util::singleton_view;
// the default ball and stick has one stimulus at the terminal end of the dendrite
auto cell = make_cell_ball_and_stick();
@@ -191,10 +194,10 @@ TEST(fvm_multi, stimulus)
// as during the stimulus windows.
std::vector<fvm_cell::target_handle> targets;
- std::vector<fvm_cell::probe_handle> probes;
+ probe_association_map<fvm_cell::probe_handle> probe_map;
fvm_cell fvcell;
- fvcell.initialize(singleton_view(cell), targets, probes);
+ fvcell.initialize({0}, cable1d_recipe(cell), targets, probe_map);
auto ref = fvcell.find_mechanism("stimulus");
ASSERT_TRUE(ref) << "no stimuli retrieved from lowered fvm cell: expected 2";
@@ -283,10 +286,10 @@ TEST(fvm_multi, mechanism_indexes)
// generate the lowered fvm cell
std::vector<fvm_cell::target_handle> targets;
- std::vector<fvm_cell::probe_handle> probes;
+ probe_association_map<fvm_cell::probe_handle> probe_map;
fvm_cell fvcell;
- fvcell.initialize(util::singleton_view(c), targets, probes);
+ fvcell.initialize({0}, cable1d_recipe(c), targets, probe_map);
// make vectors with the expected CV indexes for each mechanism
std::vector<unsigned> hh_index = {0u, 4u, 5u, 6u, 7u, 8u};
@@ -380,11 +383,11 @@ void run_target_handle_test(std::vector<handle_info> all_handles) {
}
auto n = all_handles.size();
- std::vector<fvm_cell::target_handle> targets(n);
- std::vector<fvm_cell::probe_handle> probes;
+ std::vector<fvm_cell::target_handle> targets;
+ probe_association_map<fvm_cell::probe_handle> probe_map;
fvm_cell fvcell;
- fvcell.initialize(cells, targets, probes);
+ fvcell.initialize({0, 1}, cable1d_recipe(cells), targets, probe_map);
ASSERT_EQ(n, util::size(targets));
unsigned i = 0;
diff --git a/tests/unit/test_mc_cell_group.cpp b/tests/unit/test_mc_cell_group.cpp
index 4c157cb1e71578d7cc2f9982bde90c4ccd58f334..4d5b14692e3150d59c70fa40624a5aface3b91b4 100644
--- a/tests/unit/test_mc_cell_group.cpp
+++ b/tests/unit/test_mc_cell_group.cpp
@@ -7,7 +7,8 @@
#include <util/rangeutil.hpp>
#include "common.hpp"
-#include "../test_common_cells.hpp"
+#include "../common_cells.hpp"
+#include "../simple_recipes.hpp"
using namespace nest::mc;
using fvm_cell = fvm::fvm_multicell<nest::mc::multicore::backend>;
@@ -21,16 +22,15 @@ cell make_cell() {
return c;
}
-
TEST(mc_cell_group, get_kind) {
- mc_cell_group<fvm_cell> group{ {0}, util::singleton_view(make_cell()) };
+ mc_cell_group<fvm_cell> group{{0}, cable1d_recipe(make_cell()) };
// we are generating a mc_cell_group which should be of the correct type
EXPECT_EQ(cell_kind::cable1d_neuron, group.get_cell_kind());
}
TEST(mc_cell_group, test) {
- mc_cell_group<fvm_cell> group{{0}, util::singleton_view(make_cell())};
+ mc_cell_group<fvm_cell> group{{0}, cable1d_recipe(make_cell()) };
group.advance(50, 0.01);
@@ -40,28 +40,34 @@ TEST(mc_cell_group, test) {
}
TEST(mc_cell_group, sources) {
- using cell_group_type = mc_cell_group<fvm_cell>;
+ // Make twenty cells, with an extra detector on gids 0, 3 and 17
+ // to make things more interesting.
+ std::vector<cell> cells;
+
+ for (int i=0; i<20; ++i) {
+ cells.push_back(make_cell());
+ if (i==0 || i==3 || i==17) {
+ cells.back().add_detector({1, 0.3}, 2.3);
+ }
- auto cell = make_cell();
- EXPECT_EQ(cell.detectors().size(), 1u);
- // add another detector on the cell to make things more interesting
- cell.add_detector({1, 0.3}, 2.3);
+ EXPECT_EQ(1u + (i==0 || i==3 || i==17), cells.back().detectors().size());
+ }
- cell_gid_type first_gid = 37u;
- auto group = cell_group_type{{first_gid}, util::singleton_view(cell)};
+ std::vector<cell_gid_type> gids = {3u, 4u, 10u, 16u, 17u, 18u};
+ mc_cell_group<fvm_cell> group{gids, cable1d_recipe(cells)};
- // expect group sources to be lexicographically sorted by source id
- // with gids in cell group's range and indices starting from zero
+ // Expect group sources to be lexicographically sorted by source id
+ // with gids in cell group's range and indices starting from zero.
const auto& sources = group.spike_sources();
- for (unsigned i = 0; i<sources.size(); ++i) {
- auto id = sources[i];
- if (i==0) {
- EXPECT_EQ(id.gid, first_gid);
+ for (unsigned j = 0; j<sources.size(); ++j) {
+ auto id = sources[j];
+ if (j==0) {
+ EXPECT_EQ(id.gid, gids[0]);
EXPECT_EQ(id.index, 0u);
}
else {
- auto prev = sources[i-1];
+ auto prev = sources[j-1];
EXPECT_GT(id, prev);
EXPECT_EQ(id.index, id.gid==prev.gid? prev.index+1: 0u);
}
diff --git a/tests/unit/test_probe.cpp b/tests/unit/test_probe.cpp
index 8ea355eaf261e2610d8ba736244fe49232b29f1d..1caa2a53b0bbb6a52d9de8f7cb53e07a5ee04431 100644
--- a/tests/unit/test_probe.cpp
+++ b/tests/unit/test_probe.cpp
@@ -6,79 +6,60 @@
#include <fvm_multicell.hpp>
#include <util/rangeutil.hpp>
-TEST(probe, instantiation)
-{
- using namespace nest::mc;
-
- cell c1;
-
- segment_location loc1{0, 0};
- segment_location loc2{1, 0.6};
-
- auto p1 = c1.add_probe({loc1, probeKind::membrane_voltage});
- auto p2 = c1.add_probe({loc2, probeKind::membrane_current});
-
- // expect locally provided probe ids to be numbered sequentially from zero.
-
- EXPECT_EQ(0u, p1);
- EXPECT_EQ(1u, p2);
-
- // expect the probes() return to be a collection with these two probes.
-
- auto probes = c1.probes();
- EXPECT_EQ(2u, probes.size());
+#include "../common_cells.hpp"
+#include "../simple_recipes.hpp"
- EXPECT_EQ(loc1, probes[0].location);
- EXPECT_EQ(probeKind::membrane_voltage, probes[0].kind);
-
- EXPECT_EQ(loc2, probes[1].location);
- EXPECT_EQ(probeKind::membrane_current, probes[1].kind);
-}
+using namespace nest::mc;
TEST(probe, fvm_multicell)
{
- using namespace nest::mc;
+ using fvm_cell = fvm::fvm_multicell<nest::mc::multicore::backend>;
- cell bs;
+ cell bs = make_cell_ball_and_stick(false);
- // ball-and-stick model morphology
+ i_clamp stim(0, 100, 0.3);
+ bs.add_stimulus({1, 1}, stim);
- bs.add_soma(12.6157/2.0);
- bs.add_cable(0, section_kind::dendrite, 0.5, 0.5, 200);
- bs.soma()->set_compartments(5);
+ cable1d_recipe rec(bs);
segment_location loc0{0, 0};
segment_location loc1{1, 1};
- segment_location loc2{1, 0.5};
+ segment_location loc2{1, 0.7};
- bs.add_probe({loc0, probeKind::membrane_voltage});
- bs.add_probe({loc1, probeKind::membrane_voltage});
- bs.add_probe({loc2, probeKind::membrane_current});
+ rec.add_probe(0, 10, cell_probe_address{loc0, cell_probe_address::membrane_voltage});
+ rec.add_probe(0, 20, cell_probe_address{loc1, cell_probe_address::membrane_voltage});
+ rec.add_probe(0, 30, cell_probe_address{loc2, cell_probe_address::membrane_current});
- i_clamp stim(0, 100, 0.3);
- bs.add_stimulus({1, 1}, stim);
+ std::vector<fvm_cell::target_handle> targets;
+ probe_association_map<fvm_cell::probe_handle> probe_map;
+
+ fvm_cell lcell;
+ lcell.initialize({0}, rec, targets, probe_map);
+
+ EXPECT_EQ(3u, rec.num_probes(0));
+ EXPECT_EQ(3u, probe_map.size());
- using fvm_multicell = fvm::fvm_multicell<nest::mc::multicore::backend>;
- std::vector<fvm_multicell::target_handle> targets;
- std::vector<fvm_multicell::probe_handle> probes{3};
+ EXPECT_EQ(10, probe_map.at({0, 0}).tag);
+ EXPECT_EQ(20, probe_map.at({0, 1}).tag);
+ EXPECT_EQ(30, probe_map.at({0, 2}).tag);
- fvm_multicell lcell;
- lcell.initialize(util::singleton_view(bs), targets, probes);
+ fvm_cell::probe_handle p0 = probe_map.at({0, 0}).handle;
+ fvm_cell::probe_handle p1 = probe_map.at({0, 1}).handle;
+ fvm_cell::probe_handle p2 = probe_map.at({0, 2}).handle;
// Know from implementation that probe_handle.second
// is a compartment index: expect probe values and
- // direct queries of voltage and current
- // to be equal in fvm cell
+ // direct queries of voltage and current to be equal in fvm_cell.
- EXPECT_EQ(lcell.voltage()[probes[0].second], lcell.probe(probes[0]));
- EXPECT_EQ(lcell.voltage()[probes[1].second], lcell.probe(probes[1]));
- EXPECT_EQ(lcell.current()[probes[2].second], lcell.probe(probes[2]));
+ EXPECT_EQ(lcell.voltage()[p0.second], lcell.probe(p0));
+ EXPECT_EQ(lcell.voltage()[p1.second], lcell.probe(p1));
+ EXPECT_EQ(lcell.current()[p2.second], lcell.probe(p2));
lcell.setup_integration(0.05, 0.05);
lcell.step_integration();
- EXPECT_EQ(lcell.voltage()[probes[0].second], lcell.probe(probes[0]));
- EXPECT_EQ(lcell.voltage()[probes[1].second], lcell.probe(probes[1]));
- EXPECT_EQ(lcell.current()[probes[2].second], lcell.probe(probes[2]));
+ EXPECT_EQ(lcell.voltage()[p0.second], lcell.probe(p0));
+ EXPECT_EQ(lcell.voltage()[p1.second], lcell.probe(p1));
+ EXPECT_EQ(lcell.current()[p2.second], lcell.probe(p2));
}
diff --git a/tests/unit/test_range.cpp b/tests/unit/test_range.cpp
index e058604f7b31d1341328f42ee06efba97eaf550e..52e054eed48fd6de96920f19889ee0254efecc7a 100644
--- a/tests/unit/test_range.cpp
+++ b/tests/unit/test_range.cpp
@@ -22,6 +22,7 @@
using namespace nest::mc;
using testing::null_terminated;
+using testing::nocopy;
TEST(range, list_iterator) {
std::list<int> l = { 2, 4, 6, 8, 10 };
@@ -477,6 +478,43 @@ TEST(range, all_of_any_of) {
EXPECT_TRUE(util::any_of(cstr("87654x"), pred));
}
+TEST(range, keys) {
+ {
+ std::map<int, double> map = {{10, 2.0}, {3, 8.0}};
+ std::vector<int> expected = {3, 10};
+ std::vector<int> keys = util::assign_from(util::keys(map));
+ EXPECT_EQ(expected, keys);
+ }
+
+ {
+ struct cmp {
+ bool operator()(const nocopy<int>& a, const nocopy<int>& b) const {
+ return a.value<b.value;
+ }
+ };
+ std::map<nocopy<int>, double, cmp> map;
+ map.insert(std::pair<nocopy<int>, double>(11, 2.0));
+ map.insert(std::pair<nocopy<int>, double>(2, 0.3));
+ map.insert(std::pair<nocopy<int>, double>(2, 0.8));
+ map.insert(std::pair<nocopy<int>, double>(5, 0.1));
+
+ std::vector<int> expected = {2, 5, 11};
+ std::vector<int> keys;
+ for (auto& k: util::keys(map)) {
+ keys.push_back(k.value);
+ }
+ EXPECT_EQ(expected, keys);
+ }
+
+ {
+ std::unordered_multimap<int, double> map = {{3, 0.1}, {5, 0.4}, {11, 0.8}, {5, 0.2}};
+ std::vector<int> expected = {3, 5, 5, 11};
+ std::vector<int> keys = util::assign_from(util::keys(map));
+ util::sort(keys);
+ EXPECT_EQ(expected, keys);
+ }
+}
+
#ifdef NMC_HAVE_TBB
TEST(range, tbb_split) {
diff --git a/tests/unit/test_rss_cell.cpp b/tests/unit/test_rss_cell.cpp
index e12bbb472743db453a4097bfae2bcf5293326996..5812c2c2e35b888f54d5fa79d823bf0647d8c7f9 100644
--- a/tests/unit/test_rss_cell.cpp
+++ b/tests/unit/test_rss_cell.cpp
@@ -1,57 +1,68 @@
#include "../gtest.h"
-#include "rss_cell.hpp"
+#include <rss_cell.hpp>
+#include <rss_cell_group.hpp>
-using namespace nest::mc;
+#include "../simple_recipes.hpp"
-TEST(rss_cell, constructor)
-{
- rss_cell test(0.0, 0.01, 1.0);
-}
+using namespace nest::mc;
+
+using rss_recipe = homogeneous_recipe<cell_kind::regular_spike_source, rss_cell>;
TEST(rss_cell, basic_usage)
{
- rss_cell sut(0.1, 0.01, 0.2);
+ constexpr time_type dt = 0.01; // dt is ignored by rss_cell_group::advance().
- // no spikes in this time frame
- auto spikes = sut.spikes_until(0.09);
- EXPECT_EQ(size_t(0), spikes.size());
+ // Use floating point times with an exact representation in order to avoid
+ // rounding issues.
+ rss_cell desc{0.125, 0.03125, 0.5};
+ rss_cell_group sut({0}, rss_recipe(1u, desc));
- //only on in this time frame
- spikes = sut.spikes_until(0.11);
- EXPECT_EQ(size_t(1), spikes.size());
+ // No spikes in this time frame.
+ sut.advance(0.1, dt);
+ EXPECT_EQ(0u, sut.spikes().size());
- // Reset the internal state to null
+ // Only on in this time frame
+ sut.clear_spikes();
+ sut.advance(0.127, dt);
+ EXPECT_EQ(1u, sut.spikes().size());
+
+ // Reset cell group state.
sut.reset();
- // Expect 10 excluding the 0.2
- spikes = sut.spikes_until(0.2);
- EXPECT_EQ(size_t(10), spikes.size());
+ // Expect 12 spikes excluding the 0.5 end point.
+ sut.advance(0.5, dt);
+ EXPECT_EQ(12u, sut.spikes().size());
}
TEST(rss_cell, poll_time_after_end_time)
{
- rss_cell sut(0.1, 0.01, 0.2);
+ constexpr time_type dt = 0.01; // dt is ignored by rss_cell_group::advance().
+
+ rss_cell desc{0.125, 0.03125, 0.5};
+ rss_cell_group sut({0}, rss_recipe(1u, desc));
- // no spikes in this time frame
- auto spikes = sut.spikes_until(0.3);
- EXPECT_EQ(size_t(10), spikes.size());
+ // Expect 12 spikes in this time frame.
+ sut.advance(0.7, dt);
+ EXPECT_EQ(12u, sut.spikes().size());
- // now ask for spikes for a time slot already passed.
- spikes = sut.spikes_until(0.2);
+ // Now ask for spikes for a time slot already passed:
// It should result in zero spikes because of the internal state!
- EXPECT_EQ(size_t(0), spikes.size());
+ sut.clear_spikes();
+ sut.advance(0.2, dt);
+ EXPECT_EQ(0u, sut.spikes().size());
sut.reset();
- // Expect 10 excluding the 0.2
- spikes = sut.spikes_until(0.2);
- EXPECT_EQ(size_t(10), spikes.size());
+ // Expect 12 excluding the 0.5
+ sut.advance(0.5, dt);
+ EXPECT_EQ(12u, sut.spikes().size());
}
TEST(rss_cell, cell_kind_correct)
{
- rss_cell sut(0.1, 0.01, 0.2);
+ rss_cell desc{0.1, 0.01, 0.2};
+ rss_cell_group sut({0}, rss_recipe(1u, desc));
EXPECT_EQ(cell_kind::regular_spike_source, sut.get_cell_kind());
}
diff --git a/tests/unit/test_schedule.cpp b/tests/unit/test_schedule.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..055b17eebfeea64c46497481cfd36f2a7f6a3642
--- /dev/null
+++ b/tests/unit/test_schedule.cpp
@@ -0,0 +1,305 @@
+#include <algorithm>
+#include <random>
+#include <stdexcept>
+#include <vector>
+
+#include <common_types.hpp>
+#include <schedule.hpp>
+#include <util/partition.hpp>
+#include <util/rangeutil.hpp>
+
+#include "common.hpp"
+#include "stats.hpp"
+
+using namespace nest::mc;
+using namespace testing;
+
+// Pull events from n non-contiguous subintervals of [t0, t1)
+// and check for monotonicity and boundedness.
+
+void run_invariant_checks(schedule S, time_type t0, time_type t1, unsigned n, int seed=0) {
+ if (!n) return;
+
+ std::minstd_rand R(seed);
+ std::uniform_real_distribution<time_type> U(t0, t1);
+
+ std::vector<time_type> divisions = {t0, t1};
+ std::generate_n(std::back_inserter(divisions), 2*(n-1), [&] { return U(R); });
+ util::sort(divisions);
+
+ bool skip = false;
+ for (auto ival: util::partition_view(divisions)) {
+ if (!skip) {
+ auto ts = S.events(ival.first, ival.second);
+
+ EXPECT_TRUE(std::is_sorted(ts.begin(), ts.end()));
+ if (!ts.empty()) {
+ EXPECT_LE(t0, ts.front());
+ EXPECT_GT(t1, ts.back());
+ }
+ }
+ skip = !skip;
+ }
+}
+
+// Take events from n contiguous intervals comprising [t0, t1), reset, and
+// then compare with events taken from a different set of contiguous
+// intervals comprising [t0, t1).
+
+void run_reset_check(schedule S, time_type t0, time_type t1, unsigned n, int seed=0) {
+ if (!n) return;
+
+ std::minstd_rand R(seed);
+ std::uniform_real_distribution<time_type> U(t0, t1);
+
+ std::vector<time_type> first_div = {t0, t1};
+ std::generate_n(std::back_inserter(first_div), n-1, [&] { return U(R); });
+ util::sort(first_div);
+
+ std::vector<time_type> second_div = {t0, t1};
+ std::generate_n(std::back_inserter(second_div), n-1, [&] { return U(R); });
+ util::sort(second_div);
+
+ std::vector<time_type> first;
+ for (auto ival: util::partition_view(first_div)) {
+ util::append(first, S.events(ival.first, ival.second));
+ }
+
+ S.reset();
+ std::vector<time_type> second;
+ for (auto ival: util::partition_view(second_div)) {
+ util::append(second, S.events(ival.first, ival.second));
+ }
+
+ EXPECT_EQ(first, 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));
+
+ S.reset();
+ EXPECT_EQ(expected, 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));
+}
+
+TEST(schedule, regular_invariants) {
+ SCOPED_TRACE("regular_invariants");
+ run_invariant_checks(regular_schedule(0.3), 3, 12, 7);
+}
+
+TEST(schedule, regular_reset) {
+ SCOPED_TRACE("regular_reset");
+ run_reset_check(regular_schedule(0.3), 3, 12, 7);
+}
+
+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));
+
+ S.reset();
+ EXPECT_EQ(expected, 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));
+}
+
+TEST(schedule, explicit_invariants) {
+ SCOPED_TRACE("explicit_invariants");
+
+ time_type times[] = {0.1, 0.3, 0.4, 0.42, 2.1, 2.3, 6.01, 9, 9.1, 9.8, 10, 11.2, 13};
+ run_invariant_checks(explicit_schedule(times), 0.4, 10.2, 5);
+}
+
+TEST(schedule, explicit_reset) {
+ SCOPED_TRACE("explicit_reset");
+
+ time_type times[] = {0.1, 0.3, 0.4, 0.42, 2.1, 2.3, 6.01, 9, 9.1, 9.8, 10, 11.2, 13};
+ run_reset_check(explicit_schedule(times), 0.4, 10.2, 5);
+}
+
+// A Uniform Random Bit Generator[*] adaptor that deliberately
+// skews the generated numbers by raising their quantile to
+// the given power.
+//
+// [*] Not actually uniform.
+
+template <typename RNG>
+struct skew_adaptor {
+ using result_type = typename RNG::result_type;
+ static constexpr result_type min() { return RNG::min(); }
+ static constexpr result_type max() { return RNG::max(); }
+
+ explicit skew_adaptor(double power): power_(power) {}
+ result_type operator()() {
+ constexpr double scale = max()-min();
+ constexpr double ooscale = 1./scale;
+
+ double x = ooscale*(G_()-min());
+ x = std::pow(x, power_);
+ return min()+scale*x;
+ }
+
+private:
+ RNG G_;
+ double power_;
+};
+
+template <typename RNG>
+double poisson_schedule_dispersion(int nbin, double mean_dt, RNG& G) {
+ schedule S = poisson_schedule(mean_dt, G);
+
+ std::vector<int> bin(nbin);
+ for (auto t: S.events(0, nbin)) {
+ int j = (int)t;
+ if (j<0 || j>=nbin) {
+ throw std::logic_error("poisson schedule result out of bounds");
+ }
+ ++bin[j];
+ }
+
+ summary_stats stats = summarize(bin);
+ return stats.mean/stats.variance;
+}
+
+// NOTE: schedule.poisson_uniformity tests can be expected to
+// fail approximately 1% of the time, if the underlying
+// random sequence were allowed to vary freely.
+
+TEST(schedule, poisson_uniformity) {
+ // Run Poisson dispersion test for N=101 with two-sided
+ // χ²-test critical value α=0.01.
+ //
+ // Test based on: N·dispersion ~ χ²(N-1) (approximately)
+ //
+ // F(chi2_lb; N-1) = α/2
+ // F(chi2_ub; N-1) = 1-α/2
+ //
+ // Numbers taken from scipy:
+ // scipy.stats.chi2.isf(0.01/2, 1000)
+ // scipy.stats.chi2.isf(1-0.01/2, 1000)
+
+ constexpr int N = 1001;
+ //constexpr double alpha = 0.01;
+ constexpr double chi2_lb = 888.56352318146696;
+ constexpr double chi2_ub = 1118.9480663231843;
+
+ std::mt19937_64 G;
+ double dispersion = poisson_schedule_dispersion(N, 1.23, G);
+ double test_value = N*dispersion;
+ EXPECT_GT(test_value, chi2_lb);
+ EXPECT_LT(test_value, chi2_ub);
+
+ // 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);
+ int n = (int)events.size();
+ double dn = ks::dn_statistic(events);
+
+ EXPECT_LT(ks::dn_cdf(dn, n), 0.99);
+
+ // Check that these tests fail for a non-Poisson
+ // source.
+
+ skew_adaptor<std::mt19937_64> W(1.5);
+ dispersion = poisson_schedule_dispersion(N, 1.23, 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);
+ n = (int)events.size();
+ dn = ks::dn_statistic(events);
+
+ // This test is currently failing, because we can't
+ // use a sufficiently high `n` in the `dn_cdf` function
+ // to get enough discrimination from the K-S test at
+ // 1%. TODO: Fix this by implementing n>140 case in
+ // `dn_cdf`.
+
+ // EXPECT_GT(ks::dn_cdf(dn, n), 0.99);
+}
+
+TEST(schedule, poisson_rate) {
+ // Test Poisson events over an interval against
+ // corresponding Poisson distribution.
+
+ constexpr double alpha = 0.01;
+ constexpr double lambda = 123.4;
+
+ std::mt19937_64 G;
+ schedule S = poisson_schedule(1./lambda, G);
+ int n = (int)S.events(0,1).size();
+ double cdf = poisson::poisson_cdf_approx(n, lambda);
+
+ EXPECT_GT(cdf, alpha/2);
+ EXPECT_LT(cdf, 1-alpha/2);
+
+ // Check that the test fails for a non-Poisson
+ // source.
+
+ skew_adaptor<std::mt19937_64> W(1.5);
+ S = poisson_schedule(1./lambda, W);
+ n = (int)S.events(0,1).size();
+ cdf = poisson::poisson_cdf_approx(n, lambda);
+
+ EXPECT_FALSE(cdf>=alpha/2 && cdf<=1-alpha/2);
+}
+
+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);
+}
+
+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);
+}
+
+TEST(schedule, poisson_offset) {
+ // Expect Poisson schedule with an offset to give exactly the
+ // same sequence, after the offset, as a regular zero-based Poisson.
+
+ const double offset = 3.3;
+
+ std::mt19937_64 G1;
+ G1.discard(300);
+
+ std::vector<time_type> expected;
+ for (auto t: poisson_schedule(3.4, G1).events(0., 100.)) {
+ t += offset;
+ if (t<100.) {
+ expected.push_back(t);
+ }
+ }
+
+ std::mt19937_64 G2;
+ G2.discard(300);
+
+ EXPECT_TRUE(seq_almost_eq<time_type>(expected, poisson_schedule(offset, 3.4, 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);
+}
+
diff --git a/tests/unit/test_stats.cpp b/tests/unit/test_stats.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e9c84ade8005ea2762cd1d2e0ac11ce58e830de1
--- /dev/null
+++ b/tests/unit/test_stats.cpp
@@ -0,0 +1,157 @@
+#include <string>
+#include <string>
+
+
+#include "common.hpp"
+#include "stats.hpp"
+
+/* Unit tests for the test statistics routines */
+
+using namespace testing;
+
+TEST(stats, dn_statistic) {
+ // Dn should equal the maximum vertical distance between the
+ // empirical distribution and the straight line from (0,0) to
+ // (1,1).
+
+ // Expected: statistic=0.75 (empirical cdf should be 1 at 0.25).
+ double x1[] = {0.25};
+ EXPECT_DOUBLE_EQ(0.75, ks::dn_statistic(x1));
+
+ // Expected: statistic=0.55 (empirical cdf should be 0.75 at 0.2).
+ double x2[] = {0.1, 0.15, 0.2, 0.8};
+ EXPECT_DOUBLE_EQ(0.55, ks::dn_statistic(x2));
+
+ // Expected: statistic=1 (empirical cdf should be 0 at 1¯).
+ double x3[] = {1.0};
+ EXPECT_DOUBLE_EQ(1.0, ks::dn_statistic(x3));
+}
+
+TEST(stats, dn_cdf) {
+ using std::pow;
+ using std::tgamma;
+
+ // Compute cdf F(x; n) of Dn for some known values and accuracies.
+ // (See e.g. Simard and L'Ecuyer 2011.).
+
+ // Zero and one tails:
+ // 0 if x ≤ 1/2n; 1 if x ≥ 1.
+
+ EXPECT_EQ(0.0, ks::dn_cdf(0.0, 1));
+ EXPECT_EQ(0.0, ks::dn_cdf(0.0, 10));
+ EXPECT_EQ(0.0, ks::dn_cdf(0.01, 10));
+ EXPECT_EQ(0.0, ks::dn_cdf(0.04999, 10));
+ EXPECT_EQ(0.0, ks::dn_cdf(0.00004999, 10000));
+ EXPECT_EQ(1.0, ks::dn_cdf(1, 1));
+ EXPECT_EQ(1.0, ks::dn_cdf(1234.45, 1));
+ EXPECT_EQ(1.0, ks::dn_cdf(1, 10000));
+ EXPECT_EQ(1.0, ks::dn_cdf(1234.45, 10000));
+
+ // When x in [1/2n, 1/n), F(x; n) = n!(2x-1/n)^n.
+ int n = 3;
+ double x = 0.3;
+ double expected = tgamma(n+1)*pow(2*x-1./n, n);
+ EXPECT_NEAR(expected, ks::dn_cdf(x, n), expected*1e-15);
+
+ // When x in [1-1/n, 1), F(x; n) = 1-2(1-x)^n.
+ n = 5;
+ x = 0.81;
+ expected = 1-2*pow(1-x, n);
+ EXPECT_NEAR(expected, ks::dn_cdf(x, n), expected*1e-15);
+
+ // When n·x^2 > 18.37, F(x; n) should be within double epsilon of 1.
+ n = 75;
+ x = 0.5;
+ EXPECT_EQ(1., ks::dn_cdf(x, n));
+
+ // Various spots in the middle (avoiding n>140 until we have
+ // a more complete implementation).
+
+ n = 100;
+ x = 0.2;
+ expected = 1-0.000555192732802810;
+ EXPECT_NEAR(expected, ks::dn_cdf(x, n), 1e-15); // note: absolute error bound
+
+ n = 140;
+ x = 0.0464158883361278;
+ expected = 0.0902623294750042;
+ EXPECT_NEAR(expected, ks::dn_cdf(x, n), expected*1e-14); // note: larger rel tol here
+}
+
+TEST(stats, running) {
+ // Exercise simple summary statistics:
+ summary_stats S;
+
+ double x1[] = {};
+ S = summarize(x1);
+
+ EXPECT_EQ(0, S.n);
+ EXPECT_EQ(0, S.mean);
+
+ double x2[] = {3.1, 3.1, 3.1};
+ S = summarize(x2);
+
+ EXPECT_EQ(3, S.n);
+ EXPECT_EQ(3.1, S.mean);
+ EXPECT_EQ(0, S.variance);
+ EXPECT_EQ(3.1, S.min);
+ EXPECT_EQ(3.1, S.max);
+
+ constexpr double inf = std::numeric_limits<double>::infinity();
+ double x3[] = {0.3, inf, 1.};
+ S = summarize(x3);
+
+ EXPECT_EQ(3, S.n);
+ EXPECT_EQ(0.3, S.min);
+ EXPECT_EQ(inf, S.max);
+
+ double x4[] = {-1, -2, 1, 0, 2};
+ S = summarize(x4);
+
+ EXPECT_EQ(5, S.n);
+ EXPECT_DOUBLE_EQ(0, S.mean);
+ EXPECT_DOUBLE_EQ(2, S.variance);
+ EXPECT_EQ(-2, S.min);
+ EXPECT_EQ(2, S.max);
+}
+
+TEST(stats, poisson_cdf) {
+ // Wilson-Hilferty transform approximation of the Poisson CDF
+ // is expected to be within circa 1e-3 (absolute) of true for
+ // parameter μ ≥ 6.
+
+ // Exact values in table computed with scipy.stats.
+ struct point {
+ int n;
+ double mu;
+ double cdf;
+ };
+
+ point points[] = {
+ {1, 6.00, 0.017351265236665},
+ {3, 6.00, 0.151203882776648},
+ {6, 6.00, 0.606302782412592},
+ {10, 6.00, 0.957379076417462},
+ {4, 11.30, 0.012323331570747},
+ {7, 11.30, 0.124852978325848},
+ {11, 11.30, 0.543501467974833},
+ {18, 11.30, 0.977530236363240},
+ {13, 23.00, 0.017428210696087},
+ {18, 23.00, 0.174768719569196},
+ {23, 23.00, 0.555149935616771},
+ {32, 23.00, 0.971056607478885},
+ {31, 44.70, 0.019726646735143},
+ {38, 44.70, 0.177614024770857},
+ {44, 44.70, 0.498035947942948},
+ {58, 44.70, 0.976892082270190},
+ {171, 200.00, 0.019982398243966},
+ {185, 200.00, 0.152411852483996},
+ {200, 200.00, 0.518794309678716},
+ {228, 200.00, 0.976235501016406},
+ };
+
+ for (auto p: points) {
+ SCOPED_TRACE("n="+std::to_string(p.n)+"; mu="+std::to_string(p.mu));
+ EXPECT_NEAR(p.cdf, poisson::poisson_cdf_approx(p.n, p.mu), 1.e-3);
+ }
+}
diff --git a/tests/validation/convergence_test.hpp b/tests/validation/convergence_test.hpp
index 7c9f593634fc80334f4b1cb3afbdaa6673cc82e6..daf6f4b885c7f09f872e648be4a709812bfaa286 100644
--- a/tests/validation/convergence_test.hpp
+++ b/tests/validation/convergence_test.hpp
@@ -1,8 +1,12 @@
#pragma once
+#include <vector>
+
+#include <model.hpp>
+#include <sampling.hpp>
+#include <simple_sampler.hpp>
#include <util/filter.hpp>
#include <util/rangeutil.hpp>
-#include <cell.hpp>
#include <json/json.hpp>
@@ -14,10 +18,9 @@
namespace nest {
namespace mc {
-struct sampler_info {
+struct probe_label {
const char* label;
- cell_member_type probe;
- simple_sampler sampler;
+ cell_member_type probe_id;
};
/* Common functionality for testing convergence towards
@@ -34,25 +37,25 @@ private:
std::string param_name_;
bool run_validation_;
nlohmann::json meta_;
- std::vector<sampler_info> cell_samplers_;
- std::map<std::string, trace_data> ref_data_;
+ std::vector<probe_label> probe_labels_;
+ std::map<std::string, trace_data<double>> ref_data_;
std::vector<conv_entry<Param>> conv_tbl_;
public:
- template <typename SamplerInfoSeq>
+ template <typename ProbeLabelSeq>
convergence_test_runner(
const std::string& param_name,
- const SamplerInfoSeq& samplers,
+ const ProbeLabelSeq& probe_labels,
const nlohmann::json& meta
):
param_name_(param_name),
run_validation_(false),
meta_(meta)
{
- util::assign(cell_samplers_, samplers);
+ util::assign(probe_labels_, probe_labels);
}
- // allow free access to JSON meta data attached to saved traces
+ // Allow free access to JSON meta data attached to saved traces.
nlohmann::json& metadata() { return meta_; }
void load_reference_data(const util::path& ref_path) {
@@ -60,8 +63,8 @@ public:
try {
ref_data_ = g_trace_io.load_traces(ref_path);
- run_validation_ = util::all_of(cell_samplers_,
- [&](const sampler_info& se) { return ref_data_.count(se.label)>0; });
+ run_validation_ = util::all_of(probe_labels_,
+ [&](const probe_label& pl) { return ref_data_.count(pl.label)>0; });
EXPECT_TRUE(run_validation_);
}
@@ -70,34 +73,50 @@ public:
}
}
- template <typename Model>
- void run(Model& m, Param p, float t_end, float dt, const std::vector<float>& excl) {
- // reset samplers and attach to probe locations
- for (auto& se: cell_samplers_) {
- se.sampler.reset();
- m.attach_sampler(se.probe, se.sampler.sampler());
+ void run(model& m, Param p, float sample_dt, float t_end, float dt, const std::vector<float>& excl) {
+ struct sampler_state {
+ sampler_association_handle h; // Keep these for clean up at end.
+ const char* label;
+ trace_data<double> trace;
+ };
+
+ // Attach new samplers to labelled probe ids.
+ std::vector<sampler_state> samplers;
+ samplers.reserve(probe_labels_.size());
+
+ for (auto& pl: probe_labels_) {
+ samplers.push_back({});
+ auto& entry = samplers.back();
+
+ entry.label = pl.label;
+ entry.h = m.add_sampler(one_probe(pl.probe_id), regular_schedule(sample_dt), simple_sampler<double>(entry.trace));
}
m.run(t_end, dt);
- for (auto& se: cell_samplers_) {
- std::string label = se.label;
- const auto& trace = se.sampler.trace;
+ for (auto& entry: samplers) {
+ std::string label = entry.label;
+ const auto& trace = entry.trace;
- // save trace
+ // Save trace.
nlohmann::json trace_meta(meta_);
trace_meta[param_name_] = p;
g_trace_io.save_trace(label, trace, trace_meta);
- // compute metrics
+ // Compute metreics.
if (run_validation_) {
- double linf = linf_distance(trace, ref_data_[label], excl);
+ double linf = linf_distance(ref_data_[label], trace, excl);
auto pd = peak_delta(trace, ref_data_[label]);
conv_tbl_.push_back({label, p, linf, pd});
}
}
+
+ // Remove added samplers.
+ for (const auto& entry: samplers) {
+ m.remove_sampler(entry.h);
+ }
}
void report() {
@@ -109,10 +128,10 @@ public:
}
void assert_all_convergence() const {
- for (const sampler_info& se: cell_samplers_) {
- SCOPED_TRACE(se.label);
+ for (const auto& pl: probe_labels_) {
+ SCOPED_TRACE(pl.label);
assert_convergence(util::filter(conv_tbl_,
- [&](const conv_entry<Param>& e) { return e.id==se.label; }));
+ [&](const conv_entry<Param>& e) { return e.id==pl.label; }));
}
}
};
diff --git a/tests/validation/trace_analysis.cpp b/tests/validation/trace_analysis.cpp
index b83d41a2ca9ce76d79d63dd84744294366a56c18..befe39f57ed43fad6965316c73bd8f59df704d72 100644
--- a/tests/validation/trace_analysis.cpp
+++ b/tests/validation/trace_analysis.cpp
@@ -18,7 +18,7 @@ namespace nest {
namespace mc {
struct trace_interpolant {
- trace_interpolant(const trace_data& trace): trace_(trace) {}
+ trace_interpolant(const trace_data<double>& trace): trace_(trace) {}
double operator()(float t) const {
if (trace_.empty()) return std::nan("");
@@ -37,15 +37,15 @@ struct trace_interpolant {
return math::lerp(vx[i], vx[i+1], (t-p.first)/(p.second-p.first));
}
- const trace_data& trace_;
+ const trace_data<double>& trace_;
};
-double linf_distance(const trace_data& u, const trace_data& ref) {
- trace_interpolant f{ref};
+double linf_distance(const trace_data<double>& u, const trace_data<double>& r) {
+ trace_interpolant f{r};
return util::max_value(
util::transform_view(u,
- [&](trace_entry x) { return std::abs(x.v-f(x.t)); }));
+ [&](trace_entry<double> x) { return std::abs(x.v-f(x.t)); }));
}
// Compute linf distance as above, but excluding sample points that lie
@@ -53,10 +53,10 @@ double linf_distance(const trace_data& u, const trace_data& ref) {
//
// `excl` contains the times to exclude, in ascending order.
-double linf_distance(const trace_data& u, const trace_data& ref, const std::vector<float>& excl) {
+double linf_distance(const trace_data<double>& u, const trace_data<double>& ref, const std::vector<float>& excl) {
trace_interpolant f{ref};
- trace_data reduced;
+ trace_data<double> reduced;
unsigned nexcl = excl.size();
unsigned ei = 0;
@@ -85,7 +85,7 @@ double linf_distance(const trace_data& u, const trace_data& ref, const std::vect
return linf_distance(reduced, ref);
}
-std::vector<trace_peak> local_maxima(const trace_data& u) {
+std::vector<trace_peak> local_maxima(const trace_data<double>& u) {
std::vector<trace_peak> peaks;
if (u.size()<2) return peaks;
@@ -116,7 +116,7 @@ std::vector<trace_peak> local_maxima(const trace_data& u) {
return peaks;
}
-util::optional<trace_peak> peak_delta(const trace_data& a, const trace_data& b) {
+util::optional<trace_peak> peak_delta(const trace_data<double>& a, const trace_data<double>& b) {
auto p = local_maxima(a);
auto q = local_maxima(b);
diff --git a/tests/validation/trace_analysis.hpp b/tests/validation/trace_analysis.hpp
index 8b0976237bc1790dc13141ee86c625cf196fc4d4..3bed1446512cc55a44195f9a415c5af4c9473fab 100644
--- a/tests/validation/trace_analysis.hpp
+++ b/tests/validation/trace_analysis.hpp
@@ -17,14 +17,14 @@ namespace mc {
// Compute max |v_i - f(t_i)| where (t, v) is the
// first trace `u` and f is the piece-wise linear interpolant
-// of the second trace `ref`.
+// of the second trace `r`.
-double linf_distance(const trace_data& u, const trace_data& ref);
+double linf_distance(const trace_data<double>& u, const trace_data<double>& r);
-// Compute linf distance as above, excluding samples near
-// times given in `excl`, monotonically increasing.
+// Compute linf distance as above, excluding samples in the first trace
+// near times given in `excl`, monotonically increasing.
-double linf_distance(const trace_data& u, const trace_data& ref, const std::vector<float>& excl);
+double linf_distance(const trace_data<double>& u, const trace_data<double>& r, const std::vector<float>& excl);
// Find local maxima (peaks) in a trace, excluding end points.
@@ -38,14 +38,14 @@ struct trace_peak {
}
};
-std::vector<trace_peak> local_maxima(const trace_data& u);
+std::vector<trace_peak> local_maxima(const trace_data<double>& u);
// Compare differences in peak times across two traces.
// Returns largest magnitute displacement between peaks,
// together with a sampling error bound, or `nothing`
// if the number of peaks differ.
-util::optional<trace_peak> peak_delta(const trace_data& a, const trace_data& b);
+util::optional<trace_peak> peak_delta(const trace_data<double>& a, const trace_data<double>& b);
// Record for error data for convergence testing.
// Only linf and peak_delta are used for convergence testing below;
diff --git a/tests/validation/validate.cpp b/tests/validation/validate.cpp
index e186a64f9ddcab1699ef45dae9a59128ae0a7a72..911e3067f3c4392bbb07a1d5bd6f430e5b1b8781 100644
--- a/tests/validation/validate.cpp
+++ b/tests/validation/validate.cpp
@@ -21,8 +21,8 @@ const char* usage_str =
" -p, --path=DIR Look for validation reference data in DIR\n"
" -m, --max-comp=N Run convergence tests to a maximum of N\n"
" compartments per segment\n"
-" -d, --min-dt=DT Run convergence tests with or with a minimumf\n"
-" timestep DT\n"
+" -d, --min-dt=DT Run convergence tests with a minimum timestep DT [ms]\n"
+" -s, --sample-dt=DT Sample rate for simulations [ms]\n"
" -h, --help Display usage information and exit\n"
"\n"
"Validation data is searched by default in the directory specified by\n"
@@ -52,6 +52,9 @@ int main(int argc, char **argv) {
else if (auto o = parse_opt<float>(arg, 'd', "min-dt")) {
g_trace_io.set_min_dt(*o);
}
+ else if (auto o = parse_opt<float>(arg, 's', "sample-dt")) {
+ g_trace_io.set_sample_dt(*o);
+ }
else if (auto o = parse_opt(arg, 'v', "verbose")) {
g_trace_io.set_verbose(true);
}
@@ -66,7 +69,6 @@ int main(int argc, char **argv) {
return RUN_ALL_TESTS();
}
-
catch (to::parse_opt_error& e) {
to::usage(argv[0], usage_str, e.what());
return 1;
diff --git a/tests/validation/validate_ball_and_stick.hpp b/tests/validation/validate_ball_and_stick.hpp
index 21dd6469fd21b4bc167226696c9e4e654efda3a8..516d82c77560d4987c9fe30d0925bee9818c1873 100644
--- a/tests/validation/validate_ball_and_stick.hpp
+++ b/tests/validation/validate_ball_and_stick.hpp
@@ -7,29 +7,40 @@
#include <hardware/node_info.hpp>
#include <model.hpp>
#include <recipe.hpp>
+#include <segment.hpp>
#include <simple_sampler.hpp>
+#include <util/meta.hpp>
#include <util/path.hpp>
#include "../gtest.h"
-#include "../test_common_cells.hpp"
+#include "../common_cells.hpp"
+#include "../simple_recipes.hpp"
#include "convergence_test.hpp"
#include "trace_analysis.hpp"
#include "validation_data.hpp"
-template <typename SamplerInfoSeq>
+#include <iostream>
+
+struct probe_point {
+ const char* label;
+ nest::mc::segment_location where;
+};
+
+template <typename ProbePointSeq>
void run_ncomp_convergence_test(
const char* model_name,
const nest::mc::util::path& ref_data_path,
nest::mc::backend_kind backend,
const nest::mc::cell& c,
- SamplerInfoSeq& samplers,
+ ProbePointSeq& probe_points,
float t_end=100.f)
{
using namespace nest::mc;
auto max_ncomp = g_trace_io.max_ncomp();
auto dt = g_trace_io.min_dt();
+ auto sample_dt = g_trace_io.sample_dt();
nlohmann::json meta = {
{"name", "membrane voltage"},
@@ -42,7 +53,14 @@ void run_ncomp_convergence_test(
auto exclude = stimulus_ends(c);
- convergence_test_runner<int> runner("ncomp", samplers, meta);
+ auto n_probe = util::size(probe_points);
+ std::vector<probe_label> plabels;
+ plabels.reserve(n_probe);
+ for (unsigned i = 0; i<n_probe; ++i) {
+ plabels.push_back(probe_label{probe_points[i].label, {0u, i}});
+ }
+
+ convergence_test_runner<int> runner("ncomp", plabels, meta);
runner.load_reference_data(ref_data_path);
for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) {
@@ -51,11 +69,16 @@ void run_ncomp_convergence_test(
seg->set_compartments(ncomp);
}
}
+ cable1d_recipe rec{c};
+ for (const auto& p: probe_points) {
+ rec.add_probe(0, 0, cell_probe_address{p.where, cell_probe_address::membrane_voltage});
+ }
+
hw::node_info nd(1, backend==backend_kind::gpu? 1: 0);
- auto decomp = partition_load_balance(singleton_recipe{c}, nd);
- model m(singleton_recipe{c}, decomp);
+ auto decomp = partition_load_balance(rec, nd);
+ model m(rec, decomp);
- runner.run(m, ncomp, t_end, dt, exclude);
+ runner.run(m, ncomp, sample_dt, t_end, dt, exclude);
}
runner.report();
runner.assert_all_convergence();
@@ -65,13 +88,10 @@ void validate_ball_and_stick(nest::mc::backend_kind backend) {
using namespace nest::mc;
cell c = make_cell_ball_and_stick();
- add_common_voltage_probes(c);
-
- float sample_dt = 0.025f;
- sampler_info samplers[] = {
- {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)},
- {"dend.mid", {0u, 1u}, simple_sampler(sample_dt)},
- {"dend.end", {0u, 2u}, simple_sampler(sample_dt)}
+ probe_point points[] = {
+ {"soma.mid", {0u, 0.5}},
+ {"dend.mid", {1u, 0.5}},
+ {"dend.end", {1u, 1.0}}
};
run_ncomp_convergence_test(
@@ -79,20 +99,17 @@ void validate_ball_and_stick(nest::mc::backend_kind backend) {
"neuron_ball_and_stick.json",
backend,
c,
- samplers);
+ points);
}
void validate_ball_and_taper(nest::mc::backend_kind backend) {
using namespace nest::mc;
cell c = make_cell_ball_and_taper();
- add_common_voltage_probes(c);
-
- float sample_dt = 0.025f;
- sampler_info samplers[] = {
- {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)},
- {"taper.mid", {0u, 1u}, simple_sampler(sample_dt)},
- {"taper.end", {0u, 2u}, simple_sampler(sample_dt)}
+ probe_point points[] = {
+ {"soma.mid", {0u, 0.5}},
+ {"taper.mid", {1u, 0.5}},
+ {"taper.end", {1u, 1.0}}
};
run_ncomp_convergence_test(
@@ -100,24 +117,21 @@ void validate_ball_and_taper(nest::mc::backend_kind backend) {
"neuron_ball_and_taper.json",
backend,
c,
- samplers);
+ points);
}
void validate_ball_and_3stick(nest::mc::backend_kind backend) {
using namespace nest::mc;
cell c = make_cell_ball_and_3stick();
- add_common_voltage_probes(c);
-
- float sample_dt = 0.025f;
- sampler_info samplers[] = {
- {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)},
- {"dend1.mid", {0u, 1u}, simple_sampler(sample_dt)},
- {"dend1.end", {0u, 2u}, simple_sampler(sample_dt)},
- {"dend2.mid", {0u, 3u}, simple_sampler(sample_dt)},
- {"dend2.end", {0u, 4u}, simple_sampler(sample_dt)},
- {"dend3.mid", {0u, 5u}, simple_sampler(sample_dt)},
- {"dend3.end", {0u, 6u}, simple_sampler(sample_dt)}
+ probe_point points[] = {
+ {"soma.mid", {0u, 0.5}},
+ {"dend1.mid", {1u, 0.5}},
+ {"dend1.end", {1u, 1.0}},
+ {"dend2.mid", {2u, 0.5}},
+ {"dend2.end", {2u, 1.0}},
+ {"dend3.mid", {3u, 0.5}},
+ {"dend3.end", {3u, 1.0}}
};
run_ncomp_convergence_test(
@@ -125,24 +139,17 @@ void validate_ball_and_3stick(nest::mc::backend_kind backend) {
"neuron_ball_and_3stick.json",
backend,
c,
- samplers);
+ points);
}
void validate_rallpack1(nest::mc::backend_kind backend) {
using namespace nest::mc;
cell c = make_cell_simple_cable();
-
- // three probes: left end, 30% along, right end.
- c.add_probe({{1, 0.0}, probeKind::membrane_voltage});
- c.add_probe({{1, 0.3}, probeKind::membrane_voltage});
- c.add_probe({{1, 1.0}, probeKind::membrane_voltage});
-
- float sample_dt = 0.025f;
- sampler_info samplers[] = {
- {"cable.x0.0", {0u, 0u}, simple_sampler(sample_dt)},
- {"cable.x0.3", {0u, 1u}, simple_sampler(sample_dt)},
- {"cable.x1.0", {0u, 2u}, simple_sampler(sample_dt)},
+ probe_point points[] = {
+ {"cable.x0.0", {1u, 0.0}},
+ {"cable.x0.3", {1u, 0.3}},
+ {"cable.x1.0", {1u, 1.0}}
};
run_ncomp_convergence_test(
@@ -150,7 +157,7 @@ void validate_rallpack1(nest::mc::backend_kind backend) {
"numeric_rallpack1.json",
backend,
c,
- samplers,
+ points,
250.f);
}
@@ -158,13 +165,10 @@ void validate_ball_and_squiggle(nest::mc::backend_kind backend) {
using namespace nest::mc;
cell c = make_cell_ball_and_squiggle();
- add_common_voltage_probes(c);
-
- float sample_dt = 0.025f;
- sampler_info samplers[] = {
- {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)},
- {"dend.mid", {0u, 1u}, simple_sampler(sample_dt)},
- {"dend.end", {0u, 2u}, simple_sampler(sample_dt)}
+ probe_point points[] = {
+ {"soma.mid", {0u, 0.5}},
+ {"dend.mid", {1u, 0.5}},
+ {"dend.end", {1u, 1.0}}
};
#if 0
@@ -183,5 +187,5 @@ void validate_ball_and_squiggle(nest::mc::backend_kind backend) {
"neuron_ball_and_squiggle.json",
backend,
c,
- samplers);
+ points);
}
diff --git a/tests/validation/validate_compartment_policy.cpp b/tests/validation/validate_compartment_policy.cpp
index 4848ea12a14c9da2bce2dff0fd28ff4f042b4bd5..3c78c3a2e415be7aa1ce51c1e1a77b0d4a87883e 100644
--- a/tests/validation/validate_compartment_policy.cpp
+++ b/tests/validation/validate_compartment_policy.cpp
@@ -13,8 +13,10 @@
#include "../gtest.h"
-#include "../test_common_cells.hpp"
+#include "../common_cells.hpp"
+#include "../simple_recipes.hpp"
#include "../test_util.hpp"
+
#include "trace_analysis.hpp"
#include "validation_data.hpp"
diff --git a/tests/validation/validate_kinetic.hpp b/tests/validation/validate_kinetic.hpp
index 75a4d5efe4d235e1ae90f1d9a24d71bf5d3f0bfb..919324932ed5b00afb7de3e873c57f1c19fc0c70 100644
--- a/tests/validation/validate_kinetic.hpp
+++ b/tests/validation/validate_kinetic.hpp
@@ -10,7 +10,9 @@
#include <simple_sampler.hpp>
#include <util/rangeutil.hpp>
-#include "../test_common_cells.hpp"
+#include "../common_cells.hpp"
+#include "../simple_recipes.hpp"
+
#include "convergence_test.hpp"
#include "trace_analysis.hpp"
#include "validation_data.hpp"
@@ -24,19 +26,21 @@ void run_kinetic_dt(
{
using namespace nest::mc;
- float sample_dt = .025f;
- sampler_info samplers[] = {
- {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)}
- };
+ float sample_dt = g_trace_io.sample_dt();
+
+ cable1d_recipe rec{c};
+ rec.add_probe(0, 0, cell_probe_address{{0, 0.5}, cell_probe_address::membrane_voltage});
+ probe_label plabels[1] = {"soma.mid", {0u, 0u}};
meta["sim"] = "nestmc";
meta["backend_kind"] = to_string(backend);
- convergence_test_runner<float> runner("dt", samplers, meta);
+
+ convergence_test_runner<float> runner("dt", plabels, meta);
runner.load_reference_data(ref_file);
hw::node_info nd(1, backend==backend_kind::gpu? 1: 0);
- auto decomp = partition_load_balance(singleton_recipe{c}, nd);
- model model(singleton_recipe{c}, decomp);
+ auto decomp = partition_load_balance(rec, nd);
+ model model(rec, decomp);
auto exclude = stimulus_ends(c);
@@ -49,7 +53,7 @@ void run_kinetic_dt(
model.reset();
float dt = float(1./oo_dt);
- runner.run(model, dt, t_end, dt, exclude);
+ runner.run(model, dt, sample_dt, t_end, dt, exclude);
}
}
@@ -64,7 +68,6 @@ void validate_kinetic_kin1(nest::mc::backend_kind backend) {
// 20 µm diameter soma with single mechanism, current probe
cell c;
auto soma = c.add_soma(10);
- c.add_probe({{0, 0.5}, probeKind::membrane_current});
soma->add_mechanism(std::string("test_kin1"));
nlohmann::json meta = {
@@ -82,7 +85,6 @@ void validate_kinetic_kinlva(nest::mc::backend_kind backend) {
// 20 µm diameter soma with single mechanism, current probe
cell c;
auto soma = c.add_soma(10);
- c.add_probe({{0, 0.5}, probeKind::membrane_voltage});
c.add_stimulus({0,0.5}, {20., 130., -0.025});
soma->add_mechanism(std::string("test_kinlva"));
diff --git a/tests/validation/validate_soma.hpp b/tests/validation/validate_soma.hpp
index b33fa7b80959555a728c57d393fe42915d1f74bc..edc985b84ca734dbccf6024216c91f3c1538ff98 100644
--- a/tests/validation/validate_soma.hpp
+++ b/tests/validation/validate_soma.hpp
@@ -10,7 +10,9 @@
#include <simple_sampler.hpp>
#include <util/rangeutil.hpp>
-#include "../test_common_cells.hpp"
+#include "../common_cells.hpp"
+#include "../simple_recipes.hpp"
+
#include "convergence_test.hpp"
#include "trace_analysis.hpp"
#include "validation_data.hpp"
@@ -18,15 +20,17 @@
void validate_soma(nest::mc::backend_kind backend) {
using namespace nest::mc;
+ float sample_dt = g_trace_io.sample_dt();
+
cell c = make_cell_soma_only();
- add_common_voltage_probes(c);
- hw::node_info nd(1, backend==backend_kind::gpu? 1: 0);
- auto decomp = partition_load_balance(singleton_recipe{c}, nd);
- model m(singleton_recipe{c}, decomp);
+ cable1d_recipe rec{c};
+ rec.add_probe(0, 0, cell_probe_address{{0, 0.5}, cell_probe_address::membrane_voltage});
+ probe_label plabels[1] = {"soma.mid", {0u, 0u}};
- float sample_dt = .025f;
- sampler_info samplers[] = {{"soma.mid", {0u, 0u}, simple_sampler(sample_dt)}};
+ hw::node_info nd(1, backend==backend_kind::gpu? 1: 0);
+ auto decomp = partition_load_balance(rec, nd);
+ model m(rec, decomp);
nlohmann::json meta = {
{"name", "membrane voltage"},
@@ -36,7 +40,7 @@ void validate_soma(nest::mc::backend_kind backend) {
{"backend_kind", to_string(backend)}
};
- convergence_test_runner<float> runner("dt", samplers, meta);
+ convergence_test_runner<float> runner("dt", plabels, meta);
runner.load_reference_data("numeric_soma.json");
float t_end = 100.f;
@@ -50,7 +54,7 @@ void validate_soma(nest::mc::backend_kind backend) {
m.reset();
float dt = float(1./oo_dt);
- runner.run(m, dt, t_end, dt, {});
+ runner.run(m, dt, sample_dt, t_end, dt, {});
}
}
end:
diff --git a/tests/validation/validate_synapses.hpp b/tests/validation/validate_synapses.hpp
index 3e814add278ccb91cad9c8722c1c23261a6d8c6b..429da370fb0ed8cc164cc8e3c953e46591c2162c 100644
--- a/tests/validation/validate_synapses.hpp
+++ b/tests/validation/validate_synapses.hpp
@@ -12,7 +12,9 @@
#include "../gtest.h"
-#include "../test_common_cells.hpp"
+#include "../common_cells.hpp"
+#include "../simple_recipes.hpp"
+
#include "convergence_test.hpp"
#include "trace_analysis.hpp"
#include "validation_data.hpp"
@@ -38,7 +40,6 @@ void run_synapse_test(
cell c = make_cell_ball_and_stick(false); // no stimuli
parameter_list syn_default(syn_type);
c.add_synapse({1, 0.5}, syn_default);
- add_common_voltage_probes(c);
// injected spike events
std::vector<postsynaptic_spike_event> synthetic_events = {
@@ -50,24 +51,33 @@ void run_synapse_test(
// exclude points of discontinuity from linf analysis
std::vector<float> exclude = {10.f, 20.f, 40.f};
- float sample_dt = 0.025f;
- sampler_info samplers[] = {
- {"soma.mid", {0u, 0u}, simple_sampler(sample_dt)},
- {"dend.mid", {0u, 1u}, simple_sampler(sample_dt)},
- {"dend.end", {0u, 2u}, simple_sampler(sample_dt)}
+ float sample_dt = g_trace_io.sample_dt();
+ probe_label plabels[3] = {
+ {"soma.mid", {0u, 0u}},
+ {"dend.mid", {0u, 1u}},
+ {"dend.end", {0u, 2u}}
};
- convergence_test_runner<int> runner("ncomp", samplers, meta);
+ convergence_test_runner<int> runner("ncomp", plabels, meta);
runner.load_reference_data(ref_data_path);
hw::node_info nd(1, backend==backend_kind::gpu? 1: 0);
for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) {
c.cable(1)->set_compartments(ncomp);
- auto decomp = partition_load_balance(singleton_recipe{c}, nd);
- model m(singleton_recipe{c}, decomp);
+
+ cable1d_recipe rec{c};
+ // soma.mid
+ rec.add_probe(0, 0, cell_probe_address{{0, 0.5}, cell_probe_address::membrane_voltage});
+ // dend.mid
+ rec.add_probe(0, 0, cell_probe_address{{1, 0.5}, cell_probe_address::membrane_voltage});
+ // dend.end
+ rec.add_probe(0, 0, cell_probe_address{{1, 1.0}, cell_probe_address::membrane_voltage});
+
+ auto decomp = partition_load_balance(rec, nd);
+ model m(rec, decomp);
m.group(0).enqueue_events(synthetic_events);
- runner.run(m, ncomp, t_end, dt, exclude);
+ runner.run(m, ncomp, sample_dt, t_end, dt, exclude);
}
runner.report();
runner.assert_all_convergence();
diff --git a/tests/validation/validation_data.cpp b/tests/validation/validation_data.cpp
index 7b44713bf96a6b7e03c054612f398fe35dade33b..7bb9d86152fd5a580e231bae5a86174bd59a6bdf 100644
--- a/tests/validation/validation_data.cpp
+++ b/tests/validation/validation_data.cpp
@@ -48,11 +48,11 @@ util::path trace_io::find_datadir() {
return util::path();
}
-void trace_io::save_trace(const std::string& label, const trace_data& data, const nlohmann::json& meta) {
+void trace_io::save_trace(const std::string& label, const trace_data<double>& data, const nlohmann::json& meta) {
save_trace("time", label, data, meta);
}
-void trace_io::save_trace(const std::string& abscissa, const std::string& label, const trace_data& data, const nlohmann::json& meta) {
+void trace_io::save_trace(const std::string& abscissa, const std::string& label, const trace_data<double>& data, const nlohmann::json& meta) {
using namespace nest::mc;
nlohmann::json j = meta;
@@ -65,8 +65,8 @@ void trace_io::save_trace(const std::string& abscissa, const std::string& label,
}
template <typename Seq1, typename Seq2>
-static trace_data zip_trace_data(const Seq1& ts, const Seq2& vs) {
- trace_data trace;
+static trace_data<double> zip_trace_data(const Seq1& ts, const Seq2& vs) {
+ trace_data<double> trace;
auto ti = std::begin(ts);
auto te = std::end(ts);
@@ -79,7 +79,7 @@ static trace_data zip_trace_data(const Seq1& ts, const Seq2& vs) {
return trace;
}
-static void parse_trace_json(const nlohmann::json& j, std::map<std::string, trace_data>& traces) {
+static void parse_trace_json(const nlohmann::json& j, std::map<std::string, trace_data<double>>& traces) {
if (j.is_array()) {
for (auto& i: j) parse_trace_json(i, traces);
}
@@ -94,7 +94,7 @@ static void parse_trace_json(const nlohmann::json& j, std::map<std::string, trac
}
}
-std::map<std::string, trace_data> trace_io::load_traces(const util::path& name) {
+std::map<std::string, trace_data<double>> trace_io::load_traces(const util::path& name) {
util::path file = datadir_/name;
std::ifstream fid(file);
if (!fid) {
@@ -104,7 +104,7 @@ std::map<std::string, trace_data> trace_io::load_traces(const util::path& name)
nlohmann::json data;
fid >> data;
- std::map<std::string, trace_data> traces;
+ std::map<std::string, trace_data<double>> traces;
parse_trace_json(data, traces);
return traces;
}
diff --git a/tests/validation/validation_data.hpp b/tests/validation/validation_data.hpp
index a470ebf47f8c47c5dec1d9088d80c435814eae04..e462b5d8b2fbecc69652f44d481c9f45953d9a83 100644
--- a/tests/validation/validation_data.hpp
+++ b/tests/validation/validation_data.hpp
@@ -23,6 +23,9 @@ namespace mc {
* output JSON file, if specified.
*/
+// TODO: split simulation options from trace_io structure; rename
+// file to e.g. 'trace_io.hpp'
+
class trace_io {
public:
// Try to find the data directory on construction.
@@ -42,9 +45,9 @@ public:
}
}
- void save_trace(const std::string& label, const trace_data& data, const nlohmann::json& meta);
- void save_trace(const std::string& abscissa, const std::string& label, const trace_data& data, const nlohmann::json& meta);
- std::map<std::string, trace_data> load_traces(const util::path& name);
+ void save_trace(const std::string& label, const trace_data<double>& data, const nlohmann::json& meta);
+ void save_trace(const std::string& abscissa, const std::string& label, const trace_data<double>& data, const nlohmann::json& meta);
+ std::map<std::string, trace_data<double>> load_traces(const util::path& name);
// common flags, options set by driver
@@ -57,6 +60,9 @@ public:
void set_min_dt(float dt) { min_dt_ = dt; }
float min_dt() const { return min_dt_; }
+ void set_sample_dt(float dt) { sample_dt_ = dt; }
+ float sample_dt() const { return sample_dt_; }
+
void set_datadir(const util::path& dir) { datadir_ = dir; }
void set_output(const util::path& file) {
@@ -81,6 +87,7 @@ private:
bool verbose_flag_ = false;
int max_ncomp_ = 100;
float min_dt_ = 0.001f;
+ float sample_dt_ = 0.005f;
// Returns value of NMC_DATADIR environment variable if set,
// otherwise make a 'best-effort' search for the data directory,