diff --git a/.ycm_extra_conf.py b/.ycm_extra_conf.py
index 6f707b356d9c22203c04167a980c93c24056c4a5..ea97fe481c903d1e5319bf4a9e0174069c75ddc5 100644
--- a/.ycm_extra_conf.py
+++ b/.ycm_extra_conf.py
@@ -41,16 +41,16 @@ flags = [
'c++',
'-I',
'.',
- '-I',
- '/usr/include/c++/4.9.2',
'-isystem',
- '/usr/lib/gcc/x86_64-unknown-linux-gnu/4.9.2/include'
+ '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1'
'-isystem',
- '/usr/local/include',
- '-I',
- '/home/bcumming/software/github/modparser/tests',
- '-I',
- '/home/bcumming/software/msolve/json',
+ '/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.10.sdk/usr/include/c++/4.2.1/bits'
+# '-I',
+# '/usr/include/c++/4.9.2',
+# '-isystem',
+# '/usr/lib/gcc/x86_64-unknown-linux-gnu/4.9.2/include'
+# '-isystem',
+# '/usr/local/include',
]
diff --git a/cell_tree.hpp b/cell_tree.hpp
index 9073dd18b35ecceab1e1fb8922dc29a8d708c0f9..bcacb25c171f31439f0c4711649fccfd36ab8642 100644
--- a/cell_tree.hpp
+++ b/cell_tree.hpp
@@ -1,34 +1,34 @@
#pragma once
#include <algorithm>
+#include <iomanip>
#include <iostream>
#include <fstream>
-#include <iomanip>
+#include <memory>
#include <numeric>
#include <ostream>
#include <vector>
-#include <memory>
#include "vector/include/Vector.hpp"
#include "tree.hpp"
#include "util.hpp"
-// The tree data structure that describes the branches 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 branches,
-// 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.
+/// 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 :
// use a signed 16-bit integer for storage of indexes, which is reasonable given
- // that typical cells have at most 1000-2000 nodes
+ // that typical cells have at most 1000-2000 segments
using int_type = int16_t;
using index_type = memory::HostVector<int_type>;
using index_view = index_type::view_type;
@@ -37,140 +37,70 @@ class cell_tree {
cell_tree(std::vector<int> const& parent_index)
{
// handle the case of an empty parent list, which implies a single-compartment model
- std::vector<int> branch_index;
+ std::vector<int> segment_index;
if(parent_index.size()>0) {
- branch_index = tree_.init_from_parent_index(parent_index);
+ segment_index = tree_.init_from_parent_index(parent_index);
}
else {
- branch_index = tree_.init_from_parent_index(std::vector<int>({0}));
+ segment_index = tree_.init_from_parent_index(std::vector<int>({0}));
}
// if needed, calculate meta-data like length[] and end[] arrays for data
}
- /// Minimize the depth of the tree.
- /// Pick a root node that minimizes the depth of the tree.
- /*
- int_type balance() {
- if(num_branches()==1) {
- return 0;
- }
+ /// construct from a tree
+ cell_tree(tree const& t)
+ : tree_(t)
+ {}
- // calculate the depth of each branch from the root
- // pre-order traversal of the tree
- index_type depth(num_branches());
- auto depth_from_root = [this, &depth] (int_type b) -> void
- {
- auto d = depth[tree.parent(b)] + 1;
- depth[b] = d;
- for(auto c : tree.children(b)) {
- depth_from_root(c);
- }
- }
- depth[0]=0;
- depth_from_root(0);
-
- auto max = std::max_element(depth.begin(), depth.end());
- auto max_leaf = std::distance(depth.begin(), max);
- auto original_depth = *max;
+ /// construct from a tree
+ cell_tree(tree&& t)
+ : tree_(std::move(t))
+ {}
- // Calculate the depth of each compartment as the maximum distance
- // from a child leaf
- // post-order traversal of the tree
- auto depth_from_leaf = [this, &depth] (int_type b)
- {
- int_type max_depth = 0;
- for(auto c : children(branch)) {
- max_depth = std::max(max_depth, depth_from_leaf(c));
- }
- depth[b] = max_depth;
- return max_depth+1;
- }
- depth_from_leaf(0);
-
- // 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);
- auto distance_from_max_leaf = 1;
- auto parent = max_leaf;
- auto pos = parent(max_leaf);
- while(pos != -1) {
- for(auto c : children(pos)) {
- if(c!=parent) {
- 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);
- }
- }
- }
- parent = pos;
- pos = parent(pos);
- ++distance_from_max_leaf;
- }
-
- std::cout << memory::util::green(std::to_string(max_sub_tree.root)) << " ";
-
- // 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) {
- std::cout << " root " << 0 << " depth " << original_depth << std::endl;
- return *max;
- }
-
- // 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;
- }
+ /// 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
- tree.change_root(new_root);
+ tree_.change_root(new_root);
- return new_depth;
+ 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 branches in the cell
- size_t num_branches() const {
+ /// returns the number of segments in the cell
+ size_t num_segments() const {
return tree_.num_nodes();
}
- /// returns the number of child branches of branch b
+ /// returns the number of child segments of segment b
size_t num_children(size_t b) const {
return tree_.num_children(b);
}
- /// returns a list of the children of branch b
+ /// returns a list of the children of segment b
const index_view children(size_t b) const {
return tree_.children(b);
}
- /// returns the parent of branch b
+ /// returns the parent index of segment b
int_type parent(size_t b) const {
return tree_.parent(b);
}
- /// generates a graphviz .dot file that visualizes cell branch structure
+ /// 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 {" << std::endl;
- for(auto b : range(0,num_branches())) {
+ for(auto b : range(0,num_segments())) {
if(children(b).size()) {
for(auto c : children(b)) {
fid << " " << b << " -- " << c << ";" << std::endl;
@@ -180,6 +110,21 @@ class cell_tree {
fid << "}" << std::endl;
}
+ index_type depth_from_leaf()
+ {
+ tree::index_type depth(num_segments());
+ depth_from_leaf(depth, 0);
+ return depth;
+ }
+
+ index_type depth_from_root()
+ {
+ tree::index_type depth(num_segments());
+ depth[0] = 0;
+ depth_from_root(depth, 1);
+ return depth;
+ }
+
private :
@@ -211,7 +156,89 @@ class cell_tree {
int depth;
};
- // storage for the tree structure of cell branches
+ /// 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);
+ auto original_depth = *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);
+ auto distance_from_max_leaf = 1;
+ auto pnt = max_leaf;
+ auto pos = parent(max_leaf);
+ while(pos != -1) {
+ 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(index_type& depth, int 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(index_type& depth, int segment)
+ {
+ auto d = depth[parent(segment)] + 1;
+ depth[segment] = d;
+ for(auto c : children(segment)) {
+ depth_from_root(depth, c);
+ }
+ }
+
+ // storage for the tree structure of cell segments
tree tree_;
};
-
diff --git a/main.cpp b/main.cpp
index cd7a910ee3ac2c8c8b1cd7053ad10af1be46de6a..e15fb88707d58ee6b0b72670e7559bcf0d5eb360 100644
--- a/main.cpp
+++ b/main.cpp
@@ -17,22 +17,22 @@ TEST(cell_tree, from_parent_index) {
{
std::vector<int> parent_index = {0};
cell_tree tree(parent_index);
- EXPECT_EQ(tree.num_branches(), 1);
+ EXPECT_EQ(tree.num_segments(), 1);
EXPECT_EQ(tree.num_children(0), 0);
}
// CASE 2 : empty parent_index
{
std::vector<int> parent_index;
cell_tree tree(parent_index);
- EXPECT_EQ(tree.num_branches(), 1);
+ EXPECT_EQ(tree.num_segments(), 1);
EXPECT_EQ(tree.num_children(0), 0);
}
- // tree with two branches off the root node
+ // tree with two segments off the root node
{
std::vector<int> parent_index =
{0, 0, 1, 2, 0, 4};
cell_tree tree(parent_index);
- EXPECT_EQ(tree.num_branches(), 3);
+ EXPECT_EQ(tree.num_segments(), 3);
// the root has 2 children
EXPECT_EQ(tree.num_children(0), 2);
// the children are leaves
@@ -40,11 +40,11 @@ TEST(cell_tree, from_parent_index) {
EXPECT_EQ(tree.num_children(2), 0);
}
{
- // tree with three branches off the root node
+ // tree with three segments off the root node
std::vector<int> parent_index =
{0, 0, 1, 2, 0, 4, 0, 6, 7, 8};
cell_tree tree(parent_index);
- EXPECT_EQ(tree.num_branches(), 4);
+ EXPECT_EQ(tree.num_segments(), 4);
// the root has 3 children
EXPECT_EQ(tree.num_children(0), 3);
// the children are leaves
@@ -53,11 +53,11 @@ TEST(cell_tree, from_parent_index) {
EXPECT_EQ(tree.num_children(3), 0);
}
{
- // tree with three branches off the root node, and another 2 branches off of the third branch from the root node
+ // tree with three segments off the root node, and another 2 segments off of the third branch from the root node
std::vector<int> parent_index =
{0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12};
cell_tree tree(parent_index);
- EXPECT_EQ(tree.num_branches(), 6);
+ EXPECT_EQ(tree.num_segments(), 6);
// the root has 3 children
EXPECT_EQ(tree.num_children(0), 3);
// one of the chilren has 2 children ...
@@ -78,7 +78,7 @@ TEST(cell_tree, from_parent_index) {
std::vector<int> parent_index = {0,0,1,1};
cell_tree tree(parent_index);
- EXPECT_EQ(tree.num_branches(), 4);
+ EXPECT_EQ(tree.num_segments(), 4);
EXPECT_EQ(tree.num_children(0), 1);
EXPECT_EQ(tree.num_children(1), 2);
@@ -95,7 +95,7 @@ TEST(cell_tree, from_parent_index) {
std::vector<int> parent_index = {0,0,0,1,1};
cell_tree tree(parent_index);
- EXPECT_EQ(tree.num_branches(), 5);
+ EXPECT_EQ(tree.num_segments(), 5);
EXPECT_EQ(tree.num_children(0), 2);
EXPECT_EQ(tree.num_children(1), 2);
@@ -114,7 +114,7 @@ TEST(cell_tree, from_parent_index) {
std::vector<int> parent_index = {0,0,0,1,1,4,4};
cell_tree tree(parent_index);
- EXPECT_EQ(tree.num_branches(), 7);
+ EXPECT_EQ(tree.num_segments(), 7);
EXPECT_EQ(tree.num_children(0), 2);
EXPECT_EQ(tree.num_children(1), 2);
@@ -129,12 +129,54 @@ TEST(cell_tree, from_parent_index) {
TEST(tree, change_root) {
{
// a cell with the following structure
- // should be rebalanced around node 1
- // 0
- // / \
- // 1 2
- // / \
- // 3 4
+ // make 1 the new root
+ // 0 0
+ // / \ |
+ // 1 2 -> 1
+ // |
+ // 2
+ std::vector<int> parent_index = {0,0,0};
+ tree t;
+ t.init_from_parent_index(parent_index);
+ auto new_tree = t.change_root(1);
+
+ EXPECT_EQ(new_tree.num_nodes(), 3);
+
+ EXPECT_EQ(new_tree.num_children(0), 1);
+ EXPECT_EQ(new_tree.num_children(1), 1);
+ EXPECT_EQ(new_tree.num_children(2), 0);
+ }
+ {
+ // a cell with the following structure
+ // make 1 the new root
+ // 0 0
+ // / \ /|\
+ // 1 2 -> 1 2 3
+ // / \ |
+ // 3 4 4
+ std::vector<int> parent_index = {0,0,0,1,1};
+ tree t;
+ t.init_from_parent_index(parent_index);
+ auto new_tree = t.change_root(1);
+
+ EXPECT_EQ(new_tree.num_nodes(), 5);
+
+ EXPECT_EQ(new_tree.num_children(0), 3);
+ EXPECT_EQ(new_tree.num_children(1), 0);
+ EXPECT_EQ(new_tree.num_children(2), 0);
+ EXPECT_EQ(new_tree.num_children(3), 1);
+ EXPECT_EQ(new_tree.num_children(4), 0);
+ }
+ {
+ // 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> parent_index = {0,0,0,1,1,4,4};
@@ -142,6 +184,40 @@ TEST(tree, change_root) {
t.init_from_parent_index(parent_index);
auto new_tree = t.change_root(1);
+
+ EXPECT_EQ(new_tree.num_nodes(), 7);
+
+ EXPECT_EQ(new_tree.num_children(0), 3);
+ EXPECT_EQ(new_tree.num_children(1), 0);
+ EXPECT_EQ(new_tree.num_children(2), 2);
+ EXPECT_EQ(new_tree.num_children(3), 0);
+ EXPECT_EQ(new_tree.num_children(4), 0);
+ EXPECT_EQ(new_tree.num_children(5), 1);
+ EXPECT_EQ(new_tree.num_children(6), 0);
+
+ auto T = cell_tree(new_tree);
+ }
+}
+
+TEST(cell_tree, balance) {
+ {
+ // 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> parent_index = {0,0,0,1,1,4,4};
+ cell_tree t(parent_index);
+
+ t.balance();
+
+ EXPECT_EQ(t.num_segments(), 7);
+
+ t.to_graphviz("cell.dot");
}
}
diff --git a/tree.hpp b/tree.hpp
index a887234a4cd76bbbfab377ed8fcf938792b19a15..6edc61e77ae2d12b60d77b0d61ac2499fdd40f1d 100644
--- a/tree.hpp
+++ b/tree.hpp
@@ -1,8 +1,11 @@
#pragma once
#include <algorithm>
+#include <numeric>
#include <vector>
+#include <cassert>
+
#include "vector/include/Vector.hpp"
#include "util.hpp"
@@ -23,17 +26,14 @@ class tree {
}
tree(tree&& other)
+ : data_(std::move(other.data_))
{
- data_ = std::move(other.data_);
set_ranges(other.num_nodes());
}
tree() = default;
/// create the tree from a parent_index
- /// maybe this should be an initializer, not a constructor, because
- /// we currently return the branch index in place of the parent_index
- /// which feels like bad design.
template <typename I>
std::vector<I>
init_from_parent_index(std::vector<I> const& parent_index)
@@ -187,18 +187,28 @@ class tree {
tree new_tree;
new_tree.init(num_nodes());
+ // mark all nodes
+ new_tree.parents_(memory::all) = -1;
+ new_tree.child_index_(memory::all) = -1;
+ new_tree.children_(memory::all) = -1;
+
// add the root node
- //new_tree.parents_[0] = -1;
+ new_tree.parents_[0] = -1;
new_tree.child_index_[0] = 0;
// allocate space for the permutation vector that
// will represent the permutation performed on the branches
// during the rebalancing
- index_type p(num_nodes()+1);
- std::cout << "allocated for array of size " << p.size() << std::endl;
+ index_type p(num_nodes(), -1);
// recersively rebalance the tree
- new_tree.add_children(0, b, p, *this, true);
+ new_tree.add_children(0, b, 0, p, *this);
+
+ // renumber the child indexes
+ std::transform(
+ new_tree.children_.begin(), new_tree.children_.end(),
+ new_tree.children_.begin(), [&p] (int i) {return p[i];}
+ );
return new_tree;
}
@@ -236,58 +246,67 @@ class tree {
/// will be applied recursively until the old root has been processed,
/// which indicates that the renumbering is finished.
///
- /// precondition - the node new node has already been placed in the tree
+ /// precondition - the node new_node has already been placed in the tree
+ /// precondition - all of new_node's children have been added to the tree
+ /// new_node : the new index of the node whose children are to be added
+ /// to the tree
+ /// old_node : the index of new_node in the original tree
+ /// parent_node : equals index of old_node's parent in the original tree
+ /// should be a child of new_node
+ /// : equals -1 if the old_node's parent is not a child of
+ /// new_node
+ /// p : permutation vector, p[i] is the new index of node i in the old
+ /// tree
int add_children(
int new_node,
int old_node,
+ int parent_node,
index_view p,
- tree const& old_tree,
- bool parent_as_child
+ tree const& old_tree
)
{
- std::cout << "add_children(old " << old_node << ", new " << new_node << ", " << parent_as_child << ")" << std::endl;
-
// check for the senitel that indicates that the old root has
// been processed
if(old_node==-1) {
- assert(parent_as_child); // sanity check
+ //assert(parent_node<0); // sanity check
return new_node;
}
+ p[old_node] = new_node;
+
+ // the list of the children of the original node
auto old_children = old_tree.children(old_node);
- auto pos = child_index_[new_node];
- auto n_children = old_children.size() + (parent_as_child ? 1 : 0);
- child_index_[new_node+1] = pos + n_children;
- std::cout << " inserting " << n_children << " into " << child_index_(new_node, new_node+2) << std::endl;
+ auto this_node = new_node;
+ auto pos = child_index_[this_node];
- std::cout << " inserting";
- // first renumber the children
+ auto add_parent_as_child = parent_node>=0 && old_node>0;
+ //
+ // STEP 1 : add the child indexes for this_node
+ //
+ // first add the children of the node
for(auto b : old_children) {
- std::cout << " " << b;
- new_node++;
- p[new_node] = b;
- children_[pos++] = new_node;
+ if(b != parent_node) {
+ children_[pos++] = b;
+ }
}
- // then add and renumber the parent as a child
- if(parent_as_child) {
- std::cout << " " << yellow(std::to_string(old_tree.parent(old_node)));
- new_node++;
- p[new_node] = old_tree.parent(old_node);
- children_[pos++] = new_node;
+ // then add the node's parent as a child if applicable
+ if(add_parent_as_child) {
+ children_[pos++] = old_tree.parent(old_node);
}
- std::cout << std::endl;
+ child_index_[this_node+1] = pos;
- // then visit the sub-tree of each child recursively
- // - traverse _down_ the tree
- //auto child_range = range();
+ //
+ // STEP 2 : recursively add each child's children
+ //
+ new_node++;
for(auto b : old_children) {
- new_node = add_children(new_node, b, p, old_tree, false);
+ if(b != parent_node) {
+ new_node = add_children(new_node, b, -1, p, old_tree);
+ }
}
- // finally visit the parent recursively
- // - traverse _up_ the tree towards the old root
- if(parent_as_child) {
- new_node = add_children(new_node, parents_[old_node], p, old_tree, true);
+ if(add_parent_as_child) {
+ new_node = add_children(new_node, old_tree.parent(old_node), old_node, p, old_tree);
}
return new_node;
diff --git a/util.hpp b/util.hpp
index efc28a7d7b150d884ba69aa3340aa72e4ac8b9ed..d37c7f65427f1bd7e12edc8ee10b0f95fcb47ba8 100644
--- a/util.hpp
+++ b/util.hpp
@@ -9,8 +9,8 @@ using memory::util::white;
using memory::util::blue;
using memory::util::cyan;
-#include <vector>
#include <ostream>
+#include <vector>
template <typename T>
std::ostream&