From f91d0b3ccd62ff6826679977a41203ef63735e6d Mon Sep 17 00:00:00 2001
From: Sam Yates <yates@cscs.ch>
Date: Thu, 9 Mar 2017 16:48:26 +0100
Subject: [PATCH] Morphologies in miniapp (#178)
Fix morphology section ctor bug.
Add morphology pools for miniapp from which morphologies are drawn in the recipe.
Add command-line options to expose above.
Add option --report-compartments to (slowly) check the min, mean, and max number of compartments in the generated cells across the simulation.
Use morphology::add_section to do all the heavy lifting; no need to make section_geometry objects by hand, unless you really, really want to.
---
lmorpho/lsystem.cpp | 2 +-
miniapp/CMakeLists.txt | 2 +
miniapp/io.cpp | 35 +++++++++++++--
miniapp/io.hpp | 3 ++
miniapp/miniapp.cpp | 28 ++++++++++++
miniapp/miniapp_recipes.cpp | 87 +++++++++++++++++++++++++------------
miniapp/miniapp_recipes.hpp | 16 ++++++-
miniapp/morphology_pool.cpp | 56 ++++++++++++++++++++++++
miniapp/morphology_pool.hpp | 42 ++++++++++++++++++
src/CMakeLists.txt | 4 +-
src/morphology.cpp | 11 +++--
src/morphology.hpp | 10 +++++
src/swcio.cpp | 3 +-
src/swcio.hpp | 24 +++++-----
src/util/meta.hpp | 2 +-
src/util/path.cpp | 43 ++++++++++++++++++
src/util/path.hpp | 6 ++-
src/util/scope_exit.hpp | 54 +++++++++++++++++++++++
18 files changed, 373 insertions(+), 55 deletions(-)
create mode 100644 miniapp/morphology_pool.cpp
create mode 100644 miniapp/morphology_pool.hpp
create mode 100644 src/util/path.cpp
create mode 100644 src/util/scope_exit.hpp
diff --git a/lmorpho/lsystem.cpp b/lmorpho/lsystem.cpp
index f88ba290..24f92179 100644
--- a/lmorpho/lsystem.cpp
+++ b/lmorpho/lsystem.cpp
@@ -303,7 +303,7 @@ nest::mc::morphology generate_morphology(const lsys_param& P, lsys_generator &g)
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};
+ 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});
diff --git a/miniapp/CMakeLists.txt b/miniapp/CMakeLists.txt
index 65973201..801a48d2 100644
--- a/miniapp/CMakeLists.txt
+++ b/miniapp/CMakeLists.txt
@@ -4,11 +4,13 @@ set(MINIAPP_SOURCES
miniapp.cpp
io.cpp
miniapp_recipes.cpp
+ morphology_pool.cpp
)
set(MINIAPP_SOURCES_CUDA
miniapp.cu
io.cpp
miniapp_recipes.cpp
+ morphology_pool.cpp
)
if(NMC_WITH_CUDA)
diff --git a/miniapp/io.cpp b/miniapp/io.cpp
index 9771de35..7a6e9c67 100644
--- a/miniapp/io.cpp
+++ b/miniapp/io.cpp
@@ -104,7 +104,7 @@ static void update_option(util::optional<T>& opt, const nlohmann::json& j, const
opt = util::nothing;
}
else {
- opt = value;
+ opt = value.get<T>();
}
}
}
@@ -127,6 +127,9 @@ cl_options read_options(int argc, char** argv, bool allow_write) {
0.0, // probe_ratio
"trace_", // trace_prefix
util::nothing, // trace_max_gid
+ util::nothing, // morphologies
+ false, // morph_rr;
+ false, // report_compartments;
// spike_output_parameters:
false, // spike output
@@ -194,13 +197,18 @@ cl_options read_options(int argc, char** argv, bool allow_write) {
TCLAP::ValueArg<util::optional<unsigned>> trace_max_gid_arg(
"T", "trace-max-gid", "only trace probes on cells up to and including <gid>",
false, defopts.trace_max_gid, "gid", cmd);
+ TCLAP::ValueArg<util::optional<std::string>> morphologies_arg(
+ "M", "morphologies", "load morphologies from SWC files matching <glob>",
+ false, defopts.morphologies, "glob", cmd);
+ TCLAP::SwitchArg morph_rr_arg(
+ "", "morph-rr", "Serial rather than random morphology selection from pool", cmd, false);
+ TCLAP::SwitchArg report_compartments_arg(
+ "", "report-compartments", "Count compartments in cells before simulation", cmd, false);
TCLAP::SwitchArg spike_output_arg(
"f","spike-file-output","save spikes to file", cmd, false);
-
TCLAP::ValueArg<unsigned> dry_run_ranks_arg(
"D","dry-run-ranks","number of ranks in dry run mode",
false, defopts.dry_run_ranks, "positive integer", cmd);
-
TCLAP::SwitchArg profile_only_zero_arg(
"z", "profile-only-zero", "Only output profile information for rank 0", cmd, false);
@@ -232,6 +240,9 @@ cl_options read_options(int argc, char** argv, bool allow_write) {
update_option(options.probe_soma_only, fopts, "probe_soma_only");
update_option(options.trace_prefix, fopts, "trace_prefix");
update_option(options.trace_max_gid, fopts, "trace_max_gid");
+ update_option(options.morphologies, fopts, "morphologies");
+ update_option(options.morph_rr, fopts, "morph_rr");
+ update_option(options.report_compartments, fopts, "report_compartments");
// Parameters for spike output
update_option(options.spike_file_output, fopts, "spike_file_output");
@@ -271,6 +282,9 @@ cl_options read_options(int argc, char** argv, bool allow_write) {
update_option(options.probe_soma_only, probe_soma_only_arg);
update_option(options.trace_prefix, trace_prefix_arg);
update_option(options.trace_max_gid, trace_max_gid_arg);
+ update_option(options.morphologies, morphologies_arg);
+ update_option(options.morph_rr, morph_rr_arg);
+ update_option(options.report_compartments, report_compartments_arg);
update_option(options.spike_file_output, spike_output_arg);
update_option(options.profile_only_zero, profile_only_zero_arg);
update_option(options.dry_run_ranks, dry_run_ranks_arg);
@@ -314,6 +328,14 @@ cl_options read_options(int argc, char** argv, bool allow_write) {
else {
fopts["trace_max_gid"] = nullptr;
}
+ if (options.morphologies) {
+ fopts["morphologies"] = options.morphologies.get();
+ }
+ else {
+ fopts["morphologies"] = nullptr;
+ }
+ fopts["morph_rr"] = options.morph_rr;
+ fopts["report_compartments"] = options.report_compartments;
fid << std::setw(3) << fopts << "\n";
}
@@ -347,6 +369,13 @@ std::ostream& operator<<(std::ostream& o, const cl_options& options) {
o << *options.trace_max_gid;
}
o << "\n";
+ o << " morphologies : ";
+ if (options.morphologies) {
+ o << *options.morphologies;
+ }
+ o << "\n";
+ o << " morphology r-r : " << (options.morph_rr ? "yes" : "no") << "\n";
+ o << " report compartments : " << (options.report_compartments ? "yes" : "no") << "\n";
return o;
}
diff --git a/miniapp/io.hpp b/miniapp/io.hpp
index bf6b23ef..cd3eca66 100644
--- a/miniapp/io.hpp
+++ b/miniapp/io.hpp
@@ -27,6 +27,9 @@ struct cl_options {
double probe_ratio;
std::string trace_prefix;
util::optional<unsigned> trace_max_gid;
+ util::optional<std::string> morphologies;
+ bool morph_rr;
+ bool report_compartments;
// Parameters for spike output
bool spike_file_output;
diff --git a/miniapp/miniapp.cpp b/miniapp/miniapp.cpp
index 8d520c9f..a865f42e 100644
--- a/miniapp/miniapp.cpp
+++ b/miniapp/miniapp.cpp
@@ -45,6 +45,7 @@ using communicator_type = communication::communicator<model_type::time_type, com
using spike_type = typename communicator_type::spike_type;
void write_trace_json(const sample_trace_type& trace, const std::string& prefix = "trace_");
+void report_compartment_stats(const recipe&);
int main(int argc, char** argv) {
nest::mc::communication::global_policy_guard global_guard(argc, argv);
@@ -97,6 +98,9 @@ int main(int argc, char** argv) {
};
model_type m(*recipe, util::partition_view(group_divisions));
+ if (options.report_compartments) {
+ report_compartment_stats(*recipe);
+ }
// inject some artificial spikes, 1 per 20 neurons.
std::vector<cell_gid_type> local_sources;
@@ -205,6 +209,14 @@ void banner() {
std::unique_ptr<recipe> make_recipe(const io::cl_options& options, const probe_distribution& pdist) {
basic_recipe_param p;
+ if (options.morphologies) {
+ std::cout << "loading morphologies...\n";
+ p.morphologies.clear();
+ load_swc_morphology_glob(p.morphologies, options.morphologies.get());
+ std::cout << "loading morphologies: " << p.morphologies.size() << " loaded.\n";
+ }
+ p.morphology_round_robin = options.morph_rr;
+
p.num_compartments = options.compartments_per_segment;
p.num_synapses = options.all_to_all? options.cells-1: options.synapses_per_cell;
p.synapse_type = options.syn_type;
@@ -260,3 +272,19 @@ void write_trace_json(const sample_trace_type& trace, const std::string& prefix)
std::ofstream file(path);
file << std::setw(1) << jrep << std::endl;
}
+
+void report_compartment_stats(const recipe& rec) {
+std::size_t ncell = rec.num_cells();
+ std::size_t ncomp_total = 0;
+ std::size_t ncomp_min = -1;
+ std::size_t ncomp_max = 0;
+
+ for (std::size_t i = 0; i<ncell; ++i) {
+ std::size_t ncomp = rec.get_cell(i).num_compartments();
+ ncomp_total += ncomp;
+ ncomp_min = std::min(ncomp_min, ncomp);
+ ncomp_max = std::max(ncomp_max, ncomp);
+ }
+
+ std::cout << "compartments/cell: min=" << ncomp_min <<"; max=" << ncomp_max << "; mean=" << (double)ncomp_total/ncell << "\n";
+}
diff --git a/miniapp/miniapp_recipes.cpp b/miniapp/miniapp_recipes.cpp
index 62d4b48c..88474975 100644
--- a/miniapp/miniapp_recipes.cpp
+++ b/miniapp/miniapp_recipes.cpp
@@ -8,6 +8,7 @@
#include <util/debug.hpp>
#include "miniapp_recipes.hpp"
+#include "morphology_pool.hpp"
namespace nest {
namespace mc {
@@ -17,37 +18,51 @@ namespace mc {
template <typename RNG>
cell make_basic_cell(
+ const morphology& morph,
unsigned compartments_per_segment,
unsigned num_synapses,
const std::string& syn_type,
RNG& rng)
{
- nest::mc::cell cell;
-
- // Soma with diameter 12.6157 um and HH channel
- auto soma = cell.add_soma(12.6157/2.0);
- soma->add_mechanism(mc::hh_parameters());
-
- // 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, 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());
- d->set_compartments(compartments_per_segment);
- d->mechanism("membrane").set("r_L", 100);
+ nest::mc::cell cell = make_cell(morph, true);
+
+ for (auto& segment: cell.segments()) {
+ if (compartments_per_segment!=0) {
+ if (cable_segment* cable = segment->as_cable()) {
+ cable->set_compartments(compartments_per_segment);
+ }
+ }
+
+ if (segment->is_dendrite()) {
+ segment->add_mechanism(mc::pas_parameters());
+ segment->mechanism("membrane").set("r_L", 100);
+ }
}
+ cell.soma()->add_mechanism(mc::hh_parameters());
cell.add_detector({0,0}, 20);
auto distribution = std::uniform_real_distribution<float>(0.f, 1.0f);
- // distribute the synapses at random locations the terminal dendrites in a
- // round robin manner; the terminal dendrites in this cell have indices 2 and 3.
+
+ // Distribute the synapses at random locations the terminal dendrites in a
+ // round robin manner.
+
+ morph.assert_valid();
+ std::vector<unsigned> terminals;
+ for (const auto& section: morph.sections) {
+ // Note that morphology section ids should match up exactly with cell
+ // segment ids!
+ if (section.terminal) {
+ terminals.push_back(section.id);
+ }
+ }
+
+ EXPECTS(!terminals.empty());
+
nest::mc::parameter_list syn_default(syn_type);
for (unsigned i=0; i<num_synapses; ++i) {
- cell.add_synapse({2+(i%2), distribution(rng)}, syn_default);
+ unsigned id = terminals[i%terminals.size()];
+ cell.add_synapse({id, distribution(rng)}, syn_default);
}
return cell;
@@ -58,7 +73,8 @@ public:
basic_cell_recipe(cell_gid_type ncell, basic_recipe_param param, probe_distribution pdist):
ncell_(ncell), param_(std::move(param)), pdist_(std::move(pdist))
{
- delay_distribution_param = exp_param{param_.mean_connection_delay_ms
+ EXPECTS(param_.morphologies.size()>0);
+ delay_distribution_param_ = exp_param{param_.mean_connection_delay_ms
- param_.min_connection_delay_ms};
}
@@ -68,17 +84,20 @@ public:
auto gen = std::mt19937(i); // TODO: replace this with hashing generator...
auto cc = get_cell_count_info(i);
- auto cell = make_basic_cell(param_.num_compartments, cc.num_targets,
+ const auto& morph = get_morphology(i);
+ unsigned cell_segments = morph.components();
+
+ auto cell = make_basic_cell(morph, param_.num_compartments, cc.num_targets,
param_.synapse_type, gen);
- EXPECTS(cell.num_segments()==basic_cell_segments);
+ EXPECTS(cell.num_segments()==cell_segments);
EXPECTS(cell.probes().size()==0);
EXPECTS(cell.synapses().size()==cc.num_targets);
EXPECTS(cell.detectors().size()==cc.num_sources);
// add probes
if (cc.num_probes) {
- unsigned n_probe_segs = pdist_.all_segments? basic_cell_segments: 1u;
+ unsigned n_probe_segs = pdist_.all_segments? cell_segments: 1u;
for (unsigned i = 0; i<n_probe_segs; ++i) {
if (pdist_.membrane_voltage) {
cell.add_probe({{i, i? 0.5: 0.0}, mc::probeKind::membrane_voltage});
@@ -94,12 +113,13 @@ public:
cell_count_info get_cell_count_info(cell_gid_type i) const override {
cell_count_info cc = {1, param_.num_synapses, 0 };
+ unsigned cell_segments = get_morphology(i).components();
// probe this cell?
if (std::floor(i*pdist_.proportion)!=std::floor((i-1.0)*pdist_.proportion)) {
std::size_t np = pdist_.membrane_voltage + pdist_.membrane_current;
if (pdist_.all_segments) {
- np *= basic_cell_segments;
+ np *= cell_segments;
}
cc.num_probes = np;
@@ -111,7 +131,7 @@ public:
protected:
template <typename RNG>
cell_connection draw_connection_params(RNG& rng) const {
- std::exponential_distribution<float> delay_dist(delay_distribution_param);
+ std::exponential_distribution<float> delay_dist(delay_distribution_param_);
float delay = param_.min_connection_delay_ms + delay_dist(rng);
float weight = param_.syn_weight_per_cell/param_.num_synapses;
return cell_connection{{0, 0}, {0, 0}, weight, delay};
@@ -120,10 +140,23 @@ protected:
cell_gid_type ncell_;
basic_recipe_param param_;
probe_distribution pdist_;
- static constexpr int basic_cell_segments = 4;
using exp_param = std::exponential_distribution<float>::param_type;
- exp_param delay_distribution_param;
+ exp_param delay_distribution_param_;
+
+ const morphology& get_morphology(cell_gid_type gid) const {
+ // Allocate to gids sequentially?
+ if (param_.morphology_round_robin) {
+ return param_.morphologies[gid%param_.morphologies.size()];
+ }
+
+ // Morphologies are otherwise selected deterministically pseudo-randomly from pool.
+ std::uniform_int_distribution<unsigned> morph_select_dist_(0, param_.morphologies.size()-1);
+
+ // TODO: definitely replace this with a random hash!
+ auto gen = std::mt19937(gid+0xbad0cafe);
+ return param_.morphologies[morph_select_dist_(gen)];
+ }
};
class basic_ring_recipe: public basic_cell_recipe {
diff --git a/miniapp/miniapp_recipes.hpp b/miniapp/miniapp_recipes.hpp
index 2e519bd1..afc1a010 100644
--- a/miniapp/miniapp_recipes.hpp
+++ b/miniapp/miniapp_recipes.hpp
@@ -4,7 +4,9 @@
#include <memory>
#include <stdexcept>
-#include "recipe.hpp"
+#include <recipe.hpp>
+
+#include "morphology_pool.hpp"
// miniapp-specific recipes
@@ -19,12 +21,24 @@ struct probe_distribution {
};
struct basic_recipe_param {
+ // `num_compartments` is the number of compartments to place in each
+ // unbranched section of the morphology, A value of zero indicates that
+ // the number of compartments should equal the number of piecewise
+ // linear segments in the morphology description of that branch.
unsigned num_compartments = 1;
+
+ // Total number of synapses on each cell.
unsigned num_synapses = 1;
+
std::string synapse_type = "expsyn";
float min_connection_delay_ms = 20.0;
float mean_connection_delay_ms = 20.75;
float syn_weight_per_cell = 0.3;
+
+ morphology_pool morphologies = default_morphology_pool;
+
+ // If true, iterate through morphologies rather than select randomly.
+ bool morphology_round_robin = false;
};
std::unique_ptr<recipe> make_basic_ring_recipe(
diff --git a/miniapp/morphology_pool.cpp b/miniapp/morphology_pool.cpp
new file mode 100644
index 00000000..68f3a617
--- /dev/null
+++ b/miniapp/morphology_pool.cpp
@@ -0,0 +1,56 @@
+#include <fstream>
+#include <memory>
+#include <vector>
+
+#include <morphology.hpp>
+#include <swcio.hpp>
+#include <util/path.hpp>
+
+#include "morphology_pool.hpp"
+
+namespace nest {
+namespace mc {
+
+static morphology make_basic_y_morphology() {
+ morphology morph;
+
+ // soma of diameter 12.6157 microns.
+ // proximal section: 200 microns, radius 0.5 microns.
+ // two terminal branches, each: 100 microns, terminal radius 0.25 microns.
+ morph.soma.r = 12.6157/2;
+ double x = morph.soma.r;
+ morph.add_section({{x, 0, 0, 0.5}, {x+200, 0, 0, 0.5}});
+ x += 200;
+ morph.add_section({{x, 0, 0, 0.5}, {x+100, 0, 0, 0.25}}, 1u);
+ morph.add_section({{x, 0, 0, 0.5}, {x+100, 0, 0, 0.25}}, 1u);
+
+ morph.assert_valid();
+ return morph;
+}
+
+morphology_pool default_morphology_pool(make_basic_y_morphology());
+
+void load_swc_morphology(morphology_pool& pool, const util::path& swc_path) {
+ std::ifstream fi;
+ fi.exceptions(std::ifstream::failbit);
+
+ fi.open(swc_path.c_str());
+ pool.insert(io::swc_as_morphology(io::parse_swc_file(fi)));
+}
+
+void load_swc_morphology_glob(morphology_pool& pool, const std::string& swc_pattern) {
+ std::ifstream fi;
+ fi.exceptions(std::ifstream::failbit);
+
+ auto swc_paths = util::glob(swc_pattern);
+ for (const auto& p: swc_paths) {
+ fi.open(p.c_str());
+ pool.insert(io::swc_as_morphology(io::parse_swc_file(fi)));
+ pool[pool.size()-1].assert_valid();
+ fi.close();
+ }
+}
+
+
+} // namespace mc
+} // namespace nest
diff --git a/miniapp/morphology_pool.hpp b/miniapp/morphology_pool.hpp
new file mode 100644
index 00000000..30eadcc0
--- /dev/null
+++ b/miniapp/morphology_pool.hpp
@@ -0,0 +1,42 @@
+#pragma once
+
+// Maintain a pool of morphologies for use with miniapp recipes.
+// The default pool comprises a single ball-and-stick morphology;
+// sets of morphologies can be loaded from SWC files.
+
+#include <memory>
+#include <string>
+#include <vector>
+
+#include <morphology.hpp>
+#include <util/path.hpp>
+
+namespace nest {
+namespace mc {
+
+class morphology_pool {
+ std::shared_ptr<std::vector<morphology>> pool;
+
+public:
+ // Construct default empty pool.
+ morphology_pool(): pool(new std::vector<morphology>) {}
+
+ // Construct pool with one starting morphology.
+ explicit morphology_pool(morphology m): pool(new std::vector<morphology>) {
+ insert(std::move(m));
+ }
+
+ std::size_t size() const { return pool->size(); }
+ const morphology& operator[](std::ptrdiff_t i) const { return (*pool)[i]; }
+
+ void insert(morphology m) { (*pool).push_back(std::move(m)); }
+ void clear() { (*pool).clear(); }
+};
+
+extern morphology_pool default_morphology_pool;
+
+void load_swc_morphology(morphology_pool& pool, const util::path& swc_path);
+void load_swc_morphology_glob(morphology_pool& pool, const std::string& pattern);
+
+} // namespace mc
+} // namespace nest
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index b3faee92..4d600486 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,6 +1,3 @@
-set(HEADERS
- swcio.hpp
-)
set(BASE_SOURCES
common_types_io.cpp
cell.cpp
@@ -9,6 +6,7 @@ set(BASE_SOURCES
profiling/profiler.cpp
swcio.cpp
util/debug.cpp
+ util/path.cpp
util/unwind.cpp
backends/fvm_multicore.cpp
)
diff --git a/src/morphology.cpp b/src/morphology.cpp
index 740d4883..35b520ad 100644
--- a/src/morphology.cpp
+++ b/src/morphology.cpp
@@ -90,14 +90,17 @@ section_geometry& morphology::add_section(std::vector<section_point> points, uns
if (section.parent_id >= section.id) {
throw morphology_error("improper parent id for section");
}
- sections[section.parent_id].terminal = false;
+
+ 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(true, nsection);
+ std::vector<int> terminal(nsection, true);
for (std::size_t i=0; i<nsection; ++i) {
auto id = m.sections[i].id;
@@ -106,9 +109,9 @@ static const char* morphology_invariant_violation(const morphology& m) {
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;
+ 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) {
diff --git a/src/morphology.hpp b/src/morphology.hpp
index e57278f2..4faf7535 100644
--- a/src/morphology.hpp
+++ b/src/morphology.hpp
@@ -28,6 +28,11 @@ struct section_geometry {
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);
};
@@ -53,6 +58,11 @@ struct morphology {
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
diff --git a/src/swcio.cpp b/src/swcio.cpp
index 023bbf28..3b825d3e 100644
--- a/src/swcio.cpp
+++ b/src/swcio.cpp
@@ -118,12 +118,13 @@ std::ostream& operator<<(std::ostream& os, const swc_record& record) {
}
std::vector<swc_record> parse_swc_file(std::istream& is) {
+ constexpr auto eof = std::char_traits<char>::eof();
std::vector<swc_record> records;
std::size_t line_number = 1;
std::string line;
try {
- while (is) {
+ while (is && is.peek()!=eof) {
std::getline(is, line);
if (is_comment(line)) {
continue;
diff --git a/src/swcio.hpp b/src/swcio.hpp
index 019b26ff..a7e728fd 100644
--- a/src/swcio.hpp
+++ b/src/swcio.hpp
@@ -137,41 +137,39 @@ morphology swc_as_morphology(const RandomAccessSequence& swc_records) {
auto b_start = std::next(swc_records.begin(), branch_index[i]);
auto b_end = std::next(swc_records.begin(), branch_index[i+1]);
- section_geometry section;
- section.id = i;
- section.parent_id = parent_branch_index[i];
+ unsigned parent_id = parent_branch_index[i];
+ std::vector<section_point> points;
+ section_kind kind = section_kind::none;
- if (section.parent_id != 0) {
+ 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]]];
- section_point point{parent_record.x, parent_record.y, parent_record.z, parent_record.r};
- section.points.push_back(point);
+ 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) {
- section_point point{b->x, b->y, b->z, b->r};
- section.points.push_back(point);
+ points.push_back(section_point{b->x, b->y, b->z, b->r});
switch (b->type) {
case swc_record::kind::axon:
- section.kind = section_kind::axon;
+ kind = section_kind::axon;
break;
case swc_record::kind::dendrite:
case swc_record::kind::apical_dendrite:
- section.kind = section_kind::dendrite;
+ kind = section_kind::dendrite;
break;
case swc_record::kind::soma:
- section.kind = section_kind::soma;
+ kind = section_kind::soma;
break;
default: ; // stick with what we have
}
}
- morph.sections.push_back(section);
+ morph.add_section(std::move(points), parent_id, kind);
}
- EXPECTS(morph.check_valid());
+ morph.assert_valid();
return morph;
}
diff --git a/src/util/meta.hpp b/src/util/meta.hpp
index 3728a916..52c67669 100644
--- a/src/util/meta.hpp
+++ b/src/util/meta.hpp
@@ -210,7 +210,7 @@ struct is_sequence:
std::false_type {};
template<typename T>
-struct is_sequence<T, impl::sink<decltype(std::declval<T>())>>:
+struct is_sequence<T, impl::sink<decltype(cbegin(std::declval<T>()))>>:
std::true_type {};
template <typename T>
diff --git a/src/util/path.cpp b/src/util/path.cpp
new file mode 100644
index 00000000..53cc710b
--- /dev/null
+++ b/src/util/path.cpp
@@ -0,0 +1,43 @@
+#include <util/scope_exit.hpp>
+#include <util/path.hpp>
+
+// POSIX header
+#include <glob.h>
+
+namespace nest {
+namespace mc {
+namespace util {
+
+std::vector<path> glob(const std::string& pattern) {
+ std::vector<path> paths;
+
+ int flags = GLOB_MARK | GLOB_NOCHECK;
+#if defined(GLOB_TILDE)
+ flags |= GLOB_TILDE;
+#endif
+#if defined(GLOB_TILDE)
+ flags |= GLOB_BRACE;;
+#endif
+
+ glob_t matches;
+ auto r = ::glob(pattern.c_str(), flags, nullptr, &matches);
+ auto glob_guard = on_scope_exit([&]() { ::globfree(&matches); });
+
+ if (r==GLOB_NOSPACE) {
+ throw std::bad_alloc{};
+ }
+ else if (r==0) {
+ // success
+ paths.reserve(matches.gl_pathc);
+ for (auto pathp = matches.gl_pathv; *pathp; ++pathp) {
+ paths.push_back(path{*pathp});
+ }
+ }
+
+ return paths;
+}
+
+} // namespace util
+} // namespace mc
+} // namespace nest
+
diff --git a/src/util/path.hpp b/src/util/path.hpp
index 6df1194d..d89873ad 100644
--- a/src/util/path.hpp
+++ b/src/util/path.hpp
@@ -19,6 +19,7 @@
#include <string>
#include <iostream>
+#include <utility>
#include <util/meta.hpp>
#include <util/rangeutil.hpp>
@@ -50,7 +51,7 @@ namespace posix {
// Construct or assign from value_type string or sequence.
template <typename Source>
- path(const Source& source) { assign(source); }
+ path(Source&& source) { assign(std::forward<Source>(source)); }
template <typename Iter>
path(Iter b, Iter e) { assign(b, e); }
@@ -270,6 +271,9 @@ namespace posix {
using path = posix::path;
+// POSIX glob (3) wrapper.
+std::vector<path> glob(const std::string& pattern);
+
} // namespace util
} // namespace mc
} // namespace nest
diff --git a/src/util/scope_exit.hpp b/src/util/scope_exit.hpp
new file mode 100644
index 00000000..c5bb28b4
--- /dev/null
+++ b/src/util/scope_exit.hpp
@@ -0,0 +1,54 @@
+#pragma once
+
+#include <type_traits>
+#include <utility>
+
+// Convenience class for RAII control of resources.
+
+namespace nest {
+namespace mc {
+namespace util {
+
+// `scope_exit` guard object will call provided functional object
+// on destruction. The provided functional object must be nothrow
+// move constructible.
+
+template <
+ typename F,
+ typename = typename std::enable_if<std::is_nothrow_move_constructible<F>::value>::type
+>
+class scope_exit {
+ F on_exit;
+ bool trigger = true;
+
+public:
+ template <
+ typename F2,
+ typename = typename std::enable_if<std::is_nothrow_constructible<F, F2>::value>::type
+ >
+ explicit scope_exit(F2&& f) noexcept:
+ on_exit(std::forward<F2>(f)) {}
+
+ scope_exit(scope_exit&& other) noexcept:
+ on_exit(std::move(other.on_exit))
+ {
+ other.release();
+ }
+
+ void release() noexcept {
+ trigger = false;
+ }
+
+ ~scope_exit() noexcept(noexcept(on_exit())) {
+ if (trigger) on_exit();
+ }
+};
+
+template <typename F>
+scope_exit<typename std::decay<F>::type> on_scope_exit(F&& f) {
+ return scope_exit<typename std::decay<F>::type>(std::forward<F>(f));
+}
+
+} // namespace util
+} // namespace mc
+} // namespace nest
--
GitLab