diff --git a/CMakeLists.txt b/CMakeLists.txt
index 72ac867917e76a6b92ffe0a2073d0fd6d4f72b0a..f91e332eecdd77b327846951ed9271d740e926df 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -20,27 +20,29 @@ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib)
# TBB support
-set( WITH_TBB OFF CACHE BOOL "use TBB for on-node threading" )
-if(${WITH_TBB})
+set(WITH_TBB OFF CACHE BOOL "use TBB for on-node threading" )
+if(WITH_TBB)
find_package(TBB REQUIRED)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DWITH_TBB ${TBB_DEFINITIONS}")
- link_directories(${TBB_LIBRARY})
endif()
# MPI support
-# relies on the user specifying an MPI-aware compiler wrapper
-set( WITH_MPI OFF CACHE BOOL "use MPI for distrubuted parallelism" )
-if(${WITH_MPI})
- set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DWITH_MPI")
+set(WITH_MPI OFF CACHE BOOL "use MPI for distrubuted parallelism")
+if(WITH_MPI)
+ find_package(MPI REQUIRED)
+ include_directories(SYSTEM ${MPI_C_INCLUDE_PATH})
+ add_definitions(-DWITH_MPI)
+ # unfortunate workaround for C++ detection in system mpi.h
+ add_definitions(-DMPICH_SKIP_MPICXX=1 -DOMPI_SKIP_MPICXX=1)
+ set_property(DIRECTORY APPEND_STRING PROPERTY COMPILE_OPTIONS "${MPI_C_COMPILE_FLAGS}")
endif()
# Cray systems
-set( SYSTEM_CRAY OFF CACHE BOOL "add flags for compilation on Cray systems" )
-if( ${SYSTEM_CRAY} )
+set(SYSTEM_CRAY OFF CACHE BOOL "add flags for compilation on Cray systems")
+if(SYSTEM_CRAY)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -dynamic")
endif()
-
# targets for extermal dependencies
include(ExternalProject)
externalproject_add(modparser
diff --git a/external/vector b/external/vector
index 193df3683a3eaa2bf203bbcfa935d7e7019b1719..8284611f05b0fbe21a1f84630e2726015cb1d96d 160000
--- a/external/vector
+++ b/external/vector
@@ -1 +1 @@
-Subproject commit 193df3683a3eaa2bf203bbcfa935d7e7019b1719
+Subproject commit 8284611f05b0fbe21a1f84630e2726015cb1d96d
diff --git a/miniapp/CMakeLists.txt b/miniapp/CMakeLists.txt
index 883905caa9a5ffe2ebd6c244ac0565e620ed6a81..a5b9561198875a6515b4e8a69f7685ea75c19803 100644
--- a/miniapp/CMakeLists.txt
+++ b/miniapp/CMakeLists.txt
@@ -11,6 +11,11 @@ add_executable(miniapp.exe ${MINIAPP_SOURCES} ${HEADERS})
target_link_libraries(miniapp.exe LINK_PUBLIC cellalgo)
target_link_libraries(miniapp.exe LINK_PUBLIC ${TBB_LIBRARIES})
+if(WITH_MPI)
+ target_link_libraries(miniapp.exe LINK_PUBLIC ${MPI_C_LIBRARIES})
+ set_property(TARGET miniapp.exe APPEND_STRING PROPERTY LINK_FLAGS "${MPI_C_LINK_FLAGS}")
+endif()
+
set_target_properties(miniapp.exe
PROPERTIES
RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/miniapp"
diff --git a/src/algorithms.hpp b/src/algorithms.hpp
index 83eeaf8dd7e32f776384d0c47fec5e5d0b43d1a8..738f163fa8381363b099c981d57df36923bc921d 100644
--- a/src/algorithms.hpp
+++ b/src/algorithms.hpp
@@ -19,230 +19,302 @@
namespace nest {
namespace mc {
-namespace algorithms{
-
- template <typename C>
- typename C::value_type
- sum(C const& c)
- {
- using value_type = typename C::value_type;
- return std::accumulate(c.begin(), c.end(), value_type{0});
- }
-
- template <typename C>
- typename C::value_type
- mean(C const& c)
- {
- return sum(c)/c.size();
+namespace algorithms {
+
+template <typename C>
+typename C::value_type
+sum(C const& c)
+{
+ using value_type = typename C::value_type;
+ return std::accumulate(c.begin(), c.end(), value_type{0});
+}
+
+template <typename C>
+typename C::value_type
+mean(C const& c)
+{
+ return sum(c)/c.size();
+}
+
+template <typename C>
+C make_index(C const& c)
+{
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "make_index only applies to integral types"
+ );
+
+ C out(c.size()+1);
+ out[0] = 0;
+ std::partial_sum(c.begin(), c.end(), out.begin()+1);
+ return out;
+}
+
+/// works like std::is_sorted(), but with stronger condition that succesive
+/// elements must be greater than those before them
+template <typename C>
+bool is_strictly_monotonic_increasing(C const& c)
+{
+ using value_type = typename C::value_type;
+ return std::is_sorted(
+ c.begin(),
+ c.end(),
+ [] (value_type const& lhs, value_type const& rhs) {
+ return lhs <= rhs;
+ }
+ );
+}
+
+template <typename C>
+bool is_strictly_monotonic_decreasing(C const& c)
+{
+ using value_type = typename C::value_type;
+ return std::is_sorted(
+ c.begin(),
+ c.end(),
+ [] (value_type const& lhs, value_type const& rhs) {
+ return lhs >= rhs;
+ }
+ );
+}
+
+template <
+ typename C,
+ typename = typename std::enable_if<std::is_integral<typename C::value_type>::value>
+>
+bool is_minimal_degree(C const& c)
+{
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "is_minimal_degree only applies to integral types"
+ );
+
+ if(c.size()==0u) {
+ return true;
}
- template <typename C>
- C make_index(C const& c)
- {
- static_assert(
- std::is_integral<typename C::value_type>::value,
- "make_index only applies to integral types"
- );
-
- C out(c.size()+1);
- out[0] = 0;
- std::partial_sum(c.begin(), c.end(), out.begin()+1);
- return out;
+ using value_type = typename C::value_type;
+ if(c[0] != value_type(0)) {
+ return false;
}
-
- /// works like std::is_sorted(), but with stronger condition that succesive
- /// elements must be greater than those before them
- template <typename C>
- bool is_strictly_monotonic_increasing(C const& c)
- {
- using value_type = typename C::value_type;
- return std::is_sorted(
- c.begin(),
- c.end(),
- [] (value_type const& lhs, value_type const& rhs) {
- return lhs <= rhs;
- }
- );
+ auto i = value_type(1);
+ auto it = std::find_if(
+ c.begin()+1, c.end(), [&i](value_type v) { return v>=(i++); }
+ );
+ return it==c.end();
+}
+
+template <typename C>
+bool is_positive(C const& c)
+{
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "is_positive only applies to integral types"
+ );
+ for(auto v : c) {
+ if(v<1) {
+ return false;
+ }
}
-
- template <typename C>
- bool is_strictly_monotonic_decreasing(C const& c)
- {
- using value_type = typename C::value_type;
- return std::is_sorted(
- c.begin(),
- c.end(),
- [] (value_type const& lhs, value_type const& rhs) {
- return lhs >= rhs;
- }
- );
+ return true;
+}
+
+template<typename C>
+bool has_contiguous_segments(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;
}
- template <
- typename C,
- typename = typename std::enable_if<std::is_integral<typename C::value_type>::value>
- >
- bool is_minimal_degree(C const& c)
- {
- static_assert(
- std::is_integral<typename C::value_type>::value,
- "is_minimal_degree only applies to integral types"
- );
-
- if(c.size()==0u) {
- return true;
- }
+ int n = parent_index.size();
+ std::vector<bool> is_leaf(n, false);
- using value_type = typename C::value_type;
- if(c[0] != value_type(0)) {
+ for(auto i=1; i<n; ++i) {
+ auto p = parent_index[i];
+ if(is_leaf[p]) {
return false;
}
- auto i = value_type(1);
- auto it = std::find_if(
- c.begin()+1, c.end(), [&i](value_type v) { return v>=(i++); }
- );
- return it==c.end();
- }
- template <typename C>
- bool is_positive(C const& c)
- {
- static_assert(
- std::is_integral<typename C::value_type>::value,
- "is_positive only applies to integral types"
- );
- for(auto v : c) {
- if(v<1) {
- return false;
- }
+ if(p != i-1) {
+ // we have a branch and i-1 is a leaf node
+ is_leaf[i-1] = true;
}
- return true;
}
- template<typename C>
- bool has_contiguous_segments(const C &parent_index)
- {
- static_assert(
- std::is_integral<typename C::value_type>::value,
- "integral type required"
- );
+ return true;
+}
- if (!is_minimal_degree(parent_index)) {
- return false;
- }
+template<typename C>
+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"
+ );
- int n = parent_index.size();
- std::vector<bool> is_leaf(n, false);
+ 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]];
+ }
- for(auto i=1; i<n; ++i) {
- auto p = parent_index[i];
- if(is_leaf[p]) {
- return false;
- }
+ return count;
+}
- if(p != i-1) {
- // we have a branch and i-1 is a leaf node
- is_leaf[i-1] = true;
- }
- }
+template<typename C>
+std::vector<typename C::value_type> branches(const C& parent_index)
+{
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "integral type required"
+ );
- return true;
+ EXPECTS(has_contiguous_segments(parent_index));
+
+ std::vector<typename C::value_type> branch_index;
+ if (parent_index.empty()) {
+ return branch_index;
}
- template<typename C>
- 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"
- );
-
- 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]];
+ auto num_child = child_count(parent_index);
+ branch_index.push_back(0);
+ for (std::size_t i = 1; i < parent_index.size(); ++i) {
+ auto p = parent_index[i];
+ if (num_child[p] > 1 || parent_index[p] == p) {
+ // parent_index[p] == p -> parent_index[i] is the soma
+ branch_index.push_back(i);
}
-
- return count;
}
- template<typename C>
- std::vector<typename C::value_type> branches(const C &parent_index)
- {
- static_assert(
- std::is_integral<typename C::value_type>::value,
- "integral type required"
- );
-
- EXPECTS(has_contiguous_segments(parent_index));
-
- auto num_child = child_count(parent_index);
- std::vector<typename C::value_type> branch_runs(
- parent_index.size(), 0
- );
-
- std::size_t num_branches = (num_child[0] == 1) ? 1 : 0;
- for (std::size_t i = 1; i < parent_index.size(); ++i) {
- auto p = parent_index[i];
- if (num_child[p] > 1) {
- ++num_branches;
- }
-
- branch_runs[i] = num_branches;
- }
+ branch_index.push_back(parent_index.size());
+ return branch_index;
+}
+
- return branch_runs;
+template<typename C>
+std::vector<typename C::value_type> expand_branches(const C& branch_index)
+{
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "integral type required"
+ );
+
+ if (branch_index.empty())
+ return {};
+
+ std::vector<typename C::value_type> expanded(branch_index.back());
+ for (std::size_t i = 0; i < branch_index.size()-1; ++i) {
+ for (auto j = branch_index[i]; j < branch_index[i+1]; ++j) {
+ expanded[j] = i;
+ }
}
- template<typename C>
- bool is_sorted(const C& c)
- {
- return std::is_sorted(c.begin(), c.end());
+ return expanded;
+}
+
+template<typename C>
+typename C::value_type find_branch(const C& branch_index,
+ typename C::value_type nid)
+{
+ using value_type = typename C::value_type;
+ static_assert(
+ std::is_integral<value_type>::value,
+ "integral type required"
+ );
+
+ auto it = std::find_if(
+ branch_index.begin(), branch_index.end(),
+ [nid](const value_type &v) { return v > nid; }
+ );
+
+ return it - branch_index.begin() - 1;
+}
+
+
+template<typename C>
+std::vector<typename C::value_type> make_parent_index(
+ const C& parent_index, const C& branch_index)
+{
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "integral type required"
+ );
+
+ if (parent_index.empty() && branch_index.empty()) {
+ return {};
}
- template<typename C>
- bool is_unique(const C& c)
- {
- return std::adjacent_find(c.begin(), c.end()) == c.end();
+ EXPECTS(parent_index.size() == branch_index.back());
+ EXPECTS(has_contiguous_segments(parent_index));
+ EXPECTS(is_strictly_monotonic_increasing(branch_index));
+
+ // expand the branch index
+ auto expanded_branch = expand_branches(branch_index);
+
+ std::vector<typename C::value_type> new_parent_index;
+ for (std::size_t i = 0; i < branch_index.size()-1; ++i) {
+ auto p = parent_index[branch_index[i]];
+ new_parent_index.push_back(expanded_branch[p]);
}
- /// Return and index that maps entries in sub to their corresponding
- /// values in super, where sub is a subset of super.
- ///
- /// Both sets are sorted and have unique entries.
- /// Complexity is O(n), where n is size of super
- template<typename C>
- // C::iterator models forward_iterator
- // C::value_type is_integral
- C index_into(const C& super, const C& sub)
- {
- //EXPECTS {s \in super : \forall s \in sub};
- EXPECTS(is_unique(super) && is_unique(sub));
- EXPECTS(is_sorted(super) && is_sorted(sub));
- EXPECTS(sub.size() <= super.size());
-
- static_assert(
- std::is_integral<typename C::value_type>::value,
- "index_into only applies to integral types"
- );
-
- C out(sub.size()); // out will have one entry for each index in sub
-
- auto sub_it=sub.begin();
- auto super_it=super.begin();
- auto sub_idx=0u, super_idx = 0u;
-
- while(sub_it!=sub.end() && super_it!=super.end()) {
- if(*sub_it==*super_it) {
- out[sub_idx] = super_idx;
- ++sub_it; ++sub_idx;
- }
- ++super_it; ++super_idx;
+ return new_parent_index;
+}
+
+
+template<typename C>
+bool is_sorted(const C& c)
+{
+ return std::is_sorted(c.begin(), c.end());
+}
+
+template<typename C>
+bool is_unique(const C& c)
+{
+ return std::adjacent_find(c.begin(), c.end()) == c.end();
+}
+
+/// Return and index that maps entries in sub to their corresponding
+/// values in super, where sub is a subset of super.
+///
+/// Both sets are sorted and have unique entries.
+/// Complexity is O(n), where n is size of super
+template<typename C>
+// C::iterator models forward_iterator
+// C::value_type is_integral
+C index_into(const C& super, const C& sub)
+{
+ //EXPECTS {s \in super : \forall s \in sub};
+ EXPECTS(is_unique(super) && is_unique(sub));
+ EXPECTS(is_sorted(super) && is_sorted(sub));
+ EXPECTS(sub.size() <= super.size());
+
+ static_assert(
+ std::is_integral<typename C::value_type>::value,
+ "index_into only applies to integral types"
+ );
+
+ C out(sub.size()); // out will have one entry for each index in sub
+
+ auto sub_it=sub.begin();
+ auto super_it=super.begin();
+ auto sub_idx=0u, super_idx = 0u;
+
+ while(sub_it!=sub.end() && super_it!=super.end()) {
+ if(*sub_it==*super_it) {
+ out[sub_idx] = super_idx;
+ ++sub_it; ++sub_idx;
}
+ ++super_it; ++super_idx;
+ }
- EXPECTS(sub_idx==sub.size());
+ EXPECTS(sub_idx==sub.size());
- return out;
- }
+ return out;
+}
} // namespace algorithms
} // namespace mc
diff --git a/src/cell_tree.hpp b/src/cell_tree.hpp
index 1e04fc2c3bcc4e0fe887c7048997afd272fa4bff..4d9f70b1a934ab3b5bdb840d6bb0bd1b8419dfec 100644
--- a/src/cell_tree.hpp
+++ b/src/cell_tree.hpp
@@ -32,7 +32,7 @@ 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 segments
- using int_type = int16_t;
+ using int_type = int;
using index_type = memory::HostVector<int_type>;
using view_type = index_type::view_type;
using const_view_type = index_type::const_view_type;
@@ -298,4 +298,3 @@ private :
} // namespace mc
} // namespace nest
-
diff --git a/src/swcio.cpp b/src/swcio.cpp
index 792d0fc7bdb3b4898b575dae552d0347efd85e9c..203f69eb1bdca8905145edb4edbe86d564682e9b 100644
--- a/src/swcio.cpp
+++ b/src/swcio.cpp
@@ -1,4 +1,5 @@
#include <algorithm>
+#include <functional>
#include <iomanip>
#include <map>
#include <sstream>
@@ -230,88 +231,60 @@ swc_record_range_clean::swc_record_range_clean(std::istream &is)
}
}
-//
-// Convenience functions for returning the radii and the coordinates of a series
-// of swc records
-//
-static std::vector<swc_record::coord_type>
-swc_radii(const std::vector<swc_record> &records)
-{
- std::vector<swc_record::coord_type> radii;
- for (const auto &r : records) {
- radii.push_back(r.radius());
- }
-
- return radii;
-}
-
-static std::vector<nest::mc::point<swc_record::coord_type> >
-swc_points(const std::vector<swc_record> &records)
-{
- std::vector<nest::mc::point<swc_record::coord_type> > points;
- for (const auto &r : records) {
- points.push_back(r.coord());
- }
-
- return points;
-}
-
-static void make_cable(cell &cell,
- const std::vector<swc_record::id_type> &branch_index,
- const std::vector<swc_record> &branch_run)
-{
- auto new_parent = branch_index[branch_run.back().id()] - 1;
- cell.add_cable(new_parent, nest::mc::segmentKind::dendrite,
- swc_radii(branch_run), swc_points(branch_run));
-}
-
cell swc_read_cell(std::istream &is)
{
+ using namespace nest::mc;
+
cell newcell;
- std::vector<swc_record::id_type> parent_list;
+ std::vector<swc_record::id_type> parent_index;
std::vector<swc_record> swc_records;
for (const auto &r : swc_get_records<swc_io_clean>(is)) {
swc_records.push_back(r);
- parent_list.push_back(r.parent());
+ parent_index.push_back(r.parent());
}
- // The parent of soma must be 0
- if (!parent_list.empty()) {
- parent_list[0] = 0;
+ if (parent_index.empty()) {
+ return newcell;
}
- auto branch_index = nest::mc::algorithms::branches(parent_list);
- std::vector<swc_record> branch_run;
-
- branch_run.reserve(parent_list.size());
- auto last_branch_point = branch_index[0];
- for (auto i = 0u; i < swc_records.size(); ++i) {
- if (branch_index[i] != last_branch_point) {
- // New branch encountered; add to cell the current one
- const auto &p = branch_run.back();
- if (p.parent() == -1) {
- // This is a soma
- newcell.add_soma(p.radius(), p.coord());
- last_branch_point = i;
- } else {
- last_branch_point = i - 1;
- make_cable(newcell, branch_index, branch_run);
- }
-
- // Reset the branch run
- branch_run.clear();
- if (p.parent() != -1) {
- // Add parent of the current cell to the branch,
- // if not branching from soma
- branch_run.push_back(swc_records[parent_list[i]]);
- }
+ // The parent of soma must be 0, while in SWC files is -1
+ parent_index[0] = 0;
+ auto branch_index = algorithms::branches(parent_index);
+ auto new_parent_index = algorithms::make_parent_index(parent_index,
+ branch_index);
+
+ // sanity check
+ EXPECTS(new_parent_index.size() == branch_index.size() - 1);
+
+ // Add the soma first; then the segments
+ newcell.add_soma(swc_records[0].radius(), swc_records[0].coord());
+ for (std::size_t i = 1; i < new_parent_index.size(); ++i) {
+ auto b_start = std::next(swc_records.begin(), branch_index[i]);
+ auto b_end = std::next(swc_records.begin(), branch_index[i+1]);
+
+ std::vector<swc_record::coord_type> radii;
+ std::vector<nest::mc::point<swc_record::coord_type>> points;
+ if (new_parent_index[i] != 0) {
+ // include the parent of current record if not branching from soma
+ auto p = parent_index[branch_index[i]];
+ radii.push_back(swc_records[p].radius());
+ points.push_back(swc_records[p].coord());
}
- branch_run.push_back(swc_records[i]);
- }
+ // extract the radii and the points
+ std::for_each(b_start, b_end,
+ [&radii](const swc_record& r) {
+ radii.push_back(r.radius());
+ });
+
+ std::for_each(b_start, b_end,
+ [&points](const swc_record& r) {
+ points.push_back(r.coord());
+ });
- if (!branch_run.empty()) {
- make_cable(newcell, branch_index, branch_run);
+ // add the new cable
+ newcell.add_cable(new_parent_index[i],
+ nest::mc::segmentKind::dendrite, radii, points);
}
return newcell;
diff --git a/src/tree.hpp b/src/tree.hpp
index 3d76a9dc1e8171f62026320a98f71078e91f1dc5..98ff9aa43656f9d68b879297a6ffd5399348b6f0 100644
--- a/src/tree.hpp
+++ b/src/tree.hpp
@@ -17,7 +17,7 @@ class tree {
public :
- using int_type = int16_t;
+ using int_type = int;
using index_type = memory::HostVector<int_type>;
using view_type = index_type::view_type;
@@ -62,92 +62,22 @@ class tree {
);
}
- // n = number of compartment in cell
- auto n = parent_index.size();
+ auto new_parent_index = algorithms::make_parent_index(
+ parent_index, algorithms::branches(parent_index));
- // On completion of this loop child_count[i] is the number of children
- // of compartment i compensate count for compartment 0, which has itself
- // as its own parent
- index_type child_count(n, 0);
- child_count[0] = -1;
- for(auto i : parent_index) {
- ++child_count[i];
- }
+ init(new_parent_index.size());
+ parents_(memory::all) = new_parent_index;
+ parents_[0] = -1;
- // Find the number of branches by summing the number of children of all
- // compartments with more than 1 child.
- auto nbranches = 1 + child_count[0]; // children of the root node are counted differently
- for(auto i : range(1,n)) {
- if(child_count[i]>1) {
- nbranches += child_count[i];
- }
- }
+ child_index_(memory::all) =
+ algorithms::make_index(algorithms::child_count(parents_));
- // allocate memory for storing the tree
- init(nbranches);
-
- // index of the branch for each compartment
- std::vector<I> branch_index(n);
-
- // Mark the parent of the root node as -1.
- // This simplifies the implementation of some algorithms on the tree
- parents_[0]=-1;
-
- // bcount records how many branches have been created in the loop
- auto bcount=0;
-
- for(auto i : range(1,n)) {
- // the branch index of the parent of compartment i
- auto parent_node = parent_index[i];
- // index of the parent of compartment i
- auto parent_branch = branch_index[parent_node];
-
- // if this compartments's parent
- // - has more than one child
- // - or is the root branch
- // this is the first compartment in a branch, so mark it as such
- if(child_count[parent_node]>1 || parent_node==0) {
- bcount++;
- branch_index[i] = bcount;
- parents_[bcount] = parent_branch;
- }
- // not the first compartment in a branch,
- // so inherit the parent's branch number
- else {
- branch_index[i] = parent_branch;
- }
- }
-
- child_index_(memory::all) = 0;
- // the root node has to be handled separately: all of its children must
- // be counted
- child_index_[1] = child_count[0];
- for(auto i : range(1, n)) {
- if(child_count[i]>1) {
- child_index_[branch_index[i]+1] = child_count[i];
- }
- }
- std::partial_sum(child_index_.begin(), child_index_.end(), child_index_.begin());
-
- // Fill in the list of children of each branch.
- // Requires some additional book keeping to keep track of how many
- // children have already been filled in for each branch.
- for(auto i : range(1, nbranches)) {
- // parents_[i] is the parent of branch i, for which i must be added
- // as a child, and child_index_[p] is the index into which the next
- // child of p is to be stored
+ std::vector<int> pos(parents_.size(), 0);
+ for (auto i = 1u; i < parents_.size(); ++i) {
auto p = parents_[i];
- children_[child_index_[p]] = i;
- ++child_index_[p];
- }
-
- // The child index has already been calculated as a side-effect of the
- // loop above, but is shifted one value to the left, so perform a
- // rotation to the right.
- for(auto i=nbranches-1; i>=00; --i) {
- child_index_[i+1] = child_index_[i];
+ children_[child_index_[p] + pos[p]] = i;
+ ++pos[p];
}
- child_index_[0] = 0;
}
size_t num_children() const {
diff --git a/src/util/meta.hpp b/src/util/meta.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2e22e31ceda0ff83588625486ec7a36161569733
--- /dev/null
+++ b/src/util/meta.hpp
@@ -0,0 +1,36 @@
+#pragma once
+
+/* Type utilities and convenience expressions. */
+
+#include <type_traits>
+
+namespace nest {
+namespace mc {
+namespace util {
+
+// Until C++14 ...
+
+template <typename T>
+using result_of_t = typename std::result_of<T>::type;
+
+template <bool V>
+using enable_if_t = typename std::enable_if<V>::type;
+
+// Convenience short cuts
+
+template <typename T>
+using enable_if_copy_constructible_t =
+ enable_if_t<std::is_copy_constructible<T>::value>;
+
+template <typename T>
+using enable_if_move_constructible_t =
+ enable_if_t<std::is_move_constructible<T>::value>;
+
+template <typename... T>
+using enable_if_constructible_t =
+ enable_if_t<std::is_constructible<T...>::value>;
+
+
+} // namespace util
+} // namespace mc
+} // namespace nest
diff --git a/src/util/optional.hpp b/src/util/optional.hpp
index a3c2da998d468c8564894f2db61b0340c056c70c..0d71c4a5d564b4a85db457a930f803388ee0861b 100644
--- a/src/util/optional.hpp
+++ b/src/util/optional.hpp
@@ -1,25 +1,25 @@
#pragma once
-/*! \file optional.h
- * \brief An option class with a monadic interface.
+/* An option class with a monadic interface.
*
- * The std::option<T> class was proposed for inclusion into C++14, but was
- * ultimately rejected. (See N3672 proposal for details.) This class offers
- * similar functionality, namely a class that can represent a value (or
- * reference), or nothing at all.
+ * The std::option<T> class was proposed for inclusion into C++14, but was
+ * ultimately rejected. (See N3672 proposal for details.) This class offers
+ * similar functionality, namely a class that can represent a value (or
+ * reference), or nothing at all.
*
- * In addition, this class offers monadic and monoidal bindings, allowing
- * the chaining of operations any one of which might represent failure with
- * an unset optional value.
+ * In addition, this class offers monadic and monoidal bindings, allowing
+ * the chaining of operations any one of which might represent failure with
+ * an unset optional value.
*
- * One point of difference between the proposal N3672 and this implementation
- * is the lack of constexpr versions of the methods and constructors.
+ * One point of difference between the proposal N3672 and this implementation
+ * is the lack of constexpr versions of the methods and constructors.
*/
#include <type_traits>
#include <stdexcept>
#include <utility>
+#include "util/meta.hpp"
#include "util/uninitialized.hpp"
namespace nest {
@@ -29,135 +29,208 @@ namespace util {
template <typename X> struct optional;
struct optional_unset_error: std::runtime_error {
- explicit optional_unset_error(const std::string &what_str): std::runtime_error(what_str) {}
- optional_unset_error(): std::runtime_error("optional value unset") {}
+ explicit optional_unset_error(const std::string &what_str)
+ : std::runtime_error(what_str)
+ {}
+
+ optional_unset_error()
+ : std::runtime_error("optional value unset")
+ {}
};
struct optional_invalid_dereference: std::runtime_error {
- explicit optional_invalid_dereference(const std::string &what_str): std::runtime_error(what_str) {}
- optional_invalid_dereference(): std::runtime_error("derefernce of optional<void> value") {}
+ explicit optional_invalid_dereference(const std::string &what_str)
+ : std::runtime_error(what_str)
+ {}
+
+ optional_invalid_dereference()
+ : std::runtime_error("derefernce of optional<void> value")
+ {}
};
struct nothing_t {};
constexpr nothing_t nothing{};
namespace detail {
- template <typename Y> struct lift_type { using type=optional<Y>; };
- template <typename Y> struct lift_type<optional<Y>> { using type=optional<Y>; };
+ template <typename Y>
+ struct lift_type {
+ using type = optional<Y>;
+ };
+
+ template <typename Y>
+ struct lift_type<optional<Y>> {
+ using type = optional<Y>;
+ };
+
+ template <typename Y>
+ using lift_type_t = typename lift_type<Y>::type;
struct optional_tag {};
- template <typename X> struct is_optional {
- enum {value=std::is_base_of<optional_tag,typename std::decay<X>::type>::value };
+ template <typename X>
+ using is_optional = std::is_base_of<optional_tag, typename std::decay<X>::type>;
+
+ template <typename D,typename X>
+ struct wrapped_type_impl {
+ using type = X;
};
- template <typename D,typename X> struct wrapped_type_ { using type=X; };
- template <typename D,typename X> struct wrapped_type_<optional<D>,X> { using type=D; };
+ template <typename D,typename X>
+ struct wrapped_type_impl<optional<D>,X> {
+ using type = D;
+ };
- template <typename X> struct wrapped_type { using type=typename wrapped_type_<typename std::decay<X>::type,X>::type; };
+ template <typename X>
+ struct wrapped_type {
+ using type = typename wrapped_type_impl<typename std::decay<X>::type,X>::type;
+ };
+
+ template <typename X>
+ using wrapped_type_t = typename wrapped_type<X>::type;
template <typename X>
struct optional_base: detail::optional_tag {
template <typename Y> friend struct optional;
protected:
- using D=util::uninitialized<X>;
+ using data_type = util::uninitialized<X>;
public:
- using reference_type=typename D::reference_type;
- using const_reference_type=typename D::const_reference_type;
- using pointer_type=typename D::pointer_type;
- using const_pointer_type=typename D::const_pointer_type;
+ using reference_type = typename data_type::reference_type;
+ using const_reference_type = typename data_type::const_reference_type;
+ using pointer_type = typename data_type::pointer_type;
+ using const_pointer_type = typename data_type::const_pointer_type;
protected:
bool set;
- D data;
+ data_type data;
- optional_base(): set(false) {}
+ optional_base() : set(false) {}
template <typename T>
- optional_base(bool set_,T&& init): set(set_) { if (set) data.construct(std::forward<T>(init)); }
+ optional_base(bool set_, T&& init) : set(set_) {
+ if (set) {
+ data.construct(std::forward<T>(init));
+ }
+ }
- reference_type ref() { return data.ref(); }
+ reference_type ref() { return data.ref(); }
const_reference_type ref() const { return data.cref(); }
public:
- ~optional_base() { if (set) data.destruct(); }
+ ~optional_base() {
+ if (set) {
+ data.destruct();
+ }
+ }
const_pointer_type operator->() const { return data.ptr(); }
- pointer_type operator->() { return data.ptr(); }
-
+ pointer_type operator->() { return data.ptr(); }
+
const_reference_type operator*() const { return ref(); }
- reference_type operator*() { return ref(); }
-
+ reference_type operator*() { return ref(); }
+
reference_type get() {
- if (set) return ref();
- else throw optional_unset_error();
+ if (set) {
+ return ref();
+ }
+ else {
+ throw optional_unset_error();
+ }
}
const_reference_type get() const {
- if (set) return ref();
- else throw optional_unset_error();
+ if (set) {
+ return ref();
+ }
+ else {
+ throw optional_unset_error();
+ }
}
explicit operator bool() const { return set; }
template <typename Y>
- bool operator==(const Y &y) const { return set && ref()==y; }
+ bool operator==(const Y& y) const {
+ return set && ref()==y;
+ }
template <typename Y>
- bool operator==(const optional<Y> &o) const {
+ bool operator==(const optional<Y>& o) const {
return (set && o.set && ref()==o.ref()) || (!set && !o.set);
}
void reset() {
- if (set) data.destruct();
- set=false;
+ if (set) {
+ data.destruct();
+ }
+ set = false;
}
template <typename F>
- auto bind(F &&f) -> typename lift_type<decltype(data.apply(std::forward<F>(f)))>::type {
- using F_result_type=decltype(data.apply(std::forward<F>(f)));
- using result_type=typename lift_type<F_result_type>::type;
+ auto bind(F&& f) -> lift_type_t<decltype(data.apply(std::forward<F>(f)))> {
+ using F_result_type = decltype(data.apply(std::forward<F>(f)));
+ using result_type = lift_type_t<F_result_type>;
+
+ if (!set) {
+ return result_type();
+ }
- if (!set) return result_type();
- else return bind_impl<result_type,std::is_same<F_result_type,void>::value>::bind(data,std::forward<F>(f));
+ return bind_impl<result_type, std::is_void<F_result_type>::value>::
+ bind(data, std::forward<F>(f));
}
template <typename F>
- auto bind(F &&f) const -> typename lift_type<decltype(data.apply(std::forward<F>(f)))>::type {
- using F_result_type=decltype(data.apply(std::forward<F>(f)));
- using result_type=typename lift_type<F_result_type>::type;
+ auto bind(F&& f) const -> lift_type_t<decltype(data.apply(std::forward<F>(f)))> {
+ using F_result_type = decltype(data.apply(std::forward<F>(f)));
+ using result_type = lift_type_t<F_result_type>;
- if (!set) return result_type();
- else return bind_impl<result_type,std::is_same<F_result_type,void>::value>::bind(data,std::forward<F>(f));
+ if (!set) {
+ return result_type();
+ }
+
+ return bind_impl<result_type, std::is_void<F_result_type>::value>::
+ bind(data, std::forward<F>(f));
}
template <typename F>
- auto operator>>(F &&f) -> decltype(this->bind(std::forward<F>(f))) { return bind(std::forward<F>(f)); }
+ auto operator>>(F&& f) -> decltype(this->bind(std::forward<F>(f))) {
+ return bind(std::forward<F>(f));
+ }
template <typename F>
- auto operator>>(F &&f) const -> decltype(this->bind(std::forward<F>(f))) { return bind(std::forward<F>(f)); }
+ auto operator>>(F&& f) const -> decltype(this->bind(std::forward<F>(f))) {
+ return bind(std::forward<F>(f));
+ }
private:
- template <typename R,bool F_void_return>
+ template <typename R, bool F_void_return>
struct bind_impl {
template <typename DT,typename F>
- static R bind(DT &d,F &&f) { return R(d.apply(std::forward<F>(f))); }
+ static R bind(DT& d,F&& f) {
+ return R(d.apply(std::forward<F>(f)));
+ }
};
template <typename R>
struct bind_impl<R,true> {
template <typename DT,typename F>
- static R bind(DT &d,F &&f) { d.apply(std::forward<F>(f)); return R(true); }
+ static R bind(DT& d,F&& f) {
+ d.apply(std::forward<F>(f));
+ return R(true);
+ }
};
-
};
-}
+
+ // type utilities
+ template <typename T>
+ using enable_unless_optional_t = enable_if_t<!is_optional<T>::value>;
+
+} // namespace detail
template <typename X>
struct optional: detail::optional_base<X> {
- using base=detail::optional_base<X>;
+ using base = detail::optional_base<X>;
using base::set;
using base::ref;
using base::reset;
@@ -166,85 +239,109 @@ struct optional: detail::optional_base<X> {
optional(): base() {}
optional(nothing_t): base() {}
- template <typename Y=X,typename = typename std::enable_if<std::is_copy_constructible<Y>::value>::type>
- optional(const X &x): base(true,x) {}
+ template <
+ typename Y = X,
+ typename = enable_if_copy_constructible_t<Y>
+ >
+ optional(const X& x): base(true, x) {}
- template <typename Y=X,typename = typename std::enable_if<std::is_move_constructible<Y>::value>::type>
- optional(X &&x): base(true,std::move(x)) {}
+ template <
+ typename Y = X,
+ typename = enable_if_move_constructible_t<Y>
+ >
+ optional(X&& x): base(true, std::move(x)) {}
- optional(const optional &ot): base(ot.set,ot.ref()) {}
+ optional(const optional& ot): base(ot.set, ot.ref()) {}
template <typename T>
- optional(const optional<T> &ot): base(ot.set,ot.ref()) {}
+ optional(const optional<T>& ot): base(ot.set, ot.ref()) {}
- optional(optional &&ot): base(ot.set,std::move(ot.ref())) {}
+ optional(optional&& ot): base(ot.set, std::move(ot.ref())) {}
template <typename T>
- optional(optional<T> &&ot): base(ot.set,std::move(ot.ref())) {}
+ optional(optional<T>&& ot): base(ot.set, std::move(ot.ref())) {}
- template <typename Y,typename = typename std::enable_if<!detail::is_optional<Y>::value>::type>
- optional &operator=(Y &&y) {
- if (set) ref()=std::forward<Y>(y);
+ template <
+ typename Y,
+ typename = detail::enable_unless_optional_t<Y>
+ >
+ optional& operator=(Y&& y) {
+ if (set) {
+ ref() = std::forward<Y>(y);
+ }
else {
- set=true;
+ set = true;
data.construct(std::forward<Y>(y));
}
return *this;
}
- optional &operator=(const optional &o) {
+ optional& operator=(const optional& o) {
if (set) {
- if (o.set) ref()=o.ref();
- else reset();
+ if (o.set) {
+ ref() = o.ref();
+ }
+ else {
+ reset();
+ }
}
else {
- set=o.set;
- if (set) data.construct(o.ref());
+ set = o.set;
+ if (set) {
+ data.construct(o.ref());
+ }
}
return *this;
}
- template <typename Y=X, typename =typename std::enable_if<
- std::is_move_assignable<Y>::value &&
- std::is_move_constructible<Y>::value
- >::type>
- optional &operator=(optional &&o) {
+ template <
+ typename Y = X,
+ typename = typename std::enable_if<
+ std::is_move_assignable<Y>::value &&
+ std::is_move_constructible<Y>::value
+ >::type
+ >
+ optional& operator=(optional&& o) {
if (set) {
- if (o.set) ref()=std::move(o.ref());
+ if (o.set) {
+ ref() = std::move(o.ref());
+ }
else reset();
}
else {
- set=o.set;
- if (set) data.construct(std::move(o.ref()));
+ set = o.set;
+ if (set) {
+ data.construct(std::move(o.ref()));
+ }
}
return *this;
}
};
template <typename X>
-struct optional<X &>: detail::optional_base<X &> {
- using base=detail::optional_base<X &>;
+struct optional<X&>: detail::optional_base<X&> {
+ using base=detail::optional_base<X&>;
using base::set;
using base::ref;
using base::data;
optional(): base() {}
optional(nothing_t): base() {}
- optional(X &x): base(true,x) {}
+ optional(X&x): base(true,x) {}
template <typename T>
- optional(optional<T &> &ot): base(ot.set,ot.ref()) {}
+ optional(optional<T&>& ot): base(ot.set,ot.ref()) {}
- template <typename Y,typename = typename std::enable_if<!detail::is_optional<Y>::value>::type>
- optional &operator=(Y &y) {
- set=true;
- ref()=y;
+ template <typename Y,typename = typename std::enable_if<!detail::is_optional<Y>()>::type>
+ optional& operator=(Y& y) {
+ set = true;
+ ref() = y;
return *this;
}
template <typename Y>
- optional &operator=(optional<Y &> &o) {
- set=o.set;
+ optional& operator=(optional<Y&>& o) {
+ set = o.set;
data.construct(o);
return *this;
}
@@ -256,7 +353,7 @@ struct optional<X &>: detail::optional_base<X &> {
template <>
struct optional<void>: detail::optional_base<void> {
- using base=detail::optional_base<void>;
+ using base = detail::optional_base<void>;
using base::set;
optional(): base() {}
@@ -265,19 +362,24 @@ struct optional<void>: detail::optional_base<void> {
optional(T): base(true,true) {}
template <typename T>
- optional(const optional<T> &o): base(o.set,true) {}
+ optional(const optional<T>& o): base(o.set,true) {}
template <typename T>
- optional &operator=(T) { set=true; return *this; }
+ optional& operator=(T) {
+ set = true;
+ return *this;
+ }
template <typename T>
- optional &operator=(const optional<T> &o) { set=o.set; return *this; }
+ optional& operator=(const optional<T>& o) {
+ set = o.set;
+ return *this;
+ }
- // override equality operators
template <typename Y>
- bool operator==(const Y &y) const { return false; }
+ bool operator==(const Y& y) const { return false; }
- bool operator==(const optional<void> &o) const {
+ bool operator==(const optional<void>& o) const {
return (set && o.set) || (!set && !o.set);
}
};
@@ -286,25 +388,36 @@ struct optional<void>: detail::optional_base<void> {
template <typename A,typename B>
typename std::enable_if<
detail::is_optional<A>::value || detail::is_optional<B>::value,
- optional<typename std::common_type<typename detail::wrapped_type<A>::type,typename detail::wrapped_type<B>::type>::type>
+ optional<
+ typename std::common_type<
+ detail::wrapped_type_t<A>,
+ detail::wrapped_type_t<B>
+ >::type
+ >
>::type
-operator|(A &&a,B &&b) {
- return a?a:b;
+operator|(A&& a,B&& b) {
+ return a ? a : b;
}
template <typename A,typename B>
typename std::enable_if<
detail::is_optional<A>::value || detail::is_optional<B>::value,
- optional<typename detail::wrapped_type<B>::type>
+ optional<detail::wrapped_type_t<B>>
>::type
-operator&(A &&a,B &&b) {
- using result_type=optional<typename detail::wrapped_type<B>::type>;
- return a?b:result_type();
+operator&(A&& a,B&& b) {
+ using result_type = optional<detail::wrapped_type_t<B>>;
+ return a ? b: result_type();
}
-inline optional<void> provided(bool condition) { return condition?optional<void>(true):optional<void>(); }
+inline optional<void> provided(bool condition) {
+ return condition ? optional<void>(true) : optional<void>();
+}
template <typename X>
-optional<X> just(X &&x) { return optional<X>(std::forward<X>(x)); }
+optional<X> just(X&& x) {
+ return optional<X>(std::forward<X>(x));
+}
-}}} // namespace nest::mc::util
+} // namespace util
+} // namespace mc
+} // namespace nest
diff --git a/src/util/uninitialized.hpp b/src/util/uninitialized.hpp
index 26ebcf90ef2c94ea0645ccb7d3a94f6a01fb2510..846e46d5a2451fe1535f026ee4b801e5303a6f0c 100644
--- a/src/util/uninitialized.hpp
+++ b/src/util/uninitialized.hpp
@@ -5,10 +5,15 @@
* The uninitialized<X> structure holds space for an item of
* type X, leaving its construction or destruction to the user.
*
- * Specialisations for reference types X & and for the void type
+ * Specialisations for reference types X& and for the void type
* allow for the handling of non-value types in a uniform manner.
*/
+#include <type_traits>
+#include <utility>
+
+#include "util/meta.hpp"
+
namespace nest {
namespace mc {
namespace util {
@@ -19,39 +24,47 @@ namespace util {
template <typename X>
struct uninitialized {
private:
- typename std::aligned_storage<sizeof(X),alignof(X)>::type data;
+ typename std::aligned_storage<sizeof(X), alignof(X)>::type data;
public:
- using pointer_type=X *;
- using const_pointer_type=const X *;
- using reference_type=X &;
- using const_reference_type=const X &;
+ using pointer_type = X*;
+ using const_pointer_type = const X*;
+ using reference_type = X&;
+ using const_reference_type = const X&;
- pointer_type ptr() { return reinterpret_cast<X *>(&data); }
- const_pointer_type cptr() const { return reinterpret_cast<const X *>(&data); }
+ pointer_type ptr() { return reinterpret_cast<X*>(&data); }
+ const_pointer_type cptr() const { return reinterpret_cast<const X*>(&data); }
- reference_type ref() { return *reinterpret_cast<X *>(&data); }
- const_reference_type cref() const { return *reinterpret_cast<const X *>(&data); }
+ reference_type ref() { return *reinterpret_cast<X*>(&data); }
+ const_reference_type cref() const { return *reinterpret_cast<const X*>(&data); }
// Copy construct the value.
- template <typename Y=X,
- typename =typename std::enable_if<std::is_copy_constructible<Y>::value>::type>
- void construct(const X &x) { new(&data) X(x); }
+ template <
+ typename Y = X,
+ typename = enable_if_copy_constructible_t<Y>
+ >
+ void construct(const X& x) {
+ new(&data) X(x);
+ }
// General constructor for X, forwarding arguments.
- template <typename... Y,
- typename =typename std::enable_if<std::is_constructible<X,Y...>::value>::type>
- void construct(Y&& ...args) { new(&data) X(std::forward<Y>(args)...); }
+ template <
+ typename... Y,
+ typename = enable_if_constructible_t<X, Y...>
+ >
+ void construct(Y&& ...args) {
+ new(&data) X(std::forward<Y>(args)...);
+ }
void destruct() { ptr()->~X(); }
// Apply the one-parameter functor F to the value by reference.
template <typename F>
- typename std::result_of<F(reference_type)>::type apply(F &&f) { return f(ref()); }
+ result_of_t<F(reference_type)> apply(F&& f) { return f(ref()); }
// Apply the one-parameter functor F to the value by const reference.
template <typename F>
- typename std::result_of<F(const_reference_type)>::type apply(F &&f) const { return f(cref()); }
+ result_of_t<F(const_reference_type)> apply(F&& f) const { return f(cref()); }
};
/* Maintains storage for a pointer of type X, representing
@@ -63,10 +76,10 @@ private:
X *data;
public:
- using pointer_type=X *;
- using const_pointer_type=const X *;
- using reference_type=X &;
- using const_reference_type=const X &;
+ using pointer_type = X*;
+ using const_pointer_type = const X*;
+ using reference_type = X&;
+ using const_reference_type = const X&;
pointer_type ptr() { return data; }
const_pointer_type cptr() const { return data; }
@@ -74,15 +87,20 @@ public:
reference_type ref() { return *data; }
const_reference_type cref() const { return *data; }
- void construct(X &x) { data=&x; }
+ void construct(X& x) { data = &x; }
void destruct() {}
// Apply the one-parameter functor F to the value by reference.
template <typename F>
- typename std::result_of<F(reference_type)>::type apply(F &&f) { return f(ref()); }
+ result_of_t<F(reference_type)> apply(F&& f) {
+ return f(ref());
+ }
+
// Apply the one-parameter functor F to the value by const reference.
template <typename F>
- typename std::result_of<F(const_reference_type)>::type apply(F &&f) const { return f(cref()); }
+ result_of_t<F(const_reference_type)> apply(F&& f) const {
+ return f(cref());
+ }
};
/* Wrap a void type in an uninitialized template.
@@ -91,10 +109,10 @@ public:
*/
template <>
struct uninitialized<void> {
- using pointer_type=void *;
- using const_pointer_type=const void *;
- using reference_type=void;
- using const_reference_type=void;
+ using pointer_type = void*;
+ using const_pointer_type = const void*;
+ using reference_type = void;
+ using const_reference_type = void;
pointer_type ptr() { return nullptr; }
const_pointer_type cptr() const { return nullptr; }
@@ -109,20 +127,9 @@ struct uninitialized<void> {
// Equivalent to f()
template <typename F>
- typename std::result_of<F()>::type apply(F &&f) const { return f(); }
+ result_of_t<F()> apply(F&& f) const { return f(); }
};
-template <typename...>
-struct uninitialized_can_construct: std::false_type {};
-
-template <typename X,typename... Y>
-struct uninitialized_can_construct<X,Y...>: std::integral_constant<bool,std::is_constructible<X,Y...>::value> {};
-
-template <typename X,typename Y>
-struct uninitialized_can_construct<X &,Y>: std::integral_constant<bool,std::is_convertible<X &,Y>::value> {};
-
-template <typename... Y>
-struct uninitialized_can_construct<void,Y...>: std::true_type {};
-
-}}} // namespace nest::mc::util
-
+} // namespace util
+} // namespace mc
+} // namespace nest
diff --git a/tests/test_algorithms.cpp b/tests/test_algorithms.cpp
index 584bdf8ad10db170b44604d3db8f8047d3c29860..ec25bd42c0941b8324fbea3288ed5b4f1d1f7a59 100644
--- a/tests/test_algorithms.cpp
+++ b/tests/test_algorithms.cpp
@@ -369,11 +369,11 @@ TEST(algorithms, branches)
{
//
- // 0
- // /|\.
- // 1 4 6
- // / | \.
- // 2 5 7
+ // 0 0
+ // /|\. /|\.
+ // 1 4 6 1 2 3
+ // / | \. => / \.
+ // 2 5 7 4 5
// / \.
// 3 8
// / \.
@@ -386,58 +386,100 @@ TEST(algorithms, branches)
std::vector<int> parent_index =
{ 0, 0, 1, 2, 0, 4, 0, 6, 7, 8, 9, 8, 11, 12 };
std::vector<int> expected_branches =
- { 0, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5 };
+ { 0, 1, 4, 6, 9, 11, 14 };
+ std::vector<int> expected_parent_index =
+ { 0, 0, 0, 0, 3, 3 };
auto actual_branches = algorithms::branches(parent_index);
EXPECT_EQ(expected_branches, actual_branches);
+
+ auto actual_parent_index =
+ algorithms::make_parent_index(parent_index, actual_branches);
+ EXPECT_EQ(expected_parent_index, actual_parent_index);
+
+ // Check find_branch
+ EXPECT_EQ(0, algorithms::find_branch(actual_branches, 0));
+ EXPECT_EQ(1, algorithms::find_branch(actual_branches, 1));
+ EXPECT_EQ(1, algorithms::find_branch(actual_branches, 2));
+ EXPECT_EQ(1, algorithms::find_branch(actual_branches, 3));
+ EXPECT_EQ(2, algorithms::find_branch(actual_branches, 4));
+ EXPECT_EQ(2, algorithms::find_branch(actual_branches, 5));
+ EXPECT_EQ(3, algorithms::find_branch(actual_branches, 6));
+ EXPECT_EQ(3, algorithms::find_branch(actual_branches, 7));
+ EXPECT_EQ(3, algorithms::find_branch(actual_branches, 8));
+ EXPECT_EQ(4, algorithms::find_branch(actual_branches, 9));
+ EXPECT_EQ(4, algorithms::find_branch(actual_branches, 10));
+ EXPECT_EQ(5, algorithms::find_branch(actual_branches, 11));
+ EXPECT_EQ(5, algorithms::find_branch(actual_branches, 12));
+ EXPECT_EQ(5, algorithms::find_branch(actual_branches, 13));
+ EXPECT_EQ(6, algorithms::find_branch(actual_branches, 55));
+
+ // Check expand_branches
+ auto expanded = algorithms::expand_branches(actual_branches);
+ EXPECT_EQ(parent_index.size(), expanded.size());
+ for (std::size_t i = 0; i < parent_index.size(); ++i) {
+ EXPECT_EQ(algorithms::find_branch(actual_branches, i),
+ expanded[i]);
+ }
}
{
//
- // 0
- // |
- // 1
+ // 0 0
+ // | |
+ // 1 => 1
// |
// 2
// |
// 3
//
- std::vector<int> parent_index =
- { 0, 0, 1, 2 };
- std::vector<int> expected_branches =
- { 0, 1, 1, 1 };
+ std::vector<int> parent_index = { 0, 0, 1, 2 };
+ std::vector<int> expected_branches = { 0, 1, 4 };
+ std::vector<int> expected_parent_index = { 0, 0 };
auto actual_branches = algorithms::branches(parent_index);
EXPECT_EQ(expected_branches, actual_branches);
+
+ auto actual_parent_index =
+ algorithms::make_parent_index(parent_index, actual_branches);
+ EXPECT_EQ(expected_parent_index, actual_parent_index);
}
{
//
- // 0
- // |
- // 1
- // |
- // 2
+ // 0 0
+ // | |
+ // 1 => 1
+ // | / \.
+ // 2 2 3
// / \.
// 3 4
// \.
// 5
//
- std::vector<int> parent_index =
- { 0, 0, 1, 2, 2, 4 };
- std::vector<int> expected_branches =
- { 0, 1, 1, 2, 3, 3 };
+ std::vector<int> parent_index = { 0, 0, 1, 2, 2, 4 };
+ std::vector<int> expected_branches = { 0, 1, 3, 4, 6 };
+ std::vector<int> expected_parent_index = { 0, 0, 1, 1 };
auto actual_branches = algorithms::branches(parent_index);
EXPECT_EQ(expected_branches, actual_branches);
+
+ auto actual_parent_index =
+ algorithms::make_parent_index(parent_index, actual_branches);
+ EXPECT_EQ(expected_parent_index, actual_parent_index);
}
{
- std::vector<int> parent_index = { 0 };
- std::vector<int> expected_branches = { 0 };
+ std::vector<int> parent_index = { 0 };
+ std::vector<int> expected_branches = { 0, 1 };
+ std::vector<int> expected_parent_index = { 0 };
auto actual_branches = algorithms::branches(parent_index);
EXPECT_EQ(expected_branches, actual_branches);
+
+ auto actual_parent_index =
+ algorithms::make_parent_index(parent_index, actual_branches);
+ EXPECT_EQ(expected_parent_index, actual_parent_index);
}
}
diff --git a/tests/test_swcio.cpp b/tests/test_swcio.cpp
index 62c52b099e6142f9bc409ca461568dbce7ffa4d7..5ba0b991a865dfb397e5c1dbbb1e8813e71b628c 100644
--- a/tests/test_swcio.cpp
+++ b/tests/test_swcio.cpp
@@ -521,4 +521,3 @@ TEST(swc_parser, from_file_ball_and_stick)
EXPECT_TRUE(nest::mc::cell_basic_equality(local_cell, cell));
}
-
diff --git a/tests/test_tree.cpp b/tests/test_tree.cpp
index 1e06b3a90fdf92f96e43e2c93036d37788c3efa6..ca9c6d73b1d9e4bde0bb90fc9061ff34f69f0503 100644
--- a/tests/test_tree.cpp
+++ b/tests/test_tree.cpp
@@ -142,20 +142,27 @@ TEST(cell_tree, from_parent_index) {
{
//
// 0
- // / \.
- // 1 2
+ // /|\.
+ // 1 4 5
// / \.
- // 3 4
- std::vector<int> parent_index = {0,0,0,1,1};
+ // 2 3
+ std::vector<int> parent_index = {0,0,1,1,0,0};
cell_tree tree(parent_index);
- EXPECT_EQ(tree.num_segments(), 5u);
+ EXPECT_EQ(tree.num_segments(), 6u);
- EXPECT_EQ(tree.num_children(0), 2u);
+ EXPECT_EQ(tree.num_children(0), 3u);
EXPECT_EQ(tree.num_children(1), 2u);
EXPECT_EQ(tree.num_children(2), 0u);
EXPECT_EQ(tree.num_children(3), 0u);
EXPECT_EQ(tree.num_children(4), 0u);
+
+ // Check children
+ EXPECT_EQ(1, tree.children(0)[0]);
+ EXPECT_EQ(4, tree.children(0)[1]);
+ EXPECT_EQ(5, tree.children(0)[2]);
+ EXPECT_EQ(2, tree.children(1)[0]);
+ EXPECT_EQ(3, tree.children(1)[1]);
}
{
// 0