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

@@ -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});

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)

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;
 }

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;

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";
+}

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 {

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(

miniapp/morphology_pool.cpp

+#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

miniapp/morphology_pool.hpp

+#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

src/CMakeLists.txt

-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
 )

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) {

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

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;

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;
 }
 

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>

src/util/path.cpp

+#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
+

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

src/util/scope_exit.hpp

+#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