diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0ffce50ffc32f9e136adc25ef1fb4309667eab04..c89d2d979e30142e96da464c9d9d46ec9aec9f2e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -25,15 +25,6 @@ option(ARB_VECTORIZE "use explicit SIMD code in generated mechanisms" OFF)
set(ARB_MODCC "" CACHE STRING "path to external modcc NMODL compiler")
-# Generate validation data for validation tests?
-
-option(ARB_BUILD_VALIDATION_DATA "generate validation data" OFF)
-
-# Where to generate and find validation data?
-
-set(ARB_VALIDATION_DATA_DIR "${PROJECT_SOURCE_DIR}/validation/data" CACHE PATH
- "location of generated validation data")
-
# For sup::glop, just wrap POSIX glob(3)?
# Turn off for platforms without glob(3) in libc, e.g. Android Bionic.
@@ -372,7 +363,7 @@ add_subdirectory(arbor)
# arborenv, arborenv-public-headers:
add_subdirectory(arborenv)
-# unit, unit-mpi, unit-local, unit-modcc, validate
+# unit, unit-mpi, unit-local, unit-modcc
add_subdirectory(test)
# self contained examples:
@@ -381,11 +372,6 @@ add_subdirectory(example)
# html:
add_subdirectory(doc)
-# validation-data:
-if(ARB_BUILD_VALIDATION_DATA)
- add_subdirectory(validation)
-endif()
-
# python interface:
if (ARB_WITH_PYTHON)
add_subdirectory(python)
diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp
index affac000ba94f30a01317cbd23a21eaecd3c661d..bf41a4b53111ad44ebd3b2a482de6aea3d8b9f3e 100644
--- a/arbor/cable_cell.cpp
+++ b/arbor/cable_cell.cpp
@@ -30,11 +30,85 @@ struct cable_cell_impl {
using region_map = cable_cell::region_map;
using locset_map = cable_cell::locset_map;
- cable_cell_impl() {
- segments.push_back(make_segment<placeholder_segment>());
- parents.push_back(0);
+ cable_cell_impl(const arb::morphology& m,
+ const label_dict& dictionary,
+ bool compartments_from_discretization):
+ morph(m)
+ {
+ using point = cable_cell::point_type;
+ if (!m.num_branches()) {
+ segments.push_back(make_segment<placeholder_segment>());
+ parents.push_back(0);
+ return;
+ }
+
+ // Add the soma.
+ auto loc = m.samples()[0].loc; // location of soma.
+
+ // If there is no spherical root/soma use a zero-radius soma.
+ double srad = m.spherical_root()? loc.radius: 0.;
+ segments.push_back(make_segment<soma_segment>(srad, point(loc.x, loc.y, loc.z)));
+ parents.push_back(-1);
+
+ auto& samples = m.samples();
+ auto& props = m.sample_props();
+ for (auto i: util::make_span(1, m.num_branches())) {
+ auto index = util::make_range(m.branch_indexes(i));
+
+ // find kind for the branch. Use the tag of the last sample in the branch.
+ int tag = samples[index.back()].tag;
+ section_kind kind;
+ switch (tag) {
+ case 1: // soma
+ throw cable_cell_error("No support for complex somata (yet)");
+ case 2: // axon
+ kind = section_kind::axon;
+ case 3: // dendrite
+ case 4: // apical dendrite
+ default: // just take dendrite as default
+ kind = section_kind::dendrite;
+ }
+
+ std::vector<value_type> radii;
+ std::vector<cable_cell::point_type> points;
+
+ // The current discretization code does not handle collocated points correctly,
+ // particularly if they lie at the start of a branch, so we have to skip the first
+ // point on a branch if it is collocated with the second point.
+ bool skip_first = is_collocated(props[index[1]]);
+ for (auto j: util::make_span(skip_first, index.size())) {
+ auto& s = samples[index[j]];
+ radii.push_back(s.loc.radius);
+ points.push_back(cable_cell::point_type(s.loc.x, s.loc.y, s.loc.z));
+ }
+
+ // Find the id of this branch's parent.
+ auto pid = m.branch_parent(i);
+ // Adjust pid if a zero-radius soma was used.
+ if (!m.spherical_root()) {
+ pid = pid==mnpos? 0: pid+1;
+ }
+ segments.push_back(make_segment<cable_segment>(kind, radii, points));
+ parents.push_back(pid);
+ if (compartments_from_discretization) {
+ int ncolloc = std::count_if(index.begin(), index.end(), [&props](auto i){return is_collocated(props[i]);});
+ int ncomp = index.size()-ncolloc-1;
+ ncomp -= is_collocated(props[index[0]]);
+ segments.back()->as_cable()->set_compartments(ncomp);
+ }
+ }
+
+ for (auto& r: dictionary.regions()) {
+ regions[r.first] = thingify(r.second, morph);
+ }
+
+ for (auto& l: dictionary.locsets()) {
+ locations[l.first] = thingify(l.second, morph);
+ }
}
+ cable_cell_impl(): cable_cell_impl({},{},false) {}
+
cable_cell_impl(const cable_cell_impl& other) {
parents = other.parents;
stimuli = other.stimuli;
@@ -93,6 +167,16 @@ struct cable_cell_impl {
return lid_range(first, list.size());
}
+ template <typename Desc, typename T>
+ lid_range place(const std::string& target, const Desc& desc, std::vector<T>& list) {
+ const auto first = list.size();
+
+ const auto it = locations.find(target);
+ if (it==locations.end()) return lid_range(first, first);
+
+ return place(it->second, desc, list);
+ }
+
lid_range place_gj(const mlocation_list& locs) {
const auto first = gap_junction_sites.size();
@@ -101,6 +185,15 @@ struct cable_cell_impl {
return lid_range(first, gap_junction_sites.size());
}
+ lid_range place_gj(const std::string& target) {
+ const auto first = gap_junction_sites.size();
+
+ const auto it = locations.find(target);
+ if (it==locations.end()) return lid_range(first, first);
+
+ return place_gj(it->second);
+ }
+
void assert_valid_segment(index_type i) const {
if (i>=segments.size()) {
throw cable_cell_error("no such segment");
@@ -110,6 +203,28 @@ struct cable_cell_impl {
bool valid_location(const mlocation& loc) const {
return test_invariants(loc) && loc.branch<segments.size();
}
+
+ template <typename F>
+ void paint(const std::string& target, F&& f) {
+ auto it = regions.find(target);
+
+ // Nothing to do if there are no regions that match.
+ if (it==regions.end()) return;
+
+ paint(it->second, std::forward<F>(f));
+ }
+
+ template <typename F>
+ void paint(const mcable_list& cables, F&& f) {
+ for (auto c: cables) {
+ 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]);
+ }
+ }
};
using impl_ptr = std::unique_ptr<cable_cell_impl, void (*)(cable_cell_impl*)>;
@@ -117,6 +232,12 @@ impl_ptr make_impl(cable_cell_impl* c) {
return impl_ptr(c, [](cable_cell_impl* p){delete p;});
}
+cable_cell::cable_cell(const arb::morphology& m,
+ const label_dict& dictionary,
+ bool compartments_from_discretization):
+ impl_(make_impl(new cable_cell_impl(m, dictionary, compartments_from_discretization)))
+{}
+
cable_cell::cable_cell():
impl_(make_impl(new cable_cell_impl()))
{}
@@ -126,42 +247,10 @@ cable_cell::cable_cell(const cable_cell& other):
impl_(make_impl(new cable_cell_impl(*other.impl_)))
{}
-size_type cable_cell::num_segments() const {
+size_type cable_cell::num_branches() const {
return impl_->segments.size();
}
-//
-// note: I think that we have to enforce that the soma is the first
-// segment that is added
-//
-soma_segment* cable_cell::add_soma(value_type radius, point_type center) {
- if (has_soma()) {
- throw cable_cell_error("cell already has soma");
- }
- impl_->segments[0] = make_segment<soma_segment>(radius, center);
- return impl_->segments[0]->as_soma();
-}
-
-cable_segment* cable_cell::add_cable(index_type parent, segment_ptr&& cable) {
- if (!cable->as_cable()) {
- throw cable_cell_error("segment is not a cable segment");
- }
-
- if (parent>num_segments()) {
- throw cable_cell_error("parent index out of range");
- }
-
- impl_->segments.push_back(std::move(cable));
- impl_->parents.push_back(parent);
-
- return impl_->segments.back()->as_cable();
-}
-
-segment* cable_cell::segment(index_type index) {
- impl_->assert_valid_segment(index);
- return impl_->segments[index].get();
-}
-
segment const* cable_cell::parent(index_type index) const {
impl_->assert_valid_segment(index);
return impl_->segments[impl_->parents[index]].get();
@@ -206,12 +295,8 @@ const std::vector<cable_cell::stimulus_instance>& cable_cell::stimuli() const {
return impl_->stimuli;
}
-void cable_cell::set_regions(cable_cell::region_map r) {
- impl_->regions = std::move(r);
-}
-
-void cable_cell::set_locsets(cable_cell::locset_map l) {
- impl_->locations = std::move(l);
+const em_morphology* cable_cell::morphology() const {
+ return &(impl_->morph);
}
//
@@ -221,19 +306,23 @@ void cable_cell::set_locsets(cable_cell::locset_map l) {
//
void cable_cell::paint(const std::string& target, mechanism_desc desc) {
- auto it = impl_->regions.find(target);
+ impl_->paint(target,
+ [&desc](segment_ptr& s){return s->add_mechanism(desc);});
+}
- // Nothing to do if there are no regions that match.
- if (it==impl_->regions.end()) return;
+void cable_cell::paint(const std::string& target, cable_cell_local_parameter_set params) {
+ impl_->paint(target,
+ [¶ms](segment_ptr& s){return s->parameters = params;});
+}
- for (auto c: it->second) {
- if (c.prox_pos!=0 || c.dist_pos!=1) {
- throw cable_cell_error(util::pprintf(
- "cable_cell does not support regions with partial branches: \"{}\": {}",
- target, c));
- }
- segment(c.branch)->add_mechanism(desc);
- }
+void cable_cell::paint(const region& target, mechanism_desc desc) {
+ impl_->paint(thingify(target, impl_->morph),
+ [&desc](segment_ptr& s){return s->add_mechanism(desc);});
+}
+
+void cable_cell::paint(const region& target, cable_cell_local_parameter_set params) {
+ impl_->paint(thingify(target, impl_->morph),
+ [¶ms](segment_ptr& s){return s->parameters = params;});
}
//
@@ -244,59 +333,27 @@ void cable_cell::paint(const std::string& target, mechanism_desc desc) {
//
//
-// Synapses
+// Synapses.
//
lid_range cable_cell::place(const std::string& target, const mechanism_desc& desc) {
- const auto first = impl_->synapses.size();
-
- const auto it = impl_->locations.find(target);
- if (it==impl_->locations.end()) return lid_range(first, first);
-
- return impl_->place(it->second, desc, impl_->synapses);
+ return impl_->place(target, desc, impl_->synapses);
}
-/*
lid_range cable_cell::place(const locset& ls, const mechanism_desc& desc) {
- const auto locs = thingify(ls, impl_->morph);
- return impl_->place(locs, desc, impl_->synapses);
-}
-*/
-
-lid_range cable_cell::place(const mlocation& loc, const mechanism_desc& desc) {
- if (!impl_->valid_location(loc)) {
- throw cable_cell_error(util::pprintf(
- "Attempt to add synapse at invalid location: \"{}\"", loc));
- }
- return impl_->place({loc}, desc, impl_->synapses);
+ return impl_->place(thingify(ls, impl_->morph), desc, impl_->synapses);
}
//
-// Stimuli
+// Stimuli.
//
lid_range cable_cell::place(const std::string& target, const i_clamp& desc) {
- const auto first = impl_->stimuli.size();
-
- const auto it = impl_->locations.find(target);
- if (it==impl_->locations.end()) return lid_range(first, first);
-
- return impl_->place(it->second, desc, impl_->stimuli);
+ return impl_->place(target, desc, impl_->stimuli);
}
-/*
lid_range cable_cell::place(const locset& ls, const i_clamp& desc) {
- const auto locs = thingify(ls, impl_->morph);
- return impl_->place(locs, desc, impl_->stimuli);
-}
-*/
-
-lid_range cable_cell::place(const mlocation& loc, const i_clamp& desc) {
- if (!impl_->valid_location(loc)) {
- throw cable_cell_error(util::pprintf(
- "Attempt to add stimulus at invalid location: {}", loc));
- }
- return impl_->place({loc}, desc, impl_->stimuli);
+ return impl_->place(thingify(ls, impl_->morph), desc, impl_->stimuli);
}
//
@@ -304,71 +361,31 @@ lid_range cable_cell::place(const mlocation& loc, const i_clamp& desc) {
//
lid_range cable_cell::place(const std::string& target, gap_junction_site) {
- const auto first = impl_->stimuli.size();
-
- const auto it = impl_->locations.find(target);
- if (it==impl_->locations.end()) return lid_range(first, first);
-
- return impl_->place_gj(it->second);
+ return impl_->place_gj(target);
}
-/*
lid_range cable_cell::place(const locset& ls, gap_junction_site) {
- const auto locs = thingify(ls, impl_->morph);
- return impl_->place_gj(locs);
-}
-*/
-
-lid_range cable_cell::place(const mlocation& loc, gap_junction_site) {
- if (!impl_->valid_location(loc)) {
- throw cable_cell_error(util::pprintf(
- "Attempt to add gap junction site at invalid location: {}", loc));
- }
- return impl_->place_gj({loc});
+ return impl_->place_gj(thingify(ls, impl_->morph));
}
//
// Spike detectors.
//
lid_range cable_cell::place(const std::string& target, const threshold_detector& desc) {
- const auto first = impl_->stimuli.size();
-
- const auto it = impl_->locations.find(target);
- if (it==impl_->locations.end()) return lid_range(first, first);
-
- return impl_->place(it->second, desc.threshold, impl_->spike_detectors);
+ return impl_->place(target, desc.threshold, impl_->spike_detectors);
}
-/*
lid_range cable_cell::place(const locset& ls, const threshold_detector& desc) {
- const auto locs = thingify(ls, impl_->morph);
- return impl_->place(locs, desc.threshold, impl_->spike_detectors);
-}
-*/
-
-lid_range cable_cell::place(const mlocation& loc, const threshold_detector& desc) {
- if (!impl_->valid_location(loc)) {
- throw cable_cell_error(util::pprintf(
- "Attempt to add spike detector at invalid location: {}", loc));
- }
- return impl_->place({loc}, desc.threshold, impl_->spike_detectors);
-}
-
-
-soma_segment* cable_cell::soma() {
- return has_soma()? segment(0)->as_soma(): nullptr;
+ return impl_->place(thingify(ls, impl_->morph), desc.threshold, impl_->spike_detectors);
}
+//
+// TODO: deprectate the following as soon as discretization code catches up with em_morphology
+//
const soma_segment* cable_cell::soma() const {
return has_soma()? segment(0)->as_soma(): nullptr;
}
-cable_segment* cable_cell::cable(index_type index) {
- impl_->assert_valid_segment(index);
- auto cable = segment(index)->as_cable();
- return cable? cable: throw cable_cell_error("segment is not a cable segment");
-}
-
const cable_segment* cable_cell::cable(index_type index) const {
impl_->assert_valid_segment(index);
auto cable = segment(index)->as_cable();
@@ -377,21 +394,13 @@ const cable_segment* cable_cell::cable(index_type index) const {
std::vector<size_type> cable_cell::compartment_counts() const {
std::vector<size_type> comp_count;
- comp_count.reserve(num_segments());
+ comp_count.reserve(num_branches());
for (const auto& s: segments()) {
comp_count.push_back(s->num_compartments());
}
return comp_count;
}
-const em_morphology* cable_cell::morphology() const {
- return &(impl_->morph);
-}
-
-void cable_cell::set_morphology(em_morphology m) {
- impl_->morph = std::move(m);
-}
-
size_type cable_cell::num_compartments() const {
return util::sum_by(impl_->segments,
[](const segment_ptr& s) { return s->num_compartments(); });
@@ -448,82 +457,4 @@ value_type cable_cell::segment_mean_attenuation(
return 2*std::sqrt(math::pi<double>*tau_per_um*frequency)*length_factor; // [1/µm]
}
-cable_cell make_cable_cell(const morphology& m,
- const label_dict& dictionary,
- bool compartments_from_discretization)
-{
- using point3d = cable_cell::point_type;
- cable_cell cell;
-
- if (!m.num_branches()) {
- return cell;
- }
-
- // Add the soma.
- auto loc = m.samples()[0].loc; // location of soma.
-
- // If there is no spherical root/soma use a zero-radius soma.
- double srad = m.spherical_root()? loc.radius: 0.;
- cell.add_soma(srad, point3d(loc.x, loc.y, loc.z));
-
- auto& samples = m.samples();
- for (auto i: util::make_span(1, m.num_branches())) {
- auto index = util::make_range(m.branch_indexes(i));
-
- // find kind for the branch. Use the tag of the last sample in the branch.
- int tag = samples[index.back()].tag;
- section_kind kind;
- switch (tag) {
- case 1: // soma
- throw cable_cell_error("No support for complex somata (yet)");
- case 2: // axon
- kind = section_kind::axon;
- case 3: // dendrite
- case 4: // apical dendrite
- default: // just take dendrite as default
- kind = section_kind::dendrite;
- }
-
- std::vector<value_type> radii;
- std::vector<point3d> points;
-
- for (auto i: index) {
- auto& s = samples[i];
- radii.push_back(s.loc.radius);
- points.push_back(point3d(s.loc.x, s.loc.y, s.loc.z));
- }
-
- // Find the id of this branch's parent.
- auto pid = m.branch_parent(i);
- // Adjust pid if a zero-radius soma was used.
- if (!m.spherical_root()) {
- pid = pid==mnpos? 0: pid+1;
- }
- auto cable = cell.add_cable(pid, make_segment<cable_segment>(kind, radii, points));
- if (compartments_from_discretization) {
- cable->as_cable()->set_compartments(radii.size()-1);
- }
- }
-
- // Construct concrete regions.
- // Ignores the pointsets, for now.
- auto em = em_morphology(m); // for converting labels to "segments"
-
- std::unordered_map<std::string, mcable_list> regions;
- for (auto r: dictionary.regions()) {
- regions[r.first] = thingify(r.second, em);
- }
- cell.set_regions(std::move(regions));
-
- std::unordered_map<std::string, mlocation_list> locsets;
- for (auto l: dictionary.locsets()) {
- locsets[l.first] = thingify(l.second, em);
- }
- cell.set_locsets(std::move(locsets));
-
- cell.set_morphology(std::move(em));
-
- return cell;
-}
-
} // namespace arb
diff --git a/arbor/fvm_layout.cpp b/arbor/fvm_layout.cpp
index fd293f8434685cef4e49a256a44edf4bdfa52f36..b10449e17b41a2b59124cfe3a096f4e8135668c5 100644
--- a/arbor/fvm_layout.cpp
+++ b/arbor/fvm_layout.cpp
@@ -210,7 +210,7 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells, const ca
fvm_discretization D;
util::make_partition(D.cell_segment_bounds,
- transform_view(cells, [](const cable_cell& c) { return c.num_segments(); }));
+ transform_view(cells, [](const cable_cell& c) { return c.num_branches(); }));
std::vector<index_type> cell_cv_bounds;
auto cell_cv_part = make_partition(cell_cv_bounds,
@@ -261,7 +261,7 @@ fvm_discretization fvm_discretize(const std::vector<cable_cell>& cells, const ca
seg_cv_bounds.clear();
auto seg_cv_part = make_partition(
seg_cv_bounds,
- transform_view(make_span(c.num_segments()), [&c](const unsigned s) {
+ transform_view(make_span(c.num_branches()), [&c](const unsigned s) {
if (!c.segment(s)->is_soma() && c.parent(s)->is_soma()) {
return c.segment(s)->num_compartments() + 1;
}
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index 2e35ded6c276c9a441c756b71b1134cc26d3b034..d6c8ff4f2129d7fdbac954f74ba832a58a8bfb90 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -79,39 +79,39 @@ public:
/// Move constructor
cable_cell(cable_cell&& other) = default;
- /// add a soma to the cell
- /// radius must be specified
- soma_segment* add_soma(value_type radius, point_type center=point_type());
+ /// construct from morphology
+ cable_cell(const class morphology& m,
+ const label_dict& dictionary={},
+ bool compartments_from_discretization=false);
- /// add a cable
- /// parent is the index of the parent segment for the cable section
- /// cable is the segment that will be moved into the cell
- cable_segment* add_cable(index_type parent, segment_ptr&& cable);
+ // the number of branches in the cell
+ size_type num_branches() const;
+ // All of the members marked with LEGACY below will be removed once
+ // the discretization code has moved from consuming segments to em_morphology.
+
+ // LEGACY
bool has_soma() const;
- class segment* segment(index_type index);
+ // LEGACY
const class segment* parent(index_type index) const;
+ // LEGACY
const class segment* segment(index_type index) const;
// access pointer to the soma
// returns nullptr if the cell has no soma
// LEGACY
- soma_segment* soma();
const soma_segment* soma() const;
// access pointer to a cable segment
// will throw an cable_cell_error exception if
// the cable index is not valid
// LEGACY
- cable_segment* cable(index_type index);
const cable_segment* cable(index_type index) const;
+ // LEGACY
const std::vector<segment_ptr>& segments() const;
- // the number of segments in the cell
- size_type num_segments() const;
-
// return a vector with the compartment count for each segment in the cell
// LEGACY
std::vector<size_type> compartment_counts() const;
@@ -128,27 +128,28 @@ public:
//
// Density channels.
- void paint(const std::string& target, mechanism_desc);
+ void paint(const std::string&, mechanism_desc);
+ void paint(const region&, mechanism_desc);
+
+ // Properties.
+ void paint(const std::string&, cable_cell_local_parameter_set);
+ void paint(const region&, cable_cell_local_parameter_set);
// Synapses.
- lid_range place(const std::string& target, const mechanism_desc&);
- lid_range place(const mlocation&, const mechanism_desc&); // LEGACY
- //lid_range place(const locset&, const mechanism_desc&);
+ lid_range place(const std::string&, const mechanism_desc&);
+ lid_range place(const locset&, const mechanism_desc&);
// Stimuli.
- lid_range place(const std::string& target, const i_clamp&);
- lid_range place(const mlocation&, const i_clamp&); // LEGACY
- //lid_range place(const locset&, const i_clamp&);
+ lid_range place(const std::string&, const i_clamp&);
+ lid_range place(const locset&, const i_clamp&);
// Gap junctions.
lid_range place(const std::string&, gap_junction_site);
- lid_range place(const mlocation& loc, gap_junction_site); // LEGACY
- //lid_range place(const locset&, gap_junction_site);
+ lid_range place(const locset&, gap_junction_site);
// spike detectors
lid_range place(const std::string&, const threshold_detector&);
- lid_range place(const mlocation&, const threshold_detector&); // LEGACY
- //lid_range place(const locset&, const threshold_detector&);
+ lid_range place(const locset&, const threshold_detector&);
//
// access to placed items
@@ -159,13 +160,6 @@ public:
const std::vector<detector_instance>& detectors() const;
const std::vector<stimulus_instance>& stimuli() const;
- // These setters are temporary, for "side-loading" in make_cable_cell.
- // In the regions, locset and morphology descriptions will be passed directly
- // to the cable_cell constructor.
- void set_regions(region_map r);
- void set_locsets(locset_map l);
- void set_morphology(em_morphology m);
-
const em_morphology* morphology() const;
// Checks that two cells have the same
@@ -194,12 +188,4 @@ private:
std::unique_ptr<cable_cell_impl, void (*)(cable_cell_impl*)> impl_;
};
-// Create a cable cell from a morphology specification.
-// If compartments_from_discretization is true, set number of compartments
-// in each segment to be the number of piecewise linear sections in the
-// corresponding section of the morphology.
-cable_cell make_cable_cell(const morphology& morph,
- const label_dict& labels={},
- bool compartments_from_discretization=false);
-
} // namespace arb
diff --git a/example/dryrun/branch_cell.hpp b/example/dryrun/branch_cell.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..bbb818188a3694f4c3df060c6744c5a7246fc697
--- /dev/null
+++ b/example/dryrun/branch_cell.hpp
@@ -0,0 +1,124 @@
+#pragma once
+
+#include <array>
+#include <random>
+
+#include <nlohmann/json.hpp>
+
+#include <arbor/common_types.hpp>
+#include <arbor/cable_cell.hpp>
+
+#include <sup/json_params.hpp>
+
+// Parameters used to generate the random cell morphologies.
+struct cell_parameters {
+ cell_parameters() = default;
+
+ // Maximum number of levels in the cell (not including the soma)
+ unsigned max_depth = 5;
+
+ // The following parameters are described as ranges.
+ // The first value is at the soma, and the last value is used on the last level.
+ // Values at levels in between are found by linear interpolation.
+ std::array<double,2> branch_probs = {1.0, 0.5}; // Probability of a branch occuring.
+ std::array<unsigned,2> compartments = {20, 2}; // Compartment count on a branch.
+ std::array<double,2> lengths = {200, 20}; // Length of branch in μm.
+
+ // The number of synapses per cell.
+ unsigned synapses = 1;
+};
+
+cell_parameters parse_cell_parameters(nlohmann::json& json) {
+ cell_parameters params;
+ sup::param_from_json(params.max_depth, "depth", json);
+ sup::param_from_json(params.branch_probs, "branch-probs", json);
+ sup::param_from_json(params.compartments, "compartments", json);
+ sup::param_from_json(params.lengths, "lengths", json);
+ sup::param_from_json(params.synapses, "synapses", json);
+
+ return params;
+}
+
+// Helper used to interpolate in branch_cell.
+template <typename T>
+double interp(const std::array<T,2>& r, unsigned i, unsigned n) {
+ double p = i * 1./(n-1);
+ double r0 = r[0];
+ double r1 = r[1];
+ return r[0] + p*(r1-r0);
+}
+
+arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params) {
+ arb::sample_tree tree;
+
+ // Add soma.
+ double soma_radius = 12.6157/2.0;
+ tree.append(arb::mnpos, {{0,0,0,soma_radius}, 1}); // For area of 500 μm².
+
+ std::vector<std::vector<unsigned>> levels;
+ levels.push_back({0});
+
+ // Standard mersenne_twister_engine seeded with gid.
+ std::mt19937 gen(gid);
+ std::uniform_real_distribution<double> dis(0, 1);
+
+ double dend_radius = 0.5; // Diameter of 1 μm for each cable.
+
+ double dist_from_soma = soma_radius;
+ for (unsigned i=0; i<params.max_depth; ++i) {
+ // Branch prob at this level.
+ double bp = interp(params.branch_probs, i, params.max_depth);
+ // Length at this level.
+ double l = interp(params.lengths, i, params.max_depth);
+ // Number of compartments at this level.
+ unsigned nc = std::round(interp(params.compartments, i, params.max_depth));
+
+ std::vector<unsigned> sec_ids;
+ for (unsigned sec: levels[i]) {
+ for (unsigned j=0; j<2; ++j) {
+ if (dis(gen)<bp) {
+ auto z = dist_from_soma;
+ auto p = tree.append(sec, {{0,0,z,dend_radius}, 3});
+ if (nc>1) {
+ auto dz = l/nc;
+ for (unsigned k=1; k<nc; ++k) {
+ p = tree.append(p, {{0,0,z+k*dz, dend_radius}, 3});
+ }
+ }
+ sec_ids.push_back(tree.append(p, {{0,0,z+l,dend_radius}, 3}));
+ }
+ }
+ }
+ if (sec_ids.empty()) {
+ break;
+ }
+ levels.push_back(sec_ids);
+
+ dist_from_soma += l;
+ }
+
+ arb::label_dict d;
+
+ using arb::reg::tagged;
+ d.set("soma", tagged(1));
+ d.set("dendrites", join(tagged(3), tagged(4)));
+
+ arb::cable_cell cell(arb::morphology(tree, true), d, true);
+
+ cell.paint("soma", "hh");
+ cell.paint("dendrites", "pas");
+ cell.default_parameters.axial_resistivity = 100; // [Ω·cm]
+
+ // Add spike threshold detector at the soma.
+ cell.place(arb::mlocation{0,0}, arb::threshold_detector{10});
+
+ // Add a synapse to the mid point of the first dendrite.
+ cell.place(arb::mlocation{1, 0.5}, "expsyn");
+
+ // Add additional synapses that will not be connected to anything.
+ for (unsigned i=1u; i<params.synapses; ++i) {
+ cell.place(arb::mlocation{1, 0.5}, "expsyn");
+ }
+
+ return cell;
+}
diff --git a/example/dryrun/dryrun.cpp b/example/dryrun/dryrun.cpp
index ff300d93951699d26a4557598d91ffccafde49cb..ea21ea4ff3a1417672ca1ca78d381355cfc23fc4 100644
--- a/example/dryrun/dryrun.cpp
+++ b/example/dryrun/dryrun.cpp
@@ -27,24 +27,13 @@
#include <sup/json_meter.hpp>
#include <sup/json_params.hpp>
+#include "branch_cell.hpp"
+
#ifdef ARB_MPI_ENABLED
#include <mpi.h>
#include <arborenv/with_mpi.hpp>
#endif
-// Parameters used to generate the random cell morphologies.
-struct cell_parameters {
- // Maximum number of levels in the cell (not including the soma)
- unsigned max_depth = 5;
-
- // The following parameters are described as ranges.
- // The first value is at the soma, and the last value is used on the last level.
- // Values at levels in between are found by linear interpolation.
- std::array<double,2> branch_probs = {1.0, 0.5}; // Probability of a branch occuring.
- std::array<unsigned,2> compartments = {20, 2}; // Compartment count on a branch.
- std::array<double,2> lengths = {200, 20}; // Length of branch in μm.
-};
-
struct run_params {
std::string name = "default";
bool dry_run = false;
@@ -175,7 +164,7 @@ struct cell_stats {
size_type ncomp_tmp = 0;
for (size_type i=b; i<e; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
- nsegs_tmp += c.num_segments();
+ nsegs_tmp += c.num_branches();
ncomp_tmp += c.num_compartments();
}
MPI_Allreduce(&nsegs_tmp, &nsegs, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
@@ -187,7 +176,7 @@ struct cell_stats {
ncells = r.num_cells();
for (size_type i = 0; i < ncells; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
- nsegs += c.num_segments();
+ nsegs += c.num_branches();
ncomp += c.num_compartments();
}
}
@@ -198,7 +187,7 @@ struct cell_stats {
for (size_type i = 0; i < params.num_cells_per_rank; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
- nsegs += c.num_segments();
+ nsegs += c.num_branches();
ncomp += c.num_compartments();
}
@@ -211,12 +200,11 @@ struct cell_stats {
return o << "cell stats: "
<< s.nranks << " ranks; "
<< s.ncells << " cells; "
- << s.nsegs << " segments; "
+ << s.nsegs << " branches; "
<< s.ncomp << " compartments.";
}
};
-
int main(int argc, char** argv) {
try {
#ifdef ARB_MPI_ENABLED
@@ -324,67 +312,6 @@ int main(int argc, char** argv) {
return 0;
}
-// Helper used to interpolate in branch_cell.
-template <typename T>
-double interp(const std::array<T,2>& r, unsigned i, unsigned n) {
- double p = i * 1./(n-1);
- double r0 = r[0];
- double r1 = r[1];
- return r[0] + p*(r1-r0);
-}
-
-arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params) {
- arb::cable_cell cell;
- cell.default_parameters.axial_resistivity = 100; // [Ω·cm]
-
- // Add soma.
- auto soma = cell.add_soma(12.6157/2.0); // For area of 500 μm².
- soma->add_mechanism("hh");
-
- std::vector<std::vector<unsigned>> levels;
- levels.push_back({0});
-
- // Standard mersenne_twister_engine seeded with gid.
- std::mt19937 gen(gid);
- std::uniform_real_distribution<double> dis(0, 1);
-
- double dend_radius = 0.5; // Diameter of 1 μm for each cable.
-
- unsigned nsec = 1;
- for (unsigned i=0; i<params.max_depth; ++i) {
- // Branch prob at this level.
- double bp = interp(params.branch_probs, i, params.max_depth);
- // Length at this level.
- double l = interp(params.lengths, i, params.max_depth);
- // Number of compartments at this level.
- unsigned nc = std::round(interp(params.compartments, i, params.max_depth));
-
- std::vector<unsigned> sec_ids;
- for (unsigned sec: levels[i]) {
- for (unsigned j=0; j<2; ++j) {
- if (dis(gen)<bp) {
- sec_ids.push_back(nsec++);
- auto dend = cell.add_cable(sec, arb::make_segment<arb::cable_segment>(arb::section_kind::dendrite, dend_radius, dend_radius, l));
- dend->set_compartments(nc);
- dend->add_mechanism("pas");
- }
- }
- }
- if (sec_ids.empty()) {
- break;
- }
- levels.push_back(sec_ids);
- }
-
- // Add spike threshold detector at the soma.
- cell.place(arb::mlocation{0,0}, arb::threshold_detector{10});
-
- // Add a synapse to the mid point of the first dendrite.
- cell.place(arb::mlocation{1, 0.5}, "expsyn");
-
- return cell;
-}
-
run_params read_options(int argc, char** argv) {
using sup::param_from_json;
@@ -414,10 +341,7 @@ run_params read_options(int argc, char** argv) {
param_from_json(params.num_ranks, "num-ranks", json);
param_from_json(params.duration, "duration", json);
param_from_json(params.min_delay, "min-delay", json);
- param_from_json(params.cell.max_depth, "depth", json);
- param_from_json(params.cell.branch_probs, "branch-probs", json);
- param_from_json(params.cell.compartments, "compartments", json);
- param_from_json(params.cell.lengths, "lengths", json);
+ params.cell = parse_cell_parameters(json);
if (!json.empty()) {
for (auto it=json.begin(); it!=json.end(); ++it) {
diff --git a/example/gap_junctions/gap_junctions.cpp b/example/gap_junctions/gap_junctions.cpp
index 4e3aa846f75176670c1060b23abab22c5812b617..32c789d9950c31ea81f491fdc11808a4110e94ad 100644
--- a/example/gap_junctions/gap_junctions.cpp
+++ b/example/gap_junctions/gap_junctions.cpp
@@ -161,7 +161,7 @@ struct cell_stats {
size_type ncomp_tmp = 0;
for (size_type i=b; i<e; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
- nsegs_tmp += c.num_segments();
+ nsegs_tmp += c.num_branches();
ncomp_tmp += c.num_compartments();
}
MPI_Allreduce(&nsegs_tmp, &nsegs, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
@@ -170,7 +170,7 @@ struct cell_stats {
ncells = r.num_cells();
for (size_type i=0; i<ncells; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
- nsegs += c.num_segments();
+ nsegs += c.num_branches();
ncomp += c.num_compartments();
}
#endif
@@ -179,7 +179,7 @@ struct cell_stats {
friend std::ostream& operator<<(std::ostream& o, const cell_stats& s) {
return o << "cell stats: "
<< s.ncells << " cells; "
- << s.nsegs << " segments; "
+ << s.nsegs << " branchess; "
<< s.ncomp << " compartments.";
}
};
@@ -334,45 +334,52 @@ void write_trace_json(const std::vector<arb::trace_data<double>>& trace, unsigne
}
arb::cable_cell gj_cell(cell_gid_type gid, unsigned ncell, double stim_duration) {
- arb::cable_cell cell;
+ // Create the sample tree that defines the morphology.
+ arb::sample_tree tree;
+ double soma_rad = 22.360679775/2.0; // convert diameter to radius in μm
+ tree.append({{0,0,0,soma_rad}, 1}); // soma is a single sample point
+ double dend_rad = 3./2;
+ tree.append(0, {{0,0,soma_rad, dend_rad}, 3}); // proximal point of the dendrite
+ tree.append(1, {{0,0,soma_rad+300, dend_rad}, 3}); // distal end of the dendrite
+
+ // Create a label dictionary that creates a single region that covers the whole cell.
+ arb::label_dict d;
+ d.set("all", arb::reg::all());
+
+ // Create the cell and set its electrical properties.
+ arb::cable_cell cell(tree, d);
cell.default_parameters.axial_resistivity = 100; // [Ω·cm]
cell.default_parameters.membrane_capacitance = 0.018; // [F/m²]
+ // Define the density channels and their parameters.
arb::mechanism_desc nax("nax");
+ nax["gbar"] = 0.04;
+ nax["sh"] = 10;
+
arb::mechanism_desc kdrmt("kdrmt");
+ kdrmt["gbar"] = 0.0001;
+
arb::mechanism_desc kamt("kamt");
+ kamt["gbar"] = 0.004;
+
arb::mechanism_desc pas("pas");
+ pas["g"] = 1.0/12000.0;
+ pas["e"] = -65;
- auto set_reg_params = [&]() {
- nax["gbar"] = 0.04;
- nax["sh"] = 10;
- kdrmt["gbar"] = 0.0001;
- kamt["gbar"] = 0.004;
- pas["g"] = 1.0/12000.0;
- pas["e"] = -65;
- };
-
- auto setup_seg = [&](auto seg) {
- seg->add_mechanism(nax);
- seg->add_mechanism(kdrmt);
- seg->add_mechanism(kamt);
- seg->add_mechanism(pas);
- };
-
- auto soma = cell.add_soma(22.360679775/2.0);
- set_reg_params();
- setup_seg(soma);
-
- auto dend = cell.add_cable(0, arb::make_segment<arb::cable_segment>(arb::section_kind::dendrite, 3.0/2.0, 3.0/2.0, 300)); //cable 1
- dend->set_compartments(1);
- set_reg_params();
- setup_seg(dend);
+ // Paint density channels on all parts of the cell
+ cell.paint("all", nax);
+ cell.paint("all", kdrmt);
+ cell.paint("all", kamt);
+ cell.paint("all", pas);
+ // Add a spike detector to the soma.
cell.place(arb::mlocation{0,0}, arb::threshold_detector{10});
+ // Add two gap junction sites.
cell.place(arb::mlocation{0, 1}, arb::gap_junction_site{});
cell.place(arb::mlocation{1, 1}, arb::gap_junction_site{});
+ // Attach a stimulus to the second cell.
if (!gid) {
arb::i_clamp stim(0, stim_duration, 0.4);
cell.place(arb::mlocation{0, 0.5}, stim);
diff --git a/example/generators/generators.cpp b/example/generators/generators.cpp
index 9d48f0c6bb1f999a2785eb564219c3720b44b13a..e5234ec96980dae656103e1697d7ed7389ba58d5 100644
--- a/example/generators/generators.cpp
+++ b/example/generators/generators.cpp
@@ -47,10 +47,15 @@ public:
// capacitance: 0.01 F/m² [default]
// synapses: 1 * expsyn
arb::util::unique_any get_cell_description(cell_gid_type gid) const override {
- arb::cable_cell c;
+ arb::sample_tree tree;
+ double r = 18.8/2.0; // convert 18.8 μm diameter to radius
+ tree.append({{0,0,0,r}, 1});
- c.add_soma(18.8/2.0); // convert 18.8 μm diameter to radius
- c.soma()->add_mechanism("pas");
+ arb::label_dict d;
+ d.set("soma", arb::reg::tagged(1));
+
+ arb::cable_cell c(tree, d);
+ c.paint("soma", "pas");
// Add one synapse at the soma.
// This synapse will be the target for all events, from both
diff --git a/example/ring/branch_cell.hpp b/example/ring/branch_cell.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..bbb818188a3694f4c3df060c6744c5a7246fc697
--- /dev/null
+++ b/example/ring/branch_cell.hpp
@@ -0,0 +1,124 @@
+#pragma once
+
+#include <array>
+#include <random>
+
+#include <nlohmann/json.hpp>
+
+#include <arbor/common_types.hpp>
+#include <arbor/cable_cell.hpp>
+
+#include <sup/json_params.hpp>
+
+// Parameters used to generate the random cell morphologies.
+struct cell_parameters {
+ cell_parameters() = default;
+
+ // Maximum number of levels in the cell (not including the soma)
+ unsigned max_depth = 5;
+
+ // The following parameters are described as ranges.
+ // The first value is at the soma, and the last value is used on the last level.
+ // Values at levels in between are found by linear interpolation.
+ std::array<double,2> branch_probs = {1.0, 0.5}; // Probability of a branch occuring.
+ std::array<unsigned,2> compartments = {20, 2}; // Compartment count on a branch.
+ std::array<double,2> lengths = {200, 20}; // Length of branch in μm.
+
+ // The number of synapses per cell.
+ unsigned synapses = 1;
+};
+
+cell_parameters parse_cell_parameters(nlohmann::json& json) {
+ cell_parameters params;
+ sup::param_from_json(params.max_depth, "depth", json);
+ sup::param_from_json(params.branch_probs, "branch-probs", json);
+ sup::param_from_json(params.compartments, "compartments", json);
+ sup::param_from_json(params.lengths, "lengths", json);
+ sup::param_from_json(params.synapses, "synapses", json);
+
+ return params;
+}
+
+// Helper used to interpolate in branch_cell.
+template <typename T>
+double interp(const std::array<T,2>& r, unsigned i, unsigned n) {
+ double p = i * 1./(n-1);
+ double r0 = r[0];
+ double r1 = r[1];
+ return r[0] + p*(r1-r0);
+}
+
+arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params) {
+ arb::sample_tree tree;
+
+ // Add soma.
+ double soma_radius = 12.6157/2.0;
+ tree.append(arb::mnpos, {{0,0,0,soma_radius}, 1}); // For area of 500 μm².
+
+ std::vector<std::vector<unsigned>> levels;
+ levels.push_back({0});
+
+ // Standard mersenne_twister_engine seeded with gid.
+ std::mt19937 gen(gid);
+ std::uniform_real_distribution<double> dis(0, 1);
+
+ double dend_radius = 0.5; // Diameter of 1 μm for each cable.
+
+ double dist_from_soma = soma_radius;
+ for (unsigned i=0; i<params.max_depth; ++i) {
+ // Branch prob at this level.
+ double bp = interp(params.branch_probs, i, params.max_depth);
+ // Length at this level.
+ double l = interp(params.lengths, i, params.max_depth);
+ // Number of compartments at this level.
+ unsigned nc = std::round(interp(params.compartments, i, params.max_depth));
+
+ std::vector<unsigned> sec_ids;
+ for (unsigned sec: levels[i]) {
+ for (unsigned j=0; j<2; ++j) {
+ if (dis(gen)<bp) {
+ auto z = dist_from_soma;
+ auto p = tree.append(sec, {{0,0,z,dend_radius}, 3});
+ if (nc>1) {
+ auto dz = l/nc;
+ for (unsigned k=1; k<nc; ++k) {
+ p = tree.append(p, {{0,0,z+k*dz, dend_radius}, 3});
+ }
+ }
+ sec_ids.push_back(tree.append(p, {{0,0,z+l,dend_radius}, 3}));
+ }
+ }
+ }
+ if (sec_ids.empty()) {
+ break;
+ }
+ levels.push_back(sec_ids);
+
+ dist_from_soma += l;
+ }
+
+ arb::label_dict d;
+
+ using arb::reg::tagged;
+ d.set("soma", tagged(1));
+ d.set("dendrites", join(tagged(3), tagged(4)));
+
+ arb::cable_cell cell(arb::morphology(tree, true), d, true);
+
+ cell.paint("soma", "hh");
+ cell.paint("dendrites", "pas");
+ cell.default_parameters.axial_resistivity = 100; // [Ω·cm]
+
+ // Add spike threshold detector at the soma.
+ cell.place(arb::mlocation{0,0}, arb::threshold_detector{10});
+
+ // Add a synapse to the mid point of the first dendrite.
+ cell.place(arb::mlocation{1, 0.5}, "expsyn");
+
+ // Add additional synapses that will not be connected to anything.
+ for (unsigned i=1u; i<params.synapses; ++i) {
+ cell.place(arb::mlocation{1, 0.5}, "expsyn");
+ }
+
+ return cell;
+}
diff --git a/example/ring/ring.cpp b/example/ring/ring.cpp
index ac7d344f304f569a6a7dbba52a08c200db0b748d..2220d6f4eb0ac84d8da60dcf2578b2f27b29bf3e 100644
--- a/example/ring/ring.cpp
+++ b/example/ring/ring.cpp
@@ -29,29 +29,13 @@
#include <sup/json_meter.hpp>
#include <sup/json_params.hpp>
+#include "branch_cell.hpp"
+
#ifdef ARB_MPI_ENABLED
#include <mpi.h>
#include <arborenv/with_mpi.hpp>
#endif
-// Parameters used to generate the random cell morphologies.
-struct cell_parameters {
- cell_parameters() = default;
-
- // Maximum number of levels in the cell (not including the soma)
- unsigned max_depth = 5;
-
- // The following parameters are described as ranges.
- // The first value is at the soma, and the last value is used on the last level.
- // Values at levels in between are found by linear interpolation.
- std::array<double,2> branch_probs = {1.0, 0.5}; // Probability of a branch occuring.
- std::array<unsigned,2> compartments = {20, 2}; // Compartment count on a branch.
- std::array<double,2> lengths = {200, 20}; // Length of branch in μm.
-
- // The number of synapses per cell.
- unsigned synapses = 1;
-};
-
struct ring_params {
ring_params() = default;
@@ -171,7 +155,7 @@ struct cell_stats {
size_type ncomp_tmp = 0;
for (size_type i=b; i<e; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
- nsegs_tmp += c.num_segments();
+ nsegs_tmp += c.num_branches();
ncomp_tmp += c.num_compartments();
}
MPI_Allreduce(&nsegs_tmp, &nsegs, 1, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD);
@@ -180,7 +164,7 @@ struct cell_stats {
ncells = r.num_cells();
for (size_type i=0; i<ncells; ++i) {
auto c = arb::util::any_cast<arb::cable_cell>(r.get_cell_description(i));
- nsegs += c.num_segments();
+ nsegs += c.num_branches();
ncomp += c.num_compartments();
}
#endif
@@ -189,7 +173,7 @@ struct cell_stats {
friend std::ostream& operator<<(std::ostream& o, const cell_stats& s) {
return o << "cell stats: "
<< s.ncells << " cells; "
- << s.nsegs << " segments; "
+ << s.nsegs << " branches; "
<< s.ncomp << " compartments.";
}
};
@@ -327,72 +311,6 @@ void write_trace_json(const arb::trace_data<double>& trace) {
file << std::setw(1) << json << "\n";
}
-// Helper used to interpolate in branch_cell.
-template <typename T>
-double interp(const std::array<T,2>& r, unsigned i, unsigned n) {
- double p = i * 1./(n-1);
- double r0 = r[0];
- double r1 = r[1];
- return r[0] + p*(r1-r0);
-}
-
-arb::cable_cell branch_cell(arb::cell_gid_type gid, const cell_parameters& params) {
- arb::cable_cell cell;
- cell.default_parameters.axial_resistivity = 100; // [Ω·cm]
-
- // Add soma.
- auto soma = cell.add_soma(12.6157/2.0); // For area of 500 μm².
- soma->add_mechanism("hh");
-
- std::vector<std::vector<unsigned>> levels;
- levels.push_back({0});
-
- // Standard mersenne_twister_engine seeded with gid.
- std::mt19937 gen(gid);
- std::uniform_real_distribution<double> dis(0, 1);
-
- double dend_radius = 0.5; // Diameter of 1 μm for each cable.
-
- unsigned nsec = 1;
- for (unsigned i=0; i<params.max_depth; ++i) {
- // Branch prob at this level.
- double bp = interp(params.branch_probs, i, params.max_depth);
- // Length at this level.
- double l = interp(params.lengths, i, params.max_depth);
- // Number of compartments at this level.
- unsigned nc = std::round(interp(params.compartments, i, params.max_depth));
-
- std::vector<unsigned> sec_ids;
- for (unsigned sec: levels[i]) {
- for (unsigned j=0; j<2; ++j) {
- if (dis(gen)<bp) {
- sec_ids.push_back(nsec++);
- auto dend = cell.add_cable(sec, arb::make_segment<arb::cable_segment>(arb::section_kind::dendrite, dend_radius, dend_radius, l));
- dend->set_compartments(nc);
- dend->add_mechanism("pas");
- }
- }
- }
- if (sec_ids.empty()) {
- break;
- }
- levels.push_back(sec_ids);
- }
-
- // Add spike threshold detector at the soma.
- cell.place(arb::mlocation{0,0}, arb::threshold_detector{10});
-
- // Add a synapse to the mid point of the first dendrite.
- cell.place(arb::mlocation{1, 0.5}, "expsyn");
-
- // Add additional synapses that will not be connected to anything.
- for (unsigned i=1u; i<params.synapses; ++i) {
- cell.place(arb::mlocation{1, 0.5}, "expsyn");
- }
-
- return cell;
-}
-
ring_params read_options(int argc, char** argv) {
using sup::param_from_json;
@@ -420,11 +338,7 @@ ring_params read_options(int argc, char** argv) {
param_from_json(params.num_cells, "num-cells", json);
param_from_json(params.duration, "duration", json);
param_from_json(params.min_delay, "min-delay", json);
- param_from_json(params.cell.max_depth, "depth", json);
- param_from_json(params.cell.branch_probs, "branch-probs", json);
- param_from_json(params.cell.compartments, "compartments", json);
- param_from_json(params.cell.lengths, "lengths", json);
- param_from_json(params.cell.synapses, "synapses", json);
+ params.cell = parse_cell_parameters(json);
if (!json.empty()) {
for (auto it=json.begin(); it!=json.end(); ++it) {
diff --git a/example/single/single.cpp b/example/single/single.cpp
index 501969f5de77ca9aafbe9fc7853cfd0945e248a5..3e7eb01758e13d17f68d09d7560be4b86f2f97e3 100644
--- a/example/single/single.cpp
+++ b/example/single/single.cpp
@@ -58,26 +58,28 @@ struct single_recipe: public arb::recipe {
using arb::reg::tagged;
dict.set("soma", tagged(1));
dict.set("dend", join(tagged(3), tagged(4), tagged(42)));
- arb::cable_cell c = make_cable_cell(morpho, dict, false);
+ arb::cable_cell c(morpho, dict, false);
// Add HH mechanism to soma, passive channels to dendrites.
c.paint("soma", "hh");
c.paint("dend", "pas");
// Discretize dendrites according to the NEURON d-lambda rule.
- for (std::size_t i=1; i<c.num_segments(); ++i) {
+ /* skip this during refactoring of the morphology interface
+ for (std::size_t i=1; i<c.num_branches(); ++i) {
arb::cable_segment* branch = c.cable(i);
double dx = c.segment_length_constant(100., i, gprop.default_parameters)*0.3;
unsigned n = std::ceil(branch->length()/dx);
branch->set_compartments(n);
}
+ */
// Add synapse to last branch.
- arb::cell_lid_type last_segment = c.num_segments()-1;
- arb::mlocation end_last_segment = { last_segment, 1. };
- c.place(end_last_segment, "exp2syn");
+ arb::cell_lid_type last_branch = c.num_branches()-1;
+ arb::mlocation end_last_branch = { last_branch, 1. };
+ c.place(end_last_branch, "exp2syn");
return c;
}
diff --git a/python/cells.cpp b/python/cells.cpp
index 21a93db2d8df0b9f6f7c7fc0c38d002e1df48d0f..8150dffe46313ff7f01c2d56ffb39bb81a4512d1 100644
--- a/python/cells.cpp
+++ b/python/cells.cpp
@@ -84,13 +84,11 @@ double interp(const std::array<T,2>& r, unsigned i, unsigned n) {
}
arb::cable_cell make_cable_cell(arb::cell_gid_type gid, const cell_parameters& params) {
- arb::cable_cell cell;
-
- cell.default_parameters.axial_resistivity = 100;
+ arb::sample_tree tree;
// Add soma.
- auto soma = cell.add_soma(12.6157/2.0); // For area of 500 μm².
- soma->add_mechanism("hh");
+ double soma_radius = 12.6157/2.0;
+ tree.append(arb::mnpos, {{0,0,0,soma_radius}, 1}); // For area of 500 μm².
std::vector<std::vector<unsigned>> levels;
levels.push_back({0});
@@ -101,7 +99,7 @@ arb::cable_cell make_cable_cell(arb::cell_gid_type gid, const cell_parameters& p
double dend_radius = 0.5; // Diameter of 1 μm for each cable.
- unsigned nsec = 1;
+ double dist_from_soma = soma_radius;
for (unsigned i=0; i<params.max_depth; ++i) {
// Branch prob at this level.
double bp = interp(params.branch_probs, i, params.max_depth);
@@ -114,10 +112,15 @@ arb::cable_cell make_cable_cell(arb::cell_gid_type gid, const cell_parameters& p
for (unsigned sec: levels[i]) {
for (unsigned j=0; j<2; ++j) {
if (dis(gen)<bp) {
- sec_ids.push_back(nsec++);
- auto dend = cell.add_cable(sec, arb::make_segment<arb::cable_segment>(arb::section_kind::dendrite, dend_radius, dend_radius, l));
- dend->set_compartments(nc);
- dend->add_mechanism("pas");
+ auto z = dist_from_soma;
+ auto p = tree.append(sec, {{0,0,z,dend_radius}, 3});
+ if (nc>1) {
+ auto dz = l/nc;
+ for (unsigned k=1; k<nc; ++k) {
+ p = tree.append(p, {{0,0,z+k*dz, dend_radius}, 3});
+ }
+ }
+ sec_ids.push_back(tree.append(p, {{0,0,z+l,dend_radius}, 3}));
}
}
}
@@ -125,15 +128,29 @@ arb::cable_cell make_cable_cell(arb::cell_gid_type gid, const cell_parameters& p
break;
}
levels.push_back(sec_ids);
+
+ dist_from_soma += l;
}
+ arb::label_dict d;
+
+ using arb::reg::tagged;
+ d.set("soma", tagged(1));
+ d.set("dendrites", join(tagged(3), tagged(4)));
+
+ arb::cable_cell cell(arb::morphology(tree, true), d, true);
+
+ cell.paint("soma", "hh");
+ cell.paint("dendrites", "pas");
+ cell.default_parameters.axial_resistivity = 100; // [Ω·cm]
+
// Add spike threshold detector at the soma.
cell.place(arb::mlocation{0,0}, arb::threshold_detector{10});
// Add a synapse to the mid point of the first dendrite.
cell.place(arb::mlocation{1, 0.5}, "expsyn");
- // Add additional synapses.
+ // Add additional synapses that will not be connected to anything.
for (unsigned i=1u; i<params.synapses; ++i) {
cell.place(arb::mlocation{1, 0.5}, "expsyn");
}
diff --git a/sup/include/sup/json_params.hpp b/sup/include/sup/json_params.hpp
index 7903b6cd6478a61af7a667b588519925721c6d48..fdcf0a20e2be4df7bffa7ec8ca06a5b457743c33 100644
--- a/sup/include/sup/json_params.hpp
+++ b/sup/include/sup/json_params.hpp
@@ -1,3 +1,5 @@
+#pragma once
+
#include <array>
#include <exception>
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index ef725caab660680f7b503a439d1f29f76f817178..7010a1a71af252a4183740d2f1dc1ee6dbf984e8 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -13,10 +13,6 @@ add_custom_target(tests)
# Builds: unit.
add_subdirectory(unit)
-# Model validation.
-# Builds: validate.
-add_subdirectory(validation)
-
# Test MPI wrappers and distribution operations.
# Builds: unit-local and unit-mpi (if MPI enabled).
add_subdirectory(unit-distributed)
diff --git a/test/common_cells.hpp b/test/common_cells.hpp
index dd40bc91d0add02c8f101651d685d1bc6ab8a316..02821167bc0c550413c8b763bf4c2af1fbc5468f 100644
--- a/test/common_cells.hpp
+++ b/test/common_cells.hpp
@@ -3,10 +3,88 @@
#include <arbor/cable_cell.hpp>
#include <arbor/segment.hpp>
#include <arbor/mechinfo.hpp>
+#include <arbor/morph/label_dict.hpp>
#include <arbor/recipe.hpp>
namespace arb {
+class soma_cell_builder {
+ double soma_rad;
+ sample_tree tree;
+ std::vector<msize_t> branch_distal_id;
+ std::unordered_map<std::string, int> tag_map;
+ int tag_count = 0;
+
+ // Get tag id of region.
+ // Add a new tag if region with that name has not already had a tag associated with it.
+ int get_tag(const std::string& name) {
+ auto it = tag_map.find(name);
+ // If the name is not in the map, make a unique tag.
+ // Tags start from 1.
+ if (it==tag_map.end()) {
+ tag_map[name] = ++tag_count;
+ return tag_count;
+ }
+ return it->second;
+ }
+
+public:
+ soma_cell_builder(double r): soma_rad(r) {
+ tree.append({{0,0,0,r}, get_tag("soma")});
+ branch_distal_id.push_back(0);
+ }
+
+ // Add a new branch that is attached to parent_branch.
+ // Returns the id of the new branch.
+ msize_t add_branch(msize_t parent_branch, double len, double r1, double r2, int ncomp,
+ const std::string& region)
+ {
+ // Get tag id of region (add a new tag if region does not already exist).
+ int tag = get_tag(region);
+
+ msize_t p = branch_distal_id[parent_branch];
+ double z = parent_branch? tree.samples()[p].loc.z: soma_rad;
+
+ p = tree.append(p, {{0,0,z,r1}, tag});
+ if (ncomp>1) {
+ double dz = len/ncomp;
+ double dr = (r2-r1)/ncomp;
+ for (auto i=1; i<ncomp; ++i) {
+ p = tree.append(p, {{0,0,z+i*dz, r1+i*dr}, tag});
+ }
+ }
+ p = tree.append(p, {{0,0,z+len,r2}, tag});
+ branch_distal_id.push_back(p);
+
+ return branch_distal_id.size()-1;
+ }
+
+ cable_cell make_cell() const {
+ // Test that a valid tree was generated, that is, every branch has
+ // either 0 children, or at least 2 children.
+ for (auto i: branch_distal_id) {
+ if (i==0) continue;
+ auto prop = tree.properties()[i];
+ if (!is_fork(prop) && !is_terminal(prop)) {
+ throw cable_cell_error(
+ "attempt to construct a cable_cell from a soma_cell_builder "
+ "where a branch has only one child branch.");
+ }
+ }
+
+ // Make label dictionary with one entry for each tag.
+ label_dict dict;
+ for (auto& tag: tag_map) {
+ dict.set(tag.first, reg::tagged(tag.second));
+ }
+
+ // Make cable_cell from sample tree and dictionary.
+ // The true flag is used to force the discretization to make compartments
+ // at sample points.
+ return cable_cell(tree, dict, true);
+ }
+};
+
/*
* Create cell with just a soma:
*
@@ -21,11 +99,10 @@ namespace arb {
*/
inline cable_cell make_cell_soma_only(bool with_stim = true) {
- cable_cell c;
-
- auto soma = c.add_soma(18.8/2.0);
- soma->add_mechanism("hh");
+ soma_cell_builder builder(18.8/2.0);
+ auto c = builder.make_cell();
+ c.paint("soma", "hh");
if (with_stim) {
c.place(mlocation{0,0.5}, i_clamp{10., 100., 0.1});
}
@@ -55,125 +132,16 @@ inline cable_cell make_cell_soma_only(bool with_stim = true) {
*/
inline cable_cell make_cell_ball_and_stick(bool with_stim = true) {
- cable_cell c;
-
- auto soma = c.add_soma(12.6157/2.0);
- soma->add_mechanism("hh");
-
- c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 1.0/2, 1.0/2, 200.0));
-
- for (auto& seg: c.segments()) {
- if (seg->is_dendrite()) {
- seg->add_mechanism("pas");
- seg->set_compartments(4);
- }
- }
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 4, "dend");
+ auto c = builder.make_cell();
+ c.paint("soma", "hh");
+ c.paint("dend", "pas");
if (with_stim) {
- c.place(mlocation{1,1}, i_clamp{5., 80., 0.3});
- }
- return c;
-}
-
-/*
- * Create cell with a soma and unbranched tapered dendrite:
- *
- * Common properties:
- * bulk resistivity: 100 Ω·cm [default]
- * capacitance: 0.01 F/m² [default]
- *
- * Soma:
- * mechanisms: HH (default params)
- * diameter: 12.6157 µm
- *
- * Dendrite:
- * mechanisms: passive (default params)
- * diameter proximal: 1 µm
- * diameter distal: 0.4 µm
- * length: 200 µm
- * bulk resistivity: 100 Ω·cm
- * compartments: 4
- *
- * Stimulus:
- * end of dendrite, t=[5 ms, 85 ms), 0.3 nA
- */
-
-inline cable_cell make_cell_ball_and_taper(bool with_stim = true) {
- cable_cell c;
-
- auto soma = c.add_soma(12.6157/2.0);
- soma->add_mechanism("hh");
-
- c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 1.0/2, 0.4/2, 200.0));
-
- for (auto& seg: c.segments()) {
- if (seg->is_dendrite()) {
- seg->add_mechanism("pas");
- seg->set_compartments(4);
- }
+ c.place(mlocation{1,1}, i_clamp{5, 80, 0.3});
}
- if (with_stim) {
- c.place(mlocation{1,1}, i_clamp{5., 80., 0.3});
- }
- return c;
-}
-
-/*
- * Create cell with a soma and unbranched dendrite with varying diameter:
- *
- * Common properties:
- * bulk resistivity: 100 Ω·cm [default]
- * capacitance: 0.01 F/m² [default]
- *
- * Soma:
- * mechanisms: HH
- * diameter: 12.6157 µm
- *
- * Dendrite:
- * mechanisms: passive (default params)
- * length: 100 µm
- * membrane resistance: 100 Ω·cm
- * compartments: 4
- *
- * Stimulus:
- * end of dendrite, t=[5 ms, 85 ms), 0.3 nA
- */
-
-inline cable_cell make_cell_ball_and_squiggle(bool with_stim = true) {
- cable_cell c;
-
- auto soma = c.add_soma(12.6157/2.0);
- soma->add_mechanism("hh");
-
- std::vector<cable_cell::value_type> radii;
- std::vector<cable_cell::point_type> points;
-
- double length = 100.0;
- int npoints = 200;
-
- for (int i=0; i<npoints; ++i) {
- double x = i*(1.0/(npoints-1));
- double r = std::exp(-x)*(std::sin(40*x)*0.05+0.1)+0.1;
-
- radii.push_back(r);
- points.push_back({x*length, 0., 0.});
- };
-
- auto dendrite =
- make_segment<cable_segment>(section_kind::dendrite, radii, points);
- c.add_cable(0, std::move(dendrite));
-
- for (auto& seg: c.segments()) {
- if (seg->is_dendrite()) {
- seg->add_mechanism("pas");
- seg->set_compartments(4);
- }
- }
-
- if (with_stim) {
- c.place(mlocation{1,1}, i_clamp{5., 80., 0.3});
- }
return c;
}
@@ -202,78 +170,20 @@ inline cable_cell make_cell_ball_and_squiggle(bool with_stim = true) {
*/
inline cable_cell make_cell_ball_and_3stick(bool with_stim = true) {
- cable_cell c;
-
- auto soma = c.add_soma(12.6157/2.0);
- soma->add_mechanism("hh");
-
- c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 100));
- c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 100));
- c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 100));
-
- for (auto& seg: c.segments()) {
- if (seg->is_dendrite()) {
- seg->add_mechanism("pas");
- seg->set_compartments(4);
- }
- }
-
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 100, 0.5, 0.5, 4, "dend");
+ builder.add_branch(1, 100, 0.5, 0.5, 4, "dend");
+ builder.add_branch(1, 100, 0.5, 0.5, 4, "dend");
+
+ auto c = builder.make_cell();
+ c.paint("soma", "hh");
+ c.paint("dend", "pas");
if (with_stim) {
c.place(mlocation{2,1}, i_clamp{5., 80., 0.45});
c.place(mlocation{3,1}, i_clamp{40., 10.,-0.2});
}
- return c;
-}
-
-/*
- * Create 'soma-less' cell with single cable, with physical
- * parameters from Rallpack 1 model.
- *
- * Common properties:
- * mechanisms: passive
- * membrane conductance: 0.000025 S/cm² ( = 1/(4Ω·m²) )
- * membrane reversal potential: -65 mV (default)
- * diameter: 1 µm
- * length: 1000 µm
- * bulk resistivity: 100 Ω·cm [default]
- * capacitance: 0.01 F/m² [default]
- * compartments: 4
- *
- * Stimulus:
- * end of dendrite, t=[0 ms, inf), 0.1 nA
- *
- * Note: zero-volume soma added with same mechanisms, as
- * work-around for some existing fvm modelling issues.
- */
-
-inline cable_cell make_cell_simple_cable(bool with_stim = true) {
- cable_cell c;
-
- c.default_parameters.axial_resistivity = 100;
- c.default_parameters.membrane_capacitance = 0.01;
-
- c.add_soma(0);
- c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 1000));
-
- double gbar = 0.000025;
- double I = 0.1;
-
- mechanism_desc pas("pas");
- pas["g"] = gbar;
-
- for (auto& seg: c.segments()) {
- seg->add_mechanism(pas);
-
- if (seg->is_dendrite()) {
- seg->set_compartments(4);
- }
- }
- if (with_stim) {
- // stimulus in the middle of our zero-volume 'soma'
- // corresponds to proximal end of cable.
- c.place(mlocation{0,0.5}, i_clamp{0., INFINITY, I});
- }
return c;
}
+
} // namespace arb
diff --git a/test/unit/test_cable_cell.cpp b/test/unit/test_cable_cell.cpp
index 49d16a9e74b4e82d79f4a7f4b7350d0e8423db96..62c54706886dee11cf1b5ba829e9c9ede03fb1b1 100644
--- a/test/unit/test_cable_cell.cpp
+++ b/test/unit/test_cable_cell.cpp
@@ -4,6 +4,8 @@
#include <arbor/math.hpp>
#include <arbor/morph/locset.hpp>
+#include "io/sepval.hpp"
+
#include "tree.hpp"
using namespace arb;
@@ -11,226 +13,69 @@ using ::arb::math::pi;
TEST(cable_cell, soma) {
// test that insertion of a soma works
- // define with no centre point
- {
- cable_cell c;
- auto soma_radius = 2.1;
+ // define with centre point @ (0,0,1)
+ double soma_radius = 3.2;
- EXPECT_EQ(c.has_soma(), false);
- c.add_soma(soma_radius);
- EXPECT_EQ(c.has_soma(), true);
+ arb::sample_tree samples;
+ samples.append({0,0,0,soma_radius,1});
+ auto c = cable_cell(arb::morphology(samples));
- auto s = c.soma();
- EXPECT_EQ(s->radius(), soma_radius);
- EXPECT_EQ(s->center().is_set(), false);
- }
+ EXPECT_EQ(c.has_soma(), true);
- // test that insertion of a soma works
- // define with centre point @ (0,0,1)
- {
- cable_cell c;
- auto soma_radius = 3.2;
-
- EXPECT_EQ(c.has_soma(), false);
- c.add_soma(soma_radius, {0,0,1});
- EXPECT_EQ(c.has_soma(), true);
-
- auto s = c.soma();
- EXPECT_EQ(s->radius(), soma_radius);
- EXPECT_EQ(s->center().is_set(), true);
- }
-}
-
-TEST(cable_cell, add_segment) {
- // add a pre-defined segment
- {
- cable_cell c;
-
- auto soma_radius = 2.1;
- auto cable_radius = 0.1;
- auto cable_length = 8.3;
-
- // add a soma because we need something to attach the first dendrite to
- c.add_soma(soma_radius, {0,0,1});
-
- auto seg =
- make_segment<cable_segment>(
- section_kind::dendrite,
- cable_radius, cable_radius, cable_length
- );
- c.add_cable(0, std::move(seg));
-
- EXPECT_EQ(c.num_segments(), 2u);
- }
-
- // add segment on the fly
- {
- cable_cell c;
-
- auto soma_radius = 2.1;
- auto cable_radius = 0.1;
- auto cable_length = 8.3;
-
- // add a soma because we need something to attach the first dendrite to
- c.add_soma(soma_radius, {0,0,1});
-
- c.add_cable(
- 0,
- make_segment<cable_segment>(section_kind::dendrite, cable_radius, cable_radius, cable_length)
- );
-
- EXPECT_EQ(c.num_segments(), 2u);
- }
- {
- cable_cell c;
-
- auto soma_radius = 2.1;
- auto cable_radius = 0.1;
- auto cable_length = 8.3;
-
- // add a soma because we need something to attach the first dendrite to
- c.add_soma(soma_radius, {0,0,1});
-
- c.add_cable(
- 0,
- make_segment<cable_segment>(section_kind::dendrite,
- std::vector<double>{cable_radius, cable_radius, cable_radius, cable_radius},
- std::vector<double>{cable_length, cable_length, cable_length})
- );
-
- EXPECT_EQ(c.num_segments(), 2u);
- }
+ auto s = c.soma();
+ EXPECT_EQ(s->radius(), soma_radius);
}
TEST(cable_cell, multiple_cables) {
- // generate a cylindrical cable segment of length 1/pi and radius 1
+ double soma_radius = std::pow(3./(4.*pi<double>), 1./3.);
+
+ // Generate a cylindrical cable segment of length 1/pi and radius 1
// volume = 1
// area = 2
- auto seg = [](section_kind k) {
- return make_segment<cable_segment>( k, 1.0, 1.0, 1./pi<double> );
+ // Returns the distal point of the added cable.
+ auto append_branch = [soma_radius](sample_tree& stree, int proximal) {
+ constexpr int tag = 2;
+ if (!proximal) {
+ double z = soma_radius;
+ proximal = stree.append(0, {0,0,z, 1/pi<double>, tag});
+ }
+ return stree.append(proximal, msample{0, 0, stree.samples()[proximal].loc.z+1, 1/pi<double>, tag});
};
- // add a pre-defined segment
- {
- cable_cell c;
-
- auto soma_radius = std::pow(3./(4.*pi<double>), 1./3.);
-
- // cell strucure as follows
- // left : segment numbering
- // right : segment type (soma, axon, dendrite)
- //
- // 0 s
- // / \ / \.
- // 1 2 d a
- // / \ / \.
- // 3 4 d d
-
- // add a soma
- c.add_soma(soma_radius, {0,0,1});
-
- // hook the dendrite and axons
- c.add_cable(0, seg(section_kind::dendrite));
- c.add_cable(0, seg(section_kind::axon));
- c.add_cable(1, seg(section_kind::dendrite));
- c.add_cable(1, seg(section_kind::dendrite));
-
- EXPECT_EQ(c.num_segments(), 5u);
-
- // construct the graph
- tree con(c.parents());
-
- auto no_parent = tree::no_parent;
- EXPECT_EQ(con.num_segments(), 5u);
- EXPECT_EQ(con.parent(0), no_parent);
- EXPECT_EQ(con.parent(1), 0u);
- EXPECT_EQ(con.parent(2), 0u);
- EXPECT_EQ(con.parent(3), 1u);
- EXPECT_EQ(con.parent(4), 1u);
- EXPECT_EQ(con.num_children(0), 2u);
- EXPECT_EQ(con.num_children(1), 2u);
- EXPECT_EQ(con.num_children(2), 0u);
- EXPECT_EQ(con.num_children(3), 0u);
- EXPECT_EQ(con.num_children(4), 0u);
-
- }
-}
-
-TEST(cable_cell, unbranched_chain) {
- cable_cell c;
-
- auto soma_radius = std::pow(3./(4.*pi<double>), 1./3.);
-
- // Cell strucure that looks like a centipede: i.e. each segment has only one child
+ // cell strucure with branch numbers
//
- // | |
- // 0|1-2-3-4|5-6-7-8
- // | |
+ // 0
+ // / \.
+ // 1 2
+ // / \.
+ // 3 4
- // add a soma
- c.add_soma(soma_radius, {0,0,1});
+ arb::sample_tree samples;
+ samples.append({0,0,-soma_radius,soma_radius,1});
// hook the dendrite and axons
- c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 1.0, 1.0, 1./pi<double>));
- c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 1.0, 1.0, 1./pi<double>));
+ append_branch(samples, 0);
+ append_branch(samples, 0);
+ append_branch(samples, 2);
+ append_branch(samples, 2);
+
+ auto c = cable_cell(arb::morphology(samples, true));
- EXPECT_EQ(c.num_segments(), 3u);
+ EXPECT_EQ(c.num_branches(), 5u);
// construct the graph
tree con(c.parents());
auto no_parent = tree::no_parent;
- EXPECT_EQ(con.num_segments(), 3u);
+ EXPECT_EQ(con.num_segments(), 5u);
EXPECT_EQ(con.parent(0), no_parent);
EXPECT_EQ(con.parent(1), 0u);
- EXPECT_EQ(con.parent(2), 1u);
- EXPECT_EQ(con.num_children(0), 1u);
- EXPECT_EQ(con.num_children(1), 1u);
+ EXPECT_EQ(con.parent(2), 0u);
+ EXPECT_EQ(con.parent(3), 1u);
+ EXPECT_EQ(con.parent(4), 1u);
+ EXPECT_EQ(con.num_children(0), 2u);
+ EXPECT_EQ(con.num_children(1), 2u);
EXPECT_EQ(con.num_children(2), 0u);
-}
-
-TEST(cable_cell, clone) {
- // make simple cell with multiple segments
-
- cable_cell c;
- c.add_soma(2.1);
- c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.3, 0.2, 10));
- c.segment(1)->set_compartments(3);
- c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.2, 0.15, 20));
- c.segment(2)->set_compartments(5);
-
- c.place(mlocation{1, 0.3}, "expsyn");
- c.place(mlocation{0, 0.5}, threshold_detector{10.0});
-
- // make clone
-
- cable_cell d(c);
-
- // check equality
-
- ASSERT_EQ(c.num_segments(), d.num_segments());
- EXPECT_EQ(c.soma()->radius(), d.soma()->radius());
- EXPECT_EQ(c.segment(1)->as_cable()->length(), d.segment(1)->as_cable()->length());
- {
- const auto& csyns = c.synapses();
- const auto& dsyns = d.synapses();
-
- ASSERT_EQ(csyns.size(), dsyns.size());
- for (unsigned i=0; i<csyns.size(); ++i) {
- ASSERT_EQ(csyns[i].location, dsyns[i].location);
- }
- }
-
- ASSERT_EQ(1u, c.detectors().size());
- ASSERT_EQ(1u, d.detectors().size());
- EXPECT_EQ(c.detectors()[0].threshold, d.detectors()[0].threshold);
-
- // check clone is independent
-
- c.add_cable(2, make_segment<cable_segment>(section_kind::dendrite, 0.15, 0.1, 20));
- EXPECT_NE(c.num_segments(), d.num_segments());
-
- c.segment(1)->set_compartments(7);
- EXPECT_NE(c.segment(1)->num_compartments(), d.segment(1)->num_compartments());
- EXPECT_EQ(c.segment(2)->num_compartments(), d.segment(2)->num_compartments());
+ EXPECT_EQ(con.num_children(3), 0u);
+ EXPECT_EQ(con.num_children(4), 0u);
}
diff --git a/test/unit/test_cv_policy.cpp b/test/unit/test_cv_policy.cpp
index a7a099682dcf54b38f57f660a5d87c343c590b43..273a413ef8423372ff8a63a9987bc616de4a67db 100644
--- a/test/unit/test_cv_policy.cpp
+++ b/test/unit/test_cv_policy.cpp
@@ -69,7 +69,7 @@ TEST(cv_policy, explicit_policy) {
cv_policy pol = cv_policy_explicit(lset);
for (auto& m: {m_sph_b6, m_reg_b6, m_mlt_b6}) {
- cable_cell cell = make_cable_cell(m);
+ cable_cell cell(m);
locset result = pol.cv_boundary_points(cell);
EXPECT_EQ(thingify(lset, *cell.morphology()), thingify(result, *cell.morphology()));
@@ -89,7 +89,7 @@ TEST(cv_policy, empty_morphology) {
cv_policy_max_extent(0.234, single_root_cv|interior_forks)
};
- cable_cell cell = make_cable_cell(m_empty);
+ cable_cell cell(m_empty);
auto empty_loclist = thingify(ls::nil(), *cell.morphology());
for (auto& pol: policies) {
@@ -112,7 +112,7 @@ TEST(cv_policy, single_root_cv) {
};
for (auto& morph: {m_sph_b1, m_reg_b1, m_sph_b6, m_reg_b6, m_mlt_b6}) {
- cable_cell cell = make_cable_cell(morph);
+ cable_cell cell(morph);
for (auto& polpair: policy_pairs) {
mlocation_list p1 = thingify(polpair.first.cv_boundary_points(cell), *cell.morphology());
@@ -133,7 +133,7 @@ TEST(cv_policy, fixed_per_branch) {
// root branch only
for (auto& morph: {m_sph_b1, m_reg_b1}) {
- cable_cell cell = make_cable_cell(morph);
+ cable_cell cell(morph);
{
// boundary fork points
cv_policy pol = cv_policy_fixed_per_branch(4);
@@ -153,7 +153,7 @@ TEST(cv_policy, fixed_per_branch) {
// spherical root, six branches and multiple top level branches cases:
// expected points are the same.
for (auto& morph: {m_sph_b6, m_mlt_b6}) {
- cable_cell cell = make_cable_cell(morph);
+ cable_cell cell(morph);
{
// boundary fork points
@@ -184,7 +184,7 @@ TEST(cv_policy, max_extent) {
// root branch only
for (auto& morph: {m_sph_b1, m_reg_b1}) {
- cable_cell cell = make_cable_cell(morph);
+ cable_cell cell(morph);
ASSERT_EQ(1.0, cell.morphology()->branch_length(0));
{
@@ -205,7 +205,7 @@ TEST(cv_policy, max_extent) {
// cell with varying branch lengths; extent not exact fraction.
{
- cable_cell cell = make_cable_cell(m_mlt_b6);
+ cable_cell cell(m_mlt_b6);
ASSERT_EQ(1.0, cell.morphology()->branch_length(0));
ASSERT_EQ(1.0, cell.morphology()->branch_length(1));
ASSERT_EQ(2.0, cell.morphology()->branch_length(2));
diff --git a/test/unit/test_domain_decomposition.cpp b/test/unit/test_domain_decomposition.cpp
index 30524dec9ff0aac23ea52dc61f86c753227b50ae..d941780eff13bc2b32584d46aa884e2ccd0bb873 100644
--- a/test/unit/test_domain_decomposition.cpp
+++ b/test/unit/test_domain_decomposition.cpp
@@ -8,6 +8,7 @@
#include "util/span.hpp"
+#include "../common_cells.hpp"
#include "../simple_recipes.hpp"
using namespace arb;
@@ -64,8 +65,7 @@ namespace {
}
arb::util::unique_any get_cell_description(cell_gid_type) const override {
- cable_cell c;
- c.add_soma(20);
+ auto c = arb::make_cell_soma_only(false);
c.place(mlocation{0,1}, gap_junction_site{});
return {std::move(c)};
}
diff --git a/test/unit/test_event_delivery.cpp b/test/unit/test_event_delivery.cpp
index a4d4b632c7edca338f6719c31a6c2181aa3c9713..148289e3c4c1d9cd7bda9d10f8f975bc290a8c04 100644
--- a/test/unit/test_event_delivery.cpp
+++ b/test/unit/test_event_delivery.cpp
@@ -27,11 +27,17 @@ struct test_recipe: public n_cable_cell_recipe {
explicit test_recipe(int n): n_cable_cell_recipe(n, test_cell()) {}
static cable_cell test_cell() {
- cable_cell c;
- c.add_soma(10.)->add_mechanism("pas");
+ sample_tree st;
+ st.append({0,0,0,10,1});
+
+ label_dict d;
+ d.set("soma", arb::reg::tagged(1));
+
+ cable_cell c(st, d, true);
c.place(mlocation{0, 0.5}, "expsyn");
c.place(mlocation{0, 0.5}, threshold_detector{-64});
c.place(mlocation{0, 0.5}, gap_junction_site{});
+
return c;
}
diff --git a/test/unit/test_fvm_layout.cpp b/test/unit/test_fvm_layout.cpp
index eb16024e46ce3575d087bdbac2b8dffffff6e56d..2808d3ea8fb379791b938989fb8e411d9a67bd8f 100644
--- a/test/unit/test_fvm_layout.cpp
+++ b/test/unit/test_fvm_layout.cpp
@@ -11,6 +11,7 @@
#include "util/rangeutil.hpp"
#include "util/span.hpp"
#include "io/sepval.hpp"
+#include "morph/em_morphology.hpp"
#include "common.hpp"
#include "unit_test_catalogue.hpp"
@@ -92,40 +93,40 @@ namespace {
//
// All dendrite segments with 4 compartments.
- cable_cell c2;
- segment* s;
+ soma_cell_builder builder(7.);
+ auto b1 = builder.add_branch(0, 200, 0.5, 0.5, 4, "dend");
+ auto b2 = builder.add_branch(1, 300, 0.4, 0.4, 4, "dend");
+ auto b3 = builder.add_branch(1, 180, 0.35, 0.35, 4, "dend");
+ auto cell = builder.make_cell();
- c2.default_parameters.axial_resistivity = 90;
+ cell.paint("soma", "hh");
+ cell.paint("dend", "pas");
- s = c2.add_soma(14./2);
- s->add_mechanism("hh");
+ 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;
- s = c2.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 1.0/2, 1.0/2, 200));
- s->parameters.membrane_capacitance = 0.017;
+ using ::arb::reg::branch;
+ cell.paint(branch(b1), p1);
+ cell.paint(branch(b2), p2);
+ cell.paint(branch(b3), p3);
- s = c2.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.8/2, 0.8/2, 300));
- s->parameters.membrane_capacitance = 0.013;
+ cell.place(mlocation{2,1}, i_clamp{5., 80., 0.45});
+ cell.place(mlocation{3,1}, i_clamp{40., 10.,-0.2});
- s = c2.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.7/2, 0.7/2, 180));
- s->parameters.membrane_capacitance = 0.018;
+ cell.default_parameters.axial_resistivity = 90;
- c2.place(mlocation{2,1}, i_clamp{5., 80., 0.45});
- c2.place(mlocation{3,1}, i_clamp{40., 10.,-0.2});
-
- for (auto& seg: c2.segments()) {
- if (seg->is_dendrite()) {
- seg->add_mechanism("pas");
- seg->set_compartments(4);
- }
- }
- cells.push_back(std::move(c2));
+ cells.push_back(std::move(cell));
return cells;
}
void check_two_cell_system(std::vector<cable_cell>& cells) {
- ASSERT_EQ(2u, cells[0].num_segments());
+ ASSERT_EQ(2u, cells[0].num_branches());
ASSERT_EQ(cells[0].segment(1)->num_compartments(), 4u);
- ASSERT_EQ(cells[1].num_segments(), 4u);
+ ASSERT_EQ(cells[1].num_branches(), 4u);
ASSERT_EQ(cells[1].segment(1)->num_compartments(), 4u);
ASSERT_EQ(cells[1].segment(2)->num_compartments(), 4u);
ASSERT_EQ(cells[1].segment(3)->num_compartments(), 4u);
@@ -186,7 +187,6 @@ TEST(fvm_layout, topology) {
EXPECT_EQ(ivec({0,0,1,2,3,4,6,6,7,8,9,10,11,12,13,14,11,16,17,18}), D.parent_cv);
-
EXPECT_FALSE(D.segments[0].has_parent());
EXPECT_EQ(1, D.segments[1].parent_cv);
@@ -246,7 +246,7 @@ TEST(fvm_layout, diam_and_area) {
std::vector<double> A;
for (auto ci: make_span(D.ncell)) {
- for (auto si: make_span(cells[ci].num_segments())) {
+ for (auto si: make_span(cells[ci].num_branches())) {
A.push_back(area(cells[ci].segment(si)));
}
}
@@ -285,7 +285,6 @@ TEST(fvm_layout, diam_and_area) {
double c = A[3]/(2*n)*cm1+A[4]/(2*n)*cm2+A[5]/(2*n)*cm3;
EXPECT_FLOAT_EQ(c, D.cv_capacitance[11]);
- //double cm0 = cells[1].soma()->cm;
double cm0 = neuron_parameter_defaults.membrane_capacitance.value();
c = A[2]*cm0;
EXPECT_FLOAT_EQ(c, D.cv_capacitance[6]);
@@ -295,7 +294,7 @@ TEST(fvm_layout, diam_and_area) {
// area, and h is the compartment length (given the
// regular discretization).
- cable_segment* cable = cells[1].segment(2)->as_cable();
+ auto cable = cells[1].segment(2)->as_cable();
double a = volume(cable)/cable->length();
EXPECT_FLOAT_EQ(math::pi<double>*0.8*0.8/4, a);
@@ -662,7 +661,6 @@ TEST(fvm_layout, synapse_targets) {
}
}
-
namespace {
double wm_impl(double wa, double xa) {
return wa? xa/wa: 0;
@@ -706,15 +704,12 @@ TEST(fvm_layout, density_norm_area) {
//
// Use divided compartment view on segments to compute area contributions.
- std::vector<cable_cell> cells(1);
- cable_cell& c = cells[0];
- auto soma = c.add_soma(12.6157/2.0);
-
- c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 100));
- c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.1, 200));
- c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.4, 0.4, 150));
+ soma_cell_builder builder(12.6157/2.0);
- auto& segs = c.segments();
+ // p len r1 r2 ncomp tag
+ builder.add_branch(0, 100, 0.5, 0.5, 3, "reg1");
+ builder.add_branch(1, 200, 0.5, 0.1, 3, "reg2");
+ builder.add_branch(1, 150, 0.4, 0.4, 3, "reg3");
double dflt_gkbar = .036;
double dflt_gl = 0.0003;
@@ -724,26 +719,25 @@ TEST(fvm_layout, density_norm_area) {
double seg3_gkbar = .0004;
double seg3_gl = .0004;
- for (int i = 0; i<4; ++i) {
- segment& seg = *segs[i];
- seg.set_compartments(3);
-
- mechanism_desc hh("hh");
- switch (i) {
- case 1:
- hh["gl"] = seg1_gl;
- break;
- case 2:
- hh["gkbar"] = seg2_gkbar;
- break;
- case 3:
- hh["gkbar"] = seg3_gkbar;
- hh["gl"] = seg3_gl;
- break;
- default: ;
- }
- seg.add_mechanism(hh);
- }
+ auto hh_0 = mechanism_desc("hh");
+
+ auto hh_1 = mechanism_desc("hh");
+ hh_1["gl"] = seg1_gl;
+
+ auto hh_2 = mechanism_desc("hh");
+ hh_2["gkbar"] = seg2_gkbar;
+
+ auto hh_3 = mechanism_desc("hh");
+ hh_3["gkbar"] = seg3_gkbar;
+ hh_3["gl"] = seg3_gl;
+
+ auto cell = builder.make_cell();
+ cell.paint("soma", std::move(hh_0));
+ cell.paint("reg1", std::move(hh_1));
+ cell.paint("reg2", std::move(hh_2));
+ cell.paint("reg3", std::move(hh_3));
+
+ std::vector<cable_cell> cells{std::move(cell)};
int ncv = 11; //ncomp + 1
std::vector<double> expected_gkbar(ncv, dflt_gkbar);
@@ -752,7 +746,8 @@ TEST(fvm_layout, density_norm_area) {
auto div_by_ends = [](const cable_segment* cable) {
return div_compartment_by_ends(cable->num_compartments(), cable->radii(), cable->lengths());
};
- double soma_area = area(soma);
+ auto& segs = cells[0].segments();
+ double soma_area = area(segs[0].get());
auto seg1_divs = div_by_ends(segs[1]->as_cable());
auto seg2_divs = div_by_ends(segs[2]->as_cable());
auto seg3_divs = div_by_ends(segs[3]->as_cable());
@@ -819,11 +814,9 @@ TEST(fvm_layout, density_norm_area) {
}
TEST(fvm_layout, valence_verify) {
- std::vector<cable_cell> cells(1);
- cable_cell& c = cells[0];
- auto soma = c.add_soma(6);
-
- soma->add_mechanism("test_cl_valence");
+ auto cell = soma_cell_builder(6).make_cell();
+ cell.paint("soma", "test_cl_valence");
+ std::vector<cable_cell> cells{std::move(cell)};
cable_cell_global_properties gprop;
gprop.default_parameters = neuron_parameter_defaults;
@@ -863,29 +856,27 @@ TEST(fvm_layout, ion_weights) {
//
// Geometry:
// soma 0: radius 5 µm
- // dend 1: 100 µm long, 1 µm diameter cynlinder
- // dend 2: 200 µm long, 1 µm diameter cynlinder
- // dend 3: 100 µm long, 1 µm diameter cynlinder
+ // dend 1: 100 µm long, 1 µm diameter cylinder, tag 2
+ // dend 2: 200 µm long, 1 µm diameter cylinder, tag 3
+ // dend 3: 100 µm long, 1 µm diameter cylinder, tag 4
//
// The radius of the soma is chosen such that the surface area of soma is
// the same as a 100µm dendrite, which makes it easier to describe the
// expected weights.
- auto construct_cell = [](cable_cell& c) {
- c.add_soma(5);
-
- c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 100));
- c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 200));
- c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 100));
-
- for (auto& s: c.segments()) s->set_compartments(1);
+ auto construct_cell = []() {
+ soma_cell_builder builder(5);
+ builder.add_branch(0, 100, 0.5, 0.5, 1, "dend");
+ builder.add_branch(1, 200, 0.5, 0.5, 1, "dend");
+ builder.add_branch(1, 100, 0.5, 0.5, 1, "dend");
+ return builder.make_cell();
};
using uvec = std::vector<fvm_size_type>;
using ivec = std::vector<fvm_index_type>;
using fvec = std::vector<fvm_value_type>;
- uvec mech_segs[] = {
+ uvec mech_branches[] = {
{0}, {0,2}, {2, 3}, {0, 1, 2, 3}, {3}
};
@@ -909,16 +900,16 @@ TEST(fvm_layout, ion_weights) {
}
}
- for (auto run: count_along(mech_segs)) {
+ for (auto run: count_along(mech_branches)) {
SCOPED_TRACE("run "+std::to_string(run));
- std::vector<cable_cell> cells(1);
- cable_cell& c = cells[0];
- construct_cell(c);
+ auto c = construct_cell();
- for (auto i: mech_segs[run]) {
- c.segments()[i]->add_mechanism("test_ca");
+ for (auto i: mech_branches[run]) {
+ c.paint(reg::branch(i), "test_ca");
}
+ std::vector<cable_cell> cells{std::move(c)};
+
fvm_discretization D = fvm_discretize(cells, gprop.default_parameters);
fvm_mechanism_data M = fvm_build_mechanism_data(gprop, cells, D);
@@ -943,24 +934,18 @@ TEST(fvm_layout, revpot) {
// * Reversal potential mechanisms are only extended where there exists another
// mechanism that reads them.
- auto construct_cell = [](cable_cell& c) {
- c.add_soma(5);
-
- c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 100));
- c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 200));
- c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 100));
-
- for (auto& s: c.segments()) s->set_compartments(1);
-
- // Read ea everywhere, ec only on soma.
- for (auto& s: c.segments()) s->add_mechanism("read_eX/a");
- c.soma()->add_mechanism("read_eX/c");
- };
-
mechanism_catalogue testcat = make_unit_test_catalogue();
- std::vector<cable_cell> cells(2);
- for (auto& c: cells) construct_cell(c);
+ soma_cell_builder builder(5);
+ builder.add_branch(0, 100, 0.5, 0.5, 1, "dend");
+ builder.add_branch(1, 200, 0.5, 0.5, 1, "dend");
+ builder.add_branch(1, 100, 0.5, 0.5, 1, "dend");
+ auto cell = builder.make_cell();
+ cell.paint("soma", "read_eX/c");
+ cell.paint("soma", "read_eX/a");
+ cell.paint("dend", "read_eX/a");
+
+ std::vector<cable_cell> cells{cell, cell};
cable_cell_global_properties gprop;
gprop.default_parameters = neuron_parameter_defaults;
diff --git a/test/unit/test_fvm_lowered.cpp b/test/unit/test_fvm_lowered.cpp
index a82132f497579e916fee647bd3a7fab7c84b6798..191249958991bbd5eb60680beeca7d8086e1fa20 100644
--- a/test/unit/test_fvm_lowered.cpp
+++ b/test/unit/test_fvm_lowered.cpp
@@ -92,8 +92,7 @@ public:
}
arb::util::unique_any get_cell_description(cell_gid_type gid) const override {
- cable_cell c;
- c.add_soma(20);
+ cable_cell c = soma_cell_builder(20).make_cell();
c.place(mlocation{0, 1}, gap_junction_site{});
return {std::move(c)};
}
@@ -138,9 +137,7 @@ public:
}
arb::util::unique_any get_cell_description(cell_gid_type) const override {
- cable_cell c;
- c.add_soma(20);
- return {std::move(c)};
+ return soma_cell_builder(20).make_cell();
}
cell_kind get_cell_kind(cell_gid_type gid) const override {
@@ -160,8 +157,7 @@ public:
}
arb::util::unique_any get_cell_description(cell_gid_type) const override {
- cable_cell c;
- c.add_soma(20);
+ cable_cell c = soma_cell_builder(20).make_cell();
c.place(mlocation{0,1}, gap_junction_site{});
return {std::move(c)};
}
@@ -213,10 +209,9 @@ TEST(fvm_lowered, matrix_init)
auto ispos = [](auto v) { return v>0; };
auto isneg = [](auto v) { return v<0; };
- cable_cell cell = make_cell_ball_and_stick();
-
- ASSERT_EQ(2u, cell.num_segments());
- cell.segment(1)->set_compartments(10);
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 200, 1.0/2, 1.0/2, 10, "dend"); // 10 compartments
+ cable_cell cell = builder.make_cell();
std::vector<target_handle> targets;
std::vector<fvm_index_type> cell_to_intdom;
@@ -253,8 +248,8 @@ TEST(fvm_lowered, target_handles) {
make_cell_ball_and_3stick()
};
- EXPECT_EQ(cells[0].num_segments(), 2u);
- EXPECT_EQ(cells[1].num_segments(), 4u);
+ EXPECT_EQ(cells[0].num_branches(), 2u);
+ EXPECT_EQ(cells[1].num_branches(), 4u);
// (in increasing target order)
cells[0].place(mlocation{1, 0.4}, "expsyn");
@@ -397,30 +392,26 @@ TEST(fvm_lowered, derived_mechs) {
//
// 3. Cell with both test_kin1 and custom_kin1.
- std::vector<cable_cell> cells(3);
+ std::vector<cable_cell> cells;
+ cells.reserve(3);
for (int i = 0; i<3; ++i) {
- cable_cell& c = cells[i];
- c.add_soma(6.0);
- c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 100));
-
- c.segment(1)->set_compartments(4);
- for (auto& seg: c.segments()) {
- if (!seg->is_soma()) {
- seg->as_cable()->set_compartments(4);
- }
- switch (i) {
+ soma_cell_builder builder(6);
+ builder.add_branch(0, 100, 0.5, 0.5, 4, "dend");
+ auto cell = builder.make_cell();
+
+ switch (i) {
case 0:
- seg->add_mechanism("test_kin1");
+ cell.paint(reg::all(), "test_kin1");
break;
case 1:
- seg->add_mechanism("custom_kin1");
+ cell.paint(reg::all(), "custom_kin1");
break;
case 2:
- seg->add_mechanism("test_kin1");
- seg->add_mechanism("custom_kin1");
+ cell.paint(reg::all(), "test_kin1");
+ cell.paint(reg::all(), "custom_kin1");
break;
- }
}
+ cells.push_back(std::move(cell));
}
cable1d_recipe rec(cells);
@@ -506,11 +497,10 @@ TEST(fvm_lowered, read_valence) {
{
std::vector<cable_cell> cells(1);
- cable_cell& c = cells[0];
- auto soma = c.add_soma(6.0);
- soma->add_mechanism("test_ca_read_valence");
-
- cable1d_recipe rec(cells);
+ soma_cell_builder builder(6);
+ auto cell = builder.make_cell();
+ cell.paint("soma", "test_ca_read_valence");
+ cable1d_recipe rec({std::move(cell)});
rec.catalogue() = make_unit_test_catalogue();
fvm_cell fvcell(context);
@@ -529,13 +519,11 @@ TEST(fvm_lowered, read_valence) {
{
// Check ion renaming.
- std::vector<cable_cell> cells(1);
-
- cable_cell& c = cells[0];
- auto soma = c.add_soma(6.0);
- soma->add_mechanism("cr_read_valence");
-
- cable1d_recipe rec(cells);
+ soma_cell_builder builder(6);
+ auto cell = builder.make_cell();
+ cell.paint("soma", "cr_read_valence");
+ cable1d_recipe rec({std::move(cell)});
+ rec.catalogue() = make_unit_test_catalogue();
rec.catalogue() = make_unit_test_catalogue();
rec.catalogue().derive("na_read_valence", "test_ca_read_valence", {}, {{"ca", "na"}});
@@ -557,13 +545,13 @@ TEST(fvm_lowered, read_valence) {
// Test correct scaling of ionic currents in reading and writing
TEST(fvm_lowered, ionic_currents) {
- cable_cell c;
- auto soma = c.add_soma(6.0);
+ soma_cell_builder b(6);
// Mechanism parameter is in NMODL units, i.e. mA/cm².
const double jca = 1.5;
- soma->add_mechanism(mechanism_desc("fixed_ica_current").set("ica_density", jca));
+ mechanism_desc m1("fixed_ica_current");
+ m1["ica_density"] = jca;
// Mechanism models a well-mixed fixed-depth volume without replenishment,
// giving a linear response to ica over time.
@@ -573,9 +561,14 @@ TEST(fvm_lowered, ionic_currents) {
// with NMODL units: cai' [mM/ms]; ica [mA/cm²], giving coeff in [mol/cm/C].
const double coeff = 0.5;
- soma->add_mechanism(mechanism_desc("linear_ca_conc").set("coeff", coeff));
+ mechanism_desc m2("linear_ca_conc");
+ m2["coeff"] = coeff;
- cable1d_recipe rec(c);
+ auto c = b.make_cell();
+ c.paint("soma", m1);
+ c.paint("soma", m2);
+
+ cable1d_recipe rec(std::move(c));
rec.catalogue() = make_unit_test_catalogue();
execution_context context;
@@ -604,10 +597,9 @@ TEST(fvm_lowered, ionic_currents) {
// Test correct scaling of an ionic current updated via a point mechanism
TEST(fvm_lowered, point_ionic_current) {
- cable_cell c;
-
double r = 6.0; // [µm]
- c.add_soma(6.0);
+ soma_cell_builder b(r);
+ cable_cell c = b.make_cell();
double soma_area_m2 = 4*math::pi<double>*r*r*1e-12; // [m²]
@@ -673,14 +665,12 @@ TEST(fvm_lowered, weighted_write_ion) {
execution_context context;
- cable_cell c;
- c.add_soma(5);
+ soma_cell_builder b(5);
+ b.add_branch(0, 100, 0.5, 0.5, 1, "dend");
+ b.add_branch(1, 200, 0.5, 0.5, 1, "dend");
+ b.add_branch(1, 100, 0.5, 0.5, 1, "dend");
- c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 100));
- c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 200));
- c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.5, 0.5, 100));
-
- for (auto& s: c.segments()) s->set_compartments(1);
+ cable_cell c = b.make_cell();
const double con_int = 80;
const double con_ext = 120;
@@ -768,18 +758,21 @@ TEST(fvm_lowered, gj_coords_simple) {
gap_recipe rec;
std::vector<cable_cell> cells;
- cable_cell c, d;
- c.add_soma(2.1);
- c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.3, 0.2, 10));
- c.segment(1)->set_compartments(5);
- c.place(mlocation{1, 0.8}, gap_junction_site{});
- cells.push_back(std::move(c));
-
- d.add_soma(2.4);
- d.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.3, 0.2, 10));
- d.segment(1)->set_compartments(2);
- d.place(mlocation{1, 1}, gap_junction_site{});
- cells.push_back(std::move(d));
+ {
+ soma_cell_builder b(2.1);
+ b.add_branch(0, 10, 0.3, 0.2, 5, "dend");
+ auto c = b.make_cell();
+ c.place(mlocation{1, 0.8}, gap_junction_site{});
+ cells.push_back(std::move(c));
+ }
+
+ {
+ soma_cell_builder b(2.4);
+ b.add_branch(0, 10, 0.3, 0.2, 2, "dend");
+ auto c = b.make_cell();
+ c.place(mlocation{1, 1}, gap_junction_site{});
+ cells.push_back(std::move(c));
+ }
fvm_discretization D = fvm_discretize(cells, neuron_parameter_defaults);
@@ -840,56 +833,44 @@ TEST(fvm_lowered, gj_coords_complex) {
execution_context context;
fvm_cell fvcell(context);
- gap_recipe rec;
- cable_cell c0, c1, c2;
- std::vector<cable_cell> cells;
-
- // Make 3 cells
- c0.add_soma(2.1);
- c0.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.3, 0.2, 8));
- c0.segment(1)->set_compartments(4);
-
- c1.add_soma(1.4);
- c1.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.3, 0.5, 12));
- c1.segment(1)->set_compartments(6);
- c1.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.3, 0.2, 9));
- c1.segment(2)->set_compartments(3);
- c1.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.2, 0.2, 15));
- c1.segment(3)->set_compartments(5);
-
- c2.add_soma(2.9);
- c2.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 0.3, 0.5, 4));
- c2.segment(1)->set_compartments(2);
- c2.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.4, 0.2, 6));
- c2.segment(2)->set_compartments(2);
- c2.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 0.1, 0.2, 8));
- c2.segment(3)->set_compartments(2);
- c2.add_cable(2, make_segment<cable_segment>(section_kind::dendrite, 0.2, 0.2, 4));
- c2.segment(4)->set_compartments(2);
- c2.add_cable(2, make_segment<cable_segment>(section_kind::dendrite, 0.2, 0.2, 4));
- c2.segment(5)->set_compartments(2);
-
// Add 5 gap junctions
+ soma_cell_builder b0(2.1);
+ b0.add_branch(0, 8, 0.3, 0.2, 4, "dend");
+
+ auto c0 = b0.make_cell();
c0.place(mlocation{1, 1}, gap_junction_site{});
c0.place(mlocation{1, 0.5}, gap_junction_site{});
+ soma_cell_builder b1(1.4);
+ b1.add_branch(0, 12, 0.3, 0.5, 6, "dend");
+ b1.add_branch(1, 9, 0.3, 0.2, 3, "dend");
+ b1.add_branch(1, 5, 0.2, 0.2, 5, "dend");
+
+ auto c1 = b1.make_cell();
c1.place(mlocation{2, 1}, gap_junction_site{});
c1.place(mlocation{1, 1}, gap_junction_site{});
c1.place(mlocation{1, 0.45}, gap_junction_site{});
c1.place(mlocation{1, 0.1}, gap_junction_site{});
+ soma_cell_builder b2(2.9);
+ b2.add_branch(0, 4, 0.3, 0.5, 2, "dend");
+ b2.add_branch(1, 6, 0.4, 0.2, 2, "dend");
+ b2.add_branch(1, 8, 0.1, 0.2, 2, "dend");
+ b2.add_branch(2, 4, 0.2, 0.2, 2, "dend");
+ b2.add_branch(2, 4, 0.2, 0.2, 2, "dend");
+
+ auto c2 = b2.make_cell();
c2.place(mlocation{1, 0.5}, gap_junction_site{});
c2.place(mlocation{4, 1}, gap_junction_site{});
c2.place(mlocation{2, 1}, gap_junction_site{});
- cells.push_back(std::move(c0));
- cells.push_back(std::move(c1));
- cells.push_back(std::move(c2));
+ std::vector<cable_cell> cells{std::move(c0), std::move(c1), std::move(c2)};
std::vector<fvm_index_type> cell_to_intdom;
std::vector<cell_gid_type> gids = {0, 1, 2};
+ gap_recipe rec;
fvcell.fvm_intdom(rec, gids, cell_to_intdom);
fvm_discretization D = fvm_discretize(cells, neuron_parameter_defaults);
@@ -955,8 +936,7 @@ TEST(fvm_lowered, cell_group_gj) {
// Make 20 cells
for (unsigned i = 0; i < 20; i++) {
- cable_cell c;
- c.add_soma(2.1);
+ cable_cell c = soma_cell_builder(2.1).make_cell();
if (i % 2 == 0) {
c.place(mlocation{0, 1}, gap_junction_site{});
}
diff --git a/test/unit/test_mc_cell_group.cpp b/test/unit/test_mc_cell_group.cpp
index 25415273f1d19c78a1fc0ec5710ff3eeb94f0c87..1be195302555b8cf917c8bc1c78031c92f489d42 100644
--- a/test/unit/test_mc_cell_group.cpp
+++ b/test/unit/test_mc_cell_group.cpp
@@ -21,11 +21,13 @@ namespace {
}
cable_cell make_cell() {
- auto c = make_cell_ball_and_stick();
-
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 200, 0.5, 0.5, 101, "dend");
+ cable_cell c = builder.make_cell();
+ c.paint("soma", "hh");
+ c.paint("dend", "pas");
+ c.place(mlocation{1,1}, i_clamp{5, 80, 0.3});
c.place(mlocation{0, 0}, threshold_detector{0});
- c.segment(1)->set_compartments(101);
-
return c;
}
}
@@ -36,6 +38,7 @@ ACCESS_BIND(
&mc_cell_group::spike_sources_)
TEST(mc_cell_group, get_kind) {
+ auto x = make_cell();
mc_cell_group group{{0}, cable1d_recipe(make_cell()), lowered_cell()};
EXPECT_EQ(cell_kind::cable, group.get_cell_kind());
diff --git a/test/unit/test_mc_cell_group_gpu.cpp b/test/unit/test_mc_cell_group_gpu.cpp
index 05b3d68193969ea63461e956fc3be763948b26cc..1053b0b854b07eb9c51e162f26866efb9c503607 100644
--- a/test/unit/test_mc_cell_group_gpu.cpp
+++ b/test/unit/test_mc_cell_group_gpu.cpp
@@ -19,11 +19,13 @@ namespace {
}
cable_cell make_cell() {
- auto c = make_cell_ball_and_stick();
-
+ soma_cell_builder builder(12.6157/2.0);
+ builder.add_branch(0, 200, 0.5, 0.5, 101, "dend");
+ cable_cell c = builder.make_cell();
+ c.paint("soma", "hh");
+ c.paint("dend", "pas");
+ c.place(mlocation{1,1}, i_clamp{5, 80, 0.3});
c.place(mlocation{0, 0}, threshold_detector{0});
- c.segment(1)->set_compartments(101);
-
return c;
}
}
diff --git a/test/unit/test_morphology.cpp b/test/unit/test_morphology.cpp
index fddbd74434a83d13bf1981d9a664170c36b33f03..eb34f28939ec7574e02bbd9fa51d3cb4f6cf3c49 100644
--- a/test/unit/test_morphology.cpp
+++ b/test/unit/test_morphology.cpp
@@ -570,7 +570,7 @@ TEST(morphology, swc) {
EXPECT_EQ(31u, m.num_branches());
// Confirm that converting to a cable_cell generates the same number of branches.
- auto c = arb::make_cable_cell(m, {}, false);
- EXPECT_EQ(31u, c.num_segments());
+ auto c = arb::cable_cell(m, {}, false);
+ EXPECT_EQ(31u, c.num_branches());
}
diff --git a/test/unit/test_synapses.cpp b/test/unit/test_synapses.cpp
index 5101732d2edd991386c2ac9784191c99228d407f..0c10474c08d58fe8d2bd7b3af51c3f6f800fc52f 100644
--- a/test/unit/test_synapses.cpp
+++ b/test/unit/test_synapses.cpp
@@ -14,6 +14,7 @@
#include "util/maputil.hpp"
#include "util/range.hpp"
+#include "../common_cells.hpp"
#include "common.hpp"
#include "mech_private_field_access.hpp"
@@ -32,11 +33,7 @@ ACCESS_BIND(value_type* multicore::mechanism::*, vec_i_ptr, &multicore::mechanis
TEST(synapses, add_to_cell) {
using namespace arb;
- cable_cell cell;
-
- // Soma with diameter 12.6157 um and HH channel
- auto soma = cell.add_soma(12.6157/2.0);
- soma->add_mechanism("hh");
+ auto cell = make_cell_soma_only(false);
cell.place(mlocation{0, 0.1}, "expsyn");
cell.place(mlocation{0, 0.2}, "exp2syn");
@@ -58,7 +55,7 @@ TEST(synapses, add_to_cell) {
EXPECT_EQ(syns[2].mechanism.name(), "expsyn");
// adding a synapse to an invalid branch location should throw.
- EXPECT_THROW(cell.place(mlocation{1, 0.3}, "expsyn"), arb::cable_cell_error);
+ EXPECT_THROW(cell.place(mlocation{1, 0.3}, "expsyn"), std::runtime_error);
}
template <typename Seq>
diff --git a/test/validation/CMakeLists.txt b/test/validation/CMakeLists.txt
deleted file mode 100644
index 7d885e07b8964a9f74810f6e5db398f4d871ad62..0000000000000000000000000000000000000000
--- a/test/validation/CMakeLists.txt
+++ /dev/null
@@ -1,25 +0,0 @@
-set(validation_sources
- # unit tests
- validate_ball_and_stick.cpp
- validate_soma.cpp
- validate_kinetic.cpp
- validate_synapses.cpp
-
- # support code
- validation_data.cpp
- trace_analysis.cpp
-
- # unit test driver
- validate.cpp
-)
-
-add_executable(validate EXCLUDE_FROM_ALL ${validation_sources})
-add_dependencies(tests validate)
-
-target_compile_options(validate PRIVATE ${ARB_CXXOPT_ARCH})
-target_compile_definitions(validate PRIVATE "ARB_DATADIR=\"${ARB_VALIDATION_DATA_DIR}\"")
-target_link_libraries(validate PRIVATE gtest arbor arbor-sup ext-json)
-
-if(ARB_BUILD_VALIDATION_DATA)
- add_dependencies(validate validation_data)
-endif()
diff --git a/test/validation/convergence_test.hpp b/test/validation/convergence_test.hpp
deleted file mode 100644
index 189615bbc551637e1a362e9f0fdb803ee3e52967..0000000000000000000000000000000000000000
--- a/test/validation/convergence_test.hpp
+++ /dev/null
@@ -1,176 +0,0 @@
-#pragma once
-
-#include <iterator>
-#include <vector>
-
-#include <nlohmann/json.hpp>
-
-#include <arbor/sampling.hpp>
-#include <arbor/simple_sampler.hpp>
-#include <arbor/simulation.hpp>
-#include <arbor/schedule.hpp>
-#include <sup/path.hpp>
-
-#include "../gtest.h"
-
-#include "trace_analysis.hpp"
-#include "validation_data.hpp"
-
-namespace arb {
-
-struct probe_label {
- const char* label;
- cell_member_type probe_id;
-};
-
-/* Common functionality for testing convergence towards
- * a reference data set as some parameter of the model
- * is changed.
- *
- * Type parameter Param is the type of the parameter that
- * is changed between runs.
- */
-
-template <typename Param>
-class convergence_test_runner {
-private:
- std::string param_name_;
- bool run_validation_;
- nlohmann::json meta_;
- std::vector<probe_label> probe_labels_;
- std::map<std::string, trace_data<double>> ref_data_;
- std::vector<conv_entry<Param>> conv_tbl_;
-
-public:
- template <typename ProbeLabelSeq>
- convergence_test_runner(
- const std::string& param_name,
- const ProbeLabelSeq& probe_labels,
- const nlohmann::json& meta
- ):
- param_name_(param_name),
- run_validation_(false),
- meta_(meta)
- {
- using std::begin;
- using std::end;
-
- probe_labels_.assign(begin(probe_labels), end(probe_labels));
- }
-
- // Allow free access to JSON meta data attached to saved traces.
- nlohmann::json& metadata() { return meta_; }
-
- void load_reference_data(const sup::path& ref_path) {
- run_validation_ = false;
- try {
- ref_data_ = g_trace_io.load_traces(ref_path);
-
- run_validation_ = true;
- for (const auto& pl: probe_labels_) {
- if (!(ref_data_.count(pl.label)>0)) {
- run_validation_ = false;
- break;
- }
- }
-
- EXPECT_TRUE(run_validation_);
- }
- catch (std::runtime_error&) {
- ADD_FAILURE() << "failure loading reference data: " << ref_path;
- }
- }
-
- void run(simulation& sim, Param p, float sample_dt, float t_end, float dt, const std::vector<float>& excl) {
- struct sampler_state {
- sampler_association_handle h; // Keep these for clean up at end.
- const char* label;
- trace_data<double> trace;
- };
-
- // Attach new samplers to labelled probe ids.
- std::vector<sampler_state> samplers;
- samplers.reserve(probe_labels_.size());
-
- for (auto& pl: probe_labels_) {
- samplers.push_back({});
- auto& entry = samplers.back();
-
- entry.label = pl.label;
- entry.h = sim.add_sampler(one_probe(pl.probe_id), regular_schedule(sample_dt), simple_sampler<double>(entry.trace));
- }
-
- sim.run(t_end, dt);
-
- for (auto& entry: samplers) {
- std::string label = entry.label;
- const auto& trace = entry.trace;
-
- // Save trace.
- nlohmann::json trace_meta(meta_);
- trace_meta[param_name_] = p;
-
- g_trace_io.save_trace(label, trace, trace_meta);
-
- // Compute metreics.
- if (run_validation_) {
- double linf = linf_distance(ref_data_[label], trace, excl);
- auto pd = peak_delta(trace, ref_data_[label]);
-
- conv_tbl_.push_back({label, p, linf, pd});
- }
- }
-
- // Remove added samplers.
- for (const auto& entry: samplers) {
- sim.remove_sampler(entry.h);
- }
- }
-
- void report() {
- if (run_validation_ && g_trace_io.verbose()) {
- // reorder to group by id
- std::stable_sort(conv_tbl_.begin(), conv_tbl_.end(),
- [](const auto& a, const auto& b) { return a.id<b.id; });
-
- report_conv_table(std::cout, conv_tbl_, param_name_);
- }
- }
-
- void assert_all_convergence() const {
- std::vector<conv_entry<Param>> with_label;
-
- for (const auto& pl: probe_labels_) {
- SCOPED_TRACE(pl.label);
-
- with_label.clear();
- for (const auto& e: conv_tbl_) {
- if (e.id==pl.label) {
- with_label.push_back(e);
- }
- }
-
- assert_convergence(with_label);
- }
- }
-};
-
-/*
- * Extract time points to exclude from current stimulus end-points.
- */
-
-inline std::vector<float> stimulus_ends(const cable_cell& c) {
- std::vector<float> ts;
-
- for (const auto& stimulus: c.stimuli()) {
- float t0 = stimulus.clamp.delay;
- float t1 = t0+stimulus.clamp.duration;
- ts.push_back(t0);
- ts.push_back(t1);
- }
-
- std::sort(ts.begin(), ts.end());
- return ts;
-}
-
-} // namespace arb
diff --git a/test/validation/interpolate.hpp b/test/validation/interpolate.hpp
deleted file mode 100644
index bb6ec2423099e39737c175ed945cb42452507431..0000000000000000000000000000000000000000
--- a/test/validation/interpolate.hpp
+++ /dev/null
@@ -1,52 +0,0 @@
-#pragma once
-
-#include <cmath>
-
-#include <arbor/util/compat.hpp>
-
-template <typename T, typename U>
-inline T lerp(T a, T b, U u) {
- return compat::fma(T(u), b, compat::fma(T(-u), a, a));
-}
-
-// Piece-wise linear interpolation across a sequence of points (u_i, x_i),
-// monotonically increasing in u.
-//
-// Parameters get_u and get_x provide the accessors for the point sequence;
-// consider moving to structured bindings in C++17 instead.
-
-template <typename U, typename Seq, typename GetU, typename GetX>
-auto pw_linear_interpolate(U u, const Seq& seq, GetU get_u, GetX get_x) {
- using std::begin;
- using std::end;
- using value_type = decltype(get_x(*begin(seq)));
-
- auto i = begin(seq);
- auto e = end(seq);
-
- if (i==e) {
- return value_type(NAN);
- }
-
- auto u0 = get_u(*i);
- auto x0 = get_x(*i);
-
- if (u<u0) {
- return x0;
- }
-
- while (++i!=e) {
- auto u1 = get_u(*i);
- auto x1 = get_x(*i);
-
- if (u<u1) {
- return lerp(x0, x1, (u-u0)/(u1-u0));
- }
-
- u0 = u1;
- x0 = x1;
- }
-
- return x0;
-}
-
diff --git a/test/validation/make_image.sh b/test/validation/make_image.sh
deleted file mode 100644
index c3b4a91096bba267810067fceebfe73c8ee271d8..0000000000000000000000000000000000000000
--- a/test/validation/make_image.sh
+++ /dev/null
@@ -1,5 +0,0 @@
-dot cell.dot -Tpdf -o cell.pdf
-# uses display from imagemagik
-display cell.pdf
-# or just call open on os x
-#open cell.pdf
diff --git a/test/validation/plot.py b/test/validation/plot.py
deleted file mode 100644
index 754eb6f2cfcf336215d5b1c1ba440138833ad3c8..0000000000000000000000000000000000000000
--- a/test/validation/plot.py
+++ /dev/null
@@ -1,20 +0,0 @@
-from matplotlib import pyplot
-import numpy as np
-
-raw = np.fromfile("den_syn.txt", sep=" ")
-n = raw.size/3
-data = raw.reshape(n,3)
-
-t = data[:, 0]
-soma = data[:, 1]
-dend = data[:, 2]
-
-pyplot.plot(t, soma)
-pyplot.plot(t, dend)
-
-pyplot.xlabel('time (ms)')
-pyplot.ylabel('mV')
-pyplot.xlim([t[0], t[n-1]])
-pyplot.grid()
-pyplot.show()
-
diff --git a/test/validation/trace_analysis.cpp b/test/validation/trace_analysis.cpp
deleted file mode 100644
index ccf71a4284842d46bf70cf553800c4390b389bc2..0000000000000000000000000000000000000000
--- a/test/validation/trace_analysis.cpp
+++ /dev/null
@@ -1,125 +0,0 @@
-#include <algorithm>
-#include <fstream>
-#include <string>
-
-#include <nlohmann/json.hpp>
-
-#include "../gtest.h"
-
-#include <arbor/math.hpp>
-#include <arbor/simple_sampler.hpp>
-#include <arbor/util/optional.hpp>
-
-#include "interpolate.hpp"
-#include "trace_analysis.hpp"
-
-namespace arb {
-
-struct trace_interpolant {
- trace_interpolant(const trace_data<double>& trace): trace_(trace) {}
-
- double operator()(float t) const {
- return pw_linear_interpolate(t, trace_,
- [](auto& entry) { return entry.t; },
- [](auto& entry) { return entry.v; });
- }
-
- const trace_data<double>& trace_;
-};
-
-double linf_distance(const trace_data<double>& u, const trace_data<double>& r) {
- trace_interpolant f{r};
-
- double linf = 0;
- for (auto entry: u) {
- linf = std::max(linf, std::abs(entry.v-f(entry.t)));
- }
-
- return linf;
-}
-
-// Compute linf distance as above, but excluding sample points that lie
-// near points in `excl`.
-//
-// `excl` contains the times to exclude, in ascending order.
-
-double linf_distance(const trace_data<double>& u, const trace_data<double>& ref, const std::vector<float>& excl) {
- trace_interpolant f{ref};
-
- trace_data<double> reduced;
- unsigned nexcl = excl.size();
- unsigned ei = 0;
-
- unsigned nu = u.size();
- unsigned ui = 0;
-
- while (ei<nexcl && ui<nu) {
- float t = excl[ei++];
-
- unsigned uj = ui;
- while (uj<nu && u[uj].t<t) ++uj;
-
- // include points up to and including uj-2, and then proceed from point uj+1,
- // excluding the two points closest to the discontinuity.
-
- for (unsigned k = ui; k+1<uj; ++k) {
- reduced.push_back(u[k]);
- }
- ui = uj+1;
- }
-
- for (auto k = ui; k<nu; ++k) {
- reduced.push_back(u[k]);
- }
-
- return linf_distance(reduced, ref);
-}
-
-std::vector<trace_peak> local_maxima(const trace_data<double>& u) {
- std::vector<trace_peak> peaks;
- if (u.size()<2) return peaks;
-
- int s_prev = math::signum(u[1].v-u[0].v);
- std::size_t i_start = 0;
-
- for (std::size_t i = 2; i<u.size()-1; ++i) {
- int s = math::signum(u[i].v-u[i-1].v);
- if (s_prev==1 && s==-1) {
- // found peak between i_start and i,
- // observerd peak value at i-1.
- float t0 = u[i_start].t;
- float t1 = u[i].t;
-
- peaks.push_back({(t0+t1)/2, u[i-1].v, (t1-t0)/2});
- }
-
- if (s!=0) {
- s_prev = s;
- if (s_prev>0) {
- i_start = i-1;
- }
- }
- }
- return peaks;
-}
-
-util::optional<trace_peak> peak_delta(const trace_data<double>& a, const trace_data<double>& b) {
- auto p = local_maxima(a);
- auto q = local_maxima(b);
-
- if (p.size()!=q.size() || p.empty()) return util::nullopt;
-
- auto max_delta = p[0]-q[0];
-
- for (std::size_t i = 1u; i<p.size(); ++i) {
- auto delta = p[i]-q[i];
- // pick maximum delta by greatest lower bound on dt
- if (std::abs(delta.t)-delta.t_err>std::abs(max_delta.t)-max_delta.t_err) {
- max_delta = delta;
- }
- }
- return max_delta;
-}
-
-} // namespace arb
-
diff --git a/test/validation/trace_analysis.hpp b/test/validation/trace_analysis.hpp
deleted file mode 100644
index 92b8d04811f5666bd7ec5d932b0b7544944b3545..0000000000000000000000000000000000000000
--- a/test/validation/trace_analysis.hpp
+++ /dev/null
@@ -1,109 +0,0 @@
-#pragma once
-
-#include <cfloat>
-#include <vector>
-
-#include "../gtest.h"
-
-#include <arbor/simple_sampler.hpp>
-#include <arbor/util/optional.hpp>
-
-#include "util.hpp"
-
-namespace arb {
-
-/* Trace data comparison */
-
-// Compute max |v_i - f(t_i)| where (t, v) is the
-// first trace `u` and f is the piece-wise linear interpolant
-// of the second trace `r`.
-
-double linf_distance(const trace_data<double>& u, const trace_data<double>& r);
-
-// Compute linf distance as above, excluding samples in the first trace
-// near times given in `excl`, monotonically increasing.
-
-double linf_distance(const trace_data<double>& u, const trace_data<double>& r, const std::vector<float>& excl);
-
-// Find local maxima (peaks) in a trace, excluding end points.
-
-struct trace_peak {
- float t;
- double v;
- float t_err;
-
- friend trace_peak operator-(trace_peak x, trace_peak y) {
- return {x.t-y.t, x.v-y.v, x.t_err+y.t_err};
- }
-};
-
-std::vector<trace_peak> local_maxima(const trace_data<double>& u);
-
-// Compare differences in peak times across two traces.
-// Returns largest magnitute displacement between peaks,
-// together with a sampling error bound, or `nothing`
-// if the number of peaks differ.
-
-util::optional<trace_peak> peak_delta(const trace_data<double>& a, const trace_data<double>& b);
-
-// Record for error data for convergence testing.
-// Only linf and peak_delta are used for convergence testing below;
-// if and param are for record keeping in the validation test itself.
-
-template <typename Param>
-struct conv_entry {
- std::string id;
- Param param;
- double linf;
- util::optional<trace_peak> peak_delta;
-};
-
-template <typename Param>
-using conv_data = std::vector<conv_entry<Param>>;
-
-// Assert error convergence (gtest).
-
-template <typename ConvEntrySeq>
-void assert_convergence(const ConvEntrySeq& cs) {
- if (size(cs)==0) return;
-
- auto tbound = [](trace_peak p) { return std::abs(p.t)+p.t_err; };
- float peak_dt_bound = INFINITY;
-
- for (auto pi = std::begin(cs); std::next(pi)!=std::end(cs); ++pi) {
- const auto& p = *pi;
- const auto& c = *std::next(pi);
-
- EXPECT_LE(c.linf, p.linf) << "L∞ error increase";
-
- if (!c.peak_delta) {
- EXPECT_FALSE(p.peak_delta) << "divergence in peak count";
- }
- else {
- double t = std::abs(c.peak_delta->t);
- double t_limit = c.peak_delta->t_err+peak_dt_bound;
-
- EXPECT_LE(t, t_limit) << "divergence in max peak displacement";
-
- peak_dt_bound = std::min(peak_dt_bound, tbound(*c.peak_delta));
- }
- }
-}
-
-// Report table of convergence results.
-
-template <typename ConvEntrySeq>
-void report_conv_table(std::ostream& out, const ConvEntrySeq& tbl, const std::string& param_name) {
- out << "id," << param_name << ",linf,peak_dt,peak_dt_err\n";
- for (const auto& c: tbl) {
- out << c.id << "," << c.param << "," << c.linf << ",";
- if (c.peak_delta) {
- out << c.peak_delta->t << "," << c.peak_delta->t_err << "\n";
- }
- else {
- out << "NA,NA\n";
- }
- }
-}
-
-} // namespace arb
diff --git a/test/validation/util.hpp b/test/validation/util.hpp
deleted file mode 100644
index af02639276c7c58ec8ec0ed7215efd8455c28d78..0000000000000000000000000000000000000000
--- a/test/validation/util.hpp
+++ /dev/null
@@ -1,22 +0,0 @@
-#pragma once
-
-// Simple helper utilities for validation tests.
-
-#include <sstream>
-#include <string>
-
-#include <arbor/common_types.hpp>
-
-template <typename T, std::size_t N>
-constexpr std::size_t size(T (&)[N]) noexcept {
- return N;
-}
-
-template <typename X>
-constexpr std::size_t size(const X& x) { return x.size(); }
-
-inline std::string to_string(arb::backend_kind kind) {
- std::stringstream out;
- out << kind;
- return out.str();
-}
diff --git a/test/validation/validate.cpp b/test/validation/validate.cpp
deleted file mode 100644
index 082937079028043123ac89a32707531c028a4415..0000000000000000000000000000000000000000
--- a/test/validation/validate.cpp
+++ /dev/null
@@ -1,80 +0,0 @@
-#include <cstring>
-#include <iostream>
-#include <sstream>
-#include <string>
-#include <exception>
-
-#include <sup/tinyopt.hpp>
-
-#include "../gtest.h"
-
-#include "validation_data.hpp"
-
-using namespace arb;
-
-const char* usage_str =
-"[OPTION]...\n"
-"\n"
-" -v, --verbose Print results to stdout\n"
-" -o, --output=FILE Save traces from simulations to FILE\n"
-" -p, --path=DIR Look for validation reference data in DIR\n"
-" -m, --max-comp=N Run convergence tests to a maximum of N\n"
-" compartments per segment\n"
-" -d, --min-dt=DT Run convergence tests with a minimum timestep DT [ms]\n"
-" -s, --sample-dt=DT Sample rate for simulations [ms]\n"
-" -h, --help Display usage information and exit\n"
-"\n"
-"Validation data is searched by default in the directory specified by\n"
-"ARB_DATADIR at compile time. If ARB_DATADIR does not correspond to a\n"
-"directory, the tests will try the paths './validation/data' and\n"
-"'../validation/data'. This default path can be overridden with the\n"
-"ARB_DATADIR environment variable, or with the -p command-line option.\n";
-
-int main(int argc, char **argv) {
- using to::parse_opt;
-
- ::testing::InitGoogleTest(&argc, argv);
-
- try {
- auto arg = argv+1;
- while (*arg) {
- if (auto o = parse_opt<std::string>(arg, 'p', "path")) {
- g_trace_io.set_datadir(*o);
- }
- else if (auto o = parse_opt<std::string>(arg, 'o', "output")) {
- g_trace_io.set_output(*o);
- }
- else if (auto o = parse_opt<int>(arg, 'm', "max-comp")) {
- g_trace_io.set_max_ncomp(*o);
- }
- else if (auto o = parse_opt<float>(arg, 'd', "min-dt")) {
- g_trace_io.set_min_dt(*o);
- }
- else if (auto o = parse_opt<float>(arg, 's', "sample-dt")) {
- g_trace_io.set_sample_dt(*o);
- }
- else if (auto o = parse_opt(arg, 'v', "verbose")) {
- g_trace_io.set_verbose(true);
- }
- else if (auto o = parse_opt(arg, 'h', "help")) {
- to::usage(argv[0], usage_str);
- return 0;
- }
- else {
- throw to::parse_opt_error(*arg, "unrecognized option");
- }
- }
-
- return RUN_ALL_TESTS();
- }
- catch (to::parse_opt_error& e) {
- to::usage(argv[0], usage_str, e.what());
- return 1;
- }
- catch (std::exception& e) {
- std::cerr << "caught exception: " << e.what() << "\n";
- return 2;
- }
-
- return 0;
-}
diff --git a/test/validation/validate_ball_and_stick.cpp b/test/validation/validate_ball_and_stick.cpp
deleted file mode 100644
index 4f60de9e945b0b568d5d4fbb1cef29dbe32c563b..0000000000000000000000000000000000000000
--- a/test/validation/validate_ball_and_stick.cpp
+++ /dev/null
@@ -1,258 +0,0 @@
-#include <iostream>
-
-#include <nlohmann/json.hpp>
-
-#include <arbor/common_types.hpp>
-#include <arbor/context.hpp>
-#include <arbor/domain_decomposition.hpp>
-#include <arbor/context.hpp>
-#include <arbor/load_balance.hpp>
-#include <arbor/cable_cell.hpp>
-#include <arbor/recipe.hpp>
-#include <arbor/simple_sampler.hpp>
-#include <arbor/simulation.hpp>
-#include <sup/path.hpp>
-
-#include "../common_cells.hpp"
-#include "../simple_recipes.hpp"
-
-#include "convergence_test.hpp"
-#include "trace_analysis.hpp"
-#include "util.hpp"
-#include "validation_data.hpp"
-
-#include "../gtest.h"
-
-using namespace arb;
-
-struct probe_point {
- const char* label;
- mlocation where;
-};
-
-template <typename ProbePointSeq>
-void run_ncomp_convergence_test(
- const char* model_name,
- const sup::path& ref_data_path,
- context& context,
- const cable_cell& c,
- ProbePointSeq& probe_points,
- float t_end=100.f)
-{
- using namespace arb;
-
- auto max_ncomp = g_trace_io.max_ncomp();
- auto dt = g_trace_io.min_dt();
- auto sample_dt = g_trace_io.sample_dt();
-
- nlohmann::json meta = {
- {"name", "membrane voltage"},
- {"model", model_name},
- {"dt", dt},
- {"sim", "arbor"},
- {"units", "mV"},
- {"backend_kind", (has_gpu(context)? "gpu": "multicore")}
- };
-
- auto exclude = stimulus_ends(c);
-
- auto n_probe = size(probe_points);
- std::vector<probe_label> plabels;
- plabels.reserve(n_probe);
- for (unsigned i = 0; i<n_probe; ++i) {
- plabels.push_back(probe_label{probe_points[i].label, {0u, i}});
- }
-
- convergence_test_runner<int> runner("ncomp", plabels, meta);
- runner.load_reference_data(ref_data_path);
-
- for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) {
- for (auto& seg: c.segments()) {
- if (!seg->is_soma()) {
- seg->set_compartments(ncomp);
- }
- }
- cable1d_recipe rec{c};
- for (const auto& p: probe_points) {
- rec.add_probe(0, 0, cell_probe_address{p.where, cell_probe_address::membrane_voltage});
- }
-
- auto decomp = partition_load_balance(rec, context);
- simulation sim(rec, decomp, context);
-
- runner.run(sim, ncomp, sample_dt, t_end, dt, exclude);
- }
- runner.report();
- runner.assert_all_convergence();
-}
-
-void validate_ball_and_stick(context& ctx) {
- using namespace arb;
-
- cable_cell c = make_cell_ball_and_stick();
- probe_point points[] = {
- {"soma.mid", {0u, 0.5}},
- {"dend.mid", {1u, 0.5}},
- {"dend.end", {1u, 1.0}}
- };
-
- run_ncomp_convergence_test(
- "ball_and_stick",
- "neuron_ball_and_stick.json",
- ctx,
- c,
- points);
-}
-
-void validate_ball_and_taper(context& ctx) {
- using namespace arb;
-
- cable_cell c = make_cell_ball_and_taper();
- probe_point points[] = {
- {"soma.mid", {0u, 0.5}},
- {"taper.mid", {1u, 0.5}},
- {"taper.end", {1u, 1.0}}
- };
-
- run_ncomp_convergence_test(
- "ball_and_taper",
- "neuron_ball_and_taper.json",
- ctx,
- c,
- points);
-}
-
-void validate_ball_and_3stick(context& ctx) {
- using namespace arb;
-
- cable_cell c = make_cell_ball_and_3stick();
- probe_point points[] = {
- {"soma.mid", {0u, 0.5}},
- {"dend1.mid", {1u, 0.5}},
- {"dend1.end", {1u, 1.0}},
- {"dend2.mid", {2u, 0.5}},
- {"dend2.end", {2u, 1.0}},
- {"dend3.mid", {3u, 0.5}},
- {"dend3.end", {3u, 1.0}}
- };
-
- run_ncomp_convergence_test(
- "ball_and_3stick",
- "neuron_ball_and_3stick.json",
- ctx,
- c,
- points);
-}
-
-void validate_rallpack1(context& ctx) {
- using namespace arb;
-
- cable_cell c = make_cell_simple_cable();
- probe_point points[] = {
- {"cable.x0.0", {1u, 0.0}},
- {"cable.x0.3", {1u, 0.3}},
- {"cable.x1.0", {1u, 1.0}}
- };
-
- run_ncomp_convergence_test(
- "rallpack1",
- "numeric_rallpack1.json",
- ctx,
- c,
- points,
- 250.f);
-}
-
-void validate_ball_and_squiggle(context& ctx) {
- using namespace arb;
-
- cable_cell c = make_cell_ball_and_squiggle();
- probe_point points[] = {
- {"soma.mid", {0u, 0.5}},
- {"dend.mid", {1u, 0.5}},
- {"dend.end", {1u, 1.0}}
- };
-
-#if 0
- // *temporarily* disabled: compartment division policy will
- // be moved into backend policy classes.
-
- run_ncomp_convergence_test<lowered_cell_div<div_compartment_sampler>>(
- "ball_and_squiggle_sampler",
- "neuron_ball_and_squiggle.json",
- c,
- samplers);
-#endif
-
- run_ncomp_convergence_test(
- "ball_and_squiggle_integrator",
- "neuron_ball_and_squiggle.json",
- ctx,
- c,
- points);
-}
-
-TEST(ball_and_stick, neuron_ref) {
- proc_allocation resources;
- {
- auto ctx = make_context(resources);
- validate_ball_and_stick(ctx);
- }
- if (resources.has_gpu()) {
- resources.gpu_id = -1;
- auto ctx = make_context(resources);
- validate_ball_and_stick(ctx);
- }
-}
-
-TEST(ball_and_taper, neuron_ref) {
- proc_allocation resources;
- {
- auto ctx = make_context(resources);
- validate_ball_and_taper(ctx);
- }
- if (resources.has_gpu()) {
- resources.gpu_id = -1;
- auto ctx = make_context(resources);
- validate_ball_and_taper(ctx);
- }
-}
-
-TEST(ball_and_3stick, neuron_ref) {
- proc_allocation resources;
- {
- auto ctx = make_context(resources);
- validate_ball_and_3stick(ctx);
- }
- if (resources.has_gpu()) {
- resources.gpu_id = -1;
- auto ctx = make_context(resources);
- validate_ball_and_3stick(ctx);
- }
-}
-
-TEST(rallpack1, numeric_ref) {
- proc_allocation resources;
- {
- auto ctx = make_context(resources);
- validate_rallpack1(ctx);
- }
- if (resources.has_gpu()) {
- resources.gpu_id = -1;
- auto ctx = make_context(resources);
- validate_rallpack1(ctx);
- }
-}
-
-TEST(ball_and_squiggle, neuron_ref) {
- proc_allocation resources;
- {
- auto ctx = make_context(resources);
- validate_ball_and_squiggle(ctx);
- }
- if (resources.has_gpu()) {
- resources.gpu_id = -1;
- auto ctx = make_context(resources);
- validate_ball_and_squiggle(ctx);
- }
-}
diff --git a/test/validation/validate_kinetic.cpp b/test/validation/validate_kinetic.cpp
deleted file mode 100644
index e855d02187d5f203dc74a6108255112fe37b8f0b..0000000000000000000000000000000000000000
--- a/test/validation/validate_kinetic.cpp
+++ /dev/null
@@ -1,134 +0,0 @@
-#include "../gtest.h"
-
-#include <string>
-
-#include <nlohmann/json.hpp>
-
-#include <arbor/context.hpp>
-#include <arbor/common_types.hpp>
-#include <arbor/domain_decomposition.hpp>
-#include <arbor/load_balance.hpp>
-#include <arbor/cable_cell.hpp>
-#include <arbor/recipe.hpp>
-#include <arbor/simple_sampler.hpp>
-#include <arbor/simulation.hpp>
-
-#include "../common_cells.hpp"
-#include "../simple_recipes.hpp"
-
-#include "convergence_test.hpp"
-#include "trace_analysis.hpp"
-#include "util.hpp"
-#include "validation_data.hpp"
-
-void run_kinetic_dt(
- const arb::context& context,
- arb::cable_cell& c,
- arb::cell_probe_address probe,
- float t_end,
- nlohmann::json meta,
- const std::string& ref_file)
-{
- using namespace arb;
-
- float sample_dt = g_trace_io.sample_dt();
-
- cable1d_recipe rec{c};
- rec.add_probe(0, 0, probe);
-
- probe_label plabels[1] = {{"soma.mid", {0u, 0u}}};
-
- meta["sim"] = "arbor";
- meta["backend_kind"] = arb::has_gpu(context)? "gpu": "multicore";
-
- convergence_test_runner<float> runner("dt", plabels, meta);
- runner.load_reference_data(ref_file);
-
- auto decomp = partition_load_balance(rec, context);
- simulation sim(rec, decomp, context);
-
- auto exclude = stimulus_ends(c);
-
- // use dt = 0.05, 0.02, 0.01, 0.005, 0.002, ...
- double max_oo_dt = std::round(1.0/g_trace_io.min_dt());
- for (double base = 100; ; base *= 10) {
- for (double multiple: {5., 2., 1.}) {
- double oo_dt = base/multiple;
- if (oo_dt>max_oo_dt) goto end;
-
- sim.reset();
- float dt = float(1./oo_dt);
- runner.run(sim, dt, sample_dt, t_end, dt, exclude);
- }
- }
-
-end:
- runner.report();
- runner.assert_all_convergence();
-}
-
-void validate_kinetic_kin1(const arb::context& ctx) {
- using namespace arb;
-
- // 20 µm diameter soma with single mechanism, current probe
- cable_cell c;
- auto soma = c.add_soma(10);
- soma->add_mechanism("test_kin1");
- cell_probe_address probe{{0, 0.5}, cell_probe_address::membrane_current};
-
- nlohmann::json meta = {
- {"model", "test_kin1"},
- {"name", "membrane current"},
- {"units", "nA"}
- };
-
- run_kinetic_dt(ctx, c, probe, 100.f, meta, "numeric_kin1.json");
-}
-
-void validate_kinetic_kinlva(const arb::context& ctx) {
- using namespace arb;
-
- // 20 µm diameter soma with single mechanism, current probe
- cable_cell c;
- auto soma = c.add_soma(10);
- c.place(mlocation{0,0.5}, {20., 130., -0.025});
- soma->add_mechanism("test_kinlva");
- cell_probe_address probe{{0, 0.5}, cell_probe_address::membrane_voltage};
-
- nlohmann::json meta = {
- {"model", "test_kinlva"},
- {"name", "membrane voltage"},
- {"units", "mV"}
- };
-
- run_kinetic_dt(ctx, c, probe, 300.f, meta, "numeric_kinlva.json");
-}
-
-
-using namespace arb;
-
-TEST(kinetic, kin1_numeric_ref) {
- proc_allocation resources;
- {
- auto ctx = make_context(resources);
- validate_kinetic_kin1(ctx);
- }
- if (resources.has_gpu()) {
- resources.gpu_id = -1;
- auto ctx = make_context(resources);
- validate_kinetic_kin1(ctx);
- }
-}
-
-TEST(kinetic, kinlva_numeric_ref) {
- proc_allocation resources;
- {
- auto ctx = make_context(resources);
- validate_kinetic_kinlva(ctx);
- }
- if (resources.has_gpu()) {
- resources.gpu_id = -1;
- auto ctx = make_context(resources);
- validate_kinetic_kinlva(ctx);
- }
-}
diff --git a/test/validation/validate_soma.cpp b/test/validation/validate_soma.cpp
deleted file mode 100644
index 6d5c9c6db7e4697b1b8b3cbb03fcb96d1baf6409..0000000000000000000000000000000000000000
--- a/test/validation/validate_soma.cpp
+++ /dev/null
@@ -1,78 +0,0 @@
-#include <nlohmann/json.hpp>
-
-#include <arbor/common_types.hpp>
-#include <arbor/context.hpp>
-#include <arbor/domain_decomposition.hpp>
-#include <arbor/load_balance.hpp>
-#include <arbor/cable_cell.hpp>
-#include <arbor/recipe.hpp>
-#include <arbor/simple_sampler.hpp>
-#include <arbor/simulation.hpp>
-
-#include "../common_cells.hpp"
-#include "../simple_recipes.hpp"
-
-#include "convergence_test.hpp"
-#include "trace_analysis.hpp"
-#include "util.hpp"
-#include "validation_data.hpp"
-
-#include "../gtest.h"
-
-using namespace arb;
-
-void validate_soma(const context& context) {
- float sample_dt = g_trace_io.sample_dt();
-
- cable_cell c = make_cell_soma_only();
-
- cable1d_recipe rec{c};
- rec.add_probe(0, 0, cell_probe_address{{0, 0.5}, cell_probe_address::membrane_voltage});
- probe_label plabels[1] = {{"soma.mid", {0u, 0u}}};
-
- auto decomp = partition_load_balance(rec, context);
- simulation sim(rec, decomp, context);
-
- nlohmann::json meta = {
- {"name", "membrane voltage"},
- {"model", "soma"},
- {"sim", "arbor"},
- {"units", "mV"},
- {"backend_kind", has_gpu(context)? "gpu": "multicore"}
- };
-
- convergence_test_runner<float> runner("dt", plabels, meta);
- runner.load_reference_data("numeric_soma.json");
-
- float t_end = 100.f;
-
- // use dt = 0.05, 0.02, 0.01, 0.005, 0.002, ...
- double max_oo_dt = std::round(1.0/g_trace_io.min_dt());
- for (double base = 100; ; base *= 10) {
- for (double multiple: {5., 2., 1.}) {
- double oo_dt = base/multiple;
- if (oo_dt>max_oo_dt) goto end;
-
- sim.reset();
- float dt = float(1./oo_dt);
- runner.run(sim, dt, sample_dt, t_end, dt, {});
- }
- }
-end:
-
- runner.report();
- runner.assert_all_convergence();
-}
-
-TEST(soma, numeric_ref) {
- proc_allocation resources;
- {
- auto ctx = make_context(resources);
- validate_soma(ctx);
- }
- if (resources.has_gpu()) {
- resources.gpu_id = -1;
- auto ctx = make_context(resources);
- validate_soma(ctx);
- }
-}
diff --git a/test/validation/validate_synapses.cpp b/test/validation/validate_synapses.cpp
deleted file mode 100644
index daa3a209ff9fb98a53b88884b58aa4f6d98ff885..0000000000000000000000000000000000000000
--- a/test/validation/validate_synapses.cpp
+++ /dev/null
@@ -1,115 +0,0 @@
-#include <nlohmann/json.hpp>
-
-#include <arbor/context.hpp>
-#include <arbor/domain_decomposition.hpp>
-#include <arbor/load_balance.hpp>
-#include <arbor/cable_cell.hpp>
-#include <arbor/recipe.hpp>
-#include <arbor/simple_sampler.hpp>
-#include <arbor/simulation.hpp>
-#include <sup/path.hpp>
-
-
-#include "../gtest.h"
-
-#include "../common_cells.hpp"
-#include "../simple_recipes.hpp"
-
-#include "convergence_test.hpp"
-#include "trace_analysis.hpp"
-#include "util.hpp"
-#include "validation_data.hpp"
-
-using namespace arb;
-
-void run_synapse_test(
- const char* syn_type,
- const sup::path& ref_data_path,
- const context& context,
- float t_end=70.f,
- float dt=0.001)
-{
- auto max_ncomp = g_trace_io.max_ncomp();
- nlohmann::json meta = {
- {"name", "membrane voltage"},
- {"model", syn_type},
- {"sim", "arbor"},
- {"units", "mV"},
- {"backend_kind", arb::has_gpu(context)? "gpu": "multicore"}
- };
-
- cable_cell c = make_cell_ball_and_stick(false); // no stimuli
- mechanism_desc syn_default(syn_type);
- c.place(mlocation{1, 0.5}, syn_default);
-
- // injected spike events
- std::vector<spike_event> synthetic_events = {
- {{0u, 0u}, 10.0, 0.04},
- {{0u, 0u}, 20.0, 0.04},
- {{0u, 0u}, 40.0, 0.04}
- };
-
- // exclude points of discontinuity from linf analysis
- std::vector<float> exclude = {10.f, 20.f, 40.f};
-
- float sample_dt = g_trace_io.sample_dt();
- probe_label plabels[3] = {
- {"soma.mid", {0u, 0u}},
- {"dend.mid", {0u, 1u}},
- {"dend.end", {0u, 2u}}
- };
-
- convergence_test_runner<int> runner("ncomp", plabels, meta);
- runner.load_reference_data(ref_data_path);
-
- for (int ncomp = 10; ncomp<max_ncomp; ncomp*=2) {
- c.cable(1)->set_compartments(ncomp);
-
- cable1d_recipe rec{c};
- // soma.mid
- rec.add_probe(0, 0, cell_probe_address{{0, 0.5}, cell_probe_address::membrane_voltage});
- // dend.mid
- rec.add_probe(0, 0, cell_probe_address{{1, 0.5}, cell_probe_address::membrane_voltage});
- // dend.end
- rec.add_probe(0, 0, cell_probe_address{{1, 1.0}, cell_probe_address::membrane_voltage});
-
- auto decomp = partition_load_balance(rec, context);
- simulation sim(rec, decomp, context);
-
- sim.inject_events(synthetic_events);
-
- runner.run(sim, ncomp, sample_dt, t_end, dt, exclude);
- }
- runner.report();
- runner.assert_all_convergence();
-}
-
-TEST(simple_synapse, expsyn_neuron_ref) {
- proc_allocation resources;
- {
- auto ctx = make_context(resources);
- SCOPED_TRACE("expsyn-multicore");
- run_synapse_test("expsyn", "neuron_simple_exp_synapse.json", ctx);
- }
- if (resources.has_gpu()) {
- resources.gpu_id = -1;
- auto ctx = make_context(resources);
- SCOPED_TRACE("expsyn-gpu");
- run_synapse_test("expsyn", "neuron_simple_exp_synapse.json", ctx);
- }
-}
-
-TEST(simple_synapse, exp2syn_neuron_ref) {
- proc_allocation resources;
- {
- auto ctx = make_context(resources);
- SCOPED_TRACE("exp2syn-multicore");
- run_synapse_test("exp2syn", "neuron_simple_exp_synapse.json", ctx);
- }
- if (resources.has_gpu()) {
- resources.gpu_id = -1;
- auto ctx = make_context(resources);
- SCOPED_TRACE("exp2syn-gpu");
- run_synapse_test("exp2syn", "neuron_simple_exp_synapse.json", ctx);
- }
-}
diff --git a/test/validation/validation_data.cpp b/test/validation/validation_data.cpp
deleted file mode 100644
index 3ddbaaa7aadc002389203361de2db902648a8088..0000000000000000000000000000000000000000
--- a/test/validation/validation_data.cpp
+++ /dev/null
@@ -1,116 +0,0 @@
-#include <algorithm>
-#include <cstdlib>
-#include <fstream>
-#include <stdexcept>
-#include <string>
-
-#include <nlohmann/json.hpp>
-
-#include <arbor/simple_sampler.hpp>
-#include <sup/path.hpp>
-
-#include "trace_analysis.hpp"
-#include "validation_data.hpp"
-
-namespace arb {
-
-trace_io g_trace_io;
-
-#ifndef ARB_DATADIR
-#define ARB_DATADIR ""
-#endif
-
-sup::path trace_io::find_datadir() {
- // If environment variable is set, use that in preference.
-
- if (const char* env_path = std::getenv("ARB_DATADIR")) {
- return env_path;
- }
-
- // Otherwise try compile-time path ARB_DATADIR and hard-coded
- // relative paths below in turn, returning the first that
- // corresponds to an existing directory.
-
- const char* paths[] = {
- ARB_DATADIR,
- "./validation/data",
- "../validation/data"
- };
-
- std::error_code ec;
- for (auto p: paths) {
- if (sup::is_directory(p, ec)) {
- return p;
- }
- }
-
- // Otherwise set to empty path, and rely on command-line option.
- return "";
-}
-
-void trace_io::save_trace(const std::string& label, const trace_data<double>& data, const nlohmann::json& meta) {
- save_trace("time", label, data, meta);
-}
-
-void trace_io::save_trace(const std::string& abscissa, const std::string& label, const trace_data<double>& data, const nlohmann::json& meta) {
- using nlohmann::json;
-
- json j = meta;
- json& times = j["data"][abscissa];
- json& values = j["data"][label];
-
- for (const auto& e: data) {
- times.push_back(e.t);
- values.push_back(e.v);
- }
-
- jtraces_ += std::move(j);
-}
-
-template <typename Seq1, typename Seq2>
-static trace_data<double> zip_trace_data(const Seq1& ts, const Seq2& vs) {
- trace_data<double> trace;
-
- auto ti = std::begin(ts);
- auto te = std::end(ts);
- auto vi = std::begin(vs);
- auto ve = std::end(vs);
-
- while (ti!=te && vi!=ve) {
- trace.push_back({*ti++, *vi++});
- }
- return trace;
-}
-
-static void parse_trace_json(const nlohmann::json& j, std::map<std::string, trace_data<double>>& traces) {
- if (j.is_array()) {
- for (auto& i: j) parse_trace_json(i, traces);
- }
- else if (j.is_object() && j.count("data")>0 && j["data"].count("time")>0) {
- auto data = j["data"];
- auto time = data["time"].get<std::vector<float>>();
- for (const auto& p: nlohmann::json::iterator_wrapper(data)) {
- if (p.key()=="time") continue;
-
- traces[p.key()] = zip_trace_data(time, p.value().get<std::vector<double>>());
- }
- }
-}
-
-std::map<std::string, trace_data<double>> trace_io::load_traces(const sup::path& name) {
- sup::path file = datadir_/name;
- std::ifstream fid(file);
- if (!fid) {
- throw std::runtime_error("unable to load validation data: "+file.native());
- }
-
- nlohmann::json data;
- fid >> data;
-
- std::map<std::string, trace_data<double>> traces;
- parse_trace_json(data, traces);
- return traces;
-}
-
-} // namespace arb
-
diff --git a/test/validation/validation_data.hpp b/test/validation/validation_data.hpp
deleted file mode 100644
index bf827bf1198dd2a8272759ad53d4baa871834a5a..0000000000000000000000000000000000000000
--- a/test/validation/validation_data.hpp
+++ /dev/null
@@ -1,101 +0,0 @@
-#pragma once
-
-#include <fstream>
-#include <stdexcept>
-#include <string>
-#include <utility>
-
-#include <nlohmann/json.hpp>
-
-#include <arbor/simple_sampler.hpp>
-#include <sup/path.hpp>
-
-namespace arb {
-
-/*
- * Class manages input (loading and parsing) of JSON
- * reference traces, and output of saved traces from
- * the validation tests.
- *
- * Global object has paths initalised by the test
- * main() and will write all saved traces to a single
- * output JSON file, if specified.
- */
-
-// TODO: split simulation options from trace_io structure; rename
-// file to e.g. 'trace_io.hpp'
-
-class trace_io {
-public:
- // Try to find the data directory on construction.
-
- trace_io() {
- datadir_ = find_datadir();
- }
-
- void clear_traces() {
- jtraces_ = nlohmann::json::array();
- }
-
- void write_traces() {
- if (out_) {
- out_ << jtraces_;
- out_.close();
- }
- }
-
- void save_trace(const std::string& label, const trace_data<double>& data, const nlohmann::json& meta);
- void save_trace(const std::string& abscissa, const std::string& label, const trace_data<double>& data, const nlohmann::json& meta);
- std::map<std::string, trace_data<double>> load_traces(const sup::path& name);
-
- // common flags, options set by driver
-
- void set_verbose(bool v) { verbose_flag_ = v; }
- bool verbose() const { return verbose_flag_; }
-
- void set_max_ncomp(int ncomp) { max_ncomp_ = ncomp; }
- int max_ncomp() const { return max_ncomp_; }
-
- void set_min_dt(float dt) { min_dt_ = dt; }
- float min_dt() const { return min_dt_; }
-
- void set_sample_dt(float dt) { sample_dt_ = dt; }
- float sample_dt() const { return sample_dt_; }
-
- void set_datadir(const sup::path& dir) { datadir_ = dir; }
-
- void set_output(const sup::path& file) {
- out_.open(file);
- if (!out_) {
- throw std::runtime_error("unable to open file for writing");
- }
- }
-
- // Write traces on exit.
-
- ~trace_io() {
- if (out_) {
- write_traces();
- }
- }
-
-private:
- sup::path datadir_;
- std::ofstream out_;
- nlohmann::json jtraces_ = nlohmann::json::array();
- bool verbose_flag_ = false;
- int max_ncomp_ = 100;
- float min_dt_ = 0.001f;
- float sample_dt_ = 0.005f;
-
- // Returns value of ARB_DATADIR environment variable if set,
- // otherwise make a 'best-effort' search for the data directory,
- // starting with ARB_DATADIR preprocessor define if defined and
- // if the directory exists, or else try './validation/data'
- // and '../validation/data'.
- static sup::path find_datadir();
-};
-
-extern trace_io g_trace_io;
-
-} // namespace arb
diff --git a/todo.md b/todo.md
new file mode 100644
index 0000000000000000000000000000000000000000..150c42584f1639783a2b4284bb09aaaeed5e1fec
--- /dev/null
+++ b/todo.md
@@ -0,0 +1,5 @@
+Things to do in morophology part 3.
+
+* add ability to set the number of compartments per branch (removed by allowing only sample_tree descriptions)
+* remove `cable_cell::place(mlocation, ...)` with `cable_cell::place(locset, ...)`
+ * ensure consistent error messages in new code