Skip to content
Snippets Groups Projects
Commit 338be38f authored by Benjamin Cumming's avatar Benjamin Cumming
Browse files

balancing of cell trees now works

parent 42431d4c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
.ycm_extra_conf.py
+ 8
− 8
View file @ 338be38f
......@@ -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',
]
......
cell_tree.hpp
+ 137
− 110
View file @ 338be38f
#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_;
};
main.cpp
+ 93
− 17
View file @ 338be38f
......@@ -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");
}
}
......
tree.hpp
+ 58
− 39
View file @ 338be38f
#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;
......
util.hpp
+ 1
− 1
View file @ 338be38f
......@@ -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&
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment