diff --git a/src/algorithms.hpp b/src/algorithms.hpp index da3de60f6d7489b028b1cd1903d58c3095365650..88f6ff94cc40dc7bf66c2cfab87a73789647f9c6 100644 --- a/src/algorithms.hpp +++ b/src/algorithms.hpp @@ -242,9 +242,33 @@ typename C::value_type find_branch(const C& branch_index, return it - branch_index.begin() - 1; } - +// Find the reduced form of a tree represented as a parent index. +// The operation of transforming a tree into a homeomorphic pre-image is called +// 'reduction'. The 'reduced' tree is the smallest tree that maps into the tree +// homeomorphically (i.e. preserving labels and child relationships). +// For example, the tree represented by the following index: +// {0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12} +// Has reduced form +// {0, 0, 0, 0, 3, 3} +// +// This transformation can be represented graphically: +// +// 0 0 +// /|\ /|\. +// 1 4 6 => 1 2 3 +// / | \ / \. +// 2 5 7 4 5 +// / \. +// 3 8 +// / \. +// 9 11 +// / \. +// 10 12 +// \. +// 13 +// template<typename C> -std::vector<typename C::value_type> make_parent_index( +std::vector<typename C::value_type> tree_reduce( const C& parent_index, const C& branch_index) { static_assert( diff --git a/src/cell.cpp b/src/cell.cpp index 5dfb1920b96cb534c5fa26a2d6ae2a91f952d654..12380329636119a67801372158e450b4d2dc07a5 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -169,9 +169,9 @@ compartment_model cell::model() const { compartment_model m; - m.tree = cell_tree(parents_); + m.tree = tree(parents_); auto counts = compartment_counts(); - m.parent_index = make_parent_index(m.tree.graph(), counts); + m.parent_index = make_parent_index(m.tree, counts); m.segment_index = algorithms::make_index(counts); return m; diff --git a/src/cell.hpp b/src/cell.hpp index 28149edf897c096e990ab341a443d3e13f1b39e2..200691b73bb1d148e09e1755aebe6ce7083f24d3 100644 --- a/src/cell.hpp +++ b/src/cell.hpp @@ -7,7 +7,7 @@ #include <vector> #include <common_types.hpp> -#include <cell_tree.hpp> +#include <tree.hpp> #include <morphology.hpp> #include <segment.hpp> #include <stimulus.hpp> @@ -20,9 +20,9 @@ namespace arb { /// wrapper around compartment layout information derived from a high level cell /// description struct compartment_model { - cell_tree tree; - std::vector<cell_tree::int_type> parent_index; - std::vector<cell_tree::int_type> segment_index; + arb::tree tree; + std::vector<tree::int_type> parent_index; + std::vector<tree::int_type> segment_index; }; int find_compartment_index( diff --git a/src/cell_tree.hpp b/src/cell_tree.hpp deleted file mode 100644 index 14d232d0711a124a4bcd5a9b3e5b903e09b2a866..0000000000000000000000000000000000000000 --- a/src/cell_tree.hpp +++ /dev/null @@ -1,302 +0,0 @@ -#pragma once - -#include <algorithm> -#include <iomanip> -#include <iostream> -#include <fstream> -#include <memory> -#include <numeric> -#include <ostream> -#include <vector> - -#include <memory/memory.hpp> -#include <util/span.hpp> -#include <common_types.hpp> -#include <tree.hpp> - -namespace arb { - -/// The tree data structure that describes the segments of a cell tree. -/// A cell is represented as a tree where each node may have any number of -/// children. Typically in a cell only the soma has more than two segments, -/// however this rule does not appear to be strictly followed in the PCP data -/// sets. -/// -/// To optimize some computations it is important the the tree be balanced, -/// in the sense that the depth of the tree be minimized. This means that it is -/// necessary that any node in the tree may be used as the root. In the PCP data -/// sets it appears that the soma was always index 0, however we need more -/// flexibility in choosing the root. -class cell_tree { -public: - using int_type = cell_lid_type; - using size_type = cell_local_size_type; - - using iarray = memory::host_vector<int_type>; - using view_type = iarray::view_type; - using const_view_type = iarray::const_view_type; - - using tree = arb::tree<int_type, size_type>; - static constexpr int_type no_parent = tree::no_parent; - - /// default empty constructor - cell_tree() = default; - - /// construct from a parent index - cell_tree(std::vector<int_type> const& parent_index) - { - // handle the case of an empty parent list, which implies a single-compartment model - if(parent_index.size()>0) { - tree_ = tree(parent_index); - } - else { - tree_ = tree(std::vector<int_type>({0})); - } - } - - /// construct from a tree - // copy constructor - cell_tree(tree const& t, int_type s) - : tree_(t), - soma_(s) - { } - - // move constructor - cell_tree(tree&& t, int_type s) - : tree_(std::move(t)), - soma_(s) - { } - - /// construct from a cell tree - // copy constructor - cell_tree(cell_tree const& other) - : tree_(other.tree_), - soma_(other.soma()) - { } - - // assignment from rvalue - cell_tree& operator=(cell_tree&& other) - { - std::swap(other.tree_, tree_); - std::swap(other.soma_, soma_); - return *this; - } - - // assignment - cell_tree& operator=(cell_tree const& other) - { - tree_ = other.tree_; - soma_ = other.soma_; - return *this; - } - - // move constructor - cell_tree(cell_tree&& other) - { - *this = std::move(other); - } - - tree const& graph() const { - return tree_; - } - - int_type soma() const { - return soma_; - } - - /// Minimize the depth of the tree. - int_type balance() { - // find the new root - auto new_root = find_minimum_root(); - - // change the root on the tree - auto p = tree_.change_root(new_root); - - // keep track of the soma_ - if(p.size()) { - soma_ = p[soma_]; - } - - return new_root; - } - - /// memory used to store cell tree (in bytes) - size_t memory() const { - return tree_.memory() + sizeof(cell_tree) - sizeof(tree); - } - - /// returns the number of segments in the cell - size_t num_segments() const { - return tree_.num_nodes(); - } - - /// returns the number of child segments of segment b - size_type num_children(int_type b) const { - return tree_.num_children(b); - } - - /// returns a list of the children of segment b - const_view_type children(int_type b) const { - return tree_.children(b); - } - - /// returns the parent index of segment b - int_type parent(size_t b) const { - return tree_.parent(b); - } - - /// writes a graphviz .dot file that visualizes cell segment structure - void to_graphviz(std::string const& fname) const { - - std::ofstream fid(fname); - - fid << "graph cell {" << '\n'; - for(auto b : util::make_span(0,num_segments())) { - if(children(b).size()) { - for(auto c : children(b)) { - fid << " " << b << " -- " << c << ";" << '\n'; - } - } - } - fid << "}" << std::endl; // flush at end of output? - } - - iarray depth_from_leaf() - { - tree::iarray depth(num_segments()); - depth_from_leaf(depth, int_type{0}); - return depth; - } - - iarray depth_from_root() - { - tree::iarray depth(num_segments()); - depth[0] = 0; - for (auto c: children(0)) { - depth_from_root(depth, c); - } - return depth; - } - -private : - - /// helper type for sub-tree computation - /// use in balance() - struct sub_tree { - sub_tree(int_type r, int_type diam, int_type dpth): - root(r), diameter(diam), depth(dpth) - {} - - void set(int r, int diam, int dpth) { - root = r; - diameter = diam; - depth = dpth; - } - - std::string to_string() const { - return - "[" + std::to_string(root) + "," - + std::to_string(diameter) + "," - + std::to_string(depth) + - "]"; - } - - int_type root; - int_type diameter; - int_type depth; - }; - - /// returns the index of the segment that would minimise the depth of the - /// tree if used as the root segment - int_type find_minimum_root() { - if (num_segments()==1) { - return 0; - } - - // calculate the depth of each segment from the root - // pre-order traversal of the tree - auto depth = depth_from_root(); - - auto max = std::max_element(depth.begin(), depth.end()); - auto max_leaf = std::distance(depth.begin(), max); - - // Calculate the depth of each compartment as the maximum distance - // from a child leaf - // post-order traversal of the tree - depth = depth_from_leaf(); - - // Walk back from the deepest leaf towards the root. - // At each node test to find the deepest sub-tree that doesn't include - // node max_leaf. Then check if the total diameter of this sub-tree and - // the sub-tree and node max_leaf is the largest found so far. When the - // walk has been completed to the root node, the node that has been - // selected will be the root of the sub-tree with the largest diameter. - sub_tree max_sub_tree(0, 0, 0); - int_type distance_from_max_leaf = 1; - auto pnt = max_leaf; - auto pos = parent(max_leaf); - while(pos != no_parent) { - for(auto c : children(pos)) { - if(c!=pnt) { - auto diameter = depth[c] + 1 + distance_from_max_leaf; - if (diameter>max_sub_tree.diameter) { - max_sub_tree.set(pos, diameter, distance_from_max_leaf); - } - } - } - pnt = pos; - pos = parent(pos); - ++distance_from_max_leaf; - } - - // calculate the depth of the balanced tree - auto new_depth = (max_sub_tree.diameter+1) / 2; - - // nothing to do if the current root is also the root of the - // balanced tree - if(max_sub_tree.root==0 && max_sub_tree.depth==new_depth) { - return 0; - } - - // perform another walk from max leaf towards max_sub_tree.root - auto count = new_depth; - auto new_root = max_leaf; - while(count) { - new_root = parent(new_root); - --count; - } - - return new_root; - } - - int_type depth_from_leaf(iarray& depth, int_type segment) - { - int_type max_depth = 0; - for(auto c : children(segment)) { - max_depth = std::max(max_depth, depth_from_leaf(depth, c)); - } - depth[segment] = max_depth; - return max_depth+1; - } - - void depth_from_root(iarray& depth, int_type segment) - { - auto d = depth[parent(segment)] + 1; - depth[segment] = d; - for(auto c : children(segment)) { - depth_from_root(depth, c); - } - } - - ////////////////////////////////////////////////// - // state - ////////////////////////////////////////////////// - - /// storage for the tree structure of cell segments - tree tree_; - - /// index of the soma - int_type soma_ = 0; -}; - -} // namespace arb diff --git a/src/swcio.hpp b/src/swcio.hpp index 2ba4b3943d963085841d49097e389b14a1850e73..ee347995e896ed66f235d959f4f6e9f08a49bf4c 100644 --- a/src/swcio.hpp +++ b/src/swcio.hpp @@ -122,7 +122,7 @@ morphology swc_as_morphology(const RandomAccessSequence& swc_records) { // The parent of soma must be 0, while in SWC files is -1 swc_parent_index[0] = 0; auto branch_index = algorithms::branches(swc_parent_index); // partitions [0, #records] by branch. - auto parent_branch_index = algorithms::make_parent_index(swc_parent_index, branch_index); + auto parent_branch_index = algorithms::tree_reduce(swc_parent_index, branch_index); // sanity check EXPECTS(parent_branch_index.size() == branch_index.size() - 1); diff --git a/src/tree.hpp b/src/tree.hpp index 7d9eab05ffb8f32c784bdb151b95feacca822f37..96f90043bc81d57afd503fc1244927233770957f 100644 --- a/src/tree.hpp +++ b/src/tree.hpp @@ -4,18 +4,18 @@ #include <cassert> #include <numeric> #include <vector> -#include <memory/memory.hpp> -#include <util/span.hpp> #include <algorithms.hpp> +#include <common_types.hpp> +#include <memory/memory.hpp> +#include <util/span.hpp> namespace arb { -template <typename Int, typename Size = std::size_t> class tree { public: - using int_type = Int; - using size_type = Size; + using int_type = cell_lid_type; + using size_type = cell_local_size_type; using iarray = memory::host_vector<int_type>; using view_type = typename iarray::view_type; @@ -34,42 +34,37 @@ public: tree& operator=(tree const& other) { data_ = other.data_; - set_ranges(other.num_nodes()); + set_ranges(other.num_segments()); return *this; } // copy constructors take advantage of the assignment operators // defined above - tree(tree const& other) - { + tree(tree const& other) { *this = other; } - tree(tree&& other) - { + tree(tree&& other) { *this = std::move(other); } - /// create the tree from a parent_index - template <typename I> - tree(std::vector<I> const& parent_index) - { - // validate the inputs + /// Create the tree from a parent index that lists the parent segment + /// of each segment in a cell tree. + tree(std::vector<int_type> parent_index) { + // validate the input if(!algorithms::is_minimal_degree(parent_index)) { throw std::domain_error( "parent index used to build a tree did not satisfy minimal degree ordering" ); } - auto new_parent_index = algorithms::make_parent_index( - parent_index, algorithms::branches(parent_index)); + // an empty parent_index implies a single-compartment/segment cell + EXPECTS(parent_index.size()!=0u); - init(new_parent_index.size()); - //parents_(memory::all) = new_parent_index; - memory::copy(new_parent_index, parents_); + init(parent_index.size()); + memory::copy(parent_index, parents_); parents_[0] = no_parent; - //child_index_(memory::all) = algorithms::make_index(algorithms::child_count(parents_)); memory::copy(algorithms::make_index(algorithms::child_count(parents_)), child_index_); std::vector<int_type> pos(parents_.size(), 0); @@ -88,8 +83,8 @@ public: return child_index_[b+1] - child_index_[b]; } - size_type num_nodes() const { - // the number of nodes is the size of the child index minus 1 + size_type num_segments() const { + // the number of segments/nodes is the size of the child index minus 1 // ... except for the case of an empty tree auto sz = static_cast<size_type>(child_index_.size()); return sz ? sz - 1 : 0; @@ -129,7 +124,7 @@ public: } iarray change_root(size_t b) { - assert(b<num_nodes()); + assert(b<num_segments()); // no need to rebalance if the root node has been requested if(b==0) { @@ -138,7 +133,7 @@ public: // create new tree with memory allocated tree new_tree; - new_tree.init(num_nodes()); + new_tree.init(num_segments()); // add the root node new_tree.parents_[0] = no_parent; @@ -147,7 +142,7 @@ public: // allocate space for the permutation vector that // will represent the permutation performed on the branches // during the rebalancing - iarray p(num_nodes(), -1); + iarray p(num_segments(), -1); // recersively rebalance the tree new_tree.add_children(0, b, 0, p, *this); @@ -161,7 +156,7 @@ public: // copy in new data with a move because the information in // new_tree is not kept std::swap(data_, new_tree.data_); - set_ranges(new_tree.num_nodes()); + set_ranges(new_tree.num_segments()); return p; } @@ -287,14 +282,14 @@ private: view_type parents_ = data_(0, 0); }; -template <typename IntT, typename SizeT, typename C> -std::vector<IntT> make_parent_index(tree<IntT, SizeT> const& t, C const& counts) +template <typename C> +std::vector<tree::int_type> make_parent_index(tree const& t, C const& counts) { using util::make_span; - using int_type = typename tree<IntT, SizeT>::int_type; - constexpr auto no_parent = tree<IntT, SizeT>::no_parent; + using int_type = tree::int_type; + constexpr auto no_parent = tree::no_parent; - if (!algorithms::is_positive(counts) || counts.size() != t.num_nodes()) { + if (!algorithms::is_positive(counts) || counts.size() != t.num_segments()) { throw std::domain_error( "make_parent_index requires one non-zero count per segment" ); @@ -303,7 +298,7 @@ std::vector<IntT> make_parent_index(tree<IntT, SizeT> const& t, C const& counts) auto num_compartments = index.back(); std::vector<int_type> parent_index(num_compartments); int_type pos = 0; - for (int_type i : make_span(0, t.num_nodes())) { + for (int_type i : make_span(0, t.num_segments())) { // get the parent of this segment // taking care for the case where the root node has -1 as its parent auto parent = t.parent(i); diff --git a/tests/unit/test_algorithms.cpp b/tests/unit/test_algorithms.cpp index ac15a43ca50eedf5fc91116387cdb9c45c54637c..c3fe55baf47ea2a2e8d88ba7c665be550bbd03e9 100644 --- a/tests/unit/test_algorithms.cpp +++ b/tests/unit/test_algorithms.cpp @@ -437,7 +437,7 @@ TEST(algorithms, branches) EXPECT_EQ(expected_branches, actual_branches); auto actual_parent_index = - algorithms::make_parent_index(parent_index, actual_branches); + algorithms::tree_reduce(parent_index, actual_branches); EXPECT_EQ(expected_parent_index, actual_parent_index); // Check find_branch @@ -484,7 +484,7 @@ TEST(algorithms, branches) EXPECT_EQ(expected_branches, actual_branches); auto actual_parent_index = - algorithms::make_parent_index(parent_index, actual_branches); + algorithms::tree_reduce(parent_index, actual_branches); EXPECT_EQ(expected_parent_index, actual_parent_index); } @@ -508,7 +508,7 @@ TEST(algorithms, branches) EXPECT_EQ(expected_branches, actual_branches); auto actual_parent_index = - algorithms::make_parent_index(parent_index, actual_branches); + algorithms::tree_reduce(parent_index, actual_branches); EXPECT_EQ(expected_parent_index, actual_parent_index); } @@ -521,7 +521,7 @@ TEST(algorithms, branches) EXPECT_EQ(expected_branches, actual_branches); auto actual_parent_index = - algorithms::make_parent_index(parent_index, actual_branches); + algorithms::tree_reduce(parent_index, actual_branches); EXPECT_EQ(expected_parent_index, actual_parent_index); } } diff --git a/tests/unit/test_cell.cpp b/tests/unit/test_cell.cpp index 5c197809cbf0f9598aef476e54e993af8e0c18ce..587f8fa95919b1548c3b9101604c5f8f2bfc68a3 100644 --- a/tests/unit/test_cell.cpp +++ b/tests/unit/test_cell.cpp @@ -2,7 +2,7 @@ #include "cell.hpp" -TEST(cell_type, soma) +TEST(cell, soma) { // test that insertion of a soma works // define with no centre point @@ -38,7 +38,7 @@ TEST(cell_type, soma) } } -TEST(cell_type, add_segment) +TEST(cell, add_segment) { using namespace arb; // add a pre-defined segment @@ -101,7 +101,7 @@ TEST(cell_type, add_segment) } } -TEST(cell_type, multiple_cables) +TEST(cell, multiple_cables) { using namespace arb; @@ -148,8 +148,7 @@ TEST(cell_type, multiple_cables) const auto model = c.model(); auto const& con = model.tree; - auto no_parent = cell_tree::no_parent; - + auto no_parent = tree::no_parent; EXPECT_EQ(con.num_segments(), 5u); EXPECT_EQ(con.parent(0), no_parent); EXPECT_EQ(con.parent(1), 0u); @@ -165,7 +164,48 @@ TEST(cell_type, multiple_cables) } } -TEST(cell_type, clone) +TEST(cell, unbranched_chain) +{ + using namespace arb; + + cell c; + + auto soma_radius = std::pow(3./(4.*math::pi<double>()), 1./3.); + + // Cell strucure that looks like a centipede: i.e. each segment has only one child + // + // | | + // 0|1-2-3-4|5-6-7-8 + // | | + + // add a soma + c.add_soma(soma_radius, {0,0,1}); + + // hook the dendrite and axons + c.add_cable(0, make_segment<cable_segment>(section_kind::dendrite, 1.0, 1.0, 1./math::pi<double>())); + c.add_cable(1, make_segment<cable_segment>(section_kind::dendrite, 1.0, 1.0, 1./math::pi<double>())); + + EXPECT_EQ(c.num_segments(), 3u); + // each of the 3 segments has volume 1 by design + EXPECT_EQ(c.volume(), 3.); + // each of the 2 cables has volume 2., and the soma has an awkward area + // that isn't a round number + EXPECT_EQ(c.area(), 4. + math::area_sphere(soma_radius)); + + // construct the graph + const auto tree = c.model().tree; + + auto no_parent = tree::no_parent; + EXPECT_EQ(tree.num_segments(), 3u); + EXPECT_EQ(tree.parent(0), no_parent); + EXPECT_EQ(tree.parent(1), 0u); + EXPECT_EQ(tree.parent(2), 1u); + EXPECT_EQ(tree.num_children(0), 1u); + EXPECT_EQ(tree.num_children(1), 1u); + EXPECT_EQ(tree.num_children(2), 0u); +} + +TEST(cell, clone) { using namespace arb; @@ -220,7 +260,7 @@ TEST(cell_type, clone) EXPECT_EQ(c.segment(2)->num_compartments(), d.segment(2)->num_compartments()); } -TEST(cell_type, get_kind) +TEST(cell, get_kind) { using namespace arb; diff --git a/tests/unit/test_tree.cpp b/tests/unit/test_tree.cpp index a424d9fab50f1d062ced8e733ca3d2f742b5dc43..a6426ce964373cce04bdec6f87453ff51f13213d 100644 --- a/tests/unit/test_tree.cpp +++ b/tests/unit/test_tree.cpp @@ -6,7 +6,7 @@ #include "../gtest.h" -#include <cell_tree.hpp> +#include <tree.hpp> #include <util/debug.hpp> // Path to data directory can be overriden at compile time. @@ -17,42 +17,30 @@ using json = nlohmann::json; using namespace arb; -using int_type = cell_tree::int_type; +using int_type = tree::int_type; -TEST(cell_tree, from_parent_index) { - auto no_parent = cell_tree::no_parent; +TEST(tree, from_segment_index) { + auto no_parent = tree::no_parent; // tree with single branch corresponding to the root node // this is equivalent to a single compartment model // CASE 1 : single root node in parent_index { std::vector<int_type> parent_index = {0}; - cell_tree tree(parent_index); - EXPECT_EQ(tree.num_segments(), 1u); - EXPECT_EQ(tree.num_children(0), 0u); - } - // CASE 2 : empty parent_index - { - std::vector<int_type> parent_index; - cell_tree tree(parent_index); + tree tree(parent_index); EXPECT_EQ(tree.num_segments(), 1u); EXPECT_EQ(tree.num_children(0), 0u); } { // - // 0 0 - // / \ / \. - // 1 4 => 1 2 - // / \. - // 2 5 - // / - // 3 + // 0 + // / \. + // 1 2 // - std::vector<int_type> parent_index = - {0, 0, 1, 2, 0, 4}; - cell_tree tree(parent_index); + std::vector<int_type> parent_index = {0, 0, 0}; + tree tree(parent_index); EXPECT_EQ(tree.num_segments(), 3u); // the root has 2 children EXPECT_EQ(tree.num_children(0), 2u); @@ -62,18 +50,25 @@ TEST(cell_tree, from_parent_index) { } { // - // 0 0 - // /|\ /|\. - // 1 4 6 => 1 2 3 - // / | \. - // 2 5 7 - // / \. - // 3 8 + // 0-1-2-3 // - std::vector<int_type> parent_index = - {0, 0, 1, 2, 0, 4, 0, 6, 7, 8}; - - cell_tree tree(parent_index); + std::vector<int_type> parent_index = {0, 0, 1, 2}; + tree tree(parent_index); + EXPECT_EQ(tree.num_segments(), 4u); + // all non-leaf nodes have 1 child + EXPECT_EQ(tree.num_children(0), 1u); + EXPECT_EQ(tree.num_children(1), 1u); + EXPECT_EQ(tree.num_children(2), 1u); + EXPECT_EQ(tree.num_children(3), 0u); + } + { + // + // 0 + // /|\. + // 1 2 3 + // + std::vector<int_type> parent_index = {0, 0, 0, 0}; + tree tree(parent_index); EXPECT_EQ(tree.num_segments(), 4u); // the root has 3 children EXPECT_EQ(tree.num_children(0), 3u); @@ -90,23 +85,14 @@ TEST(cell_tree, from_parent_index) { } { // - // 0 0 - // /|\ /|\. - // 1 4 6 => 1 2 3 - // / | \ / \. - // 2 5 7 4 5 - // / \. - // 3 8 - // / \. - // 9 11 - // / \. - // 10 12 - // \. - // 13 + // 0 + // /|\. + // 1 2 3 + // / \. + // 4 5 // - std::vector<int_type> parent_index = - {0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12}; - cell_tree tree(parent_index); + std::vector<int_type> parent_index = {0, 0, 0, 0, 3, 3}; + tree tree(parent_index); EXPECT_EQ(tree.num_segments(), 6u); // the root has 3 children EXPECT_EQ(tree.num_children(0), 3u); @@ -134,7 +120,7 @@ TEST(cell_tree, from_parent_index) { // / \. // 2 3 std::vector<int_type> parent_index = {0,0,1,1}; - cell_tree tree(parent_index); + tree tree(parent_index); EXPECT_EQ(tree.num_segments(), 4u); @@ -151,7 +137,7 @@ TEST(cell_tree, from_parent_index) { // / \. // 2 3 std::vector<int_type> parent_index = {0,0,1,1,0,0}; - cell_tree tree(parent_index); + tree tree(parent_index); EXPECT_EQ(tree.num_segments(), 6u); @@ -177,7 +163,7 @@ TEST(cell_tree, from_parent_index) { // / \. // 5 6 std::vector<int_type> parent_index = {0,0,0,1,1,4,4}; - cell_tree tree(parent_index); + tree tree(parent_index); EXPECT_EQ(tree.num_segments(), 7u); @@ -191,242 +177,3 @@ TEST(cell_tree, from_parent_index) { } } -TEST(cell_tree, depth_from_root) { - { - // 0 - // / \. - // 1 2 - // / \. - // 3 4 - // / \. - // 5 6 - std::vector<int_type> parent_index = {0,0,0,1,1,4,4}; - cell_tree tree(parent_index); - auto depth = tree.depth_from_root(); - - EXPECT_EQ(depth[0], 0u); - EXPECT_EQ(depth[1], 1u); - EXPECT_EQ(depth[2], 1u); - EXPECT_EQ(depth[3], 2u); - EXPECT_EQ(depth[4], 2u); - EXPECT_EQ(depth[5], 3u); - EXPECT_EQ(depth[6], 3u); - } - { - // 0 - // / \. - // 1 2 - // / \. - // 3 4 - // / \. - // 5 6 - std::vector<int_type> parent_index = {0,0,0,2,2,4,4}; - cell_tree tree(parent_index); - auto depth = tree.depth_from_root(); - - EXPECT_EQ(depth[0], 0u); - EXPECT_EQ(depth[1], 1u); - EXPECT_EQ(depth[2], 1u); - EXPECT_EQ(depth[3], 2u); - EXPECT_EQ(depth[4], 2u); - EXPECT_EQ(depth[5], 3u); - EXPECT_EQ(depth[6], 3u); - } -} - -TEST(tree, change_root) { - { - // a cell with the following structure - // make 1 the new root - // 0 0 - // / \ | - // 1 2 -> 1 - // | - // 2 - std::vector<int_type> parent_index = {0,0,0}; - tree<int_type> t(parent_index); - t.change_root(1); - - EXPECT_EQ(t.num_nodes(), 3u); - - EXPECT_EQ(t.num_children(0), 1u); - EXPECT_EQ(t.num_children(1), 1u); - EXPECT_EQ(t.num_children(2), 0u); - } - { - // a cell with the following structure - // make 1 the new root - // 0 0 - // / \ /|\. - // 1 2 -> 1 2 3 - // / \ | - // 3 4 4 - std::vector<int_type> parent_index = {0,0,0,1,1}; - tree<int_type> t(parent_index); - t.change_root(1u); - - EXPECT_EQ(t.num_nodes(), 5u); - - EXPECT_EQ(t.num_children(0), 3u); - EXPECT_EQ(t.num_children(1), 0u); - EXPECT_EQ(t.num_children(2), 0u); - EXPECT_EQ(t.num_children(3), 1u); - EXPECT_EQ(t.num_children(4), 0u); - } - { - // a cell with the following structure - // make 1 the new root - // unlike earlier tests, this decreases the depth - // of the tree - // 0 0 - // / \ /|\. - // 1 2 -> 1 2 5 - // / \ / \ \. - // 3 4 3 4 6 - // / \. - // 5 6 - std::vector<int_type> parent_index = {0,0,0,1,1,4,4}; - tree<int_type> t(parent_index); - - t.change_root(1); - - EXPECT_EQ(t.num_nodes(), 7u); - - EXPECT_EQ(t.num_children(0), 3u); - EXPECT_EQ(t.num_children(1), 0u); - EXPECT_EQ(t.num_children(2), 2u); - EXPECT_EQ(t.num_children(3), 0u); - EXPECT_EQ(t.num_children(4), 0u); - EXPECT_EQ(t.num_children(5), 1u); - EXPECT_EQ(t.num_children(6), 0u); - } -} - -TEST(cell_tree, balance) { - auto no_parent = cell_tree::no_parent; - - { - // a cell with the following structure - // will balance around 1 - // 0 0 - // / \ /|\. - // 1 2 -> 1 2 5 - // / \ / \ \. - // 3 4 3 4 6 - // / \. - // 5 6 - std::vector<int_type> parent_index = {0,0,0,1,1,4,4}; - cell_tree t(parent_index); - - t.balance(); - - // the soma (original root) has moved to 5 in the new tree - EXPECT_EQ(t.soma(), 5u); - - EXPECT_EQ(t.num_segments(), 7u); - EXPECT_EQ(t.num_children(0),3u); - EXPECT_EQ(t.num_children(1),0u); - EXPECT_EQ(t.num_children(2),2u); - EXPECT_EQ(t.num_children(3),0u); - EXPECT_EQ(t.num_children(4),0u); - EXPECT_EQ(t.num_children(5),1u); - EXPECT_EQ(t.num_children(6),0u); - EXPECT_EQ(t.parent(0), no_parent); - EXPECT_EQ(t.parent(1), 0u); - EXPECT_EQ(t.parent(2), 0u); - EXPECT_EQ(t.parent(3), 0u); - EXPECT_EQ(t.parent(4), 2u); - EXPECT_EQ(t.parent(5), 2u); - EXPECT_EQ(t.parent(6), 5u); - - //t.to_graphviz("cell.dot"); - } -} - -// this test doesn't test anything yet... it just loads each cell in turn -// from a json file and creates a .dot file for it -TEST(cell_tree, json_load) -{ - json cell_data; - std::string path{DATADIR}; - - path += "/cells_small.json"; - std::ifstream(path) >> cell_data; - - for(auto c : util::make_span(0,cell_data.size())) { - std::vector<int_type> parent_index = cell_data[c]["parent_index"]; - cell_tree tree(parent_index); - //tree.to_graphviz("cell" + std::to_string(c) + ".dot"); - } -} - -TEST(tree, make_parent_index) -{ - /* - // just the soma - { - std::vector<int> parent_index = {0}; - std::vector<int> counts = {1}; - arb::tree t(parent_index); - auto new_parent_index = make_parent_index(t, counts); - EXPECT_EQ(parent_index.size(), new_parent_index.size()); - } - // just a soma with 5 compartments - { - std::vector<int> parent_index = {0}; - std::vector<int> counts = {5}; - arb::tree t(parent_index); - auto new_parent_index = make_parent_index(t, counts); - EXPECT_EQ(new_parent_index.size(), (unsigned)counts[0]); - EXPECT_EQ(new_parent_index[0], 0); - for(auto i=1u; i<new_parent_index.size(); ++i) { - EXPECT_EQ((unsigned)new_parent_index[i], i-1); - } - } - // some trees with single compartment per segment - { - auto trees = { - // 0 - // | - // 1 - std::vector<int>{0,0}, - // 0 - // / \. - // 1 2 - std::vector<int>{0,0,0}, - // 0 - // / \. - // 1 4 - // / \ |\. - // 2 3 5 6 - std::vector<int>{0,0,0,1,1,2,2} - }; - for(auto &parent_index : trees) { - std::vector<int> counts(parent_index.size(), 1); - arb::tree t(parent_index); - auto new_parent_index = make_parent_index(t, counts); - EXPECT_EQ(parent_index, new_parent_index); - } - } - // a tree with multiple compartments per segment - // - // 0 - // / \. - // 1 8 - // / \. - // 2 9 - // /. - // 3 - // / \. - // 4 6 - // / \. - // 5 7 - { - std::vector<int> parent_index = {0,0,1,2,3,4,3,6,0,8}; - std::vector<int> counts = {1,3,2,2,2}; - arb::tree t(parent_index); - auto new_parent_index = make_parent_index(t, counts); - EXPECT_EQ(parent_index, new_parent_index); - } - */ -}