diff --git a/miniapp/miniapp.cpp b/miniapp/miniapp.cpp
index d99e43e5bdeba33439f044cd2c3aa2efe657c16c..4a619cf9a087dc7fa151ee9613adc6d5ff6357c4 100644
--- a/miniapp/miniapp.cpp
+++ b/miniapp/miniapp.cpp
@@ -1,4 +1,6 @@
#include <iostream>
+#include <fstream>
+#include <sstream>
#include <cell.hpp>
#include <cell_group.hpp>
@@ -10,6 +12,7 @@
#include "profiling/profiler.hpp"
#include "communication/communicator.hpp"
#include "communication/global_policy.hpp"
+#include "util/optional.hpp"
using namespace nest;
@@ -23,6 +26,8 @@ using global_policy = nest::mc::communication::global_policy;
using communicator_type =
mc::communication::communicator<global_policy>;
+using nest::mc::util::optional;
+
struct model {
communicator_type communicator;
std::vector<cell_group> cell_groups;
@@ -92,6 +97,63 @@ struct model {
// TODO : only stored here because init_communicator() and update_gids() are split
std::vector<id_type> source_map;
std::vector<id_type> target_map;
+
+ // traces from probes
+ struct trace_data {
+ using sample_type = std::pair<float,double>;
+ std::string name;
+ index_type id;
+ std::vector<sample_type> samples;
+ };
+
+ // different traces may be written to by different threads;
+ // during simulation, each trace_sampler will be responsible for its
+ // corresponding element in the traces vector.
+
+ std::vector<trace_data> traces;
+
+ // make a sampler that records to traces
+ struct simple_sampler_functor {
+ std::vector<trace_data> &traces_;
+ size_t trace_index_ = 0;
+ float requested_sample_time_ = 0;
+ float dt_ = 0;
+
+ simple_sampler_functor(std::vector<trace_data> &traces, size_t index, float dt) :
+ traces_(traces), trace_index_(index), dt_(dt) {}
+
+ optional<float> operator()(float t, double v) {
+ traces_[trace_index_].samples.push_back({t,v});
+ return requested_sample_time_ += dt_;
+ }
+ };
+
+ mc::sampler make_simple_sampler(index_type probe_gid, const std::string name,
+ index_type id, float dt)
+ {
+ traces.push_back(trace_data{name, id});
+ return {probe_gid, simple_sampler_functor(traces, traces.size()-1, dt)};
+ }
+
+ void reset_traces() {
+ // do not call during simulation: thread-unsafe access to traces.
+ traces.clear();
+ }
+
+ void dump_traces() {
+ // do not call during simulation: thread-unsafe access to traces.
+ for (const auto& trace: traces) {
+ std::stringstream path;
+ path << "trace_" << trace.id << "_" << trace.name << ".dat";
+
+ std::ofstream file(path.str());
+ file << "time\t" << trace.name << "\n";
+ for (const auto& sample: trace.samples) {
+ file << sample.first << "\t" << sample.second << "\n";
+ }
+ file.close();
+ }
+ }
};
// define some global model parameters
@@ -166,6 +228,7 @@ int main(int argc, char** argv) {
if (!id) {
std::cout << "there were " << m.communicator.num_spikes() << " spikes\n";
}
+ m.dump_traces();
#ifdef SPLAT
if (!global_policy::id()) {
@@ -254,7 +317,23 @@ void all_to_all_model(nest::mc::io::options& opt, model& m) {
//
m.init_communicator();
- mc::util::profiler_enter("setup", "connections");
+ // monitor soma and dendrite on a few cells
+ float sample_dt = 0.1;
+ index_type monitor_group_gids[] = { 0, 1, 2 };
+ for (auto gid : monitor_group_gids) {
+ if (!m.communicator.is_local_group(gid)) {
+ continue;
+ }
+
+ auto lid = m.communicator.group_lid(gid);
+ auto probe_soma = m.cell_groups[lid].probe_gid_range().first;
+ auto probe_dend = probe_soma+1;
+
+ m.cell_groups[lid].add_sampler(m.make_simple_sampler(probe_soma, "vsoma", gid, sample_dt));
+ m.cell_groups[lid].add_sampler(m.make_simple_sampler(probe_dend, "vdend", gid, sample_dt));
+ }
+
+ mc::util::profiler_enter("setup", "connections");
// lid is local cell/group id
for (auto lid=0u; lid<ncell_local; ++lid) {
auto target = m.communicator.target_gid_from_group_lid(lid);
@@ -271,7 +350,7 @@ void all_to_all_model(nest::mc::io::options& opt, model& m) {
}
m.communicator.construct();
- mc::util::profiler_leave(2);
+ mc::util::profiler_leave(2);
m.update_gids();
}
@@ -323,6 +402,14 @@ mc::cell make_cell(int compartments_per_segment, int num_synapses) {
cell.add_synapse({1, 0.5});
}
+ // add probes:
+ auto probe_soma = cell.add_probe({0, 0}, mc::cell::membrane_voltage);
+ auto probe_dendrite = cell.add_probe({1, 0.5}, mc::cell::membrane_voltage);
+
+ EXPECTS(probe_soma==0);
+ EXPECTS(probe_dendrite==1);
+ (void)probe_soma, (void)probe_dendrite;
+
return cell;
}
diff --git a/src/cell.hpp b/src/cell.hpp
index 0e13e588d42f7d213a7aa6f4928f1d043f95e07d..f53fb0592c94d76c3dfcd7552814345fde28f874 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -40,13 +40,15 @@ int find_compartment_index(
/// high-level abstract representation of a cell and its segments
class cell {
- public:
+public:
// types
using index_type = int;
using value_type = double;
using point_type = point<value_type>;
+ enum probe_sort { membrane_voltage, membrane_current };
+
// constructor
cell();
@@ -139,7 +141,18 @@ class cell {
return spike_detectors_;
}
- private:
+ //////////////////
+ // probes
+ //////////////////
+ index_type add_probe(segment_location loc, enum probe_sort sort) {
+ probes_.push_back({loc, sort});
+ return probes_.size()-1;
+ }
+
+ const std::vector<std::pair<segment_location, enum probe_sort>>&
+ probes() const { return probes_; }
+
+private:
// storage for connections
std::vector<index_type> parents_;
@@ -155,6 +168,9 @@ class cell {
// the sensors
std::vector<std::pair<segment_location, double>> spike_detectors_;
+
+ // the probes
+ std::vector<std::pair<segment_location, enum probe_sort>> probes_;
};
// Checks that two cells have the same
diff --git a/src/cell_group.hpp b/src/cell_group.hpp
index 4964ad10a2328c96593f9e6b07b136baa7ca7fcd..922b7b197091af6cf4ba7ecc345030bcff3e380e 100644
--- a/src/cell_group.hpp
+++ b/src/cell_group.hpp
@@ -1,5 +1,6 @@
#pragma once
+#include <cstdint>
#include <vector>
#include <cell.hpp>
@@ -12,10 +13,21 @@
namespace nest {
namespace mc {
+// samplers take a time and sample value, and return an optional time
+// for the next desired sample.
+
+struct sampler {
+ using index_type = int;
+ using time_type = float;
+ using value_type = double;
+
+ index_type probe_gid; // samplers are attached to probes
+ std::function<util::optional<time_type>(time_type, value_type)> sample;
+};
+
template <typename Cell>
class cell_group {
- public :
-
+public:
using index_type = uint32_t;
using cell_type = Cell;
using value_type = typename cell_type::value_type;
@@ -52,6 +64,18 @@ class cell_group {
first_target_gid_ = lid;
}
+ index_type num_probes() const {
+ return cell_.num_probes();
+ }
+
+ void set_probe_gids(index_type gid) {
+ first_probe_gid_ = gid;
+ }
+
+ std::pair<index_type, index_type> probe_gid_range() const {
+ return { first_probe_gid_, first_probe_gid_+cell_.num_probes() };
+ }
+
#ifdef SPLAT
void splat(std::string fname) {
char buffer[128];
@@ -65,6 +89,20 @@ class cell_group {
void advance(double tfinal, double dt) {
while (cell_.time()<tfinal) {
+ // take any pending samples
+ float cell_time = cell_.time();
+ while (auto m = sample_events_.pop_if_before(cell_time)) {
+ auto &sampler = samplers_[m->sampler_index];
+ EXPECTS((bool)sampler.sample);
+
+ index_type probe_index = sampler.probe_gid-first_probe_gid_;
+ auto next = sampler.sample(cell_.time(), cell_.probe(probe_index));
+ if (next) {
+ m->time = std::max(*next, cell_time);
+ sample_events_.push(*m);
+ }
+ }
+
#ifdef SPLAT
tt.push_back(cell_.time());
vs.push_back(cell_.voltage({0,0.0}));
@@ -126,7 +164,13 @@ class cell_group {
spikes_.clear();
}
- private :
+ void add_sampler(const sampler &s, float start_time = 0) {
+ unsigned sampler_index = samplers_.size();
+ samplers_.push_back(s);
+ sample_events_.push({sampler_index, start_time});
+ }
+
+private:
#ifdef SPLAT
// REMOVE as soon as we have a better way to probe cell state
@@ -145,10 +189,19 @@ class cell_group {
std::vector<communication::spike<index_type>> spikes_;
/// pending events to be delivered
- event_queue events_;
+ event_queue<postsynaptic_spike_event> events_;
+
+ /// pending samples to be taken
+ event_queue<sample_event> sample_events_;
/// the global id of the first target (e.g. a synapse) in this group
index_type first_target_gid_;
+
+ /// the global id of the first probe in this group
+ index_type first_probe_gid_;
+
+ /// collection of samplers to be run against probes in this group
+ std::vector<sampler> samplers_;
};
} // namespace mc
diff --git a/src/communication/communicator.hpp b/src/communication/communicator.hpp
index 517a2f6c3884c66df22914776849fcfd2da7e136..34185ac0fbd795294211ced7190d8b84a3e70874 100644
--- a/src/communication/communicator.hpp
+++ b/src/communication/communicator.hpp
@@ -90,11 +90,17 @@ public:
}
id_type target_lid(id_type gid) {
- EXPECTS(is_local_group(gid));
+ EXPECTS(is_local_group(gid));
return gid - target_gid_map_[domain_id()];
}
+ id_type group_lid(id_type gid) {
+ EXPECTS(is_local_group(gid));
+
+ return gid - group_gid_map_[domain_id()];
+ }
+
// builds the optimized data structure
void construct() {
if (!std::is_sorted(connections_.begin(), connections_.end())) {
@@ -188,7 +194,7 @@ public:
return communication_policy_.size();
}
- const std::vector<local_event>& queue(int i) const {
+ const std::vector<postsynaptic_spike_event>& queue(int i) const {
return events_[i];
}
@@ -234,7 +240,7 @@ private:
local_spike_store_type thread_spikes_;
std::vector<connection> connections_;
- std::vector<std::vector<nest::mc::local_event>> events_;
+ std::vector<std::vector<postsynaptic_spike_event>> events_;
// local target group i has targets in the half open range
// [target_map_[i], target_map_[i+1])
diff --git a/src/communication/connection.hpp b/src/communication/connection.hpp
index 528c57b59fddbfc4b7f3c96772c5c0c04685758b..ee1acfd72dc21cffe01f28239bca29526366df31 100644
--- a/src/communication/connection.hpp
+++ b/src/communication/connection.hpp
@@ -25,7 +25,7 @@ public:
id_type source() const { return source_; }
id_type destination() const { return destination_; }
- local_event make_event(spike<id_type> s) {
+ postsynaptic_spike_event make_event(spike<id_type> s) {
return {destination_, s.time + delay_, weight_};
}
diff --git a/src/event_queue.hpp b/src/event_queue.hpp
index 989e34873a1f3ff851d0a3a9f85276cea9d9bf7a..7f3034f3f08bbb2235bf457c95e8eedfe8d0895d 100644
--- a/src/event_queue.hpp
+++ b/src/event_queue.hpp
@@ -9,22 +9,31 @@
namespace nest {
namespace mc {
-struct local_event {
+struct postsynaptic_spike_event {
uint32_t target;
float time;
float weight;
};
-inline bool operator < (local_event const& l, local_event const& r) {
- return l.time < r.time;
-}
+inline float event_time(const postsynaptic_spike_event &ev) { return ev.time; }
-inline bool operator > (local_event const& l, local_event const& r) {
- return l.time > r.time;
-}
+struct sample_event {
+ uint32_t sampler_index;
+ float time;
+};
+inline float event_time(const sample_event &ev) { return ev.time; }
+
+/* Event objects must have a method event_time() which returns a value
+ * from a type with a total ordering with respect to <, >, etc.
+ */
+
+template <typename Event>
class event_queue {
public :
+ using value_type = Event;
+ using time_type = decltype(event_time(std::declval<Event>()));
+
// create
event_queue() {}
@@ -37,7 +46,7 @@ public :
}
// push thing
- void push(local_event e) {
+ void push(const value_type &e) {
queue_.push(e);
}
@@ -46,8 +55,8 @@ public :
}
// pop until
- util::optional<local_event> pop_if_before(float t_until) {
- if (!queue_.empty() && queue_.top().time < t_until) {
+ util::optional<value_type> pop_if_before(time_type t_until) {
+ if (!queue_.empty() && event_time(queue_.top()) < t_until) {
auto ev = queue_.top();
queue_.pop();
return ev;
@@ -58,10 +67,16 @@ public :
}
private:
+ struct event_greater {
+ bool operator()(const Event &a, const Event &b) {
+ return event_time(a) > event_time(b);
+ }
+ };
+
std::priority_queue<
- local_event,
- std::vector<local_event>,
- std::greater<local_event>
+ Event,
+ std::vector<Event>,
+ event_greater
> queue_;
};
@@ -69,7 +84,7 @@ private:
} // namespace mc
inline
-std::ostream& operator<< (std::ostream& o, nest::mc::local_event e)
+std::ostream& operator<< (std::ostream& o, const nest::mc::postsynaptic_spike_event& e)
{
return o << "event[" << e.target << "," << e.time << "," << e.weight << "]";
}
diff --git a/src/fvm_cell.hpp b/src/fvm_cell.hpp
index c89dcdf0c3e920b4aa370f71a48767f4e20d1040..947b7a7ba233340274f5113c12d075d7851a1ea1 100644
--- a/src/fvm_cell.hpp
+++ b/src/fvm_cell.hpp
@@ -29,7 +29,7 @@ namespace fvm {
template <typename T, typename I>
class fvm_cell {
- public :
+public:
fvm_cell() = default;
@@ -116,7 +116,7 @@ class fvm_cell {
void advance(value_type dt);
/// pass an event to the appropriate synapse and call net_receive
- void apply_event(local_event e) {
+ void apply_event(postsynaptic_spike_event e) {
mechanisms_[synapse_index_]->net_receive(e.target, e.weight);
}
@@ -135,7 +135,14 @@ class fvm_cell {
value_type time() const { return t_; }
- private:
+ value_type probe(uint32_t i) const {
+ auto p = probes_[i];
+ return (this->*p.first)[p.second];
+ }
+
+ std::size_t num_probes() const { return probes_.size(); }
+
+private:
/// current time
value_type t_ = value_type{0};
@@ -179,8 +186,7 @@ class fvm_cell {
std::vector<std::pair<uint32_t, i_clamp>> stimulii_;
- /// event queue
- event_queue events_;
+ std::vector<std::pair<const vector_type fvm_cell::*, uint32_t>> probes_;
};
////////////////////////////////////////////////////////////////////////////////
@@ -397,6 +403,21 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
synapse_index_ = mechanisms_.size()-1;
// don't forget to give point processes access to cv_areas_
mechanisms_[synapse_index_]->set_areas(cv_areas_);
+
+ // record probe locations by index into corresponding state vector
+ for (auto probe : cell.probes()) {
+ uint32_t comp = find_compartment_index(probe.first, graph);
+ switch (probe.second) {
+ case nest::mc::cell::membrane_voltage:
+ probes_.push_back({&fvm_cell::voltage_, comp});
+ break;
+ case nest::mc::cell::membrane_current:
+ probes_.push_back({&fvm_cell::current_, comp});
+ break;
+ default:
+ throw std::logic_error("unrecognized probe sort");
+ }
+ }
}
template <typename T, typename I>
diff --git a/tests/test_event_queue.cpp b/tests/test_event_queue.cpp
index e67ade0bea22a417ca417a78c550b67bb0c3e45b..5ac4c75747d761638b0bbd4b90a9b6fed6fa3c18 100644
--- a/tests/test_event_queue.cpp
+++ b/tests/test_event_queue.cpp
@@ -7,8 +7,9 @@
TEST(event_queue, push)
{
using namespace nest::mc;
+ using ps_event_queue = event_queue<postsynaptic_spike_event>;
- event_queue q;
+ ps_event_queue q;
q.push({1u, 2.f, 2.f});
q.push({4u, 1.f, 2.f});
@@ -30,15 +31,16 @@ TEST(event_queue, push)
TEST(event_queue, push_range)
{
using namespace nest::mc;
+ using ps_event_queue = event_queue<postsynaptic_spike_event>;
- local_event events[] = {
+ postsynaptic_spike_event events[] = {
{1u, 2.f, 2.f},
{4u, 1.f, 2.f},
{8u, 20.f, 2.f},
{2u, 8.f, 2.f}
};
- event_queue q;
+ ps_event_queue q;
q.push(std::begin(events), std::end(events));
std::vector<float> times;
@@ -54,15 +56,16 @@ TEST(event_queue, push_range)
TEST(event_queue, pop_if_before)
{
using namespace nest::mc;
+ using ps_event_queue = event_queue<postsynaptic_spike_event>;
- local_event events[] = {
+ postsynaptic_spike_event events[] = {
{1u, 1.f, 2.f},
{2u, 2.f, 2.f},
{3u, 3.f, 2.f},
{4u, 4.f, 2.f}
};
- event_queue q;
+ ps_event_queue q;
q.push(std::begin(events), std::end(events));
EXPECT_EQ(q.size(), 4u);
diff --git a/tests/test_optional.cpp b/tests/test_optional.cpp
index 6755525af3cfb70d610f29dc5fc0c525f68856be..9c6f35cfa658ab1c45cee53250eaa747d89dea25 100644
--- a/tests/test_optional.cpp
+++ b/tests/test_optional.cpp
@@ -282,7 +282,7 @@ TEST(optionalm,conversion) {
TEST(optionalm,or_operator) {
optional<const char *> default_msg="default";
- auto x=nullptr | default_msg;
+ auto x=(char *)0 | default_msg;
EXPECT_TRUE((bool)x);
EXPECT_STREQ("default",x.get());