diff --git a/lmorpho/CMakeLists.txt b/lmorpho/CMakeLists.txt
index 1bb074d891494098c73d140e4b63538e98364cc6..126b93f8f280d888ddb4575a48dcc7a37d7ea316 100644
--- a/lmorpho/CMakeLists.txt
+++ b/lmorpho/CMakeLists.txt
@@ -1,4 +1,4 @@
-add_executable(lmorpho lmorpho.cpp lsystem.cpp lsys_models.cpp morphology.cpp morphio.cpp)
+add_executable(lmorpho lmorpho.cpp lsystem.cpp lsys_models.cpp morphio.cpp)
target_link_libraries(lmorpho LINK_PUBLIC nestmc)
target_link_libraries(lmorpho LINK_PUBLIC ${EXTERNAL_LIBRARIES})
diff --git a/lmorpho/lsystem.cpp b/lmorpho/lsystem.cpp
index 4cf824f9994833f56232ff6fa0df358ab5b040c6..f88ba290190d247ae2eb3e614f5b55f1ff7c936b 100644
--- a/lmorpho/lsystem.cpp
+++ b/lmorpho/lsystem.cpp
@@ -5,12 +5,16 @@
#include <vector>
#include <math.hpp>
+#include <morphology.hpp>
#include "lsystem.hpp"
-#include "morphology.hpp"
using namespace nest::mc::math;
+using nest::mc::section_geometry;
+using nest::mc::section_point;
+using nest::mc::morphology;
+
// L-system implementation.
// Random distribution used in the specification of L-system parameters.
@@ -263,9 +267,9 @@ grow_result grow(section_tip tip, const lsys_sampler& S, Gen &g) {
}
}
-morphology generate_morphology(const lsys_param& P, lsys_generator &g) {
+nest::mc::morphology generate_morphology(const lsys_param& P, lsys_generator &g) {
constexpr quaternion xaxis = {0, 1, 0, 0};
- morphology morph;
+ nest::mc::morphology morph;
lsys_sampler S(P);
double soma_radius = 0.5*S.diam_soma(g);
diff --git a/lmorpho/lsystem.hpp b/lmorpho/lsystem.hpp
index f3ce4deb3ca578df970153800da40c4cdd436acc..6568e279d5e859fece0a30e047e1e8d7be566f1a 100644
--- a/lmorpho/lsystem.hpp
+++ b/lmorpho/lsystem.hpp
@@ -2,13 +2,13 @@
#include <random>
-#include "morphology.hpp"
+#include <morphology.hpp>
struct lsys_param;
using lsys_generator = std::minstd_rand;
-morphology generate_morphology(const lsys_param& P, lsys_generator& g);
+nest::mc::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.
diff --git a/lmorpho/morphio.cpp b/lmorpho/morphio.cpp
index e6259be54a8d68b12f13f069b337057836c29f76..4c50fb6b77178e3466ae2937dd40b13b18eaaccd 100644
--- a/lmorpho/morphio.cpp
+++ b/lmorpho/morphio.cpp
@@ -4,13 +4,14 @@
#include <string>
#include <vector>
+#include <morphology.hpp>
#include <swcio.hpp>
#include "morphio.hpp"
using nest::mc::io::swc_record;
-std::vector<swc_record> as_swc(const morphology& morph);
+std::vector<swc_record> as_swc(const nest::mc::morphology& morph);
// printf wrappers.
@@ -78,7 +79,8 @@ void multi_file::open(unsigned n) {
using nest::mc::io::swc_record;
-std::vector<swc_record> as_swc(const morphology& morph) {
+// TODO: Move this functionality to nestmc library.
+std::vector<swc_record> as_swc(const nest::mc::morphology& morph) {
using kind = swc_record::kind;
std::map<int, int> parent_end_id;
std::vector<swc_record> swc;
@@ -127,7 +129,7 @@ std::vector<swc_record> as_swc(const morphology& morph) {
// SWC emitter implementation.
-void swc_emitter::operator()(unsigned index, const morphology& m) {
+void swc_emitter::operator()(unsigned index, const nest::mc::morphology& m) {
file_.open(index);
auto& stream = file_.stream();
@@ -138,7 +140,7 @@ void swc_emitter::operator()(unsigned index, const morphology& m) {
// pvector emitter implementation.
-std::vector<int> as_pvector(const morphology& morph, unsigned offset) {
+std::vector<int> as_pvector(const nest::mc::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
@@ -168,7 +170,7 @@ std::vector<int> as_pvector(const morphology& morph, unsigned offset) {
return pvec;
}
-void pvector_emitter::operator()(unsigned index, const morphology& m) {
+void pvector_emitter::operator()(unsigned index, const nest::mc::morphology& m) {
auto pvec = as_pvector(m, offset_);
if (coalesce_) offset_ += pvec.size();
diff --git a/lmorpho/morphio.hpp b/lmorpho/morphio.hpp
index 37ec56ae69e06caed6b9bb96078a94f41f20df8a..87a80db2f85488d3c9196fa6b4320c3860827658 100644
--- a/lmorpho/morphio.hpp
+++ b/lmorpho/morphio.hpp
@@ -4,7 +4,7 @@
#include <iostream>
#include <vector>
-#include "morphology.hpp"
+#include <morphology.hpp>
// Manage access to a single file, std::cout, or an indexed
// sequence of files.
@@ -59,7 +59,7 @@ public:
swc_emitter(swc_emitter&&) = default;
// write `index`th morphology as SWC.
- void operator()(unsigned index, const morphology& m);
+ void operator()(unsigned index, const nest::mc::morphology& m);
void close() { file_.close(); }
~swc_emitter() { close(); }
@@ -85,7 +85,7 @@ public:
pvector_emitter(pvector_emitter&&) = default;
// write pvector for `index`th morphology.
- void operator()(unsigned index, const morphology& m);
+ void operator()(unsigned index, const nest::mc::morphology& m);
void close() { file_.close(); }
~pvector_emitter() { close(); }
diff --git a/lmorpho/morphology.cpp b/lmorpho/morphology.cpp
deleted file mode 100644
index c06c76d873a41c4a97f4f9be8f28101d897c0edc..0000000000000000000000000000000000000000
--- a/lmorpho/morphology.cpp
+++ /dev/null
@@ -1,65 +0,0 @@
-#include <cmath>
-#include <vector>
-
-#include <math.hpp>
-
-#include "morphology.hpp"
-
-using nest::mc::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;
-}
diff --git a/lmorpho/morphology.hpp b/lmorpho/morphology.hpp
deleted file mode 100644
index d1adaa6fdf6a3c0966cb827a0f4894d0db8f3b51..0000000000000000000000000000000000000000
--- a/lmorpho/morphology.hpp
+++ /dev/null
@@ -1,31 +0,0 @@
-#pragma once
-
-// Representation of 3-d embedded cable morphology.
-
-#include <vector>
-
-struct section_point {
- double x, y, z, r; // [µm], r is radius.
-};
-
-struct section_geometry {
- unsigned id;
- unsigned parent_id;
- bool terminal;
- std::vector<section_point> points;
- double length; // µm
-
- // re-discretize the section into ceil(length/dx) segments.
- void segment(double dx);
-};
-
-struct morphology {
- section_point soma; // origin + spherical radius
- std::vector<section_geometry> sections;
-
- // re-discretize all sections.
- void segment(double dx) {
- for (auto& s: sections) s.segment(dx);
- }
-};
-
diff --git a/miniapp/miniapp_recipes.cpp b/miniapp/miniapp_recipes.cpp
index 742951dff6df83cb13c2f13af6e2d1ace2ae97eb..62d4b48cd1411e30b4589a1e8cab19c5f82f58f5 100644
--- a/miniapp/miniapp_recipes.cpp
+++ b/miniapp/miniapp_recipes.cpp
@@ -4,6 +4,7 @@
#include <utility>
#include <cell.hpp>
+#include <morphology.hpp>
#include <util/debug.hpp>
#include "miniapp_recipes.hpp"
@@ -29,9 +30,9 @@ cell make_basic_cell(
// add dendrite of length 200 um and diameter 1 um with passive channel
std::vector<mc::cable_segment*> dendrites;
- dendrites.push_back(cell.add_cable(0, mc::segmentKind::dendrite, 0.5, 0.5, 200));
- dendrites.push_back(cell.add_cable(1, mc::segmentKind::dendrite, 0.5, 0.25,100));
- dendrites.push_back(cell.add_cable(1, mc::segmentKind::dendrite, 0.5, 0.25,100));
+ dendrites.push_back(cell.add_cable(0, section_kind::dendrite, 0.5, 0.5, 200));
+ dendrites.push_back(cell.add_cable(1, section_kind::dendrite, 0.5, 0.25,100));
+ dendrites.push_back(cell.add_cable(1, section_kind::dendrite, 0.5, 0.25,100));
for (auto d : dendrites) {
d->add_mechanism(mc::pas_parameters());
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ca4032a3a9b4d161de99d409a90c177da40760b9..b3faee92ef7df76c17935e6806ea7eb8dd2006cd 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -4,6 +4,7 @@ set(HEADERS
set(BASE_SOURCES
common_types_io.cpp
cell.cpp
+ morphology.cpp
parameter_list.cpp
profiling/profiler.cpp
swcio.cpp
diff --git a/src/cell.cpp b/src/cell.cpp
index 1828ff83ca6ad2dd6ce7900a18ce4061082c9c22..04766c84a4f80e8b20ca539fee1f0effc08beb20 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -1,4 +1,5 @@
#include <cell.hpp>
+#include <morphology.hpp>
#include <tree.hpp>
#include <util/debug.hpp>
@@ -235,5 +236,49 @@ bool cell_basic_equality(cell const& lhs, cell const& rhs)
return true;
}
+// Construct cell from flat morphology specification.
+
+cell make_cell(const morphology& morph, bool compartments_from_discretization) {
+ using point3d = cell::point_type;
+ cell newcell;
+
+ if (!morph) {
+ return newcell;
+ }
+
+ EXPECTS(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 std::invalid_argument("no support for complex somata");
+ break;
+ default: ;
+ }
+
+ std::vector<cell::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));
+ }
+
+ auto cable = newcell.add_cable(section.parent_id, kind, radii, points);
+ if (compartments_from_discretization) {
+ cable->as_cable()->set_compartments(radii.size()-1);
+ }
+ }
+
+ return newcell;
+}
+
} // namespace mc
} // namespace nest
diff --git a/src/cell.hpp b/src/cell.hpp
index 4b9052f2c2698ee6dd66192b98f5cd3c268d18ff..0c2bdd13a095e724062a430ad8ca24a8fdba1c7e 100644
--- a/src/cell.hpp
+++ b/src/cell.hpp
@@ -7,6 +7,7 @@
#include <common_types.hpp>
#include <cell_tree.hpp>
+#include <morphology.hpp>
#include <segment.hpp>
#include <stimulus.hpp>
#include <util/debug.hpp>
@@ -245,5 +246,11 @@ cable_segment* cell::add_cable(cell::index_type parent, Args&&... args)
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.
+cell make_cell(const morphology&, bool compartments_from_discretization=false);
+
} // namespace mc
} // namespace nest
diff --git a/src/morphology.cpp b/src/morphology.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..740d488386ea6ea407d43c2dd96d7ab26fedd512
--- /dev/null
+++ b/src/morphology.cpp
@@ -0,0 +1,135 @@
+#include <cmath>
+#include <vector>
+
+#include <math.hpp>
+
+#include <morphology.hpp>
+
+namespace nest {
+namespace mc {
+
+using ::nest::mc::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");
+ }
+ sections[section.parent_id].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(true, nsection);
+
+ 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) {
+ auto parent_index = parent_id-1;
+ terminal[parent_index] = false;
+ }
+ }
+
+ 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 mc
+} // namespace nest
diff --git a/src/morphology.hpp b/src/morphology.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e57278f28ed89db6249c8ca5b231df71f3f41d65
--- /dev/null
+++ b/src/morphology.hpp
@@ -0,0 +1,75 @@
+#pragma once
+
+// Representation of 3-d embedded cable morphology, independent of other
+// cell information.
+
+#include <stdexcept>
+#include <vector>
+
+namespace nest {
+namespace mc {
+
+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;
+
+ // 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(); }
+
+ // 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 mc
+} // namespace nest
diff --git a/src/segment.hpp b/src/segment.hpp
index dc518170485c03b4616462e18132f32782026eb6..53a486c2742b06375687a0eea22dc969dbfc3b71 100644
--- a/src/segment.hpp
+++ b/src/segment.hpp
@@ -7,6 +7,7 @@
#include <common_types.hpp>
#include <compartment.hpp>
#include <math.hpp>
+#include <morphology.hpp>
#include <parameter_list.hpp>
#include <point.hpp>
#include <util/make_unique.hpp>
@@ -21,13 +22,6 @@ struct segment_properties {
T cm = 0.01; // capacitance [F/m^2] : 10 nF/mm^2 = 0.01 F/m^2
};
-enum class segmentKind {
- soma,
- dendrite,
- axon,
- none
-};
-
// forward declarations of segment specializations
class soma_segment;
class cable_segment;
@@ -42,23 +36,23 @@ public:
// (Yet more motivation for a separate morphology description class!)
virtual std::unique_ptr<segment> clone() const = 0;
- segmentKind kind() const {
+ section_kind kind() const {
return kind_;
}
bool is_soma() const
{
- return kind_==segmentKind::soma;
+ return kind_==section_kind::soma;
}
bool is_dendrite() const
{
- return kind_==segmentKind::dendrite;
+ return kind_==section_kind::dendrite;
}
bool is_axon() const
{
- return kind_==segmentKind::axon;
+ return kind_==section_kind::axon;
}
virtual size_type num_compartments() const = 0;
@@ -152,7 +146,7 @@ public:
}
protected:
- segmentKind kind_;
+ section_kind kind_;
std::vector<parameter_list> mechanisms_;
};
@@ -165,7 +159,7 @@ public:
placeholder_segment()
{
- kind_ = segmentKind::none;
+ kind_ = section_kind::none;
}
std::unique_ptr<segment> clone() const override {
@@ -209,7 +203,7 @@ public:
soma_segment(value_type r):
radius_{r}
{
- kind_ = segmentKind::soma;
+ kind_ = section_kind::soma;
mechanisms_.push_back(membrane_parameters());
}
@@ -281,12 +275,12 @@ public:
// constructors for a cable with no location information
cable_segment(
- segmentKind k,
+ section_kind k,
std::vector<value_type> r,
std::vector<value_type> lens
) {
kind_ = k;
- assert(k==segmentKind::dendrite || k==segmentKind::axon);
+ assert(k==section_kind::dendrite || k==section_kind::axon);
radii_ = std::move(r);
lengths_ = std::move(lens);
@@ -296,7 +290,7 @@ public:
}
cable_segment(
- segmentKind k,
+ section_kind k,
value_type r1,
value_type r2,
value_type len
@@ -307,12 +301,12 @@ public:
// constructor that lets the user describe the cable as a
// seriew of radii and locations
cable_segment(
- segmentKind k,
+ section_kind k,
std::vector<value_type> r,
std::vector<point_type> p
) {
kind_ = k;
- assert(k==segmentKind::dendrite || k==segmentKind::axon);
+ assert(k==section_kind::dendrite || k==section_kind::axon);
radii_ = std::move(r);
locations_ = std::move(p);
@@ -326,7 +320,7 @@ public:
// of just the end points of the cable
// i.e. describing the cable as a single frustrum
cable_segment(
- segmentKind k,
+ section_kind k,
value_type r1,
value_type r2,
point_type const& p1,
@@ -463,7 +457,7 @@ using segment_ptr = std::unique_ptr<segment>;
/// Helper for constructing segments in a segment_ptr unique pointer wrapper.
/// Forwards the supplied arguments to construct a segment of type SegmentType.
-/// e.g. auto my_cable = make_segment<cable>(segmentKind::dendrite, ... );
+/// e.g. auto my_cable = make_segment<cable>(section_kind::dendrite, ... );
template <typename SegmentType, typename... Args>
segment_ptr make_segment(Args&&... args) {
return segment_ptr(new SegmentType(std::forward<Args>(args)...));
diff --git a/src/swcio.cpp b/src/swcio.cpp
index 1c5ee0bb0fb1345d627e68f6596cde8459145035..023bbf2880c44f97cee788356c616fc4525d7302 100644
--- a/src/swcio.cpp
+++ b/src/swcio.cpp
@@ -6,6 +6,8 @@
#include <unordered_set>
#include <algorithms.hpp>
+#include <cell.hpp>
+#include <morphology.hpp>
#include <point.hpp>
#include <swcio.hpp>
#include <util/debug.hpp>
@@ -14,314 +16,143 @@ namespace nest {
namespace mc {
namespace io {
-//
// swc_record implementation
-//
-void swc_record::renumber(id_type new_id, std::map<id_type, id_type>& idmap)
-{
- auto old_id = id_;
- id_ = new_id;
- // Obtain parent_id from the map
- auto new_parent_id = idmap.find(parent_id_);
- if (new_parent_id != idmap.end()) {
- parent_id_ = new_parent_id->second;
- }
- check_consistency();
- idmap.insert(std::make_pair(old_id, new_id));
-}
+// 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);
-void swc_record::check_consistency() const
-{
- // Check record type as well; enum's do not offer complete type safety,
- // since you can cast anything that fits to its underlying type
- if (static_cast<int>(type_) < 0 ||
- static_cast<int>(type_) > static_cast<int>(kind::custom)) {
- throw std::invalid_argument("unknown record type");
+ if (static_cast<int>(r.type) < 0 || static_cast<int>(r.type) > max_type) {
+ return "unknown record type";
}
- if (id_ < 0) {
- throw std::invalid_argument("negative ids not allowed");
+ if (r.id < 0) {
+ return "negative ids not allowed";
}
- if (parent_id_ < -1) {
- throw std::invalid_argument("parent_id < -1 not allowed");
+ if (r.parent_id < -1) {
+ return "parent_id < -1 not allowed";
}
- if (parent_id_ >= id_) {
- throw std::invalid_argument("parent_id >= id is not allowed");
+ if (r.parent_id >= r.id) {
+ return "parent_id >= id is not allowed";
}
- if (r_ < 0) {
- throw std::invalid_argument("negative radii are not allowed");
+ if (r.r < 0) {
+ return "negative radii are not allowed";
}
-}
-std::istream& operator>>(std::istream& is, swc_record& record)
-{
- swc_parser parser;
- parser.parse_record(is, record);
- return is;
+ return nullptr;
}
-std::ostream& operator<<(std::ostream& os, const swc_record& record)
-{
- // output in one-based indexing
- os << record.id_+1 << " "
- << static_cast<int>(record.type_) << " "
- << std::setprecision(7) << record.x_ << " "
- << std::setprecision(7) << record.y_ << " "
- << std::setprecision(7) << record.z_ << " "
- << std::setprecision(7) << record.r_ << " "
- << ((record.parent_id_ == -1) ? record.parent_id_ : record.parent_id_+1);
-
- return os;
-}
-
-
-//
-// Utility functions
-//
-
-std::string::size_type find_first_non_whitespace(const std::string& str)
-{
- return str.find_first_not_of(" \f\n\r\t\v");
-}
-
-bool starts_with(const std::string& str, const std::string& prefix)
-{
- // ignore leading whitespace
- auto pos = find_first_non_whitespace(str);
- if (pos == std::string::npos) {
- return false;
- }
-
- return str.find(prefix, pos) == pos;
-}
-
-bool is_space(const std::string& str)
-{
- return find_first_non_whitespace(str) == std::string::npos;
+bool swc_record::is_consistent() const {
+ return swc_record_error(*this)==nullptr;
}
-void check_parse_status(const std::istream& is, const swc_parser& parser)
-{
- if (is.fail()) {
- // If we try to read past the eof; fail bit will also be set
- throw swc_parse_error("could not parse value", parser.lineno());
+void swc_record::assert_consistent() const {
+ const char* error = swc_record_error(*this);
+ if (error) {
+ throw swc_error(error);
}
}
-nest::mc::segmentKind convert_kind(const swc_record::kind& kind)
-{
- switch (kind) {
- case swc_record::kind::soma:
- return segmentKind::soma;
- case swc_record::kind::dendrite:
- return segmentKind::dendrite;
- case swc_record::kind::axon:
- return segmentKind::axon;
- default:
- throw swc_parse_error("no known conversion for swc record type", 0);
+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);
+ if (is) {
+ // Convert to zero-based, leaving parent_id as-is if -1
+ --r.id;
+ if (r.parent_id>=0) {
+ --r.parent_id;
+ }
+ record = r;
+ return true;
}
+ return false;
}
-
-template<typename T>
-T parse_value_strict(std::istream& is, const swc_parser& parser)
-{
- T val;
- check_parse_status(is >> val, parser);
-
- // everything's fine
- return val;
+bool is_comment(const std::string& line) {
+ auto pos = line.find_first_not_of(" \f\n\r\t\v");
+ return pos==std::string::npos || line[pos]=='#';
}
-// specialize parsing for record types
-template<>
-swc_record::kind parse_value_strict(std::istream& is, const swc_parser& parser)
-{
- swc_record::id_type val;
- check_parse_status(is >> val, parser);
-
- // Let swc_record's constructor check for the type validity
- return static_cast<swc_record::kind>(val);
-}
+std::istream& operator>>(std::istream& is, swc_record& record) {
+ std::string line;
-//
-// swc_parser implementation
-//
-
-std::istream& swc_parser::parse_record(std::istream& is, swc_record& record)
-{
- while (!is.eof() && !is.bad()) {
- // consume empty and comment lines first
- std::getline(is, linebuff_);
-
- ++lineno_;
- if (!is_space(linebuff_) &&
- !starts_with(linebuff_, comment_prefix_)) {
+ while (is) {
+ std::getline(is, line);
+ if (!is) {
break;
}
- }
- if (is.bad()) {
- // let the caller check for such events
- return is;
- }
-
- if (is.eof() &&
- (is_space(linebuff_) || starts_with(linebuff_, comment_prefix_))) {
- // last line is either empty or a comment; don't parse anything
- return is;
- }
-
- if (is.fail()) {
- throw swc_parse_error("too long line detected", lineno_);
- }
+ if (is_comment(line)) {
+ continue;
+ }
- std::istringstream line(linebuff_);
- try {
- record = parse_record(line);
- } catch (std::invalid_argument& e) {
- // Rethrow as a parse error
- throw swc_parse_error(e.what(), lineno_);
+ bool ok = parse_record(line, record);
+ if (!ok) {
+ is.setstate(std::ios::failbit);
+ }
+ break;
}
-
return is;
}
-swc_record swc_parser::parse_record(std::istringstream& is)
-{
- auto id = parse_value_strict<int>(is, *this);
- auto type = parse_value_strict<swc_record::kind>(is, *this);
- auto x = parse_value_strict<swc_record::coord_type>(is, *this);
- auto y = parse_value_strict<swc_record::coord_type>(is, *this);
- auto z = parse_value_strict<swc_record::coord_type>(is, *this);
- auto r = parse_value_strict<swc_record::coord_type>(is, *this);
- auto parent_id = parse_value_strict<swc_record::id_type>(is, *this);
-
- // Convert to zero-based, leaving parent_id as-is if -1
- if (parent_id != -1) {
- parent_id--;
- }
+std::ostream& operator<<(std::ostream& os, const swc_record& record) {
+ // output in one-based indexing
+ os << record.id+1 << " "
+ << static_cast<int>(record.type) << " "
+ << std::setprecision(7) << record.x << " "
+ << std::setprecision(7) << record.y << " "
+ << std::setprecision(7) << record.z << " "
+ << std::setprecision(7) << record.r << " "
+ << ((record.parent_id == -1) ? record.parent_id : record.parent_id+1);
- return swc_record(type, id-1, x, y, z, r, parent_id);
+ return os;
}
+std::vector<swc_record> parse_swc_file(std::istream& is) {
+ std::vector<swc_record> records;
+ std::size_t line_number = 1;
+ std::string line;
-swc_record_range_clean::swc_record_range_clean(std::istream& is)
-{
- std::unordered_set<swc_record::id_type> ids;
-
- std::size_t num_trees = 0;
- swc_record::id_type last_id = -1;
- bool needsort = false;
-
- swc_record curr_record;
- for (auto r : swc_get_records<swc_io_raw>(is)) {
- if (r.parent() == -1 && ++num_trees > 1) {
- // only a single tree is allowed
- break;
- }
-
- auto inserted = ids.insert(r.id());
- if (inserted.second) {
- // not a duplicate; insert record
- records_.push_back(r);
- if (!needsort && r.id() < last_id) {
- needsort = true;
+ try {
+ while (is) {
+ std::getline(is, line);
+ if (is_comment(line)) {
+ continue;
}
- last_id = r.id();
- }
- }
-
- if (needsort) {
- std::sort(records_.begin(), records_.end());
- }
-
- // Renumber records if necessary
- std::map<swc_record::id_type, swc_record::id_type> idmap;
- swc_record::id_type next_id = 0;
- for (auto& r : records_) {
- if (r.id() != next_id) {
- r.renumber(next_id, idmap);
+ swc_record record;
+ bool ok = parse_record(line, record);
+ if (!ok) {
+ is.setstate(std::ios::failbit);
+ }
+ else {
+ record.assert_consistent();
+ records.push_back(std::move(record));
+ }
+ ++line_number;
}
-
- ++next_id;
}
-
- // Reject if branches are not contiguously numbered
- std::vector<swc_record::id_type> parent_list = { 0 };
- for (std::size_t i = 1; i < records_.size(); ++i) {
- parent_list.push_back(records_[i].parent());
+ catch (swc_error& e) {
+ e.line_number = line_number; // rethrow with saved line number
+ throw e;
}
- if (!nest::mc::algorithms::has_contiguous_compartments(parent_list)) {
- throw swc_parse_error("branches are not contiguously numbered", 0);
- }
-}
-
-cell swc_read_cell(std::istream& is)
-{
- using namespace nest::mc;
-
- cell newcell;
- std::vector<swc_record::id_type> parent_index;
- std::vector<swc_record> swc_records;
- for (const auto& r : swc_get_records<swc_io_clean>(is)) {
- swc_records.push_back(r);
- parent_index.push_back(r.parent());
- }
-
- if (parent_index.empty()) {
- return newcell;
- }
-
- // The parent of soma must be 0, while in SWC files is -1
- parent_index[0] = 0;
- auto branch_index = algorithms::branches(parent_index);
- auto new_parent_index = algorithms::make_parent_index(parent_index,
- branch_index);
-
- // sanity check
- EXPECTS(new_parent_index.size() == branch_index.size() - 1);
-
- // Add the soma first; then the segments
- newcell.add_soma(swc_records[0].radius(), swc_records[0].coord());
- for (std::size_t i = 1; i < new_parent_index.size(); ++i) {
- auto b_start = std::next(swc_records.begin(), branch_index[i]);
- auto b_end = std::next(swc_records.begin(), branch_index[i+1]);
-
- std::vector<swc_record::coord_type> radii;
- std::vector<nest::mc::point<swc_record::coord_type>> points;
- if (new_parent_index[i] != 0) {
- // include the parent of current record if not branching from soma
- auto p = parent_index[branch_index[i]];
- radii.push_back(swc_records[p].radius());
- points.push_back(swc_records[p].coord());
- }
-
- // extract the radii and the points
- std::for_each(b_start, b_end,
- [&radii](const swc_record& r) {
- radii.push_back(r.radius());
- });
-
- std::for_each(b_start, b_end,
- [&points](const swc_record& r) {
- points.push_back(r.coord());
- });
-
- // add the new cable
- newcell.add_cable(new_parent_index[i],
- convert_kind(b_start->type()), radii, points);
+ if (!is.eof()) {
+ // parse error, so throw exception
+ throw swc_error("SWC parse error", line_number);
}
- return newcell;
+ swc_canonicalize_sequence(records);
+ return records;
}
} // namespace io
diff --git a/src/swcio.hpp b/src/swcio.hpp
index 6530cba35190049d58be20dcebdcf9d81570053e..019b26ff384f6ac627cef07576054a09004d234b 100644
--- a/src/swcio.hpp
+++ b/src/swcio.hpp
@@ -3,11 +3,15 @@
#include <exception>
#include <iostream>
#include <iterator>
+#include <map>
#include <string>
+#include <unordered_set>
#include <vector>
-#include "cell.hpp"
-#include "point.hpp"
+#include <algorithms.hpp>
+#include <morphology.hpp>
+#include <point.hpp>
+#include <util/debug.hpp>
namespace nest {
namespace mc {
@@ -30,384 +34,213 @@ public:
custom
};
+ kind type = kind::undefined; // 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
+
// swc records assume zero-based indexing; root's parent remains -1
swc_record(swc_record::kind type, 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)
- {
- check_consistency();
- }
-
- swc_record()
- : type_(swc_record::kind::undefined)
- , id_(0)
- , x_(0)
- , y_(0)
- , z_(0)
- , r_(0)
- , parent_id_(-1)
- { }
+ int parent_id):
+ type(type), id(id), x(x), y(y), z(z), r(r), parent_id(parent_id)
+ {}
+ swc_record() = default;
swc_record(const swc_record& other) = default;
swc_record& operator=(const swc_record& other) = default;
- bool strict_equals(const swc_record& other) const
- {
- return id_ == other.id_ &&
- x_ == other.x_ &&
- y_ == other.y_ &&
- z_ == other.z_ &&
- r_ == other.r_ &&
- parent_id_ == other.parent_id_;
- }
-
- // Equality and comparison operators
- friend bool operator==(const swc_record& lhs,
- const swc_record& rhs)
- {
- return lhs.id_ == rhs.id_;
- }
-
- friend bool operator<(const swc_record& lhs,
- const swc_record& rhs)
- {
- return lhs.id_ < rhs.id_;
- }
-
- friend bool operator<=(const swc_record& lhs,
- const swc_record& rhs)
- {
- return (lhs < rhs) || (lhs == rhs);
+ bool operator==(const swc_record& other) const {
+ return id == other.id &&
+ x == other.x &&
+ y == other.y &&
+ z == other.z &&
+ r == other.r &&
+ parent_id == other.parent_id;
}
- friend bool operator!=(const swc_record& lhs,
- const swc_record& rhs)
- {
+ friend bool operator!=(const swc_record& lhs, const swc_record& rhs) {
return !(lhs == rhs);
}
- friend bool operator>(const swc_record& lhs,
- const swc_record& rhs)
- {
- return !(lhs < rhs) && (lhs != rhs);
- }
-
- friend bool operator>=(const swc_record& lhs,
- const swc_record& rhs)
- {
- return !(lhs < rhs);
- }
-
friend std::ostream& operator<<(std::ostream& os, const swc_record& record);
- kind type() const
- {
- return type_;
- }
-
- id_type id() const
- {
- return id_;
- }
-
- id_type parent() const
- {
- return parent_id_;
+ coord_type diameter() const {
+ return 2*r;
}
- coord_type x() const
- {
- return x_;
+ nest::mc::point<coord_type> coord() const {
+ return nest::mc::point<coord_type>(x, y, z);
}
- coord_type y() const
- {
- return y_;
+ nest::mc::section_point as_section_point() const {
+ return nest::mc::section_point{x, y, z, r};
}
- coord_type z() const
- {
- return z_;
- }
+ // validity checks
+ bool is_consistent() const;
+ void assert_consistent() const; // throw swc_error if inconsistent.
+};
- coord_type radius() const
- {
- return r_;
- }
- coord_type diameter() const
- {
- return 2*r_;
- }
+class swc_error: public std::runtime_error {
+public:
+ explicit swc_error(const char* msg, std::size_t lineno = 0):
+ std::runtime_error(msg), line_number(lineno)
+ {}
- nest::mc::point<coord_type> coord() const
- {
- return nest::mc::point<coord_type>(x_, y_, z_);
- }
+ explicit swc_error(const std::string& msg, std::size_t lineno = 0):
+ std::runtime_error(msg), line_number(lineno)
+ {}
- void renumber(id_type new_id, std::map<id_type, id_type>& idmap);
+ std::size_t line_number;
+};
-private:
- void check_consistency() const;
+// Parse one record, skipping comments and blank lines.
+std::istream& operator>>(std::istream& is, swc_record& record);
- kind type_; // record type
- id_type id_; // record id
- coord_type x_, y_, z_; // record coordinates
- coord_type r_; // record radius
- id_type parent_id_; // record parent's id
-};
+// Parse and canonicalize an EOF-terminated sequence of records.
+// Throw on parsing failure.
+std::vector<swc_record> parse_swc_file(std::istream& is);
+// Convert a canonical (see below) sequence of SWC records to a morphology object.
+template <typename RandomAccessSequence>
+morphology swc_as_morphology(const RandomAccessSequence& swc_records) {
+ morphology morph;
-class swc_parse_error : public std::runtime_error
-{
-public:
- explicit swc_parse_error(const char* msg, std::size_t lineno)
- : std::runtime_error(msg)
- , lineno_(lineno)
- { }
-
- explicit swc_parse_error(const std::string& msg, std::size_t lineno)
- : std::runtime_error(msg)
- , lineno_(lineno)
- { }
-
- std::size_t lineno() const
- {
- return lineno_;
+ std::vector<swc_record::id_type> swc_parent_index;
+ for (const auto& r: swc_records) {
+ swc_parent_index.push_back(r.parent_id);
}
-private:
- std::size_t lineno_;
-};
-
-class swc_parser {
-public:
- swc_parser(const std::string& delim, std::string comment_prefix)
- : delim_(delim)
- , comment_prefix_(comment_prefix)
- , lineno_(0)
- { }
-
- swc_parser()
- : delim_(" ")
- , comment_prefix_("#")
- , lineno_(0)
- { }
-
- std::size_t lineno() const
- {
- return lineno_;
+ if (swc_parent_index.empty()) {
+ return morph;
}
- std::istream& parse_record(std::istream& is, swc_record& record);
+ // 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::make_parent_index(swc_parent_index, branch_index);
-private:
- // Read the record from a string stream; will be treated like a single line
- swc_record parse_record(std::istringstream& is);
+ // sanity check
+ EXPECTS(parent_branch_index.size() == branch_index.size() - 1);
- std::string delim_;
- std::string comment_prefix_;
- std::string linebuff_;
- std::size_t lineno_;
-};
+ // Add the soma first; then the segments
+ const auto& soma = swc_records[0];
+ morph.soma = { soma.x, soma.y, soma.z, soma.r };
+ auto n_branches = parent_branch_index.size();
+ for (std::size_t i = 1; i < n_branches; ++i) {
+ auto b_start = std::next(swc_records.begin(), branch_index[i]);
+ auto b_end = std::next(swc_records.begin(), branch_index[i+1]);
-std::istream& operator>>(std::istream& is, swc_record& record);
+ section_geometry section;
+ section.id = i;
+ section.parent_id = parent_branch_index[i];
-class swc_record_stream_iterator :
- public std::iterator<std::forward_iterator_tag, swc_record> {
-public:
- struct eof_tag { };
-
- swc_record_stream_iterator(std::istream& is)
- : is_(is)
- , eof_(false)
- {
- is_.clear();
- is_.seekg(std::ios_base::beg);
- read_next_record();
- }
+ if (section.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]]];
- swc_record_stream_iterator(std::istream& is, eof_tag)
- : is_(is)
- , eof_(true)
- { }
-
- swc_record_stream_iterator(const swc_record_stream_iterator& other)
- : is_(other.is_)
- , parser_(other.parser_)
- , curr_record_(other.curr_record_)
- , eof_(other.eof_)
- { }
-
- swc_record_stream_iterator& operator++()
- {
- if (eof_) {
- throw std::out_of_range("attempt to read past eof");
+ section_point point{parent_record.x, parent_record.y, parent_record.z, parent_record.r};
+ section.points.push_back(point);
}
- read_next_record();
- return *this;
- }
-
- swc_record_stream_iterator operator++(int)
- {
- swc_record_stream_iterator ret(*this);
- operator++();
- return ret;
- }
-
- value_type operator*() const
- {
- if (eof_) {
- throw std::out_of_range("attempt to read past eof");
+ for (auto b = b_start; b!=b_end; ++b) {
+ section_point point{b->x, b->y, b->z, b->r};
+ section.points.push_back(point);
+
+ switch (b->type) {
+ case swc_record::kind::axon:
+ section.kind = section_kind::axon;
+ break;
+ case swc_record::kind::dendrite:
+ case swc_record::kind::apical_dendrite:
+ section.kind = section_kind::dendrite;
+ break;
+ case swc_record::kind::soma:
+ section.kind = section_kind::soma;
+ break;
+ default: ; // stick with what we have
+ }
}
- return curr_record_;
+ morph.sections.push_back(section);
}
- bool operator==(const swc_record_stream_iterator& other) const
- {
- if (eof_ && other.eof_) {
- return true;
- } else {
- return curr_record_.strict_equals(other.curr_record_);
- }
- }
+ EXPECTS(morph.check_valid());
+ return morph;
+}
- bool operator!=(const swc_record_stream_iterator& other) const
- {
- return !(*this == other);
- }
+// Given a 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
+// are ordered with repect to parent indices.
+template <typename RandomAccessSequence>
+void swc_canonicalize_sequence(RandomAccessSequence& swc_records) {
+ std::unordered_set<swc_record::id_type> ids;
- friend std::ostream& operator<<(std::ostream& os,
- const swc_record_stream_iterator& iter)
- {
- os << "{ is_.tellg(): " << iter.is_.tellg() << ", "
- << "curr_record_: " << iter.curr_record_ << ", "
- << "eof_: " << iter.eof_ << "}";
+ std::size_t num_trees = 0;
+ swc_record::id_type last_id = -1;
+ bool needsort = false;
- return os;
- }
+ for (const auto& r: swc_records) {
+ r.assert_consistent();
-private:
- void read_next_record()
- {
- parser_.parse_record(is_, curr_record_);
- if (is_.eof()) {
- eof_ = true;
+ if (r.parent_id == -1 && ++num_trees > 1) {
+ // only a single tree is allowed
+ throw swc_error("multiple trees found in SWC record sequence");
+ }
+ if (ids.count(r.id)) {
+ throw swc_error("records with duplicated ids in SWC record sequence");
}
- }
-
- std::istream& is_;
- swc_parser parser_;
- swc_record curr_record_;
-
- // indicator of eof; we need a way to define an end() iterator without
- // seeking to the end of file
- bool eof_;
-};
-
-class swc_record_range_raw {
-public:
- using value_type = swc_record;
- using reference = value_type&;
- using const_reference = const value_type&;
- using iterator = swc_record_stream_iterator;
- using const_iterator = const swc_record_stream_iterator;
-
- swc_record_range_raw(std::istream& is)
- : is_(is)
- { }
-
- iterator begin() const
- {
- return swc_record_stream_iterator(is_);
- }
+ if (!needsort && r.id < last_id) {
+ needsort = true;
+ }
- iterator end() const
- {
- iterator::eof_tag eof;
- return swc_record_stream_iterator(is_, eof);
+ last_id = r.id;
+ ids.insert(r.id);
}
- bool empty() const
- {
- return begin() == end();
+ if (needsort) {
+ std::sort(std::begin(swc_records), std::end(swc_records),
+ [](const swc_record& a, const swc_record& b) { return a.id<b.id; });
}
-private:
- std::istream& is_;
-};
-
-//
-// Reads records from an input stream until an eof is encountered and returns a
-// cleaned sequence of swc records.
-//
-// For more information check here:
-// https://github.com/eth-cscs/cell_algorithms/wiki/SWC-file-parsing
-//
-
-class swc_record_range_clean {
-public:
- using value_type = swc_record;
- using reference = value_type&;
- using const_referene = const value_type&;
- using iterator = std::vector<swc_record>::iterator;
- using const_iterator = std::vector<swc_record>::const_iterator;
-
- swc_record_range_clean(std::istream& is);
+ // Renumber records if necessary
+ std::map<swc_record::id_type, swc_record::id_type> idmap;
+ swc_record::id_type next_id = 0;
+ for (auto& r: swc_records) {
+ if (r.id != next_id) {
+ auto old_id = r.id;
+ r.id = next_id;
- iterator begin()
- {
- return records_.begin();
- }
+ auto new_parent_id = idmap.find(r.parent_id);
+ if (new_parent_id != idmap.end()) {
+ r.parent_id = new_parent_id->second;
+ }
- iterator end()
- {
- return records_.end();
+ r.assert_consistent();
+ idmap.insert(std::make_pair(old_id, next_id));
+ }
+ ++next_id;
}
- std::size_t size()
- {
- return records_.size();
+ // Reject if branches are not contiguously numbered
+ std::vector<swc_record::id_type> parent_list = { 0 };
+ for (std::size_t i = 1; i < swc_records.size(); ++i) {
+ parent_list.push_back(swc_records[i].parent_id);
}
- bool empty() const
- {
- return records_.empty();
+ if (!nest::mc::algorithms::has_contiguous_compartments(parent_list)) {
+ throw swc_error("branches are not contiguously numbered", 0);
}
-
-private:
- std::vector<swc_record> records_;
-};
-
-struct swc_io_raw
-{
- using record_range_type = swc_record_range_raw;
-};
-
-struct swc_io_clean
-{
- using record_range_type = swc_record_range_clean;
-};
-
-template<typename T = swc_io_clean>
-typename T::record_range_type swc_get_records(std::istream& is)
-{
- return typename T::record_range_type(is);
}
-cell swc_read_cell(std::istream& is);
-
} // namespace io
} // namespace mc
} // namespace nest
diff --git a/tests/test_common_cells.hpp b/tests/test_common_cells.hpp
index 43601ee9956e87ec74c39df990ece68ec6bf6261..e5f1a49bc40147822ad9e22e18c26bb87d63ecbe 100644
--- a/tests/test_common_cells.hpp
+++ b/tests/test_common_cells.hpp
@@ -63,7 +63,7 @@ inline cell make_cell_ball_and_stick(bool with_stim = true) {
auto soma = c.add_soma(12.6157/2.0);
soma->add_mechanism(hh_parameters());
- c.add_cable(0, segmentKind::dendrite, 1.0/2, 1.0/2, 200.0);
+ c.add_cable(0, section_kind::dendrite, 1.0/2, 1.0/2, 200.0);
for (auto& seg: c.segments()) {
seg->mechanism("membrane").set("r_L", 100);
@@ -108,7 +108,7 @@ inline cell make_cell_ball_and_taper(bool with_stim = true) {
auto soma = c.add_soma(12.6157/2.0);
soma->add_mechanism(hh_parameters());
- c.add_cable(0, segmentKind::dendrite, 1.0/2, 0.4/2, 200.0);
+ c.add_cable(0, section_kind::dendrite, 1.0/2, 0.4/2, 200.0);
for (auto& seg: c.segments()) {
seg->mechanism("membrane").set("r_L", 100);
@@ -162,7 +162,7 @@ inline cell make_cell_ball_and_squiggle(bool with_stim = true) {
};
auto dendrite =
- make_segment<cable_segment>(segmentKind::dendrite, radii, points);
+ make_segment<cable_segment>(section_kind::dendrite, radii, points);
c.add_cable(0, std::move(dendrite));
for (auto& seg: c.segments()) {
@@ -209,9 +209,9 @@ inline cell make_cell_ball_and_3stick(bool with_stim = true) {
auto soma = c.add_soma(12.6157/2.0);
soma->add_mechanism(hh_parameters());
- c.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 100);
- c.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100);
- c.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100);
+ c.add_cable(0, section_kind::dendrite, 0.5, 0.5, 100);
+ c.add_cable(1, section_kind::dendrite, 0.5, 0.5, 100);
+ c.add_cable(1, section_kind::dendrite, 0.5, 0.5, 100);
for (auto& seg: c.segments()) {
seg->mechanism("membrane").set("r_L", 100);
@@ -260,7 +260,7 @@ inline cell make_cell_simple_cable(bool with_stim = true) {
cell c;
c.add_soma(0);
- c.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 1000);
+ c.add_cable(0, section_kind::dendrite, 0.5, 0.5, 1000);
double r_L = 100;
double c_m = 0.01;
diff --git a/tests/unit/test_cell.cpp b/tests/unit/test_cell.cpp
index 143128fd43865df24daca89f9ea4a632e7181c7d..e08c0650b62e3e20383d683bb6cac56ea3307325 100644
--- a/tests/unit/test_cell.cpp
+++ b/tests/unit/test_cell.cpp
@@ -54,7 +54,7 @@ TEST(cell_type, add_segment)
auto seg =
make_segment<cable_segment>(
- segmentKind::dendrite,
+ section_kind::dendrite,
cable_radius, cable_radius, cable_length
);
c.add_cable(0, std::move(seg));
@@ -75,7 +75,7 @@ TEST(cell_type, add_segment)
c.add_cable(
0,
- segmentKind::dendrite, cable_radius, cable_radius, cable_length
+ section_kind::dendrite, cable_radius, cable_radius, cable_length
);
EXPECT_EQ(c.num_segments(), 2u);
@@ -92,7 +92,7 @@ TEST(cell_type, add_segment)
c.add_cable(
0,
- segmentKind::dendrite,
+ section_kind::dendrite,
std::vector<double>{cable_radius, cable_radius, cable_radius, cable_radius},
std::vector<double>{cable_length, cable_length, cable_length}
);
@@ -108,7 +108,7 @@ TEST(cell_type, multiple_cables)
// generate a cylindrical cable segment of length 1/pi and radius 1
// volume = 1
// area = 2
- auto seg = [](segmentKind k) {
+ auto seg = [](section_kind k) {
return make_segment<cable_segment>( k, 1.0, 1.0, 1./math::pi<double>() );
};
@@ -132,10 +132,10 @@ TEST(cell_type, multiple_cables)
c.add_soma(soma_radius, {0,0,1});
// hook the dendrite and axons
- c.add_cable(0, seg(segmentKind::dendrite));
- c.add_cable(0, seg(segmentKind::axon));
- c.add_cable(1, seg(segmentKind::dendrite));
- c.add_cable(1, seg(segmentKind::dendrite));
+ c.add_cable(0, seg(section_kind::dendrite));
+ c.add_cable(0, seg(section_kind::axon));
+ c.add_cable(1, seg(section_kind::dendrite));
+ c.add_cable(1, seg(section_kind::dendrite));
EXPECT_EQ(c.num_segments(), 5u);
// each of the 5 segments has volume 1 by design
@@ -173,9 +173,9 @@ TEST(cell_type, clone)
cell c;
c.add_soma(2.1);
- c.add_cable(0, segmentKind::dendrite, 0.3, 0.2, 10);
+ c.add_cable(0, section_kind::dendrite, 0.3, 0.2, 10);
c.segment(1)->set_compartments(3);
- c.add_cable(1, segmentKind::dendrite, 0.2, 0.15, 20);
+ c.add_cable(1, section_kind::dendrite, 0.2, 0.15, 20);
c.segment(2)->set_compartments(5);
parameter_list exp_default("expsyn");
@@ -208,7 +208,7 @@ TEST(cell_type, clone)
// check clone is independent
- c.add_cable(2, segmentKind::dendrite, 0.15, 0.1, 20);
+ c.add_cable(2, section_kind::dendrite, 0.15, 0.1, 20);
EXPECT_NE(c.num_segments(), d.num_segments());
d.detectors()[0].threshold = 13.0;
diff --git a/tests/unit/test_fvm_multi.cpp b/tests/unit/test_fvm_multi.cpp
index 8c08b29470c1c81d4a3351a3b7d770dedd360572..7eb0606487239d2372a699bf23b95855dec7ebd5 100644
--- a/tests/unit/test_fvm_multi.cpp
+++ b/tests/unit/test_fvm_multi.cpp
@@ -253,9 +253,9 @@ TEST(fvm_multi, mechanism_indexes)
soma->add_mechanism(hh_parameters());
// add dendrite of length 200 um and diameter 1 um with passive channel
- c.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 100);
- c.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100);
- c.add_cable(1, segmentKind::dendrite, 0.5, 0.5, 100);
+ c.add_cable(0, section_kind::dendrite, 0.5, 0.5, 100);
+ c.add_cable(1, section_kind::dendrite, 0.5, 0.5, 100);
+ c.add_cable(1, section_kind::dendrite, 0.5, 0.5, 100);
auto& segs = c.segments();
segs[1]->add_mechanism(pas_parameters());
diff --git a/tests/unit/test_probe.cpp b/tests/unit/test_probe.cpp
index 7ebbefc7adfd0c81de952479cc94d8d2f03d2e88..8027fc2fef77e3dfb194a9312889452066db534c 100644
--- a/tests/unit/test_probe.cpp
+++ b/tests/unit/test_probe.cpp
@@ -43,7 +43,7 @@ TEST(probe, fvm_multicell)
// ball-and-stick model morphology
bs.add_soma(12.6157/2.0);
- bs.add_cable(0, segmentKind::dendrite, 0.5, 0.5, 200);
+ bs.add_cable(0, section_kind::dendrite, 0.5, 0.5, 200);
bs.soma()->set_compartments(5);
segment_location loc0{0, 0};
diff --git a/tests/unit/test_segment.cpp b/tests/unit/test_segment.cpp
index d2c52959020945e022f2606933133ea51b6be50f..97f0adc98544bb0a5239ccc0075c76cef69c79bd 100644
--- a/tests/unit/test_segment.cpp
+++ b/tests/unit/test_segment.cpp
@@ -14,7 +14,7 @@ TEST(segments, soma)
EXPECT_EQ(s->volume(), pi<double>()*4./3.);
EXPECT_EQ(s->area(), pi<double>()*4.);
- EXPECT_EQ(s->kind(), segmentKind::soma);
+ EXPECT_EQ(s->kind(), section_kind::soma);
}
{
@@ -22,7 +22,7 @@ TEST(segments, soma)
EXPECT_EQ(s->volume(), pi<double>()*4./3.);
EXPECT_EQ(s->area(), pi<double>()*4.);
- EXPECT_EQ(s->kind(), segmentKind::soma);
+ EXPECT_EQ(s->kind(), section_kind::soma);
}
}
@@ -38,25 +38,25 @@ TEST(segments, cable)
// single cylindrical frustrum
{
- auto s = make_segment<cable_segment>(segmentKind::dendrite, radius, radius, length);
+ auto s = make_segment<cable_segment>(section_kind::dendrite, radius, radius, length);
EXPECT_EQ(s->volume(), 1.0);
EXPECT_EQ(s->area(), 2.0);
- EXPECT_EQ(s->kind(), segmentKind::dendrite);
+ EXPECT_EQ(s->kind(), section_kind::dendrite);
}
// cable made up of three identical cylindrical frustrums
{
auto s =
make_segment<cable_segment>(
- segmentKind::axon,
+ section_kind::axon,
std::vector<double>{radius, radius, radius, radius},
std::vector<double>{length, length, length}
);
EXPECT_EQ(s->volume(), 3.0);
EXPECT_EQ(s->area(), 6.0);
- EXPECT_EQ(s->kind(), segmentKind::axon);
+ EXPECT_EQ(s->kind(), section_kind::axon);
}
}
@@ -70,14 +70,14 @@ TEST(segments, cable_positions)
{
auto s =
make_segment<cable_segment>(
- segmentKind::dendrite,
+ section_kind::dendrite,
1, 2,
point<double>(0,0,0), point<double>(0,1,0)
);
EXPECT_EQ(s->volume(), math::volume_frustrum(1., 1., 2.));
EXPECT_EQ(s->area(), math::area_frustrum (1., 1., 2.));
- EXPECT_EQ(s->kind(), segmentKind::dendrite);
+ EXPECT_EQ(s->kind(), section_kind::dendrite);
}
// cable made up of three frustrums
@@ -85,13 +85,13 @@ TEST(segments, cable_positions)
{
auto s =
make_segment<cable_segment>(
- segmentKind::axon,
+ section_kind::axon,
std::vector<double>{1, 1.5, 2},
std::vector<point<double>>{ {0,0,0}, {0,0.5,0}, {0,1,0} }
);
EXPECT_EQ(s->volume(), math::volume_frustrum(1., 1., 2.));
EXPECT_EQ(s->area(), math::area_frustrum(1., 1., 2.));
- EXPECT_EQ(s->kind(), segmentKind::axon);
+ EXPECT_EQ(s->kind(), section_kind::axon);
}
}
diff --git a/tests/unit/test_swcio.cpp b/tests/unit/test_swcio.cpp
index 99fd5750c09fd3aacd8eb59724ae6c81463e91df..9917bbfe12bfd2f1e11561f667607ae57cb872fa 100644
--- a/tests/unit/test_swcio.cpp
+++ b/tests/unit/test_swcio.cpp
@@ -1,6 +1,7 @@
#include <array>
#include <exception>
#include <iostream>
+#include <iterator>
#include <fstream>
#include <numeric>
#include <type_traits>
@@ -16,17 +17,19 @@
# define DATADIR "../data"
#endif
+using namespace nest::mc;
+
// SWC tests
-void expect_record_equals(const nest::mc::io::swc_record& expected,
- const nest::mc::io::swc_record& actual)
+void expect_record_equals(const io::swc_record& expected,
+ const io::swc_record& actual)
{
- EXPECT_EQ(expected.id(), actual.id());
- EXPECT_EQ(expected.type(), actual.type());
- EXPECT_FLOAT_EQ(expected.x(), actual.x());
- EXPECT_FLOAT_EQ(expected.y(), actual.y());
- EXPECT_FLOAT_EQ(expected.z(), actual.z());
- EXPECT_FLOAT_EQ(expected.radius(), actual.radius());
- EXPECT_EQ(expected.parent(), actual.parent());
+ EXPECT_EQ(expected.id, actual.id);
+ EXPECT_EQ(expected.type, actual.type);
+ EXPECT_FLOAT_EQ(expected.x, actual.x);
+ EXPECT_FLOAT_EQ(expected.y, actual.y);
+ EXPECT_FLOAT_EQ(expected.z, actual.z);
+ EXPECT_FLOAT_EQ(expected.r, actual.r);
+ EXPECT_EQ(expected.parent_id, actual.parent_id);
}
TEST(swc_record, construction)
@@ -36,56 +39,57 @@ TEST(swc_record, construction)
{
// force an invalid type
swc_record::kind invalid_type = static_cast<swc_record::kind>(100);
- EXPECT_THROW(swc_record record(invalid_type, 7, 1., 1., 1., 1., 5),
- std::invalid_argument);
+ EXPECT_THROW(swc_record(invalid_type, 7, 1., 1., 1., 1., 5).assert_consistent(),
+ swc_error);
}
{
// invalid id
- EXPECT_THROW(swc_record record(
- swc_record::kind::custom, -3, 1., 1., 1., 1., 5),
- std::invalid_argument);
+ EXPECT_THROW(swc_record(
+ swc_record::kind::custom, -3, 1., 1., 1., 1., 5).assert_consistent(),
+ swc_error);
}
{
// invalid parent id
- EXPECT_THROW(swc_record record(
- swc_record::kind::custom, 0, 1., 1., 1., 1., -5),
- std::invalid_argument);
+ EXPECT_THROW(swc_record(
+ swc_record::kind::custom, 0, 1., 1., 1., 1., -5).assert_consistent(),
+ swc_error);
}
{
// invalid radius
- EXPECT_THROW(swc_record record(
- swc_record::kind::custom, 0, 1., 1., 1., -1., -1),
- std::invalid_argument);
+ EXPECT_THROW(swc_record(
+ swc_record::kind::custom, 0, 1., 1., 1., -1., -1).assert_consistent(),
+ swc_error);
}
{
// parent_id > id
- EXPECT_THROW(swc_record record(
- swc_record::kind::custom, 0, 1., 1., 1., 1., 2),
- std::invalid_argument);
+ EXPECT_THROW(swc_record(
+ swc_record::kind::custom, 0, 1., 1., 1., 1., 2).assert_consistent(),
+ swc_error);
}
{
// parent_id == id
- EXPECT_THROW(swc_record record(
- swc_record::kind::custom, 0, 1., 1., 1., 1., 0),
- std::invalid_argument);
+ EXPECT_THROW(swc_record(
+ swc_record::kind::custom, 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);
- EXPECT_EQ(record.id(), 0);
- EXPECT_EQ(record.type(), swc_record::kind::custom);
- EXPECT_EQ(record.x(), 1.);
- EXPECT_EQ(record.y(), 1.);
- EXPECT_EQ(record.z(), 1.);
- EXPECT_EQ(record.radius(), 1.);
+ EXPECT_TRUE(record.is_consistent());
+ EXPECT_EQ(record.id, 0);
+ EXPECT_EQ(record.type, swc_record::kind::custom);
+ EXPECT_EQ(record.x, 1.);
+ EXPECT_EQ(record.y, 1.);
+ EXPECT_EQ(record.z, 1.);
+ EXPECT_EQ(record.r, 1.);
EXPECT_EQ(record.diameter(), 2*1.);
- EXPECT_EQ(record.parent(), -1);
+ EXPECT_EQ(record.parent_id, -1);
}
{
@@ -96,47 +100,51 @@ TEST(swc_record, construction)
}
}
-TEST(swc_record, comparison)
+
+TEST(swc_parser, invalid_input_istream)
{
using namespace nest::mc::io;
{
- // check comparison operators
- swc_record record0(swc_record::kind::custom, 0, 1., 1., 1., 1., -1);
- swc_record record1(swc_record::kind::custom, 0, 2., 3., 4., 5., -1);
- swc_record record2(swc_record::kind::custom, 1, 2., 3., 4., 5., -1);
- EXPECT_EQ(record0, record1);
- EXPECT_LT(record0, record2);
- EXPECT_GT(record2, record1);
+ // check incomplete lines; missing parent
+ std::istringstream is("1 1 14.566132 34.873772 7.857000 0.717830\n");
+ swc_record record;
+ is >> record;
+ EXPECT_TRUE(is.fail());
}
+ {
+ // Check non-parsable values
+ std::istringstream is(
+ "1a 1 14.566132 34.873772 7.857000 0.717830 -1\n");
+ swc_record record;
+ is >> record;
+ EXPECT_TRUE(is.fail());
+ }
}
-TEST(swc_parser, invalid_input)
+TEST(swc_parser, invalid_input_parse)
{
using namespace nest::mc::io;
{
// check incomplete lines; missing parent
std::istringstream is("1 1 14.566132 34.873772 7.857000 0.717830\n");
- swc_record record;
- EXPECT_THROW(is >> record, swc_parse_error);
+ EXPECT_THROW(parse_swc_file(is), swc_error);
}
{
// Check non-parsable values
std::istringstream is(
"1a 1 14.566132 34.873772 7.857000 0.717830 -1\n");
- swc_record record;
- EXPECT_THROW(is >> record, swc_parse_error);
+ 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");
- swc_record record;
- EXPECT_THROW(is >> record, swc_parse_error);
+ EXPECT_THROW(parse_swc_file(is), swc_error);
}
{
@@ -152,20 +160,10 @@ TEST(swc_parser, invalid_input)
is << "3 1 14.566132 34.873772 7.857000 0.717830 1\n";
is << "4 1 14.566132 34.873772 7.857000 0.717830 2\n";
- std::vector<swc_record> records;
- try {
- for (auto c : swc_get_records<swc_io_clean>(is)) {
- records.push_back(std::move(c));
- }
-
- FAIL() << "expected swc_parse_error, none was thrown\n";
- } catch (const swc_parse_error& e) {
- SUCCEED();
- }
+ EXPECT_THROW(parse_swc_file(is), swc_error);
}
}
-
TEST(swc_parser, valid_input)
{
using namespace nest::mc::io;
@@ -174,7 +172,11 @@ TEST(swc_parser, valid_input)
// check empty file; no record may be parsed
swc_record record, record_orig;
std::istringstream is("");
- EXPECT_NO_THROW(is >> record);
+ is >> record;
+
+ EXPECT_TRUE(is.eof());
+ EXPECT_TRUE(is.fail());
+ EXPECT_FALSE(is.bad());
expect_record_equals(record_orig, record);
}
@@ -183,15 +185,23 @@ TEST(swc_parser, valid_input)
// no record may be parsed
swc_record record, record_orig;
std::istringstream is("#comment\n#comment");
- EXPECT_NO_THROW(is >> record);
+ is >> record;
+
+ EXPECT_TRUE(is.eof());
+ EXPECT_TRUE(is.fail());
+ EXPECT_FALSE(is.bad());
expect_record_equals(record_orig, record);
}
{
// check comment not starting at first character
swc_record record, record_orig;
- std::istringstream is(" #comment");
- EXPECT_NO_THROW(is >> record);
+ std::istringstream is(" #comment\n");
+ is >> record;
+
+ EXPECT_TRUE(is.eof());
+ EXPECT_TRUE(is.fail());
+ EXPECT_FALSE(is.bad());
expect_record_equals(record_orig, record);
}
@@ -203,11 +213,12 @@ TEST(swc_parser, valid_input)
is << " \t\n";
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
- EXPECT_NO_THROW(is >> record);
+ is >> record;
+ EXPECT_TRUE(is);
+
swc_record record_expected(
swc_record::kind::soma,
0, 14.566132, 34.873772, 7.857000, 0.717830, -1);
-
expect_record_equals(record_expected, record);
}
@@ -219,39 +230,44 @@ TEST(swc_parser, valid_input)
is << "\r\n";
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\r\n";
- EXPECT_NO_THROW(is >> record);
+ is >> record;
+ EXPECT_TRUE(is);
+
swc_record record_expected(
swc_record::kind::soma,
0, 14.566132, 34.873772, 7.857000, 0.717830, -1);
-
expect_record_equals(record_expected, record);
}
{
// check old-style mac eol; these eol are treated as simple whitespace
- // characters, so in the following case no parse error shall be thrown
- // and no record shall be read
- swc_record record, record_expected;
+ // characters, so should look line a long comment.
+ swc_record record;
std::stringstream is;
is << "#comment\r";
is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\r";
- EXPECT_NO_THROW(is >> record);
- expect_record_equals(record_expected, record);
+ is >> record;
+ EXPECT_TRUE(is.eof());
+ EXPECT_TRUE(is.fail());
+ EXPECT_FALSE(is.bad());
}
{
// check last line case (no newline at the end)
std::istringstream is("1 1 14.566132 34.873772 7.857000 0.717830 -1");
swc_record record;
- EXPECT_NO_THROW(is >> record);
- EXPECT_EQ(0, record.id()); // zero-based indexing
- EXPECT_EQ(swc_record::kind::soma, record.type());
- EXPECT_FLOAT_EQ(14.566132, record.x());
- EXPECT_FLOAT_EQ(34.873772, record.y());
- EXPECT_FLOAT_EQ( 7.857000, record.z());
- EXPECT_FLOAT_EQ( 0.717830, record.radius());
- EXPECT_FLOAT_EQ( -1, record.parent());
+ is >> record;
+ EXPECT_TRUE(is.eof());
+ 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_FLOAT_EQ(14.566132, record.x);
+ EXPECT_FLOAT_EQ(34.873772, record.y);
+ EXPECT_FLOAT_EQ( 7.857000, record.z);
+ EXPECT_FLOAT_EQ( 0.717830, record.r);
+ EXPECT_FLOAT_EQ( -1, record.parent_id);
}
{
@@ -271,18 +287,22 @@ TEST(swc_parser, valid_input)
swc_input << "# this is a final comment\n";
+ using swc_iter = std::istream_iterator<swc_record>;
+ swc_iter end;
+
std::size_t nr_records = 0;
- for (auto record : swc_get_records<swc_io_raw>(swc_input)) {
+ for (swc_iter i = swc_iter(swc_input); i!=end; ++i) {
ASSERT_LT(nr_records, records_orig.size());
- expect_record_equals(records_orig[nr_records], record);
+ expect_record_equals(records_orig[nr_records], *i);
++nr_records;
}
+ EXPECT_EQ(2u, nr_records);
}
}
TEST(swc_parser, from_allen_db)
{
- using namespace nest::mc;
+ using namespace nest::mc::io;
std::string datadir{DATADIR};
auto fname = datadir + "/example.swc";
@@ -293,13 +313,10 @@ TEST(swc_parser, from_allen_db)
}
// load the record records into a std::vector
- std::vector<io::swc_record> nodes;
- for (auto node : io::swc_get_records<io::swc_io_raw>(fid)) {
- nodes.push_back(std::move(node));
- }
+ std::vector<swc_record> nodes = parse_swc_file(fid);
// verify that the correct number of nodes was read
- EXPECT_EQ(nodes.size(), 1058u);
+ EXPECT_EQ(1058u, nodes.size());
}
TEST(swc_parser, input_cleaning)
@@ -314,7 +331,7 @@ TEST(swc_parser, input_cleaning)
is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
is << "2 2 14.566132 34.873772 7.857000 0.717830 1\n";
- EXPECT_EQ(2u, swc_get_records<swc_io_clean>(is).size());
+ EXPECT_THROW(parse_swc_file(is), swc_error);
}
{
@@ -325,8 +342,7 @@ TEST(swc_parser, input_cleaning)
is << "3 1 14.566132 34.873772 7.857000 0.717830 -1\n";
is << "4 2 14.566132 34.873772 7.857000 0.717830 1\n";
- auto records = swc_get_records<swc_io_clean>(is);
- EXPECT_EQ(2u, records.size());
+ EXPECT_THROW(parse_swc_file(is), swc_error);
}
{
@@ -339,14 +355,12 @@ TEST(swc_parser, input_cleaning)
std::array<swc_record::id_type, 4> expected_id_list = {{ 0, 1, 2, 3 }};
- auto expected_id = expected_id_list.cbegin();
- for (auto c : swc_get_records<swc_io_clean>(is)) {
- EXPECT_EQ(*expected_id, c.id());
- ++expected_id;
- }
+ auto records = parse_swc_file(is);
+ ASSERT_EQ(expected_id_list.size(), records.size());
- // Check that we have read through the whole input
- EXPECT_EQ(expected_id_list.end(), expected_id);
+ for (unsigned i = 0; i< expected_id_list.size(); ++i) {
+ EXPECT_EQ(expected_id_list[i], records[i].id);
+ }
}
{
@@ -364,22 +378,16 @@ TEST(swc_parser, input_cleaning)
std::array<swc_record::id_type, 6> expected_parent_list =
{{ -1, 0, 1, 1, 0, 4 }};
- auto expected_id = expected_id_list.cbegin();
- auto expected_parent = expected_parent_list.cbegin();
- for (auto c : swc_get_records<swc_io_clean>(is)) {
- EXPECT_EQ(*expected_id, c.id());
- EXPECT_EQ(*expected_parent, c.parent());
- ++expected_id;
- ++expected_parent;
+ auto records = parse_swc_file(is);
+ ASSERT_EQ(expected_id_list.size(), records.size());
+ for (unsigned i = 0; i< expected_id_list.size(); ++i) {
+ EXPECT_EQ(expected_id_list[i], records[i].id);
+ EXPECT_EQ(expected_parent_list[i], records[i].parent_id);
}
-
- // Check that we have read through the whole input
- EXPECT_EQ(expected_id_list.end(), expected_id);
- EXPECT_EQ(expected_parent_list.end(), expected_parent);
}
}
-TEST(swc_record_ranges, raw)
+TEST(swc_parser, raw)
{
using namespace nest::mc::io;
@@ -391,50 +399,11 @@ TEST(swc_record_ranges, raw)
is << "3 2 14.566132 34.873772 7.857000 0.717830 1\n";
is << "4 2 14.566132 34.873772 7.857000 0.717830 1\n";
- std::vector<swc_record> records;
- for (auto c : swc_get_records<swc_io_raw>(is)) {
- records.push_back(c);
- }
+ using swc_iter = std::istream_iterator<swc_record>;
+ std::vector<swc_record> records{swc_iter(is), swc_iter()};
EXPECT_EQ(4u, records.size());
-
- bool entered = false;
- auto citer = records.begin();
- for (auto c : swc_get_records<swc_io_raw>(is)) {
- expect_record_equals(c, *citer++);
- entered = true;
- }
-
- EXPECT_TRUE(entered);
- }
-
- {
- // Check out of bounds reads
- std::stringstream is;
- is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
-
- auto ibegin = swc_get_records<swc_io_raw>(is).begin();
-
- EXPECT_NO_THROW(++ibegin);
- EXPECT_THROW(*ibegin, std::out_of_range);
-
- }
-
- {
- // Check iterator increments
- std::stringstream is;
- is << "1 1 14.566132 34.873772 7.857000 0.717830 -1\n";
-
- auto iter = swc_get_records<swc_io_raw>(is).begin();
- auto iend = swc_get_records<swc_io_raw>(is).end();
-
- swc_record c;
- EXPECT_NO_THROW(c = *iter++);
- EXPECT_EQ(-1, c.parent());
- EXPECT_EQ(iend, iter);
-
- // Try to read past eof
- EXPECT_THROW(*iter, std::out_of_range);
+ EXPECT_EQ(3, records.back().id);
}
{
@@ -445,105 +414,108 @@ TEST(swc_record_ranges, raw)
is << "3 10 14.566132 34.873772 7.857000 0.717830 1\n";
is << "4 2 14.566132 34.873772 7.857000 0.717830 1\n";
- std::vector<swc_record> records;
try {
- for (auto c : swc_get_records<swc_io_raw>(is)) {
- records.push_back(c);
- }
-
+ parse_swc_file(is);
ADD_FAILURE() << "expected an exception\n";
- } catch (const swc_parse_error& e) {
- EXPECT_EQ(3u, e.lineno());
+ }
+ catch (const swc_error& e) {
+ EXPECT_EQ(3u, e.line_number);
}
}
{
// Test empty range
std::stringstream is("");
- EXPECT_TRUE(swc_get_records<swc_io_raw>(is).empty());
- EXPECT_TRUE(swc_get_records<swc_io_clean>(is).empty());
+ using swc_iter = std::istream_iterator<swc_record>;
+ std::vector<swc_record> records{swc_iter(is), swc_iter()};
+
+ EXPECT_TRUE(records.empty());
}
}
-TEST(swc_io, cell_construction)
-{
+TEST(swc_io, cell_construction) {
using namespace nest::mc;
- {
- //
- // 0
- // |
- // 1
- // |
- // 2
- // / \.
- // 3 4
- // \.
- // 5
- //
-
- std::stringstream is;
- is << "1 1 0 0 0 2.1 -1\n";
- is << "2 3 0.1 1.2 1.2 1.3 1\n";
- is << "3 3 1.0 2.0 2.2 1.1 2\n";
- is << "4 3 1.5 3.3 1.3 2.2 3\n";
- is << "5 3 2.5 5.3 2.5 0.7 3\n";
- is << "6 3 3.5 2.3 3.7 3.4 5\n";
-
- using point_type = point<double>;
- std::vector<point_type> points = {
- { 0.0, 0.0, 0.0 },
- { 0.1, 1.2, 1.2 },
- { 1.0, 2.0, 2.2 },
- { 1.5, 3.3, 1.3 },
- { 2.5, 5.3, 2.5 },
- { 3.5, 2.3, 3.7 },
- };
-
- cell cell = io::swc_read_cell(is);
- EXPECT_TRUE(cell.has_soma());
- EXPECT_EQ(4u, cell.num_segments());
-
- EXPECT_EQ(norm(points[1]-points[2]), cell.cable(1)->length());
- EXPECT_EQ(norm(points[2]-points[3]), cell.cable(2)->length());
- EXPECT_EQ(norm(points[2]-points[4]) + norm(points[4]-points[5]),
- cell.cable(3)->length());
-
-
- // Check each segment separately
- EXPECT_TRUE(cell.segment(0)->is_soma());
- EXPECT_EQ(2.1, cell.soma()->radius());
- EXPECT_EQ(point_type(0, 0, 0), cell.soma()->center());
-
- for (auto i = 1u; i < cell.num_segments(); ++i) {
- EXPECT_TRUE(cell.segment(i)->is_dendrite());
- }
-
- EXPECT_EQ(1u, cell.cable(1)->num_sub_segments());
- EXPECT_EQ(1u, cell.cable(2)->num_sub_segments());
- EXPECT_EQ(2u, cell.cable(3)->num_sub_segments());
-
-
- // Check the radii
- EXPECT_EQ(1.3, cell.cable(1)->radius(0));
- EXPECT_EQ(1.1, cell.cable(1)->radius(1));
-
- EXPECT_EQ(1.1, cell.cable(2)->radius(0));
- EXPECT_EQ(2.2, cell.cable(2)->radius(1));
-
- EXPECT_EQ(1.1, cell.cable(3)->radius(0));
- EXPECT_EQ(3.4, cell.cable(3)->radius(1));
-
- auto len_ratio = norm(points[2]-points[4]) / cell.cable(3)->length();
- EXPECT_NEAR(.7, cell.cable(3)->radius(len_ratio), 1e-15);
-
- // Double-check radii at joins are equal
- EXPECT_EQ(cell.cable(1)->radius(1),
- cell.cable(2)->radius(0));
-
- EXPECT_EQ(cell.cable(1)->radius(1),
- cell.cable(3)->radius(0));
- }
+ //
+ // 0
+ // |
+ // 1
+ // |
+ // 2
+ // / \.
+ // 3 4
+ // \.
+ // 5
+ //
+
+ std::stringstream is;
+ is << "1 1 0 0 0 2.1 -1\n";
+ is << "2 3 0.1 1.2 1.2 1.3 1\n";
+ is << "3 3 1.0 2.0 2.2 1.1 2\n";
+ is << "4 3 1.5 3.3 1.3 2.2 3\n";
+ is << "5 3 2.5 5.3 2.5 0.7 3\n";
+ is << "6 3 3.5 2.3 3.7 3.4 5\n";
+
+ using point_type = point<double>;
+ std::vector<point_type> points = {
+ { 0.0, 0.0, 0.0 },
+ { 0.1, 1.2, 1.2 },
+ { 1.0, 2.0, 2.2 },
+ { 1.5, 3.3, 1.3 },
+ { 2.5, 5.3, 2.5 },
+ { 3.5, 2.3, 3.7 },
+ };
+
+ // swc -> morphology
+ auto morph = io::swc_as_morphology(io::parse_swc_file(is));
+
+ cell cell = make_cell(morph, true);
+ EXPECT_TRUE(cell.has_soma());
+ EXPECT_EQ(4u, cell.num_segments());
+
+ EXPECT_EQ(norm(points[1]-points[2]), cell.cable(1)->length());
+ EXPECT_EQ(norm(points[2]-points[3]), cell.cable(2)->length());
+ EXPECT_EQ(norm(points[2]-points[4]) + norm(points[4]-points[5]),
+ cell.cable(3)->length());
+
+
+ // Check each segment separately
+ EXPECT_TRUE(cell.segment(0)->is_soma());
+ EXPECT_EQ(2.1, cell.soma()->radius());
+ EXPECT_EQ(point_type(0, 0, 0), cell.soma()->center());
+
+ for (auto i = 1u; i < cell.num_segments(); ++i) {
+ EXPECT_TRUE(cell.segment(i)->is_dendrite());
+ }
+
+ EXPECT_EQ(1u, cell.cable(1)->num_sub_segments());
+ EXPECT_EQ(1u, cell.cable(2)->num_sub_segments());
+ EXPECT_EQ(2u, cell.cable(3)->num_sub_segments());
+
+ // We asked to use the same discretization as in the SWC, so check number of compartments too.
+ EXPECT_EQ(1u, cell.cable(1)->num_compartments());
+ EXPECT_EQ(1u, cell.cable(2)->num_compartments());
+ EXPECT_EQ(2u, cell.cable(3)->num_compartments());
+
+ // Check the radii
+ EXPECT_EQ(1.3, cell.cable(1)->radius(0));
+ EXPECT_EQ(1.1, cell.cable(1)->radius(1));
+
+ EXPECT_EQ(1.1, cell.cable(2)->radius(0));
+ EXPECT_EQ(2.2, cell.cable(2)->radius(1));
+
+ EXPECT_EQ(1.1, cell.cable(3)->radius(0));
+ EXPECT_EQ(3.4, cell.cable(3)->radius(1));
+
+ auto len_ratio = norm(points[2]-points[4]) / cell.cable(3)->length();
+ EXPECT_NEAR(.7, cell.cable(3)->radius(len_ratio), 1e-15);
+
+ // Double-check radii at joins are equal
+ EXPECT_EQ(cell.cable(1)->radius(1),
+ cell.cable(2)->radius(0));
+
+ EXPECT_EQ(cell.cable(1)->radius(1),
+ cell.cable(3)->radius(0));
}
// check that simple ball and stick model with one dendrite attached to a soma
@@ -560,16 +532,16 @@ TEST(swc_parser, from_file_ball_and_stick)
}
// read the file into a cell object
- auto cell = nest::mc::io::swc_read_cell(fid);
+ auto bas_cell = make_cell(io::swc_as_morphology(io::parse_swc_file(fid)));
// verify that the correct number of nodes was read
- EXPECT_EQ(cell.num_segments(), 2u);
- EXPECT_EQ(cell.num_compartments(), 2u);
+ EXPECT_EQ(2u, bas_cell.num_segments());
+ EXPECT_EQ(2u, bas_cell.num_compartments());
// make an equivalent cell via C++ interface
- nest::mc::cell local_cell;
+ cell local_cell;
local_cell.add_soma(6.30785);
- local_cell.add_cable(0, nest::mc::segmentKind::dendrite, 0.5, 0.5, 200);
+ local_cell.add_cable(0, section_kind::dendrite, 0.5, 0.5, 200);
- EXPECT_TRUE(nest::mc::cell_basic_equality(local_cell, cell));
+ EXPECT_TRUE(cell_basic_equality(local_cell, bas_cell));
}