diff --git a/CMakeLists.txt b/CMakeLists.txt
index 83163769b4dfae29a3b0eea18ae5e314918029f8..94e856778a4ffd74f7c6c757f60f8a6f330b60aa 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -372,9 +372,6 @@ add_subdirectory(test)
# self contained examples:
add_subdirectory(example)
-# lmorpho:
-add_subdirectory(lmorpho)
-
# html:
add_subdirectory(doc)
diff --git a/arbor/CMakeLists.txt b/arbor/CMakeLists.txt
index 978d3a30514d3d189d0d51e5f0f944a0a4034fcf..e7d9740df549268667da5a7c4b21c5fcf6b9c12a 100644
--- a/arbor/CMakeLists.txt
+++ b/arbor/CMakeLists.txt
@@ -28,9 +28,11 @@ set(arbor_sources
mechcat.cpp
memory/cuda_wrappers.cpp
memory/util.cpp
+ morph/morphology.cpp
+ morph/sample_tree.cpp
+ morph/primitives.cpp
merge_events.cpp
simulation.cpp
- morphology.cpp
partition_load_balance.cpp
profile/clock.cpp
profile/memory_meter.cpp
diff --git a/arbor/cable_cell.cpp b/arbor/cable_cell.cpp
index bea58ad8b7e86b35ba31dfa9908e7e379e307063..50ca10c05c968e6200e51d7823b402c618797608 100644
--- a/arbor/cable_cell.cpp
+++ b/arbor/cable_cell.cpp
@@ -1,8 +1,9 @@
#include <arbor/cable_cell.hpp>
+#include <arbor/morph/morphology.hpp>
#include <arbor/segment.hpp>
-#include <arbor/morphology.hpp>
#include "util/rangeutil.hpp"
+#include "util/span.hpp"
namespace arb {
@@ -163,43 +164,55 @@ value_type cable_cell::segment_mean_attenuation(
return 2*std::sqrt(math::pi<double>*tau_per_um*frequency)*length_factor; // [1/µm]
}
-
-// Construct cell from flat morphology specification.
-
cable_cell make_cable_cell(const morphology& morph, bool compartments_from_discretization) {
using point3d = cable_cell::point_type;
cable_cell newcell;
- if (!morph) {
+ if (!morph.num_branches()) {
return newcell;
}
- arb_assert(morph.check_valid());
-
- // (not supporting soma-less cells yet)
- newcell.add_soma(morph.soma.r, point3d(morph.soma.x, morph.soma.y, morph.soma.z));
-
- for (const section_geometry& section: morph.sections) {
- auto kind = section.kind;
- switch (kind) {
- case section_kind::none: // default to dendrite
- kind = section_kind::dendrite;
- break;
- case section_kind::soma:
- throw cable_cell_error("no support for complex somata");
- break;
- default: ;
+ // Add the soma.
+ auto loc = morph.samples()[0].loc; // location of soma.
+
+ // If there is no spherical root/soma use a zero-radius soma.
+ double srad = morph.spherical_root()? loc.radius: 0.;
+ newcell.add_soma(srad, point3d(loc.x, loc.y, loc.z));
+
+ auto& samples = morph.samples();
+ for (auto i: util::make_span(1, morph.num_branches())) {
+ auto index = util::make_range(morph.branch_sample_span(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 (const section_point& p: section.points) {
- radii.push_back(p.r);
- points.push_back(point3d(p.x, p.y, p.z));
+ 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));
}
- auto cable = newcell.add_cable(section.parent_id, kind, radii, points);
+ // Find the id of this branch's parent.
+ auto pid = morph.branch_parent(i);
+ // Adjust pid if a zero-radius soma was used.
+ if (!morph.spherical_root()) {
+ pid = pid==mnpos? 0: pid+1;
+ }
+ auto cable = newcell.add_cable(pid, kind, radii, points);
if (compartments_from_discretization) {
cable->as_cable()->set_compartments(radii.size()-1);
}
diff --git a/arbor/include/arbor/cable_cell.hpp b/arbor/include/arbor/cable_cell.hpp
index f790627a60bf5a515fe7c12dd9a5847c22e1775a..09baa5de0507dc6a8e28f51e3f59112892b45ba0 100644
--- a/arbor/include/arbor/cable_cell.hpp
+++ b/arbor/include/arbor/cable_cell.hpp
@@ -9,7 +9,7 @@
#include <arbor/common_types.hpp>
#include <arbor/constants.hpp>
#include <arbor/mechcat.hpp>
-#include <arbor/morphology.hpp>
+#include <arbor/morph/morphology.hpp>
#include <arbor/segment.hpp>
namespace arb {
@@ -269,10 +269,10 @@ cable_segment* cable_cell::add_cable(cable_cell::index_type parent, Args&&... ar
return segments_.back()->as_cable();
}
-// Create a 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 morphologu.
-cable_cell make_cable_cell(const morphology&, bool compartments_from_discretization=false);
+// 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, bool compartments_from_discretization);
} // namespace arb
diff --git a/arbor/include/arbor/morph/error.hpp b/arbor/include/arbor/morph/error.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..14c40be496cc6be0b8873072a3e18433876c3a4d
--- /dev/null
+++ b/arbor/include/arbor/morph/error.hpp
@@ -0,0 +1,15 @@
+#pragma once
+
+#include <stdexcept>
+
+#include <arbor/arbexcept.hpp>
+
+namespace arb {
+
+struct morphology_error: public arbor_exception {
+ morphology_error(const char* what): arbor_exception(what) {}
+ morphology_error(const std::string& what): arbor_exception(what) {}
+};
+
+} // namespace arb
+
diff --git a/arbor/include/arbor/morph/morphology.hpp b/arbor/include/arbor/morph/morphology.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4ad07955f68af5ea976bdb8d89ecdf74fbffa280
--- /dev/null
+++ b/arbor/include/arbor/morph/morphology.hpp
@@ -0,0 +1,73 @@
+#pragma once
+
+#include <ostream>
+#include <vector>
+
+#include <arbor/util/lexcmp_def.hpp>
+#include <arbor/morph/primitives.hpp>
+#include <arbor/morph/sample_tree.hpp>
+
+namespace arb {
+
+// An unbranched cable segment that has root, terminal or fork point at each end.
+struct mbranch {
+ std::vector<msize_t> index; // sample index
+ msize_t parent_id = mnpos; // branch index
+
+ mbranch() = default;
+ mbranch(std::vector<msize_t> idx, msize_t parent):
+ index(std::move(idx)), parent_id(parent) {}
+
+ bool is_sphere() const { return size()==1u; }
+ msize_t size() const { return index.size(); }
+ bool has_parent() const { return parent_id!=mnpos;}
+
+ friend bool operator==(const mbranch& l, const mbranch& r);
+ friend std::ostream& operator<<(std::ostream& o, const mbranch& b);
+};
+
+class morphology {
+ // The sample tree of sample points and their parent-child relationships.
+ sample_tree samples_;
+
+ // Indicates whether the soma is a sphere.
+ bool spherical_root_ = false;
+
+ // Branch state.
+ std::vector<mbranch> branches_;
+ std::vector<msize_t> branch_parents_;
+ std::vector<std::vector<msize_t>> branch_children_;
+
+ using index_range = std::pair<const msize_t*, const msize_t*>;
+
+ void init();
+
+public:
+ morphology(sample_tree m, bool use_spherical_root);
+ morphology(sample_tree m);
+
+ // Whether the root of the morphology is spherical.
+ bool spherical_root() const;
+
+ // The number of branches in the morphology.
+ msize_t num_branches() const;
+
+ // The parent sample of sample i.
+ const std::vector<msize_t>& sample_parents() const;
+
+ // The parent branch of branch b.
+ msize_t branch_parent(msize_t b) const;
+
+ // The child branches of branch b.
+ const std::vector<msize_t>& branch_children(msize_t b) const;
+
+ // Range of indexes into the sample points in branch b.
+ index_range branch_sample_span(msize_t b) const;
+
+ // Range of the samples in branch b.
+ const std::vector<msample>& samples() const;
+
+ friend std::ostream& operator<<(std::ostream&, const morphology&);
+};
+
+} // namespace arb
diff --git a/arbor/include/arbor/morph/primitives.hpp b/arbor/include/arbor/morph/primitives.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..42baa9b7f46c6d24d28df259fd9478cfb8daf60b
--- /dev/null
+++ b/arbor/include/arbor/morph/primitives.hpp
@@ -0,0 +1,76 @@
+#pragma once
+
+#include <algorithm>
+#include <cstdlib>
+#include <ostream>
+#include <vector>
+
+#include <arbor/util/lexcmp_def.hpp>
+
+//
+// Types used to identify concrete locations.
+//
+
+namespace arb {
+
+using msize_t = std::uint32_t;
+constexpr msize_t mnpos = msize_t(-1);
+
+// a morphology sample point: a 3D location and radius.
+struct mpoint {
+ double x, y, z; // [µm]
+ double radius; // [μm]
+
+ friend std::ostream& operator<<(std::ostream&, const mpoint&);
+};
+
+mpoint lerp(const mpoint& a, const mpoint& b, double u);
+bool is_collocated(const mpoint& a, const mpoint& b);
+double distance(const mpoint& a, const mpoint& b);
+
+// A morphology sample consists of a location and an integer tag.
+// When loaded from an SWC file, the tag will correspond to the SWC label,
+// which are standardised as follows:
+// 1 - soma
+// 2 - axon
+// 3 - (basal) dendrite
+// 4 - apical dendrite
+// http://www.neuronland.org/NLMorphologyConverter/MorphologyFormats/SWC/Spec.html
+// However, any positive integer tag can be provided and labelled dynamically.
+
+struct msample {
+ mpoint loc;
+ int tag;
+
+ friend std::ostream& operator<<(std::ostream&, const msample&);
+};
+
+bool is_collocated(const msample& a, const msample& b);
+double distance(const msample& a, const msample& b);
+
+using point_prop = std::uint8_t;
+enum point_prop_mask: point_prop {
+ point_prop_mask_none = 0,
+ point_prop_mask_root = 1,
+ point_prop_mask_fork = 2,
+ point_prop_mask_terminal = 4,
+ point_prop_mask_collocated = 8
+};
+
+#define ARB_PROP(prop) \
+constexpr bool is_##prop(point_prop p) {\
+ return p&point_prop_mask_##prop;\
+} \
+inline void set_##prop(point_prop& p) {\
+ p |= point_prop_mask_##prop;\
+} \
+inline void unset_##prop(point_prop& p) {\
+ p &= ~point_prop_mask_##prop;\
+}
+
+ARB_PROP(root)
+ARB_PROP(fork)
+ARB_PROP(terminal)
+ARB_PROP(collocated)
+
+} // namespace arb
diff --git a/arbor/include/arbor/morph/sample_tree.hpp b/arbor/include/arbor/morph/sample_tree.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2ae690ebdc13cc9e95e440c763515de9e696c3a1
--- /dev/null
+++ b/arbor/include/arbor/morph/sample_tree.hpp
@@ -0,0 +1,56 @@
+#pragma once
+
+#include <cassert>
+#include <functional>
+#include <vector>
+#include <string>
+
+#include <arbor/swcio.hpp>
+#include <arbor/morph/primitives.hpp>
+
+namespace arb {
+
+/// Morphology composed of samples.
+class sample_tree {
+ std::vector<msample> samples_;
+ std::vector<msize_t> parents_;
+ std::vector<point_prop> props_;
+
+public:
+ sample_tree() = default;
+ sample_tree(std::vector<msample>, std::vector<msize_t>);
+
+ // Reserve space for n samples.
+ void reserve(msize_t n);
+
+ // The append functions return a handle to the last sample appended by the call.
+
+ // Append a single sample.
+ msize_t append(msize_t p, const msample& s); // to sample p.
+ msize_t append(const msample& s); // to the last sample in the tree.
+
+ // Append a sequence of samples.
+ msize_t append(msize_t p, const std::vector<msample>& slist); // to sample p.
+ msize_t append(const std::vector<msample>& slist); // to the last sample in the tree.
+
+ // The number of samples in the tree.
+ msize_t size() const;
+ bool empty() const;
+
+ // The samples in the tree.
+ const std::vector<msample>& samples() const;
+
+ // The parent index of the samples.
+ const std::vector<msize_t>& parents() const;
+
+ // The properties of the samples.
+ const std::vector<point_prop>& properties() const;
+
+ friend std::ostream& operator<<(std::ostream&, const sample_tree&);
+};
+
+/// Build a sample tree from a sequence of swc records.
+sample_tree swc_as_sample_tree(const std::vector<swc_record>& swc_records);
+
+} // namesapce arb
+
diff --git a/arbor/include/arbor/morphology.hpp b/arbor/include/arbor/morphology.hpp
deleted file mode 100644
index 44ac1212d1f7313b5f68bcffa98cf74ff087e5d1..0000000000000000000000000000000000000000
--- a/arbor/include/arbor/morphology.hpp
+++ /dev/null
@@ -1,83 +0,0 @@
-#pragma once
-
-// Representation of 3-d embedded cable morphology, independent of other
-// cell information.
-
-#include <stdexcept>
-#include <vector>
-
-namespace arb {
-
-struct section_point {
- double x, y, z, r; // [µm], r is radius.
-};
-
-enum class section_kind {
- soma,
- dendrite,
- axon,
- none
-};
-
-struct section_geometry {
- unsigned id = 0; // ids should be contigously numbered from 1 in the morphology.
- unsigned parent_id = 0;
- bool terminal = false;
- std::vector<section_point> points;
- double length = 0; // µm
- section_kind kind = section_kind::none;
-
- section_geometry() = default;
- section_geometry(unsigned id, unsigned parent_id, bool terminal, std::vector<section_point> points, double length, section_kind kind = section_kind::none):
- id(id), parent_id(parent_id), terminal(terminal), points(std::move(points)), length(length), kind(kind)
- {}
-
- // Re-discretize the section into ceil(length/dx) segments.
- void segment(double dx);
-};
-
-struct morphology_error: public std::runtime_error {
- morphology_error(const char* what): std::runtime_error(what) {}
- morphology_error(const std::string& what): std::runtime_error(what) {}
-};
-
-struct morphology {
- // origin + spherical radius; convention: r==0 => no soma
- section_point soma = { 0, 0, 0, 0};
-
- std::vector<section_geometry> sections;
-
- bool has_soma() const {
- return soma.r!=0;
- }
-
- bool empty() const {
- return sections.empty() && !has_soma();
- }
-
- operator bool() const { return !empty(); }
-
- // Return number of sections plus soma
- std::size_t components() const {
- return has_soma()+sections.size();
- }
-
- // Check invariants:
- // 1. sections[i].id = i+1 (id 0 corresponds to soma)
- // 2. sections[i].parent_id < sections[i].id
- // 3. sections[i].terminal iff !exists j s.t. sections[j].parent_id = sections[i].id
- bool check_valid() const;
-
- // Throw morphology_error if invariants violated.
- void assert_valid() const;
-
- // Re-discretize all sections.
- void segment(double dx) {
- for (auto& s: sections) s.segment(dx);
- }
-
- // Add new section from sequence of section points. Return reference to new section.
- section_geometry& add_section(std::vector<section_point> points, unsigned parent_id = 0, section_kind kind = section_kind::none);
-};
-
-} // namespace arb
diff --git a/arbor/include/arbor/segment.hpp b/arbor/include/arbor/segment.hpp
index ece282fbe0f122e2cf6c080201709950f0a9be55..10264eb525079b5ba51f78f3f19feee77bbb669d 100644
--- a/arbor/include/arbor/segment.hpp
+++ b/arbor/include/arbor/segment.hpp
@@ -13,13 +13,19 @@
#include <arbor/cable_cell_param.hpp>
#include <arbor/common_types.hpp>
#include <arbor/math.hpp>
-#include <arbor/morphology.hpp>
#include <arbor/mechinfo.hpp>
#include <arbor/point.hpp>
#include <arbor/util/optional.hpp>
namespace arb {
+enum class section_kind {
+ soma,
+ dendrite,
+ axon,
+ none
+};
+
// forward declarations of segment specializations
class soma_segment;
class cable_segment;
diff --git a/arbor/include/arbor/swcio.hpp b/arbor/include/arbor/swcio.hpp
index 50003bd1180db7023ba1a50000a67b5e0107907b..368a44653e85229fb39bd25308657fc745435365 100644
--- a/arbor/include/arbor/swcio.hpp
+++ b/arbor/include/arbor/swcio.hpp
@@ -10,7 +10,6 @@
#include <arbor/assert.hpp>
#include <arbor/arbexcept.hpp>
-#include <arbor/morphology.hpp>
#include <arbor/point.hpp>
namespace arb {
@@ -24,34 +23,25 @@ struct swc_error: public arbor_exception {
class swc_record {
public:
+ static constexpr int undefined_tag = 0;
+
+ using tag_type = int;
using id_type = int;
using coord_type = double;
- // More on SWC files: http://research.mssm.edu/cnic/swc.html
- enum class kind {
- undefined = 0,
- soma,
- axon,
- dendrite,
- apical_dendrite,
- fork_point,
- end_point,
- custom
- };
-
- kind type = kind::undefined; // record type
+ tag_type tag = undefined_tag; // record type
id_type id = 0; // record id
coord_type x = 0; // record coordinates
coord_type y = 0;
coord_type z = 0;
coord_type r = 0; // record radius
- id_type parent_id= -1; // record parent's id
+ id_type parent_id= -1; // record parent's id
// swc records assume zero-based indexing; root's parent remains -1
- swc_record(swc_record::kind type, int id,
+ swc_record(int tag, int id,
coord_type x, coord_type y, coord_type z, coord_type r,
int parent_id):
- type(type), id(id), x(x), y(y), z(z), r(r), parent_id(parent_id)
+ tag(tag), id(id), x(x), y(y), z(z), r(r), parent_id(parent_id)
{}
swc_record() = default;
@@ -81,10 +71,6 @@ public:
return arb::point<coord_type>(x, y, z);
}
- arb::section_point as_section_point() const {
- return arb::section_point{x, y, z, r};
- }
-
// validity checks
bool is_consistent() const;
void assert_consistent() const; // throw swc_error if inconsistent.
@@ -98,9 +84,6 @@ std::istream& operator>>(std::istream& is, swc_record& record);
// Throw on parsing failure.
std::vector<swc_record> parse_swc_file(std::istream& is);
-// Convert a canonical (see below) vector of SWC records to a morphology object.
-morphology swc_as_morphology(const std::vector<swc_record>& swc_records);
-
// Given a vector of random-access mutable sequence of `swc_record` describing
// a single morphology, check for consistency and renumber records
// so that ids are contiguous within branches, have no gaps, and
diff --git a/arbor/morph/morphology.cpp b/arbor/morph/morphology.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..23d38c1d2cb7448d863828f0ba83ee6124b42c33
--- /dev/null
+++ b/arbor/morph/morphology.cpp
@@ -0,0 +1,189 @@
+#include <cmath>
+#include <iostream>
+#include <unordered_map>
+#include <utility>
+
+#include <arbor/math.hpp>
+#include <arbor/morph/error.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/sample_tree.hpp>
+#include <arbor/morph/primitives.hpp>
+
+#include "algorithms.hpp"
+#include "io/sepval.hpp"
+#include "util/span.hpp"
+#include "util/strprintf.hpp"
+
+namespace arb {
+
+namespace impl{
+
+std::vector<mbranch> branches_from_parent_index(const std::vector<msize_t>& parents,
+ const std::vector<point_prop>& props,
+ bool spherical_root)
+{
+ using util::make_span;
+
+ const char* errstr_single_sample_root =
+ "A morphology with only one sample must have a spherical root";
+ const char* errstr_incomplete_cable =
+ "A branch must contain at least two samples";
+
+ if (parents.empty()) return {};
+
+ auto nsamp = parents.size();
+
+ // Enforce that a morphology with one sample has a spherical root.
+ if (!spherical_root && nsamp==1u) {
+ throw morphology_error(errstr_single_sample_root);
+ }
+
+ std::vector<int> bids(nsamp);
+ int nbranches = spherical_root? 1: 0;
+ for (auto i: make_span(1, nsamp)) {
+ auto p = parents[i];
+ bool first = is_root(props[p]) || is_fork(props[p]);
+ bids[i] = first? nbranches++: bids[p];
+ }
+
+ std::vector<mbranch> branches(nbranches);
+ for (auto i: make_span(nsamp)) {
+ auto p = parents[i];
+ auto& branch = branches[bids[i]];
+ if (!branch.size()) {
+ bool null_root = spherical_root? p==mnpos: p==mnpos||p==0;
+ branch.parent_id = null_root? mnpos: bids[p];
+
+ // Add the distal sample from the parent branch if the parent is not
+ // a spherical root or mnpos.
+ if (p!=mnpos && !(p==0 && spherical_root)) {
+ branch.index.push_back(p);
+ }
+ }
+
+ branch.index.push_back(i);
+ }
+
+ // Enforce that all cable branches that are potentially connected to a spherical
+ // root contain at least two samples.
+ if (spherical_root) {
+ for (auto i: make_span(1, nbranches)) { // skip the root.
+ if (branches[i].size()<2u) {
+ throw morphology_error(errstr_incomplete_cable);
+ }
+ }
+ }
+
+ return branches;
+}
+
+// Returns false if one of the root's children has the same tag as the root.
+bool root_sample_has_same_tag_as_child(const sample_tree& st) {
+ if (st.empty()) return false;
+ auto& P = st.parents();
+ auto& S = st.samples();
+ auto root_tag = S.front().tag;
+ for (auto i: util::make_span(1, st.size())) {
+ if (0u==P[i] && S[i].tag==root_tag) {
+ return false;
+ }
+ }
+ return true;
+}
+
+} // namespace impl
+
+bool operator==(const mbranch& l, const mbranch& r) {
+ return l.parent_id==r.parent_id && l.index==r.index;
+}
+
+std::ostream& operator<<(std::ostream& o, const mbranch& b) {
+ o <<"mbranch([" << io::csv(b.index) << "], ";
+ if (b.parent_id==mnpos) o << "none)";
+ else o << b.parent_id << ")";
+ return o;
+}
+
+morphology::morphology(sample_tree m, bool use_spherical_root):
+ samples_(std::move(m)),
+ spherical_root_(use_spherical_root)
+{
+ init();
+}
+
+morphology::morphology(sample_tree m):
+ samples_(std::move(m)),
+ spherical_root_(impl::root_sample_has_same_tag_as_child(samples_))
+{
+ init();
+}
+
+void morphology::init() {
+ using util::make_span;
+ using util::count_along;
+
+ auto nsamp = samples_.size();
+
+ if (!nsamp) return;
+
+ // Generate branches.
+ branches_ = impl::branches_from_parent_index(samples_.parents(), samples_.properties(), spherical_root_);
+ auto nbranch = branches_.size();
+
+ // Generate branch tree.
+ branch_children_.resize(nbranch);
+ branch_parents_.reserve(nbranch);
+ for (auto i: make_span(nbranch)) {
+ auto id = branches_[i].parent_id;
+ branch_parents_.push_back(id);
+ if (id!=mnpos) {
+ branch_children_[id].push_back(i);
+ }
+ }
+}
+
+// The parent branch of branch b.
+msize_t morphology::branch_parent(msize_t b) const {
+ return branch_parents_[b];
+}
+
+// The parent sample of sample i.
+const std::vector<msize_t>& morphology::sample_parents() const {
+ return samples_.parents();
+}
+
+// The child branches of branch b.
+const std::vector<msize_t>& morphology::branch_children(msize_t b) const {
+ return branch_children_[b];
+}
+
+// Whether the root of the morphology is spherical.
+bool morphology::spherical_root() const {
+ return spherical_root_;
+}
+
+morphology::index_range morphology::branch_sample_span(msize_t b) const {
+ const auto& idx = branches_[b].index;
+ return std::make_pair(idx.data(), idx.data()+idx.size());
+}
+
+const std::vector<msample>& morphology::samples() const {
+ return samples_.samples();
+}
+
+msize_t morphology::num_branches() const {
+ return branches_.size();
+}
+
+std::ostream& operator<<(std::ostream& o, const morphology& m) {
+ o << "morphology: "
+ << m.samples_.size() << " samples, "
+ << m.num_branches() << " branches.";
+ for (auto i: util::make_span(m.num_branches()))
+ o << "\n branch " << i << ": " << m.branches_[i];
+
+ return o;
+}
+
+} // namespace arb
+
diff --git a/arbor/morph/primitives.cpp b/arbor/morph/primitives.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6decddbe1c93fa5c2ee35da12d59a9001fff77e8
--- /dev/null
+++ b/arbor/morph/primitives.cpp
@@ -0,0 +1,50 @@
+#include <ostream>
+
+#include <arbor/math.hpp>
+#include <arbor/morph/primitives.hpp>
+
+#include "io/sepval.hpp"
+#include "util/span.hpp"
+#include "util/rangeutil.hpp"
+
+namespace arb {
+
+// interpolate between two points.
+mpoint lerp(const mpoint& a, const mpoint& b, double u) {
+ return { math::lerp(a.x, b.x, u),
+ math::lerp(a.y, b.y, u),
+ math::lerp(a.z, b.z, u),
+ math::lerp(a.radius, b.radius, u) };
+}
+
+// test if two morphology sample points share the same location.
+bool is_collocated(const mpoint& a, const mpoint& b) {
+ return a.x==b.x && a.y==b.y && a.z==b.z;
+}
+
+// calculate the distance between two morphology sample points.
+double distance(const mpoint& a, const mpoint& b) {
+ double dx = a.x - b.x;
+ double dy = a.y - b.y;
+ double dz = a.z - b.z;
+
+ return std::sqrt(dx*dx + dy*dy * dz*dz);
+}
+
+bool is_collocated(const msample& a, const msample& b) {
+ return is_collocated(a.loc, b.loc);
+}
+
+double distance(const msample& a, const msample& b) {
+ return distance(a.loc, b.loc);
+}
+
+std::ostream& operator<<(std::ostream& o, const mpoint& p) {
+ return o << "mpoint(" << p.x << "," << p.y << "," << p.z << "," << p.radius << ")";
+}
+
+std::ostream& operator<<(std::ostream& o, const msample& s) {
+ return o << "msample(" << s.loc << ", " << s.tag << ")";
+}
+
+} // namespace arb
diff --git a/arbor/morph/sample_tree.cpp b/arbor/morph/sample_tree.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7a687e52ca9eb966470fbdd23f2c9d14ff010c93
--- /dev/null
+++ b/arbor/morph/sample_tree.cpp
@@ -0,0 +1,132 @@
+#include <cmath>
+#include <vector>
+
+#include <arbor/math.hpp>
+
+#include <arbor/morph/error.hpp>
+#include <arbor/morph/sample_tree.hpp>
+
+#include "algorithms.hpp"
+#include "io/sepval.hpp"
+#include "util/span.hpp"
+
+namespace arb {
+
+sample_tree::sample_tree(std::vector<msample> samples, std::vector<msize_t> parents) {
+ if (samples.size()!=parents.size()) {
+ throw std::runtime_error(
+ "The same number of samples and parent indices used to create a sample morphology");
+ }
+ reserve(samples.size());
+ for (auto i: util::make_span(samples.size())) {
+ append(parents[i], samples[i]);
+ }
+}
+
+void sample_tree::reserve(msize_t n) {
+ samples_.reserve(n);
+ parents_.reserve(n);
+ props_.reserve(n);
+}
+
+msize_t sample_tree::append(msize_t p, const msample& s) {
+ if (empty()) {
+ if (p!=mnpos) throw morphology_error("Parent id of root sample must be mnpos");
+ }
+ else if (p>=size()) {
+ throw morphology_error("Parent id of a sample must be less than the sample id");
+ }
+ auto id = size();
+
+ samples_.push_back(s);
+ parents_.push_back(p);
+
+ // Set the point properties for the new point, and update those of its parent as needed.
+ point_prop prop = point_prop_mask_none;
+ if (!id) {
+ // if adding the first sample mark it as root
+ set_root(prop);
+ }
+ else {
+ // Mark the new node as terminal, and unset the parent sample's terminal bit.
+ set_terminal(prop);
+ const bool term_parent = is_terminal(props_[p]); // track if the parent was terminal.
+ unset_terminal(props_[p]);
+
+ // Mark if the new sample is collocated with its parent.
+ if (is_collocated(s, samples_[p])) {
+ set_collocated(prop);
+ }
+
+ // Set parent to be a fork if it was not a terminal point before the
+ // new sample was added (and if it isn't the root).
+ if (p && !term_parent) {
+ set_fork(props_[p]);
+ }
+ }
+ props_.push_back(prop);
+
+ return id;
+}
+
+msize_t sample_tree::append(const msample& s) {
+ std::cout << "decision : " << empty() << " " << (empty()? mnpos: size()-1) << "\n";
+ return append(empty()? mnpos: size()-1, s);
+}
+
+msize_t sample_tree::append(msize_t p, const std::vector<msample>& slist) {
+ if (!slist.size()) return size();
+
+ for (auto& s: slist) {
+ p = append(p, s);
+ }
+
+ return p;
+}
+
+msize_t sample_tree::append(const std::vector<msample>& slist) {
+ return append(empty()? mnpos: size()-1, slist);
+}
+
+msize_t sample_tree::size() const {
+ return samples_.size();
+}
+
+bool sample_tree::empty() const {
+ return samples_.empty();
+}
+
+const std::vector<msample>& sample_tree::samples() const {
+ return samples_;
+}
+
+const std::vector<msize_t>& sample_tree::parents() const {
+ return parents_;
+}
+
+const std::vector<point_prop>& sample_tree::properties() const {
+ return props_;
+}
+
+std::ostream& operator<<(std::ostream& o, const sample_tree& m) {
+ o << "sample_tree:"
+ << "\n " << m.size() << " samples"
+ << "\n samples [" << io::csv(m.samples_) << "]"
+ << "\n parents [" << io::csv(m.parents_) << "]";
+ return o;
+}
+
+sample_tree swc_as_sample_tree(const std::vector<swc_record>& swc_records) {
+ sample_tree m;
+ m.reserve(swc_records.size());
+
+ for (auto i: util::count_along(swc_records)) {
+ auto& r = swc_records[i];
+ // The parent of soma must be mnpos, while in SWC files is -1
+ msize_t p = i==0? mnpos: r.parent_id;
+ m.append(p, msample{{r.x, r.y, r.z, r.r}, r.tag});
+ }
+ return m;
+}
+
+} // namespace arb
diff --git a/arbor/morphology.cpp b/arbor/morphology.cpp
deleted file mode 100644
index 2243772b9cfe99f89460080cf9fd8e5b7664a349..0000000000000000000000000000000000000000
--- a/arbor/morphology.cpp
+++ /dev/null
@@ -1,135 +0,0 @@
-#include <cmath>
-#include <vector>
-
-#include <arbor/morphology.hpp>
-#include <arbor/math.hpp>
-
-namespace arb {
-
-using ::arb::math::lerp;
-
-static section_point lerp(const section_point& a, const section_point& b, double u) {
- return { lerp(a.x, b.x, u), lerp(a.y, b.y, u), lerp(a.z, b.z, u), lerp(a.r, b.r, u) };
-}
-
-static double distance(const section_point& a, const section_point& b) {
- double dx = b.x-a.x;
- double dy = b.y-a.y;
- double dz = b.z-a.z;
-
- return std::sqrt(dx*dx+dy*dy+dz*dz);
-}
-
-void section_geometry::segment(double dx) {
- unsigned npoint = points.size();
- if (dx<=0 || npoint<2) return;
-
- // Re-discretize into nseg segments (nseg+1 points).
- unsigned nseg = static_cast<unsigned>(std::ceil(length/dx));
-
- std::vector<section_point> sampled;
- sampled.push_back(points.front());
- double sampled_length = 0;
-
- // [left, right) is the path-length interval for successive
- // linear segments in the section.
- double left = 0;
- double right = left+distance(points[1], points[0]);
-
- // x is the next sample point (in path-length).
- double x = length/nseg;
-
- // Scan segments for sample points.
- unsigned i = 1;
- for (;;) {
- if (right>x) {
- double u = (x-left)/(right-left);
- sampled.push_back(lerp(points[i-1], points[i], u));
- unsigned k = sampled.size();
- sampled_length += distance(sampled[k-2], sampled[k-1]);
- x = k*length/nseg;
- }
- else {
- ++i;
- if (i>=npoint) break;
-
- left = right;
- right = left+distance(points[i-1], points[i]);
- }
- }
- if (sampled.size()<=nseg) {
- sampled.push_back(points.back());
- }
-
- points = std::move(sampled);
- length = sampled_length;
-}
-
-static double compute_length(const std::vector<section_point>& points) {
- double length = 0;
- std::size_t npoint = points.size();
-
- for (std::size_t i =1; i<npoint; ++i) {
- length += distance(points[i], points[i-1]);
- }
-
- return length;
-}
-
-section_geometry& morphology::add_section(std::vector<section_point> points, unsigned parent_id, section_kind kind) {
- section_geometry section;
- section.id = sections.size()+1;
- section.parent_id = parent_id;
- section.terminal = true;
- section.points = std::move(points);
- section.kind = kind;
- section.length = compute_length(section.points);
-
- if (section.parent_id >= section.id) {
- throw morphology_error("improper parent id for section");
- }
-
- if (section.parent_id>0) {
- sections[section.parent_id-1].terminal = false;
- }
- sections.push_back(std::move(section));
- return sections.back();
-}
-
-static const char* morphology_invariant_violation(const morphology& m) {
- std::size_t nsection = m.sections.size();
- std::vector<int> terminal(nsection, true);
-
- for (std::size_t i=0; i<nsection; ++i) {
- auto id = m.sections[i].id;
- auto parent_id = m.sections[i].parent_id;
-
- if (id!=i+1) return "section id does not correspond to index";
- if (parent_id>=id) return "section parent id not less than section id";
- if (parent_id>0) {
- terminal[parent_id-1] = false;
- }
- if (parent_id==0 && !m.has_soma()) return "section has parent 0 but morphology has no (spherical) soma";
- }
-
- for (std::size_t i=0; i<nsection; ++i) {
- if (terminal[i] && !m.sections[i].terminal) return "non-terminal section is marked terminal";
- if (!terminal[i] && m.sections[i].terminal) return "terminal section is marked non-terminal";
- }
-
- return nullptr;
-}
-
-bool morphology::check_valid() const {
- return morphology_invariant_violation(*this)==nullptr;
-}
-
-void morphology::assert_valid() const {
- auto error = morphology_invariant_violation(*this);
- if (error) throw morphology_error(error);
-}
-
-
-
-
-} // namespace arb
diff --git a/arbor/swcio.cpp b/arbor/swcio.cpp
index ef18c38619eb3466ecff5d83302b31d1cf4449b7..8c0b14b5d4a075c4b318ec88b753069e7faeff9c 100644
--- a/arbor/swcio.cpp
+++ b/arbor/swcio.cpp
@@ -6,7 +6,6 @@
#include <unordered_set>
#include <arbor/assert.hpp>
-#include <arbor/morphology.hpp>
#include <arbor/point.hpp>
#include <arbor/swcio.hpp>
@@ -19,10 +18,8 @@ namespace arb {
// helper function: return error message if inconsistent, or nullptr if ok.
const char* swc_record_error(const swc_record& r) {
- constexpr int max_type = static_cast<int>(swc_record::kind::custom);
-
- if (static_cast<int>(r.type) < 0 || static_cast<int>(r.type) > max_type) {
- return "unknown record type";
+ if (r.tag<0) {
+ return "unknown record tag";
}
if (r.id < 0) {
@@ -60,9 +57,7 @@ static bool parse_record(const std::string& line, swc_record& record) {
std::istringstream is(line);
swc_record r;
- int type_as_int;
- is >> r.id >> type_as_int >> r.x >> r.y >> r.z >> r.r >> r.parent_id;
- r.type = static_cast<swc_record::kind>(type_as_int);
+ is >> r.id >> r.tag >> r.x >> r.y >> r.z >> r.r >> r.parent_id;
if (is) {
// Convert to zero-based, leaving parent_id as-is if -1
--r.id;
@@ -105,7 +100,7 @@ std::istream& operator>>(std::istream& is, swc_record& record) {
std::ostream& operator<<(std::ostream& os, const swc_record& record) {
// output in one-based indexing
os << record.id+1 << " "
- << static_cast<int>(record.type) << " "
+ << record.tag << " "
<< std::setprecision(7) << record.x << " "
<< std::setprecision(7) << record.y << " "
<< std::setprecision(7) << record.z << " "
@@ -154,70 +149,6 @@ std::vector<swc_record> parse_swc_file(std::istream& is) {
return records;
}
-morphology swc_as_morphology(const std::vector<swc_record>& swc_records) {
- morphology morph;
-
- std::vector<swc_record::id_type> swc_parent_index;
- for (const auto& r: swc_records) {
- swc_parent_index.push_back(r.parent_id);
- }
-
- if (swc_parent_index.empty()) {
- return morph;
- }
-
- // The parent of soma must be 0, while in SWC files is -1
- swc_parent_index[0] = 0;
- auto branch_index = algorithms::branches(swc_parent_index); // partitions [0, #records] by branch.
- auto parent_branch_index = algorithms::tree_reduce(swc_parent_index, branch_index);
-
- // sanity check
- arb_assert(parent_branch_index.size() == branch_index.size() - 1);
-
- // Add the soma first; then the segments
- const auto& soma = swc_records[0];
- morph.soma = { soma.x, soma.y, soma.z, soma.r };
-
- for (auto i: util::make_span(1, parent_branch_index.size())) {
- auto b_start = swc_records.begin() + branch_index[i];
- auto b_end = swc_records.begin() + branch_index[i+1];
-
- unsigned parent_id = parent_branch_index[i];
- std::vector<section_point> points;
- section_kind kind = section_kind::none;
-
- if (parent_id != 0) {
- // include the parent of current record if not branching from soma
- auto parent_record = swc_records[swc_parent_index[branch_index[i]]];
-
- points.push_back(section_point{parent_record.x, parent_record.y, parent_record.z, parent_record.r});
- }
-
- for (auto b = b_start; b!=b_end; ++b) {
- points.push_back(section_point{b->x, b->y, b->z, b->r});
-
- switch (b->type) {
- case swc_record::kind::axon:
- kind = section_kind::axon;
- break;
- case swc_record::kind::dendrite:
- case swc_record::kind::apical_dendrite:
- kind = section_kind::dendrite;
- break;
- case swc_record::kind::soma:
- kind = section_kind::soma;
- break;
- default: ; // stick with what we have
- }
- }
-
- morph.add_section(std::move(points), parent_id, kind);
- }
-
- morph.assert_valid();
- return morph;
-}
-
void swc_canonicalize(std::vector<swc_record>& swc_records) {
std::unordered_set<swc_record::id_type> ids;
diff --git a/arbor/util/range.hpp b/arbor/util/range.hpp
index 154f8c35969b3a614c3987ace87717cc00ed09ef..9f175abf98e06e971e8280e28755bf6d8ffe41a7 100644
--- a/arbor/util/range.hpp
+++ b/arbor/util/range.hpp
@@ -28,11 +28,11 @@
#include <arbor/assert.hpp>
-#include <util/counter.hpp>
-#include <util/either.hpp>
-#include <util/iterutil.hpp>
-#include <util/meta.hpp>
-#include <util/sentinel.hpp>
+#include "util/either.hpp"
+#include "util/counter.hpp"
+#include "util/iterutil.hpp"
+#include "util/meta.hpp"
+#include "util/sentinel.hpp"
namespace arb {
namespace util {
diff --git a/example/single/single.cpp b/example/single/single.cpp
index 3855f9de4bbed6b618914773abdf67fd178ad232..df614950c3cf6b2060e48593dc71c92a283c8aca 100644
--- a/example/single/single.cpp
+++ b/example/single/single.cpp
@@ -7,7 +7,8 @@
#include <arbor/load_balance.hpp>
#include <arbor/cable_cell.hpp>
-#include <arbor/morphology.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/sample_tree.hpp>
#include <arbor/swcio.hpp>
#include <arbor/simulation.hpp>
#include <arbor/simple_sampler.hpp>
@@ -26,7 +27,7 @@ arb::morphology default_morphology();
arb::morphology read_swc(const std::string& path);
struct single_recipe: public arb::recipe {
- explicit single_recipe(arb::morphology m): morpho(m) {
+ explicit single_recipe(arb::morphology m): morpho(std::move(m)) {
gprop.default_parameters = arb::neuron_parameter_defaults;
}
@@ -53,7 +54,7 @@ struct single_recipe: public arb::recipe {
}
arb::util::unique_any get_cell_description(arb::cell_gid_type) const override {
- arb::cable_cell c = make_cable_cell(morpho);
+ arb::cable_cell c = make_cable_cell(morpho, false);
// Add HH mechanism to soma, passive channels to dendrites.
// Discretize dendrites according to the NEURON d-lambda rule.
@@ -70,9 +71,9 @@ struct single_recipe: public arb::recipe {
seg->set_compartments(n);
}
- // Add synapse to last segment.
+ // Add synapse to last branch.
- arb::cell_lid_type last_segment = morpho.components()-1;
+ arb::cell_lid_type last_segment = morpho.num_branches()-1;
arb::segment_location end_last_segment = { last_segment, 1. };
c.add_synapse(end_last_segment, "exp2syn");
@@ -147,25 +148,19 @@ options parse_options(int argc, char** argv) {
// to 0.2 µm.
arb::morphology default_morphology() {
- arb::morphology m;
+ arb::sample_tree samples;
- // Points in a morphology are expressed as (x, y, z, r),
- // with r being the radius. Units are µm.
+ auto p = samples.append(arb::msample{{ 0.0, 0.0, 0.0, 6.3}, 1});
+ std::cout << "p is " << p << "\n";
+ p = samples.append(p, arb::msample{{ 6.3, 0.0, 0.0, 0.5}, 3});
+ p = samples.append(p, arb::msample{{206.3, 0.0, 0.0, 0.2}, 3});
- m.soma = {0., 0., 0., 6.3};
-
- std::vector<arb::section_point> dendrite = {
- {6.3, 0., 0., 0.5},
- {206.3, 0., 0., 0.2}
- };
- m.add_section(dendrite, 0);
-
- return m;
+ return arb::morphology(std::move(samples));
}
arb::morphology read_swc(const std::string& path) {
std::ifstream f(path);
if (!f) throw std::runtime_error("unable to open SWC file: "+path);
- return arb::swc_as_morphology(arb::parse_swc_file(f));
+ return arb::morphology(arb::swc_as_sample_tree(arb::parse_swc_file(f)));
}
diff --git a/lmorpho/CMakeLists.txt b/lmorpho/CMakeLists.txt
deleted file mode 100644
index c198d458fd6135fe56c1e6fd2ed0e9a1e2572aa6..0000000000000000000000000000000000000000
--- a/lmorpho/CMakeLists.txt
+++ /dev/null
@@ -1,8 +0,0 @@
-add_executable(lmorpho lmorpho.cpp lsystem.cpp lsys_models.cpp morphio.cpp)
-
-target_link_libraries(lmorpho PRIVATE arbor arbor-sup)
-
-# TODO: resolve public headers
-target_link_libraries(lmorpho PRIVATE arbor-private-headers)
-
-install(TARGETS lmorpho RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
diff --git a/lmorpho/README.md b/lmorpho/README.md
deleted file mode 100644
index 561bbebd863f218a44fddb4affb65203ca9d0775..0000000000000000000000000000000000000000
--- a/lmorpho/README.md
+++ /dev/null
@@ -1,115 +0,0 @@
-# `lmporho` — generate random neuron morphologies
-
-`lmorpho` implements an L-system geometry generator for the purpose of
-generating a series of artificial neuron morphologies as described by a set of
-stasticial parameters.
-
-## Method
-
-The basic algorithm used is that of [Burke (1992)][burke1992], with extensions
-for embedding in 3-d taken from [Ascoli (2001)][ascoli2001].
-
-A dendritic tree is represented by a piece-wise linear embedding of a rooted
-tree into R³×Râº, describing the position in space and non-negative radius.
-Morphologies are represented as a sequence of these trees together with a point
-and radius describing a (spherical) soma.
-
-The generation process for a tree successively grows the terminal points by
-some length ΔL each step, changing its orientation each step according to a
-random distribution. After growing, there is a radius-dependent probability
-that it will bifurcate or terminate. Parameters describing these various
-probability distributions constitute the L-system model.
-
-The distributions used in [Ascoli (2001)][ascoli2001] can be constant (i.e. not
-random), uniform over some interval, truncated normal, or a mixture
-distribution of these. In the present version of the code, mixture
-distributions are not supported.
-
-## Output
-
-Generated morphologies can be output either in
-[SWC format](http://research.mssm.edu/cnic/swc.html), or as 'parent vectors',
-which describe the topology of the associated tree.
-
-Entry _i_ of the (zero-indexed) parent vector gives the index of the proximal
-unbranched section to which it connects. Somata have a parent index of -1.
-
-Morphology information can be written one each to separate files or
-concatenated. If concatenated, the parent vector indices will be shifted as
-required to maintain consistency.
-
-## Models
-
-Two L-system parameter sets are provided as 'built-in'; if neither model is
-selected, a simple, non-biological default model is used.
-
-As mixture distributions are not yet supported in `lmorpho`, these models are
-simplified approximations to those in the literature as described below.
-
-### Motor neuron model
-
-The 'motoneuron' model corresponds to the adult cat spinal alpha-motoneuron
-model described in [Ascoli (2001)][ascoli2001], based in turn on the model of
-[Burke (1992)][burke1992].
-
-The model parameters in Ascoli are not complete; missing parameters were taken
-from the corresponding ArborVitae model parameters in the same paper, and
-somatic diameters were taken from [Cullheim (1987)][cullheim1987].
-
-### Purkinke neuron model
-
-The 'purkinke' model corresponds to the guinea pig Purkinje cell model from
-[Ascoli (2001)][ascoli2001], which was fit to data from [Rapp
-(1994)][rapp1994]. Somatic diameters are a simple fit to the three measurements
-in the original source.
-
-Some required parameters are not recorded in [Ascoli (2001)][ascoli2001];
-the correlation component `diam_child_a` and the `length_step` ΔL are
-taken from the corresponding L-Neuron parameter description in the
-[L-Neuron database](l-neuron-database). These should be verified by
-fitting against the digitized morphologies as found in the database.
-
-Produced neurons from this model do not look especially realistic.
-
-### Other models
-
-There is not yet support for loading in arbitrary parameter sets, or for
-translating or approximating the parameters that can be found in the
-[L-Neuron database][l-neuron-database]. Support for parameter estimation
-from a population of discretized morphologies would be useful.
-
-## References
-
-1. Ascoli, G. A., Krichmar, J. L., Scorcioni, R., Nasuto, S. J., Senft, S. L.
- and Krichmar, G. L. (2001). Computer generation and quantitative morphometric
- analysis of virtual neurons. _Anatomy and Embryology_ _204_(4), 283–301.
- [DOI: 10.1007/s004290100201][ascoli2001]
-
-2. Burke, R. E., Marks, W. B. and Ulfhake, B. (1992). A parsimonious description
- of motoneuron dendritic morphology using computer simulation.
- _The Journal of Neuroscience_ _12_(6), 2403–2416. [online][burke1992]
-
-3. Cullheim, S., Fleshman, J. W., Glenn, L. L., and Burke, R. E. (1987).
- Membrane area and dendritic structure in type-identified triceps surae alpha
- motoneurons. _The Journal of Comparitive Neurology_ _255_(1), 68–81.
- [DOI: 10.1002/cne.902550106][cullheim1987]
-
-4. Rapp, M., Segev, I. and Yaom, Y. (1994). Physiology, morphology and detailed
- passive models of guinea-pig cerebellar Purkinje cells.
- _The Journal of Physiology_, _474_, 101–118.
- [DOI: 10.1113/jphysiol.1994.sp020006][rapp1994].
-
-
-[ascoli2001]: http://dx.doi.org/10.1007/s004290100201
- "Ascoli et al. (2001). Computer generation […] of virtual neurons."
-
-[burke1992]: http://www.jneurosci.org/content/12/6/2403
- "Burke et al. (1992). A parsimonious description of motoneuron dendritic morphology […]"
-
-[cullheim1987]: http://dx.doi.org/10.1002/cne.902550106
- "Cullheim et al. (1987). Membrane area and dendritic structure in […] alpha motoneurons."
-
-[rapp1994]: http://dx.doi.org/10.1113/jphysiol.1994.sp020006
- "Rapp et al. (1994). Physiology, morphology […] of guinea-pig cerebellar Purkinje cells."
-
-[l-neuron-database]: http://krasnow1.gmu.edu/cn3/L-Neuron/database/index.html
diff --git a/lmorpho/lmorpho.cpp b/lmorpho/lmorpho.cpp
deleted file mode 100644
index 56385a925e9e88cf2daec26646249eb0e2c1a5a8..0000000000000000000000000000000000000000
--- a/lmorpho/lmorpho.cpp
+++ /dev/null
@@ -1,114 +0,0 @@
-#include <algorithm>
-#include <iostream>
-#include <fstream>
-#include <map>
-#include <sstream>
-#include <vector>
-
-#include <arbor/morphology.hpp>
-#include <arbor/util/optional.hpp>
-#include <sup/tinyopt.hpp>
-
-#include "morphio.hpp"
-#include "lsystem.hpp"
-#include "lsys_models.hpp"
-
-using arb::util::optional;
-using arb::util::nullopt;
-using arb::util::just;
-
-const char* usage_str =
-"[OPTION]...\n"
-"\n"
-" -n, --count=N Number of morphologies to generate.\n"
-" -m, --model=MODEL Use L-system MODEL for generation (see below).\n"
-" -g, --segment=DX Segment model into compartments of max size DX µm.\n"
-" --swc=FILE Output morphologies as SWC to FILE (see below).\n"
-" --pvec=FILE Output 'parent vector' structural representation\n"
-" to FILE.\n"
-" -h, --help Emit this message and exit.\n"
-"\n"
-"Generate artificial neuron morphologies based on L-system descriptions.\n"
-"\n"
-"If a FILE argument contains a '%', then one file will be written for\n"
-"each generated morphology, with the '%' replaced by the index of the\n"
-"morphology, starting from zero. Output for each morphology will otherwise\n"
-"be concatenated: SWC files will be headed by a comment line with the\n"
-"index of the morphology; parent vectors will be merged into one long\n"
-"vector. A FILE argument of '-' corresponds to standard output.\n"
-"\n"
-"Currently supported MODELs:\n"
-" motoneuron Adult cat spinal alpha-motoneurons, based on models\n"
-" and data in Burke 1992 and Ascoli 2001.\n"
-" purkinje Guinea pig Purkinje cells, basd on models and data\n"
-" from Rapp 1994 and Ascoli 2001.\n";
-
-int main(int argc, char** argv) {
- // options
- int n_morph = 1;
- optional<unsigned> rng_seed;
- optional<std::string> swc_file;
- optional<std::string> pvector_file;
- double segment_dx = 0;
-
- std::pair<const char*, const lsys_param*> models[] = {
- {"motoneuron", &alpha_motoneuron_lsys},
- {"purkinje", &purkinje_lsys}
- };
- lsys_param P;
-
- try {
- auto arg = argv+1;
- while (*arg) {
- if (auto o = to::parse_opt<int>(arg, 'n', "count")) {
- n_morph = *o;
- }
- if (auto o = to::parse_opt<int>(arg, 's', "seed")) {
- rng_seed = *o;
- }
- else if (auto o = to::parse_opt<std::string>(arg, 0, "swc")) {
- swc_file = *o;
- }
- else if (auto o = to::parse_opt<double>(arg, 'g', "segment")) {
- segment_dx = *o;
- }
- else if (auto o = to::parse_opt<std::string>(arg, 'p', "pvec")) {
- pvector_file = *o;
- }
- else if (auto o = to::parse_opt<const lsys_param*>(arg, 'm', "model", to::keywords(models))) {
- P = **o;
- }
- else if (to::parse_opt(arg, 'h', "help")) {
- std::cout << "Usage: " << argv[0] << " " << usage_str;
- return 0;
- }
- else {
- throw to::parse_opt_error(*arg, "unrecognized option");
- }
- }
-
- std::minstd_rand g;
- if (rng_seed) g.seed(rng_seed.value());
-
- auto emit_swc = swc_file? just(swc_emitter(*swc_file, n_morph)): nullopt;
- auto emit_pvec = pvector_file? just(pvector_emitter(*pvector_file, n_morph)): nullopt;
-
- for (int i=0; i<n_morph; ++i) {
- auto morph = generate_morphology(P, g);
- morph.segment(segment_dx);
-
- if (emit_swc) (*emit_swc)(i, morph);
- if (emit_pvec) (*emit_pvec)(i, morph);
- }
- }
- catch (to::parse_opt_error& e) {
- std::cerr << argv[0] << ": " << e.what() << "\n";
- std::cerr << "Try '" << argv[0] << " --help' for more information.\n";
- std::exit(2);
- }
- catch (std::exception& e) {
- std::cerr << "caught exception: " << e.what() << "\n";
- std::exit(1);
- }
-}
-
diff --git a/lmorpho/lsys_models.cpp b/lmorpho/lsys_models.cpp
deleted file mode 100644
index 41caf7e19bdb70606191f304343018c16aa92464..0000000000000000000000000000000000000000
--- a/lmorpho/lsys_models.cpp
+++ /dev/null
@@ -1,131 +0,0 @@
-#include <cmath>
-
-#include "lsystem.hpp"
-#include "lsys_models.hpp"
-
-static constexpr double inf = INFINITY;
-
-// Predefined parameters for two classes of neurons. Numbers taken primarily
-// from Ascoli et al. 2001, but some details (soma diameters for example)
-// taken from earlier literature.
-//
-// Without mixture distribution support, these aren't entirely faithful to
-// the original sources.
-//
-// Refer to `README.md` for references and more information.
-
-lsys_param make_alpha_motoneuron_lsys() {
- lsys_param L;
-
- // Soma diameter [µm].
- L.diam_soma = { 47.0, 65.5, 57.6, 5.0 };
-
- // Number of dendritic trees (rounded to nearest int).
- L.n_tree = { 8.0, 16.0 };
-
- // Initial dendrite diameter [µm]. (Dstem)
- L.diam_initial = { 3.0, 12.0 };
-
- // Dendrite step length [µm]. (ΔL)
- L.length_step = { 25 };
-
- // Initial roll (intrinsic rotation about x-axis) [degrees].
- L.roll_initial = { -180.0, 180.0 };
-
- // Initial pitch (intrinsic rotation about y-axis) [degrees].
- L.pitch_initial = { 0.0, 180.0 };
-
- // Tortuousness: roll within section over ΔL [degrees].
- L.roll_section = { -180.0, 180.0 };
-
- // Tortuousness: pitch within section over ΔL [degrees].
- L.pitch_section = { -15.0, 15.0, 0.0, 5.0 };
-
- // Taper rate: diameter decrease per unit length.
- L.taper = { -1.25e-3 };
-
- // Branching torsion: roll at branch point [degrees].
- L.roll_at_branch = { 0.0, 180.0 };
-
- // Branch angle between siblings [degrees].
- L.branch_angle = { 1.0, 172.0, 45.0, 20.0 };
-
- // Correlated child branch radius distribution parameters.
- L.diam_child_a = -0.2087;
- L.diam_child_r = { 0.2, 1.8, 0.8255, 0.2125 };
-
- // Termination probability parameters [1/µm].
- L.pterm_k1 = 2.62e-2;
- L.pterm_k2 = -2.955;
-
- // Bifurcation probability parameters [1/μm].
- L.pbranch_ov_k1 = 2.58e-5;
- L.pbranch_ov_k2 = 2.219;
- L.pbranch_nov_k1 = 2.34e-3;
- L.pbranch_nov_k2 = 0.194;
-
- // Absolute maximum dendritic extent [µm].
- L.max_extent = 2000;
-
- return L;
-}
-
-lsys_param alpha_motoneuron_lsys = make_alpha_motoneuron_lsys();
-
-// Some parameters missing from the literature are taken from
-// the `purkburk.prm` L-Neuron parameter file:
-// http://krasnow1.gmu.edu/cn3/L-Neuron/database/purk/bu/purkburk.prm
-// Note 'Pov' and 'Pterm' numbers are swapped in Ascoli fig 3B.
-
-lsys_param make_purkinje_lsys() {
- lsys_param L;
-
- // Soma diameter [µm].
- // (Gaussian fit to 3 samples rom Rapp 1994); truncate to 2σ.
- L.diam_soma = { 22.12, 27.68, 24.9, 1.39 };
-
- // Number of dendritic trees (rounded to nearest int).
- L.n_tree = { 1 };
-
- // Initial dendrite diameter [µm]. (Dstem)
- L.diam_initial = { 4.8, 7.6, 6.167, 1.069 };
-
- // Dendrite step length [µm]. (ΔL)
- L.length_step = { 2.3 }; // from `purkburk.prm`
-
- // Tortuousness: roll within section over ΔL [degrees].
- L.roll_section = { 0 };
-
- // Tortuousness: pitch within section over ΔL [degrees].
- L.pitch_section = { -5, 5 }; // from `purkburk.prm`
-
- // Taper rate: diameter decrease per unit length.
- L.taper = { -inf, inf, -0.010, 0.022 };
-
- // Branching torsion: roll at branch point [degrees].
- L.roll_at_branch = { -inf, inf, 0.0, 6.5 }; // from `purkburj.prm`
-
- // Branch angle between siblings [degrees].
- L.branch_angle = { 1.0, 179.0, 71.0, 40.0 };
-
- // Correlated child branch radius distribution parameters.
- L.diam_child_a = 5.47078; // from `purkburk.prm`
- L.diam_child_r = { 0, inf, 0.112, 0.035 };
-
- // Termination probability parameters [1/µm].
- L.pterm_k1 = 0.1973; // from `purkburj.prm`; 0.4539 in Ascoli.
- L.pterm_k2 = -1.1533;
-
- // Bifurcation probability parameters [1/μm].
- L.pbranch_ov_k1 = 0.1021; // from `purkburk.prm`; 0.0355 in Ascoli.
- L.pbranch_ov_k2 = 0.7237;
- L.pbranch_nov_k1 = 0.1021; // from `purkburk.prm`; 0.2349 in Ascoli.
- L.pbranch_nov_k2 = -0.0116;
-
- // Absolute maximum dendritic extent [µm].
- L.max_extent = 20000;
-
- return L;
-}
-
-lsys_param purkinje_lsys = make_purkinje_lsys();
diff --git a/lmorpho/lsys_models.hpp b/lmorpho/lsys_models.hpp
deleted file mode 100644
index 96638e5ad0322cd830b4bf283a15fee152650592..0000000000000000000000000000000000000000
--- a/lmorpho/lsys_models.hpp
+++ /dev/null
@@ -1,9 +0,0 @@
-#pragma once
-
-#include "lsystem.hpp"
-
-// Predefined parameters for two classes of neurons.
-// See lsys_models.cc for details.
-
-extern lsys_param alpha_motoneuron_lsys;
-extern lsys_param purkinje_lsys;
diff --git a/lmorpho/lsystem.cpp b/lmorpho/lsystem.cpp
deleted file mode 100644
index ad3cf88ba0c359e5513dfe7fc00934982c71a8c2..0000000000000000000000000000000000000000
--- a/lmorpho/lsystem.cpp
+++ /dev/null
@@ -1,316 +0,0 @@
-#include <cmath>
-#include <iostream>
-#include <random>
-#include <stack>
-#include <vector>
-
-#include <arbor/morphology.hpp>
-#include <arbor/math.hpp>
-
-#include "lsystem.hpp"
-
-using namespace arb::math;
-
-using arb::section_geometry;
-using arb::section_point;
-using arb::morphology;
-
-// L-system implementation.
-
-// Random distribution used in the specification of L-system parameters.
-// It can be a constant, uniform over an interval, or truncated normal.
-// Note that mixture distributions of the above should be, but isn't yet,
-// implemented.
-
-struct lsys_distribution {
- using result_type = double;
- using param_type = lsys_distribution_param;
-
- lsys_distribution(const param_type& p): p_(p) { init(); }
-
- lsys_distribution(): p_() { init(); }
- lsys_distribution(result_type x): p_(x) { init(); }
- lsys_distribution(result_type lb, result_type ub): p_(lb, ub) { init(); }
- lsys_distribution(result_type lb, result_type ub, result_type mu, result_type sigma):
- p_(lb, ub, mu, sigma)
- {
- init();
- }
-
- param_type param() const { return p_; }
- void param(const param_type& pbis) { p_=pbis; }
-
- result_type min() const { return p_.lb; }
- result_type max() const {
- return p_.kind==param_type::constant? min(): p_.ub;
- }
-
- bool operator==(const lsys_distribution& y) const {
- return p_==y.p_;
- }
-
- bool operator!=(const lsys_distribution& y) const {
- return p_!=y.p_;
- }
-
- // omitting out ostream, istream methods
- // (...)
-
- template <typename Gen>
- result_type operator()(Gen& g, const param_type &p) const {
- lsys_distribution dist(p);
- return dist(g);
- }
-
- template <typename Gen>
- result_type operator()(Gen& g) const {
- switch (p_.kind) {
- case param_type::constant:
- return p_.lb;
- case param_type::uniform:
- return uniform_dist(g);
- case param_type::normal:
- // TODO: replace with a robust truncated normal generator!
- for (;;) {
- result_type r = normal_dist(g);
- if (r>=p_.lb && r<=p_.ub) return r;
- }
- }
- return 0;
- }
-
-private:
- param_type p_;
- mutable std::uniform_real_distribution<double> uniform_dist;
- mutable std::normal_distribution<double> normal_dist;
-
- void init() {
- switch (p_.kind) {
- case param_type::uniform:
- uniform_dist = std::uniform_real_distribution<double>(p_.lb, p_.ub);
- break;
- case param_type::normal:
- normal_dist = std::normal_distribution<result_type>(p_.mu, p_.sigma);
- break;
- default: ;
- }
- }
-};
-
-class burke_correlated_r {
- double a_;
- lsys_distribution r_;
-
-public:
- struct result_type { double rho1, rho2; };
-
- burke_correlated_r(double a, lsys_distribution_param r): a_(a), r_(r) {}
-
- template <typename Gen>
- result_type operator()(Gen& g) const {
- double r1 = r_(g);
- double r2 = r_(g);
- return { r1+a_*r2, r2+a_*r1 };
- }
-};
-
-class burke_exp_test {
- double k1_;
- double k2_;
-
-public:
- using result_type = int;
-
- burke_exp_test(double k1, double k2): k1_(k1), k2_(k2) {}
-
- template <typename Gen>
- result_type operator()(double delta_l, double radius, Gen& g) const {
- std::bernoulli_distribution B(delta_l*k1_*std::exp(k2_*radius*2.0));
- return B(g);
- }
-};
-
-class burke_exp2_test {
- double ak1_;
- double ak2_;
- double bk1_;
- double bk2_;
-
-public:
- using result_type = int;
-
- burke_exp2_test(double ak1, double ak2, double bk1, double bk2):
- ak1_(ak1), ak2_(ak2), bk1_(bk1), bk2_(bk2) {}
-
- template <typename Gen>
- result_type operator()(double delta_l, double radius, Gen& g) const {
- double diam = 2.*radius;
- std::bernoulli_distribution B(
- delta_l*std::min(ak1_*std::exp(ak2_*diam), bk1_*std::exp(bk2_*diam)));
-
- return B(g);
- }
-};
-
-// L-system parameter set with instantiated distributions.
-struct lsys_sampler {
- // Create distributions once from supplied parameter set.
- explicit lsys_sampler(const lsys_param& P):
- diam_soma(P.diam_soma),
- n_tree(P.n_tree),
- diam_initial(P.diam_initial),
- length_step(P.length_step),
- roll_initial(P.roll_initial),
- pitch_initial(P.pitch_initial),
- roll_section(P.roll_section),
- pitch_section(P.pitch_section),
- taper(P.taper),
- roll_at_branch(P.roll_at_branch),
- branch_angle(P.branch_angle),
- diam_child(P.diam_child_a, P.diam_child_r),
- term_test(P.pterm_k1, P.pterm_k2),
- branch_test(P.pbranch_ov_k1, P.pbranch_ov_k2,
- P.pbranch_nov_k1, P.pbranch_nov_k2),
- max_extent(P.max_extent)
- {}
-
- lsys_distribution diam_soma;
- lsys_distribution n_tree;
- lsys_distribution diam_initial;
- lsys_distribution length_step;
- lsys_distribution roll_initial;
- lsys_distribution pitch_initial;
- lsys_distribution roll_section;
- lsys_distribution pitch_section;
- lsys_distribution taper;
- lsys_distribution roll_at_branch;
- lsys_distribution branch_angle;
- burke_correlated_r diam_child;
- burke_exp_test term_test;
- burke_exp2_test branch_test;
- double max_extent;
-};
-
-struct section_tip {
- section_point point = {0., 0., 0., 0.};
- quaternion rotation = {1., 0., 0., 0.};
- double somatic_distance = 0.;
-};
-
-struct grow_result {
- std::vector<section_point> points;
- std::vector<section_tip> children;
- double length = 0.;
-};
-
-constexpr double deg_to_rad = 3.1415926535897932384626433832795/180.0;
-
-template <typename Gen>
-grow_result grow(section_tip tip, const lsys_sampler& S, Gen &g) {
- constexpr quaternion xaxis = {0, 1, 0, 0};
- static std::uniform_real_distribution<double> U;
-
- grow_result result;
- std::vector<section_point>& points = result.points;
-
- points.push_back(tip.point);
- for (;;) {
- quaternion step = xaxis^tip.rotation;
- double dl = S.length_step(g);
- tip.point.x += dl*step.x;
- tip.point.y += dl*step.y;
- tip.point.z += dl*step.z;
- tip.point.r += dl*0.5*S.taper(g);
- tip.somatic_distance += dl;
- result.length += dl;
-
- if (tip.point.r<0) tip.point.r = 0;
-
- double phi = S.roll_section(g);
- double theta = S.pitch_section(g);
- tip.rotation *= rotation_x(deg_to_rad*phi);
- tip.rotation *= rotation_y(deg_to_rad*theta);
-
- points.push_back(tip.point);
-
- if (tip.point.r==0 || tip.somatic_distance>=S.max_extent) {
- return result;
- }
-
- if (S.branch_test(dl, tip.point.r, g)) {
- auto r = S.diam_child(g);
- // ignore branch if we get non-positive radii
- if (r.rho1>0 && r.rho2>0) {
- double branch_phi = S.roll_at_branch(g);
- double branch_angle = S.branch_angle(g);
- double branch_theta1 = U(g)*branch_angle;
- double branch_theta2 = branch_theta1-branch_angle;
-
- tip.rotation *= rotation_x(deg_to_rad*branch_phi);
-
- section_tip t1 = tip;
- t1.point.r = t1.point.r * r.rho1;
- t1.rotation *= rotation_y(deg_to_rad*branch_theta1);
-
- section_tip t2 = tip;
- t2.point.r = t2.point.r * r.rho2;
- t2.rotation *= rotation_y(deg_to_rad*branch_theta2);
-
- result.children = {t1, t2};
- return result;
- }
- }
-
- if (S.term_test(dl, tip.point.r, g)) {
- return result;
- }
- }
-}
-
-arb::morphology generate_morphology(const lsys_param& P, lsys_generator &g) {
- constexpr quaternion xaxis = {0, 1, 0, 0};
- arb::morphology morph;
-
- lsys_sampler S(P);
- double soma_radius = 0.5*S.diam_soma(g);
- morph.soma = {0, 0, 0, soma_radius};
-
- struct section_start {
- section_tip tip;
- unsigned parent_id;
- };
- std::stack<section_start> starts;
- unsigned next_id = 1u;
-
- int n = (int)std::round(S.n_tree(g));
- for (int i=0; i<n; ++i) {
- double phi = S.roll_initial(g);
- double theta = S.pitch_initial(g);
- double radius = 0.5*S.diam_initial(g);
-
- section_tip tip;
- tip.rotation = rotation_x(deg_to_rad*phi)*rotation_y(deg_to_rad*theta);
- tip.somatic_distance = 0.0;
-
- auto p = (soma_radius*xaxis)^tip.rotation;
- tip.point = {p.x, p.y, p.z, radius};
-
- starts.push({tip, 0u});
- }
-
- while (!starts.empty() && morph.sections.size()<P.max_sections) {
- auto start = starts.top();
- starts.pop();
-
- auto branch = grow(start.tip, S, g);
- section_geometry section{next_id++, start.parent_id, branch.children.empty(), std::move(branch.points), branch.length};
-
- for (auto child: branch.children) {
- starts.push({child, section.id});
- }
- morph.sections.push_back(std::move(section));
- }
-
- return morph;
-}
-
diff --git a/lmorpho/lsystem.hpp b/lmorpho/lsystem.hpp
deleted file mode 100644
index 55e5beda31a9b909a9db6aec58765d179e97c45e..0000000000000000000000000000000000000000
--- a/lmorpho/lsystem.hpp
+++ /dev/null
@@ -1,131 +0,0 @@
-#pragma once
-
-#include <random>
-
-#include <arbor/morphology.hpp>
-
-struct lsys_param;
-
-using lsys_generator = std::minstd_rand;
-
-arb::morphology generate_morphology(const lsys_param& P, lsys_generator& g);
-
-// The distribution parameters used in the specification of the L-system parameters.
-// The distribution can be a constant, uniform over an interval, or truncated normal.
-// Note that mixture distributions of the above should be, but isn't yet,
-// implemented.
-
-struct lsys_distribution_param {
- using result_type = double;
- enum lsys_kind {
- constant, uniform, normal
- };
-
- lsys_kind kind;
- result_type lb; // lower bound
- result_type ub; // upper bound (used for only non-constant)
- result_type mu; // mean (used only for truncated Gaussian)
- result_type sigma; // s.d. (used only for truncated Gaussian)
-
- bool operator==(const lsys_distribution_param &p) const {
- if (kind!=p.kind) return false;
- switch (kind) {
- case constant:
- return lb == p.lb;
- case uniform:
- return lb == p.lb && ub == p.ub;
- case normal:
- return lb == p.lb && ub == p.ub && mu == p.mu && sigma == p.sigma;
- }
- return false;
- }
-
- bool operator!=(const lsys_distribution_param &p) const {
- return !(*this==p);
- }
-
- // zero or one parameter => constant
- lsys_distribution_param():
- kind(constant), lb(0.0) {}
-
- lsys_distribution_param(result_type x):
- kind(constant), lb(x) {}
-
- // two paramter => uniform
- lsys_distribution_param(result_type lb, result_type ub):
- kind(uniform), lb(lb), ub(ub) {}
-
- // four paramter => truncated normal
- lsys_distribution_param(result_type lb, result_type ub, result_type mu, result_type sigma):
- kind(normal), lb(lb), ub(ub), mu(mu), sigma(sigma) {}
-};
-
-
-// Parameters as referenced in Ascoli 2001 are given in parentheses below.
-// Defaults below are chosen to be safe, boring and not biological.
-struct lsys_param {
- // Soma diameter [µm].
- lsys_distribution_param diam_soma = { 20 };
-
- // Number of dendritic trees (rounded to nearest int). (Ntree)
- lsys_distribution_param n_tree = { 1 };
-
- // Initial dendrite diameter [µm]. (Dstem)
- lsys_distribution_param diam_initial = { 2 };
-
- // Dendrite step length [µm]. (ΔL)
- lsys_distribution_param length_step = { 25 };
-
- // Dendrites grow along x-axis after application of (instrinsic) rotations;
- // roll (about x-axis) is applied first, followed by pitch (about y-axis).
-
- // Initial roll (intrinsic rotation about x-axis) [degrees]. (Taz)
- lsys_distribution_param roll_initial = { 0 };
-
- // Initial pitch (intrinsic rotation about y-axis) [degrees]. (Tel)
- lsys_distribution_param pitch_initial = { 0 };
-
- // Tortuousness: roll within section over ΔL [degrees]. (Eaz)
- lsys_distribution_param roll_section = { 0 };
-
- // Tortuousness: pitch within section over ΔL [degrees]. (Eel)
- lsys_distribution_param pitch_section = { 0 };
-
- // Taper rate: diameter decrease per unit length. (TPRB)
- lsys_distribution_param taper = { 0 };
-
- // Branching torsion: roll at branch point [degrees]. (Btor)
- lsys_distribution_param roll_at_branch = { 0 };
-
- // Branch angle between siblings [degrees]. (Bamp)
- lsys_distribution_param branch_angle = { 45 };
-
- // Child branch diameter ratios are typically anticorrelated.
- // The Burke 1992 parameterization describes the two child ratios
- // as d1 = r1 + a*r2 and d2 = r2 + a*r1 where a is a constant
- // determined by the correlation, and r1 and r2 are drawn from
- // a common distribution (`d_child_r` below).
- double diam_child_a = 0;
- lsys_distribution_param diam_child_r = { 0.5 };
-
- // Bifurcation and termination probabilities per unit-length below are given by
- // P = k1 * exp(k2*diameter), described by the paremeters k1 and k2 [1/µm].
-
- // Termination probability parameters [1/µm]. (k1trm, k2trm)
- double pterm_k1 = 0.05;
- double pterm_k2 = -2;
-
- // Bifurcation probability is given by minimum of two probabilities, `pbranch_ov' and
- // `pbranch_nov' below.
- double pbranch_ov_k1 = 0.01;
- double pbranch_ov_k2 = 0;
- double pbranch_nov_k1 = 0.01;
- double pbranch_nov_k2 = 0;
-
- // Absolute maximum dendritic extent [µm]. (Forces termination of algorithm)
- double max_extent = 2000;
-
- // Absolute maximum number of unbranched sections. (Forces termination of algorithm)
- unsigned max_sections = 10000;
-};
-
diff --git a/lmorpho/morphio.cpp b/lmorpho/morphio.cpp
deleted file mode 100644
index e556f172c54e4c7d2654a7374e9b8ac970504009..0000000000000000000000000000000000000000
--- a/lmorpho/morphio.cpp
+++ /dev/null
@@ -1,167 +0,0 @@
-#include <fstream>
-#include <iomanip>
-#include <iterator>
-#include <map>
-#include <ostream>
-#include <sstream>
-#include <string>
-#include <vector>
-
-#include <arbor/morphology.hpp>
-#include <arbor/swcio.hpp>
-
-#include "morphio.hpp"
-
-using arb::swc_record;
-
-std::vector<swc_record> as_swc(const arb::morphology& morph);
-
-// Multi-file manager implementation.
-multi_file::multi_file(const std::string& pattern, int digits) {
- auto npos = std::string::npos;
-
- file_.exceptions(std::ofstream::failbit);
- concat_ = (pattern.find("%")==npos);
- use_stdout_ = pattern.empty() || pattern=="-";
-
- if (!concat_) {
- auto p = pattern.find("%");
- fmt_prefix_ = pattern.substr(0, p);
- fmt_suffix_ = pattern.substr(p+1);
- fmt_digits_ = digits;
- }
- else {
- filename_ = pattern;
- }
-}
-
-void multi_file::open(unsigned n) {
- if (use_stdout_ || (file_.is_open() && (concat_ || n==current_n_))) {
- return;
- }
-
- if (file_.is_open()) file_.close();
-
- std::string fname;
- if (concat_) {
- fname = filename_;
- }
- else {
- std::stringstream ss;
- ss << fmt_prefix_ << std::setfill('0') << std::setw(fmt_digits_) << n << fmt_suffix_;
- fname = ss.str();
- }
-
- file_.open(fname);
-
- current_n_ = n;
-}
-
-static std::string short_cable_message(int id, unsigned sz) {
- std::stringstream ss;
- ss << "surprisingly short cable: id=" << id << ", size=" << sz;
- return ss.str();
-}
-
-// SWC transform
-
-// TODO: Move this functionality to arbor library.
-std::vector<swc_record> as_swc(const arb::morphology& morph) {
- using kind = swc_record::kind;
- std::map<int, int> parent_end_id;
- std::vector<swc_record> swc;
-
- // soma
- const auto &p = morph.soma;
- int id = 0;
- parent_end_id[0] = 0;
- swc.emplace_back(kind::soma, id, p.x, p.y, p.z, p.r, -1);
-
- // dendrites:
- for (auto& sec: morph.sections) {
- int parent = parent_end_id[sec.parent_id];
-
- const auto& points = sec.points;
- auto n = points.size();
- if (n<2) {
- throw std::runtime_error(short_cable_message(sec.id, n));
- }
-
- // Include first point only for dendrites segments attached to soma.
- if (sec.parent_id==0) {
- const auto& p = points[0];
- ++id;
- swc.emplace_back(kind::fork_point, id, p.x, p.y, p.z, p.r, parent);
- parent = id;
- }
-
- // Interior points.
- for (unsigned i = 1; i<n-1; ++i) {
- const auto& p = points[i];
- ++id;
- swc.emplace_back(kind::dendrite, id, p.x, p.y, p.z, p.r, parent);
- parent = id;
- }
-
- // Final point (fork or terminal).
- const auto& p = points.back();
- ++id;
- swc.emplace_back(sec.terminal? kind::end_point: kind::fork_point, id, p.x, p.y, p.z, p.r, parent);
- parent_end_id[sec.id] = id;
- }
-
- return swc;
-}
-
-// SWC emitter implementation.
-
-void swc_emitter::operator()(unsigned index, const arb::morphology& m) {
- file_.open(index);
- auto& stream = file_.stream();
-
- auto swc = as_swc(m);
- stream << "# lmorpho generated morphology\n# index: " << index << "\n";
- std::copy(swc.begin(), swc.end(), std::ostream_iterator<swc_record>(stream, "\n"));
-}
-
-// pvector emitter implementation.
-
-std::vector<int> as_pvector(const arb::morphology& morph, unsigned offset) {
- std::map<int, unsigned> parent_index; // section id to segment index
- std::vector<int> pvec;
- unsigned index = offset; // starting segment index
-
- // soma
- parent_index[0] = index;
- pvec.push_back(-1);
- ++index;
-
- // dendrites:
- for (auto& sec: morph.sections) {
- int parent = parent_index[sec.parent_id];
-
- auto n = sec.points.size();
- if (n<2) {
- throw std::runtime_error(short_cable_message(sec.id, n));
- }
-
- for (unsigned i = 1; i<n; ++i) {
- pvec.push_back(parent);
- parent = index++;
- }
-
- parent_index[sec.id] = parent;
- }
-
- return pvec;
-}
-
-void pvector_emitter::operator()(unsigned index, const arb::morphology& m) {
- auto pvec = as_pvector(m, offset_);
- if (coalesce_) offset_ += pvec.size();
-
- file_.open(index);
- auto& stream = file_.stream();
- std::copy(pvec.begin(), pvec.end(), std::ostream_iterator<int>(stream, "\n"));
-}
-
diff --git a/lmorpho/morphio.hpp b/lmorpho/morphio.hpp
deleted file mode 100644
index 23a2854bead20fa2f796c4b1b340abcf9fd62dcc..0000000000000000000000000000000000000000
--- a/lmorpho/morphio.hpp
+++ /dev/null
@@ -1,97 +0,0 @@
-#pragma once
-
-#include <fstream>
-#include <iostream>
-#include <vector>
-
-#include <arbor/morphology.hpp>
-
-// Manage access to a single file, std::cout, or an indexed
-// sequence of files.
-
-class multi_file {
-private:
- std::ofstream file_;
- bool concat_ = false;
- bool use_stdout_ = false;
- // use if not concat_:
- std::string fmt_prefix_;
- std::string fmt_suffix_;
- int fmt_digits_ = 0;
- // use if concat_:
- std::string filename_; // use if concat_
- unsigned current_n_ = 0;
-
-public:
- // If pattern is '-' or empty, represent std::cout.
- // If pattern contains '%', use a different file for each index.
- // Otherwise, use a single file given by pattern.
- explicit multi_file(const std::string& pattern, int digits=0);
-
- multi_file(multi_file&&) = default;
-
- // Open file for index n, closing previously open file if different.
- void open(unsigned n);
-
- // Close currently open file.
- void close() { if (!use_stdout_ && file_.is_open()) file_.close(); }
-
- // True if writing to a single file or std::cout.
- bool single_file() const { return concat_; }
-
- // Return std::cout or file stream as appropriate.
- std::ostream& stream() { return use_stdout_? std::cout: file_; }
-};
-
-static constexpr unsigned digits(unsigned n) {
- return n? 1+digits(n/10): 0;
-}
-
-// Write a sequence of morphologies to one or more files
-// as given my `pattern`.
-
-class swc_emitter {
- multi_file file_;
-
-public:
- // `pattern` is as for `multi_file`; number `n` optionally
- // specifies the total number of morphologies, for better
- // numeric formatting.
- explicit swc_emitter(std::string pattern, unsigned n=0):
- file_(pattern, digits(n)) {}
-
- swc_emitter(swc_emitter&&) = default;
-
- // write `index`th morphology as SWC.
- void operator()(unsigned index, const arb::morphology& m);
-
- void close() { file_.close(); }
- ~swc_emitter() { close(); }
-};
-
-// Write pvectors for a sequence of morphologies to one or more files
-// as given my `pattern`. Coalesce the vectors if writing to a single
-// file/stream.
-
-class pvector_emitter {
- multi_file file_;
- bool coalesce_ = false;
- unsigned offset_ = 0;
-
-public:
- // `pattern` is as for `multi_file`; number `n` optionally
- // specifies the total number of morphologies, for better
- // numeric formatting.
- explicit pvector_emitter(std::string pattern, unsigned n=0):
- file_(pattern, digits(n)),
- coalesce_(file_.single_file()) {}
-
- pvector_emitter(pvector_emitter&&) = default;
-
- // write pvector for `index`th morphology.
- void operator()(unsigned index, const arb::morphology& m);
-
- void close() { file_.close(); }
- ~pvector_emitter() { close(); }
-};
-
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 2fb9912b8db2bb04f4eb0af11e7ebb4271e390a1..1e042fbb902b89aa858127b349cbbd75a4208fa0 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -98,6 +98,7 @@ set(unit_sources
test_multi_event_stream.cpp
test_optional.cpp
test_mechinfo.cpp
+ test_morphology.cpp
test_padded.cpp
test_partition.cpp
test_partition_by_constraint.cpp
diff --git a/test/unit/test_counter.cpp b/test/unit/test_counter.cpp
index f58c67cefc6b7cbf5afdc0000ac6894705c23d96..83b6dc2eff602edd8087538bf9dd43ffb857cc06 100644
--- a/test/unit/test_counter.cpp
+++ b/test/unit/test_counter.cpp
@@ -3,7 +3,7 @@
#include <iterator>
#include <type_traits>
-#include <util/counter.hpp>
+#include "util/counter.hpp"
using namespace arb;
diff --git a/test/unit/test_morphology.cpp b/test/unit/test_morphology.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f0af3a4a80e78121496013f73f6009e892f0e90b
--- /dev/null
+++ b/test/unit/test_morphology.cpp
@@ -0,0 +1,561 @@
+/*
+ * Unit tests for sample_tree and morphology.
+ */
+
+#include <fstream>
+#include <cassert>
+#include <cmath>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+#include "../test/gtest.h"
+
+#include <arbor/math.hpp>
+#include <arbor/morph/error.hpp>
+#include <arbor/morph/morphology.hpp>
+#include <arbor/morph/sample_tree.hpp>
+#include <arbor/cable_cell.hpp>
+
+#include "algorithms.hpp"
+#include "io/sepval.hpp"
+#include "util/range.hpp"
+#include "util/span.hpp"
+#include "util/strprintf.hpp"
+
+// Forward declare non-public functions that are used internally to build
+// morphologies so that they can be tested.
+namespace arb { namespace impl {
+ std::vector<mbranch> branches_from_parent_index(const std::vector<arb::msize_t>&, const std::vector<point_prop>&, bool);
+}}
+
+template <typename T>
+std::ostream& operator<<(std::ostream& o, const std::vector<T>& v) {
+ return o << "[" << arb::io::csv(v) << "]";
+}
+
+TEST(morphology, point_props) {
+ arb::point_prop p = arb::point_prop_mask_none;
+
+ EXPECT_FALSE(arb::is_terminal(p));
+ EXPECT_FALSE(arb::is_fork(p));
+ EXPECT_FALSE(arb::is_root(p));
+ EXPECT_FALSE(arb::is_collocated(p));
+
+ arb::set_root(p);
+ EXPECT_FALSE(arb::is_terminal(p));
+ EXPECT_FALSE(arb::is_fork(p));
+ EXPECT_TRUE(arb::is_root(p));
+ EXPECT_FALSE(arb::is_collocated(p));
+
+ arb::set_terminal(p);
+ EXPECT_TRUE(arb::is_terminal(p));
+ EXPECT_FALSE(arb::is_fork(p));
+ EXPECT_TRUE(arb::is_root(p));
+ EXPECT_FALSE(arb::is_collocated(p));
+
+ arb::unset_root(p);
+ EXPECT_TRUE(arb::is_terminal(p));
+ EXPECT_FALSE(arb::is_fork(p));
+ EXPECT_FALSE(arb::is_root(p));
+ EXPECT_FALSE(arb::is_collocated(p));
+
+ arb::set_collocated(p);
+ EXPECT_TRUE(arb::is_terminal(p));
+ EXPECT_FALSE(arb::is_fork(p));
+ EXPECT_FALSE(arb::is_root(p));
+ EXPECT_TRUE(arb::is_collocated(p));
+
+ arb::set_fork(p);
+ EXPECT_TRUE(arb::is_terminal(p));
+ EXPECT_TRUE(arb::is_fork(p));
+ EXPECT_FALSE(arb::is_root(p));
+ EXPECT_TRUE(arb::is_collocated(p));
+
+ arb::unset_fork(p);
+ arb::unset_terminal(p);
+ arb::unset_collocated(p);
+ EXPECT_FALSE(arb::is_terminal(p));
+ EXPECT_FALSE(arb::is_fork(p));
+ EXPECT_FALSE(arb::is_root(p));
+ EXPECT_FALSE(arb::is_collocated(p));
+}
+
+// TODO: test sample_tree marks properties correctly
+
+// Test internal function that parses a parent list, and marks
+// each node as either root, sequential, fork or terminal.
+TEST(sample_tree, properties) {
+ const auto npos = arb::mnpos;
+ using arb::sample_tree;
+ using pp = arb::point_prop;
+ using pvec = std::vector<arb::msize_t>;
+
+ pp c = arb::point_prop_mask_collocated;
+ pp r = arb::point_prop_mask_root;
+ pp t = arb::point_prop_mask_terminal;
+ pp s = arb::point_prop_mask_none;
+ pp f = arb::point_prop_mask_fork;
+ pp tc = t+c;
+ pp sc = s+c;
+ pp fc = f+c;
+
+ // make a sample tree from a parent vector with non-collocated points.
+ auto make_tree = [] (const pvec& parents) {
+ sample_tree st;
+ for (auto p: parents) st.append(p, {{0.,0.,double(st.size()),1.}, 1});
+ return st;
+ };
+ // make a sample tree from a parent vector with collocated points.
+ auto make_colloc_tree = [] (const pvec& parents) {
+ sample_tree st;
+ for (auto p: parents) st.append(p, {{0.,0.,0.,1.}, 1});
+ return st;
+ };
+
+ {
+ EXPECT_EQ(make_tree({npos}).properties(), std::vector<pp>{r});
+ EXPECT_EQ(make_colloc_tree({npos}).properties(), std::vector<pp>{r});
+ }
+
+ {
+ EXPECT_EQ(make_tree({npos,0}).properties(), (std::vector<pp>{r,t}));
+ EXPECT_EQ(make_colloc_tree({npos,0}).properties(), (std::vector<pp>{r,tc}));
+ }
+
+ {
+ EXPECT_EQ(make_tree({npos,0,1,2}).properties(), (std::vector<pp>{r,s,s,t}));
+ EXPECT_EQ(make_colloc_tree({npos,0,1,2}).properties(), (std::vector<pp>{r,sc,sc,tc}));
+ }
+
+ {
+ EXPECT_EQ(make_tree({npos,0,1,2,0,4,5}).properties(), (std::vector<pp>{r,s,s,t,s,s,t}));
+ EXPECT_EQ(make_colloc_tree({npos,0,1,2,0,4,5}).properties(), (std::vector<pp>{r,sc,sc,tc,sc,sc,tc}));
+ }
+
+ {
+ EXPECT_EQ(make_tree({npos,0,1,2,3,2,4,4,7}).properties(), (std::vector<pp>{r,s,f,s,f,t,t,s,t}));
+ EXPECT_EQ(make_colloc_tree({npos,0,1,2,3,2,4,4,7}).properties(), (std::vector<pp>{r,sc,fc,sc,fc,tc,tc,sc,tc}));
+ }
+}
+
+TEST(morphology, branches_from_parent_index) {
+ const auto npos = arb::mnpos;
+ using pvec = std::vector<arb::msize_t>;
+ using mb = arb::mbranch;
+
+ // make a sample tree from a parent vector with non-collocated points.
+ auto make_tree = [] (const pvec& parents) {
+ arb::sample_tree st;
+ for (auto p: parents) st.append(p, {{0.,0.,double(st.size()),1.}, 1});
+ return st;
+ };
+
+ { // single sample: can only be used to build a morphology with one spherical branch
+ pvec parents{npos};
+ auto tree = make_tree(parents);
+ auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), true);
+ EXPECT_EQ(1u, bc.size());
+ EXPECT_EQ(mb({0}, npos), bc[0]);
+
+ // A cable morphology can't be constructed from a single sample.
+ EXPECT_THROW(arb::impl::branches_from_parent_index(parents, tree.properties(), false),
+ arb::morphology_error);
+ }
+ {
+ pvec parents = {npos, 0};
+ auto tree = make_tree(parents);
+ auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
+ EXPECT_EQ(1u, bc.size());
+ EXPECT_EQ(mb({0,1}, npos), bc[0]);
+
+ // A morphology can't be constructed with a spherical soma from two samples.
+ EXPECT_THROW(arb::impl::branches_from_parent_index(parents, tree.properties(), true),
+ arb::morphology_error);
+ }
+
+ {
+ pvec parents{npos, 0, 1};
+
+ // With cable soma: one cable with 3 samples.
+ auto tree = make_tree(parents);
+ auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
+ EXPECT_EQ(1u, bc.size());
+ EXPECT_EQ(mb({0,1,2},npos), bc[0]);
+
+ // With spherical soma: one sphere and a 2-segment cable.
+ // The cable branch is attached to the sphere (i.e. the sphere is the parent branch).
+ auto bs = arb::impl::branches_from_parent_index(parents, tree.properties(), true);
+ EXPECT_EQ(2u, bs.size());
+ EXPECT_EQ(mb({0},npos), bs[0]);
+ EXPECT_EQ(mb({1,2},0), bs[1]);
+ }
+
+ {
+ pvec parents{npos, 0, 0};
+ auto tree = make_tree(parents);
+
+ // A spherical root is not valid: each cable branch would have only one sample.
+ EXPECT_THROW(arb::impl::branches_from_parent_index(parents, tree.properties(), true), arb::morphology_error);
+
+ // Two cables, with two samples each, with the first sample in each being the root
+ auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
+ EXPECT_EQ(2u, bc.size());
+ EXPECT_EQ(mb({0,1},npos), bc[0]);
+ EXPECT_EQ(mb({0,2},npos), bc[1]);
+ }
+
+ {
+ pvec parents{npos, 0, 1, 2};
+ auto tree = make_tree(parents);
+
+ // With cable soma: one cable with 4 samples.
+ auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
+ EXPECT_EQ(1u, bc.size());
+ EXPECT_EQ(mb({0,1,2,3},npos), bc[0]);
+
+ // With spherical soma: one sphere and on 3-segment cable.
+ // The cable branch is attached to the sphere (i.e. the sphere is the parent branch).
+ auto bs = arb::impl::branches_from_parent_index(parents, tree.properties(), true);
+ EXPECT_EQ(2u, bs.size());
+ EXPECT_EQ(mb({0},npos), bs[0]);
+ EXPECT_EQ(mb({1,2,3},0), bs[1]);
+ }
+
+ {
+ pvec parents{npos, 0, 1, 0};
+ auto tree = make_tree(parents);
+
+ // With cable soma: two cables with 3 and 2 samples respectively.
+ auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
+ EXPECT_EQ(2u, bc.size());
+ EXPECT_EQ(mb({0,1,2},npos), bc[0]);
+ EXPECT_EQ(mb({0,3},npos), bc[1]);
+
+ // A spherical root is not valid: the second cable branch would have only one sample.
+ EXPECT_THROW(arb::impl::branches_from_parent_index(parents, tree.properties(), true), arb::morphology_error);
+ }
+
+ {
+ pvec parents{npos, 0, 1, 0, 3};
+ auto tree = make_tree(parents);
+
+ // With cable soma: two cables with 3 samples each [0,1,2] and [0,3,4]
+ auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
+ EXPECT_EQ(2u, bc.size());
+ EXPECT_EQ(mb({0,1,2},npos), bc[0]);
+ EXPECT_EQ(mb({0,3,4},npos), bc[1]);
+
+ // With spherical soma: one sphere and 2 2-sample cables.
+ // each of the cable branches is attached to the sphere (i.e. the sphere is the parent branch).
+ auto bs = arb::impl::branches_from_parent_index(parents, tree.properties(), true);
+ EXPECT_EQ(3u, bs.size());
+ EXPECT_EQ(mb({0},npos), bs[0]);
+ EXPECT_EQ(mb({1,2},0), bs[1]);
+ EXPECT_EQ(mb({3,4},0), bs[2]);
+ }
+
+ {
+ pvec parents{npos, 0, 1, 0, 3, 4, 4, 6};
+ auto tree = make_tree(parents);
+
+ // With cable soma: 4 cables: [0,1,2] [0,3,4] [4,5] [4,6,7]
+ auto bc = arb::impl::branches_from_parent_index(parents, tree.properties(), false);
+ EXPECT_EQ(4u, bc.size());
+ EXPECT_EQ(mb({0,1,2},npos), bc[0]);
+ EXPECT_EQ(mb({0,3,4},npos), bc[1]);
+ EXPECT_EQ(mb({4,5}, 1), bc[2]);
+ EXPECT_EQ(mb({4,6,7},1), bc[3]);
+
+ // With spherical soma: 1 sphere and 4 cables: [1,2] [3,4] [4,5] [4,6,7]
+ auto bs = arb::impl::branches_from_parent_index(parents, tree.properties(), true);
+ EXPECT_EQ(5u, bs.size());
+ EXPECT_EQ(mb({0}, npos), bs[0]);
+ EXPECT_EQ(mb({1,2}, 0), bs[1]);
+ EXPECT_EQ(mb({3,4}, 0), bs[2]);
+ EXPECT_EQ(mb({4,5}, 2), bs[3]);
+ EXPECT_EQ(mb({4,6,7}, 2), bs[4]);
+ }
+}
+
+// For different parent index vectors, attempt multiple valid and invalid sample sets.
+TEST(morphology, construction) {
+ constexpr auto npos = arb::mnpos;
+ using arb::util::make_span;
+ using ms = arb::msample;
+ using pvec = std::vector<arb::msize_t>;
+ {
+ pvec p = {npos, 0};
+ std::vector<ms> s = {
+ {{0.0, 0.0, 0.0, 1.0}, 1},
+ {{0.0, 0.0, 1.0, 1.0}, 1} };
+
+ arb::sample_tree sm(s, p);
+ auto m = arb::morphology(sm);
+
+ EXPECT_EQ(1u, m.num_branches());
+ }
+ {
+ pvec p = {npos, 0, 1};
+ { // 2-segment cable (1 seg soma, 1 seg dendrite)
+ std::vector<ms> s = {
+ {{0.0, 0.0, 0.0, 5.0}, 1},
+ {{0.0, 0.0, 5.0, 1.0}, 1},
+ {{0.0, 0.0, 8.0, 1.0}, 2} };
+
+ arb::sample_tree sm(s, p);
+ auto m = arb::morphology(sm);
+
+ EXPECT_EQ(1u, m.num_branches());
+ }
+ { // spherical soma and single-segment cable
+ std::vector<ms> s = {
+ {{0.0, 0.0, 0.0, 5.0}, 1},
+ {{0.0, 0.0, 1.0, 1.0}, 2},
+ {{0.0, 0.0, 8.0, 1.0}, 2} };
+
+ arb::sample_tree sm(s, p);
+ auto m = arb::morphology(sm);
+
+ EXPECT_EQ(2u, m.num_branches());
+ }
+ }
+ {
+ // 0 |
+ // 1 3 |
+ // 2 |
+ pvec p = {npos, 0, 1, 0};
+ {
+ // two cables: 1x2 segments, 1x1 segment.
+ std::vector<ms> s = {
+ {{0.0, 0.0, 0.0, 5.0}, 1},
+ {{0.0, 0.0, 5.0, 1.0}, 1},
+ {{0.0, 0.0, 6.0, 1.0}, 2},
+ {{0.0, 4.0, 0.0, 1.0}, 1}};
+
+ arb::sample_tree sm(s, p);
+ auto m = arb::morphology(sm);
+
+ EXPECT_EQ(2u, m.num_branches());
+ }
+ {
+ // error: spherical soma with a single point cable attached via sample 3
+ std::vector<ms> s = {
+ {{0.0, 0.0, 0.0, 5.0}, 1},
+ {{0.0, 0.0, 5.0, 1.0}, 2},
+ {{0.0, 0.0, 8.0, 1.0}, 2},
+ {{0.0, 5.0, 0.0, 1.0}, 2}};
+
+ arb::sample_tree sm(s, p);
+ EXPECT_THROW((arb::morphology(sm)), arb::morphology_error);
+ }
+ }
+ {
+ // 0 |
+ // 1 3 |
+ // 2 4 |
+ pvec p = {npos, 0, 1, 0, 3};
+ {
+ // two cables: 1x2 segments, 1x1 segment.
+ std::vector<ms> s = {
+ {{0.0, 0.0, 0.0, 5.0}, 1},
+ {{0.0, 0.0, 5.0, 1.0}, 2},
+ {{0.0, 0.0, 8.0, 1.0}, 2},
+ {{0.0, 5.0, 0.0, 1.0}, 2},
+ {{0.0, 8.0, 0.0, 1.0}, 2}};
+
+ arb::sample_tree sm(s, p);
+ auto m = arb::morphology(sm);
+
+ EXPECT_EQ(3u, m.num_branches());
+ }
+ }
+}
+
+// test that morphology generates branch child-parent structure correctly.
+TEST(morphology, branches) {
+ using pvec = std::vector<arb::msize_t>;
+ using svec = std::vector<arb::msample>;
+ auto npos = arb::mnpos;
+ {
+ // 0
+ pvec parents = {npos};
+ svec samples = {
+ {{ 0,0,0,3}, 1},
+ };
+ arb::sample_tree sm(samples, parents);
+ arb::morphology m(sm);
+
+ EXPECT_EQ(1u, m.num_branches());
+ EXPECT_EQ(npos, m.branch_parent(0));
+ EXPECT_EQ(pvec{}, m.branch_children(0));
+ }
+ {
+ // 0 - 1
+ pvec parents = {npos, 0};
+ svec samples = {
+ {{ 0,0,0,3}, 1},
+ {{ 10,0,0,3}, 1},
+ };
+ arb::sample_tree sm(samples, parents);
+ arb::morphology m(sm);
+
+ EXPECT_EQ(1u, m.num_branches());
+ EXPECT_EQ(npos, m.branch_parent(0));
+ EXPECT_EQ(pvec{}, m.branch_children(0));
+ }
+ {
+ // 0 - 1 - 2
+ pvec parents = {npos, 0, 1};
+ {
+ // All samples have same tag -> the morphology is a single unbranched cable.
+ svec samples = {
+ {{ 0,0,0,3}, 1},
+ {{10,0,0,3}, 1},
+ {{100,0,0,3}, 1},
+ };
+ arb::sample_tree sm(samples, parents);
+ arb::morphology m(sm);
+
+ EXPECT_EQ(1u, m.num_branches());
+ EXPECT_EQ(npos, m.branch_parent(0));
+ EXPECT_EQ(pvec{}, m.branch_children(0));
+ }
+ {
+ // First sample has unique tag -> spherical soma attached to a single-segment cable.
+ svec samples = {
+ {{ 0,0,0,10}, 1},
+ {{ 10,0,0, 3}, 3},
+ {{100,0,0, 3}, 3},
+ };
+ arb::sample_tree sm(samples, parents);
+ arb::morphology m(sm);
+
+ EXPECT_EQ(2u, m.num_branches());
+ EXPECT_EQ(npos, m.branch_parent(0));
+ EXPECT_EQ(0u, m.branch_parent(1));
+ EXPECT_EQ(pvec{1}, m.branch_children(0));
+ EXPECT_EQ(pvec{}, m.branch_children(1));
+ }
+ }
+ {
+ // 2 - 0 - 1
+ pvec parents = {npos, 0, 0};
+
+ svec samples = {
+ {{ 0, 0,0, 5}, 3},
+ {{ 10, 0,0, 5}, 3},
+ {{ 0,10,0, 5}, 3},
+ };
+ arb::sample_tree sm(samples, parents);
+ arb::morphology m(sm);
+
+ EXPECT_EQ(2u, m.num_branches());
+ EXPECT_EQ(npos, m.branch_parent(0));
+ EXPECT_EQ(npos, m.branch_parent(1));
+ EXPECT_EQ(pvec{}, m.branch_children(0));
+ EXPECT_EQ(pvec{}, m.branch_children(1));
+ }
+ {
+ // 0 - 1 - 2 - 3
+ pvec parents = {npos, 0, 1, 2};
+ }
+ {
+ // 0 |
+ // / \ |
+ // 1 3 |
+ // / |
+ // 2 |
+ pvec parents = {npos, 0, 1, 0};
+ }
+ {
+ // Eight samples
+ //
+ // 0 |
+ // / \ |
+ // 1 3 |
+ // / \ |
+ // 2 4 |
+ // / \ |
+ // 5 6 |
+ // \ |
+ // 7 |
+ pvec parents = {npos, 0, 1, 0, 3, 4, 4, 6};
+ {
+ svec samples = {
+ {{ 0, 0, 0, 10}, 1},
+ {{ 10, 0, 0, 2}, 3},
+ {{100, 0, 0, 2}, 3},
+ {{ 0, 10, 0, 2}, 3},
+ {{ 0,100, 0, 2}, 3},
+ {{100,100, 0, 2}, 3},
+ {{ 0,200, 0, 2}, 3},
+ {{ 0,300, 0, 2}, 3},
+ };
+ arb::sample_tree sm(samples, parents);
+ arb::morphology m(sm);
+
+ EXPECT_EQ(5u, m.num_branches());
+ EXPECT_EQ(npos, m.branch_parent(0));
+ EXPECT_EQ(0u, m.branch_parent(1));
+ EXPECT_EQ(0u, m.branch_parent(2));
+ EXPECT_EQ(2u, m.branch_parent(3));
+ EXPECT_EQ(2u, m.branch_parent(4));
+ EXPECT_EQ((pvec{1,2}), m.branch_children(0));
+ EXPECT_EQ((pvec{}), m.branch_children(1));
+ EXPECT_EQ((pvec{3,4}), m.branch_children(2));
+ EXPECT_EQ((pvec{}), m.branch_children(3));
+ EXPECT_EQ((pvec{}), m.branch_children(4));
+ }
+ {
+ svec samples = {
+ {{ 0, 0, 0, 10}, 3},
+ {{ 10, 0, 0, 2}, 3},
+ {{100, 0, 0, 2}, 3},
+ {{ 0, 10, 0, 2}, 3},
+ {{ 0,100, 0, 2}, 3},
+ {{100,100, 0, 2}, 3},
+ {{ 0,200, 0, 2}, 3},
+ {{ 0,300, 0, 2}, 3},
+ };
+ arb::sample_tree sm(samples, parents);
+ arb::morphology m(sm);
+
+ EXPECT_EQ(4u, m.num_branches());
+ EXPECT_EQ(npos, m.branch_parent(0));
+ EXPECT_EQ(npos, m.branch_parent(1));
+ EXPECT_EQ(1u, m.branch_parent(2));
+ EXPECT_EQ(1u, m.branch_parent(3));
+ EXPECT_EQ((pvec{}), m.branch_children(0));
+ EXPECT_EQ((pvec{2,3}), m.branch_children(1));
+ EXPECT_EQ((pvec{}), m.branch_children(2));
+ EXPECT_EQ((pvec{}), m.branch_children(3));
+ }
+ }
+}
+
+TEST(morphology, swc) {
+ std::string datadir{DATADIR};
+ auto fname = datadir + "/example.swc";
+ std::ifstream fid(fname);
+ if (!fid.is_open()) {
+ std::cerr << "unable to open file " << fname << "... skipping test\n";
+ return;
+ }
+
+ // Load swc samples from file.
+ auto swc_samples = arb::parse_swc_file(fid);
+
+ // Build a sample_tree from swc samples.
+ auto sm = arb::swc_as_sample_tree(swc_samples);
+ EXPECT_EQ(1058u, sm.size()); // file contains 195 samples
+
+ // Test that the morphology contains the expected number of branches.
+ auto m = arb::morphology(sm);
+ 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());
+}
+
diff --git a/test/unit/test_range.cpp b/test/unit/test_range.cpp
index c8c0a82e4747bf7c17013de3dd7a269c6fdbae7d..b15c86abb54523ed7c38a262c822d1d2752c3fad 100644
--- a/test/unit/test_range.cpp
+++ b/test/unit/test_range.cpp
@@ -9,12 +9,12 @@
#include <type_traits>
#include <unordered_map>
-#include <util/counter.hpp>
-#include <util/meta.hpp>
-#include <util/range.hpp>
-#include <util/rangeutil.hpp>
-#include <util/sentinel.hpp>
-#include <util/transform.hpp>
+#include "util/counter.hpp"
+#include "util/meta.hpp"
+#include "util/range.hpp"
+#include "util/rangeutil.hpp"
+#include "util/sentinel.hpp"
+#include "util/transform.hpp"
#include "common.hpp"
diff --git a/test/unit/test_swcio.cpp b/test/unit/test_swcio.cpp
index 3a22d257cc3fdc29f310fc7d31f98e140d9ea7eb..c12fd36cc4c4e61c2e04438c354e3d1128867083 100644
--- a/test/unit/test_swcio.cpp
+++ b/test/unit/test_swcio.cpp
@@ -8,7 +8,6 @@
#include <vector>
#include <arbor/cable_cell.hpp>
-#include <arbor/morphology.hpp>
#include <arbor/swcio.hpp>
#include "../gtest.h"
@@ -26,7 +25,7 @@ void expect_record_equals(const swc_record& expected,
const swc_record& actual)
{
EXPECT_EQ(expected.id, actual.id);
- EXPECT_EQ(expected.type, actual.type);
+ EXPECT_EQ(expected.tag, actual.tag);
EXPECT_FLOAT_EQ(expected.x, actual.x);
EXPECT_FLOAT_EQ(expected.y, actual.y);
EXPECT_FLOAT_EQ(expected.z, actual.z);
@@ -36,54 +35,44 @@ void expect_record_equals(const swc_record& expected,
TEST(swc_record, construction)
{
- {
- // force an invalid type
- swc_record::kind invalid_type = static_cast<swc_record::kind>(100);
- EXPECT_THROW(swc_record(invalid_type, 7, 1., 1., 1., 1., 5).assert_consistent(),
- swc_error);
- }
+ int soma_tag = 1;
{
// invalid id
- EXPECT_THROW(swc_record(
- swc_record::kind::custom, -3, 1., 1., 1., 1., 5).assert_consistent(),
+ EXPECT_THROW(swc_record(soma_tag, -3, 1., 1., 1., 1., 5).assert_consistent(),
swc_error);
}
{
// invalid parent id
- EXPECT_THROW(swc_record(
- swc_record::kind::custom, 0, 1., 1., 1., 1., -5).assert_consistent(),
+ EXPECT_THROW(swc_record(soma_tag, 0, 1., 1., 1., 1., -5).assert_consistent(),
swc_error);
}
{
// invalid radius
- EXPECT_THROW(swc_record(
- swc_record::kind::custom, 0, 1., 1., 1., -1., -1).assert_consistent(),
+ EXPECT_THROW(swc_record(soma_tag, 0, 1., 1., 1., -1., -1).assert_consistent(),
swc_error);
}
{
// parent_id > id
- EXPECT_THROW(swc_record(
- swc_record::kind::custom, 0, 1., 1., 1., 1., 2).assert_consistent(),
+ EXPECT_THROW(swc_record(soma_tag, 0, 1., 1., 1., 1., 2).assert_consistent(),
swc_error);
}
{
// parent_id == id
- EXPECT_THROW(swc_record(
- swc_record::kind::custom, 0, 1., 1., 1., 1., 0).assert_consistent(),
+ EXPECT_THROW(swc_record(soma_tag, 0, 1., 1., 1., 1., 0).assert_consistent(),
swc_error);
}
{
// check standard construction by value
- swc_record record(swc_record::kind::custom, 0, 1., 1., 1., 1., -1);
+ swc_record record(soma_tag, 0, 1., 1., 1., 1., -1);
EXPECT_TRUE(record.is_consistent());
EXPECT_EQ(record.id, 0);
- EXPECT_EQ(record.type, swc_record::kind::custom);
+ EXPECT_EQ(record.tag, soma_tag);
EXPECT_EQ(record.x, 1.);
EXPECT_EQ(record.y, 1.);
EXPECT_EQ(record.z, 1.);
@@ -94,7 +83,7 @@ TEST(swc_record, construction)
{
// check copy constructor
- swc_record record_orig(swc_record::kind::custom, 0, 1., 1., 1., 1., -1);
+ swc_record record_orig(soma_tag, 0, 1., 1., 1., 1., -1);
swc_record record(record_orig);
expect_record_equals(record_orig, record);
}
@@ -136,13 +125,6 @@ TEST(swc_parser, invalid_input_parse)
EXPECT_THROW(parse_swc_file(is), swc_error);
}
- {
- // Check invalid record type
- std::istringstream is(
- "1 10 14.566132 34.873772 7.857000 0.717830 -1\n");
- EXPECT_THROW(parse_swc_file(is), swc_error);
- }
-
{
// Non-contiguous numbering in branches is considered invalid
// 1
@@ -210,8 +192,9 @@ TEST(swc_parser, valid_input)
is >> record;
EXPECT_TRUE(is);
+ int soma_tag = 1;
swc_record record_expected(
- swc_record::kind::soma,
+ soma_tag,
0, 14.566132, 34.873772, 7.857000, 0.717830, -1);
expect_record_equals(record_expected, record);
}
@@ -227,8 +210,9 @@ TEST(swc_parser, valid_input)
is >> record;
EXPECT_TRUE(is);
+ int soma_tag = 1;
swc_record record_expected(
- swc_record::kind::soma,
+ soma_tag,
0, 14.566132, 34.873772, 7.857000, 0.717830, -1);
expect_record_equals(record_expected, record);
}
@@ -256,7 +240,7 @@ TEST(swc_parser, valid_input)
EXPECT_FALSE(is.fail());
EXPECT_FALSE(is.bad());
EXPECT_EQ(0, record.id); // zero-based indexing
- EXPECT_EQ(swc_record::kind::soma, record.type);
+ EXPECT_EQ(1, record.tag);
EXPECT_FLOAT_EQ(14.566132, record.x);
EXPECT_FLOAT_EQ(34.873772, record.y);
EXPECT_FLOAT_EQ( 7.857000, record.z);
@@ -265,11 +249,12 @@ TEST(swc_parser, valid_input)
}
{
+
// check valid input with a series of records
std::vector<swc_record> records_orig = {
- swc_record(swc_record::kind::soma, 0,
+ swc_record(1, 0,
14.566132, 34.873772, 7.857000, 0.717830, -1),
- swc_record(swc_record::kind::dendrite, 1,
+ swc_record(3, 1,
14.566132+1, 34.873772+1, 7.857000+1, 0.717830+1, -1)
};
@@ -399,7 +384,7 @@ TEST(swc_parser, raw)
std::stringstream is;
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
- is << "3 10 14.566132 34.873772 7.857000 0.717830 1\n";
+ is << "-3 2 14.566132 34.873772 7.857000 0.717830 1\n"; // invalid sample identifier -3
is << "4 2 14.566132 34.873772 7.857000 0.717830 1\n";
try {
@@ -421,6 +406,10 @@ TEST(swc_parser, raw)
}
}
+/***
+ * If you have proper unit tests for the steps from
+ * {sample_tree -> morphology -> cable_cell}
+ * then these tests can be simplified to check that {swc -> sample_tree} works
TEST(swc_io, cell_construction) {
using namespace arb;
@@ -557,3 +546,4 @@ TEST(swc_parser, from_file_ball_and_stick) {
expect_morph_eq(expected, bas_morph);
}
+*/