diff --git a/miniapp/io.cpp b/miniapp/io.cpp index 7a6e9c67ae4fee503a0230baeab2e18dbcafaf79..3dd5e96dfc3432581fc26103275644a6ca1fb86e 100644 --- a/miniapp/io.cpp +++ b/miniapp/io.cpp @@ -111,46 +111,13 @@ static void update_option(util::optional<T>& opt, const nlohmann::json& j, const // Read options from (optional) json file and command line arguments. cl_options read_options(int argc, char** argv, bool allow_write) { - - // Default options: - const cl_options defopts{ - 1000, // number of cells - 500, // synapses_per_cell - "expsyn", // synapse type - 100, // compartments_per_segment - 100., // tfinal - 0.025, // dt - false, // all_to_all - false, // ring - 1, // group_size - false, // probe_soma_only - 0.0, // probe_ratio - "trace_", // trace_prefix - util::nothing, // trace_max_gid - util::nothing, // morphologies - false, // morph_rr; - false, // report_compartments; - - // spike_output_parameters: - false, // spike output - false, // single_file_per_simulation - true, // Overwrite outputfile if exists - "./", // output path - "spikes", // file name - "gdf", // file extension - - // dry run parameters: - 1, // default dry run size - - // Turn on/off profiling output for all ranks - false - }; - cl_options options; std::string save_file = ""; // Parse command line arguments. try { + cl_options defopts; + CustomCmdLine cmd("nest mc miniapp harness", "0.1"); TCLAP::ValueArg<std::string> ifile_arg( @@ -179,6 +146,11 @@ cl_options read_options(int argc, char** argv, bool allow_write) { TCLAP::ValueArg<double> dt_arg( "d", "dt", "set simulation time step to <time> ms", false, defopts.dt, "time", cmd); + TCLAP::ValueArg<double> bin_dt_arg( + "", "bin-dt", "set event binning interval to <time> ms", + false, defopts.bin_dt, "time", cmd); + TCLAP::SwitchArg bin_regular_arg( + "","bin-regular","use 'regular' binning policy instead of 'following'", cmd, false); TCLAP::SwitchArg all_to_all_arg( "m","alltoall","all to all network", cmd, false); TCLAP::SwitchArg ring_arg( @@ -215,8 +187,6 @@ cl_options read_options(int argc, char** argv, bool allow_write) { cmd.reorder_arguments(); cmd.parse(argc, argv); - options = defopts; - std::string ifile_name = ifile_arg.getValue(); if (ifile_name != "") { // Read parameters from specified JSON file first, to allow @@ -232,6 +202,8 @@ cl_options read_options(int argc, char** argv, bool allow_write) { update_option(options.syn_type, fopts, "syn_type"); update_option(options.compartments_per_segment, fopts, "compartments"); update_option(options.dt, fopts, "dt"); + update_option(options.bin_dt, fopts, "bin_dt"); + update_option(options.bin_regular, fopts, "bin_regular"); update_option(options.tfinal, fopts, "tfinal"); update_option(options.all_to_all, fopts, "all_to_all"); update_option(options.ring, fopts, "ring"); @@ -275,6 +247,8 @@ cl_options read_options(int argc, char** argv, bool allow_write) { update_option(options.compartments_per_segment, ncompartments_arg); update_option(options.tfinal, tfinal_arg); update_option(options.dt, dt_arg); + update_option(options.bin_dt, bin_dt_arg); + update_option(options.bin_regular, bin_regular_arg); update_option(options.all_to_all, all_to_all_arg); update_option(options.ring, ring_arg); update_option(options.group_size, group_size_arg); @@ -315,6 +289,8 @@ cl_options read_options(int argc, char** argv, bool allow_write) { fopts["syn_type"] = options.syn_type; fopts["compartments"] = options.compartments_per_segment; fopts["dt"] = options.dt; + fopts["bin_dt"] = options.bin_dt; + fopts["bin_regular"] = options.bin_regular; fopts["tfinal"] = options.tfinal; fopts["all_to_all"] = options.all_to_all; fopts["ring"] = options.ring; @@ -358,6 +334,9 @@ std::ostream& operator<<(std::ostream& o, const cl_options& options) { o << " synapses/cell : " << options.synapses_per_cell << "\n"; o << " simulation time : " << options.tfinal << "\n"; o << " dt : " << options.dt << "\n"; + o << " binning dt : " << options.bin_dt << "\n"; + o << " binning policy : " << + (options.bin_dt==0? "none": options.bin_regular? "regular": "following") << "\n"; o << " all to all network : " << (options.all_to_all ? "yes" : "no") << "\n"; o << " ring network : " << (options.ring ? "yes" : "no") << "\n"; o << " group size : " << options.group_size << "\n"; diff --git a/miniapp/io.hpp b/miniapp/io.hpp index cd3eca6676d6d9e4d4f7a1b14ffe2167df4dadac..a005cd290daf065fd4fc52121a04b08699c3e3e4 100644 --- a/miniapp/io.hpp +++ b/miniapp/io.hpp @@ -12,38 +12,51 @@ namespace nest { namespace mc { namespace io { -// holds the options for a simulation run +// Holds the options for a simulation run. +// Default constructor gives default options. + struct cl_options { - uint32_t cells; - uint32_t synapses_per_cell; - std::string syn_type; - uint32_t compartments_per_segment; - double tfinal; - double dt; - bool all_to_all; - bool ring; - uint32_t group_size; - bool probe_soma_only; - double probe_ratio; - std::string trace_prefix; - util::optional<unsigned> trace_max_gid; + // Cell parameters: + uint32_t cells = 1000; + uint32_t synapses_per_cell = 500; + std::string syn_type = "expsyn"; + uint32_t compartments_per_segment = 100; util::optional<std::string> morphologies; - bool morph_rr; - bool report_compartments; - - // Parameters for spike output - bool spike_file_output; - bool single_file_per_rank; - bool over_write; - std::string output_path; - std::string file_name; - std::string file_extension; - - // dry run parameters - int dry_run_ranks; - - // Turn on/off profiling output for all ranks - bool profile_only_zero; + bool morph_rr = false; // False => pick morphologies randomly, true => pick morphologies round-robin. + + // Network type (default is rgraph): + bool all_to_all = false; + bool ring = false; + + // Simulation running parameters: + double tfinal = 100.; + double dt = 0.025; + uint32_t group_size = 1; + bool bin_regular = false; // False => use 'following' instead of 'regular'. + double bin_dt = 0.0025; // 0 => no binning. + + // Probe/sampling specification. + bool probe_soma_only = false; + double probe_ratio = 0; // Proportion of cells to probe. + std::string trace_prefix = "trace_"; + util::optional<unsigned> trace_max_gid; // Only make traces up to this gid. + + // Parameters for spike output. + bool spike_file_output = false; + bool single_file_per_rank = false; + bool over_write = true; + std::string output_path = "./"; + std::string file_name = "spikes"; + std::string file_extension = "gdf"; + + // Dry run parameters (pertinent only when built with 'dryrun' distrib model). + int dry_run_ranks = 1; + + // Turn on/off profiling output for all ranks. + bool profile_only_zero = false; + + // Report (inefficiently) on number of cell compartments in sim. + bool report_compartments = false; }; class usage_error: public std::runtime_error { diff --git a/miniapp/miniapp.cpp b/miniapp/miniapp.cpp index c4add5583f21a69707998d6cf1a6591a78b58d39..2a6ef9b740a2f48921529de2fa90daaf403ce002 100644 --- a/miniapp/miniapp.cpp +++ b/miniapp/miniapp.cpp @@ -101,7 +101,15 @@ int main(int argc, char** argv) { report_compartment_stats(*recipe); } - // inject some artificial spikes, 1 per 20 neurons. + // Specify event binning/coalescing. + auto binning_policy = + options.bin_dt==0? binning_kind::none: + options.bin_regular? binning_kind::regular: + binning_kind::following; + + m.set_binning_policy(binning_policy, options.bin_dt); + + // Inject some artificial spikes, 1 per 20 neurons. std::vector<cell_gid_type> local_sources; cell_gid_type first_spike_cell = 20*((cell_range.first+19)/20); for (auto c=first_spike_cell; c<cell_range.second; c+=20) { @@ -109,7 +117,7 @@ int main(int argc, char** argv) { m.add_artificial_spike({c, 0}); } - // attach samplers to all probes + // Attach samplers to all probes std::vector<std::unique_ptr<sample_trace_type>> traces; const model_type::time_type sample_dt = 0.1; for (auto probe: m.probes()) { diff --git a/src/cell_group.hpp b/src/cell_group.hpp index 55cf7dd9d4bb9f50b76e49fdc6f99bdd0d867681..50c696c3c157c2523a946c058da3649ae919ad8c 100644 --- a/src/cell_group.hpp +++ b/src/cell_group.hpp @@ -143,51 +143,61 @@ public: } void advance(time_type tfinal, time_type dt) { - while (lowered_.time()<tfinal) { - // take any pending samples - time_type cell_time = lowered_.time(); + EXPECTS(lowered_.state_synchronized()); + + // Bin pending events and enqueue on lowered state. + time_type ev_min_time = lowered_.max_time(); // (but we're synchronized here) + while (auto ev = events_.pop_if_before(tfinal)) { + auto handle = get_target_handle(ev->target); + auto binned_ev_time = binner_.bin(ev->target.gid, ev->time, ev_min_time); + lowered_.add_event(binned_ev_time, handle, ev->weight); + } + + lowered_.setup_integration(tfinal, dt); + + std::vector<sample_event<time_type>> requeue_sample_events; + while (!lowered_.integration_complete()) { + // Take any pending samples. + // TODO: Placeholder: this will be replaced by a backend polling implementation. PE("sampling"); - while (auto m = sample_events_.pop_if_before(cell_time)) { + time_type cell_max_time = lowered_.max_time(); + + requeue_sample_events.clear(); + while (auto m = sample_events_.pop_if_before(cell_max_time)) { auto& s = samplers_[m->sampler_index]; EXPECTS((bool)s.sampler); - auto next = s.sampler(lowered_.time(), lowered_.probe(s.handle)); - if (next) { - m->time = std::max(*next, cell_time); - sample_events_.push(*m); + time_type cell_time = lowered_.time(s.cell_gid-gid_base_); + 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); + } } } + for (auto& ev: requeue_sample_events) { + sample_events_.push(std::move(ev)); + } PL(); - // look for events in the next time step - time_type tstep = lowered_.time()+dt; - tstep = std::min(tstep, tfinal); - auto next = events_.pop_if_before(tstep); + // Ask lowered_ cell to integrate 'one step', delivering any + // events accordingly. + // TODO: Placeholder: with backend polling for samplers, we will + // request that the lowered cell perform the integration all the + // way to tfinal. - // apply events that are due within the smallest allowed time step. - while (next && (next->time-lowered_.time()) < 0.1*(dt)) { - auto handle = get_target_handle(next->target); - lowered_.deliver_event(handle, next->weight); - next = events_.pop_if_before(tstep); - } - - // integrate cell state - time_type tnext = next ? next->time: tstep; - lowered_.advance(tnext - lowered_.time()); + lowered_.step_integration(); if (util::is_debug_mode() && !lowered_.is_physical_solution()) { std::cerr << "warning: solution out of bounds for cell " - << gid_base_ << " at t " << lowered_.time() << " ms\n"; + << gid_base_ << " at (max) t " << lowered_.max_time() << " ms\n"; } - - // apply events - PE("events"); - if (next) { - auto handle = get_target_handle(next->target); - lowered_.deliver_event(handle, next->weight); - } - PL(); } // Copy out spike voltage threshold crossings from the back end, then @@ -206,7 +216,7 @@ public: template <typename R> void enqueue_events(const R& events) { - for (auto e : events) { + for (auto& e: events) { events_.push(e); } } @@ -231,7 +241,7 @@ public: auto handle = get_probe_handle(probe_id); auto sampler_index = uint32_t(samplers_.size()); - samplers_.push_back({handle, s}); + samplers_.push_back({handle, probe_id.gid, s}); sampler_start_times_.push_back(start_time); sample_events_.push({sampler_index, start_time}); } @@ -286,6 +296,7 @@ private: struct sampler_entry { typename lowered_cell_type::probe_handle handle; + cell_gid_type cell_gid; sampler_function sampler; }; diff --git a/src/event_queue.hpp b/src/event_queue.hpp index 2a4c7ec7a21ce58460bb8d121c4b0c50df52789b..d6573e0fee63c0cf84d3f3d52ef06b996e4850d3 100644 --- a/src/event_queue.hpp +++ b/src/event_queue.hpp @@ -4,6 +4,7 @@ #include <ostream> #include <queue> #include <type_traits> +#include <utility> #include "common_types.hpp" #include "util/meta.hpp" @@ -64,6 +65,10 @@ public : queue_.push(e); } + void push(value_type&& e) { + queue_.push(std::move(e)); + } + bool empty() const { return size()==0; } @@ -72,10 +77,23 @@ public : return queue_.size(); } - // Pop and return top event `ev` of queue if `t_until` > `event_time(ev)`. - util::optional<value_type> pop_if_before(const time_type& t_until) { + // Return time t of head of queue if `t_until` > `t`. + util::optional<time_type> time_if_before(const time_type& t_until) { + if (queue_.empty()) { + return util::nothing; + } + + using ::nest::mc::event_time; + auto t = event_time(queue_.top()); + return t_until > t? util::just(t): util::nothing; + } + + // Generic conditional pop: pop and return head of queue if + // queue non-empty and the head satisfies predicate. + template <typename Pred> + util::optional<value_type> pop_if(Pred&& pred) { using ::nest::mc::event_time; - if (!queue_.empty() && t_until > event_time(queue_.top())) { + if (!queue_.empty() && pred(queue_.top())) { auto ev = queue_.top(); queue_.pop(); return ev; @@ -85,6 +103,22 @@ public : } } + // Pop and return top event `ev` of queue if `t_until` > `event_time(ev)`. + util::optional<value_type> pop_if_before(const time_type& t_until) { + using ::nest::mc::event_time; + return pop_if( + [&t_until](const value_type& ev) { return t_until > event_time(ev); } + ); + } + + // Pop and return top event `ev` of queue unless `event_time(ev)` > `t_until` + util::optional<value_type> pop_if_not_after(const time_type& t_until) { + using ::nest::mc::event_time; + return pop_if( + [&t_until](const value_type& ev) { return !(event_time(ev) > t_until); } + ); + } + // Clear queue and free storage. void clear() { queue_ = decltype(queue_){}; diff --git a/src/fvm_multicell.hpp b/src/fvm_multicell.hpp index 5cc34aa8f9d01087fb7b0e1477ba5204fd819bdb..78cb84288d5c8fb0ab28a811404dfef686465d69 100644 --- a/src/fvm_multicell.hpp +++ b/src/fvm_multicell.hpp @@ -85,8 +85,56 @@ public: return (this->*h.first)[h.second]; } - /// integrate all cell state forward in time - void advance(double dt); + /// Initialize state prior to a sequence of integration steps. + void setup_integration(value_type tfinal, value_type dt_max) { + EXPECTS(tfinal>t_); + EXPECTS(dt_max>0); + + tfinal_ = tfinal; + dt_max_ = dt_max; + integration_running_ = true; + + // TODO: Placeholder; construct backend event delivery + // data structure from events added. + EXPECTS(events_.empty()); + for (auto& ev: staged_events_) { + EXPECTS(ev.time<tfinal); + events_.push(std::move(ev)); + } + + staged_events_.clear(); + } + + /// Advance one integration step, up to `dt_max_` in each cell. + void step_integration(); + + /// Query integration completion state. + bool integration_complete() const { + return !integration_running_; + } + + /// Query per-cell time state. (Placeholders) + value_type time(size_type cell_index) const { + return t_; + } + + value_type max_time() const { + return t_; + } + + bool state_synchronized() const { + return true; + } + + /// Add an event for processing in next integration stage. + void add_event(value_type ev_time, target_handle h, value_type weight) { + EXPECTS(!integration_running_); + + // TODO: Placeholder; either add to backend event list structure + // incrementally here, or store by cell gid offset for construction + // all at once in `start_integration()`. + staged_events_.push_back({ev_time, h, weight}); + } /// Following types and methods are public only for testing: @@ -167,8 +215,6 @@ public: return it==mechanisms_.end() ? util::nothing: util::just(*it); } - value_type time() const { return t_; } - std::size_t num_probes() const { return probes_.size(); } // @@ -194,15 +240,38 @@ public: } private: - threshold_watcher threshold_watcher_; /// current time [ms] - value_type t_ = value_type{0}; + value_type t_ = 0; /// resting potential (initial voltage condition) value_type resting_potential_ = -65; + /// final time in integration round [ms] + value_type tfinal_ = 0; + + /// max time step for integration [ms] + value_type dt_max_ = 0; + + /// flag: true after a call to `setup_integration()`; reset + /// once integration to `tfinal_` is complete. + bool integration_running_ = false; + + struct deliverable_event { + double time; + target_handle handle; + double weight; + }; + + /// events staged for upcoming integration stage + std::vector<deliverable_event> staged_events_; + + /// event queue for integration period + /// TODO: Placeholder; this will change when event lists are moved to + /// a backend structure. + event_queue<deliverable_event> events_; + /// the linear system for implicit time stepping of cell state matrix_type matrix_; @@ -753,10 +822,35 @@ void fvm_multicell<Backend>::reset() { // NOTE: this has to come after the voltage_ values have been reinitialized, // because these values are used by the watchers to set their initial state. threshold_watcher_.reset(t_); + + // Reset integration state. + tfinal_ = 0; + dt_max_ = 0; + integration_running_ = false; + staged_events_.clear(); + events_.clear(); } + template <typename Backend> -void fvm_multicell<Backend>::advance(double dt) { +void fvm_multicell<Backend>::step_integration() { + // Integrate cell states from `t_` by `dt_max_` if permissible, + // or otherwise until the next event time or `t_final`. + EXPECTS(integration_running_); + + while (auto ev = events_.pop_if_not_after(t_)) { + deliver_event(ev->handle, ev->weight); + } + + value_type t_to = std::min(t_+dt_max_, tfinal_); + + if (auto t = events_.time_if_before(t_to)) { + t_to = *t; + } + + value_type dt = t_to-t_; + EXPECTS(dt>0); + PE("current"); memory::fill(current_, 0.); @@ -788,10 +882,15 @@ void fvm_multicell<Backend>::advance(double dt) { } PL(); - t_ += dt; + t_ = t_to; // update spike detector thresholds threshold_watcher_.test(t_); + + // are we there yet? + integration_running_ = t_<tfinal_; + + EXPECTS(integration_running_ || events_.empty()); } } // namespace fvm diff --git a/src/model.hpp b/src/model.hpp index c411b4a8092551aa974b052d3235a22094f33fab..ee0a12eb00ce18b2ada44459a5897c731188fe49 100644 --- a/src/model.hpp +++ b/src/model.hpp @@ -231,6 +231,13 @@ public: return bounds.second-bounds.first; } + // Set event binning policy on all our groups. + void set_binning_policy(binning_kind policy, time_type bin_interval) { + for (auto& group: cell_groups_) { + group.set_binning_policy(policy, bin_interval); + } + } + // access cell_group directly cell_group_type& group(int i) { return cell_groups_[i]; diff --git a/tests/unit/test_cell_group.cpp b/tests/unit/test_cell_group.cpp index afd2d20354b22f49f461a6427d50a64805bd02f0..af5382d740c05c0b48fb73fd74e78b2cf286e1d1 100644 --- a/tests/unit/test_cell_group.cpp +++ b/tests/unit/test_cell_group.cpp @@ -23,8 +23,7 @@ nest::mc::cell make_cell() { } TEST(cell_group, test) { - using cell_group_type = cell_group<fvm_cell>; - auto group = cell_group_type{0, util::singleton_view(make_cell())}; + cell_group<fvm_cell> group{0, util::singleton_view(make_cell())}; group.advance(50, 0.01); diff --git a/tests/unit/test_event_queue.cpp b/tests/unit/test_event_queue.cpp index 6a6752b651bbd1538682a7f48f1535d9a06e1757..60a8da476221f7595657be304fc6c8e67f6fb111 100644 --- a/tests/unit/test_event_queue.cpp +++ b/tests/unit/test_event_queue.cpp @@ -5,9 +5,9 @@ #include <event_queue.hpp> -TEST(event_queue, push) -{ - using namespace nest::mc; +using namespace nest::mc; + +TEST(event_queue, push) { using ps_event_queue = event_queue<postsynaptic_spike_event<float>>; ps_event_queue q; @@ -28,9 +28,7 @@ TEST(event_queue, push) EXPECT_TRUE(std::is_sorted(times.begin(), times.end())); } -TEST(event_queue, pop_if_before) -{ - using namespace nest::mc; +TEST(event_queue, pop_if_before) { using ps_event_queue = event_queue<postsynaptic_spike_event<float>>; cell_member_type target[4] = { @@ -86,6 +84,34 @@ TEST(event_queue, pop_if_before) EXPECT_FALSE(e6); } +TEST(event_queue, pop_if_not_after) { + struct event { + int time; + + event(int t): time(t) {} + }; + + event_queue<event> queue; + + queue.push(1); + queue.push(3); + queue.push(5); + + auto e1 = queue.pop_if_not_after(2); + EXPECT_TRUE(e1); + EXPECT_EQ(1, e1->time); + + auto e2 = queue.pop_if_before(3); + EXPECT_FALSE(e2); + + auto e3 = queue.pop_if_not_after(3); + EXPECT_TRUE(e3); + EXPECT_EQ(3, e3->time); + + auto e4 = queue.pop_if_not_after(4); + EXPECT_FALSE(e4); +} + // Event queues can be defined for arbitrary copy-constructible events // for which `event_time(ev)` returns the corresponding time. Time values just // need to be well-ordered on '>'. @@ -105,10 +131,7 @@ struct minimal_event { const wrapped_float& event_time(const minimal_event& ev) { return ev.value; } -TEST(event_queue, minimal_event_impl) -{ - using nest::mc::event_queue; - +TEST(event_queue, minimal_event_impl) { minimal_event events[] = { minimal_event(3.f), minimal_event(2.f), diff --git a/tests/unit/test_fvm_multi.cpp b/tests/unit/test_fvm_multi.cpp index 5b678bb48c347620157be8b3ee80f32287bf4840..66226d9e9464df8fe014a6ff12c41d8817331bbc 100644 --- a/tests/unit/test_fvm_multi.cpp +++ b/tests/unit/test_fvm_multi.cpp @@ -77,7 +77,9 @@ TEST(fvm_multi, init) // test that the matrix is initialized with sensible values //J.build_matrix(0.01); - fvcell.advance(0.01); + fvcell.setup_integration(0.01, 0.01); + fvcell.step_integration(); + auto& mat = J.state_; auto test_nan = [](decltype(mat.u) v) { for(auto val : v) if(val != val) return false; diff --git a/tests/unit/test_probe.cpp b/tests/unit/test_probe.cpp index 8027fc2fef77e3dfb194a9312889452066db534c..005eead2b61f2e4eb526edbe96ed2bddd47d611a 100644 --- a/tests/unit/test_probe.cpp +++ b/tests/unit/test_probe.cpp @@ -73,7 +73,8 @@ TEST(probe, fvm_multicell) EXPECT_EQ(lcell.voltage()[probes[1].second], lcell.probe(probes[1])); EXPECT_EQ(lcell.current()[probes[2].second], lcell.probe(probes[2])); - lcell.advance(0.05); + 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]));