diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp
index 6dea2d3991bb90191d8470dd16b64cc46dbe3d02..65088e389a70454a0fcfcb3e7ae806eda8527ca6 100644
--- a/arbor/cable_cell.cpp
+++ b/arbor/cable_cell.cpp
@@ -7,6 +7,7 @@
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/mprovider.hpp>
#include <arbor/segment.hpp>
+#include <arbor/util/pp_util.hpp>
#include "util/piecewise.hpp"
#include "util/rangeutil.hpp"
@@ -22,16 +23,15 @@ using value_type = cable_cell::value_type;
using index_type = cable_cell::index_type;
using size_type = cable_cell::size_type;
+template <typename T> struct constant_type {
+ template <typename> using type = T;
+};
+
struct cable_cell_impl {
using value_type = cable_cell::value_type;
using index_type = cable_cell::index_type;
using size_type = cable_cell::size_type;
- using stimulus_instance = cable_cell::stimulus_instance;
- using synapse_instance = cable_cell::synapse_instance;
- using gap_junction_instance = cable_cell::gap_junction_instance;
- using detector_instance = cable_cell::detector_instance;
-
cable_cell_impl(const arb::morphology& m,
const label_dict& dictionary,
bool compartments_from_discretization):
@@ -105,11 +105,9 @@ struct cable_cell_impl {
cable_cell_impl(const cable_cell_impl& other):
parents(other.parents),
- stimuli(other.stimuli),
- synapses(other.synapses),
- gap_junction_sites(other.gap_junction_sites),
- spike_detectors(other.spike_detectors),
- provider(other.provider)
+ provider(other.provider),
+ region_map(other.region_map),
+ location_map(other.location_map)
{
// unique_ptr's cannot be copy constructed, do a manual assignment
segments.reserve(other.segments.size());
@@ -126,75 +124,101 @@ struct cable_cell_impl {
// the segments
std::vector<segment_ptr> segments;
- // the stimuli
- std::vector<stimulus_instance> stimuli;
+ // Embedded morphology and labelled region/locset lookup.
+ mprovider provider;
- // the synapses
- std::vector<synapse_instance> synapses;
+ // Regional assignments.
+ cable_cell_region_map region_map;
- // the gap_junctions
- std::vector<gap_junction_instance> gap_junction_sites;
+ // Point assignments.
+ cable_cell_location_map location_map;
- // the sensors
- std::vector<detector_instance> spike_detectors;
+ // Track number of point assignments by type for lid/target numbers.
+ dynamic_typed_map<constant_type<cell_lid_type>::type> placed_count;
- // Embedded morphology and labelled region/locset lookup.
- mprovider provider;
+ template <typename T>
+ mlocation_map<T>& get_location_map(const T&) {
+ return location_map.get<T>();
+ }
- template <typename Desc, typename T>
- lid_range place(const mlocation_list& locs, const Desc& desc, std::vector<T>& list) {
- const auto first = list.size();
+ mlocation_map<mechanism_desc>& get_location_map(const mechanism_desc& desc) {
+ return location_map.get<mechanism_desc>()[desc.name()];
+ }
+
+ template <typename Item>
+ lid_range place(const locset& ls, const Item& item) {
+ auto& mm = get_location_map(item);
+ cell_lid_type& lid = placed_count.get<Item>();
+ cell_lid_type first = lid;
- for (auto loc: locs) {
- list.push_back({loc, desc});
+ for (auto l: thingify(ls, provider)) {
+ placed<Item> p{l, lid++, item};
+ mm.push_back(p);
}
+ return lid_range(first, lid);
+ }
- return lid_range(first, list.size());
+ void assert_valid_segment(index_type i) const {
+ if (i>=segments.size()) {
+ throw cable_cell_error("no such segment");
+ }
}
- template <typename Desc, typename T>
- lid_range place(const locset& locs, const Desc& desc, std::vector<T>& list) {
- return place(thingify(locs, provider), desc, list);
+ void paint_segment(segment_ptr& s, const mechanism_desc& p) {
+ s->add_mechanism(p);
}
- lid_range place_gj(const mlocation_list& locs) {
- const auto first = gap_junction_sites.size();
+ void paint_segment(segment_ptr& s, init_membrane_potential p) {
+ s->parameters.init_membrane_potential = p.value;
+ }
- gap_junction_sites.insert(gap_junction_sites.end(), locs.begin(), locs.end());
+ void paint_segment(segment_ptr& s, axial_resistivity p) {
+ s->parameters.axial_resistivity = p.value;
+ }
- return lid_range(first, gap_junction_sites.size());
+ void paint_segment(segment_ptr& s, temperature_K p) {
+ s->parameters.temperature_K = p.value;
}
- lid_range place_gj(const locset& locs) {
- return place_gj(thingify(locs, provider));
+ void paint_segment(segment_ptr& s, membrane_capacitance p) {
+ s->parameters.membrane_capacitance = p.value;
}
- void assert_valid_segment(index_type i) const {
- if (i>=segments.size()) {
- throw cable_cell_error("no such segment");
- }
+ void paint_segment(segment_ptr& s, initial_ion_data p) {
+ s->parameters.ion_data[p.ion] = p.initial;
+ }
+
+ template <typename T>
+ mcable_map<T>& get_region_map(const T&) {
+ return region_map.get<T>();
}
- bool valid_location(const mlocation& loc) const {
- return test_invariants(loc) && loc.branch<segments.size();
+ mcable_map<mechanism_desc>& get_region_map(const mechanism_desc& desc) {
+ return region_map.get<mechanism_desc>()[desc.name()];
}
- template <typename F>
- void paint(const mcable_list& cables, F&& f) {
+ mcable_map<initial_ion_data>& get_region_map(const initial_ion_data& init) {
+ return region_map.get<initial_ion_data>()[init.ion];
+ }
+
+ template <typename Property>
+ void paint(const region& reg, const Property& prop) {
+ mcable_list cables = thingify(reg, provider);
+ auto& mm = get_region_map(prop);
+
for (auto c: cables) {
+ if (!mm.insert(c, prop)) {
+ throw cable_cell_error(util::pprintf("cable {} overpaints", c));
+ }
+
if (c.prox_pos!=0 || c.dist_pos!=1) {
throw cable_cell_error(util::pprintf(
"cable_cell does not support regions with partial branches: {}", c));
}
assert_valid_segment(c.branch);
- f(segments[c.branch]);
+ paint_segment(segments[c.branch], prop);
}
}
-
- template <typename F>
- void paint(const region& reg, F&& f) {
- paint(thingify(reg, provider), std::forward<F>(f));
- }
};
using impl_ptr = std::unique_ptr<cable_cell_impl, void (*)(cable_cell_impl*)>;
@@ -249,22 +273,6 @@ bool cable_cell::has_soma() const {
return !segment(0)->is_placeholder();
}
-const std::vector<cable_cell::gap_junction_instance>& cable_cell::gap_junction_sites() const {
- return impl_->gap_junction_sites;
-}
-
-const std::vector<cable_cell::synapse_instance>& cable_cell::synapses() const {
- return impl_->synapses;
-}
-
-const std::vector<cable_cell::detector_instance>& cable_cell::detectors() const {
- return impl_->spike_detectors;
-}
-
-const std::vector<cable_cell::stimulus_instance>& cable_cell::stimuli() const {
- return impl_->stimuli;
-}
-
const concrete_embedding& cable_cell::embedding() const {
return impl_->provider.embedding();
}
@@ -277,61 +285,32 @@ const mprovider& cable_cell::provider() const {
return impl_->provider;
}
-
-//
-// Painters.
-//
-// Implementation of user API for painting density channel and electrical properties on cells.
-//
-
-void cable_cell::paint(const region& target, mechanism_desc desc) {
- impl_->paint(target,
- [&desc](segment_ptr& s){return s->add_mechanism(desc);});
-}
-
-void cable_cell::paint(const region& target, cable_cell_local_parameter_set params) {
- impl_->paint(target,
- [¶ms](segment_ptr& s){return s->parameters = params;});
-}
-
-//
-// Placers.
-//
-// Implementation of user API for placing discrete items on cell morphology,
-// such as synapses, spike detectors and stimuli.
-//
-
-//
-// Synapses.
-//
-
-lid_range cable_cell::place(const locset& ls, const mechanism_desc& desc) {
- return impl_->place(ls, desc, impl_->synapses);
+const cable_cell_location_map& cable_cell::location_assignments() const {
+ return impl_->location_map;
}
-//
-// Stimuli.
-//
-
-lid_range cable_cell::place(const locset& ls, const i_clamp& desc) {
- return impl_->place(ls, desc, impl_->stimuli);
+const cable_cell_region_map& cable_cell::region_assignments() const {
+ return impl_->region_map;
}
-//
-// Gap junctions.
-//
+// Forward paint methods to implementation class.
-lid_range cable_cell::place(const locset& ls, gap_junction_site) {
- return impl_->place_gj(ls);
+#define FWD_PAINT(proptype)\
+void cable_cell::paint(const region& target, proptype prop) {\
+ impl_->paint(target, prop);\
}
+ARB_PP_FOREACH(FWD_PAINT,\
+ mechanism_desc, init_membrane_potential, axial_resistivity,\
+ temperature_K, membrane_capacitance, initial_ion_data)
-//
-// Spike detectors.
-//
+// Forward place methods to implementation class.
-lid_range cable_cell::place(const locset& ls, const threshold_detector& desc) {
- return impl_->place(ls, desc.threshold, impl_->spike_detectors);
+#define FWD_PLACE(proptype)\
+lid_range cable_cell::place(const locset& target, proptype prop) {\
+ return impl_->place(target, prop);\
}
+ARB_PP_FOREACH(FWD_PLACE,\
+ mechanism_desc, i_clamp, gap_junction_site, threshold_detector)
//
// TODO: deprectate the following as soon as discretization code catches up with em_morphology
diff --git a/arbor/cable_cell_param.cpp b/arbor/cable_cell_param.cpp
index f6843b80b462111474cfbc4cdb6f78934c1112b8..bc7e26152b127f40b26f15514b5669fc36c6ea4c 100644
--- a/arbor/cable_cell_param.cpp
+++ b/arbor/cable_cell_param.cpp
@@ -50,14 +50,12 @@ void check_global_properties(const cable_cell_global_properties& G) {
}
}
}
+ util::optional<double> init_membrane_potential; // [mV]
+ util::optional<double> temperature_K; // [K]
+ util::optional<double> axial_resistivity; // [Ω·cm]
+ util::optional<double> membrane_capacitance; // [F/m²]
-cable_cell_local_parameter_set neuron_parameter_defaults = {
- // ion defaults: internal concentration [mM], external concentration [mM], reversal potential [mV]
- {
- {"na", {10.0, 140.0, 115 - 65.}},
- {"k", {54.4, 2.5, -12 - 65.}},
- {"ca", {5e-5, 2.0, 12.5*std::log(2.0/5e-5)}}
- },
+cable_cell_parameter_set neuron_parameter_defaults = {
// initial membrane potential [mV]
-65.0,
// temperatue [K]
@@ -65,7 +63,14 @@ cable_cell_local_parameter_set neuron_parameter_defaults = {
// axial resistivity [Ω·cm]
35.4,
// membrane capacitance [F/m²]
- 0.01
+ 0.01,
+ // ion defaults:
+ // internal concentration [mM], external concentration [mM], reversal potential [mV]
+ {
+ {"na", {10.0, 140.0, 115 - 65.}},
+ {"k", {54.4, 2.5, -12 - 65.}},
+ {"ca", {5e-5, 2.0, 12.5*std::log(2.0/5e-5)}}
+ },
};
// Discretization policy implementations:
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index c0c5a1d8a46af8db9372f4a4096a1d7a01a403ef..deb1d7806c7bfa9f4edcf355cb0df59009b90e54 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -5,6 +5,7 @@
#include <arbor/arbexcept.hpp>
#include <arbor/cable_cell.hpp>
+#include <arbor/morph/mcable_map.hpp>
#include <arbor/util/optional.hpp>
#include "algorithms.hpp"
@@ -14,6 +15,7 @@
#include "util/maputil.hpp"
#include "util/meta.hpp"
#include "util/partition.hpp"
+#include "util/piecewise.hpp"
#include "util/rangeutil.hpp"
#include "util/transform.hpp"
@@ -549,15 +551,16 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
};
auto cell_segment_part = D.cell_segment_part();
- size_type target_id = 0;
+ size_type target_count = 0;
for (auto cell_idx: make_span(0, D.ncell)) {
auto& cell = cells[cell_idx];
auto seg_range = cell_segment_part[cell_idx];
+ size_type target_offset = target_count;
auto add_ion_segment =
[&gparam, &cell, &ion_segments]
- (const std::string& ion_name, size_type segment_idx, const cable_cell_local_parameter_set& seg_param, const ion_dependency* iondep = nullptr)
+ (const std::string& ion_name, size_type segment_idx, const cable_cell_parameter_set& seg_param, const ion_dependency* iondep = nullptr)
{
const auto& global_ion_data = gparam.ion_data;
const auto& cell_ion_data = cell.default_parameters.ion_data;
@@ -614,33 +617,40 @@ fvm_mechanism_data fvm_build_mechanism_data(const cable_cell_global_properties&
}
}
- for (const auto& cellsyn: cell.synapses()) {
- const mechanism_desc& desc = cellsyn.mechanism;
- size_type cv = D.branch_location_cv(cell_idx, cellsyn.location);
- const auto& name = desc.name();
-
+ for (const auto& cellsyn_by_name: cell.synapses()) {
+ const auto& name = cellsyn_by_name.first;
point_mech_data& entry = point_mech_table[name];
- update_paramset_and_validate(desc, entry.info, entry.paramset);
- entry.points.push_back({cv, target_id++, &desc});
- const segment_ptr& seg = cell.segments()[cellsyn.location.branch];
- size_type segment_idx = D.cell_segment_bounds[cell_idx]+cellsyn.location.branch;
+ for (const auto& placed_synapse: cellsyn_by_name.second) {
+ mlocation loc = placed_synapse.loc;
+ cell_lid_type target = target_offset + placed_synapse.lid;
+ ++target_count;
+ const mechanism_desc& desc = placed_synapse.item;
+
+ update_paramset_and_validate(desc, entry.info, entry.paramset);
+
+ size_type cv = D.branch_location_cv(cell_idx, loc);
+ entry.points.push_back({cv, target, &desc});
+
+ const segment_ptr& seg = cell.segments()[loc.branch];
+ size_type segment_idx = D.cell_segment_bounds[cell_idx]+loc.branch;
- for (const auto& ion_entry: entry.info->ions) {
- const std::string& ion_name = ion_entry.first;
- const ion_dependency& iondep = ion_entry.second;
+ for (const auto& ion_entry: entry.info->ions) {
+ const std::string& ion_name = ion_entry.first;
+ const ion_dependency& iondep = ion_entry.second;
- add_ion_segment(ion_name, segment_idx, seg->parameters, &iondep);
+ add_ion_segment(ion_name, segment_idx, seg->parameters, &iondep);
- if (ion_entry.second.read_reversal_potential) {
- ion_revpot_segments[ion_name][cell_idx].insert(segment_idx);
+ if (ion_entry.second.read_reversal_potential) {
+ ion_revpot_segments[ion_name][cell_idx].insert(segment_idx);
+ }
}
}
}
- for (const auto& stimulus: cell.stimuli()) {
- size_type cv = D.branch_location_cv(cell_idx, stimulus.location);
- stimuli.push_back({cv, stimulus.clamp});
+ for (const auto& loc_clamp: cell.stimuli()) {
+ size_type cv = D.branch_location_cv(cell_idx, loc_clamp.loc);
+ stimuli.push_back({cv, loc_clamp.item});
}
// Add segments to ion_segments map which intersect with existing segments, so
diff --git a/arbor/fvm_lowered_cell_impl.hpp b/arbor/fvm_lowered_cell_impl.hpp
index b8dc2f5a14fa688ebd12f8c67633aa1bd34628a6..32642764db8f1696963ce028a0d89a6c6f136246 100644
--- a/arbor/fvm_lowered_cell_impl.hpp
+++ b/arbor/fvm_lowered_cell_impl.hpp
@@ -516,9 +516,9 @@ void fvm_lowered_cell_impl<B>::initialize(
for (auto cell_idx: make_span(ncell)) {
cell_gid_type gid = gids[cell_idx];
- for (auto detector: cells[cell_idx].detectors()) {
- detector_cv.push_back(D.branch_location_cv(cell_idx, detector.location));
- detector_threshold.push_back(detector.threshold);
+ for (auto entry: cells[cell_idx].detectors()) {
+ detector_cv.push_back(D.branch_location_cv(cell_idx, entry.loc));
+ detector_threshold.push_back(entry.item.threshold);
}
for (cell_lid_type j: make_span(rec.num_probes(gid))) {
@@ -563,9 +563,9 @@ std::vector<fvm_gap_junction> fvm_lowered_cell_impl<B>::fvm_gap_junctions(
if (rec.num_gap_junction_sites(gids[cell_idx])) {
gid_to_cvs[gids[cell_idx]].reserve(rec.num_gap_junction_sites(gids[cell_idx]));
- auto cell_gj = cells[cell_idx].gap_junction_sites();
+ const auto& cell_gj = cells[cell_idx].gap_junction_sites();
for (auto gj : cell_gj) {
- auto cv = D.branch_location_cv(cell_idx, gj);
+ auto cv = D.branch_location_cv(cell_idx, gj.loc);
gid_to_cvs[gids[cell_idx]].push_back(cv);
}
}
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index 9745d922d2f4f4b5c4c2504571ef639a6a79d156..b87d5e6efa9c2211b542a61f0dafd13a68bd456d 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -1,6 +1,8 @@
#pragma once
#include <string>
+#include <unordered_map>
+#include <utility>
#include <vector>
#include <arbor/arbexcept.hpp>
@@ -9,10 +11,12 @@
#include <arbor/constants.hpp>
#include <arbor/mechcat.hpp>
#include <arbor/morph/label_dict.hpp>
+#include <arbor/morph/mcable_map.hpp>
#include <arbor/morph/mprovider.hpp>
#include <arbor/morph/morphology.hpp>
#include <arbor/morph/primitives.hpp>
#include <arbor/segment.hpp>
+#include <arbor/util/typed_map.hpp>
namespace arb {
@@ -39,6 +43,42 @@ struct cell_probe_address {
// Forward declare the implementation, for PIMPL.
struct cable_cell_impl;
+// Typed maps for access to painted and placed assignments:
+//
+// Mechanisms and initial ion data are further keyed by
+// mechanism name and ion name respectively.
+
+template <typename T>
+using region_assignment =
+ std::conditional_t<
+ std::is_same<T, mechanism_desc>::value || std::is_same<T, initial_ion_data>::value,
+ std::unordered_map<std::string, mcable_map<T>>,
+ mcable_map<T>>;
+
+template <typename T>
+struct placed {
+ mlocation loc;
+ cell_lid_type lid;
+ T item;
+};
+
+template <typename T>
+using mlocation_map = std::vector<placed<T>>;
+
+template <typename T>
+using location_assignment =
+ std::conditional_t<
+ std::is_same<T, mechanism_desc>::value,
+ std::unordered_map<std::string, mlocation_map<T>>,
+ mlocation_map<T>>;
+
+using cable_cell_region_map = static_typed_map<region_assignment,
+ mechanism_desc, init_membrane_potential, axial_resistivity,
+ temperature_K, membrane_capacitance, initial_ion_data>;
+
+using cable_cell_location_map = static_typed_map<location_assignment,
+ mechanism_desc, i_clamp, gap_junction_site, threshold_detector>;
+
// High-level abstract representation of a cell and its segments
class cable_cell {
public:
@@ -49,21 +89,6 @@ public:
using gap_junction_instance = mlocation;
- struct synapse_instance {
- mlocation location;
- mechanism_desc mechanism;
- };
-
- struct stimulus_instance {
- mlocation location;
- i_clamp clamp;
- };
-
- struct detector_instance {
- mlocation location;
- double threshold;
- };
-
cable_cell_parameter_set default_parameters;
/// Default constructor
@@ -88,6 +113,32 @@ public:
// the number of branches in the cell
size_type num_branches() const;
+ // Set cell-wide default physical and ion parameters.
+
+ void set_default(init_membrane_potential prop) {
+ default_parameters.init_membrane_potential = prop.value;
+ }
+
+ void set_default(axial_resistivity prop) {
+ default_parameters.axial_resistivity = prop.value;
+ }
+
+ void set_default(temperature_K prop) {
+ default_parameters.temperature_K = prop.value;
+ }
+
+ void set_default(membrane_capacitance prop) {
+ default_parameters.membrane_capacitance = prop.value;
+ }
+
+ void set_default(initial_ion_data prop) {
+ default_parameters.ion_data[prop.ion] = prop.initial;
+ }
+
+ void set_default(ion_reversal_potential_method prop) {
+ default_parameters.reversal_potential_method[prop.ion] = prop.method;
+ }
+
// All of the members marked with LEGACY below will be removed once
// the discretization code has moved from consuming segments to em_morphology.
@@ -130,26 +181,45 @@ public:
void paint(const region&, mechanism_desc);
// Properties.
- void paint(const region&, cable_cell_local_parameter_set);
+ void paint(const region&, init_membrane_potential);
+ void paint(const region&, axial_resistivity);
+ void paint(const region&, temperature_K);
+ void paint(const region&, membrane_capacitance);
+ void paint(const region&, initial_ion_data);
// Synapses.
- lid_range place(const locset&, const mechanism_desc&);
+ lid_range place(const locset&, mechanism_desc);
// Stimuli.
- lid_range place(const locset&, const i_clamp&);
+ lid_range place(const locset&, i_clamp);
// Gap junctions.
lid_range place(const locset&, gap_junction_site);
// Spike detectors.
- lid_range place(const locset&, const threshold_detector&);
+ lid_range place(const locset&, threshold_detector);
+
+ // Convenience access to placed items.
+
+ const std::unordered_map<std::string, mlocation_map<mechanism_desc>>& synapses() const {
+ return location_assignments().get<mechanism_desc>();
+ }
+
+ const mlocation_map<gap_junction_site>& gap_junction_sites() const {
+ return location_assignments().get<gap_junction_site>();
+ }
+
+ const mlocation_map<threshold_detector>& detectors() const {
+ return location_assignments().get<threshold_detector>();
+ }
- // Access to placed items.
+ const mlocation_map<i_clamp>& stimuli() const {
+ return location_assignments().get<i_clamp>();
+ }
- const std::vector<synapse_instance>& synapses() const;
- const std::vector<gap_junction_instance>& gap_junction_sites() const;
- const std::vector<detector_instance>& detectors() const;
- const std::vector<stimulus_instance>& stimuli() const;
+ // Generic access to painted and placed items.
+ const cable_cell_region_map& region_assignments() const;
+ const cable_cell_location_map& location_assignments() const;
// Checks that two cells have the same
// - number and type of segments
diff --git a/arbor/include/arbor/cable_cell_param.hpp b/arbor/include/arbor/cable_cell_param.hpp
index 92b1822cf2cd5b82ce1495f43373c913f31887ed..3cc08754792c3372cd70767500b1bc5d6bf789d4 100644
--- a/arbor/include/arbor/cable_cell_param.hpp
+++ b/arbor/include/arbor/cable_cell_param.hpp
@@ -18,6 +18,17 @@ struct cable_cell_error: arbor_exception {
arbor_exception("cable_cell: "+what) {}
};
+// Ion inital concentration and reversal potential
+// parameters, as used in cable_cell_parameter_set,
+// and set locally via painting initial_ion_data
+// (see below).
+
+struct cable_cell_ion_data {
+ double init_int_concentration = NAN;
+ double init_ext_concentration = NAN;
+ double init_reversal_potential = NAN;
+};
+
// Current clamp description for stimulus specification.
struct i_clamp {
using value_type = double;
@@ -41,6 +52,25 @@ struct threshold_detector {
// Tag type for dispatching cable_cell::place() calls that add gap junction sites.
struct gap_junction_site {};
+// Setter types for painting physical and ion parameters or setting
+// cell-wide default:
+
+struct init_membrane_potential {
+ double value = NAN; // [mV]
+};
+
+struct temperature_K {
+ double value = NAN; // [K]
+};
+
+struct axial_resistivity {
+ double value = NAN; // [[Ω·cm]
+};
+
+struct membrane_capacitance {
+ double value = NAN; // [F/m²]
+};
+
// Mechanism description, viz. mechanism name and
// (non-global) parameter settings. Used to assign
// density and point mechanisms to segments and
@@ -104,6 +134,16 @@ private:
std::unordered_map<std::string, double> param_;
};
+struct initial_ion_data {
+ std::string ion;
+ cable_cell_ion_data initial;
+};
+
+struct ion_reversal_potential_method {
+ std::string ion;
+ mechanism_desc method;
+};
+
// FVM discretization policies/hints.
//
// CV polices, given a cable cell, provide a locset comprising
@@ -236,45 +276,28 @@ inline cv_policy default_cv_policy() {
}
// Cable cell ion and electrical defaults.
-//
-// Parameters can be overridden with `cable_cell_local_parameter_set`
-// on unbranched segments within a cell; per-cell and global defaults
-// use `cable_cell_parameter_set`, which extends the parameter set
-// to supply per-cell or global ion reversal potential calculation
-// mechanisms.
-struct cable_cell_ion_data {
- double init_int_concentration = NAN;
- double init_ext_concentration = NAN;
- double init_reversal_potential = NAN;
-};
+// Parameters can be given as per-cell and global defaults via
+// cable_cell::default_parameters and cable_cell_global_properties::default_parameters
+// respectively.
+//
+// With the exception of `reversal_potential_method`, these properties can
+// be set locally witihin a cell using the `cable_cell::paint()`, and the
+// cell defaults can be individually set with `cable_cell:set_default()`.
-struct cable_cell_local_parameter_set {
- std::unordered_map<std::string, cable_cell_ion_data> ion_data;
+struct cable_cell_parameter_set {
util::optional<double> init_membrane_potential; // [mV]
util::optional<double> temperature_K; // [K]
util::optional<double> axial_resistivity; // [Ω·cm]
util::optional<double> membrane_capacitance; // [F/m²]
-};
-struct cable_cell_parameter_set: public cable_cell_local_parameter_set {
+ std::unordered_map<std::string, cable_cell_ion_data> ion_data;
std::unordered_map<std::string, mechanism_desc> reversal_potential_method;
- cv_policy discretization = default_cv_policy();
- // We'll need something like this until C++17, for sane initialization syntax.
- cable_cell_parameter_set() = default;
- cable_cell_parameter_set(
- cable_cell_local_parameter_set p,
- std::unordered_map<std::string, mechanism_desc> m = {},
- cv_policy d = default_cv_policy()
- ):
- cable_cell_local_parameter_set(std::move(p)),
- reversal_potential_method(std::move(m)),
- discretization(std::move(d))
- {}
+ cv_policy discretization = default_cv_policy();
};
-extern cable_cell_local_parameter_set neuron_parameter_defaults;
+extern cable_cell_parameter_set neuron_parameter_defaults;
// Global cable cell data.
diff --git a/arbor/include/arbor/morph/mcable_map.hpp b/arbor/include/arbor/morph/mcable_map.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..94b638b1e9e73db56b3440a057acca364170efa7
--- /dev/null
+++ b/arbor/include/arbor/morph/mcable_map.hpp
@@ -0,0 +1,122 @@
+#pragma once
+
+// Represent an mcable_list where each cable has an associated value.
+//
+// The only mutating operations are insert, emplace, and clear.
+
+#include <algorithm>
+#include <vector>
+
+#include <arbor/assert.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/util/optional.hpp>
+
+namespace arb {
+
+template <typename T>
+struct mcable_map {
+ // Forward subset of vector interface:
+
+ using value_type = std::pair<mcable, T>;
+ using size_type = std::size_t;
+ using difference_type = std::ptrdiff_t;
+
+ using const_iterator = typename std::vector<value_type>::const_iterator;
+ using iterator = const_iterator;
+
+ using const_reverse_iterator = typename std::vector<value_type>::const_reverse_iterator;
+ using reverse_iterator = const_reverse_iterator;
+
+ const_iterator begin() const noexcept { return elements_.cbegin(); }
+ const_iterator cbegin() const noexcept { return elements_.cbegin(); }
+
+ const_reverse_iterator rbegin() const noexcept { return elements_.crbegin(); }
+ const_reverse_iterator crbegin() const noexcept { return elements_.crbegin(); }
+
+ const_iterator end() const noexcept { return elements_.cend(); }
+ const_iterator cend() const noexcept { return elements_.cend(); }
+
+ const_reverse_iterator rend() const noexcept { return elements_.crend(); }
+ const_reverse_iterator crend() const noexcept { return elements_.crend(); }
+
+ const value_type& operator[](size_type i) const { return elements_[i]; }
+
+ decltype(auto) front() const { return *begin(); }
+ decltype(auto) back() const { return *rbegin(); }
+
+ std::size_t size() const noexcept { return elements_.size(); }
+ bool empty() const noexcept { return !size(); }
+
+ void clear() { elements_.clear(); }
+
+ // mcable_map-specific operations:
+
+ // Insertion is successful iff c intersects with cables in at most two discrete points.
+ // insert() and emplace() return true on success.
+
+ bool insert(const mcable& c, T value) {
+ auto opt_it = insertion_point(c);
+ if (!opt_it) return false;
+
+ elements_.insert(*opt_it, value_type{c, std::move(value)});
+ assert_invariants();
+ return true;
+ }
+
+ template <typename... Args>
+ bool emplace(const mcable& c, Args&&... args) {
+ auto opt_it = insertion_point(c);
+ if (!opt_it) return false;
+
+ elements_.emplace(*opt_it,
+ std::piecewise_construct,
+ std::forward_as_tuple(c),
+ std::forward_as_tuple(std::forward<Args>(args)...));
+
+ assert_invariants();
+ return true;
+ }
+
+ mcable_list support() const {
+ mcable_list s;
+ s.reserve(elements_.size());
+ std::transform(elements_.begin(), elements_.end(), std::back_inserter(s),
+ [](const auto& x) { return x.first; });
+
+ return s;
+ }
+
+private:
+ std::vector<value_type> elements_;
+
+ util::optional<typename std::vector<value_type>::iterator> insertion_point(const mcable& c) {
+ struct as_mcable {
+ mcable value;
+ as_mcable(const value_type& x): value(x.first) {}
+ as_mcable(const mcable& x): value(x) {}
+ };
+
+ auto it = std::lower_bound(elements_.begin(), elements_.end(), c,
+ [](as_mcable a, as_mcable b) { return a.value<b.value; });
+
+ if (it!=elements_.begin()) {
+ mcable prior = std::prev(it)->first;
+ if (prior.branch==c.branch && prior.dist_pos>c.prox_pos) {
+ return util::nullopt;
+ }
+ }
+ if (it!=elements_.end()) {
+ mcable next = it->first;
+ if (c.branch==next.branch && c.dist_pos>next.prox_pos) {
+ return util::nullopt;
+ }
+ }
+ return it;
+ }
+
+ void assert_invariants() {
+ arb_assert(test_invariants(support()));
+ }
+};
+
+} // namespace arb
diff --git a/arbor/include/arbor/segment.hpp b/arbor/include/arbor/segment.hpp
index 10264eb525079b5ba51f78f3f19feee77bbb669d..b87bcf36e8a54bb65d0c7fb41895e53c4f7f6031 100644
--- a/arbor/include/arbor/segment.hpp
+++ b/arbor/include/arbor/segment.hpp
@@ -113,7 +113,7 @@ public:
return mechanisms_;
}
- cable_cell_local_parameter_set parameters; // override cell and global defaults
+ cable_cell_parameter_set parameters; // override cell and global defaults
protected:
segment(section_kind kind): kind_(kind) {}
diff --git a/arbor/include/arbor/util/typed_map.hpp b/arbor/include/arbor/util/typed_map.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..785602e1e9e5352c190a85d312aee0260ca0c2ea
--- /dev/null
+++ b/arbor/include/arbor/util/typed_map.hpp
@@ -0,0 +1,83 @@
+#pragma once
+
+// Maps keyed by type, with a fixed set of keys (static_typed_map) or
+// an arbitrary set (dynamic_typed_map).
+//
+// The first template parameter maps the key type to its corresponding
+// value type.
+//
+// Example: associating a vector of ints with 'int' and a vector of
+// doubles with 'double':
+//
+// static_typed_map<std::vector, int, double> m;
+// m.get<int>() = {1, 2, 3};
+// m.get<double>() = {1.2, 2.3};
+
+#include <tuple>
+#include <typeindex>
+#include <unordered_map>
+
+#include <arbor/util/any.hpp>
+
+namespace arb {
+
+template <template <class> typename E>
+struct dynamic_typed_map {
+ // Retrieve value by reference associated with type T; create entry with
+ // default value if no entry in map for T.
+ template <typename T>
+ E<T>& get() {
+ arb::util::any& store_entry = tmap_[std::type_index(typeid(T))];
+ if (!store_entry.has_value()) {
+ store_entry = arb::util::any(E<T>{});
+ }
+
+ return arb::util::any_cast<E<T>&>(store_entry);
+ }
+
+ // Retrieve value by const reference associated with type T;
+ // throw if no entry in map for T.
+ template <typename T>
+ const E<T>& get() const {
+ return arb::util::any_cast<const E<T>&>(tmap_.at(std::type_index(typeid(T))));
+ }
+
+ // True if map has an entry for type T.
+ template <typename T>
+ bool has() const { return tmap_.count(std::type_index(typeid(T))); }
+
+private:
+ std::unordered_map<std::type_index, arb::util::any> tmap_;
+};
+
+template <template <class> typename E, typename... Keys>
+struct static_typed_map {
+ // Retrieve value by reference associated with type T.
+ template <typename T>
+ E<T>& get() {
+ return std::get<index<T, Keys...>()>(tmap_);
+ }
+
+ // Retrieve value by const reference associated with type T.
+ template <typename T>
+ const E<T>& get() const {
+ return std::get<index<T, Keys...>()>(tmap_);
+ }
+
+ // True if map has an entry for type T.
+ template <typename T>
+ constexpr bool has() const { return index<T, Keys...>()<sizeof...(Keys); }
+
+private:
+ std::tuple<E<Keys>...> tmap_;
+
+ template <typename T>
+ static constexpr int index() { return 1; }
+
+ template <typename T, typename H, typename... A>
+ static constexpr int index() {
+ return std::is_same<H, T>::value? 0: 1+index<T, A...>();
+ }
+};
+
+} // namespace arb
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 8deea24dcce1ab01983d8a176b770427ad338faa..07d3d52fcca410e6dd3fc76ce7b9336f2fd9f5b1 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -105,6 +105,7 @@ set(unit_sources
test_mask_stream.cpp
test_math.cpp
test_matrix.cpp
+ test_mcable_map.cpp
test_mc_cell_group.cpp
test_mechanisms.cpp
test_mech_temp_diam.cpp
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index e1a0a5f82ea84b4bba3156d15955b8765d53bc62..b27a7b6df5948d7991da8de8bf0df2ec419b3045 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -101,17 +101,10 @@ namespace {
cell.paint("soma", "hh");
cell.paint("dend", "pas");
- cable_cell_local_parameter_set p1;
- p1.membrane_capacitance = 0.017;
- cable_cell_local_parameter_set p2;
- p2.membrane_capacitance = 0.013;
- cable_cell_local_parameter_set p3;
- p3.membrane_capacitance = 0.018;
-
using ::arb::reg::branch;
- cell.paint(branch(b1), p1);
- cell.paint(branch(b2), p2);
- cell.paint(branch(b3), p3);
+ cell.paint(branch(b1), membrane_capacitance{0.017});
+ cell.paint(branch(b2), membrane_capacitance{0.013});
+ cell.paint(branch(b3), membrane_capacitance{0.018});
cell.place(mlocation{2,1}, i_clamp{5., 80., 0.45});
cell.place(mlocation{3,1}, i_clamp{40., 10.,-0.2});
diff --git a/test/unit/test_mcable_map.cpp b/test/unit/test_mcable_map.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f5eb215c3b312316f58457c4f2ff9eecaefd07e4
--- /dev/null
+++ b/test/unit/test_mcable_map.cpp
@@ -0,0 +1,109 @@
+#include "../test/gtest.h"
+
+#include <arbor/morph/mcable_map.hpp>
+#include <arbor/morph/primitives.hpp>
+
+#include "util/span.hpp"
+
+using namespace arb;
+
+TEST(mcable_map, insertion) {
+ mcable_list cl{{0, 0.2, 0.8}, {1, 0, 1}, {2, 0.2, 0.4}, {2, 0.4, 0.7}, {2, 0.8, 1.}, {4, 0.2, 0.4}};
+ mcable_map<int> mm;
+ bool ok = false;
+
+ int j = 10;
+ for (const mcable& c: cl) {
+ ok = mm.insert(c, j++);
+ EXPECT_TRUE(ok);
+ }
+ auto expected_size = cl.size();
+ EXPECT_EQ(expected_size, mm.size());
+
+ // insertions at front or in middle:
+ ok = mm.insert(mcable{0, 0.1, 0.2}, 1);
+ EXPECT_TRUE(ok);
+ ++expected_size;
+ EXPECT_EQ(expected_size, mm.size());
+
+ ok = mm.insert(mcable{0, 0.05, 0.2}, 1);
+ EXPECT_FALSE(ok); // overlap
+ EXPECT_EQ(expected_size, mm.size());
+
+ ok = mm.insert(mcable{1, 1., 1.}, 1);
+ EXPECT_TRUE(ok); // overlapping one point is fine
+ ++expected_size;
+ EXPECT_EQ(expected_size, mm.size());
+
+ ok = mm.insert(mcable{3, 0., 1.}, 1);
+ EXPECT_TRUE(ok);
+ ++expected_size;
+ EXPECT_EQ(expected_size, mm.size());
+
+ ok = mm.insert(mcable{2, 0.7, 0.8}, 1);
+ EXPECT_TRUE(ok);
+ ++expected_size;
+ EXPECT_EQ(expected_size, mm.size());
+
+ ok = mm.insert(mcable{2, 0., 0.3}, 1);
+ EXPECT_FALSE(ok); // overlap
+ EXPECT_EQ(expected_size, mm.size());
+}
+
+TEST(mcable_map, emplace) {
+ struct foo {
+ foo(char a, int b): sa(a+1), sb(b+2) {}
+ foo(const foo&) = delete;
+ foo(foo&&) = default;
+
+ foo& operator=(const foo&) = delete;
+ foo& operator=(foo&&) = default;
+
+ char sa;
+ int sb;
+ };
+
+ mcable_map<foo> mm;
+ mm.emplace(mcable{1, 0.2, 0.4}, 'a', 3);
+ mm.emplace(mcable{2, 0.2, 0.4}, 'x', 5);
+ ASSERT_EQ(2u, mm.size());
+
+ EXPECT_EQ('b', mm[0].second.sa);
+ EXPECT_EQ(5, mm[0].second.sb);
+
+ EXPECT_EQ('y', mm[1].second.sa);
+ EXPECT_EQ(7, mm[1].second.sb);
+}
+
+TEST(mcable_map, access) {
+ mcable_list cl{{0, 0.2, 0.8}, {1, 0, 1}, {2, 0.2, 0.4}, {2, 0.4, 0.7}, {2, 0.8, 1.}, {4, 0.2, 0.4}};
+ mcable_map<int> mm;
+ bool ok = false;
+
+ int k = 10;
+ for (const mcable& c: cl) {
+ ok = mm.insert(c, k++);
+ ASSERT_TRUE(ok);
+ }
+ auto size = cl.size();
+ ASSERT_EQ(size, mm.size());
+
+ for (std::size_t i = 0; i<size; ++i) {
+ EXPECT_EQ(cl[i], mm[i].first);
+ EXPECT_EQ(int(10+i), mm[i].second);
+ }
+
+ std::size_t j = 0;
+ for (const auto& el: mm) {
+ EXPECT_EQ(cl[j], el.first);
+ EXPECT_EQ(int(10+j), el.second);
+ ++j;
+ }
+
+ j = size;
+ for (auto r = mm.rbegin(); r!=mm.rend(); ++r) {
+ --j;
+ EXPECT_EQ(cl[j], r->first);
+ EXPECT_EQ(int(10+j), r->second);
+ }
+}
diff --git a/test/unit/test_synapses.cpp b/test/unit/test_synapses.cpp
index 0c10474c08d58fe8d2bd7b3af51c3f6f800fc52f..de1a796b88763ffe33e31b3dc96c84469d43516e 100644
--- a/test/unit/test_synapses.cpp
+++ b/test/unit/test_synapses.cpp
@@ -39,20 +39,19 @@ TEST(synapses, add_to_cell) {
cell.place(mlocation{0, 0.2}, "exp2syn");
cell.place(mlocation{0, 0.3}, "expsyn");
- EXPECT_EQ(3u, cell.synapses().size());
- const auto& syns = cell.synapses();
+ auto syns = cell.synapses();
- EXPECT_EQ(syns[0].location.branch, 0u);
- EXPECT_EQ(syns[0].location.pos, 0.1);
- EXPECT_EQ(syns[0].mechanism.name(), "expsyn");
+ ASSERT_EQ(2u, syns["expsyn"].size());
+ ASSERT_EQ(1u, syns["exp2syn"].size());
- EXPECT_EQ(syns[1].location.branch, 0u);
- EXPECT_EQ(syns[1].location.pos, 0.2);
- EXPECT_EQ(syns[1].mechanism.name(), "exp2syn");
+ EXPECT_EQ((mlocation{0, 0.1}), syns["expsyn"][0].loc);
+ EXPECT_EQ("expsyn", syns["expsyn"][0].item.name());
- EXPECT_EQ(syns[2].location.branch, 0u);
- EXPECT_EQ(syns[2].location.pos, 0.3);
- EXPECT_EQ(syns[2].mechanism.name(), "expsyn");
+ EXPECT_EQ((mlocation{0, 0.3}), syns["expsyn"][1].loc);
+ EXPECT_EQ("expsyn", syns["expsyn"][1].item.name());
+
+ EXPECT_EQ((mlocation{0, 0.2}), syns["exp2syn"][0].loc);
+ EXPECT_EQ("exp2syn", syns["exp2syn"][0].item.name());
// adding a synapse to an invalid branch location should throw.
EXPECT_THROW(cell.place(mlocation{1, 0.3}, "expsyn"), std::runtime_error);