Repository: higra/Higra (language: python)

Issue:
Memory usage of the tree structure
The tree structure now has a default overhead of ``4n * 8`` bytes to store children (with ``n`` the number of vertices). Removing this overhead would allow to process larger trees if the children information is not needed. 

A possible solution is to introduce a lazy computation strategy: children would be computed only if required and then stored as usual. In Python this could be completely transparent  to users  but an explicit call to a function requesting children in C++ might be required.


Discussion / hints:

Fix patch:
diff --git a/.travis.yml b/.travis.yml
index cc382f2c..57a9aa0f 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -78,7 +78,7 @@ env:
   - CIBW_BEFORE_BUILD_MACOS="source tools/cibuildwheel_osx.sh"
   - CIBW_ENVIRONMENT_LINUX='HG_USE_TBB=1 MFLAG="-m64"  CXXFLAGS="${MFLAG}" TBB_INCLUDE_DIR="/tbb/include/" TBB_LIBRARY="/tbb/lib/"'
   - CIBW_ENVIRONMENT_MACOS='HG_USE_TBB=1 TBB_INCLUDE_DIR="/Users/travis/tbb/include/" TBB_LIBRARY="/Users/travis/tbb/lib/"'
-  - CIBW_TEST_COMMAND="pip install scikit-learn && PYTHONMALLOC=malloc_debug python -c \"import unittest;result=unittest.TextTestRunner(verbosity=2).run(unittest.defaultTestLoader.discover('{project}/test/python/'));exit(0 if result.wasSuccessful() else 1)\""
+  - CIBW_TEST_COMMAND="pip install scikit-learn==0.23.2 && python -c \"import unittest;result=unittest.TextTestRunner(verbosity=2).run(unittest.defaultTestLoader.discover('{project}/test/python/'));exit(0 if result.wasSuccessful() else 1)\""
   - TWINE_USERNAME="bperret"
   - secure: eIQeQdE8S5hQuj2xOlfYLz2KexLd3GLN/vgj9Gyik1VxJUpTysL5Qp2Qufz57qcHBYH2bvqvRM5VfkSorlQVpjWZHweQMP8qu1MHyOh5NSjWUdaBQSMOrtHkv7IsCfGOB3rwYtHgIM+GCIxHmbmuKx3xayiwLqMDXaSg9GzMGg0FDy4M4o31zoM9eOhtHpzW0SJyN3xmZgcHHgvq7UFy32tvKSfJ3O3Tjv/ju3PqzlIvJ8zdTGl0DLep7TSWaw7XaBDUR84et0t2+7+/TZz8SbPGgB5+U+a0yNTsmiEbwE70EaXMTomyfXLGWh59+cSZEOqGk+3yp0IQPaYkxaKom69+y678I1rkOonikVzRlZjY9qW8MOnppEswGQrOOGhPzbT+PlQaX0ZN8rHfsdPeYHDimG4xcblKgRW/6JlB/3ZjU5L/uMVQssFXSZ8hjoF4ucf6u1E62KEhY7qanrXujOBlXyLIMoo/pznZny154NxQM068vmIYH7neIQlBtGQPWK5DpeXf5azOfE7ehJr2weBwvNl6nb4gkvJ3mrMHD62NL5LjKqASctus2zVWJMo7EZUb5IS6WkB8Dn+IrCBwRGeH7nyLN4zP5CsJWtCwmE/SiZc0AjUAAtxRZhHS2oW48R1Bv1zzegb/TWXguXJL9DpeqynSx0eEd0j0bSE1orQ=
 before_install:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index c81fa849..3d0b2d94 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -77,6 +77,12 @@ message("Current compiler is: ${CMAKE_CXX_COMPILER_ID}")
 
 string(TOUPPER "${CMAKE_BUILD_TYPE}" U_CMAKE_BUILD_TYPE)
 
+find_program(CCACHE_FOUND "ccache")
+if(CCACHE_FOUND)
+    message("Using ccache")
+    set_property(GLOBAL PROPERTY RULE_LAUNCH_COMPILE ccache)
+endif(CCACHE_FOUND)
+
 if (CMAKE_CXX_COMPILER_ID MATCHES "Clang" OR CMAKE_CXX_COMPILER_ID MATCHES "GNU" OR CMAKE_CXX_COMPILER_ID MATCHES "Intel")
     SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wno-unused-local-typedefs")# -Wno-unused-but-set-variable -Wno-missing-braces -Wno-mismatched-tags -Wno-unneeded-internal-declaration")
 
@@ -93,7 +99,7 @@ if (CMAKE_CXX_COMPILER_ID MATCHES "Clang" OR CMAKE_CXX_COMPILER_ID MATCHES "GNU"
 
     # Not debug => optim flag
     if (NOT ${U_CMAKE_BUILD_TYPE} MATCHES DEBUG)
-        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}  -ffast-math -Wno-strict-overflow")#-march=native
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}  -ffast-math -fno-finite-math-only -Wno-strict-overflow")#-march=native
 
         # check link time optimization -faligned-new
         CHECK_CXX_COMPILER_FLAG("-flto" HAS_LTO_FLAG)
