Arbor cable cell exporter and backend support in BluePyOpt

arbor/include/arbor/morph/morphology.hpp

@@ -41,6 +41,11 @@ class ARB_ARBOR_API morphology {
     // Range of segments in a branch.
     const std::vector<msegment>& branch_segments(msize_t b) const;
 
+    // Convert to segment tree
+    // Note: tree == arb::morphology(tree).to_segment_tree() is not guaranteed
+    // to be true.
+    segment_tree to_segment_tree() const;
+
     friend std::ostream& operator<<(std::ostream&, const morphology&);
 };
 

arbor/include/arbor/morph/segment_tree.hpp

@@ -85,4 +85,8 @@ equivalent(const segment_tree& a,
 ARB_ARBOR_API segment_tree
 apply(const segment_tree&, const isometry&);
 
-} // namesapce arb
+// Roots of regions of specific tag in segment tree
+ARB_ARBOR_API std::vector<msize_t> tag_roots(const segment_tree& in, int tag);
+
+
+} // namespace arb

arbor/morph/morphology.cpp

@@ -173,6 +173,29 @@ msize_t morphology::num_branches() const {
     return impl_->branches_.size();
 }
 
+ARB_ARBOR_API segment_tree morphology::to_segment_tree() const {
+    segment_tree st;
+    const auto& branches = impl_->branches_;
+
+    for (auto bi: make_span(branches.size()) ) {
+        const auto& branch = branches[bi];
+        for (auto si: make_span(branch.size())) {
+            const auto& seg = branch[si];
+            auto p = mnpos;
+            if (si == 0) {
+                auto bp = impl_->branch_parents_[bi];
+                p = bp == mnpos ? mnpos : branches[bp].back().id;
+            } else {
+                p = branch[si-1].id;
+            }
+            st.append(p, seg.prox, seg.dist, seg.tag);
+        }
+    }
+
+    return st;
+}
+
+
 ARB_ARBOR_API std::ostream& operator<<(std::ostream& o, const morphology& m) {
     return o << *m.impl_;
 }

arbor/morph/segment_tree.cpp

@@ -10,16 +10,18 @@
 #include "util/span.hpp"
 #include "util/transform.hpp"
 
