diff --git a/src/cell_group.hpp b/src/cell_group.hpp
index e70d4b0bb5826b6ccde449f63a7fad6dab29c783..a52dc36d953bc35a557ad0688b9d94f7be46639c 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/compartment.hpp b/src/compartment.hpp
index da6bbcf57298593a60f7510822832060c7a89846..326563ad4261bd6d3e5ce0b516c3aa8808c893a9 100644
--- a/src/compartment.hpp
+++ b/src/compartment.hpp
@@ -3,7 +3,10 @@
#include <iterator>
#include <utility>
-#include "common_types.hpp"
+#include <common_types.hpp>
+#include <util/counter.hpp>
+#include <util/span.hpp>
+#include <util/transform.hpp>
namespace nest {
namespace mc {
@@ -34,6 +37,7 @@ struct compartment {
value_type length;
};
+
/// The simplest type of compartment iterator :
/// - divide a segment into n compartments of equal length
/// - assume that the radius varies linearly from one end of the segment
@@ -161,6 +165,43 @@ public:
real_type length_;
};
+// (NB: auto type deduction and lambda in C++14 will simplify the following)
+
+template <typename size_type, typename real_type>
+class compartment_maker {
+public:
+ compartment_maker(size_type n, real_type length, real_type rL, real_type rR):
+ r0_{rL},
+ dr_{(rR-rL)/n},
+ dx_{length/n}
+ {}
+
+ compartment operator()(size_type i) const {
+ return compartment(i, dx_, r0_+i*dr_, r0_+(i+1)*dr_);
+ }
+
+private:
+ real_type r0_;
+ real_type dr_;
+ real_type dx_;
+};
+
+template <typename size_type, typename real_type>
+using compartment_iterator_bis =
+ util::transform_iterator<util::counter<size_type>, compartment_maker<size_type, real_type>>;
+
+template <typename size_type, typename real_type>
+util::range<compartment_iterator_bis<size_type, real_type>> make_compartment_range(
+ size_type num_compartments,
+ real_type radius_L,
+ real_type radius_R,
+ real_type length)
+{
+ return util::transform_view(
+ util::span<size_type>(0, num_compartments),
+ compartment_maker<size_type, real_type>(num_compartments, length, radius_L, radius_R));
+}
+
} // namespace mc
} // namespace nest
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/fvm_multicell.hpp b/src/fvm_multicell.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e7bc12bb3bcb05242438c118249daca3a31527df
--- /dev/null
+++ b/src/fvm_multicell.hpp
@@ -0,0 +1,691 @@
+#pragma once
+
+#include <algorithm>
+#include <iterator>
+#include <map>
+#include <set>
+#include <string>
+#include <vector>
+
+#include <algorithms.hpp>
+#include <cell.hpp>
+#include <event_queue.hpp>
+#include <ion.hpp>
+#include <math.hpp>
+#include <matrix.hpp>
+#include <mechanism.hpp>
+#include <mechanism_catalogue.hpp>
+#include <profiling/profiler.hpp>
+#include <segment.hpp>
+#include <stimulus.hpp>
+#include <util.hpp>
+#include <util/meta.hpp>
+
+#include <vector/include/Vector.hpp>
+
+namespace nest {
+namespace mc {
+namespace fvm {
+
+template <typename Value, typename Index>
+class fvm_cell {
+public:
+ fvm_cell() = default;
+
+ /// the real number type
+ using value_type = Value;
+
+ /// the integral index type
+ 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>;
+
+ /// mechanism type
+ using mechanism_type =
+ nest::mc::mechanisms::mechanism_ptr<value_type, size_type>;
+
+ /// ion species storage
+ using ion_type = mechanisms::ion<value_type, size_type>;
+
+ /// view into index container
+ using index_view = typename index_type::view_type;
+ using const_index_view = typename index_type::const_view_type;
+
+ /// view into value container
+ using vector_view = typename vector_type::view_type;
+ using const_vector_view = typename vector_type::const_view_type;
+
+ /// build the matrix for a given time step
+ void setup_matrix(value_type dt);
+
+ /// which requires const_view in the vector library
+ const matrix_type& jacobian() {
+ return matrix_;
+ }
+
+ /// return list of CV areas in :
+ /// um^2
+ /// 1e-6.mm^2
+ /// 1e-8.cm^2
+ const_vector_view cv_areas() const {
+ return cv_areas_;
+ }
+
+ /// return the capacitance of each CV surface
+ /// this is the total capacitance, not per unit area,
+ /// i.e. equivalent to sigma_i * c_m
+ const_vector_view cv_capacitance() const {
+ return cv_capacitance_;
+ }
+
+ /// return the voltage in each CV
+ 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
+ std::vector<mechanism_type>& mechanisms() { return mechanisms_; }
+
+ /// return reference to list of ions
+ std::map<mechanisms::ionKind, ion_type>& ions() { return ions_; }
+ std::map<mechanisms::ionKind, ion_type> const& ions() const { return ions_; }
+
+ /// return reference to sodium ion
+ ion_type& ion_na() { return ions_[mechanisms::ionKind::na]; }
+ ion_type const& ion_na() const { return ions_[mechanisms::ionKind::na]; }
+
+ /// return reference to calcium ion
+ ion_type& ion_ca() { return ions_[mechanisms::ionKind::ca]; }
+ ion_type const& ion_ca() const { return ions_[mechanisms::ionKind::ca]; }
+
+ /// return reference to pottasium ion
+ ion_type& ion_k() { return ions_[mechanisms::ionKind::k]; }
+ ion_type const& ion_k() const { return ions_[mechanisms::ionKind::k]; }
+
+ /// 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
+ /// v_soma \in (-1000mv, 1000mv)
+ bool is_physical_solution() const {
+ auto v = voltage_[0];
+ return (v>-1000.) && (v<1000.);
+ }
+
+ value_type time() const { return t_; }
+
+ std::size_t num_probes() const { return probes_.size(); }
+
+private:
+ /// 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 [µ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_; // [µ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) [F/m^2]
+ vector_type cv_capacitance_;
+
+ /// the average current density over the surface of each CV [mA/cm^2]
+ /// current_ = i_m - i_e
+ vector_type current_;
+
+ /// the potential in each CV [mV]
+ vector_type voltage_;
+
+ /// 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_;
+
+ /// the ion species
+ std::map<mechanisms::ionKind, ion_type> ions_;
+
+ std::vector<std::pair<uint32_t, i_clamp>> stimulii_;
+
+ std::vector<std::pair<const vector_type fvm_cell::*, uint32_t>> probes_;
+
+ // mechanism factory
+ using mechanism_catalogue = nest::mc::mechanisms::catalogue<value_type, size_type>;
+
+ // perform area and capacitance calculation on initialization
+ void compute_cv_area_capacitance(
+ std::pair<size_type, size_type> comps, const segment& seg, index_vector &parent);
+};
+
+////////////////////////////////////////////////////////////////////////////////
+//////////////////////////////// Implementation ////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////
+
+template <typename T, typename I>
+void fvm_cell<T, I>::compute_cv_area_capacitance(
+ std::pair<size_type, size_type> comps,
+ const segment& seg,
+ index_vector &parent)
+{
+ // precondition: group_parent_index[j] holds the correct value for
+ // j in [base_comp, base_comp+segment.num_compartments()].
+
+ if (auto soma = seg->as_soma()) {
+ // confirm assumption that there is one compartment in soma
+ if (comps.size()!=1)
+ throw std::logic_error("soma allocated more than one compartment");
+ }
+ auto i = comps.first;
+ auto area = math::area_sphere(soma->radius());
+
+ cv_areas_[i] += area;
+ cv_capacitance_[i] += area * soma->mechanism("membrane").get("c_m").value;
+ }
+ else if (auto cable = s->as_cable()) {
+ // loop over each compartment in the cable
+ // each compartment has the face between two CVs at its centre
+ // the centers of the CVs are the end points of the compartment
+ //
+ // __________________________________
+ // | ........ | .cvleft. | cv |
+ // | ........ L ........ C R
+ // |__________|__________|__________|
+ //
+ // The compartment has end points marked L and R (left and right).
+ // The left compartment is assumed to be closer to the soma
+ // (i.e. it follows the minimal degree ordering)
+ // The face is at the center, marked C.
+ // The full control volume to the left (marked with .)
+
+ auto c_m = cable->mechanism("membrane").get("c_m").value;
+ auto r_L = cable->mechanism("membrane").get("r_L").value;
+ const auto& compartments = cable->compartments();
+
+ EXPECTS(util::size(compartments)==comps.second-comps.first);
+
+ for (auto i: util::make_span(comps)) {
+ const auto c& = compartments[i-comps.first];
+ auto j = parent[i];
+
+ auto radius_center = math::mean(c.radius);
+ auto area_face = math::area_circle(radius_center);
+ face_alpha_[i] = area_face / (c_m * r_L * c.length);
+
+ auto halflen = c.length/2;
+ auto al = math::area_frustrum(halflen, left(c.radius), radius_center);
+ auto ar = math::area_frustrum(halflen, right(c.radius), radius_center);
+
+ cv_areas_[j] += al;
+ cv_areas_[i] += ar;
+ cv_capacitance_[j] += al * c_m;
+ cv_capacitance_[i] += ar * c_m;
+ }
+ }
+ else {
+ throw std::domain_error("FVM lowering encountered unsuported segment type");
+ }
+}
+
+
+template <typename T, typename I>
+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)
+{
+ using util::transform_view;
+ using util::make_partition;
+ using util::make_span;
+ using util::size;
+
+ auto cell_num_compartments =
+ transform_view(cells, [](const cell& c) { return c.num_compartments(); });
+
+ std::vector<cell_lid_type> cell_comp_bounds;
+ auto cell_comp_part = make_partition(cell_comp_bounds, cell_num_compartments);
+ auto ncomp = cell_comp_part.bounds().second;
+
+ index_vector group_parent_index{cell_comp_part.bounds().second, 0};
+
+ auto comps = cell_comp_part.begin();
+ for (const auto& c: cells) {
+ auto parent_index = c.model().parent_index;
+
+ for (auto j: make_span(*comps)) {
+ group_parent_index[j] = parent_index[j-comps->first]+comps->first;
+ }
+
+ auto seg_num_compartments =
+ transform_view(c.segments(), [](const segment& s) { return s.num_compartments(); });
+
+ std::vector<cell_lid_type> seg_comp_bounds;
+ auto seg_comp_part = make_partition(seg_comp_bounds, seg_num_compartments, comps->first);
+
+ auto seg_comps = seg_comp_part.begin();
+ for (auto j=0; j<size(seg_comp_part); ++j) {
+ calculate_cv_area_capacitance(
+ }
+
+ ++comps;
+ }
+
+
+
+
+
+ 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;
+
+ const auto graph = cell.model();
+ matrix_ = matrix_type(graph.parent_index);
+
+ auto parent_index = matrix_.p();
+ segment_index_ = graph.segment_index;
+
+ auto seg_idx = 0;
+ for(auto const& s : cell.segments()) {
+ if(auto soma = s->as_soma()) {
+ // assert the assumption that the soma is at 0
+ if(seg_idx!=0) {
+ throw std::domain_error(
+ "FVM lowering encountered soma with non-zero index"
+ );
+ }
+ auto area = math::area_sphere(soma->radius());
+ cv_areas_[0] += area;
+ cv_capacitance_[0] += area * soma->mechanism("membrane").get("c_m").value;
+ }
+ else if(auto cable = s->as_cable()) {
+ // loop over each compartment in the cable
+ // each compartment has the face between two CVs at its centre
+ // the centers of the CVs are the end points of the compartment
+ //
+ // __________________________________
+ // | ........ | .cvleft. | cv |
+ // | ........ L ........ C R
+ // |__________|__________|__________|
+ //
+ // The compartment has end points marked L and R (left and right).
+ // The left compartment is assumed to be closer to the soma
+ // (i.e. it follows the minimal degree ordering)
+ // The face is at the center, marked C.
+ // The full control volume to the left (marked with .)
+ auto c_m = cable->mechanism("membrane").get("c_m").value;
+ auto r_L = cable->mechanism("membrane").get("r_L").value;
+ for(auto c : cable->compartments()) {
+ auto i = segment_index_[seg_idx] + c.index;
+ auto j = parent_index[i];
+
+ auto radius_center = math::mean(c.radius);
+ auto area_face = math::area_circle( radius_center );
+ face_alpha_[i] = area_face / (c_m * r_L * c.length);
+
+ auto halflen = c.length/2;
+
+ auto al = math::area_frustrum(halflen, left(c.radius), radius_center);
+ auto ar = math::area_frustrum(halflen, right(c.radius), radius_center);
+ cv_areas_[j] += al;
+ cv_areas_[i] += ar;
+ cv_capacitance_[j] += al * c_m;
+ cv_capacitance_[i] += ar * c_m;
+ }
+ }
+ else {
+ throw std::domain_error("FVM lowering encountered unsuported segment type");
+ }
+ ++seg_idx;
+ }
+
+ // normalize the capacitance by cv_area
+ for(auto i=0u; i<size(); ++i) {
+ cv_capacitance_[i] /= cv_areas_[i];
+ }
+
+ /////////////////////////////////////////////
+ // create mechanisms
+ /////////////////////////////////////////////
+
+ // FIXME : candidate for a private member function
+
+ // for each mechanism in the cell record the indexes of the segments that
+ // contain the mechanism
+ std::map<std::string, std::vector<unsigned>> mech_map;
+
+ for (unsigned i=0; i<cell.num_segments(); ++i) {
+ for (const auto& mech : cell.segment(i)->mechanisms()) {
+ // FIXME : Membrane has to be a proper mechanism,
+ // because it is exposed via the public interface.
+ // This if statement is bad
+ if(mech.name() != "membrane") {
+ mech_map[mech.name()].push_back(i);
+ }
+ }
+ }
+
+ // Create the mechanism implementations with the state for each mechanism
+ // instance.
+ // 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) {
+ // calculate the number of compartments that contain the mechanism
+ auto num_comp = 0u;
+ for (auto seg : mech.second) {
+ num_comp += segment_index_[seg+1] - segment_index_[seg];
+ }
+
+ // build a vector of the indexes of the compartments that contain
+ // the mechanism
+ index_type compartment_index(num_comp);
+ auto pos = 0u;
+ for (auto seg : mech.second) {
+ auto seg_size = segment_index_[seg+1] - segment_index_[seg];
+ std::iota(
+ compartment_index.data() + pos,
+ compartment_index.data() + pos + seg_size,
+ segment_index_[seg]
+ );
+ pos += seg_size;
+ }
+
+ // instantiate the mechanism
+ mechanisms_.push_back(
+ mechanism_catalogue::make(mech.first, voltage_, current_, compartment_index)
+ );
+ }
+
+ 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();
+
+ for (const auto& syn : cell.synapses()) {
+ 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_mech_indices) {
+ const auto& mech_name = syni.first;
+
+ 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));
+ }
+
+
+ /////////////////////////////////////////////
+ // build the ion species
+ /////////////////////////////////////////////
+ for(auto ion : mechanisms::ion_kinds()) {
+ // find the compartment indexes of all compartments that have a
+ // mechanism that depends on/influences ion
+ std::set<int> index_set;
+ for(auto& mech : mechanisms_) {
+ if(mech->uses_ion(ion)) {
+ for(auto idx : mech->node_index()) {
+ index_set.insert(idx);
+ }
+ }
+ }
+ std::vector<cell_lid_type> indexes(index_set.begin(), index_set.end());
+
+ // create the ion state
+ if(indexes.size()) {
+ ions_.emplace(ion, index_type(indexes));
+ }
+
+ // join the ion reference in each mechanism into the cell-wide ion state
+ for(auto& mech : mechanisms_) {
+ if(mech->uses_ion(ion)) {
+ mech->set_ion(ion, ions_[ion]);
+ }
+ }
+ }
+
+ // FIXME: Hard code parameters for now.
+ // Take defaults for reversal potential of sodium and potassium from
+ // the default values in Neuron.
+ // Neuron's defaults are defined in the file
+ // nrn/src/nrnoc/membdef.h
+ auto all = memory::all;
+
+ constexpr value_type DEF_vrest = -65.0; // same name as #define in Neuron
+
+ ion_na().reversal_potential()(all) = 115+DEF_vrest; // mV
+ ion_na().internal_concentration()(all) = 10.0; // mM
+ ion_na().external_concentration()(all) = 140.0; // mM
+
+ ion_k().reversal_potential()(all) = -12.0+DEF_vrest;// mV
+ ion_k().internal_concentration()(all) = 54.4; // mM
+ ion_k().external_concentration()(all) = 2.5; // mM
+
+ ion_ca().reversal_potential()(all) = 12.5 * std::log(2.0/5e-5);// mV
+ ion_ca().internal_concentration()(all) = 5e-5; // mM
+ ion_ca().external_concentration()(all) = 2.0; // mM
+
+ // add the stimulii
+ for(const auto& stim : cell.stimulii()) {
+ auto idx = find_compartment_index(stim.location, graph);
+ stimulii_.push_back( {idx, stim.clamp} );
+ }
+
+ // record probe locations by index into corresponding state vector
+ auto probe_hi = probe_handles.begin();
+ for (auto probe : cell.probes()) {
+ auto comp = find_compartment_index(probe.location, graph);
+ switch (probe.kind) {
+ 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) {
+ using memory::all;
+
+ // convenience accesors to matrix storage
+ auto l = matrix_.l();
+ auto d = matrix_.d();
+ auto u = matrix_.u();
+ auto p = matrix_.p();
+ auto rhs = matrix_.rhs();
+
+ // The matrix has the following layout in memory
+ // where j is the parent index of i, i.e. i<j
+ //
+ // d[i] is the diagonal entry at a_ii
+ // u[i] is the upper triangle entry at a_ji
+ // l[i] is the lower triangle entry at a_ij
+ //
+ // d[j] . . u[i]
+ // . . .
+ // . . .
+ // l[i] . . d[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;
+ l[i] = -a;
+ u[i] = -a;
+
+ // add contribution to the diagonal of parent
+ d[p[i]] += a;
+ }
+
+ // the RHS of the linear system is
+ // V[i] - dt/cm*(im - ie)
+ auto factor = 10.*dt; // units: 10·ms/(F/m^2)·(mA/cm^2) ≡ mV
+ for(auto i=0u; i<d.size(); ++i) {
+ rhs[i] = cv_areas_[i]*(voltage_[i] - factor/cv_capacitance_[i]*current_[i]);
+ }
+}
+
+template <typename T, typename I>
+void fvm_cell<T, I>::reset() {
+ voltage_(memory::all) = resting_potential_;
+ t_ = 0.;
+ for (auto& m : mechanisms_) {
+ m->nrn_init();
+ }
+}
+
+template <typename T, typename I>
+void fvm_cell<T, I>::advance(T dt) {
+ using memory::all;
+
+ PE("current");
+ current_(all) = 0.;
+
+ // update currents from ion channels
+ for(auto& m : mechanisms_) {
+ PE(m->name().c_str());
+ m->set_params(t_, dt);
+ m->nrn_current();
+ PL();
+ }
+
+ // add current contributions from stimulii
+ for (auto& stim : stimulii_) {
+ auto ie = stim.second.amplitude(t_); // [nA]
+ auto loc = stim.first;
+
+ // 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();
+
+ // solve the linear system
+ PE("matrix", "setup");
+ setup_matrix(dt);
+ PL(); PE("solve");
+ matrix_.solve();
+ PL();
+ voltage_(all) = matrix_.rhs();
+ PL();
+
+ // integrate state of gating variables etc.
+ PE("state");
+ for(auto& m : mechanisms_) {
+ PE(m->name().c_str());
+ m->nrn_state();
+ PL();
+ }
+ PL();
+
+ t_ += dt;
+}
+
+} // namespace fvm
+} // 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 7b1fb0d6f228175377242b607a111457ed90b414..e949e251dc2f92101555cd7845931728c5b19f0b 100644
--- a/src/util/meta.hpp
+++ b/src/util/meta.hpp
@@ -3,6 +3,7 @@
/* Type utilities and convenience expressions. */
#include <iterator>
+#include <cstddef>
#include <type_traits>
namespace nest {
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_compartments.cpp b/tests/unit/test_compartments.cpp
index c167b16cce5585a43efc13fc9f136025b21f7811..176bf18177d2379d2c42d3fc31199f904f22bc52 100644
--- a/tests/unit/test_compartments.cpp
+++ b/tests/unit/test_compartments.cpp
@@ -131,3 +131,28 @@ TEST(compartments, compartment_range)
EXPECT_EQ(rng.begin(), rng.end());
}
}
+
+TEST(compartments, make_compartment_range)
+{
+ using namespace nest::mc;
+ auto rng = make_compartment_range(10, 1.0, 2.0, 10.);
+
+ EXPECT_EQ((*rng.begin()).index, 0u);
+ EXPECT_EQ((*rng.end()).index, 10u);
+ EXPECT_NE(rng.begin(), rng.end());
+
+ unsigned count = 0;
+ for (auto c : rng) {
+ EXPECT_EQ(c.index, count);
+ auto er = 1.0 + double(count)/10.;
+ EXPECT_DOUBLE_EQ(left(c.radius), er);
+ EXPECT_DOUBLE_EQ(right(c.radius), er+0.1);
+ EXPECT_EQ(c.length, 1.0);
+ ++count;
+ }
+ EXPECT_EQ(count, 10u);
+
+ // test case of zero length range
+ auto rng_empty = make_compartment_range(0, 1.0, 1.0, 0.);
+ EXPECT_EQ(rng_empty.begin(), rng_empty.end());
+}
diff --git a/tests/unit/test_fvm.cpp b/tests/unit/test_fvm.cpp
index c0ed438f20841d40a9f2653e38424aeeba4a64ea..ae7f0b5a3c15c01fc2c952c3ca84ff860bd10294 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/range.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..cf7562c94fbf9ddcf33e16308d418dbdddd0ebf7 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/range.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});