@@ -164,6 +170,7 @@ option(HG_USE_TBB
 if (HG_USE_TBB)
     find_package(TBB REQUIRED)
     message(STATUS "Found intel TBB: ${TBB_INCLUDE_DIRS}")
+    add_definitions("-DTBB_SUPPRESS_DEPRECATED_MESSAGES")
 endif ()
 
 option(HG_BUILD_WHEEL
diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt
index e239eb06..f89a7f95 100644
--- a/benchmark/CMakeLists.txt
+++ b/benchmark/CMakeLists.txt
@@ -69,7 +69,9 @@ include_directories(${GBENCHMARK_INCLUDE_DIRS})
 
 set(FILES_BENCHMARK
         main.cpp
-        benchmark_parallel_sort.cpp
+        utils.cpp
+        benchmark_accumulator.cpp
+        #benchmark_parallel_sort.cpp
         # benchmark_tree_iterator.cpp
         #benchmark_array_accessor.cpp
         #benchmark_views.cpp
diff --git a/benchmark/benchmark_accumulator.cpp b/benchmark/benchmark_accumulator.cpp
new file mode 100644
index 00000000..939c16c8
--- /dev/null
+++ b/benchmark/benchmark_accumulator.cpp
@@ -0,0 +1,49 @@
+/***************************************************************************
+* Copyright ESIEE Paris (2021)                                             *
+*                                                                          *
+* Contributor(s) : Benjamin Perret                                         *
+*                                                                          *
+* Distributed under the terms of the CECILL-B License.                     *
+*                                                                          *
+* The full license is in the file LICENSE, distributed with this software. *
+****************************************************************************/
+
+
+#include <benchmark/benchmark.h>
+#include "utils.h"
+#include "higra/structure/array.hpp"
+#include "higra/accumulator/tree_accumulator.hpp"
+#include "xtensor/xrandom.hpp"
+
+using namespace xt;
+using namespace hg;
+
+
+std::size_t min_tree_size = 10;
+std::size_t max_tree_size = 20;
+
+
+static void BM_tree_accumulator(benchmark::State &state) {
+    for (auto _ : state) {
+        state.PauseTiming();
+
+        std::size_t size = state.range(0);
+        xt::random::seed(42);
+
+        auto t = get_complete_binary_tree(size);
+        array_1d<int> area = xt::ones<int>({t.num_leaves()});
+        //auto altitude = eval(xt::arange<double>(t.num_vertices())/(t.num_vertices()*255.));
+
+        state.ResumeTiming();
+        auto res = hg::accumulate_sequential(t, area, hg::accumulator_sum());
+
+        benchmark::DoNotOptimize(res[t.root()]);
+    }
+}
+
+BENCHMARK(BM_tree_accumulator)->Range(1 << min_tree_size, 1 << max_tree_size);
+
+
+
+
+
diff --git a/benchmark/benchmark_tree_attributes.cpp b/benchmark/benchmark_tree_attributes.cpp
index c4f51890..190ad26a 100644
--- a/benchmark/benchmark_tree_attributes.cpp
+++ b/benchmark/benchmark_tree_attributes.cpp
@@ -27,21 +27,6 @@ using namespace hg;
 std::size_t min_tree_size = 10;
 std::size_t max_tree_size = 16;
 
-bool is_power_of_two(std::size_t x) {
-    return (x != 0) && ((x & (x - 1)) == 0);
-}
-
-auto get_complete_binary_tree(std::size_t num_leaves) {
-    assert(is_power_of_two(num_leaves));
-    array_1d<std::size_t> parent = array_1d<std::size_t>::from_shape({num_leaves * 2 - 1});
-    for (std::size_t i = 0, j = num_leaves; i < parent.size() - 1; j++) {
-        parent(i++) = j;
-        parent(i++) = j;
-    }
-    parent(parent.size() - 1) = parent.size() - 1;
-    return tree(parent);
-}
-
 static void BM_tree_volume_cstyle(benchmark::State &state) {
     for (auto _ : state) {
         state.PauseTiming();
diff --git a/benchmark/benchmark_tree_iterator.cpp b/benchmark/benchmark_tree_iterator.cpp
index c162c874..62c5c7b5 100644
--- a/benchmark/benchmark_tree_iterator.cpp
+++ b/benchmark/benchmark_tree_iterator.cpp
@@ -9,6 +9,7 @@
 ****************************************************************************/
 
 #include <benchmark/benchmark.h>
+#include "utils.h"
 
 #include "higra/graph.hpp"
 #include "higra/accumulator/tree_accumulator.hpp"
@@ -26,20 +27,7 @@ using namespace hg;
 std::size_t min_tree_size = 10;
 std::size_t max_tree_size = 16;
 
-bool is_power_of_two(std::size_t x) {
-    return (x != 0) && ((x & (x - 1)) == 0);
-}
 
-auto get_complete_binary_tree(std::size_t num_leaves) {
-    assert(is_power_of_two(num_leaves));
-    array_1d<std::size_t> parent = array_1d<std::size_t>::from_shape({num_leaves * 2 - 1});
-    for (std::size_t i = 0, j = num_leaves; i < parent.size() - 1; j++) {
-        parent(i++) = j;
-        parent(i++) = j;
-    }
-    parent(parent.size() - 1) = parent.size() - 1;
-    return tree(parent);
-}
 
 static void BM_tree_accumulate_parallel_scalar_cstyle(benchmark::State &state) {
     for (auto _ : state) {
diff --git a/benchmark/utils.cpp b/benchmark/utils.cpp
new file mode 100644
index 00000000..d8112d6d
--- /dev/null
+++ b/benchmark/utils.cpp
@@ -0,0 +1,28 @@
+/***************************************************************************
+* Copyright ESIEE Paris (2021)                                             *
+*                                                                          *
+* Contributor(s) : Benjamin Perret                                         *
+*                                                                          *
+* Distributed under the terms of the CECILL-B License.                     *
+*                                                                          *
+* The full license is in the file LICENSE, distributed with this software. *
+****************************************************************************/
+
+#include "utils.h"
+
+using namespace hg;
+
+bool is_power_of_two(std::size_t x) {
+    return (x != 0) && ((x & (x - 1)) == 0);
+}
+
+hg::tree get_complete_binary_tree(std::size_t num_leaves) {
+    assert(is_power_of_two(num_leaves));
+    array_1d<std::size_t> parent = array_1d<std::size_t>::from_shape({num_leaves * 2 - 1});
+    for (std::size_t i = 0, j = num_leaves; i < parent.size() - 1; j++) {
+        parent(i++) = j;
+        parent(i++) = j;
+    }
+    parent(parent.size() - 1) = parent.size() - 1;
+    return tree(std::move(parent));
+}
\ No newline at end of file
diff --git a/benchmark/utils.h b/benchmark/utils.h
new file mode 100644
index 00000000..bbd0f57a
--- /dev/null
+++ b/benchmark/utils.h
@@ -0,0 +1,14 @@
+/***************************************************************************
+* Copyright ESIEE Paris (2021)                                             *
+*                                                                          *
+* Contributor(s) : Benjamin Perret                                         *
+*                                                                          *
+* Distributed under the terms of the CECILL-B License.                     *
+*                                                                          *
+* The full license is in the file LICENSE, distributed with this software. *
+****************************************************************************/
+
+#pragma once
+#include "higra/graph.hpp"
+
+hg::tree get_complete_binary_tree(std::size_t num_leaves);
diff --git a/doc/source/graph.rst b/doc/source/graph.rst
index abbd4217..24f5bfa9 100644
--- a/doc/source/graph.rst
+++ b/doc/source/graph.rst
@@ -20,22 +20,25 @@ All functions acting on graphs have the same name in C++ and in Python, except f
 In c++, graph related methods are free functions (as in BGL), while in Python they are member functions.
 For example, the function ``num_vertices`` that returns the number of vertices in a graph, will be called:
 
+.. important::
+
+    In c++, one must call the function ``compute_children`` on a tree before using it as a graph.
 
 .. tabs::
 
-    .. tab:: c++
+    .. tab:: python
 
-        .. code-block:: cpp
+        .. code-block:: python
             :linenos:
 
-            auto n = num_vertices(graph);
+            n = graph.num_vertices()
 
-    .. tab:: python
+    .. tab:: c++
 
-        .. code-block:: python
+        .. code-block:: cpp
             :linenos:
 
-            n = graph.num_vertices()
+            auto n = num_vertices(graph);
 
 
 Vertices
@@ -69,6 +72,21 @@ Example:
 
 .. tabs::
 
+    .. tab:: python
+
+        .. code-block:: python
+            :linenos:
+
+            import higra as hg
+
+            g = hg.UndirectedGraph()
+            g.add_vertex()
+            g.add_vertex()
+
+            g.add_vertices(3)
+
+            nv = g.num_vertices() # 5
+
     .. tab:: c++
 
         .. code-block:: cpp
@@ -86,22 +104,6 @@ Example:
             auto nv = num_vertices(g); // 5
 
 
-    .. tab:: python
-
-        .. code-block:: python
-            :linenos:
-
-            import higra as hg
-
-            g = hg.UndirectedGraph()
-            g.add_vertex()
-            g.add_vertex()
-
-            g.add_vertices(3)
-
-            nv = g.num_vertices() # 5
-
-
 Iterating on vertices
 *********************
 
@@ -124,6 +126,20 @@ Iterating on vertices
 
 .. tabs::
 
+    .. tab:: python
+
+        .. code-block:: python
+            :linenos:
+
+            g = hg.UndirectedGraph()
+            ...
+
+            for v in g.vertices():
+                ... # all vertices of g
+
+            for v in g.adjacent_vertices(1):
+                ... # all vertices adjacent to vertex 1 in g
+
     .. tab:: c++
 
         .. code-block:: cpp
@@ -141,20 +157,6 @@ Iterating on vertices
             }
 
 
-    .. tab:: python
-
-        .. code-block:: python
-            :linenos:
-
-            g = hg.UndirectedGraph()
-            ...
-
-            for v in g.vertices():
-                ... # all vertices of g
-
-            for v in g.adjacent_vertices(1):
-                ... # all vertices adjacent to vertex 1 in g
-
 Edges
 -----
 
@@ -216,33 +218,6 @@ Example:
 
 .. tabs::
 
-    .. tab:: c++
-
-        .. code-block:: cpp
-            :linenos:
-
-            #include "higra/graph.hpp"
-            using namespace hg;
-
-            // create a graph with 4 vertices and no edge
-            ugraph g(4);
-
-            // add an edge, between vertex 0 and 1
-            add_edge(0, 1, g);
-            // add an edge, between vertex 0 and 1
-            auto e = add_edge(1, 2, g);
-
-            auto s = source(e, g); // 1
-            auto t = target(e, g); // 2
-            auto ei = index(e, g); // 1
-
-            // add the two edges (3, 0) and (3, 1)
-            add_edges({3, 3}, {0, 1}, g);
-
-            auto ne = num_edges(g); // 4
-
-            auto edges = edge_list(g); // edges.first = {0, 1, 0, 1}, edges.second = {1, 2, 3, 3}
-
     .. tab:: python
 
         .. code-block:: python
@@ -269,6 +244,32 @@ Example:
 
             sources, targets = g.edge_list() # sources = [0, 1, 0, 1], targets = [1, 2, 3, 3]
 
+    .. tab:: c++
+
+        .. code-block:: cpp
+            :linenos:
+
+            #include "higra/graph.hpp"
+            using namespace hg;
+
+            // create a graph with 4 vertices and no edge
+            ugraph g(4);
+
+            // add an edge, between vertex 0 and 1
+            add_edge(0, 1, g);
+            // add an edge, between vertex 0 and 1
+            auto e = add_edge(1, 2, g);
+
+            auto s = source(e, g); // 1
+            auto t = target(e, g); // 2
+            auto ei = index(e, g); // 1
+
+            // add the two edges (3, 0) and (3, 1)
+            add_edges({3, 3}, {0, 1}, g);
+
+            auto ne = num_edges(g); // 4
+
+            auto edges = edge_list(g); // edges.first = {0, 1, 0, 1}, edges.second = {1, 2, 3, 3}
 
 Iterating on edges
 ******************
@@ -297,6 +298,23 @@ Iterating on edges
 
 .. tabs::
 
+    .. tab:: python
+
+        .. code-block:: python
+            :linenos:
+
+            g = hg.UndirectedGraph()
+            ...
+
+            for e in g.edges():
+                print(g.source(e), g.target(e))
+
+            for e in g.in_edges(1):
+                ... # all edges e such that g.target(e) == 1
+
+            for e in g.out_edges(1):
+                ... # all edges e such that g.source(e) == 1
+
     .. tab:: c++
 
         .. code-block:: cpp
@@ -317,25 +335,6 @@ Iterating on edges
                 ... // all edges e such that source(e, g) == 1
             }
 
-
-    .. tab:: python
-
-        .. code-block:: python
-            :linenos:
-
-            g = hg.UndirectedGraph()
-            ...
-
-            for e in g.edges():
-                print(g.source(e), g.target(e))
-
-            for e in g.in_edges(1):
-                ... # all edges e such that g.target(e) == 1
-
-            for e in g.out_edges(1):
-                ... # all edges e such that g.source(e) == 1
-
-
 Degrees
 -------
 
@@ -366,24 +365,6 @@ Operations are done in constant time in ``ugraph``, ``tree``. Operations are don
 
 .. tabs::
 
-    .. tab:: c++
-
-        .. code-block:: cpp
-            :linenos:
-
-            ugraph g;
-            ...
-
-            // degree of vertex 1
-            auto d1 = degree(1, g);
-
-            // in degree of vertex 2
-            auto d2 = in_degree(2, g);
-
-            // out degree of vertex 3
-            auto d3 = out_degree(3, g);
-
-
     .. tab:: python
 
         .. code-block:: python
@@ -401,6 +382,23 @@ Operations are done in constant time in ``ugraph``, ``tree``. Operations are don
             # out degree of vertex 3
             d3 = g.out_degree(3)
 
+    .. tab:: c++
+
+        .. code-block:: cpp
+            :linenos:
+
+            ugraph g;
+            ...
+
+            // degree of vertex 1
+            auto d1 = degree(1, g);
+
+            // in degree of vertex 2
+            auto d2 = in_degree(2, g);
+
+            // out degree of vertex 3
+            auto d3 = out_degree(3, g);
+
 
 Weighted graph
 --------------
@@ -412,6 +410,17 @@ arrays in python.
 
 .. tabs::
 
+    .. tab:: python
+
+        .. code-block:: python
+            :linenos:
+
+            def sum_adjacent_vertices_weights(graph, vertex_weights, vertex):
+                result = 0
+                for v in g.adjacent_vertices(vertex);
+                    result += vertex_weights[v]
+                return result
+
     .. tab:: c++
 
         .. code-block:: cpp
@@ -427,15 +436,3 @@ arrays in python.
                 }
                 return result
             }
-
-
-    .. tab:: python
-
-        .. code-block:: python
-            :linenos:
-
-            def sum_adjacent_vertices_weights(graph, vertex_weights, vertex):
-                result = 0
-                for v in g.adjacent_vertices(vertex);
-                    result += vertex_weights[v]
-                return result
diff --git a/doc/source/tree.rst b/doc/source/tree.rst
index f521d48d..1ea31d46 100644
--- a/doc/source/tree.rst
+++ b/doc/source/tree.rst
@@ -42,19 +42,6 @@ Example:
 
 .. tabs::
 
-    .. tab:: c++
-
-        .. code-block:: cpp
-            :linenos:
-
-            #include "higra/graph.hpp"
-            using namespace hg;
-
-            // creates the tree shown in the figure above
-            tree t({7, 7, 8, 8, 8, 9, 9, 11, 10, 10, 11, 11});
-
-
-
     .. tab:: python
 
         .. code-block:: python
@@ -65,7 +52,16 @@ Example:
             # creates the tree shown in the figure above
             g = hg.Tree((7, 7, 8, 8, 8, 9, 9, 11, 10, 10, 11, 11))
 
+    .. tab:: c++
+
+        .. code-block:: cpp
+            :linenos:
+
+            #include "higra/graph.hpp"
+            using namespace hg;
 
+            // creates the tree shown in the figure above
+            tree t({7, 7, 8, 8, 8, 9, 9, 11, 10, 10, 11, 11});
 
 Basic functions
 ---------------
@@ -88,20 +84,35 @@ Basic functions
     *   - ``parents``
         - array of vertices
         - The parent array
-    *   - ``num_children``
-        - positive integer
-        - Number(s) of children of the given node(s)
-    *   - ``child``
-        - vertex
-        - i-th child of the given node(s)
     *   - ``is_leaf``
         - boolean
         - True if given node(s) is a leaf, False otherwise
 
+
+
 Example:
 
 .. tabs::
 
+    .. tab:: python
+
+        .. code-block:: python
+            :linenos:
+
+            # creates the tree shown in the figure above
+            t = hg.Tree((7, 7, 8, 8, 8, 9, 9, 11, 10, 10, 11, 11))
+
+            t.num_leaves()              # 7
+            t.root()                    # 11
+
+            t.parent(2)                 # 8
+            t.parent((0, 11, 8))       # array {7, 11, 10}
+            t.parents()                 # array {7, 7, 8, 8, 8, 9, 9, 11, 10, 10, 11, 11}
+
+            t.is_leaf(4)                # True
+            t.is_leaf(5)                # False
+            t.is_leaf((0, 11, 8))       # array {True, False, False}
+
     .. tab:: c++
 
         .. code-block:: cpp
@@ -121,42 +132,10 @@ Example:
             parent(vertices, t);        // array {7, 11, 10}
             parents(t);                 // array {7, 7, 8, 8, 8, 9, 9, 11, 10, 10, 11, 11}
 
-            num_children(8, t);         // 3
-            num_children(vertices, t);  // array {0, 2, 3}
-
             is_leaf(4, t);              // true
             is_leaf(5, t);              // false
             is_leaf(vertices, t);       // array {true, false, false}
 
-            child(1, 11, t);            // 10
-            child(0, vertices2, t);     // array {2, 7, 0}
-
-    .. tab:: python
-
-        .. code-block:: python
-            :linenos:
-
-            # creates the tree shown in the figure above
-            t = hg.Tree((7, 7, 8, 8, 8, 9, 9, 11, 10, 10, 11, 11))
-
-            t.num_leaves()              # 7
-            t.root()                    # 11
-
-            t.parent(2)                 # 8
-            t.parent((0, 11, 8))       # array {7, 11, 10}
-            t.parents()                 # array {7, 7, 8, 8, 8, 9, 9, 11, 10, 10, 11, 11}
-
-            t.num_children(8)           # 3
-            t.num_children((0, 11, 8))  # array {0, 2, 3}
-
-
-            t.is_leaf(4)                # True
-            t.is_leaf(5)                # False
-            t.is_leaf((0, 11, 8))       # array {True, False, False}
-
-            t.child(1, 11)              # 10
-            t.child(0, (8, 11, 7))      # array {2, 7, 0}
-
 
 Iterators
 ---------
@@ -167,9 +146,6 @@ Iterators
     *   - Function
         - Returns
         - Description
-    *   - ``children_iterator`` (cpp) ``children`` (python)
-        - a range of iterators (cpp), a list (python)
-        - iterator on the children of the given node
     *   - ``leaves_iterator`` (cpp) ``leaves`` (python)
         - a range of iterators
         - iterator on the leaves of the tree
@@ -187,6 +163,40 @@ Iterators
 
 .. tabs::
 
+    .. tab:: python
+
+        .. code-block:: python
+            :linenos:
+
+            # creates the tree shown in the figure above
+            t = hg.Tree((7, 7, 8, 8, 8, 9, 9, 11, 10, 10, 11, 11))
+
+            for n in t.leaves():
+                ... # 0, 1, 2, ..., 6
+
+            for n in t.ancestors(8):
+                ... # 8, 10, 11
+
+            for n in t.leaves_to_root_iterator(
+                    include_leaves = True, # optional: include (default) or exclude leaves from the iterator
+                    include_root = True): # optional: include (default) or exclude root from the iterator
+                ... // 0, 1, 2, ..., 11
+
+            for n in t.leaves_to_root_iterator(
+                    include_leaves = False,
+                    include_root = False):
+                ... // 7, 8, 9, 10
+
+            for n in t.root_to_leaves_iterator(
+                    include_leaves = True, # optional: include (default) or exclude leaves from the iterator
+                    include_root = True): # optional: include (default) or exclude root from the iterator
+                ... // 11, 10, 9, ..., 0
+
+            for n in t.root_to_leaves_iterator(
+                    include_leaves = False,
+                    include_root = False):
+                ... // 10, 9, 8, 7
+
     .. tab:: c++
 
         .. code-block:: cpp
@@ -195,10 +205,6 @@ Iterators
             // creates the tree shown in the figure above
             tree t({7, 7, 8, 8, 8, 9, 9, 11, 10, 10, 11, 11});
 
-            for(auto n: children_iterator(t, 8)){
-                ... // 2, 3, 4
-            }
-
             for(auto n: leaves_iterator(t)){
                 ... // 0, 1, 2, ..., 6
             }
@@ -232,6 +238,51 @@ Iterators
             }
 
 
+
+
+Children relation
+-----------------
+
+.. important::
+
+    In C++ the children relation is only available on request: one must call the function ``compute_children`` prior
+    to calling any of the following functions (otherwise the behaviour is undefined). Computing the children relation
+    is a linear time operation that will require in the order of :math:`n + 3m` words of memory where :math:`n` is the number
+    of nodes in the tree and :math:`m` is the number of non-leaf nodes (1 word is equal to 64bits on a x64 platform).
+
+    The children relation can be cleared to save space with the function ``clear_children`` and the status of the relation
+    can be checked with the function ``children_computed``.
+
+    In Python the relation is automatically computed when needed, the relation can be cleared with the function ``clear_children``.
+
+.. list-table::
+    :header-rows: 1
+
+    *   - Function
+        - Returns
+        - Description
+    *   - ``num_children``
+        - positive integer
+        - Number(s) of children of the given node(s)
+    *   - ``child``
+        - vertex
+        - i-th child of the given node(s)
+    *   - ``children_iterator`` (cpp) ``children`` (python)
+        - a range of iterators (cpp), a list (python)
+        - iterator on the children of the given node
+    *   - ``compute_children`` (cpp only)
+        -
+        - initialize the children relation (can be called several time safely)
+    *   - ``children_computed`` (cpp only)
+        -
+        - ``true`` if the children relation has already been computed
+    *   - ``clear_children``
+        -
+        - free up the space used to store the children relation
+
+
+.. tabs::
+
     .. tab:: python
 
         .. code-block:: python
@@ -240,34 +291,46 @@ Iterators
             # creates the tree shown in the figure above
             t = hg.Tree((7, 7, 8, 8, 8, 9, 9, 11, 10, 10, 11, 11))
 
+            t.num_children(8)           # 3
+            t.num_children((0, 11, 8))  # array {0, 2, 3}
+
+            t.child(1, 11)              # 10
+            t.child(0, (8, 11, 7))      # array {2, 7, 0}
+
             for n in t.children(8):
                 ... # 2, 3, 4
 
-            for n in t.leaves():
-                ... # 0, 1, 2, ..., 6
 
-            for n in t.ancestors(8):
-                ... # 8, 10, 11
+    .. tab:: c++
 
-            for n in t.leaves_to_root_iterator(
-                    include_leaves = True, # optional: include (default) or exclude leaves from the iterator
-                    include_root = True): # optional: include (default) or exclude root from the iterator
-                ... // 0, 1, 2, ..., 11
+        .. code-block:: cpp
+            :linenos:
 
-            for n in t.leaves_to_root_iterator(
-                    include_leaves = False,
-                    include_root = False):
-                ... // 7, 8, 9, 10
+            // creates the tree shown in the figure above
+            tree t({7, 7, 8, 8, 8, 9, 9, 11, 10, 10, 11, 11});
 
-            for n in t.root_to_leaves_iterator(
-                    include_leaves = True, # optional: include (default) or exclude leaves from the iterator
-                    include_root = True): # optional: include (default) or exclude root from the iterator
-                ... // 11, 10, 9, ..., 0
+            // IMPORTANT: compute the children relation first
+            t.compute_children();
+
+            t.children_computed();      // true
+
+            // two set of vertices
+            array_1d<index_t> vertices{0, 11, 8};
+            array_1d<index_t> vertices2{8, 11, 7};
+
+            num_children(8, t);         // 3
+            num_children(vertices, t);  // array {0, 2, 3}
+
+            child(1, 11, t);            // 10
+            child(0, vertices2, t);     // array {2, 7, 0}
+
+            for(auto n: children_iterator(t, 8)){
+                ... // 2, 3, 4
+            }
+
+            // only if you lack memory and if you are sure that the children relation is not needed anymore
+            t.clear_children();
 
-            for n in t.root_to_leaves_iterator(
-                    include_leaves = False,
-                    include_root = False):
-                ... // 10, 9, 8, 7
 
 Finding nodes
 -------------
@@ -302,6 +365,20 @@ linearithmic time pre-processing.
 
 .. tabs::
 
+    .. tab:: python
+
+        .. code-block:: python
+            :linenos:
+
+                # tree node altitudes
+                altitudes = numpy.asarray((0, 0, 0, 0, 0, 0, 0, 1, 1, 2, 4, 5))
+
+                lca = t.lowest_common_ancestor(2, 5)                    # 10
+                lcas = t.lowest_common_ancestor((2, 8, 0), (5, 6, 11))  # (10, 10, 11)
+
+                // vertices and altitudes
+                auto r = t.find_region((2, 8, 0), (1, 6, 2), altitudes) # (2, 11, 7)
+
     .. tab:: c++
 
         .. code-block:: cpp
@@ -323,19 +400,7 @@ linearithmic time pre-processing.
                 auto r = find_region(vertices, alts, altitudes, t);             // {2, 11, 7}
 
 
-    .. tab:: python
-
-        .. code-block:: python
-            :linenos:
 
-                # tree node altitudes
-                altitudes = numpy.asarray((0, 0, 0, 0, 0, 0, 0, 1, 1, 2, 4, 5))
-
-                lca = t.lowest_common_ancestor(2, 5)                    # 10
-                lcas = t.lowest_common_ancestor((2, 8, 0), (5, 6, 11))  # (10, 10, 11)
-
-                // vertices and altitudes
-                auto r = t.find_region((2, 8, 0), (1, 6, 2), altitudes) # (2, 11, 7)
 
 Accumulators
 ------------
@@ -360,23 +425,23 @@ Accumulators are wrapped into *factories* in C++ while the Python interface only
 
 .. tabs::
 
-    .. tab:: c++
+    .. tab:: python
 
-        .. code-block:: cpp
+        .. code-block:: python
             :linenos:
 
-            auto acc = accumulator_sum();
+            acc = hg.Accumulators.sum
 
-    .. tab:: python
+    .. tab:: c++
 
-        .. code-block:: python
+        .. code-block:: cpp
             :linenos:
 
-            acc = hg.Accumulators.sum
+            auto acc = accumulator_sum();
 
 
 Parallel accumulator
-********************
+~~~~~~~~~~~~~~~~~~~~
 
 The parallel accumulator defines the new value of a node as the accumulation of the values of its children.
 This process is done in parallel on the whole tree.
@@ -405,19 +470,6 @@ The following example demonstrates the application of a parallel sum accumulator
 
 .. tabs::
 
-    .. tab:: c++
-
-        .. code-block:: cpp
-            :linenos:
-
-            // tree in the above example
-            tree t({5, 5, 6, 6, 6, 7, 7, 7});
-            array_1d<index_t> input = xt::ones({num_vertices(t)});
-
-            auto result = accumulate_parallel(t, input, hg::accumulator_sum());
-
-            // result = {0, 0, 0, 0, 0, 2, 3, 2};
-
     .. tab:: python
 
         .. code-block:: python
@@ -431,9 +483,22 @@ The following example demonstrates the application of a parallel sum accumulator
 
             # result = (0, 0, 0, 0, 0, 2, 3, 2)
 
+    .. tab:: c++
+
+        .. code-block:: cpp
+            :linenos:
+
+            // tree in the above example
+            tree t({5, 5, 6, 6, 6, 7, 7, 7});
+            array_1d<index_t> input = xt::ones({num_vertices(t)});
+
+            auto result = accumulate_parallel(t, input, hg::accumulator_sum());
+
+            // result = {0, 0, 0, 0, 0, 2, 3, 2};
+
 
 Sequential accumulator
-**********************
+~~~~~~~~~~~~~~~~~~~~~~
 
 The sequential accumulator defines the new value of a node as the accumulation of the accumulated values of its children.
 This process is thus done sequentially from the leaves to the root of the tree.
@@ -463,19 +528,6 @@ The following example demonstrates the application of a sequential sum accumulat
 
 .. tabs::
 
-    .. tab:: c++
-
-        .. code-block:: cpp
-            :linenos:
-
-            // tree in the above example
-            tree t({5, 5, 6, 6, 6, 7, 7, 7});
-            array_1d<index_t> input = xt::ones({num_leaves(t)});
-
-            auto result = accumulate_sequential(t, input, hg::accumulator_sum());
-
-            // result = {1, 1, 1, 1, 1, 2, 3, 5};
-
     .. tab:: python
 
         .. code-block:: python
@@ -489,9 +541,22 @@ The following example demonstrates the application of a sequential sum accumulat
 
             # result = (1, 1, 1, 1, 1, 2, 3, 5)
 
+    .. tab:: c++
+
+        .. code-block:: cpp
+            :linenos:
+
+            // tree in the above example
+            tree t({5, 5, 6, 6, 6, 7, 7, 7});
+            array_1d<index_t> input = xt::ones({num_leaves(t)});
+
+            auto result = accumulate_sequential(t, input, hg::accumulator_sum());
+
+            // result = {1, 1, 1, 1, 1, 2, 3, 5};
+
 
 Sequential and combine accumulator
-**********************************
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 The sequential and combine accumulator defines the new value of a node as the accumulation of the accumulated values of its children combined with another node dependent value.
 This process is thus done sequentially from the leaves to the root of the tree.
@@ -523,6 +588,20 @@ The following example demonstrates the application of sequential max accumulator
 
 .. tabs::
 
+    .. tab:: python
+
+        .. code-block:: python
+            :linenos:
+
+            # tree in the above example
+            t = hg.Tree((5, 5, 6, 6, 6, 7, 7, 7))
+            leaf_attribute = numpy.ones((t.num_leaves(),))
+            tree_attribute = numpy.ones((t.num_vertices(),))
+
+            result = hg.accumulate_and_add_sequential(t, tree_attribute, leaf_attribute, hg.Accumulators.max)
+
+            # result = (1, 1, 1, 1, 1, 2, 2, 3)
+
     .. tab:: c++
 
         .. code-block:: cpp
@@ -533,7 +612,7 @@ The following example demonstrates the application of sequential max accumulator
             array_1d<index_t> leaf_attribute = xt::ones({num_leaves(t)});
             array_1d<index_t> tree_attribute = xt::ones({num_vertices(t)});
 
-            auto result = accumulate_and_combine_sequential(tree,
+            auto result = accumulate_and_combine_sequential(t,
                                                             tree_attribute,
                                                             leaf_attribute,
                                                             hg::accumulator_max(),
@@ -541,20 +620,6 @@ The following example demonstrates the application of sequential max accumulator
 
             // result = {1, 1, 1, 1, 1, 2, 2, 3};
 
-    .. tab:: python
-
-        .. code-block:: python
-            :linenos:
-
-            # tree in the above example
-            t = hg.Tree((5, 5, 6, 6, 6, 7, 7, 7))
-            leaf_attribute = numpy.ones((t.num_leaves(),))
-            tree_attribute = numpy.ones((t.num_vertices(),))
-
-            result = hg.accumulate_and_add_sequential( tree, tree_attribute, leaf_attribute, hg.Accumulators.max)
-
-            # result = (1, 1, 1, 1, 1, 2, 2, 3)
-
 
 Note that currently, to ease the binding of this accumulator to Python, the combining function cannot be specified at runtime
 and the library offers several statically bound functions:
@@ -573,7 +638,7 @@ They are especially important for writing efficient algorithms in Python by avoi
 Using them in C++ can also be beneficial as they are written to natively and efficiently handle n-dimensional data.
 
 Conditional parallel propagator
-*******************************
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 The conditional parallel propagator defines the new value of a node as its parent value if the condition is true and keeps its value otherwise.
 This process is done in parallel on the whole tree. The default condition (if the user does not provide one) is true for all nodes: each node takes
@@ -604,6 +669,20 @@ The following example demonstrates the application of a conditional parallel pro
 
 .. tabs::
 
+    .. tab:: python
+
+        .. code-block:: python
+            :linenos:
+
+            # tree in the above example
+            t = hg.Tree((5, 5, 6, 6, 6, 7, 7, 7))
+            input = numpy.asarray((1, 2, 3, 4, 5, 6, 7, 8))
+            condition = numpy.asarray((True, False, True, False, True, True, False, False))
+
+            result = hg.propagate_parallel(t, input, condition)
+
+            # result = (6, 2, 7, 4, 7, 8, 7, 8)
+
     .. tab:: c++
 
         .. code-block:: cpp
@@ -618,22 +697,9 @@ The following example demonstrates the application of a conditional parallel pro
 
             // result = {6, 2, 7, 4, 7, 8, 7, 8};
 
-    .. tab:: python
-
-        .. code-block:: python
-            :linenos:
-
-            # tree in the above example
-            t = hg.Tree((5, 5, 6, 6, 6, 7, 7, 7))
-            input = numpy.asarray((1, 2, 3, 4, 5, 6, 7, 8))
-            condition = numpy.asarray((True, False, True, False, True, True, False, False))
-
-            result = hg.propagate_parallel(input, condition, t)
-
-            # result = (6, 2, 7, 4, 7, 8, 7, 8)
 
 Conditional sequential propagator
-*********************************
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 The conditional sequential propagator defines the new value of a node as its parent propagated value if the condition is true and keeps its value otherwise.
 This process is thus done from the root to the leaves of the tree.
@@ -663,6 +729,20 @@ The following example demonstrates the application of a conditional sequential p
 
 .. tabs::
 
+    .. tab:: python
+
+        .. code-block:: python
+            :linenos:
+
+            # tree in the above example
+            t = hg.Tree((5, 5, 6, 6, 6, 7, 7, 7))
+            input = numpy.asarray((1, 2, 3, 4, 5, 6, 7, 8))
+            condition = numpy.asarray((True, False, True, False, True, True, False, False))
+
+            result = hg.propagate_sequential(t, input, condition)
+
+            # result = (8, 2, 7, 4, 7, 8, 7, 8)
+
     .. tab:: c++
 
         .. code-block:: cpp
@@ -677,22 +757,9 @@ The following example demonstrates the application of a conditional sequential p
 
             // result = {8, 2, 7, 4, 7, 8, 7, 8};
 
-    .. tab:: python
-
-        .. code-block:: python
-            :linenos:
-
-            # tree in the above example
-            t = hg.Tree((5, 5, 6, 6, 6, 7, 7, 7))
-            input = numpy.asarray((1, 2, 3, 4, 5, 6, 7, 8))
-            condition = numpy.asarray((True, False, True, False, True, True, False, False))
-
-            result = hg.propagate_sequential(input, condition, t)
-
-            # result = (8, 2, 7, 4, 7, 8, 7, 8)
 
 Sequential propagate and accumulate
-***********************************
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 The sequential propagate and accumulate defines the new value of a node as its parent value accumulated with its current value.
 This process is done from the root to the leaves of the tree.
@@ -721,6 +788,19 @@ The following example demonstrates the application of a propagate and accumulate
 
 .. tabs::
 
+    .. tab:: python
+
+        .. code-block:: python
+            :linenos:
+
+            # tree in the above example
+            t = hg.Tree((5, 5, 6, 6, 6, 7, 7, 7))
+            input = numpy.asarray((1, 2, 3, 4, 5, 6, 7, 8))
+
+            result = hg.propagate_sequential_and_accumulate(t, input, condition)
+
+            # result = (15, 16, 18, 19, 20, 14, 15, 8)
+
     .. tab:: c++
 
         .. code-block:: cpp
@@ -733,16 +813,3 @@ The following example demonstrates the application of a propagate and accumulate
             auto result = propagate_sequential_and_accumulate(t, input, hg.Accumulators.sum);
 
             // result = {15, 16, 18, 19, 20, 14, 15, 8};
-
-    .. tab:: python
-
-        .. code-block:: python
-            :linenos:
-
-            # tree in the above example
-            t = hg.Tree((5, 5, 6, 6, 6, 7, 7, 7))
-            input = numpy.asarray((1, 2, 3, 4, 5, 6, 7, 8))
-
-            result = hg.propagate_sequential_and_accumulate(input, condition, t)
-
-            # result = (15, 16, 18, 19, 20, 14, 15, 8)
diff --git a/higra/interop/py_scipy.cpp b/higra/interop/py_scipy.cpp
index eba5eb5a..70572fe8 100644
--- a/higra/interop/py_scipy.cpp
+++ b/higra/interop/py_scipy.cpp
@@ -30,7 +30,7 @@ auto binary_hierarchy_to_scipy_linkage_matrix(const tree_t &tree,
     auto &area = xarea.derived_cast();
     hg_assert_node_weights(tree, altitudes);
     hg_assert_node_weights(tree, area);
-
+    tree.compute_children();
     hg::index_t n_leaves = num_leaves(tree);
     hg::array_2d<value_t> M = xt::empty<value_t>({(size_t) n_leaves - 1, (size_t) 4});
     for (auto i: hg::leaves_to_root_iterator(tree, hg::leaves_it::exclude)) {
diff --git a/higra/structure/py_tree_graph.cpp b/higra/structure/py_tree_graph.cpp
index 1c634dfe..215a63ad 100644
--- a/higra/structure/py_tree_graph.cpp
+++ b/higra/structure/py_tree_graph.cpp
@@ -17,6 +17,55 @@ using graph_t = hg::tree;
 using edge_t = graph_t::edge_descriptor;
 using vertex_t = graph_t::vertex_descriptor;
 
+
+template<typename graph_t>
+struct def_out_degree_tree {
+    template<typename value_t, typename C>
+    static
+    void def(C &c, const char *doc) {
+        c.def("out_degree", [](graph_t &g,
+                               pyarray<value_t> vertices) {
+                  g.compute_children();
+                  hg_assert_vertex_indices(g, vertices);
+                  return hg::out_degree(vertices, g);
+              },
+              doc,
+              pybind11::arg("vertices_array"));
+    }
+};
+
+template<typename graph_t>
+struct def_degree_tree {
+    template<typename value_t, typename C>
+    static
+    void def(C &c, const char *doc) {
+        c.def("degree", [](graph_t &g,
+                           pyarray<value_t> vertices) {
+                  g.compute_children();
+                  hg_assert_vertex_indices(g, vertices);
+                  return hg::degree(vertices, g);
+              },
+              doc,
+              pybind11::arg("vertices_array"));
+    }
+};
+
+template<typename graph_t>
+struct def_in_degree_tree {
+    template<typename value_t, typename C>
+    static
+    void def(C &c, const char *doc) {
+        c.def("in_degree", [](graph_t &g,
+                              pyarray<value_t> vertices) {
+                  g.compute_children();
+                  hg_assert_vertex_indices(g, vertices);
+                  return hg::in_degree(vertices, g);
+              },
+              doc,
+              pybind11::arg("vertices_array"));
+    }
+};
+
 template<typename graph_t>
 struct def_tree_ctr {
     template<typename type, typename C>
@@ -90,6 +139,7 @@ struct def_num_children {
         c.def("_num_children",
               [](const graph_t &tree, const pyarray<type> &vertices) {
                   hg_assert_vertex_indices(tree, vertices);
+                  tree.compute_children();
                   return hg::num_children(vertices, tree);
               },
               doc,
@@ -105,6 +155,7 @@ struct def_child {
     void def(C &c, const char *doc) {
         c.def("_child",
               [](const graph_t &tree, hg::index_t i, const pyarray<type> &vertices) {
+                  tree.compute_children();
                   hg_assert_vertex_indices(tree, vertices);
                   hg_assert(i >= 0, "Child index cannot be negative.");
                   hg_assert(i < (index_t) (xt::amin)(num_children(vertices, tree))(),
@@ -150,9 +201,97 @@ void py_init_tree_graph(pybind11::module &m) {
     add_type_overloads<def_tree_ctr<graph_t>, HG_TEMPLATE_INTEGRAL_TYPES>
             (c, "Create a tree from the given parent relation.");
     add_edge_accessor_graph_concept<graph_t, decltype(c)>(c);
-    add_incidence_graph_concept<graph_t, decltype(c)>(c);
-    add_bidirectionnal_graph_concept<graph_t, decltype(c)>(c);
-    add_adjacency_graph_concept<graph_t, decltype(c)>(c);
+
+    ///////////////////////////////////
+    // incidence graph concepts
+
+    using iterator_transform_function = std::function<pybind11::tuple(edge_t)>;
+    using out_edge_iterator = hg::transform_forward_iterator<iterator_transform_function,
+            typename hg::graph_traits<graph_t>::out_edge_iterator,
+            pybind11::tuple
+    >;
+
+    c.def("out_edges", [](const graph_t &g,
+                          const vertex_t v) {
+              hg_assert_vertex_index(g, v);
+              g.compute_children();
+              auto it = hg::out_edges(v, g);
+              // wrapping out edge iterator to python friendly type
+              auto it1 = out_edge_iterator(it.first, cpp_edge_2_python<edge_t>);
+              auto it2 = out_edge_iterator(it.second, cpp_edge_2_python<edge_t>);
+              return pybind11::make_iterator(it1, it2);
+
+          },
+          "Iterator over all out edges from 'vertex'. An out edge is a tuple '(vertex, adjacent_vertex)'.",
+          pybind11::arg("vertex"));
+
+    c.def("out_degree", [](graph_t &g, vertex_t vertex) {
+              hg_assert_vertex_index(g, vertex);
+              g.compute_children();
+              return hg::out_degree(vertex, g);
+          },
+          "Return the out degree of the given vertex.",
+          pybind11::arg("vertex"));
+
+    add_type_overloads<def_out_degree_tree<graph_t>, int, unsigned int, long long, unsigned long long>
+            (c, "Return the out degree of the given vertices.");
+
+    ///////////////////////////////////
+    // bidirectionnal graph concepts
+    using in_edge_iterator = hg::transform_forward_iterator<
+            iterator_transform_function,
+            typename hg::graph_traits<graph_t>::in_edge_iterator,
+            pybind11::tuple>;
+
+    c.def("in_edges", [](graph_t &g,
+                         const vertex_t v) {
+              hg_assert_vertex_index(g, v);
+              g.compute_children();
+              auto it = hg::in_edges(v, g);
+              // wrapping in edge iterator to python friendly type
+              auto it1 = in_edge_iterator(it.first, cpp_edge_2_python<edge_t>);
+              auto it2 = in_edge_iterator(it.second, cpp_edge_2_python<edge_t>);
+              return pybind11::make_iterator(it1, it2);
+
+          },
+          "Iterator over all in edges from 'vertex'. An in edge is a tuple '(adjacent_vertex, vertex)'.",
+          pybind11::arg("vertex"));
+
+    c.def("degree", [](graph_t &g, vertex_t vertex) {
+              hg_assert_vertex_index(g, vertex);
+              g.compute_children();
+              return hg::degree(vertex, g);
+          },
+          "Return the degree of the given vertex.",
+          pybind11::arg("vertex"));
+
+    add_type_overloads<def_degree_tree<graph_t>, int, unsigned int, long long, unsigned long long>
+            (c, "Return the degree of the given vertices.");
+
+    c.def("in_degree", [](graph_t &g, vertex_t vertex) {
+              hg_assert_vertex_index(g, vertex);
+              g.compute_children();
+              return hg::in_degree(vertex, g);
+          },
+          "Return the in degree of the given vertex.",
+          pybind11::arg("vertex"));
+
+    add_type_overloads<def_in_degree_tree<graph_t>, int, unsigned int, long long, unsigned long long>
+            (c, "Return the in degree of the given vertices.");
+
+    ///////////////////////////////////
+    // adjacency graph concepts
+
+    c.def("adjacent_vertices", [](const graph_t &g,
+                                  const vertex_t v) {
+              g.compute_children();
+              hg_assert_vertex_index(g, v);
+              auto it = hg::adjacent_vertices(v, g);
+              return pybind11::make_iterator(it.first, it.second);
+          },
+          "Iterator over all vertices adjacent to the given vertex.",
+          pybind11::arg("vertex"));
+
     add_vertex_list_graph_concept<graph_t, decltype(c)>(c);
     add_edge_list_graph_concept<graph_t, decltype(c)>(c);
     add_edge_index_graph_concept<graph_t, decltype(c)>(c);
@@ -172,11 +311,13 @@ void py_init_tree_graph(pybind11::module &m) {
 
     c.def("_num_children", [](const graph_t &t, index_t i) {
         hg_assert_vertex_index(t, i);
+        t.compute_children();
         return t.num_children(i);
     }, "", py::arg("vertex"));
     add_type_overloads<def_num_children<graph_t>, int, unsigned int, long long, unsigned long long>(c, "");
 
     c.def("_child", [](const graph_t &t, index_t i, index_t vertex) {
+              t.compute_children();
               hg_assert_vertex_index(t, vertex);
               hg_assert(i >= 0, "Child index cannot be negative.");
               hg_assert(i < (index_t) t.num_children(vertex),
@@ -214,6 +355,7 @@ void py_init_tree_graph(pybind11::module &m) {
     c.def("children",
           [](const graph_t &g, vertex_t v) {
               hg_assert_vertex_index(g, v);
+              g.compute_children();
               hg::array_1d<hg::index_t> a = hg::array_1d<hg::index_t>::from_shape({hg::num_children(v, g)});
               auto it = hg::children(v, g);
               std::copy(it.first, it.second, a.begin());
@@ -275,5 +417,10 @@ void py_init_tree_graph(pybind11::module &m) {
           py::arg("include_leaves") = true,
           py::arg("include_root") = true);
 
+    c.def("_compute_children", &graph_t::compute_children, "Compute the children relation.");
+    c.def("_children_computed", &graph_t::children_computed,
+          "True if the children relation has already been computed.");
+    c.def("clear_children", &graph_t::clear_children, "Remove the children relation if it has already been computed. "
+                                                      "May free memory but only useful if you are sure that this relation won't be required by further processing).");
 }
 
diff --git a/include/higra/accumulator/accumulator.hpp b/include/higra/accumulator/accumulator.hpp
index e0bbf7ab..16a550b6 100644
--- a/include/higra/accumulator/accumulator.hpp
+++ b/include/higra/accumulator/accumulator.hpp
@@ -39,24 +39,19 @@ namespace hg {
         * @tparam S the storage type
         * @tparam vectorial bool: is dimension of storage > 0 (different from scalar)
         */
-        template<typename S, bool vectorial>
+        template<typename S, typename initializer, typename operation, bool vectorial>
         struct acc_marginal_impl {
         };
 
-        template<typename S>
-        struct acc_marginal_impl<S, true> {
+        template<typename S, typename initializer, typename operation>
+        struct acc_marginal_impl<S, initializer, operation, true> {
 
             static const bool is_vectorial = true;
 
-            using self_type = acc_marginal_impl<S, is_vectorial>;
+            using self_type = acc_marginal_impl<S, initializer, operation, is_vectorial>;
             using value_type = typename std::iterator_traits<S>::value_type;
-            using reducer_type = std::function<value_type(value_type, value_type)>;
 
-
-            acc_marginal_impl(const S &storage_begin, const S &storage_end, const reducer_type &reducer,
-                              const value_type &init_value) :
-                    m_init_value(init_value),
-                    m_reducer(reducer),
+            acc_marginal_impl(const S &storage_begin, const S &storage_end) :
                     m_storage_begin(storage_begin),
                     m_storage_end(storage_end) {
             }
@@ -64,7 +59,7 @@ namespace hg {
             template<typename ...Args>
             void
             initialize(Args &&...) {
-                std::fill(m_storage_begin, m_storage_end, m_init_value);
+                std::fill(m_storage_begin, m_storage_end, initializer::template init_value<value_type>());
             }
 
             template<typename T, typename ...Args>
@@ -72,7 +67,7 @@ namespace hg {
             accumulate(T value_begin, Args &&...) {
                 auto s = m_storage_begin;
                 for (; s != m_storage_end; s++, value_begin++) {
-                    *s = m_reducer(*value_begin, *s);
+                    *s = operation::template reduce<value_type>(*value_begin, *s);
                 }
             };
 
@@ -93,38 +88,32 @@ namespace hg {
             }
 
         private:
-            value_type m_init_value;
-            reducer_type m_reducer;
             S m_storage_begin;
             S m_storage_end;
         };
 
-        template<typename S>
-        struct acc_marginal_impl<S, false> {
+        template<typename S, typename initializer, typename operation>
+        struct acc_marginal_impl<S, initializer, operation, false> {
 
             static const bool is_vectorial = false;
 
-            using self_type = acc_marginal_impl<S, is_vectorial>;
+            using self_type = acc_marginal_impl<S, initializer, operation, is_vectorial>;
             using value_type = typename std::iterator_traits<S>::value_type;
-            using reducer_type = std::function<value_type(value_type, value_type)>;
 
-            acc_marginal_impl(const S &storage_begin, const S &, const reducer_type &reducer,
-                              const value_type &init_value) :
-                    m_init_value(init_value),
-                    m_reducer(reducer),
+            acc_marginal_impl(const S &storage_begin, const S &) :
                     m_storage_begin(storage_begin) {
             }
 
             template<typename ...Args>
             void
             initialize(Args &&...) {
-                *m_storage_begin = m_init_value;
+                *m_storage_begin = initializer::template init_value<value_type>();
             }
 
             template<typename T, typename ...Args>
             void
             accumulate(const T value_begin, Args &&...) {
-                *m_storage_begin = m_reducer(*value_begin, *m_storage_begin);
+                *m_storage_begin = operation::template reduce<value_type>(*value_begin, *m_storage_begin);
             };
 
             template<typename ...Args>
@@ -142,8 +131,6 @@ namespace hg {
             }
 
         private:
-            value_type m_init_value;
-            reducer_type m_reducer;
             S m_storage_begin;
         };
 
@@ -561,11 +548,9 @@ namespace hg {
         auto make_accumulator(S &storage) const {
             using value_type = typename S::value_type;
             using iterator_type = decltype(storage.begin());
-            return accumulator_detail::acc_marginal_impl<iterator_type, vectorial>(
+            return accumulator_detail::acc_marginal_impl<iterator_type, accumulator_sum, accumulator_sum, vectorial>(
                     storage.begin(),
-                    storage.end(),
-                    std::plus<value_type>(),
-                    0);
+                    storage.end());
         }
 
         template<typename shape_t>
@@ -573,6 +558,16 @@ namespace hg {
         auto get_output_shape(const shape_t &input_shape) {
             return input_shape;
         }
+
+        template <typename value_type>
+        constexpr static value_type init_value(){
+            return 0;
+        }
+
+        template <typename value_type>
+        static value_type reduce(const value_type & v1, const value_type & v2){
+            return v1 + v2;
+        }
     };
 
     struct accumulator_min {
@@ -581,11 +576,9 @@ namespace hg {
         auto make_accumulator(S &storage) const {
             using value_type = typename S::value_type;
             using iterator_type = decltype(storage.begin());
-            return accumulator_detail::acc_marginal_impl<iterator_type, vectorial>(
+            return accumulator_detail::acc_marginal_impl<iterator_type, accumulator_min, accumulator_min, vectorial>(
                     storage.begin(),
-                    storage.end(),
-                    [](value_type v1, value_type v2) { return (std::min)(v1, v2); },
-                    (std::numeric_limits<value_type>::max)());
+                    storage.end());
         }
 
         template<typename shape_t>
@@ -593,6 +586,16 @@ namespace hg {
         auto get_output_shape(const shape_t &input_shape) {
             return input_shape;
         }
+
+        template <typename value_type>
+        static value_type init_value(){
+            return (std::numeric_limits<value_type>::max)();
+        }
+
+        template <typename value_type>
+        static value_type reduce(const value_type & v1, const value_type & v2){
+            return (std::min)(v1, v2);
+        }
     };
 
     struct accumulator_max {
@@ -601,11 +604,9 @@ namespace hg {
         auto make_accumulator(S &storage) const {
             using value_type = typename S::value_type;
             using iterator_type = decltype(storage.begin());
-            return accumulator_detail::acc_marginal_impl<iterator_type, vectorial>(
+            return accumulator_detail::acc_marginal_impl<iterator_type, accumulator_max, accumulator_max, vectorial>(
                     storage.begin(),
-                    storage.end(),
-                    [](value_type v1, value_type v2) { return (std::max)(v1, v2); },
-                    std::numeric_limits<value_type>::lowest());
+                    storage.end());
         }
 
         template<typename shape_t>
@@ -613,6 +614,16 @@ namespace hg {
         auto get_output_shape(const shape_t &input_shape) {
             return input_shape;
         }
+
+        template <typename value_type>
+        static value_type init_value(){
+            return std::numeric_limits<value_type>::lowest();
+        }
+
+        template <typename value_type>
+        static value_type reduce(const value_type & v1, const value_type & v2){
+            return (std::max)(v1, v2);
+        }
     };
 
     struct accumulator_prod {
@@ -621,11 +632,9 @@ namespace hg {
         auto make_accumulator(S &storage) const {
             using value_type = typename S::value_type;
             using iterator_type = decltype(storage.begin());
-            return accumulator_detail::acc_marginal_impl<iterator_type, vectorial>(
+            return accumulator_detail::acc_marginal_impl<iterator_type, accumulator_prod, accumulator_prod, vectorial>(
                     storage.begin(),
-                    storage.end(),
-                    std::multiplies<value_type>(),
-                    1);
+                    storage.end());
         }
 
         template<typename shape_t>
@@ -633,6 +642,16 @@ namespace hg {
         auto get_output_shape(const shape_t &input_shape) {
             return input_shape;
         }
+
+        template <typename value_type>
+        static value_type init_value(){
+            return 1;
+        }
+
+        template <typename value_type>
+        static value_type reduce(const value_type & v1, const value_type & v2){
+            return v1 * v2;
+        }
     };
 
     struct accumulator_mean {
diff --git a/include/higra/accumulator/tree_accumulator.hpp b/include/higra/accumulator/tree_accumulator.hpp
index 7eeed089..ebcbc08c 100644
--- a/include/higra/accumulator/tree_accumulator.hpp
+++ b/include/higra/accumulator/tree_accumulator.hpp
@@ -39,6 +39,7 @@ namespace hg {
 
             auto input_view = make_light_axis_view<vectorial>(input);
             auto output_view = make_light_axis_view<vectorial>(output);
+
             auto acc = accumulator.template make_accumulator<vectorial>(output_view);
 
             for (auto i: leaves_iterator(tree)) {
@@ -48,15 +49,40 @@ namespace hg {
                 acc.finalize();
             }
 
-            for (auto i : leaves_to_root_iterator(tree, leaves_it::exclude)) {
-                output_view.set_position(i);
-                acc.set_storage(output_view);
-                acc.initialize();
-                for (auto c : children_iterator(i, tree)) {
-                    input_view.set_position(c);
-                    acc.accumulate(input_view.begin());
+            if (tree.children_computed()) {
+
+                for (auto i : leaves_to_root_iterator(tree, leaves_it::exclude)) {
+                    output_view.set_position(i);
+                    acc.set_storage(output_view);
+                    acc.initialize();
+                    for (auto c : children_iterator(i, tree)) {
+                        input_view.set_position(c);
+                        acc.accumulate(input_view.begin());
+                    }
+                    acc.finalize();
                 }
-                acc.finalize();
+
+            } else {
+                index_t numl = num_leaves(tree);
+                std::vector<decltype(accumulator.template make_accumulator<vectorial>(output_view))> accs;
+                accs.reserve(num_vertices(tree) - numl);
+
+                for (auto i : leaves_to_root_iterator(tree, leaves_it::exclude)) {
+                    output_view.set_position(i);
+                    accs.push_back(accumulator.template make_accumulator<vectorial>(output_view));
+                    accs.back().initialize();
+                }
+
+                for (auto i : leaves_to_root_iterator(tree, leaves_it::include, root_it::exclude)) {
+                    if (i >= numl) {
+                        accs[i - numl].finalize();
+                    }
+
+                    auto p = parent(i, tree);
+                    input_view.set_position(i);
+                    accs[p - numl].accumulate(input_view.begin());
+                }
+                accs.back().finalize();
             }
 
             return output;
@@ -83,7 +109,6 @@ namespace hg {
             auto vertex_data_view = make_light_axis_view<vectorial>(vertex_data);
             auto input_view = make_light_axis_view<vectorial>(output);
             auto output_view = make_light_axis_view<vectorial>(output);
-            auto acc = accumulator.template make_accumulator<vectorial>(output_view);
 
             for (auto i: leaves_iterator(tree)) {
                 output_view.set_position(i);
@@ -91,16 +116,42 @@ namespace hg {
                 output_view = vertex_data_view;
             }
 
-            for (auto i : leaves_to_root_iterator(tree, leaves_it::exclude)) {
-                output_view.set_position(i);
-                acc.set_storage(output_view);
-                acc.initialize();
-                for (auto c : children_iterator(i, tree)) {
-                    input_view.set_position(c);
-                    acc.accumulate(input_view.begin());
+            if (tree.children_computed()) {
+                auto acc = accumulator.template make_accumulator<vectorial>(output_view);
+
+                for (auto i : leaves_to_root_iterator(tree, leaves_it::exclude)) {
+                    output_view.set_position(i);
+                    acc.set_storage(output_view);
+                    acc.initialize();
+                    for (auto c : children_iterator(i, tree)) {
+                        input_view.set_position(c);
+                        acc.accumulate(input_view.begin());
+                    }
+                    acc.finalize();
                 }
-                acc.finalize();
+            } else {
+                index_t numl = num_leaves(tree);
+                std::vector<decltype(accumulator.template make_accumulator<vectorial>(output_view))> accs;
+                accs.reserve(num_vertices(tree) - numl);
+
+                for (auto i : leaves_to_root_iterator(tree, leaves_it::exclude)) {
+                    output_view.set_position(i);
+                    accs.push_back(accumulator.template make_accumulator<vectorial>(output_view));
+                    accs.back().initialize();
+                }
+
+                for (auto i : leaves_to_root_iterator(tree, leaves_it::include, root_it::exclude)) {
+                    if (i >= numl) {
+                        accs[i - numl].finalize();
+                    }
+
+                    auto p = parent(i, tree);
+                    input_view.set_position(i);
+                    accs[p - numl].accumulate(input_view.begin());
+                }
+                accs.back().finalize();
             }
+
             return output;
         };
 
@@ -140,7 +191,6 @@ namespace hg {
             auto input_view = make_light_axis_view<vectorial>(input);
             auto inout_view = make_light_axis_view<vectorial>(output);
             auto output_view = make_light_axis_view<vectorial>(output);
-            auto acc = accumulator.template make_accumulator<vectorial>(output_view);
 
             auto vertex_data_view = make_light_axis_view<vectorial>(vertex_data);
 
@@ -150,17 +200,48 @@ namespace hg {
                 output_view = vertex_data_view;
             }
 
-            for (auto i : leaves_to_root_iterator(tree, leaves_it::exclude)) {
-                output_view.set_position(i);
-                acc.set_storage(output_view);
-                acc.initialize();
-                for (auto c : children_iterator(i, tree)) {
+            if (tree.children_computed()) {
+                auto acc = accumulator.template make_accumulator<vectorial>(output_view);
+
+                for (auto i : leaves_to_root_iterator(tree, leaves_it::exclude)) {
+                    output_view.set_position(i);
+                    acc.set_storage(output_view);
+                    acc.initialize();
+                    for (auto c : children_iterator(i, tree)) {
+                        inout_view.set_position(c);
+                        acc.accumulate(inout_view.begin());
+                    }
+                    acc.finalize();
+                    input_view.set_position(i);
+                    output_view.combine(input_view, combine);
+                }
+            } else {
+                index_t numl = num_leaves(tree);
+                std::vector<decltype(accumulator.template make_accumulator<vectorial>(output_view))> accs;
+                accs.reserve(num_vertices(tree) - numl);
+
+                for (auto i : leaves_to_root_iterator(tree, leaves_it::exclude)) {
+                    output_view.set_position(i);
+                    accs.push_back(accumulator.template make_accumulator<vectorial>(output_view));
+                    accs.back().initialize();
+                }
 
-                    inout_view.set_position(c);
-                    acc.accumulate(inout_view.begin());
+                for (auto i : leaves_to_root_iterator(tree, leaves_it::include, root_it::exclude)) {
+                    if (i >= numl) {
+                        accs[i - numl].finalize();
+                        input_view.set_position(i);
+                        output_view.set_position(i);
+                        output_view.combine(input_view, combine);
+                    }
+
+                    auto p = parent(i, tree);
+                    inout_view.set_position(i);
+                    accs[p - numl].accumulate(inout_view.begin());
                 }
-                acc.finalize();
-                input_view.set_position(i);
+
+                accs.back().finalize();
+                input_view.set_position(root(tree));
+                output_view.set_position(root(tree));
                 output_view.combine(input_view, combine);
             }
 
@@ -404,8 +485,8 @@ namespace hg {
 
     template<typename tree_t, typename T, typename accumulator_t>
     auto propagate_sequential_and_accumulate(const tree_t &tree,
-                              const xt::xexpression<T> &xinput,
-                              const accumulator_t &accumulator) {
+                                             const xt::xexpression<T> &xinput,
+                                             const accumulator_t &accumulator) {
         auto &input = xinput.derived_cast();
 
         if (input.dimension() == 1) {
diff --git a/include/higra/algo/horizontal_cuts.hpp b/include/higra/algo/horizontal_cuts.hpp
index 92d9983f..9f43b10b 100644
--- a/include/higra/algo/horizontal_cuts.hpp
+++ b/include/higra/algo/horizontal_cuts.hpp
@@ -119,7 +119,7 @@ namespace hg {
             auto num_regions = m_num_regions_cuts[cut_index];
             array_1d<index_t> nodes = array_1d<index_t>::from_shape({(size_t) num_regions});
             const tree &ct = (m_use_node_map) ? m_sorted_tree : m_original_tree;
-
+            ct.compute_children();
             if (cut_index == 0) { // special case for single region partition
                 nodes(0) = root(ct);
             } else {
@@ -176,6 +176,7 @@ namespace hg {
         template<typename T, typename E>
         void init(const T &t, const E &a) {
             auto min_alt_children = accumulate_parallel(t, a, accumulator_min());
+            t.compute_children();
 
             // single region partition... edge case
             m_num_regions_cuts.push_back(1);
diff --git a/include/higra/algo/tree.hpp b/include/higra/algo/tree.hpp
index ff4cc38d..53a0c8fa 100644
--- a/include/higra/algo/tree.hpp
+++ b/include/higra/algo/tree.hpp
@@ -168,6 +168,7 @@ namespace hg {
         // number of nodes smaller than a supervertex
         index_t removed = 0;
 
+        tree.compute_children();
         stackv<index_t> stack;
         for (auto n: leaves_iterator(tree)) {
             auto e = n;
@@ -306,6 +307,7 @@ namespace hg {
 
         array_1d<char> attr({num_vertices(tree)}, 0);
 
+        tree.compute_children();
         for (auto i: leaves_iterator(tree)) {
             if (static_cast<bool>(background_marker(i))) {
                 attr(i) = 1;
diff --git a/include/higra/algo/tree_energy_optimization.hpp b/include/higra/algo/tree_energy_optimization.hpp
index ce821216..7a5ae965 100644
--- a/include/higra/algo/tree_energy_optimization.hpp
+++ b/include/higra/algo/tree_energy_optimization.hpp
@@ -538,6 +538,7 @@ namespace hg {
         hg_assert_node_weights(tree, energy_attribute);
         hg_assert_1d_array(energy_attribute);
 
+        tree.compute_children();
         array_1d<bool> optimal_nodes = array_1d<bool>::from_shape({num_vertices(tree)});
         array_1d<value_type> optimal_energy = array_1d<value_type>::from_shape({num_vertices(tree)});
 
@@ -630,6 +631,7 @@ namespace hg {
         using lep_t = hg::tree_energy_optimization_internal::piecewise_linear_energy_function_piece<double>;
         using lef_t = hg::tree_energy_optimization_internal::piecewise_linear_energy_function<double>;
 
+        tree.compute_children();
         std::vector<lef_t> optimal_energies{};
         array_1d<double> apparition_scales = array_1d<double>::from_shape({num_vertices(tree)});
 
diff --git a/include/higra/assessment/fragmentation_curve.hpp b/include/higra/assessment/fragmentation_curve.hpp
index be23d0a3..46fa5163 100644
--- a/include/higra/assessment/fragmentation_curve.hpp
+++ b/include/higra/assessment/fragmentation_curve.hpp
@@ -152,7 +152,9 @@ namespace hg {
             hg_assert_1d_array(ground_truth);
             hg_assert_integral_value_type(ground_truth);
 
-            max_regions = (std::min)(max_regions, num_leaves(tree));
+            m_tree.compute_children();
+
+            max_regions = (std::min)(max_regions, num_leaves(m_tree));
 
             size_t num_regions_ground_truth = xt::amax(ground_truth)() + 1;
 
@@ -171,24 +173,24 @@ namespace hg {
 #endif
 
             // for a tree node i, a gt region j: card_intersection(i, j) is the number of pixels in R_i cap R_j
-            array_2d<index_t> card_intersection_leaves{{num_leaves(tree), region_gt_areas.size()}, 0};
+            array_2d<index_t> card_intersection_leaves{{num_leaves(m_tree), region_gt_areas.size()}, 0};
             if (vertex_map.size() <= 1) { // no rag
-                hg_assert_leaf_weights(tree, ground_truth);
-                for (auto i: leaves_iterator(tree)) {
+                hg_assert_leaf_weights(m_tree, ground_truth);
+                for (auto i: leaves_iterator(m_tree)) {
                     card_intersection_leaves(i, ground_truth(i))++;
                 }
-                region_tree_area = attribute_area(tree);
+                region_tree_area = attribute_area(m_tree);
             } else { // tree on rag
                 hg_assert(vertex_map.size() == ground_truth.size(), "Vertex map and ground truth sizes do not match.");
                 for (index_t i = 0; i < (index_t) vertex_map.size(); i++) {
                     card_intersection_leaves(vertex_map(i), ground_truth(i))++;
                 }
-                region_tree_area = attribute_area(tree,
+                region_tree_area = attribute_area(m_tree,
                                                   rag_accumulate(vertex_map, xt::ones<index_t>(vertex_map.shape()),
                                                                  accumulator_counter()));
             }
 
-            array_2d<double> card_intersection = accumulate_sequential(tree, card_intersection_leaves,
+            array_2d<double> card_intersection = accumulate_sequential(m_tree, card_intersection_leaves,
                                                                        accumulator_sum());
 
             array_1d<double> scores;
@@ -212,17 +214,17 @@ namespace hg {
             }
 
             // initialize scoring for single region partitions (the node itself)
-            for (auto i: leaves_to_root_iterator(tree)) {
+            for (auto i: leaves_to_root_iterator(m_tree)) {
                 backtracking.push_back({{1, scores(i), 0, 0}});
             }
 
-            for (auto i: leaves_to_root_iterator(tree, leaves_it::exclude)) {
-                hg_assert(num_children(i, tree) == 2, "Only binary trees are supported.");
+            for (auto i: leaves_to_root_iterator(m_tree, leaves_it::exclude)) {
+                hg_assert(num_children(i, m_tree) == 2, "Only binary trees are supported.");
 
                 auto &backtrack_i = backtracking[i];
 
-                auto c1 = child(0, i, tree);
-                auto c2 = child(1, i, tree);
+                auto c1 = child(0, i, m_tree);
+                auto c2 = child(1, i, m_tree);
                 auto &backtrack_c1 = backtracking[c1];
                 auto &backtrack_c2 = backtracking[c2];
 
diff --git a/include/higra/attribute/tree_attribute.hpp b/include/higra/attribute/tree_attribute.hpp
index d2992103..7afb19c4 100644
--- a/include/higra/attribute/tree_attribute.hpp
+++ b/include/higra/attribute/tree_attribute.hpp
@@ -72,13 +72,13 @@ namespace hg {
         hg_assert_1d_array(node_area);
         hg_assert_node_weights(tree, node_altitude);
         hg_assert_1d_array(node_altitude);
-
+        tree.compute_children();
         auto &parent = tree.parents();
         array_1d<double> volume = xt::empty<double>({tree.num_vertices()});
         xt::view(volume, xt::range(0, num_leaves(tree))) = 0;
         for (auto i: leaves_to_root_iterator(tree, leaves_it::exclude)) {
             volume(i) = std::fabs(node_altitude(i) - node_altitude(parent(i))) * node_area(i);
-            for (auto c: tree.children(i)) {
+            for (auto c: children_iterator(i, tree)) {
                 volume(i) += volume(c);
             }
         }
@@ -129,6 +129,7 @@ namespace hg {
         hg_assert_1d_array(altitudes);
         using value_type = typename T::value_type;
 
+        tree.compute_children();
         if (increasing_altitudes) {
             auto min_depth = xt::empty_like(altitudes);
             xt::noalias(xt::view(min_depth, xt::range(0, num_leaves(tree)))) =
@@ -195,6 +196,7 @@ namespace hg {
         hg_assert_node_weights(tree, altitudes);
         hg_assert_1d_array(altitudes);
 
+        tree.compute_children();
         array_1d<bool> extrema = xt::zeros<bool>({num_vertices(tree)});
         for (auto n: leaves_to_root_iterator(tree, leaves_it::exclude)) {
             bool flag = true;
@@ -257,6 +259,8 @@ namespace hg {
         hg_assert_node_weights(tree, attribute);
         hg_assert_1d_array(attribute);
 
+        tree.compute_children();
+
         // identify path to the deepest extrema
         array_1d<index_t> ref_son({num_vertices(tree)}, invalid_index);
         if (increasing_altitudes) {
@@ -380,6 +384,7 @@ namespace hg {
     template<typename tree_t>
     auto attribute_sibling(const tree_t &tree, index_t skip = 1) {
 
+        tree.compute_children();
         array_1d<index_t> attribute = xt::empty<index_t>({num_vertices(tree)});
         for (auto n: leaves_to_root_iterator(tree, leaves_it::exclude, root_it::include)) {
             index_t nchs = num_children(n, tree);
@@ -430,6 +435,7 @@ namespace hg {
         array_1d<double> res = array_1d<double>::from_shape({num_vertices(tree)});
         xt::noalias(xt::view(res, xt::range(0, num_leaves(tree)))) = vertex_perimeter;
         array_1d<bool> visited({num_leaves(tree)}, false);
+        tree.compute_children();
 
         for (auto i: leaves_to_root_iterator(tree, leaves_it::exclude)) {
             res(i) = 0;
@@ -462,6 +468,7 @@ namespace hg {
      */
     template<typename tree_t>
     auto attribute_child_number(const tree_t &tree) {
+        tree.compute_children();
         array_1d<index_t> res = xt::empty<index_t>({num_vertices(tree)});
         for (auto i: leaves_to_root_iterator(tree, leaves_it::exclude)) {
             for (index_t c = 0; c < (index_t) num_children(i, tree); c++) {
@@ -526,6 +533,7 @@ namespace hg {
         auto &node_weights = xnode_weights.derived_cast();
         hg_assert_node_weights(tree, node_weights);
 
+        tree.compute_children();
         array_nd<value_type> res = xt::zeros<value_type>(node_weights.shape());
         const auto num_v = num_vertices(tree);
         const auto num_l = num_leaves(tree);
diff --git a/include/higra/hierarchy/hierarchy_core.hpp b/include/higra/hierarchy/hierarchy_core.hpp
index 0af10074..be558df1 100644
--- a/include/higra/hierarchy/hierarchy_core.hpp
+++ b/include/higra/hierarchy/hierarchy_core.hpp
@@ -159,10 +159,13 @@ namespace hg {
     auto simplify_tree(const tree &t, const criterion_t &criterion, bool process_leaves = false) {
         HG_TRACE();
 
+
         // this case is significantly harder because a reordering of the nodes
         // may append if an internal node becomes a leaf
         if (process_leaves) {
 
+            t.compute_children();
+
             // ********************************
             // identification of deleted sub-trees
             // true if all nodes below a given  node are deleted
@@ -256,48 +259,43 @@ namespace hg {
             return make_remapped_tree(tree(new_parent, t.category()), std::move(node_map));
         } else {
             auto n_nodes = num_vertices(t);
-            auto copy_parent = parents(t);
-
-            std::size_t count = 0;
-            array_1d<tree::vertex_descriptor> deleted_map = xt::zeros<tree::vertex_descriptor>(copy_parent.shape());
-            array_1d<bool> deleted = xt::zeros<bool>(copy_parent.shape());
-
-            // from root to leaves, compute the new parent relation,
-            // don't care of the  holes in the parent tab
-            for (auto i: root_to_leaves_iterator(t, leaves_it::exclude, root_it::exclude)) {
-                auto parent = copy_parent(i);
-                if (criterion(i)) {
-                    for (auto c: children_iterator(i, t)) {
-                        copy_parent(c) = parent;
-                    }
-                    count++;
+            auto n_leaves = num_leaves(t);
+
+            array_1d<index_t> new_ranks = xt::empty<index_t>({n_nodes});
+
+            xt::view(new_ranks, xt::range(0, n_leaves)) = xt::arange<index_t>(n_leaves);
+            index_t count = n_leaves;
+
+            for (auto i: leaves_to_root_iterator(t, leaves_it::exclude, root_it::exclude)){
+                if (!criterion(i)) {
+                    new_ranks(i) = count;
+                    ++count;
                 }
-                // number of deleted nodes after node i
-                deleted_map(i) = count;
             }
 
-            //correct the mapping
-            deleted_map = count - deleted_map;
+            new_ranks(root(t)) = count;
+            ++count;
 
-            array_1d<tree::vertex_descriptor> new_parent = xt::arange<tree::vertex_descriptor>(0, n_nodes - count);
-            array_1d<tree::vertex_descriptor> node_map = xt::zeros<tree::vertex_descriptor>({n_nodes - count});
+            array_1d<index_t> new_parent = xt::empty<index_t>({count});
+            array_1d<index_t> node_map = xt::empty<index_t>({count});
 
-            count = 0;
+            count = new_parent.size() - 1;
 
-            for (auto i: leaves_to_root_iterator(t, leaves_it::include, root_it::exclude)) {
+            new_parent(count) = count;
+            node_map(count) = root(t);
+            --count;
+
+            for (auto i: root_to_leaves_iterator(t, leaves_it::include, root_it::exclude)) {
                 if (!criterion(i) || is_leaf(i, t)) {
-                    auto par = copy_parent(i);
-                    auto new_par = par - deleted_map(par);
+                    new_parent(count) = new_ranks(parent(i, t));
                     node_map(count) = i;
-                    new_parent(count) = new_par;
-                    count++;
+                    --count;
+                } else {
+                    new_ranks(i) = new_ranks(parent(i, t));
                 }
             }
 
-            node_map(node_map.size() - 1) = root(t);
-
-
-            return make_remapped_tree(tree(new_parent, t.category()), std::move(node_map));
+            return make_remapped_tree(tree(std::move(new_parent), t.category()), std::move(node_map));
         }
 
     };
@@ -389,6 +387,8 @@ namespace hg {
         auto num_v = num_vertices(tree);
         auto num_l = num_leaves(tree);
         auto num_v_res = num_l * 2 - 1;
+
+        tree.compute_children();
         array_1d<index_t> node_map = array_1d<index_t>::from_shape({num_v});
         array_1d<index_t> reverse_node_map = array_1d<index_t>::from_shape({num_v_res});
         for (index_t i = 0; i < (index_t) num_l; i++) {
diff --git a/include/higra/hierarchy/watershed_hierarchy.hpp b/include/higra/hierarchy/watershed_hierarchy.hpp
index 67e240a3..482a0a70 100644
--- a/include/higra/hierarchy/watershed_hierarchy.hpp
+++ b/include/higra/hierarchy/watershed_hierarchy.hpp
@@ -25,6 +25,7 @@ namespace hg {
                                    const T1 &altitude,
                                    const T2 &attribute) {
             using value_type = typename T2::value_type;
+            tree.compute_children();
             array_1d<value_type> result = xt::empty_like(attribute);
             for (auto n: leaves_iterator(tree)) {
                 result(n) = 0;
diff --git a/include/higra/structure/lca_fast.hpp b/include/higra/structure/lca_fast.hpp
index eacad568..319ca3f8 100644
--- a/include/higra/structure/lca_fast.hpp
+++ b/include/higra/structure/lca_fast.hpp
@@ -89,6 +89,7 @@ namespace hg {
             }
 
             void LCApreprocessEuler(const tree_t &tree) {
+                tree.compute_children();
                 struct se {
                     index_t node;
                     bool first_visit;
diff --git a/include/higra/structure/tree_graph.hpp b/include/higra/structure/tree_graph.hpp
index f3c3b0b9..d20e8d78 100644
--- a/include/higra/structure/tree_graph.hpp
+++ b/include/higra/structure/tree_graph.hpp
@@ -112,31 +112,22 @@ namespace hg {
             tree(const xt::xexpression<T> &parents = xt::xarray<vertex_descriptor>({0}),
                  tree_category category = tree_category::partition_tree) :
                     _parents(parents),
-                    _children(_parents.size()),
+                    _children_computed(false),
+                    _children(0),
                     _category(category) {
                 HG_TRACE();
+                _init();
+            };
 
-                hg_assert(_parents.shape().size() == 1, "parents must be a linear (1d) array");
-                _num_vertices = _parents.size();
-                _root = _num_vertices - 1;
-                hg_assert(_parents(_root) == _root, "nodes are not in a topological order (last node is not a root)");
-
-                for (vertex_descriptor v = 0; v < _root; ++v) {
-                    vertex_descriptor parent_v = _parents(v);
-                    hg_assert(parent_v != v, "several root nodes detected");
-                    hg_assert(parent_v > v, "nodes are not in a topological order");
-                    _children[parent_v].push_back(v);
-                }
-
-                index_t num_leaves = 0;
-
-                for (vertex_descriptor v = 0; v <= _root; ++v) {
-                    if (_children[v].size() == 0) {
-                        hg_assert(num_leaves == v, "leaves nodes are not before internal nodes");
-                        num_leaves++;
-                    }
-                }
-                _num_leaves = (size_t) num_leaves;
+            template<typename T>
+            tree(xt::xexpression<T> &&parents = xt::xarray<vertex_descriptor>({0}),
+                 tree_category category = tree_category::partition_tree) :
+                    _parents(std::move(parents.derived_cast())),
+                    _children_computed(false),
+                    _children(0),
+                    _category(category) {
+                HG_TRACE();
+                _init();
             };
 
             const auto &category() const {
@@ -155,8 +146,18 @@ namespace hg {
                 return (_num_vertices == 0) ? 0 : _num_vertices - 1;
             }
 
+            const auto &children(vertex_descriptor v) const {
+                if (v < _num_leaves) {
+                    return _empty_children;
+                }
+                return _children[v - _num_leaves];
+            }
+
             size_t num_children(const vertex_descriptor v) const {
-                return _children[v].size();
+                if (v < _num_leaves) {
+                    return 0;
+                }
+                return _children[v - _num_leaves].size();
             }
 
             vertex_descriptor root() const {
@@ -164,23 +165,19 @@ namespace hg {
             }
 
             degree_size_type degree(vertex_descriptor v) const {
-                return _children[v].size() + ((v != _root) ? 1 : 0);
+                return num_children(v) + ((v != _root) ? 1 : 0);
             }
 
             auto children_cbegin(vertex_descriptor v) const {
-                return _children[v].cbegin();
+                return children(v).cbegin();
             }
 
             auto children_cend(vertex_descriptor v) const {
-                return _children[v].cend();
-            }
-
-            auto children(vertex_descriptor v) const {
-                return _children[v];
+                return children(v).cend();
             }
 
             auto child(index_t i, vertex_descriptor v) const {
-                return _children[v][i];
+                return _children[v - _num_leaves][i];
             }
 
             template<typename... Args>
@@ -215,7 +212,7 @@ namespace hg {
             }
 
             auto is_leaf(vertex_descriptor v) const {
-                return (size_t) v < _num_leaves;
+                return v < _num_leaves;
             }
 
             template<typename T>
@@ -226,14 +223,65 @@ namespace hg {
                 return v;
             }
 
+            void compute_children() const {
+                if (!_children_computed) {
+                    _children.resize(_num_vertices - _num_leaves);
+                    for (vertex_descriptor v = 0; v < _root; ++v) {
+                        vertex_descriptor parent_v = _parents(v);
+                        _children[parent_v - _num_leaves].push_back(v);
+                    }
+                    _children_computed = true;
+                }
+            }
+
+            void clear_children() const {
+                _children.clear();
+                _children_computed = false;
+            }
+
+            bool children_computed() const {
+                return _children_computed;
+            }
+
+
         private:
 
+            void _init() {
+                hg_assert(_parents.shape().size() == 1, "parents must be a linear (1d) array");
+                _num_vertices = _parents.size();
+                _root = _num_vertices - 1;
+                hg_assert(_parents(_root) == _root, "nodes are not in a topological order (last node is not a root)");
+
+                array_1d<index_t> num_children = xt::zeros_like(_parents);
+                for (vertex_descriptor v = 0; v < _root; ++v) {
+                    vertex_descriptor parent_v = _parents(v);
+                    hg_assert(parent_v != v, "several root nodes detected");
+                    hg_assert(parent_v > v, "nodes are not in a topological order");
+                    num_children(parent_v)++;
+                }
+
+                index_t num_leaves = 0;
+
+                for (vertex_descriptor v = 0; v <= _root; ++v) {
+                    if (num_children(v) == 0) {
+                        hg_assert(num_leaves == v, "leaves nodes are not before internal nodes");
+                        num_leaves++;
+                    }
+                }
+                _num_leaves = (size_t) num_leaves;
+            }
+
+
             vertex_descriptor _root;
             size_t _num_vertices;
-            size_t _num_leaves;
+            index_t _num_leaves;
             array_1d <vertex_descriptor> _parents;
-            std::vector<children_list_t> _children;
+            mutable bool _children_computed;
+            mutable std::vector<children_list_t> _children;
             tree_category _category;
+
+            // for leaf nodes
+            children_list_t _empty_children;
         };
 
 
@@ -492,7 +540,8 @@ namespace hg {
     inline
     std::pair<tree::children_iterator, tree::children_iterator>
     children(const tree::vertex_descriptor v, const tree &g) {
-        return std::make_pair(g.children_cbegin(v), g.children_cend(v));
+        auto &c = g.children(v);
+        return std::make_pair(c.cbegin(), c.cend());
     }
 
     inline
@@ -575,10 +624,11 @@ namespace hg {
     std::pair<typename hg::tree::adjacency_iterator, typename hg::tree::adjacency_iterator>
     adjacent_vertices(typename hg::tree::vertex_descriptor v, const hg::tree &g) {
         using it = typename hg::tree::adjacency_iterator;
+        auto &c = g.children(v);
         auto par = g.parent(v);
         return std::make_pair(
-                it(v, par, g.children_cbegin(v)),
-                it(par, par, g.children_cend(v)));
+                it(v, par, c.cbegin()),
+                it(par, par, c.cend()));
     }
 
 
@@ -588,12 +638,13 @@ namespace hg {
         auto fun = [v](const typename hg::tree::vertex_descriptor t) {
             return hg::tree::edge_descriptor(v, t, (std::min)(v, t));
         };
+        auto &c = g.children(v);
         using it = typename hg::tree::out_edge_iterator;
         using ita = typename hg::tree::adjacency_iterator;
         auto par = g.parent(v);
         return std::make_pair(
-                it(ita(v, par, g.children_cbegin(v)), fun),
-                it(ita(par, par, g.children_cend(v)), fun));
+                it(ita(v, par, c.cbegin()), fun),
+                it(ita(par, par, c.cend()), fun));
     }
 
     inline
@@ -602,12 +653,13 @@ namespace hg {
         auto fun = [v](const typename hg::tree::vertex_descriptor t) {
             return hg::tree::edge_descriptor(t, v, (std::min)(v, t));
         };
+        auto &c = g.children(v);
         using it = typename hg::tree::out_edge_iterator;
         using ita = typename hg::tree::adjacency_iterator;
         auto par = g.parent(v);
         return std::make_pair(
-                it(ita(v, par, g.children_cbegin(v)), fun),
-                it(ita(par, par, g.children_cend(v)), fun));
+                it(ita(v, par, c.cbegin()), fun),
+                it(ita(par, par, c.cend()), fun));
     }
 
     template<typename T>


Test patch:
diff --git a/test/cpp/accumulator/test_tree_accumulator.cpp b/test/cpp/accumulator/test_tree_accumulator.cpp
index 1ceb520c..73529c16 100644
--- a/test/cpp/accumulator/test_tree_accumulator.cpp
+++ b/test/cpp/accumulator/test_tree_accumulator.cpp
@@ -52,6 +52,65 @@ namespace tree_accumulator {
     TEST_CASE("accumulator tree vectorial", "[tree_accumulator]") {
 
         auto tree = data.t;
+        tree.compute_children();
+
+        REQUIRE(tree.children_computed());
+        array_2d<index_t> input{{1, 0},
+                                {1, 1},
+                                {1, 2},
+                                {1, 3},
+                                {1, 4},
+                                {1, 5},
+                                {1, 6},
+                                {1, 7}};
+
+        auto res1 = accumulate_parallel(tree, input, hg::accumulator_min());
+        auto longmax = std::numeric_limits<index_t>::max();
+        array_2d<index_t> ref1{{longmax, longmax},
+                               {longmax, longmax},
+                               {longmax, longmax},
+                               {longmax, longmax},
+                               {longmax, longmax},
+                               {1,       0},
+                               {1,       2},
+                               {1,       5}};
+        REQUIRE(xt::allclose(ref1, res1));
+
+        array_2d<index_t> vertex_data{{1, 0},
+                                      {1, 1},
+                                      {1, 2},
+                                      {1, 3},
+                                      {1, 4}};
+        auto res2 = accumulate_sequential(tree, vertex_data, hg::accumulator_sum());
+        array_2d<index_t> ref2{{1, 0},
+                               {1, 1},
+                               {1, 2},
+                               {1, 3},
+                               {1, 4},
+                               {2, 1},
+                               {3, 9},
+                               {5, 10}};
+        REQUIRE(xt::allclose(ref2, res2));
+
+        auto res3 = accumulate_and_combine_sequential(tree, input, vertex_data, hg::accumulator_sum(),
+                                                      std::plus<index_t>());
+        array_2d<index_t> ref3{{1, 0},
+                               {1, 1},
+                               {1, 2},
+                               {1, 3},
+                               {1, 4},
+                               {3, 6},
+                               {4, 15},
+                               {8, 28}};
+        REQUIRE(xt::allclose(ref3, res3));
+
+    }
+
+    TEST_CASE("accumulator tree vectorial no children", "[tree_accumulator]") {
+
+        auto tree = data.t;
+
+        REQUIRE(!tree.children_computed());
 
         array_2d<index_t> input{{1, 0},
                                 {1, 1},
diff --git a/test/cpp/structure/test_tree.cpp b/test/cpp/structure/test_tree.cpp
index 5560a22f..c38cc395 100644
--- a/test/cpp/structure/test_tree.cpp
+++ b/test/cpp/structure/test_tree.cpp
@@ -22,11 +22,36 @@ namespace tree {
 
         hg::tree t;
 
-        _data() : t(xt::xarray<long>{5, 5, 6, 6, 6, 7, 7, 7}) {
+        _data() : t(array_1d<index_t>{5, 5, 6, 6, 6, 7, 7, 7}) {
+            t.compute_children();
         }
 
     } data;
 
+    TEST_CASE("tree ctr", "[tree]") {
+        hg::tree t1;
+        REQUIRE(hg::num_vertices(t1) == 0);
+
+        array_1d<index_t> p2{5, 5, 6, 6, 6, 7, 7, 7};
+        hg::tree t2(p2);
+        REQUIRE(hg::num_vertices(t2) == 8);
+        REQUIRE(p2.size() == 8);
+
+        array_1d<index_t> p3{5, 5, 6, 6, 6, 7, 7, 7};
+        hg::tree t3(std::move(p3));
+        REQUIRE(hg::num_vertices(t2) == 8);
+        REQUIRE(p3.size() == 0);
+    }
+
+    TEST_CASE("tree children", "[tree]") {
+        hg::tree t(array_1d<index_t>{5, 5, 6, 6, 6, 7, 7, 7});
+        REQUIRE(t.children_computed() == false);
+        t.compute_children();
+        REQUIRE(t.children_computed() == true);
+        t.clear_children();
+        REQUIRE(t.children_computed() == false);
+    }
+
     TEST_CASE("tree sizes", "[tree]") {
         auto t = data.t;
         REQUIRE(hg::category(t) == tree_category::partition_tree);