+using arb::util::make_span;
+
 namespace arb {
 
 struct node_t {
     msize_t parent;
     msize_t id;
 };
 
-using node_p = std::function<bool(const node_t&)>;
+using id_p = std::function<bool(msize_t)>;
 
-node_p yes = [](node_t) { return true; };
+id_p yes = [](msize_t) { return true; };
 
 // invert parent <*> child relation, returns a map of parent_id -> [children_id]
 // For predictable ordering we sort the vectors.
@@ -30,52 +32,81 @@ std::map<msize_t, std::vector<msize_t>> tree_to_children(const segment_tree& tre
     for (auto& [k, v]: result) std::sort(v.begin(), v.end());
     return result;
 }
-    
-// Copy a segment tree into a new tree
+
+// Copy a segment subtree into a new tree
 // - tree to be (partially) copied
-// - start={parent, id}: start by attaching segment=`id`` from `tree` to the
-//                       output at `parent`, then its children to it recursively
+// - start={parent, id}: start by attaching segment=`id` from `tree` to the segment
+//                       `parent` of `init`, then its children to it recursively
 // - predicate: if returning false for a given node, we prune that sub-tree starting
-//              at node (inclusive). Can be used to prune trees by parent or id.
+//              at node (inclusive). Can be used to prune trees by segment id.
 // - init: initial tree to append to
-// Note: this is basically a recursive function w/ an explicit stack.
-segment_tree copy_if(const segment_tree& tree,
-                     node_t start,
-                     node_p predicate,
-                     const segment_tree& init={}) {
+// Note: this is an iterative implementation of depth-first traversal with an explicit stack.
+std::pair<segment_tree, std::vector<bool>> copy_subtree_if(const segment_tree& tree,
+                                                           const node_t& start,
+                                                           id_p predicate,
+                                                           const segment_tree& init={}) {
     auto children_of = tree_to_children(tree);
     auto& segments = tree.segments();
     segment_tree result = init;
+    std::vector<bool> copied(segments.size(), false);
     auto todo = std::vector<node_t>{start};
     while (!todo.empty()) {
         auto node = todo.back();
         todo.pop_back();
-        if (!predicate(node)) continue;
+        if (!predicate(node.id)) continue;
+        copied[node.id] = true;
         const auto& segment = segments[node.id];
         auto current = result.append(node.parent, segment.prox, segment.dist, segment.tag);
         for (auto child: children_of[node.id]) {
             todo.push_back({current, child});
         }
     }
-    return result;
+    return {result, copied};
+}
+
+// Filter a segment tree as a graph
+// - tree to be (partially) filtered
+// - predicate: if returning false for a given node, we exclude that node from the
+//              returned tree. Can be used to filter trees by segment id.
+std::pair<segment_tree, std::vector<bool>> copy_fulltree_if(const segment_tree& tree,
+                                                            id_p predicate) {
+    auto& segments = tree.segments();
+    auto& parents = tree.parents();
+
+    segment_tree result;
+    std::vector<bool> copied(segments.size(), false);
+    std::vector<msize_t> new_ids(segments.size(), mnpos);
+
+    for (msize_t id: make_span(0, segments.size())) {
+        if (!predicate(id)) continue;
+        copied[id] = true;
+        const auto& segment = segments[id];
+        auto parent = parents[id];
+        new_ids[id] = result.append(parent != mnpos ? new_ids[parent] : mnpos,
+                                    segment.prox, segment.dist, segment.tag);
+    }
+    return {result, copied};
 }
 
 std::pair<segment_tree, segment_tree>
 split_at(const segment_tree& tree, msize_t at) {
     if (at >= tree.size() || at == mnpos) throw invalid_segment_parent(at, tree.size());
     // span the sub-tree starting at the splitting node
-    segment_tree post = copy_if(tree, {mnpos, at}, yes);
-    // copy the original tree, but skip all nodes in the `post` tree
-    segment_tree pre = copy_if(tree,
-                               {mnpos, 0},
-                               [=](auto& node) { return node.id != at; });
+    const auto& post_copied = copy_subtree_if(tree, {mnpos, at}, yes);
+    auto post = post_copied.first;
+    auto copied = post_copied.second;
+
+    // copy the original segment_tree (as a graph), skipping all nodes in the `post` subtree
+    segment_tree pre = copy_fulltree_if(tree,
+                                        [&](msize_t id) { return !copied[id]; }).first;
+
     return {pre, post};
 }
 
 segment_tree
 join_at(const segment_tree& lhs, msize_t at, const segment_tree& rhs) {
     if (at >= lhs.size() && at != mnpos) throw invalid_segment_parent(at, lhs.size());
-    return copy_if(rhs, {at, 0}, yes, lhs);
+    return copy_subtree_if(rhs, {at, 0}, yes, lhs).first;
 }
 
 bool
@@ -98,7 +129,7 @@ equivalent(const segment_tree& a,
         return segs;
     };
 
-    std::vector<std::pair<arb::msize_t, arb::msize_t>> todo{{arb::mnpos, arb::mnpos}};
+    std::vector<std::pair<msize_t, msize_t>> todo{{arb::mnpos, arb::mnpos}};
     while (!todo.empty()) {
         const auto& [a_cursor, b_cursor] = todo.back();
         auto as = fetch_children(a_cursor, a.segments(), a_children_of);
@@ -200,5 +231,22 @@ ARB_ARBOR_API std::ostream& operator<<(std::ostream& o, const segment_tree& m) {
              << io::sepval(tstr, ' ') <<  "))";
 }
 
+
+ARB_ARBOR_API std::vector<msize_t> tag_roots(const segment_tree& t, int tag) {
+    const auto& segments = t.segments();
+    const auto& parents = t.parents();
+    std::vector<msize_t> tag_roots;
+
+    for (auto i: make_span(0, segments.size())) {
+        auto par = parents[i];
+        if (segments[i].tag == tag && (par == mnpos || segments[par].tag != tag)) {
+            tag_roots.push_back(i);
+        }
+    }
+
+    return tag_roots;
+}
+
+
 } // namespace arb
 

doc/concepts/decor.rst

@@ -342,13 +342,13 @@ A point mechanism (synapse) can form the target of a :term:`connection` on a cel
     expsyn = arbor.mechanism('expsyn')
 
     # Wrap the 'expsyn' mechanism in a `synapse` object and add it to the decor.
-    decor.paint('"syn_loc_0"', arbor.synapse(expsyn))
+    decor.place('"syn_loc_0"', arbor.synapse(expsyn))
 
     # Create an 'expsyn' mechanism with default parameter values as a `synapse` object, and add it to the decor.
-    decor.paint('"syn_loc_1"', arbor.synapse("expsyn"))
+    decor.place('"syn_loc_1"', arbor.synapse("expsyn"))
 
     # Create an 'expsyn' mechanism with modified 'tau' parameter as a `synapse` object, and add it to the decor.
-    decor.paint('"syn_loc_2"', arbor.synapse("expsyn", {"tau": 1.0}))
+    decor.place('"syn_loc_2"', arbor.synapse("expsyn", {"tau": 1.0}))
 
 
 .. _cablecell-threshold-detectors:
@@ -399,10 +399,10 @@ on two separate cells.
     gj = arbor.mechanism("gj", {"g": 2.0})
 
     # Wrap the 'gj' mechanism in a `junction` object and add it to the decor.
-    decor.paint('"gj_loc_0"', arbor.junction(gj))
+    decor.place('"gj_loc_0"', arbor.junction(gj))
 
     # Create a 'gj' mechanism with modified 'g' parameter as a `junction` object, and add it to the decor.
-    decor.paint('"gj_loc_1"', arbor.junction("gj", {"g": 1.5}))
+    decor.place('"gj_loc_1"', arbor.junction("gj", {"g": 1.5}))
 
 .. _cablecell-stimuli:
 
@@ -431,14 +431,14 @@ constant stimuli and constant amplitude stimuli restricted to a fixed time inter
     decor.place('(root)', arbor.iclamp(10), "iclamp0")
 
     # Constant amplitude 10 nA stimulus at 20 Hz, with initial phase of π/4 radians.
-    decor.place('(root)', arbor.iclamp(10, frequency=0.020, phasce=math.pi/4), "iclamp1")
+    decor.place('(root)', arbor.iclamp(10, frequency=0.020, phase=math.pi/4), "iclamp1")
 
     # Stimulus at 1 kHz, amplitude 10 nA, for 40 ms starting at t = 30 ms.
     decor.place('(root)', arbor.iclamp(30, 40, 10, frequency=1), "iclamp2")
 
     # Piecewise linear stimulus with amplitude ranging from 0 nA to 10 nA,
     # starting at t = 30 ms and stopping at t = 50 ms.
-    decor.place('(root)', arbor.iclamp([(30, 0), (37, 10), (43, 8), (50, 0)], "iclamp3")
+    decor.place('(root)', arbor.iclamp([(30, 0), (37, 10), (43, 8), (50, 0)]), "iclamp3")
 
 
 .. _cablecell-probes:

doc/concepts/interconnectivity.rst

@@ -68,7 +68,7 @@ full recipe.
 
    Arbor uses a lazily constructed network (from the ``recipe`` callbacks) for
    good reason; storing the full connectivity (for all ``gids``) in the
-   ``recipe`` can lead to probitivively large memory footprints. Keep this in
+   ``recipe`` can lead to prohibitively large memory footprints. Keep this in
    mind when designing your connectivity and heed the consequences of doing I/O
    in these callbacks. This is doubly important when using models with dynamic
    connectivity where the temptation to store all connections is even larger and

doc/concepts/labels.rst

@@ -300,7 +300,7 @@ Locset expressions
 .. label:: (proximal-translate ls:locset distance:real)
 
     The set of locations that correspond to moving each location in the ``ls`` in the proximal direction
-    ``distance`` μm. The locations in the output have a one to one correspondance with those in ``ls``.
+    ``distance`` μm. The locations in the output have a one to one correspondence with those in ``ls``.
 
     .. figure:: ../gen-images/proximal_translate_label.svg
       :width: 600
@@ -316,7 +316,7 @@ Locset expressions
 
     An input location will generate multiple output locations when it is translated
     past a fork point, with a new location for each child branch (see the example
-    below). For this reason there is not a one-to-one correspondance between locations
+    below). For this reason there is not a one-to-one correspondence between locations
     in the input and output sets, so the results are sorted and duplicates are removed.
 
 
@@ -862,7 +862,7 @@ the regions, locsets or iexpr defined in the label dictionary by referring to th
       'root': '(root)',   # typically the start of the soma is at the root of the cell.
       'stim_site': '(location 0 0.5)', # site for the stimulus, in the middle of branch 0.
       'axon_end': '(restrict (terminal) (region "axon"))',  # end of the axon.
-      'rad_expr': '(radius 0.5)'}  # iexpr evaluating the radius scaled by 0.5
+      'rad_expr': '(radius 0.5)'  # iexpr evaluating the radius scaled by 0.5
     })
 
 

