tpl_graph_utils.H

Go to the documentation of this file.

 
/*
                          Aleph_w
 
  Data structures & Algorithms
  version 2.0.0b
  https://github.com/lrleon/Aleph-w
 
  This file is part of Aleph-w library
 
  Copyright (c) 2002-2026 Leandro Rabindranath Leon
 
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
 
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
 
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/
 
 
# ifndef TPL_GRAPH_UTILS_H
# define TPL_GRAPH_UTILS_H
 
 
# include <cassert>
# include <cstddef>
# include <limits>
# include <memory>
# include <tuple>
# include <utility>
# include <vector>
# include <tpl_agraph.H>
# include <tpl_dynListQueue.H>
# include <ah-errors.H>
 
using namespace Aleph;
 
namespace Aleph {
 
  template <class GT> inline static bool
  __depth_first_traversal(const GT & g, typename GT::Node * node,
                          typename GT::Arc *  arc,
                          bool (*visit)(const GT & g, typename GT::Node *,
                                        typename GT::Arc *),
                          size_t & count);
 
  template <class GT> inline size_t 
  depth_first_traversal(const GT & g, typename GT::Node * start_node,
                        bool (*visit)(const GT & g, typename GT::Node *, 
                                      typename GT::Arc *))
  {
    g.reset_bit_nodes(Depth_First);
    g.reset_bit_arcs(Depth_First);
    size_t counter = 0;
 
    __depth_first_traversal(g, start_node, nullptr, visit, counter);
 
    return counter;
  }
 
  template <class GT> inline 
  size_t depth_first_traversal(const GT & g, 
                               bool (*visit)(const GT &, typename GT::Node *, 
                                             typename GT::Arc *))
  {
    return depth_first_traversal(g, g.get_first_node(), visit);
  }
 
  template <class GT>
  struct Default_Visit_Op
  {
    bool operator () (const GT &, typename GT::Node *, typename GT::Arc *)
    {
      return false;
    }
  };
 
  template <class GT, 
            class Operation = Default_Visit_Op<GT>, 
            class SA        = Dft_Show_Arc<GT>> 
  class Depth_First_Traversal
  {
    Operation * op_ptr = nullptr;
    SA          sa;
    size_t      count = 0;
    const GT *  g_ptr = nullptr;
 
  private:
 
    bool __dft(typename GT::Node * node, typename GT::Arc * arc = nullptr)
    {
      if (IS_NODE_VISITED(node, Depth_First)) 
        return false; 
 
      NODE_BITS(node).set_bit(Depth_First, true);
      count++;
 
      if ((*op_ptr)(*g_ptr, node, arc)) 
        return true;
 
      if (count == g_ptr->get_num_nodes())
        return true;
 
      // Recursively traverse arcs incident to `node`.
      for (Node_Arc_Iterator<GT, SA> it(node, sa); it.has_curr(); it.next_ne())
        {
          auto arc = it.get_current_arc_ne();
          if (IS_ARC_VISITED(arc, Depth_First)) 
            continue;
 
          ARC_BITS(arc).set_bit(Depth_First, true);
          if (__dft (it.get_tgt_node_ne(), arc))
            return true;
        }
 
      return false;
    }
 
    size_t dft(const GT & g, typename GT::Node * start_node, Operation & __op)
    {
      op_ptr = &__op;
      g_ptr  = &g;
    
      g_ptr->reset_bit_nodes(Depth_First);
      g_ptr->reset_bit_arcs(Depth_First);
 
      count = 0;
 
      __dft(start_node);
 
      return count;
    }
 
  public:
 
    Depth_First_Traversal(SA __sa = SA()) : sa(__sa) { /* empty */ }
 
    size_t operator () (const GT & g, Operation op = Operation()) 
    {
      return dft(g, g.get_first_node(), op);
    }
 
    size_t operator () (const GT & g, typename GT::Node * sn,
                        Operation op = Operation())
    {
      return dft(g, sn, op);
    }
  };
 
 
  template <class GT> inline static bool
  __depth_first_traversal(const GT & g, typename GT::Node * node, 
                          typename GT::Arc *  arc,
                          bool (*visit)(const GT & g, typename GT::Node *, 
                                        typename GT::Arc *),
                          size_t & count)
  {
    if (IS_NODE_VISITED(node, Depth_First)) 
      return false; 
 
    NODE_BITS(node).set_bit(Depth_First, true); // mark node visited
    count++;
 
    if (visit != nullptr) // invoke callback if provided
      if ((*visit)(g, node, arc)) 
        return true;
 
    if (count == g.get_num_nodes()) // all nodes discovered?
      return true;
 
    for (auto it = g.get_arc_it(node); it.has_curr(); it.next_ne())
      {
        auto arc = it.get_current_arc_ne();
        if (IS_ARC_VISITED(arc, Depth_First)) 
          continue;
 
        ARC_BITS(arc).set_bit(Depth_First, true); // mark arc visited
        if (__depth_first_traversal(g, it.get_tgt_node_ne(), arc, visit, count))
          return true;
      }
 
    return false; // keep exploring
  } 
 
 
  template <class GT> inline size_t
  breadth_first_traversal(const GT & g, typename GT::Node * start,
                          bool (*visit)(const GT &, typename GT::Node *,
                                        typename GT::Arc *) )
  {
    g.reset_bit_nodes(Breadth_First);
    g.reset_bit_arcs(Breadth_First);
    DynListQueue<typename GT::Arc*> q;
 
    for (auto it = g.get_arc_it(start); it.has_curr(); it.next_ne())
      q.put(it.get_current_arc_ne());
 
    NODE_BITS(start).set_bit(Breadth_First, true);
    size_t node_counter = 1;
 
    if (visit != nullptr)
      if ((*visit)(g, start, nullptr))
        return 1;
 
    while (not q.is_empty() and node_counter < g.get_num_nodes())
      {
        auto arc = q.get();
        ARC_BITS(arc).set_bit(Breadth_First, true);
 
        auto src = g.get_src_node(arc);
        auto tgt = g.get_tgt_node(arc);
        if (IS_NODE_VISITED(src, Breadth_First) and
            IS_NODE_VISITED(tgt, Breadth_First))
          continue;
 
        auto visit_node = IS_NODE_VISITED(src, Breadth_First) ? tgt : src;
        if (visit != nullptr)
          if ((*visit)(g, visit_node, arc))
            break;
 
        NODE_BITS(visit_node).set_bit(Breadth_First, true); 
        node_counter++;
 
        for (auto it = g.get_arc_it(visit_node); it.has_curr(); it.next_ne())
          {
            auto curr_arc = it.get_current_arc_ne();
            if (IS_ARC_VISITED(curr_arc, Breadth_First)) 
              continue; 
 
            if (IS_NODE_VISITED(g.get_src_node(curr_arc), Breadth_First) and 
                IS_NODE_VISITED(g.get_tgt_node(curr_arc), Breadth_First))
              continue; // both endpoints already visited
    
            q.put(curr_arc);
          }
      }
 
    return node_counter;
  }
 
  template <class GT> inline size_t
  breadth_first_traversal(GT & g, 
                          bool (*visit)(const GT &, typename GT::Node *, 
                                        typename GT::Arc *))
  {
    return breadth_first_traversal(g, g.get_first_node(), visit);
  }
 
      
  template <class GT, 
            class Operation = Default_Visit_Op<GT>, 
            class SA        = Dft_Show_Arc<GT>> 
  class Breadth_First_Traversal
  {
    SA     sa;
    size_t count = 0;
 
    size_t bft(const GT & g, typename GT::Node * start, Operation & op)
    {
      g.reset_bit_nodes(Breadth_First); 
      g.reset_bit_arcs(Breadth_First);  
      DynListQueue<typename GT::Arc*> q; // pending arcs queue
 
      for (Node_Arc_Iterator<GT, SA> it(start, sa); it.has_curr(); it.next_ne())
        q.put(it.get_current_arc_ne());
 
      NODE_BITS(start).set_bit(Breadth_First, true); 
      count = 1;
 
      if (op (g, start, nullptr))
        return 1;
 
      while (not q.is_empty() and count < g.get_num_nodes()) 
        {
          auto arc = q.get();
          ARC_BITS(arc).set_bit(Breadth_First, true); 
 
          auto src = g.get_src_node(arc); 
          auto tgt = g.get_tgt_node(arc);
 
          if (IS_NODE_VISITED(src, Breadth_First) and
              IS_NODE_VISITED(tgt, Breadth_First))
            continue;
 
          auto curr = IS_NODE_VISITED(src, Breadth_First) ? tgt : src;
          if (op (g, curr, arc))
            break;
 
          NODE_BITS(curr).set_bit(Breadth_First, true); 
          count++;
 
          for (Node_Arc_Iterator<GT, SA> it(curr, sa); it.has_curr(); it.next_ne())
            {
              auto curr_arc = it.get_current_arc_ne();
              if (IS_ARC_VISITED(curr_arc, Breadth_First)) 
                continue; 
 
              if (IS_NODE_VISITED(g.get_src_node(curr_arc), Breadth_First) and 
                  IS_NODE_VISITED(g.get_tgt_node(curr_arc), Breadth_First))
                continue;
 
              q.put(curr_arc);
            }
        }
 
      return count;
    }  
 
  public:
 
    Breadth_First_Traversal(SA __sa = SA()) : sa(__sa) { /* empty */ }
 
    size_t operator () (const GT & g, Operation op)
    {
      return bft (g, g.get_first_node(), op);
    }
 
    size_t operator () (const GT & g, typename GT::Node * p, 
                        Operation && op = Operation())
    {
      return bft(g, p, op);
    }
 
    size_t operator () (const GT & g, typename GT::Node * p, Operation & op)
    {
      return bft(g, p, op);
    }
  };
 
  template <class GT> inline 
  Path<GT> find_path_breadth_first(const GT & g, typename GT::Node * start, 
                                   typename GT::Node * end)
  {
    ah_invalid_argument_if(start == nullptr or end == nullptr)
      << "find_path_breadth_first(): start and end must be non-null";
 
    if (start == end)
      return Path<GT>(g, start);
 
    g.reset_nodes();
    g.reset_arcs(); 
  
    DynListQueue<typename GT::Arc*> q; 
 
    for (auto it = g.get_arc_it(start); it.has_curr(); it.next_ne())
      q.put(it.get_current_arc_ne());
        
    NODE_BITS(start).set_bit(Find_Path, true);
 
    bool path_found = false;
 
    while (not q.is_empty())
      {
        auto arc = q.get(); 
        auto src = g.get_src_node(arc); 
        auto tgt = g.get_tgt_node(arc);
 
        if (IS_NODE_VISITED(src, Find_Path) and IS_NODE_VISITED(tgt, Find_Path))
          continue;
      
        if (IS_NODE_VISITED(tgt, Find_Path))
          std::swap(src, tgt);
 
        ARC_BITS(arc).set_bit(Find_Path, true); 
        NODE_BITS(tgt).set_bit(Find_Path, true);
        NODE_COOKIE(tgt) = src;
        
        if (tgt == end)
          {
            path_found = true;
            break;
          }
 
        for (auto it = g.get_arc_it(tgt); it.has_curr(); it.next_ne())
          {
            auto a = it.get_current_arc_ne();
            if (IS_ARC_VISITED(a, Find_Path))
              continue;
 
            if (IS_NODE_VISITED(g.get_src_node(a), Find_Path) and 
                IS_NODE_VISITED(g.get_tgt_node(a), Find_Path))
              continue;
 
            q.put(a);
          }
      }
  
    if (not path_found)
      return Path<GT>(g);
 
    q.empty(); // free queue memory for eventually saving the required for path
  
    Path<GT> path(g, end);
    auto p = end;
    while (p != start)
      {
        p = (typename GT::Node *) NODE_COOKIE(p);
        path.insert(p);
      }
 
    return path;
  }
 
 
  template <class GT> inline 
  bool test_connectivity(const GT & g)
  {
    ah_domain_error_if(g.is_digraph())
      << "test_connectivity() does not work on digraphs";
 
    const auto num_nodes = g.get_num_nodes();
    if (num_nodes == 0)
      return false;
 
    if (g.get_num_arcs() + 1 < num_nodes)
      return false; 
 
    return depth_first_traversal<GT>(g, nullptr) == num_nodes;
  }
 
 
  template <class GT> inline static
  bool __test_cycle(const GT & g, typename GT::Node *, typename GT::Node *);
 
 
  template <class GT> inline 
  bool test_for_cycle(const GT & g, typename GT::Node * src)
  {
    ah_invalid_argument_if(src == nullptr)
      << "test_for_cycle(): src must be non-null";
 
    g.reset_bit_nodes(Test_Cycle);
    g.reset_bit_arcs(Test_Cycle);
    for (auto it = g.get_arc_it(src); it.has_curr(); it.next_ne())
      {
        auto arc = it.get_curr();
        if (IS_ARC_VISITED(arc, Test_Cycle)) 
          continue;
 
        ARC_BITS(arc).set_bit(Test_Cycle, true);
        if (__test_cycle(g, src, it.get_tgt_node_ne()))
          return true;
      }
 
    return false;
  }
 
 
  template <class GT> inline static bool
  __test_cycle(const GT & g, typename GT::Node * src, typename GT::Node * curr)
  {
    if (src == curr) 
      return true; // detected cycle
 
    if (IS_NODE_VISITED(curr, Test_Cycle)) 
      return false;
 
    NODE_BITS(curr).set_bit(Test_Cycle, true);
 
    for (auto it = g.get_arc_it(curr); it.has_curr(); it.next_ne())
      {
        auto arc = it.get_curr();
        if (IS_ARC_VISITED(arc, Test_Cycle)) 
          continue;
 
        ARC_BITS(arc).set_bit(Test_Cycle, true);
        if (__test_cycle(g, src, it.get_tgt_node_ne()))
          return true;
      }
 
    return false; 
  }  
 
  template <class GT> inline static
  bool __is_graph_acyclique(const GT & g, typename GT::Node * curr_node)
  {
    if (IS_NODE_VISITED(curr_node, Test_Cycle)) 
      return false;
 
    NODE_BITS(curr_node).set_bit(Test_Cycle, true);
 
    for (auto it = g.get_arc_it(curr_node); it.has_curr(); it.next_ne())
      {
        auto arc = it.get_current_arc_ne();
        if (IS_ARC_VISITED(arc, Test_Cycle)) 
          continue; 
 
        ARC_BITS(arc).set_bit(Test_Cycle, true); 
 
        if (not __is_graph_acyclique(g, it.get_tgt_node_ne())) 
          return false;
      }
    // All outgoing arcs explored without detecting a back edge.
    return true; 
  }
 
 
  template <class GT> inline 
  bool is_graph_acyclique(const GT & g, typename GT::Node * start_node)
  {
    ah_domain_error_if(g.is_digraph())
      << "is_graph_acyclique() does not work for digraphs";
 
    const auto num_nodes = g.get_num_nodes();
    if (num_nodes == 0)
      return true;
 
    ah_invalid_argument_if(start_node == nullptr)
      << "is_graph_acyclique(): start_node must be non-null";
 
    if (num_nodes == 1)
      return true;
  
    if (g.get_num_arcs() >= num_nodes) 
      return false; 
 
    g.reset_bit_arcs(Test_Cycle);
    g.reset_bit_nodes(Test_Cycle);
 
    return __is_graph_acyclique(g, start_node);
  }
 
  template <class GT> inline 
  bool is_graph_acyclique(const GT & g)
  {
    ah_domain_error_if(g.is_digraph())
      << "is_graph_acyclique() does not work for digraphs";
 
    const auto num_nodes = g.get_num_nodes();
    if (num_nodes == 0)
      return true;
 
    if (num_nodes == 1)
      return true;
 
    if (g.get_num_arcs() >= num_nodes) 
      return false; 
 
    g.reset_bit_arcs(Test_Cycle);
    g.reset_bit_nodes(Test_Cycle);
 
    for (auto it = g.get_node_it(); it.has_curr(); it.next_ne()) 
      {
        auto current_node = it.get_current_node_ne();
        if (IS_NODE_VISITED(current_node, Test_Cycle)) 
          continue; 
 
        if (not __is_graph_acyclique(g, current_node)) 
          return false;
      }
 
    return true;
  }
 
      
  template <class GT> 
  inline bool has_cycle(const GT & g)
  {
    return not is_graph_acyclique(g);
  }
 
 
  template <class GT> inline 
  bool test_for_path(const GT & g, typename GT::Node * start_node, 
                     typename GT::Node * end_node)
  { 
    ah_invalid_argument_if(start_node == nullptr or end_node == nullptr)
      << "test_for_path(): start_node and end_node must be non-null";
 
    if (start_node == end_node)
      return true;
 
    g.reset_bit_nodes(Find_Path); 
    g.reset_bit_arcs(Find_Path);
 
    for (auto it = g.get_arc_it(start_node); it.has_curr(); it.next_ne())
      {
        auto arc = it.get_current_arc_ne();
        ARC_BITS(arc).set_bit(Find_Path, true);
        if (__test_for_path(g, it.get_tgt_node_ne(), end_node)) 
          return true; 
      }
 
    return false;
  }
 
 
  template <class GT> inline static
  bool __test_for_path(const GT & g, typename GT::Node * curr_node, 
                       typename GT::Node * end_node)  
  {
    if (curr_node == end_node) 
      return true;
 
    if (IS_NODE_VISITED(curr_node, Find_Path))
      return false;
 
    NODE_BITS(curr_node).set_bit(Find_Path, true);
 
    for (auto it = g.get_arc_it(curr_node); it.has_curr(); it.next_ne())
      { 
        auto arc = it.get_current_arc_ne();
        if (IS_ARC_VISITED(arc, Find_Path)) 
          continue; 
 
        ARC_BITS(arc).set_bit(Find_Path, true);
        if (__test_for_path(g, it.get_tgt_node_ne(), end_node)) 
          return true;
      }
 
    return false;
  }
 
  template <class GT> inline 
  DynList<GT> inconnected_components(const GT & g);
 
  template <class GT> inline 
  void build_subgraph(const GT & g, GT & sg, 
                      typename GT::Node * g_src, size_t & node_count); 
    
 
  template <class GT> inline
  DynList<GT> inconnected_components(const GT & g) 
  {
    g.reset_nodes(); 
    g.reset_arcs(); 
 
    DynList<GT> list;
    size_t count = 0; // visited node counter
    for (auto it = g.get_node_it(); count < g.get_num_nodes() and it.has_curr();
         it.next_ne())
      {
        auto curr = it.get_current_node_ne();
        if (IS_NODE_VISITED(curr, Build_Subtree)) 
          continue;
 
        list.append(GT());
        GT & subgraph = list.get_last();
        build_subgraph(g, subgraph, curr, count); 
      }
 
    return list;
  }
 
 
  template <class GT> inline 
  void build_subgraph(const GT & g, GT & sg, 
                      typename GT::Node * g_src, size_t & node_count)
  {
    if (IS_NODE_VISITED(g_src, Build_Subtree)) 
      return;
 
    NODE_BITS(g_src).set_bit(Build_Subtree, true);
    ++node_count; 
 
    auto sg_src = mapped_node<GT>(g_src);
    if (sg_src == nullptr) // not mapped yet
      {
        sg_src = sg.insert_node(g_src->get_info());
        GT::map_nodes(g_src, sg_src);
      }
 
    for (auto it = g.get_arc_it(g_src); 
         node_count < g.get_num_nodes() and it.has_curr(); it.next_ne())
      {
        auto arc = it.get_current_arc_ne();
        if (IS_ARC_VISITED(arc, Build_Subtree)) 
          continue; 
 
        ARC_BITS(arc).set_bit(Build_Subtree, true);
        auto g_tgt  = it.get_tgt_node_ne();
        auto sg_tgt = mapped_node<GT>(g_tgt);
        if (sg_tgt == nullptr) // sg_tgt mapped in sg?
          {
            sg_tgt = sg.insert_node(g_tgt->get_info()); 
            GT::map_nodes(g_tgt, sg_tgt); 
          }
 
        auto sg_arc = sg.insert_arc(sg_src, sg_tgt, arc->get_info());
        GT::map_arcs(arc, sg_arc); 
 
        build_subgraph(g, sg, g_tgt, node_count);
      }
  }
 
  
  template <class GT> inline static
  bool __find_depth_first_spanning_tree(const GT &          g, 
                                        typename GT::Node * gnode, 
                                        typename GT::Arc *  garc, 
                                        GT &                tree,
                                        typename GT::Node * tnode);
 
  template <class GT> inline
  GT find_depth_first_spanning_tree(const GT & g, typename GT::Node * gnode)
  {
    g.reset_nodes();
    g.reset_arcs(); 
 
    GT tree;
 
    NODE_BITS(gnode).set_bit(Spanning_Tree, true);
  
    auto tnode = tree.insert_node(gnode->get_info()); 
    GT::map_nodes(gnode, tnode);
  
    for (auto it = g.get_arc_it(gnode); it.has_curr(); it.next_ne())
      {
        auto arc = it.get_current_arc_ne();
        if (IS_ARC_VISITED(arc, Spanning_Tree)) 
          continue;
 
        auto arc_tgt_node = it.get_tgt_node_ne();
        if (IS_NODE_VISITED(arc_tgt_node, Spanning_Tree))
          continue;
 
        if (__find_depth_first_spanning_tree(g, arc_tgt_node, arc, tree, tnode))
          return tree;
      }
 
    return tree;
  }
 
  template <class GT> inline
  GT find_depth_first_spanning_tree(const GT & g)
  {
    return find_depth_first_spanning_tree(g, g.get_node());
  }
 
  
  template <class GT> inline static
  bool __find_depth_first_spanning_tree(const GT & g, typename GT::Node * gnode, 
                                        typename GT::Arc *  garc, 
                                        GT & tree, typename GT::Node * tnode)
  {
    NODE_BITS(gnode).set_bit(Spanning_Tree, true);
    ARC_BITS(garc).set_bit(Spanning_Tree, true);
 
    auto tree_tgt_node = tree.insert_node(gnode->get_info());
    GT::map_nodes(gnode, tree_tgt_node);
 
    auto tarc = tree.insert_arc(tnode, tree_tgt_node, garc->get_info()); 
    GT::map_arcs(garc, tarc);
 
    tnode = tree_tgt_node; 
    if (tree.get_num_nodes() == g.get_num_nodes())
      return true;
 
    assert(tree.get_num_nodes() > tree.get_num_arcs()); // tree invariant
 
    for (auto it = g.get_arc_it(gnode); it.has_curr(); it.next_ne())
      {
        auto arc = it.get_current_arc_ne();
        if (IS_ARC_VISITED(arc, Spanning_Tree)) 
          continue;
 
        auto arc_tgt_node = it.get_tgt_node_ne();
        if (IS_NODE_VISITED(arc_tgt_node, Spanning_Tree))
          continue;
 
        if (__find_depth_first_spanning_tree(g, arc_tgt_node, arc, tree, tnode))
          return true; // spanning tree completed
      }
 
    return false;
  }
 
  template <class GT> inline
  GT find_breadth_first_spanning_tree(GT & g, typename GT::Node * gp)
  {
    g.reset_bit_nodes(Spanning_Tree);
    g.reset_bit_arcs(Spanning_Tree);
 
    GT tree;
 
    std::unique_ptr<typename GT::Node> tp_auto(new typename GT::Node(gp));
    tree.insert_node(tp_auto.get());
    GT::map_nodes(gp, tp_auto.release());
    NODE_BITS(gp).set_bit(Spanning_Tree, true);
 
    DynListQueue<typename GT::Arc*> q;
    for (auto it = g.get_arc_it(gp); it.has_curr(); it.next_ne())
      q.put(it.get_curr());
 
    while (not q.is_empty()) 
      {
        auto garc = q.get(); 
        ARC_BITS(garc).set_bit(Spanning_Tree, true);
        auto gsrc = g.get_src_node(garc);  
        auto gtgt = g.get_tgt_node(garc);
 
        if (IS_NODE_VISITED(gsrc, Spanning_Tree) and 
            IS_NODE_VISITED(gtgt, Spanning_Tree))
          continue;
 
        if (IS_NODE_VISITED(gtgt, Spanning_Tree)) // gtgt visited?
          std::swap(gsrc, gtgt); // make gsrc the visited endpoint
 
        auto tsrc = mapped_node<GT>(gsrc);
        NODE_BITS(gtgt).set_bit(Spanning_Tree, true);
 
        // Create and map the new tree node.
        std::unique_ptr<typename GT::Node> ttgt_auto(new typename GT::Node(gtgt));
        tree.insert_node(ttgt_auto.get());
        auto ttgt = ttgt_auto.release();
        GT::map_nodes(gtgt, ttgt);
 
        // Insert and map the new tree arc.
        auto tarc = tree.insert_arc(tsrc, ttgt, garc->get_info());
        GT::map_arcs(garc, tarc);
        if (tree.get_num_nodes() == g.get_num_nodes())
          break; 
 
        // Enqueue arcs incident to the newly visited node.
        for (auto it = g.get_arc_it(gtgt); it.has_curr(); it.next_ne())
          {
            auto current_arc = it.get_current_arc_ne();
            if (IS_ARC_VISITED(current_arc, Spanning_Tree)) 
              continue;
 
            if (IS_NODE_VISITED(g.get_src_node(current_arc),Spanning_Tree) and 
                IS_NODE_VISITED(g.get_tgt_node(current_arc),Spanning_Tree))
              continue;
            q.put(current_arc);
          }
      }
 
    return tree;
  }
 
  template <class GT>
  GT build_spanning_tree(const DynArray<typename GT::Arc*> & arcs)
  {
    using Node = typename GT::Node;
    using Arc  = typename GT::Arc;
 
    GT ret;
    DynMapTree<Node*, Node*> table;
    arcs.for_each([&table, &ret] (Arc * ga)
                  {
                    if (ga == nullptr)
                      return;
 
                    Node * gsrc = (Node*) ga->src_node;
                    Node * gtgt = (Node*) ga->tgt_node;
 
                    Node * tsrc;
                    auto * pair_ptr = table.search(gsrc);
                    if (pair_ptr)
                      tsrc = pair_ptr->second;
                    else
                      {
                        tsrc = ret.insert_node(gsrc->get_info());
                        table.insert(gsrc, tsrc);
                        NODE_COOKIE(tsrc) = gsrc;
                      }
  
                    Node * ttgt;
                    pair_ptr = table.search(gtgt);
                    if (pair_ptr)
                      ttgt = pair_ptr->second;
                    else
                      {
                        ttgt = ret.insert_node(gtgt->get_info());
                        table.insert(gtgt, ttgt);
                        NODE_COOKIE(ttgt) = gtgt;
                      }
 
                    Arc * ta = ret.insert_arc(tsrc, ttgt);
                    *ta = *ga;
                    ARC_COOKIE(ta) = ga;
                  });
 
    return ret;
  }
 
  template <class GT> inline static 
  long & df(typename GT::Node * p)
  {
    return NODE_COUNTER(p);
  }
 
  template <class GT> inline static 
  long & low(typename GT::Node * p)
  {
    return reinterpret_cast<long&>(NODE_COOKIE(p));
  }
 
  template <class GT> inline static
  void __compute_cut_nodes(const GT & g, DynList<typename GT::Node *> & list, 
                           typename GT::Node * p, typename GT::Arc * a,
                           long & curr_df);
 
  template <class GT>
  DynList<typename GT::Node*> 
  compute_cut_nodes(const GT & g, typename GT::Node * start)
  {
    DynList<typename GT::Node*> list;
 
    ah_domain_error_if(g.is_digraph())
      << "compute_cut_nodes() does not work on digraphs";
 
    if (g.get_num_nodes() == 0)
      return list;
 
    ah_invalid_argument_if(start == nullptr)
      << "compute_cut_nodes(): start must be non-null";
 
    using Node = typename GT::Node;
 
    std::vector<Node*> nodes;
    nodes.reserve(g.get_num_nodes());
    for (auto it = g.get_node_it(); it.has_curr(); it.next_ne())
      nodes.push_back(it.get_curr());
 
    std::vector<long> low_values(nodes.size(), -1);
 
    for (size_t i = 0; i < nodes.size(); ++i)
      {
        auto p = nodes[i];
        NODE_COUNTER(p) = 0;
        NODE_BITS(p).reset();
        NODE_COOKIE(p) = &low_values[i];
      }
 
    struct Cookie_Guard
    {
      std::vector<Node*> & nodes;
 
      ~Cookie_Guard()
      {
        for (auto p : nodes)
          NODE_COOKIE(p) = nullptr;
      }
    } cookie_guard{nodes};
 
    g.reset_arcs();
    long current_df = 0;
    NODE_BITS(start).set_bit(Depth_First, true);
    low<GT>(start) = df<GT>(start) = current_df++;
    int call_counter = 0;
  
    // Traverse arcs incident to `start`.
    for (auto it = g.get_arc_it(start); 
         it.has_curr() and current_df < g.get_num_nodes(); it.next_ne())
      {
        auto tgt = it.get_tgt_node_ne();
        if (IS_NODE_VISITED(tgt, Depth_First)) 
          continue; 
 
        auto arc = it.get_current_arc_ne();
        if (IS_ARC_VISITED(arc, Depth_First)) 
          continue;
 
        ARC_BITS(arc).set_bit(Depth_First, true);
        __compute_cut_nodes(g, list, tgt, arc, current_df); 
        ++call_counter;
      }
 
    // Root is an articulation point iff it has more than one DFS child.
    if (call_counter > 1)
      {
        NODE_BITS(start).set_bit(Cut, true);
        list.append(start);
      }
 
    return list;
  }
 
  template <class GT>
  DynList<typename GT::Node*> compute_cut_nodes(const GT & g)
  {
    return compute_cut_nodes(g, g.get_node());
  }
      
  template <class GT> inline static
  void __compute_cut_nodes(const GT & g, DynList<typename GT::Node *> & list, 
                           typename GT::Node * p, typename GT::Arc * a, 
                           long & curr_df)
  {
    NODE_BITS(p).set_bit(Depth_First, true);
    low<GT>(p) = df<GT>(p) = curr_df++;
 
    bool p_is_cut_node = false;
    for (auto it = g.get_arc_it(p); it.has_curr(); it.next_ne())
      {
        auto arc = it.get_current_arc_ne();
        if (arc == a) 
          continue; // parent arc
 
        auto tgt = it.get_tgt_node_ne();
        if (IS_NODE_VISITED(tgt, Depth_First)) 
          { 
            if (not IS_ARC_VISITED(arc, Depth_First))
              if (df<GT>(tgt) < low<GT>(p))
                low<GT>(p) = df<GT>(tgt);
 
            continue;
          }
 
        if (IS_ARC_VISITED(arc, Depth_First)) 
          continue;
 
        ARC_BITS(arc).set_bit(Depth_First, true);
 
        __compute_cut_nodes(g, list, tgt, arc, curr_df);
 
        if (low<GT>(tgt) < low<GT>(p)) 
          low<GT>(p) = low<GT>(tgt);
 
        if (low<GT>(tgt) >= df<GT>(p) and df<GT>(tgt) != 0)
          p_is_cut_node = true;
      }
 
    if (p_is_cut_node)
      {
        NODE_BITS(p).set_bit(Cut, true);
        list.append(p);
      }
  }
 
  const long Cross_Arc = -1;
 
  template <class GT> inline static 
  bool is_a_cross_arc(typename GT::Arc * a) 
  {
    return ARC_COUNTER(a) == Cross_Arc; 
  }
 
  template <class GT> inline static 
  bool is_a_cut_node(typename GT::Node * p)
  {
    return NODE_BITS(p).get_bit(Cut);
  }
 
  template <class GT> inline static 
  bool is_an_cut_arc(typename GT::Arc * a)
  {
    return ARC_BITS(a).get_bit(Cut);
  }
 
  template <class GT> inline static 
  bool is_node_painted(typename GT::Node * p)
  {
    return NODE_COUNTER(p) > 0;
  }
 
  template <class GT> inline static 
  bool is_arc_painted(typename GT::Arc * arc)
  {
    return ARC_COUNTER(arc) > 0;
  }
 
  template <class GT> inline static 
  void paint_node(typename GT::Node * p, const long & color)
  {
    NODE_COUNTER(p) = color;
  }
 
  template <class GT> inline static 
  void paint_arc(typename GT::Arc * a, const long & color)
  {
    ARC_COUNTER(a) = color;
  }
 
  template <class GT> inline static 
  const long & get_color(typename GT::Node * p)
  {
    return NODE_COUNTER(p);
  }
 
  template <class GT> inline static 
  const long & get_color(typename GT::Arc * a)
  {
    return ARC_COUNTER(a);
  }
 
  template <class GT> inline static
  void __paint_subgraph(const GT & g, typename GT::Node * p, long current_color)
  {
    assert(not is_a_cut_node <GT> (p));
 
    if (is_node_painted <GT> (p)) 
      return; 
 
    paint_node <GT> (p, current_color);
 
    for (auto it = g.get_arc_it(p); it.has_curr(); it.next_ne())
      {
        auto arc = it.get_current_arc_ne();
        if (is_arc_painted <GT> (arc))
          continue;
 
        auto tgt = it.get_tgt_node_ne();
        if (is_a_cut_node <GT> (tgt))
          continue;
 
        paint_arc <GT> (arc, current_color); 
 
        __paint_subgraph(g, tgt, current_color); 
      }
  }
 
  template <class GT> inline static
  void 
  __paint_from_cut_node(const GT & g, typename GT::Node * p, long & current_color)
  {
    assert(is_a_cut_node <GT> (p));
 
    // Paint each adjacent non-cut block with a new color.
    for (auto it = g.get_arc_it(p); it.has_curr(); it.next_ne())
      {
        auto arc = it.get_current_arc_ne();
 
        assert(not is_arc_painted <GT> (arc));
 
        auto tgt_node = it.get_tgt_node_ne();
        if (is_a_cut_node <GT> (tgt_node)) // cut-to-cut arc
          {
            ARC_BITS(arc).set_bit(Cut, true);
            continue;
          }
        else 
          {
            paint_arc <GT> (arc, Cross_Arc);
            if (is_node_painted <GT> (tgt_node)) 
              continue; 
          }
 
        __paint_subgraph(g, tgt_node, current_color);
 
        ++current_color;
 
        assert(not is_arc_painted <GT> (arc));
      }
  }
 
  template <class GT> inline long 
  paint_subgraphs(const GT & g, const DynList<typename GT::Node*> & cut_node_list)
  {
    g.reset_counter_nodes();
    g.reset_counter_arcs();
    long current_color = 1;
 
    for (auto it = cut_node_list.get_it(); it.has_curr(); it.next_ne())
      __paint_from_cut_node(g, it.get_curr(), current_color);
 
    return current_color;
  }
 
  template <class GT> inline static
  void __map_subgraph(const GT & g, GT & sg, typename GT::Node * gsrc, 
                      const long color)
  {
    assert(get_color <GT> (gsrc) == color);
 
    auto tsrc = mapped_node<GT>(gsrc); // image of gsrc in sg
 
    // Traverse arcs and copy those with the requested color.
    for (auto it = g.get_arc_it(gsrc); it.has_curr(); it.next_ne())
      {
        auto garc = it.get_current_arc_ne();
        if (get_color<GT>(garc) != color or IS_ARC_VISITED(garc, Build_Subtree))
          continue;
 
        ARC_BITS(garc).set_bit(Build_Subtree, true); 
 
        auto gtgt = it.get_tgt_node_ne(); 
 
        assert(get_color <GT> (gtgt) == color);
 
        typename GT::Node * ttgt = nullptr; // image of gtgt in sg
        if (IS_NODE_VISITED(gtgt, Build_Subtree))
          ttgt = mapped_node<GT> (gtgt);
        else
          {     // gtgt not in sg yet: copy & map it
            std::unique_ptr<typename GT::Node> ttgt_auto(new typename GT::Node(gtgt));
            sg.insert_node(ttgt_auto.get());
            GT::map_nodes(gtgt, ttgt_auto.get());
            NODE_BITS(gtgt).set_bit(Build_Subtree, true); 
            ttgt = ttgt_auto.release(); 
          }
 
        auto tarc = sg.insert_arc(tsrc, ttgt, garc->get_info());
        GT::map_arcs(garc, tarc);
 
        __map_subgraph(g, sg, gtgt, color);
      }
  }
 
  template <class GT>
  GT map_subgraph(const GT & g, const long color)
  {
    g.reset_bit_nodes(Build_Subtree);
    g.reset_bit_arcs(Build_Subtree);
 
    typename GT::Node * first = nullptr;
    for (auto it = g.get_node_it(); it.has_curr(); it.next_ne())
      if (get_color <GT> (it.get_current_node_ne()) == color)
        {
          first = it.get_current_node_ne();
          break;
        }
 
    ah_domain_error_if(first == nullptr)
      << "Color does not exist in the graph";
 
    GT sg;
    std::unique_ptr<typename GT::Node> auto_tsrc(new typename GT::Node(first));
    sg.insert_node(auto_tsrc.get());
    GT::map_nodes(first, auto_tsrc.release());
    NODE_BITS(first).set_bit(Build_Subtree, true);
 
    __map_subgraph(g, sg, first, color);
  
    return sg;
  }
 
  template <class GT>
  std::tuple<GT, DynList<typename GT::Arc*>>
  map_cut_graph(const GT & g, const DynList<typename GT::Node*> & cut_node_list)
  {
    GT cut_graph;
    DynList<typename GT::Arc*> cross_arc_list;
 
    for (auto it = cut_node_list.get_it(); it.has_curr(); it.next_ne())
      {
        auto gp = it.get_curr();
 
        assert(is_a_cut_node <GT> (gp));
 
        std::unique_ptr<typename GT::Node> tp_auto(new typename GT::Node(gp));
        cut_graph.insert_node(tp_auto.get());
        GT::map_nodes(gp, tp_auto.release());
      }
 
    // cut_graph contains cut-to-cut arcs; cross_arc_list collects cut-to-noncut arcs.
    for (auto it = g.get_arc_it(); it.has_curr(); it.next_ne())
      {
        auto garc = it.get_current_arc_ne();
        if (is_a_cross_arc <GT> (garc))
          {
            cross_arc_list.append(garc); 
            continue;
          }
 
        if (not is_an_cut_arc <GT> (garc)) 
          continue;
 
        auto src = mapped_node<GT>(g.get_src_node(garc));
        auto tgt = mapped_node<GT>(g.get_tgt_node(garc));
 
        assert(src != nullptr and tgt != nullptr);
 
        auto arc = cut_graph.insert_arc(src, tgt, garc->get_info());
        GT::map_arcs(garc, arc);
      }
  
    return { cut_graph, cross_arc_list };
  }
 
 
  template <class GT, class Distance> 
  struct Distance_Compare
  {
    Distance dist;
 
    Distance_Compare(Distance __dist = Distance()) : dist(__dist) { /* empty */ }
 
    bool operator () (typename GT::Arc * a1, typename GT::Arc * a2) const
    {
      return dist(a1) < dist(a2);
    }
  };
 
 
  template <class GT> 
  GT invert_digraph(const GT & g)
  {
    ah_domain_error_unless(g.is_digraph())
      << "invert_digraph() requires a digraph (or a graph temporarily treated as a digraph)";
 
    g.reset_nodes();
    g.reset_arcs();
    GT gi;
 
    // Copy all nodes first so isolated vertices are preserved.
    for (auto it = g.get_node_it(); it.has_curr(); it.next_ne())
      {
        auto gp = it.get_curr();
        auto ip = gi.insert_node(gp->get_info());
        GT::map_nodes(gp, ip);
      }
 
    for (auto it = g.get_arc_it(); it.has_curr(); it.next_ne())
      {
        auto arc   = it.get_curr();
 
        auto ssrc = g.get_src_node(arc);
        auto stgt = g.get_tgt_node(arc);
 
        auto rsrc = mapped_node<GT>(ssrc);
        auto rtgt = mapped_node<GT>(stgt);
 
        assert(rsrc != nullptr and rtgt != nullptr);
 
        typename GT::Arc * ai = gi.insert_arc(rtgt, rsrc, arc->get_info());
        GT::map_arcs(arc, ai);
      }
 
    assert(g.get_num_arcs() == gi.get_num_arcs() and 
           g.get_num_nodes() == gi.get_num_nodes());
 
    return gi;
  }
 
 
  template <class GT, class SA = Dft_Show_Arc<GT>>
  class Invert_Digraph
  {
    SA sa;
 
  public:
 
    Invert_Digraph(SA __sa) : sa(__sa) { /* empty */ }
 
    GT operator () (const GT & g) const
    {
      ah_domain_error_unless(g.is_digraph())
        << "Invert_Digraph requires a digraph (or a graph temporarily treated as a digraph)";
 
      g.reset_nodes();
      g.reset_arcs();
      GT gi;
 
      // Copy all nodes first so isolated vertices are preserved.
      for (auto it = g.get_node_it(); it.has_curr(); it.next_ne())
        {
          auto gp = it.get_curr();
          auto ip = gi.insert_node(gp->get_info());
          GT::map_nodes(gp, ip);
        }
 
      for (Arc_Iterator<GT, SA> it(g, sa); it.has_curr(); it.next_ne())
        {
          auto arc  = it.get_curr();
 
          auto ssrc = g.get_src_node(arc);
          auto stgt = g.get_tgt_node(arc);
 
          auto rsrc = mapped_node<GT>(ssrc);
          auto rtgt = mapped_node<GT>(stgt);
 
          assert(rsrc != nullptr and rtgt != nullptr);
 
          typename GT::Arc * ai = gi.insert_arc(rtgt, rsrc, arc->get_info());
          GT::map_arcs(arc, ai);
        }
 
      assert(g.get_num_nodes() == gi.get_num_nodes());
 
      return gi;
    }
  };
 
  template <class GT>
  class Dft_Dist
  {
  public:
 
    typedef typename GT::Arc_Type Distance_Type;
 
    static const Distance_Type Zero_Distance;
 
    static const Distance_Type Max_Distance;
 
    Distance_Type & operator () (typename GT::Arc * a) const
    {
      return a->get_info();
    }
 
    Distance_Type & operator () (typename GT::Arc * a, typename GT::Node*) const
    {
      return a->get_info();
    }
 
    static void set_zero(typename GT::Arc * a) { a->get_info() = 0; }
  };
 
  template <class GT>
  const typename Dft_Dist<GT>::Distance_Type Dft_Dist<GT>::Max_Distance =
    std::numeric_limits<typename Dft_Dist<GT>::Distance_Type>::max();
 
  template <class GT>
  const typename Dft_Dist<GT>::Distance_Type Dft_Dist<GT>::Zero_Distance =
    typename Dft_Dist<GT>::Distance_Type{};
  
 
  template <class GT, class SA = Dft_Show_Arc<GT>>
  typename GT::Arc *
  search_spanning_tree_arc(const GT & g,
                           typename GT::Node *src, typename GT::Node *tgt,
                           SA sa = SA()) noexcept
  {
    assert(src != nullptr and tgt != nullptr);
 
    // For efficiency, iterate over the node with fewer arcs
    if (not g.is_digraph() and tgt->num_arcs < src->num_arcs)
      std::swap(tgt, src);
 
    typename GT::Arc * first_match = nullptr;
    for (Node_Arc_Iterator<GT, SA> it(src, sa); it.has_curr(); it.next_ne())
      {
        if (it.get_tgt_node_ne() != tgt)
          continue;
 
        auto arc = it.get_current_arc_ne();
        // Prefer arc marked as spanning tree
        if (IS_ARC_VISITED(arc, Aleph::Spanning_Tree))
          return arc;
 
        // Keep first match as fallback
        if (first_match == nullptr)
          first_match = arc;
      }
 
    return first_match; // Return first match if no spanning tree arc found
  }
 
  template <class GT, class Distance = Dft_Dist<GT>>
  typename Distance::Distance_Type
  get_min_path(typename GT::Node * s, typename GT::Node * end, Path<GT> & path)
  {
    using Distance_Type = typename Distance::Distance_Type;
 
    ah_invalid_argument_if(s == nullptr or end == nullptr)
      << "get_min_path(): s and end must be non-null";
 
    Distance_Type dist{};
    path.empty();
 
    if (s == end)
      {
        path.init(end);
        return dist;
      }
 
    // Build path from end to start following cookies, collecting arcs
    Distance distance;
    DynList<std::pair<typename GT::Node*, typename GT::Arc*>> node_arc_list;
 
    auto curr = end;
    while (curr != s)
      {
        auto prev = static_cast<typename GT::Node *>(NODE_COOKIE(curr));
        ah_domain_error_if(prev == nullptr)
          << "get_min_path(): broken cookie chain (nullptr)";
 
        // Find the spanning tree arc between prev and curr
        auto arc = search_spanning_tree_arc<GT>(path.get_graph(), prev, curr);
        ah_domain_error_if(arc == nullptr)
          << "get_min_path(): no arc connecting nodes in path";
 
        node_arc_list.insert(std::make_pair(curr, arc));
        dist += distance(arc);
        curr = prev;
      }
 
    // Now build path from start to end
    path.init(s);
    for (auto it = node_arc_list.get_it(); it.has_curr(); it.next())
      {
        auto & [node, arc] = it.get_curr();
        (void)node; // suppress warning
        path.append(arc);
      }
 
    return dist;
  }
 
  template <class GT, 
            class Distance = Dft_Dist<GT>, 
            class SA       = Dft_Show_Arc<GT>>
  class Total_Cost
  {
    Distance dist;
    SA       sa;
    typename Distance::Distance_Type sum{};
 
  public:
 
    Total_Cost(Distance __dist = Distance(), SA __sa = SA()) 
      : dist(__dist), sa(__sa)
    {
      // empty
    }
 
    typename Distance::Distance_Type total_cost(GT & g)
    {
      sum = typename Distance::Distance_Type{};
 
      // Traverse all arcs and sum their weights
      for (Arc_Iterator <GT, SA> it(g, sa); it.has_curr(); it.next_ne())
        sum += dist(it.get_current_arc_ne());
 
      return sum;
    }
 
    typename Distance::Distance_Type operator () (GT & g)
    {
      return total_cost (g);
    }
 
    void reset() noexcept
    {
      sum = typename Distance::Distance_Type{};
    }
 
    typename Distance::Distance_Type value() const noexcept
    {
      return sum;
    }
 
    bool operator () (typename GT::Arc * a)
    {
      if (not sa(a))
        return true; // skip but keep traversing
 
      sum += dist(a);
      return true; // keep traversing
    }
  };
 
 
 
} // end namespace Aleph
 
# endif // TPL_GRAPH_UTILS_H