Merge pull request #11 from eth-cscs/master

merge master into my fork

.gitignore

@@ -18,6 +18,8 @@
 *.swo
 *.swn
 *.swq
+*.swm
+*.swl
 
 # tar files
 *.tar
@@ -65,6 +67,7 @@ external/tmp
 mechanisms/*.hpp
 
 # build path
-build
+build*
 
 commit.msg
+

CMakeLists.txt

@@ -51,7 +51,7 @@ if(WITH_MPI)
 endif()
 
 # Profiler support
-set(WITH_PROFILING OFF CACHE BOOL "use built in profiling of miniapp" )
+set(WITH_PROFILING OFF CACHE BOOL "use built-in profiling of miniapp" )
 if(WITH_PROFILING)
     add_definitions(-DWITH_PROFILING)
 endif()
@@ -62,18 +62,45 @@ if(SYSTEM_CRAY)
     set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -dynamic")
 endif()
 
-# targets for extermal dependencies
-include(ExternalProject)
-externalproject_add(modparser
-    PREFIX ${CMAKE_BINARY_DIR}/external
-    CMAKE_ARGS "-DCMAKE_INSTALL_PREFIX=${CMAKE_BINARY_DIR}/external"
-               "-DCMAKE_CXX_FLAGS=${SAVED_CXX_FLAGS}"
-               "-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}"
-    BINARY_DIR "${CMAKE_BINARY_DIR}/external/modparser"
-    STAMP_DIR  "${CMAKE_BINARY_DIR}/external/"
-    TMP_DIR    "${CMAKE_BINARY_DIR}/external/tmp"
-    SOURCE_DIR "${CMAKE_SOURCE_DIR}/external/modparser"
+# vectorization target
+set(VECTORIZE_TARGET "none" CACHE STRING "CPU target for vectorization {KNL,AVX,AVX2}")
+
+if(VECTORIZE_TARGET STREQUAL "KNL")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CXXOPT_KNL}")
+elseif(VECTORIZE_TARGET STREQUAL "AVX")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CXXOPT_AVX}")
+elseif(VECTORIZE_TARGET STREQUAL "AVX2")
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CXXOPT_AVX2}")
+endif()
+
+# whether to generate optimized kernels from NMODL
+set(USE_OPTIMIZED_KERNELS OFF CACHE BOOL "generate optimized code that vectorizes with the Intel compiler")
+
+# Only build modcc if it has not already been installed.
+# This is useful if cross compiling for KNL, when it is not desirable to compile
+# modcc with the same flags that are used for the KNL target.
+find_program(MODCC_BIN modcc)
+set(modcc "${MODCC_BIN}")
+set(use_external_modcc OFF BOOL)
+
+# the modcc executable was not found, so build our own copy
+if(MODCC_BIN STREQUAL "MODCC_BIN-NOTFOUND")
+    include(ExternalProject)
+    externalproject_add(modparser
+        PREFIX ${CMAKE_BINARY_DIR}/external
+        CMAKE_ARGS "-DCMAKE_INSTALL_PREFIX=${CMAKE_BINARY_DIR}/external"
+                   "-DCMAKE_CXX_FLAGS=${SAVED_CXX_FLAGS}"
+                   "-DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}"
+        BINARY_DIR "${CMAKE_BINARY_DIR}/external/modparser"
+        STAMP_DIR  "${CMAKE_BINARY_DIR}/external/"
+        TMP_DIR    "${CMAKE_BINARY_DIR}/external/tmp"
+        SOURCE_DIR "${CMAKE_SOURCE_DIR}/external/modparser"
     )
+    # Set up environment to use the version of modcc that is compiled
+    # as the ExternalProject above.
+    set(use_external_modcc ON)
+    set(modcc "${CMAKE_BINARY_DIR}/external/bin/modcc")
+endif()
 
 
 include_directories(${CMAKE_SOURCE_DIR}/external)
@@ -82,7 +109,6 @@ include_directories(${CMAKE_SOURCE_DIR}/include)
 include_directories(${CMAKE_SOURCE_DIR}/src)
 include_directories(${CMAKE_SOURCE_DIR}/miniapp)
 include_directories(${CMAKE_SOURCE_DIR})
-
 if( "${WITH_TBB}" STREQUAL "ON" )
     include_directories(${TBB_INCLUDE_DIRS})
 endif()

README.md

@@ -77,3 +77,147 @@ cmake <path to CMakeLists.txt> -DWITH_TBB=ON -DSYSTEM_CRAY=ON
 cmake <path to CMakeLists.txt> -DWITH_TBB=ON -DWITH_MPI=ON -DSYSTEM_CRAY=ON
 
 ```
+
+# targetting KNL
+
+## build modparser                                                                                       
+
+The source to source compiler "modparser" that generates the C++/CUDA kernels for the ion channels and synapses is in a separate repository.
+It is included in our project as a git submodule, and by default it will be built with the same compiler and flags that are used to build the miniapp and tests.
+
+This can cause problems if we are cross compiling, e.g. for KNL, because the modparser compiler might not be runnable on the compilation node.
+CMake will look for the source to source compiler executable, `modcc`, in the `PATH` environment variable, and will use the version if finds instead of building its own.
+
+Modparser requires a C++11 compiler, and has been tested on GCC, Intel, and Clang compilers
+  - if the default compiler on your is some ancient version of gcc you might need to load a module/set the CC and CXX environment variables.
+
+```bash
+git clone git@github.com:eth-cscs/modparser.git
+cd modparser
+
+# example of setting a C++11 compiler
+export CXX=`which gcc-4.8`
+
+cmake .
+make -j
+
+# set path and test that you can see modcc
+export PATH=`pwd`/bin:$PATH
+which modcc
+```
+
+## set up environment
+
+- source the intel compilers
+- source the TBB vars
+- I have only tested with the latest stable version from online, not the version that comes installed sometimes with the Intel compilers.
+
+## build miniapp
+
+```bash
+# clone the repo and set up the submodules
+git clone TODO
+cd cell_algorithms
+git submodule init
+git submodule update
+
+# make a path for out of source build
+mkdir build_knl
+cd build_knl
+
+## build miniapp
+
+# clone the repo and set up the submodules
+git clone https://github.com/bcumming/cell_algorithms.git
+cd cell_algorithms
+
+# checkout the knl branch
+git checkout knl
+
+# setup submodules
+git submodule init
+git submodule update
+
+# make a path for out of source build
+mkdir build_knl
+cd build_knl
+
+# run cmake with all the magic flags
+export CC=`which icc`
+export CXX=`which icpc`
+cmake .. -DCMAKE_BUILD_TYPE=release -DWITH_TBB=ON -DWITH_PROFILING=ON -DVECTORIZE_TARGET=KNL -DUSE_OPTIMIZED_KERNELS=ON
+make -j
+```
+
+The flags passed into cmake are described:
+  - `-DCMAKE_BUILD_TYPE=release` : build in release mode with `-O3`.
+  - `-WITH_TBB=ON` : use TBB for threading on multicore
+  - `-DWITH_PROFILING=ON` : use internal profilers that print profiling report at end
+  - `-DVECTORIZE_TARGET=KNL` : generate AVX512 instructions, alternatively you can use:
+    - `AVX2` for Haswell & Broadwell
+    - `AVX` for Sandy Bridge and Ivy Bridge
+  - `-DUSE_OPTIMIZED_KERNELS=ON` : tell the source to source compiler to generate optimized kernels that use Intel extensions
+    - without these vectorized code will not be generated.
+
+## run tests
+
+Run some unit tests
+```bash
+cd tests
+./test.exe
+cd ..
+```
+
+## run miniapp
+
+The miniapp is the target for benchmarking.
+First, we can run a small problem to check the build.
+For the small test run, the parameters have the following meaning
+  - `-n 1000` : 1000 cells
+  - `-s 200` : 200 synapses per cell
+  - `-t 20`  : simulated for 20ms
+  - `-p 0`   : no file output of voltage traces
+
+The number of cells is the number of discrete tasks that are distributed to the threads in each large time integration period.
+The number of synapses per cell is the amount of computational work per cell/task.
+Realistic cells have anywhere in the range of 1,000-10,000 synapses per cell.
+
+```bash
+cd miniapp
+
+# a small run to check that everything works
+./miniapp.exe -n 1000 -s 200 -t 20 -p 0
+
+# a larger run for generating meaninful benchmarks
+./miniapp.exe -n 2000 -s 2000 -t 100 -p 0
+```
+
+This generates the following profiler output (some reformatting to make the table work):
+
+```
+              ---------------------------------------
+             |       small       |       large       |
+-------------|-------------------|-------------------|
+total        |  0.791     100.0  | 38.593     100.0  |
+  stepping   |  0.738      93.3  | 36.978      95.8  |
+    matrix   |  0.406      51.3  |  6.034      15.6  |
+      solve  |  0.308      38.9  |  4.534      11.7  |
+      setup  |  0.082      10.4  |  1.260       3.3  |
+      other  |  0.016       2.0  |  0.240       0.6  |
+    state    |  0.194      24.5  | 23.235      60.2  |
+      expsyn |  0.158      20.0  | 22.679      58.8  |
+      hh     |  0.014       1.7  |  0.215       0.6  |
+      pas    |  0.003       0.4  |  0.053       0.1  |
+      other  |  0.019       2.4  |  0.287       0.7  |
+    current  |  0.107      13.5  |  7.106      18.4  |
+      expsyn |  0.047       5.9  |  6.118      15.9  |
+      pas    |  0.028       3.5  |  0.476       1.2  |
+      hh     |  0.006       0.7  |  0.096       0.2  |
+      other  |  0.026       3.3  |  0.415       1.1  |
+    events   |  0.005       0.6  |  0.125       0.3  |
+    sampling |  0.003       0.4  |  0.051       0.1  |
+    other    |  0.024       3.0  |  0.428       1.1  |
+  other      |  0.053       6.7  |  1.614       4.2  |
+-----------------------------------------------------
+```
+

cmake/CompilerOptions.cmake

@@ -13,10 +13,24 @@ if(${CMAKE_CXX_COMPILER_ID} MATCHES "Clang")
     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
+    set(CXXOPT_KNL "-march=knl")
+    set(CXXOPT_AVX "-mavx")
+    set(CXXOPT_AVX2 "-march=core-avx2")
+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
 
     set(CXXOPT_WALL "${CXXOPT_WALL} -wd488")
+
+    # compiler flags for generating KNL-specific AVX512 instructions
+    set(CXXOPT_KNL "-xMIC-AVX512")
+    set(CXXOPT_AVX "-mavx")
+    set(CXXOPT_AVX2 "-march=core-avx2")
 endif()
 

modparser @ b200bf63

-Subproject commit 588ca1a5ea28ef04d17b318e754d479e5489eb9a
+Subproject commit b200bf6376a2dc30edea98fcc2375fc9be095135

vector @ c8678e80

-Subproject commit 8284611f05b0fbe21a1f84630e2726015cb1d96d
+Subproject commit c8678e80660cd63b293ade80ece8eed2249c1b06

mechanisms/CMakeLists.txt

+# the list of built-in mechanisms to be provided by default
 set(mechanisms pas hh expsyn exp2syn)
 
-set(modcc "${CMAKE_BINARY_DIR}/external/bin/modcc")
+# set the flags for the modcc compiler that converts NMODL
+# files to C++/CUDA source.
+set(modcc_flags "-t cpu")
+if(USE_OPTIMIZED_KERNELS) # generate optimized kernels
+    set(modcc_flags ${modcc_flags} -O)
+endif()
 
+# generate source for each mechanism
 foreach(mech ${mechanisms})
     set(mod "${CMAKE_CURRENT_SOURCE_DIR}/mod/${mech}.mod")
     set(hpp "${CMAKE_CURRENT_SOURCE_DIR}/${mech}.hpp")
-    add_custom_command(OUTPUT "${hpp}"
-                       DEPENDS modparser "${mod}"
-                       WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}"
-                       COMMAND "${modcc}" -t cpu -o "${hpp}" "${mod}")
+    if(use_external_modcc)
+        add_custom_command(
+           OUTPUT "${hpp}"
+           WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}"
+           COMMAND ${modcc} ${modcc_flags} ${mod} -o ${hpp}
+       )
+    else()
+        add_custom_command(
+            OUTPUT "${hpp}"
+            DEPENDS modparser "${mod}"
+            WORKING_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}"
+            COMMAND ${modcc} ${modcc_flags} ${mod} -o ${hpp}
+        )
+    endif()
     set_source_files_properties("${hpp}" PROPERTIES GENERATED TRUE)
     list(APPEND all_mod_hpps "${hpp}")
 endforeach()

mechanisms/generate.sh

+#!/usr/bin/env bash
+
 #flags="-t cpu -O"
 flags="-t cpu"
 

miniapp/CMakeLists.txt

 set(HEADERS
 )
 set(MINIAPP_SOURCES
-    #    mpi.cpp
     io.cpp
     miniapp.cpp
+    miniapp_recipes.cpp
 )
 
 add_executable(miniapp.exe ${MINIAPP_SOURCES} ${HEADERS})

miniapp/io.cpp

-#include <fstream>
+#include <algorithm>
 #include <exception>
+#include <fstream>
+#include <istream>
+#include <type_traits>
 
 #include <tclap/CmdLine.h>
+#include <json/src/json.hpp>
+
+#include <util/meta.hpp>
+#include <util/optional.hpp>
 
 #include "io.hpp"
 
+// Let TCLAP understand value arguments that are of an optional type.
+
+namespace TCLAP {
+    template <typename V>
+    struct ArgTraits<nest::mc::util::optional<V>> {
+        using ValueCategory = ValueLike;
+    };
+} // namespace TCLAP
+
 namespace nest {
 namespace mc {
+
+namespace util {
+    // Using static here because we do not want external linkage for this operator.
+    template <typename V>
+    static std::istream& operator>>(std::istream& I, optional<V>& v) {
+        V u;
+        if (I >> u) {
+            v = u;
+        }
+        return I;
+    }
+}
+
 namespace io {
 
-/// read simulation options from json file with name fname
-/// if file name is empty or if file is not a valid json file a default
-/// set of parameters is returned :
-///      1000 cells, 500 synapses per cell, 100 compartments per segment
+// Override annoying parameters listed back-to-front behaviour.
+//
+// TCLAP argument creation _prepends_ its arguments to the internal
+// list (_argList), where standard options --help etc. are already
+// pre-inserted.
+//
+// reorder_arguments() reverses the arguments to restore ordering,
+// and moves the standard options to the end.
+class CustomCmdLine: public TCLAP::CmdLine {
+public:
+    CustomCmdLine(const std::string &message, const std::string &version = "none"):
+        TCLAP::CmdLine(message, ' ', version, true)
+    {}
+
+    void reorder_arguments() {
+        _argList.reverse();
+        for (auto opt: {"help", "version", "ignore_rest"}) {
+            auto i = std::find_if(
+                _argList.begin(), _argList.end(),
+                [&opt](TCLAP::Arg* a) { return a->getName()==opt; });
+
+            if (i!=_argList.end()) {
+                auto a = *i;
+                _argList.erase(i);
+                _argList.push_back(a);
+            }
+        }
+    }
+};
+
+// Update an option value from command line argument if set.
+template <
+    typename T,
+    typename Arg,
+    typename = util::enable_if_t<std::is_base_of<TCLAP::Arg, Arg>::value>
+>
+static void update_option(T& opt, Arg& arg) {
+    if (arg.isSet()) {
+        opt = arg.getValue();
+    }
+}
+
+// Update an option value from json object if key present.
+template <typename T>
+static void update_option(T& opt, const nlohmann::json& j, const std::string& key) {
+    if (j.count(key)) {
+        opt = j[key];
+    }
+}
+
+// --- special case for string due to ambiguous overloading in json library.
+static void update_option(std::string& opt, const nlohmann::json& j, const std::string& key) {
+    if (j.count(key)) {
+        opt = j[key].get<std::string>();
+    }
+}
+
+// --- special case for optional values.
+template <typename T>
+static void update_option(util::optional<T>& opt, const nlohmann::json& j, const std::string& key) {
+    if (j.count(key)) {
+        auto value = j[key];
+        if (value.is_null()) {
+            opt = util::nothing;
+        }
+        else {
+            opt = value;
+        }
+    }
+}
+
+// Read options from (optional) json file and command line arguments.
 cl_options read_options(int argc, char** argv) {
 
-    // set default options
-    const cl_options default_options{"", 1000, 500, "expsyn", 100, 100., 0.025, false};
+    // Default options:
+    const cl_options defopts{
+        1000,       // number of cells
+        500,        // synapses_per_cell
+        "expsyn",   // synapse type
+        100,        // compartments_per_segment
+        100.,       // tfinal
+        0.025,      // dt
+        false,      // all_to_all
+        false,      // probe_soma_only
+        0.0,        // probe_ratio
+        "trace_",   // trace_prefix
+        util::nothing,  // trace_max_gid
+
+        // spike_output_parameters:
+        false,      // spike output
+        false,      // single_file_per_simulation
+        true,       // Overwrite outputfile if exists
+        "./",       // output path
+        "spikes",   // file name
+        "gdf"       // file extension
+    };
 
     cl_options options;
-    // parse command line arguments
+    std::string save_file = "";
+
+    // Parse command line arguments.
     try {
-        TCLAP::CmdLine cmd("mod2c performance benchmark harness", ' ', "0.1");
+        CustomCmdLine cmd("nest mc miniapp harness", "0.1");
 
+        TCLAP::ValueArg<std::string> ifile_arg(
+            "i", "ifile",
+            "read parameters from json-formatted file <file name>",
+            false, "","file name", cmd);
+        TCLAP::ValueArg<std::string> ofile_arg(
+            "o", "ofile",
+            "save parameters to json-formatted file <file name>",
+            false, "","file name", cmd);
         TCLAP::ValueArg<uint32_t> ncells_arg(
             "n", "ncells", "total number of cells in the model",
-            false, 1000, "non negative integer");
+            false, defopts.cells, "integer", cmd);
         TCLAP::ValueArg<uint32_t> nsynapses_arg(
             "s", "nsynapses", "number of synapses per cell",
-            false, 500, "non negative integer");
+            false, defopts.synapses_per_cell, "integer", cmd);
         TCLAP::ValueArg<std::string> syntype_arg(
-            "S", "syntype", "type of synapse (expsyn or exp2syn)",
-            false, "expsyn", "synapse type");
+            "S", "syntype", "specify synapse type: expsyn or exp2syn",
+            false, defopts.syn_type, "string", cmd);
         TCLAP::ValueArg<uint32_t> ncompartments_arg(
             "c", "ncompartments", "number of compartments per segment",
-            false, 100, "non negative integer");
-        TCLAP::ValueArg<std::string> ifile_arg(
-            "i", "ifile", "json file with model parameters",
-            false, "","file name string");
+            false, defopts.compartments_per_segment, "integer", cmd);
         TCLAP::ValueArg<double> tfinal_arg(
-            "t", "tfinal", "time to simulate in ms",
-            false, 100., "positive real number");
+            "t", "tfinal", "run simulation to <time> ms",
+            false, defopts.tfinal, "time", cmd);
         TCLAP::ValueArg<double> dt_arg(
-            "d", "dt", "time step size in ms",
-            false, 0.025, "positive real number");
+            "d", "dt", "set simulation time step to <time> ms",
+            false, defopts.dt, "time", cmd);
         TCLAP::SwitchArg all_to_all_arg(
-            "a","alltoall","all to all network", cmd, false);
-
-        cmd.add(ncells_arg);
-        cmd.add(nsynapses_arg);
-        cmd.add(syntype_arg);
-        cmd.add(ncompartments_arg);
-        cmd.add(ifile_arg);
-        cmd.add(dt_arg);
-        cmd.add(tfinal_arg);
+            "m","alltoall","all to all network", cmd, false);
+        TCLAP::ValueArg<double> probe_ratio_arg(
+            "p", "probe-ratio", "proportion between 0 and 1 of cells to probe",
+            false, defopts.probe_ratio, "proportion", cmd);
+        TCLAP::SwitchArg probe_soma_only_arg(
+            "X", "probe-soma-only", "only probe cell somas, not dendrites", cmd, false);
+        TCLAP::ValueArg<std::string> trace_prefix_arg(
+            "P", "prefix", "write traces to files with prefix <prefix>",
+            false, defopts.trace_prefix, "stringr", cmd);
+        TCLAP::ValueArg<util::optional<unsigned>> trace_max_gid_arg(
+            "T", "trace-max-gid", "only trace probes on cells up to and including <gid>",
+            false, defopts.trace_max_gid, "gid", cmd);
+        TCLAP::SwitchArg spike_output_arg(
+            "f","spike_file_output","save spikes to file", cmd, false);
 
+        cmd.reorder_arguments();
         cmd.parse(argc, argv);
 
-        options.cells = ncells_arg.getValue();
-        options.synapses_per_cell = nsynapses_arg.getValue();
-        options.syn_type = syntype_arg.getValue();
-        options.compartments_per_segment = ncompartments_arg.getValue();
-        options.ifname = ifile_arg.getValue();
-        options.tfinal = tfinal_arg.getValue();
-        options.dt = dt_arg.getValue();
-        options.all_to_all = all_to_all_arg.getValue();
+        options = defopts;
+
+        std::string ifile_name = ifile_arg.getValue();
+        if (ifile_name != "") {
+            // Read parameters from specified JSON file first, to allow
+            // overriding arguments on the command line.
+            std::ifstream fid(ifile_name);
+            if (fid) {
+                try {
+                    nlohmann::json fopts;
+                    fid >> fopts;
+
+                    update_option(options.cells, fopts, "cells");
+                    update_option(options.synapses_per_cell, fopts, "synapses");
+                    update_option(options.syn_type, fopts, "syn_type");
+                    update_option(options.compartments_per_segment, fopts, "compartments");
+                    update_option(options.dt, fopts, "dt");
+                    update_option(options.tfinal, fopts, "tfinal");
+                    update_option(options.all_to_all, fopts, "all_to_all");
+                    update_option(options.probe_ratio, fopts, "probe_ratio");
+                    update_option(options.probe_soma_only, fopts, "probe_soma_only");
+                    update_option(options.trace_prefix, fopts, "trace_prefix");
+                    update_option(options.trace_max_gid, fopts, "trace_max_gid");
+
+                    // Parameters for spike output
+                    update_option(options.spike_file_output, fopts, "spike_file_output");
+                    if (options.spike_file_output) {
+                        update_option(options.single_file_per_rank, fopts, "single_file_per_rank");
+                        update_option(options.over_write, fopts, "over_write");
+                        update_option(options.output_path, fopts, "output_path");
+                        update_option(options.file_name, fopts, "file_name");
+                        update_option(options.file_extension, fopts, "file_extension");
+                    }
+
+                }
+                catch (std::exception& e) {
+                    throw model_description_error(
+                        "unable to parse parameters in "+ifile_name+": "+e.what());
+                }
+            }
+            else {
+                throw usage_error("unable to open model parameter file "+ifile_name);
+            }
+        }
+
+        update_option(options.cells, ncells_arg);
+        update_option(options.synapses_per_cell, nsynapses_arg);
+        update_option(options.syn_type, syntype_arg);
+        update_option(options.compartments_per_segment, ncompartments_arg);
+        update_option(options.tfinal, tfinal_arg);
+        update_option(options.dt, dt_arg);
+        update_option(options.all_to_all, all_to_all_arg);
+        update_option(options.probe_ratio, probe_ratio_arg);
+        update_option(options.probe_soma_only, probe_soma_only_arg);
+        update_option(options.trace_prefix, trace_prefix_arg);
+        update_option(options.trace_max_gid, trace_max_gid_arg);
+        update_option(options.spike_file_output, spike_output_arg);
+
+        save_file = ofile_arg.getValue();
     }
-    // catch any exceptions in command line handling
-    catch(TCLAP::ArgException &e) {
-        std::cerr << "error: parsing command line arguments:\n  "
-                  << e.error() << " for arg " << e.argId() << "\n";
-        exit(1);
+    catch (TCLAP::ArgException& e) {
+        throw usage_error("error parsing command line argument "+e.argId()+": "+e.error());
     }
 
-    if(options.ifname == "") {
-        options.check_and_normalize();
-        return options;
-    }
-    else {
-        std::ifstream fid(options.ifname, std::ifstream::in);
+    // Save option values if requested.
+    if (save_file != "") {
+        std::ofstream fid(save_file);
+        if (fid) {
+            try {
+                nlohmann::json fopts;
 
-        if(fid.is_open()) {
-            // read json data in input file
-            nlohmann::json fopts;
-            fid >> fopts;
+                fopts["cells"] = options.cells;
+                fopts["synapses"] = options.synapses_per_cell;
+                fopts["syn_type"] = options.syn_type;
+                fopts["compartments"] = options.compartments_per_segment;
+                fopts["dt"] = options.dt;
+                fopts["tfinal"] = options.tfinal;
+                fopts["all_to_all"] = options.all_to_all;
+                fopts["probe_ratio"] = options.probe_ratio;
+                fopts["probe_soma_only"] = options.probe_soma_only;
+                fopts["trace_prefix"] = options.trace_prefix;
+                if (options.trace_max_gid) {
+                    fopts["trace_max_gid"] = options.trace_max_gid.get();
+                }
+                else {
+                    fopts["trace_max_gid"] = nullptr;
+                }
+                fid << std::setw(3) << fopts << "\n";
 
-            try {
-                options.cells = fopts["cells"];
-                options.synapses_per_cell = fopts["synapses"];
-                options.compartments_per_segment = fopts["compartments"];
-                options.dt = fopts["dt"];
-                options.tfinal = fopts["tfinal"];
-                options.all_to_all = fopts["all_to_all"];
-            }
-            catch(std::domain_error e) {
-                std::cerr << "error: unable to open parameters in "
-                          << options.ifname << " : " << e.what() << "\n";
-                exit(1);
             }
-            catch(std::exception e) {
-                std::cerr << "error: unable to open parameters in "
-                          << options.ifname << "\n";
-                exit(1);
+            catch (std::exception& e) {
+                throw model_description_error(
+                    "unable to save parameters in "+save_file+": "+e.what());
             }
-            options.check_and_normalize();
-            return options;
+        }
+        else {
+            throw usage_error("unable to write to model parameter file "+save_file);
         }
     }
-
-    return default_options;
+    return options;
 }
 
 std::ostream& operator<<(std::ostream& o, const cl_options& options) {
-    o << "simultion options:\n";
+    o << "simulation options:\n";
     o << "  cells                : " << options.cells << "\n";
     o << "  compartments/segment : " << options.compartments_per_segment << "\n";
     o << "  synapses/cell        : " << options.synapses_per_cell << "\n";
     o << "  simulation time      : " << options.tfinal << "\n";
     o << "  dt                   : " << options.dt << "\n";
     o << "  all to all network   : " << (options.all_to_all ? "yes" : "no") << "\n";
-    o << "  input file name      : " << options.ifname << "\n";
+    o << "  probe ratio          : " << options.probe_ratio << "\n";
+    o << "  probe soma only      : " << (options.probe_soma_only ? "yes" : "no") << "\n";
+    o << "  trace prefix         : " << options.trace_prefix << "\n";
+    o << "  trace max gid        : ";
+    if (options.trace_max_gid) {
+       o << *options.trace_max_gid;
+    }
+    o << "\n";
 
     return o;
 }

miniapp/io.hpp

 #pragma once
 
-#include <json/src/json.hpp>
+#include <string>
+#include <cstdint>
+#include <iosfwd>
+#include <stdexcept>
+#include <utility>
+
+#include <util/optional.hpp>
 
 namespace nest {
 namespace mc {
@@ -8,7 +14,6 @@ namespace io {
 
 // holds the options for a simulation run
 struct cl_options {
-    std::string ifname;
     uint32_t cells;
     uint32_t synapses_per_cell;
     std::string syn_type;
@@ -16,20 +21,37 @@ struct cl_options {
     double tfinal;
     double dt;
     bool all_to_all;
+    bool probe_soma_only;
+    double probe_ratio;
+    std::string trace_prefix;
+    util::optional<unsigned> trace_max_gid;
+
+    // Parameters for spike output
+    bool spike_file_output;
+    bool single_file_per_rank;
+    bool over_write;
+    std::string output_path;
+    std::string file_name;
+    std::string file_extension;
+};
 
-    // TODO the normalize bit should be moved to the model_parameters when
-    // we start having more models
-    void check_and_normalize() {
-        if(all_to_all) {
-            synapses_per_cell = cells - 1;
-        }
-    }
+class usage_error: public std::runtime_error {
+public:
+    template <typename S>
+    usage_error(S&& whatmsg): std::runtime_error(std::forward<S>(whatmsg)) {}
+};
+
+class model_description_error: public std::runtime_error {
+public:
+    template <typename S>
+    model_description_error(S&& whatmsg): std::runtime_error(std::forward<S>(whatmsg)) {}
 };
 
 std::ostream& operator<<(std::ostream& o, const cl_options& opt);
 
 cl_options read_options(int argc, char** argv);
 
+
 } // namespace io
 } // namespace mc
 } // namespace nest

miniapp/miniapp_recipes.cpp

+#include <cmath>
+#include <random>
+#include <vector>
+#include <utility>
+
+#include <cell.hpp>
+#include <util/debug.hpp>
+
+#include "miniapp_recipes.hpp"
+
+namespace nest {
+namespace mc {
+
+// TODO: split cell description into separate morphology, stimulus, probes, mechanisms etc.
+// description for greater data reuse.
+
+template <typename RNG>
+cell make_basic_cell(
+    unsigned compartments_per_segment,
+    unsigned num_synapses,
+    const std::string& syn_type,
+    RNG& rng)
+{
+    nest::mc::cell cell;
+
+    // Soma with diameter 12.6157 um and HH channel
+    auto soma = cell.add_soma(12.6157/2.0);
+    soma->add_mechanism(mc::hh_parameters());
+
+    // add dendrite of length 200 um and diameter 1 um with passive channel
+    std::vector<mc::cable_segment*> dendrites;
+    dendrites.push_back(cell.add_cable(0, mc::segmentKind::dendrite, 0.5, 0.5, 200));
+    dendrites.push_back(cell.add_cable(1, mc::segmentKind::dendrite, 0.5, 0.25,100));
+    dendrites.push_back(cell.add_cable(1, mc::segmentKind::dendrite, 0.5, 0.25,100));
+
+    for (auto d : dendrites) {
+        d->add_mechanism(mc::pas_parameters());
+        d->set_compartments(compartments_per_segment);
+        d->mechanism("membrane").set("r_L", 100);
+    }
+
+    cell.add_detector({0,0}, 20);
+
+    auto distribution = std::uniform_real_distribution<float>(0.f, 1.0f);
+    // distribute the synapses at random locations the terminal dendrites in a
+    // round robin manner; the terminal dendrites in this cell have indices 2 and 3.
+    nest::mc::parameter_list syn_default(syn_type);
+    for (unsigned i=0; i<num_synapses; ++i) {
+        cell.add_synapse({2+(i%2), distribution(rng)}, syn_default);
+    }
+
+    return cell;
+}
+
+class basic_cell_recipe: public recipe {
+public:
+    basic_cell_recipe(cell_gid_type ncell, basic_recipe_param param, probe_distribution pdist):
+        ncell_(ncell), param_(std::move(param)), pdist_(std::move(pdist))
+    {
+        delay_distribution_param = exp_param{param_.mean_connection_delay_ms
+                            - param_.min_connection_delay_ms};
+    }
+
+    cell_size_type num_cells() const override { return ncell_; }
+
+    cell get_cell(cell_gid_type i) const override {
+        auto gen = std::mt19937(i); // TODO: replace this with hashing generator...
+
+        auto cc = get_cell_count_info(i);
+        auto cell = make_basic_cell(param_.num_compartments, cc.num_targets,
+                        param_.synapse_type, gen);
+
+        EXPECTS(cell.num_segments()==basic_cell_segments);
+        EXPECTS(cell.probes().size()==0);
+        EXPECTS(cell.synapses().size()==cc.num_targets);
+        EXPECTS(cell.detectors().size()==cc.num_sources);
+
+        // add probes
+        if (cc.num_probes) {
+            unsigned n_probe_segs = pdist_.all_segments? basic_cell_segments: 1u;
+            for (unsigned i = 0; i<n_probe_segs; ++i) {
+                if (pdist_.membrane_voltage) {
+                    cell.add_probe({{i, i? 0.5: 0.0}, mc::probeKind::membrane_voltage});
+                }
+                if (pdist_.membrane_current) {
+                    cell.add_probe({{i, i? 0.5: 0.0}, mc::probeKind::membrane_current});
+                }
+            }
+        }
+        EXPECTS(cell.probes().size()==cc.num_probes);
+        return cell;
+    }
+
+    cell_count_info get_cell_count_info(cell_gid_type i) const override {
+        cell_count_info cc = {1, param_.num_synapses, 0 };
+
+        // probe this cell?
+        if (std::floor(i*pdist_.proportion)!=std::floor((i-1.0)*pdist_.proportion)) {
+            std::size_t np = pdist_.membrane_voltage + pdist_.membrane_current;
+            if (pdist_.all_segments) {
+                np *= basic_cell_segments;
+            }
+
+            cc.num_probes = np;
+        }
+
+        return cc;
+    }
+
+protected:
+    template <typename RNG>
+    cell_connection draw_connection_params(RNG& rng) const {
+        std::exponential_distribution<float> delay_dist(delay_distribution_param);
+        float delay = param_.min_connection_delay_ms + delay_dist(rng);
+        float weight = param_.syn_weight_per_cell/param_.num_synapses;
+        return cell_connection{{0, 0}, {0, 0}, weight, delay};
+    }
+
+    cell_gid_type ncell_;
+    basic_recipe_param param_;
+    probe_distribution pdist_;
+    static constexpr int basic_cell_segments = 4;
+
+    using exp_param = std::exponential_distribution<float>::param_type;
+    exp_param delay_distribution_param;
+};
+
+class basic_ring_recipe: public basic_cell_recipe {
+public:
+    basic_ring_recipe(cell_gid_type ncell,
+                      basic_recipe_param param,
+                      probe_distribution pdist = probe_distribution{}):
+        basic_cell_recipe(ncell, std::move(param), std::move(pdist)) {}
+
+    std::vector<cell_connection> connections_on(cell_gid_type i) const override {
+        std::vector<cell_connection> conns;
+        auto gen = std::mt19937(i); // TODO: replace this with hashing generator...
+
+        cell_gid_type prev = i==0? ncell_-1: i-1;
+        for (unsigned t=0; t<param_.num_synapses; ++t) {
+            cell_connection cc = draw_connection_params(gen);
+            cc.source = {prev, 0};
+            cc.dest = {i, t};
+            conns.push_back(cc);
+        }
+
+        return conns;
+    }
+};
+
+std::unique_ptr<recipe> make_basic_ring_recipe(
+        cell_gid_type ncell,
+        basic_recipe_param param,
+        probe_distribution pdist)
+{
+    return std::unique_ptr<recipe>(new basic_ring_recipe(ncell, param, pdist));
+}
+
+
+class basic_rgraph_recipe: public basic_cell_recipe {
+public:
+    basic_rgraph_recipe(cell_gid_type ncell,
+                      basic_recipe_param param,
+                      probe_distribution pdist = probe_distribution{}):
+        basic_cell_recipe(ncell, std::move(param), std::move(pdist)) {}
+
+    std::vector<cell_connection> connections_on(cell_gid_type i) const override {
+        std::vector<cell_connection> conns;
+        auto conn_param_gen = std::mt19937(i); // TODO: replace this with hashing generator...
+        auto source_gen = std::mt19937(i*123+457); // ditto
+
+        std::uniform_int_distribution<cell_gid_type> source_distribution(0, ncell_-2);
+
+        for (unsigned t=0; t<param_.num_synapses; ++t) {
+            auto source = source_distribution(source_gen);
+            if (source>=i) ++source;
+
+            cell_connection cc = draw_connection_params(conn_param_gen);
+            cc.source = {source, 0};
+            cc.dest = {i, t};
+            conns.push_back(cc);
+        }
+
+        return conns;
+    }
+};
+
+std::unique_ptr<recipe> make_basic_rgraph_recipe(
+        cell_gid_type ncell,
+        basic_recipe_param param,
+        probe_distribution pdist)
+{
+    return std::unique_ptr<recipe>(new basic_rgraph_recipe(ncell, param, pdist));
+}
+
+class basic_kgraph_recipe: public basic_cell_recipe {
+public:
+    basic_kgraph_recipe(cell_gid_type ncell,
+                      basic_recipe_param param,
+                      probe_distribution pdist = probe_distribution{}):
+        basic_cell_recipe(ncell, std::move(param), std::move(pdist))
+    {
+        if (std::size_t(param.num_synapses) != ncell-1) {
+            throw invalid_recipe_error("number of synapses per cell must equal number "
+                "of cells minus one in complete graph model");
+        }
+    }
+
+    std::vector<cell_connection> connections_on(cell_gid_type i) const override {
+        std::vector<cell_connection> conns;
+        auto conn_param_gen = std::mt19937(i); // TODO: replace this with hashing generator...
+
+        for (unsigned t=0; t<param_.num_synapses; ++t) {
+            cell_gid_type source = t>=i? t+1: t;
+            EXPECTS(source<ncell_);
+
+            cell_connection cc = draw_connection_params(conn_param_gen);
+            cc.source = {source, 0};
+            cc.dest = {i, t};
+            conns.push_back(cc);
+        }
+
+        return conns;
+    }
+};
+
+std::unique_ptr<recipe> make_basic_kgraph_recipe(
+        cell_gid_type ncell,
+        basic_recipe_param param,
+        probe_distribution pdist)
+{
+    return std::unique_ptr<recipe>(new basic_kgraph_recipe(ncell, param, pdist));
+}
+
+} // namespace mc
+} // namespace nest

miniapp/miniapp_recipes.hpp

+#pragma once
+
+#include <cstddef>
+#include <memory>
+#include <stdexcept>
+
+#include "recipe.hpp"
+
+// miniapp-specific recipes
+
+namespace nest {
+namespace mc {
+
+struct probe_distribution {
+    float proportion = 1.f; // what proportion of cells should get probes?
+    bool all_segments = true;    // false => soma only
+    bool membrane_voltage = true;
+    bool membrane_current = true;
+};
+
+struct basic_recipe_param {
+    unsigned num_compartments = 1;
+    unsigned num_synapses = 1;
+    std::string synapse_type = "expsyn";
+    float min_connection_delay_ms = 20.0;
+    float mean_connection_delay_ms = 20.75;
+    float syn_weight_per_cell = 0.3;
+};
+
+std::unique_ptr<recipe> make_basic_ring_recipe(
+        cell_gid_type ncell,
+        basic_recipe_param param,
+        probe_distribution pdist = probe_distribution{});
+
+std::unique_ptr<recipe> make_basic_kgraph_recipe(
+        cell_gid_type ncell,
+        basic_recipe_param param,
+        probe_distribution pdist = probe_distribution{});
+
+std::unique_ptr<recipe> make_basic_rgraph_recipe(
+        cell_gid_type ncell,
+        basic_recipe_param param,
+        probe_distribution pdist = probe_distribution{});
+
+} // namespace mc
+} // namespace nest

miniapp/trace_sampler.hpp

+#pragma once
+
+#include <cstdlib>
+#include <vector>
+
+#include <common_types.hpp>
+#include <cell.hpp>
+#include <util/optional.hpp>
+
+#include <iostream>
+
+namespace nest {
+namespace mc {
+
+template <typename Time=float, typename Value=double>
+struct sample_trace {
+    using time_type = Time;
+    using value_type = Value;
+
+    struct sample_type {
+        time_type time;
+        value_type value;
+    };
+
+    std::string name;
+    std::string units;
+    cell_member_type probe_id;
+    std::vector<sample_type> samples;
+
+    sample_trace() =default;
+    sample_trace(cell_member_type probe_id, const std::string& name, const std::string& units):
+        name(name), units(units), probe_id(probe_id)
+    {}
+};
+
+template <typename Time=float, typename Value=double>
+struct trace_sampler {
+    using time_type = Time;
+    using value_type = Value;
+
+    float next_sample_t() const { return t_next_sample_; }
+
+    util::optional<time_type> operator()(time_type t, value_type v) {
+        if (t<t_next_sample_) {
+            return t_next_sample_;
+        }
+
+        trace_->samples.push_back({t,v});
+        return t_next_sample_+=sample_dt_;
+    }
+
+    trace_sampler(sample_trace<time_type, value_type> *trace, time_type sample_dt, time_type tfrom=0):
+       trace_(trace), sample_dt_(sample_dt), t_next_sample_(tfrom)
+    {}
+
+private:
+    sample_trace<time_type, value_type> *trace_;
+
+    time_type sample_dt_;
+    time_type t_next_sample_;
+};
+
+// with type deduction ...
+template <typename Time, typename Value>
+trace_sampler<Time, Value> make_trace_sampler(sample_trace<Time, Value> *trace, Time sample_dt, Time tfrom=0) {
+    return trace_sampler<Time, Value>(trace, sample_dt, tfrom);
+}
+
+} // namespace mc
+} // namespace nest

nrn/generate_validation.sh

-python2.7 ./soma.py
-python2.7 ./ball_and_stick.py
-python2.7 ./ball_and_3stick.py
-python2.7 ./simple_synapse.py --synapse exp2
-python2.7 ./simple_synapse.py --synapse exp
+#!/usr/bin/env bash
+
+python2 ./soma.py
+python2 ./ball_and_stick.py
+python2 ./ball_and_3stick.py
+python2 ./simple_synapse.py --synapse exp2
+python2 ./simple_synapse.py --synapse exp

src/CMakeLists.txt

@@ -2,8 +2,8 @@ set(HEADERS
     swcio.hpp
 )
 set(BASE_SOURCES
+    common_types_io.cpp
     cell.cpp
-    mechanism_interface.cpp
     parameter_list.cpp
     profiling/profiler.cpp
     swcio.cpp

src/algorithms.hpp

@@ -90,12 +90,12 @@ bool is_minimal_degree(C const& c)
         "is_minimal_degree only applies to integral types"
     );
 
-    if(c.size()==0u) {
+    if (c.size()==0u) {
         return true;
     }
 
     using value_type = typename C::value_type;
-    if(c[0] != value_type(0)) {
+    if (c[0] != value_type(0)) {
         return false;
     }
     auto i = value_type(1);
@@ -121,49 +121,43 @@ bool is_positive(C const& c)
 }
 
 template<typename C>
-bool has_contiguous_segments(const C &parent_index)
+std::vector<typename C::value_type> child_count(const C& parent_index)
 {
     static_assert(
         std::is_integral<typename C::value_type>::value,
         "integral type required"
     );
 
-    if (!is_minimal_degree(parent_index)) {
-        return false;
-    }
-
-    int n = parent_index.size();
-    std::vector<bool> is_leaf(n, false);
-
-    for(auto i=1; i<n; ++i) {
-        auto p = parent_index[i];
-        if(is_leaf[p]) {
-            return false;
-        }
-
-        if(p != i-1) {
-            // we have a branch and i-1 is a leaf node
-            is_leaf[i-1] = true;
-        }
+    std::vector<typename C::value_type> count(parent_index.size(), 0);
+    for (auto i = 1u; i < parent_index.size(); ++i) {
+        ++count[parent_index[i]];
     }
 
-    return true;
+    return count;
 }
 
 template<typename C>
-std::vector<typename C::value_type> child_count(const C &parent_index)
+bool has_contiguous_compartments(const C& parent_index)
 {
+    using value_type = typename C::value_type;
     static_assert(
-        std::is_integral<typename C::value_type>::value,
+        std::is_integral<value_type>::value,
         "integral type required"
     );
 
-    std::vector<typename C::value_type> count(parent_index.size(), 0);
-    for (std::size_t i = 1; i < parent_index.size(); ++i) {
-        ++count[parent_index[i]];
+    if (!is_minimal_degree(parent_index)) {
+        return false;
     }
 
-    return count;
+    auto num_child = child_count(parent_index);
+    for (auto i=1u; i < parent_index.size(); ++i) {
+        auto p = parent_index[i];
+        if (num_child[p]==1 && p!=value_type(i-1)) {
+            return false;
+        }
+    }
+
+    return true;
 }
 
 template<typename C>
@@ -174,7 +168,7 @@ std::vector<typename C::value_type> branches(const C& parent_index)
         "integral type required"
     );
 
-    EXPECTS(has_contiguous_segments(parent_index));
+    EXPECTS(has_contiguous_compartments(parent_index));
 
     std::vector<typename C::value_type> branch_index;
     if (parent_index.empty()) {
@@ -229,7 +223,7 @@ typename C::value_type find_branch(const C& branch_index,
 
     auto it =  std::find_if(
         branch_index.begin(), branch_index.end(),
-        [nid](const value_type &v) { return v > nid; }
+        [nid](const value_type& v) { return v > nid; }
     );
 
     return it - branch_index.begin() - 1;
@@ -249,8 +243,8 @@ std::vector<typename C::value_type> make_parent_index(
         return {};
     }
 
-    EXPECTS(parent_index.size() == unsigned(branch_index.back()));
-    EXPECTS(has_contiguous_segments(parent_index));
+    EXPECTS(parent_index.size() == branch_index.back());
+    EXPECTS(has_contiguous_compartments(parent_index));
     EXPECTS(is_strictly_monotonic_increasing(branch_index));
 
     // expand the branch index

src/cell.cpp

@@ -26,7 +26,7 @@ cell::cell()
     parents_.push_back(0);
 }
 
-int cell::num_segments() const
+cell::size_type cell::num_segments() const
 {
     return segments_.size();
 }
@@ -75,9 +75,9 @@ cable_segment* cell::add_cable(cell::index_type parent, segment_ptr&& cable)
     return segments_.back()->as_cable();
 }
 
-segment* cell::segment(int index)
+segment* cell::segment(index_type index)
 {
-    if(index<0 || index>=num_segments()) {
+    if (index>=num_segments()) {
         throw std::out_of_range(
             "attempt to access a segment with invalid index"
         );
@@ -85,9 +85,9 @@ segment* cell::segment(int index)
     return segments_[index].get();
 }
 
-segment const* cell::segment(int index) const
+segment const* cell::segment(index_type index) const
 {
-    if(index<0 || index>=num_segments()) {
+    if (index>=num_segments()) {
         throw std::out_of_range(
             "attempt to access a segment with invalid index"
         );
@@ -109,7 +109,7 @@ soma_segment* cell::soma()
     return nullptr;
 }
 
-cable_segment* cell::cable(int index)
+cable_segment* cell::cable(index_type index)
 {
     if(index>0 && index<num_segments()) {
         return segment(index)->as_cable();
@@ -146,9 +146,9 @@ std::vector<segment_ptr> const& cell::segments() const
     return segments_;
 }
 
-std::vector<int> cell::compartment_counts() const
+std::vector<cell::size_type> cell::compartment_counts() const
 {
-    std::vector<int> comp_count;
+    std::vector<size_type> comp_count;
     comp_count.reserve(num_segments());
     for(auto const& s : segments()) {
         comp_count.push_back(s->num_compartments());
@@ -156,10 +156,10 @@ std::vector<int> cell::compartment_counts() const
     return comp_count;
 }
 
-size_t cell::num_compartments() const
+cell::size_type cell::num_compartments() const
 {
     auto n = 0u;
-    for(auto &s : segments_) {
+    for(auto& s : segments_) {
         n += s->num_compartments();
     }
     return n;
@@ -196,7 +196,7 @@ void cell::add_detector(segment_location loc, double threshold)
     spike_detectors_.push_back({loc, threshold});
 }
 
-std::vector<int> const& cell::segment_parents() const
+std::vector<cell::index_type> const& cell::segment_parents() const
 {
     return parents_;
 }
@@ -212,24 +212,22 @@ std::vector<int> const& cell::segment_parents() const
 //  - number of compartments in each segment
 bool cell_basic_equality(cell const& lhs, cell const& rhs)
 {
-    if(lhs.num_segments() != rhs.num_segments()) {
+    if (lhs.num_segments() != rhs.num_segments()) {
         return false;
     }
-    if(lhs.segment_parents() != rhs.segment_parents()) {
+    if (lhs.segment_parents() != rhs.segment_parents()) {
         return false;
     }
-    for(auto i=0; i<lhs.num_segments(); ++i) {
+    for (cell::index_type i=0; i<lhs.num_segments(); ++i) {
         // a quick and dirty test
         auto& l = *lhs.segment(i);
         auto& r = *rhs.segment(i);
 
-        if(l.kind() != r.kind()) return false;
-        if(l.area() != r.area()) return false;
-        if(l.volume() != r.volume()) return false;
-        if(l.as_cable()) {
-            if(   l.as_cable()->num_compartments()
-               != r.as_cable()->num_compartments())
-            {
+        if (l.kind() != r.kind()) return false;
+        if (l.area() != r.area()) return false;
+        if (l.volume() != r.volume()) return false;
+        if (l.as_cable()) {
+            if (l.as_cable()->num_compartments() != r.as_cable()->num_compartments()) {
                 return false;
             }
         }