doc/concepts/morphology.rst

@@ -607,7 +607,7 @@ Editing morphologies
 
 While a reified morphology cannot be edited -- it is immutable by definition --
 the segment tree can be changed. If you need to make such modifications, first
-consider whethnr they should be stored in a file as this is often easier for
+consider whether they should be stored in a file as this is often easier for
 tracking provenance and version history.
 
 For the remaining cases, Arbor offers a limited suite of tools. First, most
@@ -632,6 +632,13 @@ Two common editing operations are provided
 
 Note that ``join_at`` and ``split_at`` are inverse to each other.
 
+A particular use-case for these operations is pruning a specific tag-region in the
+segment tree and replacing it with a surrogate model. This is e.g. commonly performed
+for the axon, known as axon-replacement. For this purpose, the function ``tag_roots``
+allows to obtain the IDs of root segments of a a tag region. These IDs can
+then be used with ``split_at`` to split off subtrees and ``join_at`` to attach a
+surrogate subtree (with the same parent as the one split off).
+
 .. _morph-formats:
 
 Supported file formats

doc/cpp/hardware.rst

@@ -90,7 +90,7 @@ The header ``arborenv/concurrency.hpp`` supplies lower-level functions for query
    Returns the list of logical processor ids where the calling thread has affinity,
    or an empty vector if unable to determine.
 
