diff --git a/CMakeLists.txt b/CMakeLists.txt
index fd59431143e52b4af50be8f23b6242f1acc4100e..d181f424c8e58319da6d33471fd6020d33527df0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -275,4 +275,5 @@ add_subdirectory(mechanisms)
add_subdirectory(src)
add_subdirectory(tests)
add_subdirectory(miniapp)
+add_subdirectory(lmorpho)
diff --git a/cmake/CompilerOptions.cmake b/cmake/CompilerOptions.cmake
index 3788419abfa040010fa90fc368d9a1bab23e137b..c3982b9485fc1c99a83a83bd419c09f8417eb3eb 100644
--- a/cmake/CompilerOptions.cmake
+++ b/cmake/CompilerOptions.cmake
@@ -9,7 +9,6 @@ if(${CMAKE_CXX_COMPILER_ID} MATCHES "XL")
# Disable 'missing-braces' warning: this will inappropriately
# flag initializations such as
# std::array<int,3> a={1,2,3};
-
set(CXXOPT_WALL "${CXXOPT_WALL} -Wno-missing-braces")
# CMake, bless its soul, likes to insert this unsupported flag. Hilarity ensues.
@@ -20,26 +19,28 @@ if(${CMAKE_CXX_COMPILER_ID} MATCHES "Clang")
# Disable 'missing-braces' warning: this will inappropriately
# flag initializations such as
# std::array<int,3> a={1,2,3};
-
set(CXXOPT_WALL "${CXXOPT_WALL} -Wno-missing-braces")
endif()
if(${CMAKE_CXX_COMPILER_ID} MATCHES "GNU")
-
- # compiler flags for generating KNL-specific AVX512 instructions
- # supported in gcc 4.9.x and later
+ # Compiler flags for generating KNL-specific AVX512 instructions
+ # supported in gcc 4.9.x and later.
set(CXXOPT_KNL "-march=knl")
set(CXXOPT_AVX "-mavx")
set(CXXOPT_AVX2 "-march=core-avx2")
+
+ # Disable 'maybe-uninitialized' warning: this will be raised
+ # inappropriately in some uses of util::optional<T> when T
+ # is a primitive type.
+ set(CXXOPT_WALL "${CXXOPT_WALL} -Wno-maybe-uninitialized")
endif()
if(${CMAKE_CXX_COMPILER_ID} MATCHES "Intel")
# Disable warning for unused template parameter
- # this is raised by a templated function in the json library
-
+ # this is raised by a templated function in the json library.
set(CXXOPT_WALL "${CXXOPT_WALL} -wd488")
- # compiler flags for generating KNL-specific AVX512 instructions
+ # Compiler flags for generating KNL-specific AVX512 instructions.
set(CXXOPT_KNL "-xMIC-AVX512")
set(CXXOPT_AVX "-mavx")
set(CXXOPT_AVX2 "-march=core-avx2")
diff --git a/lmorpho/CMakeLists.txt b/lmorpho/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1bb074d891494098c73d140e4b63538e98364cc6
--- /dev/null
+++ b/lmorpho/CMakeLists.txt
@@ -0,0 +1,5 @@
+add_executable(lmorpho lmorpho.cpp lsystem.cpp lsys_models.cpp morphology.cpp morphio.cpp)
+
+target_link_libraries(lmorpho LINK_PUBLIC nestmc)
+target_link_libraries(lmorpho LINK_PUBLIC ${EXTERNAL_LIBRARIES})
+
diff --git a/lmorpho/README.md b/lmorpho/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..561bbebd863f218a44fddb4affb65203ca9d0775
--- /dev/null
+++ b/lmorpho/README.md
@@ -0,0 +1,115 @@
+# `lmporho` — generate random neuron morphologies
+
+`lmorpho` implements an L-system geometry generator for the purpose of
+generating a series of artificial neuron morphologies as described by a set of
+stasticial parameters.
+
+## Method
+
+The basic algorithm used is that of [Burke (1992)][burke1992], with extensions
+for embedding in 3-d taken from [Ascoli (2001)][ascoli2001].
+
+A dendritic tree is represented by a piece-wise linear embedding of a rooted
+tree into R³×Râº, describing the position in space and non-negative radius.
+Morphologies are represented as a sequence of these trees together with a point
+and radius describing a (spherical) soma.
+
+The generation process for a tree successively grows the terminal points by
+some length ΔL each step, changing its orientation each step according to a
+random distribution. After growing, there is a radius-dependent probability
+that it will bifurcate or terminate. Parameters describing these various
+probability distributions constitute the L-system model.
+
+The distributions used in [Ascoli (2001)][ascoli2001] can be constant (i.e. not
+random), uniform over some interval, truncated normal, or a mixture
+distribution of these. In the present version of the code, mixture
+distributions are not supported.
+
+## Output
+
+Generated morphologies can be output either in
+[SWC format](http://research.mssm.edu/cnic/swc.html), or as 'parent vectors',
+which describe the topology of the associated tree.
+
+Entry _i_ of the (zero-indexed) parent vector gives the index of the proximal
+unbranched section to which it connects. Somata have a parent index of -1.
+
+Morphology information can be written one each to separate files or
+concatenated. If concatenated, the parent vector indices will be shifted as
+required to maintain consistency.
+
+## Models
+
+Two L-system parameter sets are provided as 'built-in'; if neither model is
+selected, a simple, non-biological default model is used.
+
+As mixture distributions are not yet supported in `lmorpho`, these models are
+simplified approximations to those in the literature as described below.
+
+### Motor neuron model
+
+The 'motoneuron' model corresponds to the adult cat spinal alpha-motoneuron
+model described in [Ascoli (2001)][ascoli2001], based in turn on the model of
+[Burke (1992)][burke1992].
+
+The model parameters in Ascoli are not complete; missing parameters were taken
+from the corresponding ArborVitae model parameters in the same paper, and
+somatic diameters were taken from [Cullheim (1987)][cullheim1987].
+
+### Purkinke neuron model
+
+The 'purkinke' model corresponds to the guinea pig Purkinje cell model from
+[Ascoli (2001)][ascoli2001], which was fit to data from [Rapp
+(1994)][rapp1994]. Somatic diameters are a simple fit to the three measurements
+in the original source.
+
+Some required parameters are not recorded in [Ascoli (2001)][ascoli2001];
+the correlation component `diam_child_a` and the `length_step` ΔL are
+taken from the corresponding L-Neuron parameter description in the
+[L-Neuron database](l-neuron-database). These should be verified by
+fitting against the digitized morphologies as found in the database.
+
+Produced neurons from this model do not look especially realistic.
+
+### Other models
+
+There is not yet support for loading in arbitrary parameter sets, or for
+translating or approximating the parameters that can be found in the
+[L-Neuron database][l-neuron-database]. Support for parameter estimation
+from a population of discretized morphologies would be useful.
+
+## References
+
+1. Ascoli, G. A., Krichmar, J. L., Scorcioni, R., Nasuto, S. J., Senft, S. L.
+ and Krichmar, G. L. (2001). Computer generation and quantitative morphometric
+ analysis of virtual neurons. _Anatomy and Embryology_ _204_(4), 283–301.
+ [DOI: 10.1007/s004290100201][ascoli2001]
+
+2. Burke, R. E., Marks, W. B. and Ulfhake, B. (1992). A parsimonious description
+ of motoneuron dendritic morphology using computer simulation.
+ _The Journal of Neuroscience_ _12_(6), 2403–2416. [online][burke1992]
+
+3. Cullheim, S., Fleshman, J. W., Glenn, L. L., and Burke, R. E. (1987).
+ Membrane area and dendritic structure in type-identified triceps surae alpha
+ motoneurons. _The Journal of Comparitive Neurology_ _255_(1), 68–81.
+ [DOI: 10.1002/cne.902550106][cullheim1987]
+
+4. Rapp, M., Segev, I. and Yaom, Y. (1994). Physiology, morphology and detailed
+ passive models of guinea-pig cerebellar Purkinje cells.
+ _The Journal of Physiology_, _474_, 101–118.
+ [DOI: 10.1113/jphysiol.1994.sp020006][rapp1994].
+
+
+[ascoli2001]: http://dx.doi.org/10.1007/s004290100201
+ "Ascoli et al. (2001). Computer generation […] of virtual neurons."
+
+[burke1992]: http://www.jneurosci.org/content/12/6/2403
+ "Burke et al. (1992). A parsimonious description of motoneuron dendritic morphology […]"
+
+[cullheim1987]: http://dx.doi.org/10.1002/cne.902550106
+ "Cullheim et al. (1987). Membrane area and dendritic structure in […] alpha motoneurons."
+
+[rapp1994]: http://dx.doi.org/10.1113/jphysiol.1994.sp020006
+ "Rapp et al. (1994). Physiology, morphology […] of guinea-pig cerebellar Purkinje cells."
+
+[l-neuron-database]: http://krasnow1.gmu.edu/cn3/L-Neuron/database/index.html
diff --git a/lmorpho/lmorpho.cpp b/lmorpho/lmorpho.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..65f32009b766b84209cef6954c52c848902a6bf1
--- /dev/null
+++ b/lmorpho/lmorpho.cpp
@@ -0,0 +1,115 @@
+#include <algorithm>
+#include <iostream>
+#include <fstream>
+#include <map>
+#include <sstream>
+#include <vector>
+
+#include <tinyopt.hpp>
+#include <util/optional.hpp>
+
+#include "morphology.hpp"
+#include "morphio.hpp"
+#include "lsystem.hpp"
+#include "lsys_models.hpp"
+
+namespace to = nest::mc::to;
+using nest::mc::util::optional;
+using nest::mc::util::nothing;
+using nest::mc::util::just;
+
+const char* usage_str =
+"[OPTION]...\n"
+"\n"
+" -n, --count=N Number of morphologies to generate.\n"
+" -m, --model=MODEL Use L-system MODEL for generation (see below).\n"
+" -g, --segment=DX Segment model into compartments of max size DX µm.\n"
+" --swc=FILE Output morphologies as SWC to FILE (see below).\n"
+" --pvec=FILE Output 'parent vector' structural representation\n"
+" to FILE.\n"
+" -h, --help Emit this message and exit.\n"
+"\n"
+"Generate artificial neuron morphologies based on L-system descriptions.\n"
+"\n"
+"If a FILE arrgument contains a '%', then one file will be written for\n"
+"each generated morphology, with the '%' replaced by the index of the\n"
+"morphology, starting from zero. Output for each morphology will otherwise\n"
+"be concatenated: SWC files will be headed by a comment line with the\n"
+"index of the morphology; parent vectors will be merged into one long\n"
+"vector. A FILE argument of '-' corresponds to standard output.\n"
+"\n"
+"Currently supported MODELs:\n"
+" motoneuron Adult cat spinal alpha-motoneurons, based on models\n"
+" and data in Burke 1992 and Ascoli 2001.\n"
+" purkinje Guinea pig Purkinje cells, basd on models and data\n"
+" from Rapp 1994 and Ascoli 2001.\n";
+
+int main(int argc, char** argv) {
+ // options
+ int n_morph = 1;
+ optional<unsigned> rng_seed;
+ optional<std::string> swc_file;
+ optional<std::string> pvector_file;
+ double segment_dx = 0;
+
+ std::pair<const char*, const lsys_param*> models[] = {
+ {"motoneuron", &alpha_motoneuron_lsys},
+ {"purkinje", &purkinje_lsys}
+ };
+ lsys_param P;
+
+ try {
+ auto arg = argv+1;
+ while (*arg) {
+ if (auto o = to::parse_opt<int>(arg, 'n', "count")) {
+ n_morph = *o;
+ }
+ if (auto o = to::parse_opt<int>(arg, 's', "seed")) {
+ rng_seed = *o;
+ }
+ else if (auto o = to::parse_opt<std::string>(arg, 0, "swc")) {
+ swc_file = *o;
+ }
+ else if (auto o = to::parse_opt<double>(arg, 'g', "segment")) {
+ segment_dx = *o;
+ }
+ else if (auto o = to::parse_opt<std::string>(arg, 'p', "pvec")) {
+ pvector_file = *o;
+ }
+ else if (auto o = to::parse_opt<const lsys_param*>(arg, 'm', "model", to::keywords(models))) {
+ P = **o;
+ }
+ else if (to::parse_opt(arg, 'h', "help")) {
+ std::cout << "Usage: " << argv[0] << " " << usage_str;
+ return 0;
+ }
+ else {
+ throw to::parse_opt_error(*arg, "unrecognized option");
+ }
+ }
+
+ std::minstd_rand g;
+ if (rng_seed) g.seed(rng_seed.get());
+
+ auto emit_swc = swc_file? just(swc_emitter(*swc_file, n_morph)): nothing;
+ auto emit_pvec = pvector_file? just(pvector_emitter(*pvector_file, n_morph)): nothing;
+
+ for (int i=0; i<n_morph; ++i) {
+ auto morph = generate_morphology(P, g);
+ morph.segment(segment_dx);
+
+ if (emit_swc) (*emit_swc)(i, morph);
+ if (emit_pvec) (*emit_pvec)(i, morph);
+ }
+ }
+ catch (to::parse_opt_error& e) {
+ std::cerr << argv[0] << ": " << e.what() << "\n";
+ std::cerr << "Try '" << argv[0] << " --help' for more information.\n";
+ std::exit(2);
+ }
+ catch (std::exception& e) {
+ std::cerr << "caught exception: " << e.what() << "\n";
+ std::exit(1);
+ }
+}
+
diff --git a/lmorpho/lsys_models.cpp b/lmorpho/lsys_models.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c26437d4a14b638762f4665d31d718563a89ff74
--- /dev/null
+++ b/lmorpho/lsys_models.cpp
@@ -0,0 +1,131 @@
+#include <math.hpp>
+
+#include "lsystem.hpp"
+#include "lsys_models.hpp"
+
+static constexpr double inf = nest::mc::math::infinity<double>();
+
+// Predefined parameters for two classes of neurons. Numbers taken primarily
+// from Ascoli et al. 2001, but some details (soma diameters for example)
+// taken from earlier literature.
+//
+// Without mixture distribution support, these aren't entirely faithful to
+// the original sources.
+//
+// Refer to `README.md` for references and more information.
+
+lsys_param make_alpha_motoneuron_lsys() {
+ lsys_param L;
+
+ // Soma diameter [µm].
+ L.diam_soma = { 47.0, 65.5, 57.6, 5.0 };
+
+ // Number of dendritic trees (rounded to nearest int).
+ L.n_tree = { 8.0, 16.0 };
+
+ // Initial dendrite diameter [µm]. (Dstem)
+ L.diam_initial = { 3.0, 12.0 };
+
+ // Dendrite step length [µm]. (ΔL)
+ L.length_step = { 25 };
+
+ // Initial roll (intrinsic rotation about x-axis) [degrees].
+ L.roll_initial = { -180.0, 180.0 };
+
+ // Initial pitch (intrinsic rotation about y-axis) [degrees].
+ L.pitch_initial = { 0.0, 180.0 };
+
+ // Tortuousness: roll within section over ΔL [degrees].
+ L.roll_section = { -180.0, 180.0 };
+
+ // Tortuousness: pitch within section over ΔL [degrees].
+ L.pitch_section = { -15.0, 15.0, 0.0, 5.0 };
+
+ // Taper rate: diameter decrease per unit length.
+ L.taper = { -1.25e-3 };
+
+ // Branching torsion: roll at branch point [degrees].
+ L.roll_at_branch = { 0.0, 180.0 };
+
+ // Branch angle between siblings [degrees].
+ L.branch_angle = { 1.0, 172.0, 45.0, 20.0 };
+
+ // Correlated child branch radius distribution parameters.
+ L.diam_child_a = -0.2087;
+ L.diam_child_r = { 0.2, 1.8, 0.8255, 0.2125 };
+
+ // Termination probability parameters [1/µm].
+ L.pterm_k1 = 2.62e-2;
+ L.pterm_k2 = -2.955;
+
+ // Bifurcation probability parameters [1/μm].
+ L.pbranch_ov_k1 = 2.58e-5;
+ L.pbranch_ov_k2 = 2.219;
+ L.pbranch_nov_k1 = 2.34e-3;
+ L.pbranch_nov_k2 = 0.194;
+
+ // Absolute maximum dendritic extent [µm].
+ L.max_extent = 2000;
+
+ return L;
+}
+
+lsys_param alpha_motoneuron_lsys = make_alpha_motoneuron_lsys();
+
+// Some parameters missing from the literature are taken from
+// the `purkburk.prm` L-Neuron parameter file:
+// http://krasnow1.gmu.edu/cn3/L-Neuron/database/purk/bu/purkburk.prm
+// Note 'Pov' and 'Pterm' numbers are swapped in Ascoli fig 3B.
+
+lsys_param make_purkinje_lsys() {
+ lsys_param L;
+
+ // Soma diameter [µm].
+ // (Gaussian fit to 3 samples rom Rapp 1994); truncate to 2σ.
+ L.diam_soma = { 22.12, 27.68, 24.9, 1.39 };
+
+ // Number of dendritic trees (rounded to nearest int).
+ L.n_tree = { 1 };
+
+ // Initial dendrite diameter [µm]. (Dstem)
+ L.diam_initial = { 4.8, 7.6, 6.167, 1.069 };
+
+ // Dendrite step length [µm]. (ΔL)
+ L.length_step = { 2.3 }; // from `purkburk.prm`
+
+ // Tortuousness: roll within section over ΔL [degrees].
+ L.roll_section = { 0 };
+
+ // Tortuousness: pitch within section over ΔL [degrees].
+ L.pitch_section = { -5, 5 }; // from `purkburk.prm`
+
+ // Taper rate: diameter decrease per unit length.
+ L.taper = { -inf, inf, -0.010, 0.022 };
+
+ // Branching torsion: roll at branch point [degrees].
+ L.roll_at_branch = { -inf, inf, 0.0, 6.5 }; // from `purkburj.prm`
+
+ // Branch angle between siblings [degrees].
+ L.branch_angle = { 1.0, 179.0, 71.0, 40.0 };
+
+ // Correlated child branch radius distribution parameters.
+ L.diam_child_a = 5.47078; // from `purkburk.prm`
+ L.diam_child_r = { -inf, inf, 0.112, 0.035 };
+
+ // Termination probability parameters [1/µm].
+ L.pterm_k1 = 0.1973; // from `purkburj.prm`; 0.4539 in Ascoli.
+ L.pterm_k2 = -1.1533;
+
+ // Bifurcation probability parameters [1/μm].
+ L.pbranch_ov_k1 = 0.1021; // from `purkburk.prm`; 0.0355 in Ascoli.
+ L.pbranch_ov_k2 = 0.7237;
+ L.pbranch_nov_k1 = 0.1021; // from `purkburk.prm`; 0.2349 in Ascoli.
+ L.pbranch_nov_k2 = -0.0116;
+
+ // Absolute maximum dendritic extent [µm].
+ L.max_extent = 20000;
+
+ return L;
+}
+
+lsys_param purkinje_lsys = make_purkinje_lsys();
diff --git a/lmorpho/lsys_models.hpp b/lmorpho/lsys_models.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..96638e5ad0322cd830b4bf283a15fee152650592
--- /dev/null
+++ b/lmorpho/lsys_models.hpp
@@ -0,0 +1,9 @@
+#pragma once
+
+#include "lsystem.hpp"
+
+// Predefined parameters for two classes of neurons.
+// See lsys_models.cc for details.
+
+extern lsys_param alpha_motoneuron_lsys;
+extern lsys_param purkinje_lsys;
diff --git a/lmorpho/lsystem.cpp b/lmorpho/lsystem.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7e46c913b1c02df10d7e684a0e4d22858f50b8e5
--- /dev/null
+++ b/lmorpho/lsystem.cpp
@@ -0,0 +1,309 @@
+#include <cmath>
+#include <iostream>
+#include <random>
+#include <stack>
+#include <vector>
+
+#include <math.hpp>
+
+#include "lsystem.hpp"
+#include "morphology.hpp"
+
+using namespace nest::mc::math;
+
+// L-system implementation.
+
+// Random distribution used in the specification of L-system parameters.
+// It can be a constant, uniform over an interval, or truncated normal.
+// Note that mixture distributions of the above should be, but isn't yet,
+// implemented.
+
+struct lsys_distribution {
+ using result_type = double;
+ using param_type = lsys_distribution_param;
+
+ lsys_distribution(const param_type& p): p_(p) { init(); }
+
+ lsys_distribution(): p_() { init(); }
+ lsys_distribution(result_type x): p_(x) { init(); }
+ lsys_distribution(result_type lb, result_type ub): p_(lb, ub) { init(); }
+ lsys_distribution(result_type lb, result_type ub, result_type mu, result_type sigma):
+ p_(lb, ub, mu, sigma)
+ {
+ init();
+ }
+
+ param_type param() const { return p_; }
+ void param(const param_type& pbis) { p_=pbis; }
+
+ result_type min() const { return p_.lb; }
+ result_type max() const {
+ return p_.kind==param_type::constant? min(): p_.ub;
+ }
+
+ bool operator==(const lsys_distribution& y) const {
+ return p_==y.p_;
+ }
+
+ bool operator!=(const lsys_distribution& y) const {
+ return p_!=y.p_;
+ }
+
+ // omitting out ostream, istream methods
+ // (...)
+
+ template <typename Gen>
+ result_type operator()(Gen& g, const param_type &p) const {
+ lsys_distribution dist(p);
+ return dist(g);
+ }
+
+ template <typename Gen>
+ result_type operator()(Gen& g) const {
+ switch (p_.kind) {
+ case param_type::constant:
+ return p_.lb;
+ case param_type::uniform:
+ return uniform_dist(g);
+ case param_type::normal:
+ // TODO: replace with a robust truncated normal generator!
+ for (;;) {
+ result_type r = normal_dist(g);
+ if (r>=p_.lb && r<=p_.ub) return r;
+ }
+ }
+ return 0;
+ }
+
+private:
+ param_type p_;
+ mutable std::uniform_real_distribution<double> uniform_dist;
+ mutable std::normal_distribution<double> normal_dist;
+
+ void init() {
+ switch (p_.kind) {
+ case param_type::uniform:
+ uniform_dist = std::uniform_real_distribution<double>(p_.lb, p_.ub);
+ break;
+ case param_type::normal:
+ normal_dist = std::normal_distribution<result_type>(p_.mu, p_.sigma);
+ break;
+ default: ;
+ }
+ }
+};
+
+class burke_correlated_r {
+ double a_;
+ lsys_distribution r_;
+
+public:
+ struct result_type { double rho1, rho2; };
+
+ burke_correlated_r(double a, lsys_distribution_param r): a_(a), r_(r) {}
+
+ template <typename Gen>
+ result_type operator()(Gen& g) const {
+ double r1 = r_(g);
+ double r2 = r_(g);
+ return { r1+a_*r2, r2+a_*r1 };
+ }
+};
+
+class burke_exp_test {
+ double k1_;
+ double k2_;
+
+public:
+ using result_type = int;
+
+ burke_exp_test(double k1, double k2): k1_(k1), k2_(k2) {}
+
+ template <typename Gen>
+ result_type operator()(double delta_l, double radius, Gen& g) const {
+ std::bernoulli_distribution B(delta_l*k1_*std::exp(k2_*radius*2.0));
+ return B(g);
+ }
+};
+
+class burke_exp2_test {
+ double ak1_;
+ double ak2_;
+ double bk1_;
+ double bk2_;
+
+public:
+ using result_type = int;
+
+ burke_exp2_test(double ak1, double ak2, double bk1, double bk2):
+ ak1_(ak1), ak2_(ak2), bk1_(bk1), bk2_(bk2) {}
+
+ template <typename Gen>
+ result_type operator()(double delta_l, double radius, Gen& g) const {
+ double diam = 2.*radius;
+ std::bernoulli_distribution B(
+ delta_l*std::min(ak1_*std::exp(ak2_*diam), bk1_*std::exp(bk2_*diam)));
+
+ return B(g);
+ }
+};
+
+// L-system parameter set with instantiated distributions.
+struct lsys_sampler {
+ // Create distributions once from supplied parameter set.
+ explicit lsys_sampler(const lsys_param& P):
+ diam_soma(P.diam_soma),
+ n_tree(P.n_tree),
+ diam_initial(P.diam_initial),
+ length_step(P.length_step),
+ roll_initial(P.roll_initial),
+ pitch_initial(P.pitch_initial),
+ roll_section(P.roll_section),
+ pitch_section(P.pitch_section),
+ taper(P.taper),
+ roll_at_branch(P.roll_at_branch),
+ branch_angle(P.branch_angle),
+ diam_child(P.diam_child_a, P.diam_child_r),
+ term_test(P.pterm_k1, P.pterm_k2),
+ branch_test(P.pbranch_ov_k1, P.pbranch_ov_k2,
+ P.pbranch_nov_k1, P.pbranch_nov_k2),
+ max_extent(P.max_extent)
+ {}
+
+ lsys_distribution diam_soma;
+ lsys_distribution n_tree;
+ lsys_distribution diam_initial;
+ lsys_distribution length_step;
+ lsys_distribution roll_initial;
+ lsys_distribution pitch_initial;
+ lsys_distribution roll_section;
+ lsys_distribution pitch_section;
+ lsys_distribution taper;
+ lsys_distribution roll_at_branch;
+ lsys_distribution branch_angle;
+ burke_correlated_r diam_child;
+ burke_exp_test term_test;
+ burke_exp2_test branch_test;
+ double max_extent;
+};
+
+struct section_tip {
+ section_point p = {0., 0., 0., 0.};
+ quaternion rotation = {1., 0., 0., 0.};
+ double somatic_distance = 0.;
+};
+
+struct grow_result {
+ std::vector<section_point> points;
+ std::vector<section_tip> children;
+ double length = 0.;
+};
+
+constexpr double deg_to_rad = 3.1415926535897932384626433832795l/180.0;
+
+template <typename Gen>
+grow_result grow(section_tip tip, const lsys_sampler& S, Gen &g) {
+ constexpr quaternion xaxis = {0, 1, 0, 0};
+ static std::uniform_real_distribution<double> U;
+
+ grow_result result;
+ std::vector<section_point>& points = result.points;
+
+ points.push_back(tip.p);
+ for (;;) {
+ quaternion step = xaxis^tip.rotation;
+ double dl = S.length_step(g);
+ tip.p.x += dl*step.x;
+ tip.p.y += dl*step.y;
+ tip.p.z += dl*step.z;
+ tip.p.r += dl*0.5*S.taper(g);
+ tip.somatic_distance += dl;
+ result.length += dl;
+
+ if (tip.p.r<0) tip.p.r = 0;
+
+ double phi = S.roll_section(g);
+ double theta = S.pitch_section(g);
+ tip.rotation *= rotation_x(deg_to_rad*phi);
+ tip.rotation *= rotation_y(deg_to_rad*theta);
+
+ points.push_back(tip.p);
+
+ if (tip.p.r==0 || tip.somatic_distance>=S.max_extent) {
+ return result;
+ }
+
+ if (S.branch_test(dl, tip.p.r, g)) {
+ double branch_phi = S.roll_at_branch(g);
+ double branch_angle = S.branch_angle(g);
+ double branch_theta1 = U(g)*branch_angle;
+ double branch_theta2 = branch_theta1-branch_angle;
+ auto r = S.diam_child(g);
+
+ tip.rotation *= rotation_x(deg_to_rad*branch_phi);
+
+ section_tip t1 = tip;
+ t1.p.r = t1.p.r * r.rho1;
+ t1.rotation *= rotation_y(deg_to_rad*branch_theta1);
+
+ section_tip t2 = tip;
+ t2.p.r = t2.p.r * r.rho2;
+ t2.rotation *= rotation_y(deg_to_rad*branch_theta2);
+
+ result.children = {t1, t2};
+ return result;
+ }
+
+ if (S.term_test(dl, tip.p.r, g)) {
+ return result;
+ }
+ }
+}
+
+morphology generate_morphology(const lsys_param& P, lsys_generator &g) {
+ constexpr quaternion xaxis = {0, 1, 0, 0};
+ morphology morph;
+
+ lsys_sampler S(P);
+ double soma_radius = 0.5*S.diam_soma(g);
+ morph.soma = {0, 0, 0, soma_radius};
+
+ struct section_start {
+ section_tip tip;
+ unsigned parent_id;
+ };
+ std::stack<section_start> starts;
+ unsigned next_id = 1u;
+
+ int n = (int)std::round(S.n_tree(g));
+ for (int i=0; i<n; ++i) {
+ double phi = S.roll_initial(g);
+ double theta = S.pitch_initial(g);
+ double radius = 0.5*S.diam_initial(g);
+
+ section_tip tip;
+ tip.rotation = rotation_x(deg_to_rad*phi)*rotation_y(deg_to_rad*theta);
+ tip.somatic_distance = 0.0;
+
+ auto p = (soma_radius*xaxis)^tip.rotation;
+ tip.p = {p.x, p.y, p.z, radius};
+
+ starts.push({tip, 0u});
+ }
+
+ while (!starts.empty() && morph.sections.size()<P.max_sections) {
+ auto start = starts.top();
+ starts.pop();
+
+ auto branch = grow(start.tip, S, g);
+ section_geometry section = {next_id++, start.parent_id, branch.children.empty(), std::move(branch.points), branch.length};
+
+ for (auto child: branch.children) {
+ starts.push({child, section.id});
+ }
+ morph.sections.push_back(std::move(section));
+ }
+
+ return morph;
+}
+
diff --git a/lmorpho/lsystem.hpp b/lmorpho/lsystem.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..f3ce4deb3ca578df970153800da40c4cdd436acc
--- /dev/null
+++ b/lmorpho/lsystem.hpp
@@ -0,0 +1,131 @@
+#pragma once
+
+#include <random>
+
+#include "morphology.hpp"
+
+struct lsys_param;
+
+using lsys_generator = std::minstd_rand;
+
+morphology generate_morphology(const lsys_param& P, lsys_generator& g);
+
+// The distribution parameters used in the specification of the L-system parameters.
+// The distribution can be a constant, uniform over an interval, or truncated normal.
+// Note that mixture distributions of the above should be, but isn't yet,
+// implemented.
+
+struct lsys_distribution_param {
+ using result_type = double;
+ enum lsys_kind {
+ constant, uniform, normal
+ };
+
+ lsys_kind kind;
+ result_type lb; // lower bound
+ result_type ub; // upper bound (used for only non-constant)
+ result_type mu; // mean (used only for truncated Gaussian)
+ result_type sigma; // s.d. (used only for truncated Gaussian)
+
+ bool operator==(const lsys_distribution_param &p) const {
+ if (kind!=p.kind) return false;
+ switch (kind) {
+ case constant:
+ return lb == p.lb;
+ case uniform:
+ return lb == p.lb && ub == p.ub;
+ case normal:
+ return lb == p.lb && ub == p.ub && mu == p.mu && sigma == p.sigma;
+ }
+ return false;
+ }
+
+ bool operator!=(const lsys_distribution_param &p) const {
+ return !(*this==p);
+ }
+
+ // zero or one parameter => constant
+ lsys_distribution_param():
+ kind(constant), lb(0.0) {}
+
+ lsys_distribution_param(result_type x):
+ kind(constant), lb(x) {}
+
+ // two paramter => uniform
+ lsys_distribution_param(result_type lb, result_type ub):
+ kind(uniform), lb(lb), ub(ub) {}
+
+ // four paramter => truncated normal
+ lsys_distribution_param(result_type lb, result_type ub, result_type mu, result_type sigma):
+ kind(normal), lb(lb), ub(ub), mu(mu), sigma(sigma) {}
+};
+
+
+// Parameters as referenced in Ascoli 2001 are given in parentheses below.
+// Defaults below are chosen to be safe, boring and not biological.
+struct lsys_param {
+ // Soma diameter [µm].
+ lsys_distribution_param diam_soma = { 20 };
+
+ // Number of dendritic trees (rounded to nearest int). (Ntree)
+ lsys_distribution_param n_tree = { 1 };
+
+ // Initial dendrite diameter [µm]. (Dstem)
+ lsys_distribution_param diam_initial = { 2 };
+
+ // Dendrite step length [µm]. (ΔL)
+ lsys_distribution_param length_step = { 25 };
+
+ // Dendrites grow along x-axis after application of (instrinsic) rotations;
+ // roll (about x-axis) is applied first, followed by pitch (about y-axis).
+
+ // Initial roll (intrinsic rotation about x-axis) [degrees]. (Taz)
+ lsys_distribution_param roll_initial = { 0 };
+
+ // Initial pitch (intrinsic rotation about y-axis) [degrees]. (Tel)
+ lsys_distribution_param pitch_initial = { 0 };
+
+ // Tortuousness: roll within section over ΔL [degrees]. (Eaz)
+ lsys_distribution_param roll_section = { 0 };
+
+ // Tortuousness: pitch within section over ΔL [degrees]. (Eel)
+ lsys_distribution_param pitch_section = { 0 };
+
+ // Taper rate: diameter decrease per unit length. (TPRB)
+ lsys_distribution_param taper = { 0 };
+
+ // Branching torsion: roll at branch point [degrees]. (Btor)
+ lsys_distribution_param roll_at_branch = { 0 };
+
+ // Branch angle between siblings [degrees]. (Bamp)
+ lsys_distribution_param branch_angle = { 45 };
+
+ // Child branch diameter ratios are typically anticorrelated.
+ // The Burke 1992 parameterization describes the two child ratios
+ // as d1 = r1 + a*r2 and d2 = r2 + a*r1 where a is a constant
+ // determined by the correlation, and r1 and r2 are drawn from
+ // a common distribution (`d_child_r` below).
+ double diam_child_a = 0;
+ lsys_distribution_param diam_child_r = { 0.5 };
+
+ // Bifurcation and termination probabilities per unit-length below are given by
+ // P = k1 * exp(k2*diameter), described by the paremeters k1 and k2 [1/µm].
+
+ // Termination probability parameters [1/µm]. (k1trm, k2trm)
+ double pterm_k1 = 0.05;
+ double pterm_k2 = -2;
+
+ // Bifurcation probability is given by minimum of two probabilities, `pbranch_ov' and
+ // `pbranch_nov' below.
+ double pbranch_ov_k1 = 0.01;
+ double pbranch_ov_k2 = 0;
+ double pbranch_nov_k1 = 0.01;
+ double pbranch_nov_k2 = 0;
+
+ // Absolute maximum dendritic extent [µm]. (Forces termination of algorithm)
+ double max_extent = 2000;
+
+ // Absolute maximum number of unbranched sections. (Forces termination of algorithm)
+ unsigned max_sections = 10000;
+};
+
diff --git a/lmorpho/morphio.cpp b/lmorpho/morphio.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e6259be54a8d68b12f13f069b337057836c29f76
--- /dev/null
+++ b/lmorpho/morphio.cpp
@@ -0,0 +1,179 @@
+#include <fstream>
+#include <iterator>
+#include <map>
+#include <string>
+#include <vector>
+
+#include <swcio.hpp>
+
+#include "morphio.hpp"
+
+using nest::mc::io::swc_record;
+
+std::vector<swc_record> as_swc(const morphology& morph);
+
+// printf wrappers.
+
+template <typename... Args>
+static std::string strprintf(const char* fmt, Args&&... args) {
+ thread_local static std::vector<char> buffer(1024);
+
+ for (;;) {
+ int n = std::snprintf(buffer.data(), buffer.size(), fmt, std::forward<Args>(args)...);
+ if (n<0) return ""; // error
+ if ((unsigned)n<buffer.size()) return std::string(buffer.data());
+ buffer.resize(2*n);
+ }
+}
+
+template <typename... Args>
+static std::string strprintf(std::string fmt, Args&&... args) {
+ return strprintf(fmt.c_str(), std::forward<Args>(args)...);
+}
+
+// Multi-file manager implementation.
+multi_file::multi_file(const std::string& pattern, int digits) {
+ auto npos = std::string::npos;
+
+ file_.exceptions(std::ofstream::failbit);
+ concat_ = (pattern.find("%")==npos);
+ use_stdout_ = pattern.empty() || pattern=="-";
+
+ if (!concat_) {
+ std::string nfmt = digits? "%0"+std::to_string(digits)+"d": "%d";
+ std::string::size_type i = 0;
+ for (;;) {
+ auto p = pattern.find("%", i);
+
+ if (p==npos) {
+ fmt_ += pattern.substr(i);
+ break;
+ }
+ else {
+ fmt_ += pattern.substr(i, p-i);
+ fmt_ += i==0? nfmt: "%%";
+ i = p+1;
+ }
+ }
+ }
+ else {
+ filename_ = pattern;
+ }
+}
+
+void multi_file::open(unsigned n) {
+ if (use_stdout_ || (file_.is_open() && (concat_ || n==current_n_))) {
+ return;
+ }
+
+ if (file_.is_open()) file_.close();
+
+ std::string fname = concat_? filename_: strprintf(fmt_, n);
+ file_.open(fname);
+
+ current_n_ = n;
+}
+
+// SWC transform
+
+using nest::mc::io::swc_record;
+
+std::vector<swc_record> as_swc(const morphology& morph) {
+ using kind = swc_record::kind;
+ std::map<int, int> parent_end_id;
+ std::vector<swc_record> swc;
+
+ // soma
+ const auto &p = morph.soma;
+ int id = 0;
+ parent_end_id[0] = 0;
+ swc.emplace_back(kind::soma, id, p.x, p.y, p.z, p.r, -1);
+
+ // dendrites:
+ for (auto& sec: morph.sections) {
+ int parent = parent_end_id[sec.parent_id];
+
+ const auto& points = sec.points;
+ auto n = points.size();
+ if (n<2) {
+ throw std::runtime_error(strprintf("surprisingly short cable: id=%d, size=%ul", sec.id, n));
+ }
+
+ // Include first point only for dendrites segments attached to soma.
+ if (sec.parent_id==0) {
+ const auto& p = points[0];
+ ++id;
+ swc.emplace_back(kind::fork_point, id, p.x, p.y, p.z, p.r, parent);
+ parent = id;
+ }
+
+ // Interior points.
+ for (unsigned i = 1; i<n-1; ++i) {
+ const auto& p = points[i];
+ ++id;
+ swc.emplace_back(kind::dendrite, id, p.x, p.y, p.z, p.r, parent);
+ parent = id;
+ }
+
+ // Final point (fork or terminal).
+ const auto& p = points.back();
+ ++id;
+ swc.emplace_back(sec.terminal? kind::end_point: kind::fork_point, id, p.x, p.y, p.z, p.r, parent);
+ parent_end_id[sec.id] = id;
+ }
+
+ return swc;
+}
+
+// SWC emitter implementation.
+
+void swc_emitter::operator()(unsigned index, const morphology& m) {
+ file_.open(index);
+ auto& stream = file_.stream();
+
+ auto swc = as_swc(m);
+ stream << "# lmorpho generated morphology\n# index: " << index << "\n";
+ std::copy(swc.begin(), swc.end(), std::ostream_iterator<swc_record>(stream, "\n"));
+}
+
+// pvector emitter implementation.
+
+std::vector<int> as_pvector(const morphology& morph, unsigned offset) {
+ std::map<int, unsigned> parent_index; // section id to segment index
+ std::vector<int> pvec;
+ unsigned index = offset; // starting segment index
+
+ // soma
+ parent_index[0] = index;
+ pvec.push_back(-1);
+ ++index;
+
+ // dendrites:
+ for (auto& sec: morph.sections) {
+ int parent = parent_index[sec.parent_id];
+
+ auto n = sec.points.size();
+ if (n<2) {
+ throw std::runtime_error(strprintf("surprisingly short cable: id=%d, size=%ul", sec.id, n));
+ }
+
+ for (unsigned i = 1; i<n; ++i) {
+ pvec.push_back(parent);
+ parent = index++;
+ }
+
+ parent_index[sec.id] = parent;
+ }
+
+ return pvec;
+}
+
+void pvector_emitter::operator()(unsigned index, const morphology& m) {
+ auto pvec = as_pvector(m, offset_);
+ if (coalesce_) offset_ += pvec.size();
+
+ file_.open(index);
+ auto& stream = file_.stream();
+ std::copy(pvec.begin(), pvec.end(), std::ostream_iterator<int>(stream, "\n"));
+}
+
diff --git a/lmorpho/morphio.hpp b/lmorpho/morphio.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..37ec56ae69e06caed6b9bb96078a94f41f20df8a
--- /dev/null
+++ b/lmorpho/morphio.hpp
@@ -0,0 +1,93 @@
+#pragma once
+
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+#include "morphology.hpp"
+
+// Manage access to a single file, std::cout, or an indexed
+// sequence of files.
+
+class multi_file {
+private:
+ std::ofstream file_;
+ bool concat_ = false;
+ bool use_stdout_ = false;
+ std::string fmt_; // use if not concat_
+ std::string filename_; // use if concat_
+ unsigned current_n_ = 0;
+
+public:
+ // If pattern is '-' or empty, represent std::cout.
+ // If pattern contains '%', use a different file for each index.
+ // Otherwise, use a single file given by pattern.
+ explicit multi_file(const std::string& pattern, int digits=0);
+
+ multi_file(multi_file&&) = default;
+
+ // Open file for index n, closing previously open file if different.
+ void open(unsigned n);
+
+ // Close currently open file.
+ void close() { if (!use_stdout_ && file_.is_open()) file_.close(); }
+
+ // True if writing to a single file or std::cout.
+ bool single_file() const { return concat_; }
+
+ // Return std::cout or file stream as appropriate.
+ std::ostream& stream() { return use_stdout_? std::cout: file_; }
+};
+
+static constexpr unsigned digits(unsigned n) {
+ return n? 1+digits(n/10): 0;
+}
+
+// Write a sequence of morphologies to one or more files
+// as given my `pattern`.
+
+class swc_emitter {
+ multi_file file_;
+
+public:
+ // `pattern` is as for `multi_file`; number `n` optionally
+ // specifies the total number of morphologies, for better
+ // numeric formatting.
+ explicit swc_emitter(std::string pattern, unsigned n=0):
+ file_(pattern, digits(n)) {}
+
+ swc_emitter(swc_emitter&&) = default;
+
+ // write `index`th morphology as SWC.
+ void operator()(unsigned index, const morphology& m);
+
+ void close() { file_.close(); }
+ ~swc_emitter() { close(); }
+};
+
+// Write pvectors for a sequence of morphologies to one or more files
+// as given my `pattern`. Coalesce the vectors if writing to a single
+// file/stream.
+
+class pvector_emitter {
+ multi_file file_;
+ bool coalesce_ = false;
+ unsigned offset_ = 0;
+
+public:
+ // `pattern` is as for `multi_file`; number `n` optionally
+ // specifies the total number of morphologies, for better
+ // numeric formatting.
+ explicit pvector_emitter(std::string pattern, unsigned n=0):
+ file_(pattern, digits(n)),
+ coalesce_(file_.single_file()) {}
+
+ pvector_emitter(pvector_emitter&&) = default;
+
+ // write pvector for `index`th morphology.
+ void operator()(unsigned index, const morphology& m);
+
+ void close() { file_.close(); }
+ ~pvector_emitter() { close(); }
+};
+
diff --git a/lmorpho/morphology.cpp b/lmorpho/morphology.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7e9c60ea1c69e98229fc465ce6bf7c4ecfb57b90
--- /dev/null
+++ b/lmorpho/morphology.cpp
@@ -0,0 +1,62 @@
+#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.
+ for (unsigned i = 1; i<npoint;) {
+ 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]);
+ x = k*length/nseg;
+ }
+ else {
+ ++i;
+ 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
new file mode 100644
index 0000000000000000000000000000000000000000..d1adaa6fdf6a3c0966cb827a0f4894d0db8f3b51
--- /dev/null
+++ b/lmorpho/morphology.hpp
@@ -0,0 +1,31 @@
+#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/modcc/textbuffer.hpp b/modcc/textbuffer.hpp
index 55b6da678daa5ebebc5bc373bf91864984a4686e..5ab11396c813e943ee1b9b35d705bbf8ad3f55e0 100644
--- a/modcc/textbuffer.hpp
+++ b/modcc/textbuffer.hpp
@@ -9,6 +9,13 @@ public:
TextBuffer() {
text_.precision(std::numeric_limits<double>::max_digits10);
}
+ TextBuffer(const TextBuffer& other):
+ indent_(other.indent_),
+ indentation_width_(other.indentation_width_),
+ gutter_(other.gutter_),
+ text_(other.text_.str())
+ {}
+
TextBuffer& add_gutter();
void add_line(std::string const& line);
void add_line();
diff --git a/src/math.hpp b/src/math.hpp
index 50d760effa5626c328358b45644c4ca5b2fa6959..08d60ffdfe34fa420c4aecdb65c028acb41a7f59 100644
--- a/src/math.hpp
+++ b/src/math.hpp
@@ -81,6 +81,124 @@ int signum(T x) {
return (x<T(0)) - (x>T(0));
}
+// Quaternion implementation.
+// Represents w + x.i + y.j + z.k.
+
+struct quaternion {
+ double w, x, y, z;
+
+ constexpr quaternion(): w(0), x(0), y(0), z(0) {}
+
+ // scalar
+ constexpr quaternion(double w): w(w), x(0), y(0), z(0) {}
+
+ // vector (pure imaginary)
+ constexpr quaternion(double x, double y, double z): w(0), x(x), y(y), z(z) {}
+
+ // all 4-components
+ constexpr quaternion(double w, double x, double y, double z): w(w), x(x), y(y), z(z) {}
+
+ // equality testing
+ constexpr bool operator==(quaternion q) const {
+ return w==q.w && x==q.x && y==q.y && z==q.z;
+ }
+
+ constexpr bool operator!=(quaternion q) const {
+ return !(*this==q);
+ }
+
+ constexpr quaternion conj() const {
+ return {w, -x, -y, -z};
+ }
+
+ constexpr quaternion operator*(quaternion q) const {
+ return {w*q.w-x*q.x-y*q.y-z*q.z,
+ w*q.x+x*q.w+y*q.z-z*q.y,
+ w*q.y-x*q.z+y*q.w+z*q.x,
+ w*q.z+x*q.y-y*q.x+z*q.w};
+ }
+
+ quaternion& operator*=(quaternion q) {
+ return (*this=*this*q);
+ }
+
+ constexpr quaternion operator*(double d) const {
+ return {w*d, x*d, y*d, z*d};
+ }
+
+ quaternion& operator*=(double d) {
+ return (*this=*this*d);
+ }
+
+ friend constexpr quaternion operator*(double d, quaternion q) {
+ return q*d;
+ }
+
+ constexpr quaternion operator+(quaternion q) const {
+ return {w+q.w, x+q.x, y+q.y, z+q.z};
+ }
+
+ quaternion& operator+=(quaternion q) {
+ w += q.w;
+ x += q.x;
+ y += q.y;
+ z += q.z;
+ return *this;
+ }
+
+ constexpr quaternion operator-() const {
+ return {-w, -x, -y, -z};
+ }
+
+ constexpr quaternion operator-(quaternion q) const {
+ return {w-q.w, x-q.x, y-q.y, z-q.z};
+ }
+
+ quaternion& operator-=(quaternion q) {
+ w -= q.w;
+ x -= q.x;
+ y -= q.y;
+ z -= q.z;
+ return *this;
+ }
+
+ constexpr double sqnorm() const {
+ return w*w+x*x+y*y+z*z;
+ }
+
+ double norm() const {
+ return std::sqrt(sqnorm());
+ }
+
+ // Conjugation a ^ b = b a b*.
+ constexpr quaternion operator^(quaternion b) const {
+ return b*(*this)*b.conj();
+ }
+
+ quaternion& operator^=(quaternion b) {
+ return *this = b*(*this)*b.conj();
+ }
+
+ // add more as required...
+};
+
+// Quaternionic representations of axis rotations.
+// Given a vector v = (x, y, z), and r a quaternion representing
+// a rotation as below, then then (0, x, y, z) ^ r = (0, x', y', z')
+// represents the rotated vector.
+
+inline quaternion rotation_x(double phi) {
+ return {std::cos(phi/2), std::sin(phi/2), 0, 0};
+}
+
+inline quaternion rotation_y(double theta) {
+ return {std::cos(theta/2), 0, std::sin(theta/2), 0};
+}
+
+inline quaternion rotation_z(double psi) {
+ return {std::cos(psi/2), 0, 0, std::sin(psi/2)};
+}
+
} // namespace math
} // namespace mc
} // namespace nest
diff --git a/tests/validation/tinyopt.hpp b/src/tinyopt.hpp
similarity index 58%
rename from tests/validation/tinyopt.hpp
rename to src/tinyopt.hpp
index 90eb145f312511460786fa7f5a004604c5c54b8e..1e50a926d849d4041132ad4ee1b7aa17c8b4bd96 100644
--- a/tests/validation/tinyopt.hpp
+++ b/src/tinyopt.hpp
@@ -2,15 +2,15 @@
#include <cstring>
#include <iostream>
+#include <iterator>
#include <sstream>
#include <stdexcept>
#include <string>
+#include <utility>
+#include <vector>
-#include "../gtest.h"
-
-#include <communication/global_policy.hpp>
-
-#include "util/optional.hpp"
+#include <util/meta.hpp>
+#include <util/optional.hpp>
namespace nest {
namespace mc {
@@ -37,26 +37,58 @@ void usage(const char* argv0, const std::string& usage_str, const std::string& p
std::cerr << "Usage: " << basename << " " << usage_str << "\n";
}
-template <typename V = std::string>
-util::optional<V> parse_opt(char **& argp, char shortopt, const char* longopt=nullptr) {
+template <typename V>
+struct default_parser {
+ util::optional<V> operator()(const std::string& text) const {
+ V v;
+ std::istringstream stream(text);
+ stream >> v;
+ return stream? util::just(v): util::nothing;
+ }
+};
+
+template <typename V>
+class keyword_parser {
+ std::vector<std::pair<std::string, V>> map_;
+
+public:
+ template <typename KeywordPairs>
+ keyword_parser(const KeywordPairs& pairs): map_(std::begin(pairs), std::end(pairs)) {}
+
+ util::optional<V> operator()(const std::string& text) const {
+ for (const auto& p: map_) {
+ if (text==p.first) return p.second;
+ }
+ return util::nothing;
+ }
+};
+
+template <typename KeywordPairs>
+auto keywords(const KeywordPairs& pairs) -> keyword_parser<decltype(std::begin(pairs)->second)> {
+ return keyword_parser<decltype(std::begin(pairs)->second)>(pairs);
+}
+
+template <typename V = std::string, typename P = default_parser<V>, typename = util::enable_if_t<!std::is_same<V, void>::value>>
+util::optional<V> parse_opt(char **& argp, char shortopt, const char* longopt=nullptr, const P& parse = P{}) {
const char* arg = argp[0];
if (!arg || arg[0]!='-') {
return util::nothing;
}
- std::stringstream buf;
+ std::string text;
if (arg[1]=='-' && longopt) {
const char* rest = arg+2;
const char* eq = std::strchr(rest, '=');
if (!std::strcmp(rest, longopt)) {
- buf.str(argp[1]? argp[1]: throw parse_opt_error(arg, "missing argument"));
+ if (!argp[1]) throw parse_opt_error(arg, "missing argument");
+ text = argp[1];
argp += 2;
}
else if (eq && !std::strncmp(rest, longopt, eq-rest)) {
- buf.str(eq+1);
+ text = eq+1;
argp += 1;
}
else {
@@ -64,21 +96,20 @@ util::optional<V> parse_opt(char **& argp, char shortopt, const char* longopt=nu
}
}
else if (shortopt && arg[1]==shortopt && arg[2]==0) {
- buf.str(argp[1]? argp[1]: throw parse_opt_error(arg, "missing argument"));
+ if (!argp[1]) throw parse_opt_error(arg, "missing argument");
+ text = argp[1];
argp += 2;
}
else {
return util::nothing;
}
- V v;
- if (!(buf >> v)) {
- throw parse_opt_error(arg, "failed to parse option argument");
- }
+ auto v = parse(text);
+
+ if (!v) throw parse_opt_error(arg, "failed to parse option argument");
return v;
}
-template <>
util::optional<void> parse_opt(char **& argp, char shortopt, const char* longopt) {
if (!*argp || *argp[0]!='-') {
return util::nothing;
diff --git a/tests/unit/test_math.cpp b/tests/unit/test_math.cpp
index 1a5f442f2bba5a1b04db8fd64d1c39897ac6f8f0..d281487fb97d3391d7bf2b94fc57e6aa3540bf6b 100644
--- a/tests/unit/test_math.cpp
+++ b/tests/unit/test_math.cpp
@@ -111,3 +111,148 @@ TEST(math, infinity) {
EXPECT_EQ(std::numeric_limits<double>::infinity(), check.d);
EXPECT_EQ(std::numeric_limits<long double>::infinity(), check.ld);
}
+
+TEST(quaternion, ctor) {
+ // scalar
+ quaternion q1(3.5);
+
+ EXPECT_EQ(3.5, q1.w);
+ EXPECT_EQ(0., q1.x);
+ EXPECT_EQ(0., q1.y);
+ EXPECT_EQ(0., q1.z);
+
+ // pure imaginery
+ quaternion q2(1.5, 2.5, 3.5);
+
+ EXPECT_EQ(0, q2.w);
+ EXPECT_EQ(1.5, q2.x);
+ EXPECT_EQ(2.5, q2.y);
+ EXPECT_EQ(3.5, q2.z);
+
+ // all components
+ quaternion q3(0.5, 1.5, 2.5, 3.5);
+
+ EXPECT_EQ(0.5, q3.w);
+ EXPECT_EQ(1.5, q3.x);
+ EXPECT_EQ(2.5, q3.y);
+ EXPECT_EQ(3.5, q3.z);
+
+ // copy ctor
+ quaternion q4(q3);
+
+ EXPECT_EQ(0.5, q4.w);
+ EXPECT_EQ(1.5, q4.x);
+ EXPECT_EQ(2.5, q4.y);
+ EXPECT_EQ(3.5, q4.z);
+}
+
+TEST(quaternion, assign) {
+ quaternion q1(0.5, 1.5, 2.5, 3.5);
+ quaternion q2(7.3, -2.4, 11.1, -9);
+
+ q2 = q1;
+
+ EXPECT_EQ(0.5, q2.w);
+ EXPECT_EQ(1.5, q2.x);
+ EXPECT_EQ(2.5, q2.y);
+ EXPECT_EQ(3.5, q2.z);
+
+ q2 = -11.5;
+
+ EXPECT_EQ(-11.5, q2.w);
+ EXPECT_EQ(0, q2.x);
+ EXPECT_EQ(0, q2.y);
+ EXPECT_EQ(0, q2.z);
+}
+
+TEST(quaternion, equality) {
+ quaternion q1(1., 2., 3., 5.5);
+
+ quaternion q2(q1);
+ EXPECT_EQ(q1, q2);
+
+ q2 = q1;
+ q2.w += 0.5;
+ EXPECT_NE(q1, q2);
+
+ q2 = q1;
+ q2.x += 0.5;
+ EXPECT_NE(q1, q2);
+
+ q2 = q1;
+ q2.y += 0.5;
+ EXPECT_NE(q1, q2);
+
+ q2 = q1;
+ q2.z += 0.5;
+ EXPECT_NE(q1, q2);
+}
+
+TEST(quaternion, unaryop) {
+ quaternion q(2, -3, 4.5, -5);
+
+ EXPECT_EQ(quaternion(-2, 3, -4.5, 5), -q);
+ EXPECT_EQ(quaternion(2, 3, -4.5, 5), q.conj());
+
+ quaternion r(10, 6, 4, 37);
+ EXPECT_EQ(1521., r.sqnorm());
+ EXPECT_DOUBLE_EQ(39., r.norm());
+}
+
+TEST(quaternion, binop) {
+ quaternion q1(2, -3, 4.5, -5);
+ quaternion q2(0.5, 1.5, -2.5, 3.5);
+
+ EXPECT_EQ(quaternion(2.5, -1.5, 2, -1.5), q1+q2);
+ EXPECT_EQ(quaternion(1.5, -4.5, 7, -8.5), q1-q2);
+ EXPECT_EQ(quaternion(34.25, 4.75, 0.25, 5.25), q1*q2);
+ EXPECT_EQ(quaternion(42, -41.5, 71, -131), q1^q2);
+}
+
+TEST(quaternion, assignop) {
+ quaternion q1(2, -3, 4.5, -5);
+ quaternion q2(0.5, 1.5, -2.5, 3.5);
+
+ quaternion q;
+ q = q1;
+ q += q2;
+ EXPECT_EQ(q1+q2, q);
+
+ q = q1;
+ q -= q2;
+ EXPECT_EQ(q1-q2, q);
+
+ q = q1;
+ q *= q2;
+ EXPECT_EQ(q1*q2, q);
+
+ q = q1;
+ q ^= q2;
+ EXPECT_EQ(q1^q2, q);
+}
+
+TEST(quaternion, rotate) {
+ double deg_to_rad = pi<double>()/180.;
+ double sqrt3o2 = std::sqrt(3.)/2.;
+ double eps = 1e-15;
+
+ quaternion q(0, 1, 2, 3);
+
+ auto r = q^rotation_x(deg_to_rad*30);
+ EXPECT_NEAR(0., r.w, eps);
+ EXPECT_NEAR(q.x, r.x, eps);
+ EXPECT_NEAR(q.y*sqrt3o2-q.z/2., r.y, eps);
+ EXPECT_NEAR(q.y/2.+q.z*sqrt3o2, r.z, eps);
+
+ r = q^rotation_y(deg_to_rad*30);
+ EXPECT_NEAR(0., r.w, eps);
+ EXPECT_NEAR(q.x*sqrt3o2+q.z/2., r.x, eps);
+ EXPECT_NEAR(q.y, r.y, eps);
+ EXPECT_NEAR(-q.x/2.+q.z*sqrt3o2, r.z, eps);
+
+ r = q^rotation_z(deg_to_rad*30);
+ EXPECT_NEAR(0., r.w, eps);
+ EXPECT_NEAR(q.x*sqrt3o2-q.y/2., r.x, eps);
+ EXPECT_NEAR(q.x/2.+q.y*sqrt3o2, r.y, eps);
+ EXPECT_NEAR(q.z, r.z, eps);
+}
diff --git a/tests/validation/validate.cpp b/tests/validation/validate.cpp
index a2f2f8322cc9c63e03eddf41957521e5bbdffb99..e120ecc60ffa3825d0c5e046d2eca12a2b8e62f0 100644
--- a/tests/validation/validate.cpp
+++ b/tests/validation/validate.cpp
@@ -5,10 +5,10 @@
#include <exception>
#include <communication/global_policy.hpp>
+#include <tinyopt.hpp>
#include "../gtest.h"
-#include "tinyopt.hpp"
#include "validation_data.hpp"
using namespace nest::mc;
@@ -46,10 +46,10 @@ int main(int argc, char **argv) {
else if (auto o = parse_opt<float>(arg, 'd', "min-dt")) {
g_trace_io.set_min_dt(*o);
}
- else if (auto o = parse_opt<void>(arg, 'v', "verbose")) {
+ else if (auto o = parse_opt(arg, 'v', "verbose")) {
g_trace_io.set_verbose(true);
}
- else if (auto o = parse_opt<void>(arg, 'h', "help")) {
+ else if (auto o = parse_opt(arg, 'h', "help")) {
to::usage(argv[0], usage_str);
return 0;
}