diff --git a/src/cell_group.hpp b/src/cell_group.hpp
index e70d4b0bb5826b6ccde449f63a7fad6dab29c783..b0bbccbcb46757afe72695413cbc525b0d86eb9a 100644
--- a/src/cell_group.hpp
+++ b/src/cell_group.hpp
@@ -9,20 +9,21 @@
#include <event_queue.hpp>
#include <spike.hpp>
#include <spike_source.hpp>
+#include <util/singleton.hpp>
#include <profiling/profiler.hpp>
namespace nest {
namespace mc {
-template <typename Cell>
+template <typename LoweredCell>
class cell_group {
public:
using index_type = cell_gid_type;
- using cell_type = Cell;
- using value_type = typename cell_type::value_type;
- using size_type = typename cell_type::value_type;
- using spike_detector_type = spike_detector<Cell>;
+ using lowered_cell_type = LoweredCell;
+ using value_type = typename lowered_cell_type::value_type;
+ using size_type = typename lowered_cell_type::value_type;
+ using spike_detector_type = spike_detector<lowered_cell_type>;
using source_id_type = cell_member_type;
using time_type = float;
@@ -36,9 +37,13 @@ public:
cell_group() = default;
cell_group(cell_gid_type gid, const cell& c) :
- gid_base_{gid}, cell_{c}
+ gid_base_{gid}
{
- initialize_cells();
+ detector_handles_.resize(c.detectors().size());
+ target_handles_.resize(c.synapses().size());
+ probe_handles_.resize(c.probes().size());
+
+ cell_.initialize(util::singleton_view(c), detector_handles_, target_handles_, probe_handles_);
// Create spike detectors and associate them with globally unique source ids,
// as specified by cell gid and cell-local zero-based index.
@@ -46,18 +51,19 @@ public:
cell_gid_type source_gid = gid_base_;
cell_lid_type source_lid = 0u;
+ unsigned i = 0;
for (auto& d : c.detectors()) {
cell_member_type source_id{source_gid, source_lid++};
spike_sources_.push_back({
- source_id, spike_detector_type(cell_, d.location, d.threshold, 0.f)
+ source_id, spike_detector_type(cell_, detector_handles_[i], d.threshold, 0.f)
});
}
}
void reset() {
remove_samplers();
- initialize_cells();
+ cell_.reset();
for (auto& spike_source: spike_sources_) {
spike_source.source.reset(cell_, 0.f);
}
@@ -70,11 +76,10 @@ public:
PE("sampling");
while (auto m = sample_events_.pop_if_before(cell_time)) {
- auto& sampler_spec = samplers_[m->sampler_index];
- EXPECTS((bool)sampler_spec.sampler);
+ auto& s = samplers_[m->sampler_index];
+ EXPECTS((bool)s.sampler);
+ auto next = s.sampler(cell_.time(), cell_.probe(s.handle));
- index_type probe_index = sampler_spec.probe_id.index;
- auto next = sampler_spec.sampler(cell_.time(), cell_.probe(probe_index));
if (next) {
m->time = std::max(*next, cell_time);
sample_events_.push(*m);
@@ -105,12 +110,14 @@ public:
// apply events
if (next) {
- cell_.apply_event(next.get());
+ auto handle = target_handles_[next->target.index];
+ cell_.deliver_event(handle, next->weight);
// apply events that are due within some epsilon of the current
// time step. This should be a parameter. e.g. with for variable
// order time stepping, use the minimum possible time step size.
while(auto e = events_.pop_if_before(cell_.time()+dt/10.)) {
- cell_.apply_event(e.get());
+ auto handle = target_handles_[e->target.index];
+ cell_.deliver_event(handle, e->weight);
}
}
PL();
@@ -128,9 +135,6 @@ public:
const std::vector<spike<source_id_type, time_type>>&
spikes() const { return spikes_; }
- cell_type& cell() { return cell_; }
- const cell_type& cell() const { return cell_; }
-
const std::vector<spike_source_type>&
spike_sources() const {
return spike_sources_;
@@ -141,8 +145,9 @@ public:
}
void add_sampler(cell_member_type probe_id, sampler_function s, time_type start_time = 0) {
+ EXPECTS(probe_id.gid==gid_base_);
auto sampler_index = uint32_t(samplers_.size());
- samplers_.push_back({probe_id, s});
+ samplers_.push_back({probe_handles_[probe_id.index], s});
sample_events_.push({sampler_index, start_time});
}
@@ -151,22 +156,21 @@ public:
samplers_.clear();
}
-private:
- void initialize_cells() {
- cell_.voltage()(memory::all) = -65.;
- cell_.initialize();
+ value_type probe(cell_member_type probe_id) const {
+ return cell_.probe(probe_handles_[probe_id.index]);
}
+private:
/// gid of first cell in group
cell_gid_type gid_base_;
/// the lowered cell state (e.g. FVM) of the cell
- cell_type cell_;
+ lowered_cell_type cell_;
/// spike detectors attached to the cell
std::vector<spike_source_type> spike_sources_;
- //. spikes that are generated
+ /// spikes that are generated
std::vector<spike<source_id_type, time_type>> spikes_;
/// pending events to be delivered
@@ -178,11 +182,18 @@ private:
/// 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_;
+ /// handles for accessing lowered cell
+ using detector_handle = typename lowered_cell_type::detector_handle;
+ std::vector<detector_handle> detector_handles_;
+
+ using target_handle = typename lowered_cell_type::target_handle;
+ std::vector<target_handle> target_handles_;
+
+ using probe_handle = typename lowered_cell_type::probe_handle;
+ std::vector<probe_handle> probe_handles_;
struct sampler_entry {
- cell_member_type probe_id;
+ typename lowered_cell_type::probe_handle handle;
sampler_function sampler;
};
diff --git a/src/fvm_cell.hpp b/src/fvm_cell.hpp
index 526fd22a391b477d5af3c88a4535b4721c83cbfe..5335c8129a5963d8d288094116442c431279a2f7 100644
--- a/src/fvm_cell.hpp
+++ b/src/fvm_cell.hpp
@@ -1,6 +1,7 @@
#pragma once
#include <algorithm>
+#include <iterator>
#include <map>
#include <set>
#include <string>
@@ -14,27 +15,155 @@
#include <matrix.hpp>
#include <mechanism.hpp>
#include <mechanism_catalogue.hpp>
+#include <profiling/profiler.hpp>
#include <segment.hpp>
#include <stimulus.hpp>
#include <util.hpp>
-#include <profiling/profiler.hpp>
+#include <util/meta.hpp>
#include <vector/include/Vector.hpp>
+/*
+ * Lowered cell implementation based on finite volume method.
+ *
+ * TODO: Move following description of internal API
+ * to a better location or document.
+ *
+ * Lowered cells are managed by the `cell_group` class.
+ * A `cell_group` manages one lowered cell instance, which
+ * may in turn simulate one or more cells.
+ *
+ * The outward facing interface for `cell_group` is defined
+ * in terms of cell gids and member indices; the interface
+ * between `cell_group` and a lowered cell is described below.
+ *
+ * The motivation for the following interface is to
+ * 1. Provide more flexibility in lowered cell implementation
+ * 2. Save memory and duplicated index maps by associating
+ * externally visible objects with opaque handles.
+ *
+ * In the following, `lowered_cell` represents the lowered
+ * cell class, and `lowered` an instance thereof.
+ *
+ * `lowered_cell::detector_handle`
+ * Type of handles used by spike detectors to query
+ * membrane voltage.
+ *
+ * `lowered_cell::probe_handle`
+ * Type of handle used to query the value of a probe.
+ *
+ * `lowered_cell::target_handle`
+ * Type of handle used to identify the target of a
+ * postsynaptic spike event.
+ *
+ * `lowered_cell::value_type`
+ * Floating point type used for internal states
+ * and simulation time.
+ *
+ * `lowered_cell()`
+ * Default constructor; performs no cell-specific
+ * initialization.
+ *
+ * `lowered.initialize(const Cells& cells, ...)`
+ * Allocate and initalize data structures to simulate
+ * the cells described in the collection `cells`,
+ * where each item is of type `nest::mc::cell`.
+ *
+ * Remaining arguments consist of references to
+ * collections for storing:
+ * 1. Detector handles
+ * 2. Target handles
+ * 3. Probe handles
+ *
+ * Handles are written in the same order as they
+ * appear in the provided cell descriptions.
+ *
+ * `lowered.reset()`
+ *
+ * Resets state to initial conditiions and sets
+ * internal simulation time to 0.
+ *
+ * `lowered.advance(value_type dt)`
+ *
+ * Advanece simulation state by `dt` (value in
+ * milliseconds). For `fvm_cell` at least,
+ * this corresponds to one integration step.
+ *
+ * `lowered.deliver_event(target_handle target, value_type weight)`
+ *
+ * Update target-specifc state based on the arrival
+ * of a postsynaptic spike event with given weight.
+ *
+ * `lowered.detect_voltage(detector_handle)`
+ *
+ * Return membrane voltage at detector site as specified
+ * by handle.
+ *
+ * `lowered.probe(probe_handle)`
+ *
+ * Return value of corresponding probe.
+ *
+ * `lowered.resting_potential(value_type potential)`
+ *
+ * Set the steady-state membrane voltage used for the
+ * cell initial condition. (Defaults to -65 mV currently.)
+ */
+
namespace nest {
namespace mc {
namespace fvm {
-template <typename T, typename I>
+template <typename Value, typename Index>
class fvm_cell {
public:
-
fvm_cell() = default;
/// the real number type
- using value_type = T;
+ using value_type = Value;
+
/// the integral index type
- using size_type = I;
+ using size_type = Index;
+
+ /// the container used for indexes
+ using index_type = memory::HostVector<size_type>;
+
+ /// the container used for values
+ using vector_type = memory::HostVector<value_type>;
+
+ /// API for cell_group (see above):
+
+ using detector_handle = size_type;
+ using target_handle = std::pair<size_type, size_type>;
+ using probe_handle = std::pair<const vector_type fvm_cell::*, size_type>;
+
+ void resting_potential(value_type potential_mV) {
+ resting_potential_ = potential_mV;
+ }
+
+ template <typename Cells, typename Detectors, typename Targets, typename Probes>
+ void initialize(
+ const Cells& cells, // collection of nest::mc::cell descriptions
+ Detectors& detector_handles, // (write) where to store detector handles
+ Targets& target_handles, // (write) where to store target handles
+ Probes& probe_handles); // (write) where to store probe handles
+
+ void reset();
+
+ void deliver_event(target_handle h, value_type weight) {
+ mechanisms_[synapse_base_+h.first]->net_receive(h.second, weight);
+ }
+
+ value_type detector_voltage(detector_handle h) const {
+ return voltage_[h]; // detector_handle is just the compartment index
+ }
+
+ value_type probe(probe_handle h) const {
+ return (this->*h.first)[h.second];
+ }
+
+ void advance(value_type dt);
+
+ /// Following types and methods are public only for testing:
/// the type used to store matrix information
using matrix_type = matrix<value_type, size_type>;
@@ -46,21 +175,14 @@ public:
/// ion species storage
using ion_type = mechanisms::ion<value_type, size_type>;
- /// the container used for indexes
- using index_type = memory::HostVector<size_type>;
/// view into index container
using index_view = typename index_type::view_type;
using const_index_view = typename index_type::const_view_type;
- /// the container used for values
- using vector_type = memory::HostVector<value_type>;
/// view into value container
using vector_view = typename vector_type::view_type;
using const_vector_view = typename vector_type::const_view_type;
- /// constructor
- fvm_cell(nest::mc::cell const& cell);
-
/// build the matrix for a given time step
void setup_matrix(value_type dt);
@@ -88,6 +210,10 @@ public:
vector_view voltage() { return voltage_; }
const_vector_view voltage() const { return voltage_; }
+ /// return the current in each CV
+ vector_view current() { return current_; }
+ const_vector_view current() const { return current_; }
+
std::size_t size() const { return matrix_.size(); }
/// return reference to in iterable container of the mechanisms
@@ -109,28 +235,6 @@ public:
ion_type& ion_k() { return ions_[mechanisms::ionKind::k]; }
ion_type const& ion_k() const { return ions_[mechanisms::ionKind::k]; }
- /// make a time step
- void advance(value_type dt);
-
- /// pass an event to the appropriate synapse and call net_receive
- template <typename Time>
- void apply_event(postsynaptic_spike_event<Time> e) {
- mechanisms_[synapse_index_]->net_receive(e.target.index, e.weight);
- }
-
- mechanism_type& synapses() {
- return mechanisms_[synapse_index_];
- }
-
- /// set initial states
- void initialize();
-
- /// returns the compartment index of a segment location
- int compartment_index(segment_location loc) const;
-
- /// returns voltage at a segment location
- value_type voltage(segment_location loc) const;
-
/// flags if solution is physically realistic.
/// here we define physically realistic as the voltage being within reasonable bounds.
/// use a simple test of the voltage at the soma is reasonable, i.e. in the range
@@ -140,51 +244,43 @@ public:
return (v>-1000.) && (v<1000.);
}
- /// returns current at a segment location
- value_type current(segment_location loc) const;
-
-
value_type time() const { return t_; }
- 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
+ /// current time [ms]
value_type t_ = value_type{0};
+ /// resting potential (initial voltage condition)
+ value_type resting_potential_ = -65;
+
/// the linear system for implicit time stepping of cell state
matrix_type matrix_;
/// index for fast lookup of compartment index ranges of segments
index_type segment_index_;
- /// cv_areas_[i] is the surface area of CV i
+ /// cv_areas_[i] is the surface area of CV i [µm^2]
vector_type cv_areas_;
/// alpha_[i] is the following value at the CV face between
/// CV i and its parent, required when constructing linear system
/// face_alpha_[i] = area_face / (c_m * r_L * delta_x);
- vector_type face_alpha_;
+ vector_type face_alpha_; // [µm·m^2/cm/s ≡ 10^5 µm^2/ms]
- /// cv_capacitance_[i] is the capacitance of CV i per unit area (i.e. c_m)
+ /// cv_capacitance_[i] is the capacitance of CV i per unit area (i.e. c_m) [F/m^2]
vector_type cv_capacitance_;
- /// the average current over the surface of each CV
+ /// the average current density over the surface of each CV [mA/cm^2]
/// current_ = i_m - i_e
- /// so the total current over the surface of CV i is
- /// current_[i] * cv_areas_
vector_type current_;
- /// the potential in mV in each CV
+ /// the potential in each CV [mV]
vector_type voltage_;
- std::size_t synapse_index_; // synapses at the end of mechanisms_, from here
+ /// Where point mechanisms start in the mechanisms_ list.
+ std::size_t synapse_base_;
/// the set of mechanisms present in the cell
std::vector<mechanism_type> mechanisms_;
@@ -197,7 +293,7 @@ private:
std::vector<std::pair<const vector_type fvm_cell::*, uint32_t>> probes_;
// mechanism factory
- using mechanism_catalogue = nest::mc::mechanisms::catalogue<T, I>;
+ using mechanism_catalogue = nest::mc::mechanisms::catalogue<value_type, size_type>;
};
////////////////////////////////////////////////////////////////////////////////
@@ -205,23 +301,35 @@ private:
////////////////////////////////////////////////////////////////////////////////
template <typename T, typename I>
-fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
-: cv_areas_ {cell.num_compartments(), T(0)}
-, face_alpha_ {cell.num_compartments(), T(0)}
-, cv_capacitance_{cell.num_compartments(), T(0)}
-, current_ {cell.num_compartments(), T(0)}
-, voltage_ {cell.num_compartments(), T(0)}
+template <typename Cells, typename Detectors, typename Targets, typename Probes>
+void fvm_cell<T, I>::initialize(
+ const Cells& cells,
+ Detectors& detector_handles,
+ Targets& target_handles,
+ Probes& probe_handles)
{
+ if (util::size(cells)!=1u) {
+ throw std::invalid_argument("fvm_cell accepts only one cell");
+ }
+
+ const nest::mc::cell& cell = *(std::begin(cells));
+ size_type ncomp = cell.num_compartments();
+
+ // confirm write-parameters have enough room to store handles
+ EXPECTS(util::size(detector_handles)==cell.detectors().size());
+ EXPECTS(util::size(target_handles)==cell.synapses().size());
+ EXPECTS(util::size(probe_handles)==cell.probes().size());
+
+ // initialize storage
+ cv_areas_ = vector_type{ncomp, T{0}};
+ face_alpha_ = vector_type{ncomp, T{0}};
+ cv_capacitance_ = vector_type{ncomp, T{0}};
+ current_ = vector_type{ncomp, T{0}};
+ voltage_ = vector_type{ncomp, T{0}};
+
using util::left;
using util::right;
- // TODO: potential code stink
- // matrix_ is a member, but it is not initialized with the other members
- // above because it requires the parent_index, which is calculated
- // "on the fly" by cell.model().
- // cell.model() is quite expensive, and the information it calculates is
- // used elsewhere, so we defer the intialization to inside the constructor
- // body.
const auto graph = cell.model();
matrix_ = matrix_type(graph.parent_index);
@@ -313,8 +421,6 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
// TODO : this works well for density mechanisms (e.g. ion channels), but
// does it work for point processes (e.g. synapses)?
for (auto& mech : mech_map) {
- //auto& helper = nest::mc::mechanisms::get_mechanism_helper(mech.first);
-
// calculate the number of compartments that contain the mechanism
auto num_comp = 0u;
for (auto seg : mech.second) {
@@ -338,22 +444,37 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
// instantiate the mechanism
mechanisms_.push_back(
mechanism_catalogue::make(mech.first, voltage_, current_, compartment_index)
- //helper->new_mechanism(voltage_, current_, compartment_index)
);
}
- synapse_index_ = mechanisms_.size();
+ synapse_base_ = mechanisms_.size();
+
+ // Create the synapse mechanism implementations together with the target handles
+ std::vector<std::vector<cell_lid_type>> syn_mech_map;
+ std::map<std::string, std::size_t> syn_mech_indices;
+ auto target_hi = target_handles.begin();
- std::map<std::string, std::vector<cell_lid_type>> syn_map;
for (const auto& syn : cell.synapses()) {
- syn_map[syn.mechanism.name()].push_back(find_compartment_index(syn.location, graph));
+ const auto& name = syn.mechanism.name();
+ std::size_t index = 0;
+ if (syn_mech_indices.count(name)==0) {
+ index = syn_mech_map.size();
+ syn_mech_indices[name] = index;
+ syn_mech_map.push_back(std::vector<cell_lid_type>{});
+ }
+ else {
+ index = syn_mech_indices[name];
+ }
+
+ size_type comp = find_compartment_index(syn.location, graph);
+ *target_hi++ = target_handle{index, syn_mech_map[index].size()};
+ syn_mech_map[index].push_back(comp);
}
- for (const auto& syni : syn_map) {
+ for (const auto& syni : syn_mech_indices) {
const auto& mech_name = syni.first;
- // auto& helper = nest::mc::mechanisms::get_mechanism_helper(mech_name);
- index_type compartment_index(syni.second);
+ index_type compartment_index(syn_mech_map[syni.second]);
auto mech = mechanism_catalogue::make(mech_name, voltage_, current_, compartment_index);
mech->set_areas(cv_areas_);
mechanisms_.push_back(std::move(mech));
@@ -417,24 +538,35 @@ fvm_cell<T, I>::fvm_cell(nest::mc::cell const& cell)
}
// record probe locations by index into corresponding state vector
+ auto probe_hi = probe_handles.begin();
for (auto probe : cell.probes()) {
- uint32_t comp = find_compartment_index(probe.location, graph);
+ auto comp = find_compartment_index(probe.location, graph);
switch (probe.kind) {
- case probeKind::membrane_voltage:
- probes_.push_back({&fvm_cell::voltage_, comp});
- break;
- case probeKind::membrane_current:
- probes_.push_back({&fvm_cell::current_, comp});
- break;
- default:
- throw std::logic_error("unrecognized probeKind");
+ case probeKind::membrane_voltage:
+ *probe_hi = {&fvm_cell::voltage_, comp};
+ break;
+ case probeKind::membrane_current:
+ *probe_hi = {&fvm_cell::current_, comp};
+ break;
+ default:
+ throw std::logic_error("unrecognized probeKind");
}
+ ++probe_hi;
}
+
+ // detector handles are just their corresponding compartment indices
+ auto detector_hi = detector_handles.begin();
+ for (auto detector : cell.detectors()) {
+ auto comp = find_compartment_index(detector.location, graph);
+ *detector_hi++ = comp;
+ }
+
+ // initialise mechanism and voltage state
+ reset();
}
template <typename T, typename I>
-void fvm_cell<T, I>::setup_matrix(T dt)
-{
+void fvm_cell<T, I>::setup_matrix(T dt) {
using memory::all;
// convenience accesors to matrix storage
@@ -456,9 +588,9 @@ void fvm_cell<T, I>::setup_matrix(T dt)
// . . .
// l[i] . . d[i]
//
- //d(all) = 1.0;
- d(all) = cv_areas_;
- for(auto i=1u; i<d.size(); ++i) {
+
+ d(all) = cv_areas_; // [µm^2]
+ for (auto i=1u; i<d.size(); ++i) {
auto a = 1e5*dt * face_alpha_[i];
d[i] += a;
@@ -471,50 +603,23 @@ void fvm_cell<T, I>::setup_matrix(T dt)
// the RHS of the linear system is
// V[i] - dt/cm*(im - ie)
- auto factor = 10.*dt;
+ auto factor = 10.*dt; // units: 10·ms/(F/m^2)·(mA/cm^2) ≡ mV
for(auto i=0u; i<d.size(); ++i) {
- //rhs[i] = voltage_[i] - factor/cv_capacitance_[i]*current_[i];
rhs[i] = cv_areas_[i]*(voltage_[i] - factor/cv_capacitance_[i]*current_[i]);
}
}
-template <typename T, typename I>
-int fvm_cell<T, I>::compartment_index(segment_location loc) const
-{
- EXPECTS(unsigned(loc.segment) < segment_index_.size());
-
- const auto seg = loc.segment;
-
- auto first = segment_index_[seg];
- auto n = segment_index_[seg+1] - first;
- auto index = std::floor(n*loc.position);
- return index<n ? first+index : first+n-1;
-}
template <typename T, typename I>
-T fvm_cell<T, I>::voltage(segment_location loc) const
-{
- return voltage_[compartment_index(loc)];
-}
-
-template <typename T, typename I>
-T fvm_cell<T, I>::current(segment_location loc) const
-{
- return current_[compartment_index(loc)];
-}
-
-template <typename T, typename I>
-void fvm_cell<T, I>::initialize()
-{
+void fvm_cell<T, I>::reset() {
+ voltage_(memory::all) = resting_potential_;
t_ = 0.;
-
- for(auto& m : mechanisms_) {
+ for (auto& m : mechanisms_) {
m->nrn_init();
}
}
template <typename T, typename I>
-void fvm_cell<T, I>::advance(T dt)
-{
+void fvm_cell<T, I>::advance(T dt) {
using memory::all;
PE("current");
@@ -529,17 +634,18 @@ void fvm_cell<T, I>::advance(T dt)
}
// add current contributions from stimulii
- for(auto& stim : stimulii_) {
- auto ie = stim.second.amplitude(t_);
+ for (auto& stim : stimulii_) {
+ auto ie = stim.second.amplitude(t_); // [nA]
auto loc = stim.first;
- // the factor of 100 scales the injected current to 10^2.nA
- current_[loc] -= 100.*ie/cv_areas_[loc];
+ // note: current_ in [mA/cm^2], ie in [nA], cv_areas_ in [µm^2].
+ // unit scale factor: [nA/µm^2]/[mA/cm^2] = 100
+ current_[loc] -= 100*ie/cv_areas_[loc];
}
PL();
- PE("matrix", "setup");
// solve the linear system
+ PE("matrix", "setup");
setup_matrix(dt);
PL(); PE("solve");
matrix_.solve();
@@ -547,8 +653,8 @@ void fvm_cell<T, I>::advance(T dt)
voltage_(all) = matrix_.rhs();
PL();
- PE("state");
// integrate state of gating variables etc.
+ PE("state");
for(auto& m : mechanisms_) {
PE(m->name().c_str());
m->nrn_state();
@@ -563,4 +669,3 @@ void fvm_cell<T, I>::advance(T dt)
} // namespace mc
} // namespace nest
-
diff --git a/src/spike_source.hpp b/src/spike_source.hpp
index 0537b6800a19534d23aedd7de2edae3d022b53e6..fd043a3a6252410d89fd5d64c0d9d57816f14987 100644
--- a/src/spike_source.hpp
+++ b/src/spike_source.hpp
@@ -8,13 +8,12 @@ namespace mc {
// spike detector for a lowered cell
template <typename Cell>
-class spike_detector
-{
+class spike_detector {
public:
using cell_type = Cell;
- spike_detector(const cell_type& cell, segment_location loc, double thresh, float t_init) :
- location_(loc),
+ spike_detector(const cell_type& cell, typename Cell::detector_handle h, double thresh, float t_init) :
+ handle_(h),
threshold_(thresh)
{
reset(cell, t_init);
@@ -22,7 +21,7 @@ public:
util::optional<float> test(const cell_type& cell, float t) {
util::optional<float> result = util::nothing;
- auto v = cell.voltage(location_);
+ auto v = cell.detector_voltage(handle_);
// these if statements could be simplified, but I keep them like
// this to clearly reflect the finite state machine
@@ -50,22 +49,19 @@ public:
bool is_spiking() const { return is_spiking_; }
- segment_location location() const { return location_; }
-
float t() const { return previous_t_; }
float v() const { return previous_v_; }
void reset(const cell_type& cell, float t_init) {
previous_t_ = t_init;
- previous_v_ = cell.voltage(location_);
+ previous_v_ = cell.detector_voltage(handle_);
is_spiking_ = previous_v_ >= threshold_;
}
private:
-
// parameters/data
- segment_location location_;
+ typename cell_type::detector_handle handle_;
double threshold_;
// state
diff --git a/src/stimulus.hpp b/src/stimulus.hpp
index f7c7538cf398a7154588126620d97b7c054953a4..f3c587a412eb2aba5a0dc1c127c262b86ca018e7 100644
--- a/src/stimulus.hpp
+++ b/src/stimulus.hpp
@@ -45,9 +45,9 @@ class i_clamp {
private:
- value_type delay_ = 0;
- value_type duration_ = 0;
- value_type amplitude_ = 0;
+ value_type delay_ = 0; // [ms]
+ value_type duration_ = 0; // [ms]
+ value_type amplitude_ = 0; // [nA]
};
} // namespace mc
diff --git a/src/util/meta.hpp b/src/util/meta.hpp
index 2e22e31ceda0ff83588625486ec7a36161569733..663c7a6d6e5e139aa7a25c005e4e41608dc39ede 100644
--- a/src/util/meta.hpp
+++ b/src/util/meta.hpp
@@ -2,6 +2,7 @@
/* Type utilities and convenience expressions. */
+#include <cstddef>
#include <type_traits>
namespace nest {
@@ -16,6 +17,12 @@ using result_of_t = typename std::result_of<T>::type;
template <bool V>
using enable_if_t = typename std::enable_if<V>::type;
+template <typename X>
+std::size_t size(const X& x) { return x.size(); }
+
+template <typename X, std::size_t N>
+constexpr std::size_t size(X (&)[N]) { return N; }
+
// Convenience short cuts
template <typename T>
diff --git a/src/util/singleton.hpp b/src/util/singleton.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c718eef3624af2088daa9d84abb00401ed1bd6d4
--- /dev/null
+++ b/src/util/singleton.hpp
@@ -0,0 +1,63 @@
+#pragma once
+
+/*
+ * Present a single object as a (non-owning) container with one
+ * element.
+ *
+ * (Will be subsumed by range/view code.)
+ */
+
+#include <algorithm>
+
+namespace nest {
+namespace mc {
+namespace util {
+
+template <typename X>
+struct singleton_adaptor {
+ using size_type = std::size_t;
+ using difference_type = std::ptrdiff_t;
+ using value_type = X;
+ using reference = X&;
+ using const_reference = const X&;
+ using iterator = X*;
+ using const_iterator = const X*;
+
+ X* xp;
+ singleton_adaptor(X& x): xp(&x) {}
+
+ const X* cbegin() const { return xp; }
+ const X* begin() const { return xp; }
+ X* begin() { return xp; }
+
+ const X* cend() const { return xp+1; }
+ const X* end() const { return xp+1; }
+ X* end() { return xp+1; }
+
+ std::size_t size() const { return 1u; }
+
+ bool empty() const { return false; }
+
+ const X* front() const { return *xp; }
+ X* front() { return *xp; }
+
+ const X* back() const { return *xp; }
+ X* back() { return *xp; }
+
+ const X* operator[](difference_type) const { return *xp; }
+ X* operator[](difference_type) { return *xp; }
+
+ void swap(singleton_adaptor& s) { std::swap(xp, s.xp); }
+ friend void swap(singleton_adaptor& r, singleton_adaptor& s) {
+ r.swap(s);
+ }
+};
+
+template <typename X>
+singleton_adaptor<X> singleton_view(X& x) {
+ return singleton_adaptor<X>(x);
+}
+
+} // namespace util
+} // namespace mc
+} // namespace nest
diff --git a/tests/unit/test_fvm.cpp b/tests/unit/test_fvm.cpp
index c0ed438f20841d40a9f2653e38424aeeba4a64ea..39fdb6c837f71820c3fd6981f9be2bc38cee91c8 100644
--- a/tests/unit/test_fvm.cpp
+++ b/tests/unit/test_fvm.cpp
@@ -5,6 +5,7 @@
#include <common_types.hpp>
#include <cell.hpp>
#include <fvm_cell.hpp>
+#include <util/singleton.hpp>
#include "../test_util.hpp"
@@ -42,7 +43,14 @@ TEST(fvm, cable)
cell.segment(2)->set_compartments(4);
using fvm_cell = fvm::fvm_cell<double, cell_lid_type>;
- fvm_cell fvcell(cell);
+
+ std::vector<fvm_cell::target_handle> targets;
+ std::vector<fvm_cell::detector_handle> detectors;
+ std::vector<fvm_cell::probe_handle> probes;
+
+ fvm_cell fvcell;
+ fvcell.initialize(util::singleton_view(cell), detectors, targets, probes);
+
auto& J = fvcell.jacobian();
EXPECT_EQ(cell.num_compartments(), 9u);
@@ -92,7 +100,13 @@ TEST(fvm, init)
cell.segment(1)->set_compartments(10);
using fvm_cell = fvm::fvm_cell<double, cell_lid_type>;
- fvm_cell fvcell(cell);
+ std::vector<fvm_cell::target_handle> targets;
+ std::vector<fvm_cell::detector_handle> detectors;
+ std::vector<fvm_cell::probe_handle> probes;
+
+ fvm_cell fvcell;
+ fvcell.initialize(util::singleton_view(cell), detectors, targets, probes);
+
auto& J = fvcell.jacobian();
EXPECT_EQ(J.size(), 11u);
diff --git a/tests/unit/test_probe.cpp b/tests/unit/test_probe.cpp
index 5af744d2cc5f7f9b8a403b991b22d260c9ec74d4..13a12ad6bd067ed91b37f2867f22b0f6405ba89a 100644
--- a/tests/unit/test_probe.cpp
+++ b/tests/unit/test_probe.cpp
@@ -1,8 +1,9 @@
#include "gtest.h"
-#include "common_types.hpp"
-#include "cell.hpp"
-#include "fvm_cell.hpp"
+#include <common_types.hpp>
+#include <cell.hpp>
+#include <fvm_cell.hpp>
+#include <util/singleton.hpp>
TEST(probe, instantiation)
{
@@ -49,31 +50,35 @@ TEST(probe, fvm_cell)
segment_location loc1{1, 1};
segment_location loc2{1, 0.5};
- auto pv0 = bs.add_probe({loc0, probeKind::membrane_voltage});
- auto pv1 = bs.add_probe({loc1, probeKind::membrane_voltage});
- auto pi2 = bs.add_probe({loc2, probeKind::membrane_current});
+ bs.add_probe({loc0, probeKind::membrane_voltage});
+ bs.add_probe({loc1, probeKind::membrane_voltage});
+ bs.add_probe({loc2, probeKind::membrane_current});
i_clamp stim(0, 100, 0.3);
bs.add_stimulus({1, 1}, stim);
- fvm::fvm_cell<double, cell_local_size_type> lcell(bs);
- lcell.setup_matrix(0.01);
- lcell.initialize();
+ using fvm_cell = fvm::fvm_cell<double, cell_lid_type>;
+ std::vector<fvm_cell::target_handle> targets;
+ std::vector<fvm_cell::detector_handle> detectors;
+ std::vector<fvm_cell::probe_handle> probes{3};
- EXPECT_EQ(3u, lcell.num_probes());
+ fvm_cell lcell;
+ lcell.initialize(util::singleton_view(bs), detectors, targets, probes);
- // expect probe values and direct queries of voltage and current
+ // Know from implementation that probe_handle.second
+ // is a compartment index: expect probe values and
+ // direct queries of voltage and current
// to be equal in fvm cell
- EXPECT_EQ(lcell.voltage(loc0), lcell.probe(pv0));
- EXPECT_EQ(lcell.voltage(loc1), lcell.probe(pv1));
- EXPECT_EQ(lcell.current(loc2), lcell.probe(pi2));
+ EXPECT_EQ(lcell.voltage()[probes[0].second], lcell.probe(probes[0]));
+ 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);
- EXPECT_EQ(lcell.voltage(loc0), lcell.probe(pv0));
- EXPECT_EQ(lcell.voltage(loc1), lcell.probe(pv1));
- EXPECT_EQ(lcell.current(loc2), lcell.probe(pi2));
+ EXPECT_EQ(lcell.voltage()[probes[0].second], lcell.probe(probes[0]));
+ EXPECT_EQ(lcell.voltage()[probes[1].second], lcell.probe(probes[1]));
+ EXPECT_EQ(lcell.current()[probes[2].second], lcell.probe(probes[2]));
}
diff --git a/tests/unit/test_spikes.cpp b/tests/unit/test_spikes.cpp
index 3b4862029208c1f6c5a657d83dd982f510852cd6..5546ebd30dd6021cf74d98618dd043ed0eea79ce 100644
--- a/tests/unit/test_spikes.cpp
+++ b/tests/unit/test_spikes.cpp
@@ -4,7 +4,8 @@
#include <spike_source.hpp>
struct cell_proxy {
- double voltage(nest::mc::segment_location loc) const {
+ using detector_handle = int;
+ double detector_voltage(detector_handle) const {
return v;
}
@@ -15,16 +16,17 @@ TEST(spikes, spike_detector)
{
using namespace nest::mc;
using detector_type = spike_detector<cell_proxy>;
+ using detector_handle = cell_proxy::detector_handle;
+
cell_proxy proxy;
float threshold = 10.f;
float t = 0.f;
float dt = 1.f;
- auto loc = segment_location(1, 0.1);
+ detector_handle handle{};
- auto detector = detector_type(proxy, loc, threshold, t);
+ auto detector = detector_type(proxy, handle, threshold, t);
EXPECT_FALSE(detector.is_spiking());
- EXPECT_EQ(loc, detector.location());
EXPECT_EQ(proxy.v, detector.v());
EXPECT_EQ(t, detector.t());
@@ -35,7 +37,6 @@ TEST(spikes, spike_detector)
EXPECT_FALSE(spike);
EXPECT_FALSE(detector.is_spiking());
- EXPECT_EQ(loc, detector.location());
EXPECT_EQ(proxy.v, detector.v());
EXPECT_EQ(t, detector.t());
}
@@ -49,7 +50,6 @@ TEST(spikes, spike_detector)
EXPECT_EQ(spike.get(), 1.5);
EXPECT_TRUE(detector.is_spiking());
- EXPECT_EQ(loc, detector.location());
EXPECT_EQ(proxy.v, detector.v());
EXPECT_EQ(t, detector.t());
}
@@ -62,7 +62,6 @@ TEST(spikes, spike_detector)
EXPECT_FALSE(spike);
EXPECT_FALSE(detector.is_spiking());
- EXPECT_EQ(loc, detector.location());
EXPECT_EQ(proxy.v, detector.v());
EXPECT_EQ(t, detector.t());
}
diff --git a/tests/validation/validate_ball_and_stick.cpp b/tests/validation/validate_ball_and_stick.cpp
index 62d06b4dac85885f935dd88855223551cc33fcce..59e952de58de61ae32cdbfb721ac991ec2119f8b 100644
--- a/tests/validation/validate_ball_and_stick.cpp
+++ b/tests/validation/validate_ball_and_stick.cpp
@@ -4,6 +4,7 @@
#include <common_types.hpp>
#include <cell.hpp>
#include <fvm_cell.hpp>
+#include <util/singleton.hpp>
#include "gtest.h"
#include "../test_util.hpp"
@@ -82,6 +83,11 @@ TEST(ball_and_stick, neuron_baseline)
}
};
+ using fvm_cell = fvm::fvm_cell<double, cell_local_size_type>;
+ std::vector<fvm_cell::detector_handle> detectors(cell.detectors().size());
+ std::vector<fvm_cell::target_handle> targets(cell.synapses().size());
+ std::vector<fvm_cell::probe_handle> probes(cell.probes().size());
+
std::vector<result> results;
for(auto run_index=0u; run_index<cell_data.size(); ++run_index) {
auto& run = cell_data[run_index];
@@ -90,13 +96,10 @@ TEST(ball_and_stick, neuron_baseline)
std::vector<std::vector<double>> v(3);
// make the lowered finite volume cell
- fvm::fvm_cell<double, cell_local_size_type> model(cell);
- auto graph = cell.model();
- // set initial conditions
- using memory::all;
- model.voltage()(all) = -65.;
- model.initialize(); // have to do this _after_ initial conditions are set
+ fvm_cell model;
+ model.initialize(util::singleton_view(cell), detectors, targets, probes);
+ auto graph = cell.model();
// run the simulation
auto soma_comp = nest::mc::find_compartment_index({0,0}, graph);
@@ -234,6 +237,11 @@ TEST(ball_and_3stick, neuron_baseline)
}
};
+ using fvm_cell = fvm::fvm_cell<double, cell_local_size_type>;
+ std::vector<fvm_cell::detector_handle> detectors(cell.detectors().size());
+ std::vector<fvm_cell::target_handle> targets(cell.synapses().size());
+ std::vector<fvm_cell::probe_handle> probes(cell.probes().size());
+
std::vector<result> results;
auto start = testing::tic();
for(auto run_index=0u; run_index<cell_data.size(); ++run_index) {
@@ -245,13 +253,11 @@ TEST(ball_and_3stick, neuron_baseline)
std::vector<std::vector<double>> v(3);
// make the lowered finite volume cell
- fvm::fvm_cell<double, cell_local_size_type> model(cell);
+ fvm_cell model;
+ model.initialize(util::singleton_view(cell), detectors, targets, probes);
auto graph = cell.model();
// set initial conditions
- using memory::all;
- model.voltage()(all) = -65.;
- model.initialize(); // have to do this _after_ initial conditions are set
// run the simulation
auto soma_comp = nest::mc::find_compartment_index({0,0.}, graph);
diff --git a/tests/validation/validate_soma.cpp b/tests/validation/validate_soma.cpp
index 24a0c1c59f0db24bc6f9c7832909f9f52af7beb7..af4d32826acb12bcc522daf7aa3f84c80dcc6bb2 100644
--- a/tests/validation/validate_soma.cpp
+++ b/tests/validation/validate_soma.cpp
@@ -4,6 +4,7 @@
#include <common_types.hpp>
#include <cell.hpp>
#include <fvm_cell.hpp>
+#include <util/singleton.hpp>
#include "gtest.h"
#include "../test_util.hpp"
@@ -27,7 +28,13 @@ TEST(soma, neuron_baseline)
cell.add_stimulus({0,0.5}, {10., 100., 0.1});
// make the lowered finite volume cell
- fvm::fvm_cell<double, cell_local_size_type> model(cell);
+ using fvm_cell = fvm::fvm_cell<double, cell_local_size_type>;
+ std::vector<fvm_cell::detector_handle> detectors(cell.detectors().size());
+ std::vector<fvm_cell::target_handle> targets(cell.synapses().size());
+ std::vector<fvm_cell::probe_handle> probes(cell.probes().size());
+
+ fvm_cell model;
+ model.initialize(util::singleton_view(cell), detectors, targets, probes);
// load data from file
auto cell_data = testing::g_validation_data.load("soma.json");
@@ -49,9 +56,7 @@ TEST(soma, neuron_baseline)
double dt = run["dt"];
// set initial conditions
- using memory::all;
- model.voltage()(all) = -65.;
- model.initialize(); // have to do this _after_ initial conditions are set
+ model.reset();
// run the simulation
auto tfinal = 120.; // ms
@@ -92,7 +97,13 @@ TEST(soma, convergence)
cell.add_stimulus({0,0.5}, {10., 100., 0.1});
// make the lowered finite volume cell
- fvm::fvm_cell<double, cell_local_size_type> model(cell);
+ using fvm_cell = fvm::fvm_cell<double, cell_local_size_type>;
+ std::vector<fvm_cell::detector_handle> detectors(cell.detectors().size());
+ std::vector<fvm_cell::target_handle> targets(cell.synapses().size());
+ std::vector<fvm_cell::probe_handle> probes(cell.probes().size());
+
+ fvm_cell model;
+ model.initialize(util::singleton_view(cell), detectors, targets, probes);
// generate baseline solution with small dt=0.0001
std::vector<double> baseline_spike_times;
@@ -101,9 +112,7 @@ TEST(soma, convergence)
std::vector<double> v;
// set initial conditions
- using memory::all;
- model.voltage()(all) = -65.;
- model.initialize(); // have to do this _after_ initial conditions are set
+ model.reset();
// run the simulation
auto tfinal = 120.; // ms
@@ -123,9 +132,7 @@ TEST(soma, convergence)
std::vector<double> v;
// set initial conditions
- using memory::all;
- model.voltage()(all) = -65.;
- model.initialize(); // have to do this _after_ initial conditions are set
+ model.reset();
// run the simulation
auto tfinal = 120.; // ms
diff --git a/tests/validation/validate_synapses.cpp b/tests/validation/validate_synapses.cpp
index 876b24eb55c9cad0e08c17cef9cdf16595f5bd90..9aa90abef6f6178a6285f9a1d5db2e9d1cb1d4f1 100644
--- a/tests/validation/validate_synapses.cpp
+++ b/tests/validation/validate_synapses.cpp
@@ -71,8 +71,11 @@ void run_neuron_baseline(const char* syn_type, const char* data_file)
cell.add_synapse({1, 0.5}, syn_default);
// add probes
- auto probe_soma = cell.add_probe({{0,0}, probeKind::membrane_voltage});
- auto probe_dend = cell.add_probe({{1,0.5}, probeKind::membrane_voltage});
+ auto probe_soma_idx = cell.add_probe({{0,0}, probeKind::membrane_voltage});
+ auto probe_dend_idx = cell.add_probe({{1,0.5}, probeKind::membrane_voltage});
+
+ cell_member_type probe_soma{0u, probe_soma_idx};
+ cell_member_type probe_dend{0u, probe_dend_idx};
// injected spike events
postsynaptic_spike_event<float> synthetic_events[] = {
@@ -111,15 +114,15 @@ void run_neuron_baseline(const char* syn_type, const char* data_file)
group.enqueue_events(synthetic_events);
// run the simulation
- v[0].push_back(group.cell().probe(probe_soma));
- v[1].push_back(group.cell().probe(probe_dend));
+ v[0].push_back(group.probe(probe_soma));
+ v[1].push_back(group.probe(probe_dend));
double t = 0.;
while(t < tfinal) {
t += dt;
group.advance(t, dt);
// save voltage at soma and dendrite
- v[0].push_back(group.cell().probe(probe_soma));
- v[1].push_back(group.cell().probe(probe_dend));
+ v[0].push_back(group.probe(probe_soma));
+ v[1].push_back(group.probe(probe_dend));
}
results.push_back({num_compartments, dt, v, measurements});