diff --git a/miniapp/io.cpp b/miniapp/io.cpp index e1a059ae1c072275d43211c84550d7ccb6b5c74e..5b7a1cb338999c67a6885287fdccb02195322bf6 100644 --- a/miniapp/io.cpp +++ b/miniapp/io.cpp @@ -8,7 +8,7 @@ namespace io { // for now this is just a placeholder options read_options(std::string fname) { // 10 cells, 1 synapses per cell, 10 compartments per segment - return {100, 1, 100}; + return {20, 1, 100}; } std::ostream& operator<<(std::ostream& o, const options& opt) { diff --git a/miniapp/miniapp.cpp b/miniapp/miniapp.cpp index 5f2082e1b41945ecd0452a27984b6dc34bcb20b8..fdac119896bdeadf9b9d40fffd05bb3a6345bb5b 100644 --- a/miniapp/miniapp.cpp +++ b/miniapp/miniapp.cpp @@ -101,6 +101,19 @@ 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; + + // sampling // WIP COME BACK HERE + struct simple_sampler { + index_type group_gid; + std::string name; + double dt; + + std::vector<std::pair<float,double>> tv; + util::optional<float> operator()(float t, double v) { + define + } + + }; // define some global model parameters @@ -269,6 +282,25 @@ void all_to_all_model(nest::mc::io::options& opt, model& m) { auto start_network = timer.tic(); m.init_communicator(); + // 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; + } + + 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; + + // 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) } + }; + } + // lid is local cell/group id for (auto lid=0u; lid<ncell_local; ++lid) { auto target = m.communicator.target_gid_from_group_lid(lid); @@ -348,6 +380,13 @@ 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(nest::mc::cell::membrane_potential, {0,0}); + auto probe_dendrite = cell.add_probe(nest::mc::cell::membrane_potential, {1,0.5}); + + EXPECT(probe_soma==0); + EXPECT(probe_dendrite==1); + return cell; } diff --git a/src/cell_group.hpp b/src/cell_group.hpp index 85a0832bf24a5720951b35db392133e2d877a312..f1360ffe03c6fd344191955ed7b9eec19b2d9b3f 100644 --- a/src/cell_group.hpp +++ b/src/cell_group.hpp @@ -10,10 +10,20 @@ namespace nest { namespace mc { +// samplers take a time and sample value, and return an optional time +// for the next desired sample. + +struct sampler { + 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; @@ -50,6 +60,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]; @@ -63,6 +85,20 @@ class cell_group { 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); + + 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); + } + } + #ifdef SPLAT tt.push_back(cell_.time()); vs.push_back(cell_.voltage({0,0.0})); @@ -122,7 +158,9 @@ class cell_group { spikes_.clear(); } - private : + std::vector<sampler> samplers; + +private: #ifdef SPLAT // REMOVE as soon as we have a better way to probe cell state @@ -141,10 +179,16 @@ 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_; }; } // namespace mc diff --git a/src/communication/communicator.hpp b/src/communication/communicator.hpp index 517a2f6c3884c66df22914776849fcfd2da7e136..eee922541e41d3d3d8b317660ab5b390b42d99fa 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())) { diff --git a/src/event_queue.hpp b/src/event_queue.hpp index 989e34873a1f3ff851d0a3a9f85276cea9d9bf7a..9250eb862e5a167fdf1ae0b962e0e30e3e772126 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; + + float event_time() const { return 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 bool operator > (local_event const& l, local_event const& r) { - return l.time > r.time; -} + float event_time() const { return 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> class event_queue { public : + using value_type = Event; + using time_type = template std::result_of<decltype(&Event::event_time)(Event)>::type; + // create event_queue() {} @@ -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 a.event_time() > b.event_time(); + } + }; + std::priority_queue< - local_event, - std::vector<local_event>, - std::greater<local_event> + Event, + std::vector<Event>, + event_greater > queue_; };