From 7f9288fbe2cf223d2f1426c3c59569cfa05ee4e0 Mon Sep 17 00:00:00 2001
From: Sam Yates <yates@cscs.ch>
Date: Mon, 6 Mar 2017 16:33:26 +0100
Subject: [PATCH] Morphology generation with L-systems (#162)
Adds a stand-alone program for the generation of random morphologies form a L-system description. The algorithm is that of Burke (1992), with some of the extensions provided by Ascoli et al. (2001).
Two sets of L-system parameters have been included, corresponding to alpha motoneurons and Purkinje cells, but there is certainly something wrong with the data for the latter, and more correct numbers will probably need to be synthesized from existing Purkinje cell morphological information.
Documentation for `lmorpho` is incomplete, but the command line help (`--help`) goes some way to explain the usage. In order to get output, one must specify `--swc` or `--pvec` (or both) to emit SWC files or the structural parent vectors. Coarser discretization can be obtained with the `--segment` option.
Some minor modifications have been included in other parts of the source repo:
* Added copy constructor for `TextBuffer` in `modcc/textbuffer.hpp`, required to keep clang++ 3.5 happy.
* Disabled 'maybe-uninitialized' warnings from gcc, as these get raised improperly with some uses of `util::optional`. (Seems that this is also an issue for `boost::optional`.)
* Moved `tinyopt.hpp` into the `src/` directory, and extended the functionality a bit to support keyword arguments.
* `validate.cpp` modified slightly to accommodate new `tinyopt`.
---
CMakeLists.txt | 1 +
cmake/CompilerOptions.cmake | 17 +-
lmorpho/CMakeLists.txt | 5 +
lmorpho/README.md | 115 ++++++++++
lmorpho/lmorpho.cpp | 115 ++++++++++
lmorpho/lsys_models.cpp | 131 +++++++++++
lmorpho/lsys_models.hpp | 9 +
lmorpho/lsystem.cpp | 309 ++++++++++++++++++++++++++
lmorpho/lsystem.hpp | 131 +++++++++++
lmorpho/morphio.cpp | 179 +++++++++++++++
lmorpho/morphio.hpp | 93 ++++++++
lmorpho/morphology.cpp | 62 ++++++
lmorpho/morphology.hpp | 31 +++
modcc/textbuffer.hpp | 7 +
src/math.hpp | 118 ++++++++++
{tests/validation => src}/tinyopt.hpp | 63 ++++--
tests/unit/test_math.cpp | 145 ++++++++++++
tests/validation/validate.cpp | 6 +-
18 files changed, 1510 insertions(+), 27 deletions(-)
create mode 100644 lmorpho/CMakeLists.txt
create mode 100644 lmorpho/README.md
create mode 100644 lmorpho/lmorpho.cpp
create mode 100644 lmorpho/lsys_models.cpp
create mode 100644 lmorpho/lsys_models.hpp
create mode 100644 lmorpho/lsystem.cpp
create mode 100644 lmorpho/lsystem.hpp
create mode 100644 lmorpho/morphio.cpp
create mode 100644 lmorpho/morphio.hpp
create mode 100644 lmorpho/morphology.cpp
create mode 100644 lmorpho/morphology.hpp
rename {tests/validation => src}/tinyopt.hpp (58%)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index fd594311..d181f424 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 3788419a..c3982b94 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 00000000..1bb074d8
--- /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 00000000..561bbebd
--- /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 00000000..65f32009
--- /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 00000000..c26437d4
--- /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 00000000..96638e5a
--- /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 00000000..7e46c913
--- /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 00000000..f3ce4deb
--- /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 00000000..e6259be5
--- /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 00000000..37ec56ae
--- /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 00000000..7e9c60ea
--- /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 00000000..d1adaa6f
--- /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 55b6da67..5ab11396 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 50d760ef..08d60ffd 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 90eb145f..1e50a926 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 1a5f442f..d281487f 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 a2f2f832..e120ecc6 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;
}
--
GitLab