diff --git a/miniapp/miniapp.cpp b/miniapp/miniapp.cpp
index fdac119896bdeadf9b9d40fffd05bb3a6345bb5b..549386abd307db798127666913558ecf905a212d 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>
@@ -11,6 +13,7 @@
#include "communication/communicator.hpp"
#include "communication/serial_global_policy.hpp"
#include "communication/mpi_global_policy.hpp"
+#include "util/optional.hpp"
using namespace nest;
@@ -27,6 +30,8 @@ using communicator_type =
mc::communication::communicator<mc::communication::serial_global_policy>;
#endif
+using nest::mc::util::optional;
+
struct model {
communicator_type communicator;
std::vector<cell_group> cell_groups;
@@ -102,18 +107,62 @@ struct model {
std::vector<id_type> source_map;
std::vector<id_type> target_map;
- // sampling // WIP COME BACK HERE
- struct simple_sampler {
- index_type group_gid;
+ // traces from probes
+ struct trace_data {
+ using sample_type = std::pair<float,double>;
std::string name;
- double dt;
+ 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_;
+ }
+ };
- std::vector<std::pair<float,double>> tv;
- util::optional<float> operator()(float t, double v) {
- define
+ 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
@@ -185,6 +234,7 @@ int main(int argc, char** argv) {
m.print_times();
std::cout << "there were " << m.communicator.num_spikes() << " spikes\n";
}
+ m.dump_traces();
#ifdef SPLAT
if (!mc::mpi::rank()) {
@@ -284,21 +334,18 @@ void all_to_all_model(nest::mc::io::options& opt, model& m) {
// monitor soma and dendrite on a few cells
float sample_dt = 0.1;
- index_type monitor_group_gids = { 0, 1, 2 };
+ index_type monitor_group_gids[] = { 0, 1, 2 };
for (auto gid : monitor_group_gids) {
if (!m.communicator.is_local_group(gid)) {
continue;
}
- lid = m.communicator.group_lid(gid);
- index_type probe_soma = m.cell_groups[lid].probe_gid_range().first;
- index_type probe_dend = probe_soma+1;
+ auto lid = m.communicator.group_lid(gid);
+ auto probe_soma = m.cell_groups[lid].probe_gid_range().first;
+ auto probe_dend = probe_soma+1;
- // WIP COME BACK HERE
- m.cell_groups[lid].samplers = {
- { probe_soma, m.make_simple_sampler(gid, "vsoma", sample_dt) },
- { probe_dend, m.make_simple_sampler(gid, "vdend", sample_dt) }
- };
+ 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));
}
// lid is local cell/group id
@@ -381,11 +428,12 @@ mc::cell make_cell(int compartments_per_segment, int num_synapses) {
}
// add probes:
- auto probe_soma = cell.add_probe(nest::mc::cell::membrane_potential, {0,0});
- auto probe_dendrite = cell.add_probe(nest::mc::cell::membrane_potential, {1,0.5});
+ 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);
- EXPECT(probe_soma==0);
- EXPECT(probe_dendrite==1);
+ 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 f1360ffe03c6fd344191955ed7b9eec19b2d9b3f..7b6f5bedd1b80705631a5a82cb965dc227e9768b 100644
--- a/src/cell_group.hpp
+++ b/src/cell_group.hpp
@@ -1,5 +1,6 @@
#pragma once
+#include <cstdint>
#include <vector>
#include <cell.hpp>
@@ -14,6 +15,7 @@ namespace mc {
// for the next desired sample.
struct sampler {
+ using index_type = int;
using time_type = float;
using value_type = double;
@@ -86,16 +88,16 @@ public:
void advance(double tfinal, double dt) {
while (cell_.time()<tfinal) {
// take any pending samples
- while (sample_event m = sample_events_.pop_if_before(cell_.time())) {
- auto &sampler = samplers_[m.sampler_index];
- EXPECT((bool)sampler.sample);
+ 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) {
- EXPECT(*next>m.time);
- m.time = *next;
- sample_events_.push(m);
+ m->time = std::max(*next, cell_time);
+ sample_events_.push(*m);
}
}
@@ -158,8 +160,12 @@ public:
spikes_.clear();
}
- std::vector<sampler> samplers;
-
+ 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
@@ -189,6 +195,9 @@ private:
/// 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 eee922541e41d3d3d8b317660ab5b390b42d99fa..34185ac0fbd795294211ced7190d8b84a3e70874 100644
--- a/src/communication/communicator.hpp
+++ b/src/communication/communicator.hpp
@@ -194,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];
}
@@ -240,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/communication/mpi.hpp b/src/communication/mpi.hpp
index 063e36ef3fc3ae44f952e2caeddf18663832cefd..5be71b6381bf8311bcab85d5396d64ba2c3eb26b 100644
--- a/src/communication/mpi.hpp
+++ b/src/communication/mpi.hpp
@@ -43,7 +43,7 @@ namespace mpi {
template <> \
struct mpi_traits<T> { \
constexpr static size_t count() { return 1; } \
- constexpr static MPI_Datatype mpi_type() { return M; } \
+ /* constexpr */ static MPI_Datatype mpi_type() { return M; } \
constexpr static bool is_mpi_native_type() { return true; } \
};
diff --git a/src/event_queue.hpp b/src/event_queue.hpp
index 9250eb862e5a167fdf1ae0b962e0e30e3e772126..7f3034f3f08bbb2235bf457c95e8eedfe8d0895d 100644
--- a/src/event_queue.hpp
+++ b/src/event_queue.hpp
@@ -13,26 +13,26 @@ struct postsynaptic_spike_event {
uint32_t target;
float time;
float weight;
-
- float event_time() const { return time; }
};
+inline float event_time(const postsynaptic_spike_event &ev) { return ev.time; }
+
struct sample_event {
uint32_t sampler_index;
float time;
-
- float event_time() const { return 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 <type Event>
+template <typename Event>
class event_queue {
public :
using value_type = Event;
- using time_type = template std::result_of<decltype(&Event::event_time)(Event)>::type;
+ using time_type = decltype(event_time(std::declval<Event>()));
// create
event_queue() {}
@@ -46,7 +46,7 @@ public :
}
// push thing
- void push(local_event e) {
+ void push(const value_type &e) {
queue_.push(e);
}
@@ -68,8 +68,8 @@ public :
private:
struct event_greater {
- bool operator(const Event &a, const Event &b) {
- return a.event_time() > b.event_time();
+ bool operator()(const Event &a, const Event &b) {
+ return event_time(a) > event_time(b);
}
};
@@ -84,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 4bd215716b2b77471e8c488925f3c54f6105e165..5cadd2912d79f793fa01812fc4dcb2a24aefee3d 100644
--- a/src/fvm_cell.hpp
+++ b/src/fvm_cell.hpp
@@ -27,7 +27,7 @@ namespace fvm {
template <typename T, typename I>
class fvm_cell {
- public :
+public:
fvm_cell() = default;
@@ -117,7 +117,7 @@ class fvm_cell {
void advance_to(value_type tfinal, 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);
}
@@ -136,7 +136,14 @@ class fvm_cell {
value_type time() const { return t_; }
- private:
+ value_type probe(uint32_t i) const {
+ auto p = probes_[i];
+ return (*p.first)[p.second];
+ }
+
+ std::size_t num_probes() const { return probes_.size(); }
+
+private:
/// current time
value_type t_ = value_type{0};
@@ -180,8 +187,12 @@ class fvm_cell {
std::vector<std::pair<uint32_t, i_clamp>> stimulii_;
- /// event queue
- event_queue events_;
+ std::vector<std::pair<const vector_type*, uint32_t>> probes_;
+
+ /*
+ /// spike event queue
+ event_queue<postsynaptic_spike_event> events_;
+ */
};
////////////////////////////////////////////////////////////////////////////////
@@ -398,6 +409,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({&voltage_, comp});
+ break;
+ case nest::mc::cell::membrane_current:
+ probes_.push_back({¤t_, 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);