WIP

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) {

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;
 }
 

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

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())) {

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_;
 };