-The header ``arborenv/gpu_env.hpp`` supplies lower-level functions for queruing the GPU environment.
+The header ``arborenv/gpu_env.hpp`` supplies lower-level functions for querying the GPU environment.
 
 .. cpp:function:: int find_private_gpu(MPI_Comm comm)
 

doc/cpp/morphology.rst

@@ -77,6 +77,11 @@ consistent parent-child indexing, and with ``n`` segments numbered from ``0`` to
     Join two subtrees at a given id, such that said id becomes the parent
     of the inserted sub-tree.
 
+.. cpp:function:: std::vector<msize_t> tag_roots(const segment_tree& t, int tag)
+
+    Get IDs of roots of a region with specific tag in the segment tree, i.e. segments whose
+    parent is either :data:`mnpos` or a segment with a different tag.
+
 .. cpp:function:: bool equivalent(const segment_tree& l, const segment_tree& r)
 
     Two trees are equivalent if

doc/dev/simd_api.rst

@@ -1060,7 +1060,7 @@ wrap calls to the Intel scalar vector mathematics library (SVML).
 
 Outside of this case, the functions *exp*, *log*, *expm1* and
 *exprelr* use explicit approximations as detailed below. The
-algortihms follow those used in the
+algorithms follow those used in the
 `Cephes library <http://www.netlib.org/cephes/>`_, with
 some accommodations.
 

doc/python/morphology.rst

@@ -254,6 +254,12 @@ Cable cell morphology
         ``join_at`` is inverse to ``split_at`` for a proper choice of ``id``.
         The join point ``id`` must be in ``L``.
 
+    .. method:: tag_roots(tag)
+
+        Get IDs of roots of region with a particular ``tag`` in the segment tree, i.e.
+        segments whose parent is either :data:`mnpos` or a segment with a different
+        tag.
+
     .. method:: apply_isometry(iso)
 
         Apply an :py:class:`isometry` to the segment tree, returns the transformed tree as a copy.

doc/tutorial/index.rst

@@ -24,6 +24,7 @@ Cells
    single_cell_detailed
    single_cell_cable
    single_cell_allen
+   single_cell_bluepyopt
 
 Recipes
 -------

doc/tutorial/network_ring_gpu.rst

@@ -69,7 +69,7 @@ hardware Arbor will execute on.
 Step **(14)** creates a :class:`arbor.partition_hint`, and tells it to put 1000 cells in a groups allocated to GPUs, 
 and to prefer the utilisation of the GPU if present. In fact, the default distribution strategy of 
 :class:`arbor.partition_load_balance` already spreads out cells as evenly as possible over CPUs, and groups
-(up to 1000) on GPUs, so strictly speaking it was not necesary to give that part of the hint.
+(up to 1000) on GPUs, so strictly speaking it was not necessary to give that part of the hint.
 Lastly, a dictionary is created with which hints are assigned to a particular :class:`arbor.cell_kind`.
 Different kinds may favor different execution, hence the option.
 In this simulation, there are only :class:`arbor.cell_kind.cable`, so we assign the hint to that kind.

doc/tutorial/network_ring_mpi.rst

@@ -75,7 +75,7 @@ The results
 
 Before we execute the simulation, we have to understand how Arbor distributes the computational load over the ranks.
 After executing ``mpirun``, all nodes will run the same script. In the domain decomposition step, the nodes will use
-the provided MPI communicator to divide the work. Once :py:func:`arbor.simulation.run` starts, each node wil work on
+the provided MPI communicator to divide the work. Once :py:func:`arbor.simulation.run` starts, each node will work on
 their allocated cell ``gid`` s.
 
 This is relevant for the collection of results: these are not gathered for you. Remember that in step **(15)** we

doc/tutorial/probe_lfpykit.rst

@@ -18,7 +18,7 @@ These are signals that mainly stem from transmembrane currents.
    1. Recording of transmembrane currents using :class:`arbor.cable_probe_total_current_cell`
    2. Recording of stimulus currents using :class:`arbor.cable_probe_stimulus_current_cell`
    3. Using the :class:`arbor.place_pwlin` API
-   4. Map recorded transmembrane currents to extracellular potentials by derivinging Arbor specific classes from LFPykit's :class:`lfpykit.LineSourcePotential` and :class:`lfpykit.CellGeometry`
+   4. Map recorded transmembrane currents to extracellular potentials by deriving Arbor specific classes from LFPykit's :class:`lfpykit.LineSourcePotential` and :class:`lfpykit.CellGeometry`
 
 .. _tutorial_lfpykit-linesource: