Changeset 2025:93fcadf94ab0 in lemon-0.x for lemon/min_cost_arborescence.h

Timestamp:

03/31/06 13:10:44 (18 years ago)

Author:

Balazs Dezso

Branch:

default

Phase:

public

Convert:

svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@2664

Message:

Bugfix in the minimum cost arborescence algorithm
Dual solution computation and interface for algorithm
Optimality test on random graph

File:

• lemon/min_cost_arborescence.h (modified) (28 diffs)

lemon/min_cost_arborescence.h

@@ -27 +27 @@
 
 #include <lemon/list_graph.h>
+#include <lemon/bin_heap.h>
 
 namespace lemon {
@@ -57 +58 @@
     ///
     /// The type of the map that stores which edges are in the arborescence.
-    /// It must meet the \ref concept::ReadWriteMap "ReadWriteMap" concept.
-    /// Initially it will be setted to false on each edge. The algorithm
-    /// may set each value one time to true and maybe after it to false again.
-    /// Therefore you cannot use maps like BackInserteBoolMap with this
-    /// algorithm.   
+    /// It must meet the \ref concept::WritedMap "WriteMap" concept.
+    /// Initially it will be setted to false on each edge. After it
+    /// will set all arborescence edges once.
     typedef typename Graph::template EdgeMap<bool> ArborescenceMap; 
 
@@ -73 +72 @@
     }
 
+    /// \brief The type of the PredMap
+    ///
+    /// The type of the PredMap. It is a node map with an edge value type.
+    typedef typename Graph::template NodeMap<typename Graph::Edge> PredMap;
+
+    /// \brief Instantiates a PredMap.
+    ///
+    /// This function instantiates a \ref PredMap. 
+    /// \param _graph is the graph, to which we would like to define the 
+    /// PredMap.
+    static PredMap *createPredMap(const Graph &_graph){
+      return new PredMap(_graph);
+    }
+    
   };
 
@@ -86 +99 @@
   /// given sources and spans all the nodes which are reachable from the
   /// sources. The time complexity of the algorithm is O(n^2 + e).
+  ///
+  /// The algorithm provides also an optimal dual solution to arborescence
+  /// that way the optimality of the algorithm can be proofed easily.
   ///
   /// \param _Graph The graph type the algorithm runs on. The default value
@@ -136 +152 @@
     ///The type of the costs of the edges.
     typedef typename Traits::Value Value;
+    ///The type of the predecessor map.
+    typedef typename Traits::PredMap PredMap;
     ///The type of the map that stores which edges are in the arborescence.
     typedef typename Traits::ArborescenceMap ArborescenceMap;
@@ -158 +176 @@
     };
 
-    const Graph* graph;
-    const CostMap* cost;
-
-    ArborescenceMap* _arborescence_map;
-    bool local_arborescence_map;
-
-    typedef typename Graph::template NodeMap<int> LevelMap;
-    LevelMap *_level;
+    const Graph *graph;
+    const CostMap *cost;
+
+    PredMap *_pred;
+    bool local_pred;
+
+    ArborescenceMap *_arborescence;
+    bool local_arborescence;
+
+    typedef typename Graph::template EdgeMap<int> EdgeOrder;
+    EdgeOrder *_edge_order;
+    
+    typedef typename Graph::template NodeMap<int> NodeOrder;
+    NodeOrder *_node_order;
 
     typedef typename Graph::template NodeMap<CostEdge> CostEdgeMap;
@@ -180 +204 @@
     std::vector<Node> queue;
 
-    int node_counter;
+    typedef std::vector<typename Graph::Node> DualNodeList;
+
+    DualNodeList _dual_node_list;
+
+    struct DualVariable {
+      int begin, end;
+      Value value;
+      
+      DualVariable(int _begin, int _end, Value _value)
+        : begin(_begin), end(_end), value(_value) {}
+
+    };
+
+    typedef std::vector<DualVariable> DualVariables;
+
+    DualVariables _dual_variables;
+
+    typedef typename Graph::template NodeMap<int> HeapCrossRef;
+    
+    HeapCrossRef *_heap_cross_ref;
+
+    typedef BinHeap<Node, int, HeapCrossRef> Heap;
+
+    Heap *_heap;
 
   public:
@@ -209 +256 @@
                                    DefArborescenceMapTraits<T> > Create;
     };
+
+    template <class T>
+    struct DefPredMapTraits : public Traits {
+      typedef T PredMap;
+      static PredMap *createPredMap(const Graph &)
+      {
+        throw UninitializedParameter();
+      }
+    };
+
+    /// \brief \ref named-templ-param "Named parameter" for 
+    /// setting PredMap type
+    ///
+    /// \ref named-templ-param "Named parameter" for setting 
+    /// PredMap type
+    template <class T>
+    struct DefPredMap 
+      : public MinCostArborescence<Graph, CostMap, DefPredMapTraits<T> > {
+      typedef MinCostArborescence<Graph, CostMap, DefPredMapTraits<T> > Create;
+    };
     
     /// @}
@@ -217 +284 @@
     /// \param _cost The cost map used by the algorithm.
     MinCostArborescence(const Graph& _graph, const CostMap& _cost) 
-      : graph(&_graph), cost(&_cost),
-        _arborescence_map(0), local_arborescence_map(false), 
-        _level(0), _cost_edges(0) {}
+      : graph(&_graph), cost(&_cost), _pred(0), local_pred(false),
+        _arborescence(0), local_arborescence(false), 
+        _edge_order(0), _node_order(0), _cost_edges(0), 
+        _heap_cross_ref(0), _heap(0) {}
 
     /// \brief Destructor.
@@ -231 +299 @@
     /// \return \c (*this)
     MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
-      _arborescence_map = &m;
+      if (local_arborescence) {
+        delete _arborescence;
+      }
+      local_arborescence = false;
+      _arborescence = &m;
+      return *this;
+    }
+
+    /// \brief Sets the arborescence map.
+    /// 
+    /// Sets the arborescence map.
+    /// \return \c (*this)
+    MinCostArborescence& predMap(PredMap& m) {
+      if (local_pred) {
+        delete _pred;
+      }
+      local_pred = false;
+      _pred = &m;
       return *this;
     }
@@ -247 +332 @@
     /// Returns a reference to the arborescence map.
     const ArborescenceMap& arborescenceMap() const {
-      return *_arborescence_map;
+      return *_arborescence;
     }
 
@@ -255 +340 @@
     /// \param edge The edge of the graph.
     /// \pre \ref run() must be called before using this function.
-    bool arborescenceEdge(Edge edge) const {
-      return (*_arborescence_map)[edge];
+    bool arborescence(Edge edge) const {
+      return (*_pred)[graph->target(edge)] == edge;
+    }
+
+    /// \brief Returns a reference to the pred map.
+    ///
+    /// Returns a reference to the pred map.
+    const PredMap& predMap() const {
+      return *_pred;
+    }
+
+    /// \brief Returns the predecessor edge of the given node.
+    ///
+    /// Returns the predecessor edge of the given node.
+    bool pred(Node node) const {
+      return (*_pred)[node];
     }
  
@@ -262 +361 @@
     ///
     /// Returns the cost of the arborescence.
-   Value arborescenceCost() const {
+    Value arborescenceValue() const {
       Value sum = 0;
       for (EdgeIt it(*graph); it != INVALID; ++it) {
-        if (arborescenceEdge(it)) {
+        if (arborescence(it)) {
           sum += (*cost)[it];
         }
@@ -271 +370 @@
       return sum;
     }
+
+    /// \brief Indicates that a node is reachable from the sources.
+    ///
+    /// Indicates that a node is reachable from the sources.
+    bool reached(Node node) const {
+      return (*_node_order)[node] != -3;
+    }
+
+    /// \brief Indicates that a node is processed.
+    ///
+    /// Indicates that a node is processed. The arborescence path exists 
+    /// from the source to the given node.
+    bool processed(Node node) const {
+      return (*_node_order)[node] == -1;
+    }
+
+    /// \brief Returns the number of the dual variables in basis.
+    ///
+    /// Returns the number of the dual variables in basis.
+    int dualSize() const {
+      return _dual_variables.size();
+    }
+
+    /// \brief Returns the value of the dual solution.
+    ///
+    /// Returns the value of the dual solution. It should be
+    /// equal to the arborescence value.
+    Value dualValue() const {
+      Value sum = 0;
+      for (int i = 0; i < (int)_dual_variables.size(); ++i) {
+        sum += _dual_variables[i].value;
+      }
+      return sum;
+    }
+
+    /// \brief Returns the number of the nodes in the dual variable.
+    ///
+    /// Returns the number of the nodes in the dual variable.
+    int dualSize(int k) const {
+      return _dual_variables[k].end - _dual_variables[k].begin;
+    }
+
+    /// \brief Returns the value of the dual variable.
+    ///
+    /// Returns the the value of the dual variable.
+    const Value& dualValue(int k) const {
+      return _dual_variables[k].value;
+    }
+
+    /// \brief Lemon iterator for get a dual variable.
+    ///
+    /// Lemon iterator for get a dual variable. This class provides
+    /// a common style lemon iterator which gives back a subset of
+    /// the nodes.
+    class DualIt {
+    public:
+
+      /// \brief Constructor.
+      ///
+      /// Constructor for get the nodeset of the variable. 
+      DualIt(const MinCostArborescence& algorithm, int variable) 
+        : _algorithm(&algorithm), _variable(variable) 
+      {
+        _index = _algorithm->_dual_variables[_variable].begin;
+      }
+
+      /// \brief Invalid constructor.
+      ///
+      /// Invalid constructor.
+      DualIt(Invalid) : _algorithm(0) {}
+
+      /// \brief Conversion to node.
+      ///
+      /// Conversion to node.
+      operator Node() const { 
+        return _algorithm ? _algorithm->_dual_node_list[_index] : INVALID;
+      }
+
+      /// \brief Increment operator.
+      ///
+      /// Increment operator.
+      DualIt& operator++() {
+        ++_index;
+        if (_algorithm->_dual_variables[_variable].end == _index) {
+          _algorithm = 0;
+        }
+        return *this; 
+      }
+
+      bool operator==(const DualIt& it) const { 
+        return (Node)(*this) == (Node)it; 
+      }
+      bool operator!=(const DualIt& it) const { 
+        return (Node)(*this) != (Node)it; 
+      }
+      bool operator<(const DualIt& it) const { 
+        return (Node)(*this) < (Node)it; 
+      }
+      
+    private:
+      const MinCostArborescence* _algorithm;
+      int _variable;
+      int _index;
+    };
 
     /// @}
@@ -280 +483 @@
     /// first you must call \ref init(), then you can add several 
     /// source nodes with \ref addSource().
-    /// Finally \ref start() will perform the actual path
+    /// Finally \ref start() will perform the arborescence
     /// computation.
 
@@ -291 +494 @@
     void init() {
       initStructures();
+      _heap->clear();
       for (NodeIt it(*graph); it != INVALID; ++it) {
         (*_cost_edges)[it].edge = INVALID;
-        (*_level)[it] = -3; 
+        _node_order->set(it, -3); 
+        _heap_cross_ref->set(it, Heap::PRE_HEAP);
       }
       for (EdgeIt it(*graph); it != INVALID; ++it) {
-        _arborescence_map->set(it, false);
-      }
+        _arborescence->set(it, false);
+        _edge_order->set(it, -1);
+      }
+      _dual_node_list.clear();
+      _dual_variables.clear();
     }
 
@@ -310 +518 @@
         nodes.pop_back();
         for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {
-          if ((*_level)[graph->target(it)] == -3) {
-            (*_level)[graph->target(it)] = -2;
-            nodes.push_back(graph->target(it));
-            queue.push_back(graph->target(it));
+          Node target = graph->target(it);
+          if ((*_node_order)[target] == -3) {
+            (*_node_order)[target] = -2;
+            nodes.push_back(target);
+            queue.push_back(target);
           }
         }
       }
-      (*_level)[source] = -1;
+      (*_node_order)[source] = -1;
     }
 
@@ -328 +537 @@
     /// \warning The queue must not be empty!
     Node processNextNode() {
-      node_counter = 0;
       Node node = queue.back();
       queue.pop_back();
-      if ((*_level)[node] == -2) {
+      if ((*_node_order)[node] == -2) {
         Edge edge = prepare(node);
-        while ((*_level)[graph->source(edge)] != -1) {
-          if ((*_level)[graph->source(edge)] >= 0) {
-            edge = contract(bottom((*_level)[graph->source(edge)]));
+        Node source = graph->source(edge);
+        while ((*_node_order)[source] != -1) {
+          if ((*_node_order)[source] >= 0) {
+            edge = contract(source);
           } else {
-            edge = prepare(graph->source(edge));
+            edge = prepare(source);
           }
+          source = graph->source(edge);
         }
-        finalize(graph->target(edge));
+        finalize(edge);
         level_stack.clear();        
       }
@@ -401 +611 @@
 
     void initStructures() {
-      if (!_arborescence_map) {
-        local_arborescence_map = true;
-        _arborescence_map = Traits::createArborescenceMap(*graph);
-      }
-      if (!_level) {
-        _level = new LevelMap(*graph);
+      if (!_pred) {
+        local_pred = true;
+        _pred = Traits::createPredMap(*graph);
+      }
+      if (!_arborescence) {
+        local_arborescence = true;
+        _arborescence = Traits::createArborescenceMap(*graph);
+      }
+      if (!_edge_order) {
+        _edge_order = new EdgeOrder(*graph);
+      }
+      if (!_node_order) {
+        _node_order = new NodeOrder(*graph);
       }
       if (!_cost_edges) {
         _cost_edges = new CostEdgeMap(*graph);
       }
+      if (!_heap_cross_ref) {
+        _heap_cross_ref = new HeapCrossRef(*graph, -1);
+      }
+      if (!_heap) {
+        _heap = new Heap(*_heap_cross_ref);
+      }
     }
 
     void destroyStructures() {
-      if (_level) {
-        delete _level;
+      if (local_arborescence) {
+        delete _arborescence;
+      }
+      if (local_pred) {
+        delete _pred;
+      }
+      if (!_edge_order) {
+        delete _edge_order;
+      }
+      if (_node_order) {
+        delete _node_order;
       }
       if (!_cost_edges) {
         delete _cost_edges;
       }
-      if (local_arborescence_map) {
-        delete _arborescence_map;
+      if (!_heap) {
+        delete _heap;
+      }
+      if (!_heap_cross_ref) {
+        delete _heap_cross_ref;
       }
     }
@@ -427 +662 @@
     Edge prepare(Node node) {
       std::vector<Node> nodes;
-      (*_level)[node] = node_counter;
+      (*_node_order)[node] = _dual_node_list.size();
+      StackLevel level;
+      level.node_level = _dual_node_list.size();
+      _dual_node_list.push_back(node);
       for (InEdgeIt it(*graph, node); it != INVALID; ++it) {
         Edge edge = it;
+        Node source = graph->source(edge);
         Value value = (*cost)[it];
-        if (graph->source(edge) == node || 
-            (*_level)[graph->source(edge)] == -3) continue;
-        if ((*_cost_edges)[graph->source(edge)].edge == INVALID) {
-          (*_cost_edges)[graph->source(edge)].edge = edge;
-          (*_cost_edges)[graph->source(edge)].value = value;
-          nodes.push_back(graph->source(edge));
+        if (source == node || (*_node_order)[source] == -3) continue;
+        if ((*_cost_edges)[source].edge == INVALID) {
+          (*_cost_edges)[source].edge = edge;
+          (*_cost_edges)[source].value = value;
+          nodes.push_back(source);
         } else {
-          if ((*_cost_edges)[graph->source(edge)].value > value) {
-            (*_cost_edges)[graph->source(edge)].edge = edge;
-            (*_cost_edges)[graph->source(edge)].value = value;
+          if ((*_cost_edges)[source].value > value) {
+            (*_cost_edges)[source].edge = edge;
+            (*_cost_edges)[source].value = value;
           }
         }      
@@ -450 +688 @@
         }
       }
-      StackLevel level;
-      level.node_level = node_counter;
+      _edge_order->set(minimum.edge, _dual_variables.size());
+      DualVariable var(_dual_node_list.size() - 1, 
+                       _dual_node_list.size(), minimum.value);
+      _dual_variables.push_back(var);
       for (int i = 0; i < (int)nodes.size(); ++i) {
         (*_cost_edges)[nodes[i]].value -= minimum.value;
@@ -458 +698 @@
       }
       level_stack.push_back(level);
-      ++node_counter;
-      _arborescence_map->set(minimum.edge, true);
       return minimum.edge;
     }
   
-    Edge contract(int node_bottom) {
+    Edge contract(Node node) {
+      int node_bottom = bottom(node);
       std::vector<Node> nodes;
       while (!level_stack.empty() && 
@@ -469 +708 @@
         for (int i = 0; i < (int)level_stack.back().edges.size(); ++i) {
           Edge edge = level_stack.back().edges[i].edge;
+          Node source = graph->source(edge);
           Value value = level_stack.back().edges[i].value;
-          if ((*_level)[graph->source(edge)] >= node_bottom) continue;
-          if ((*_cost_edges)[graph->source(edge)].edge == INVALID) {
-            (*_cost_edges)[graph->source(edge)].edge = edge;
-            (*_cost_edges)[graph->source(edge)].value = value;
-            nodes.push_back(graph->source(edge));
+          if ((*_node_order)[source] >= node_bottom) continue;
+          if ((*_cost_edges)[source].edge == INVALID) {
+            (*_cost_edges)[source].edge = edge;
+            (*_cost_edges)[source].value = value;
+            nodes.push_back(source);
           } else {
-            if ((*_cost_edges)[graph->source(edge)].value > value) {
-              (*_cost_edges)[graph->source(edge)].edge = edge;
-              (*_cost_edges)[graph->source(edge)].value = value;
+            if ((*_cost_edges)[source].value > value) {
+              (*_cost_edges)[source].edge = edge;
+              (*_cost_edges)[source].value = value;
             }
           }
@@ -490 +730 @@
         }
       }
+      _edge_order->set(minimum.edge, _dual_variables.size());
+      DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
+      _dual_variables.push_back(var);
       StackLevel level;
       level.node_level = node_bottom;
@@ -498 +741 @@
       }
       level_stack.push_back(level);
-      _arborescence_map->set(minimum.edge, true);
       return minimum.edge;
     }
 
-    int bottom(int level) {
+    int bottom(Node node) {
       int k = level_stack.size() - 1;
-      while (level_stack[k].node_level > level) {
+      while (level_stack[k].node_level > (*_node_order)[node]) {
         --k;
       }
@@ -510 +752 @@
     }
 
-    void finalize(Node source) {
-      std::vector<Node> nodes;
-      nodes.push_back(source);
-      while (!nodes.empty()) {
-        Node node = nodes.back();
-        nodes.pop_back();
-        for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {
-          if ((*_level)[graph->target(it)] >= 0 && (*_arborescence_map)[it]) {
-            (*_level)[graph->target(it)] = -1;
-            nodes.push_back(graph->target(it));
-          } else {
-            _arborescence_map->set(it, false);
+    void finalize(Edge edge) {
+      Node node = graph->target(edge);
+      _heap->push(node, (*_edge_order)[edge]);
+      _pred->set(node, edge);
+      while (!_heap->empty()) {
+        Node source = _heap->top();
+        _heap->pop();
+        _node_order->set(source, -1);
+        for (OutEdgeIt it(*graph, source); it != INVALID; ++it) {
+          if ((*_edge_order)[it] < 0) continue;
+          Node target = graph->target(it);
+          switch(_heap->state(target)) {
+          case Heap::PRE_HEAP:
+            _heap->push(target, (*_edge_order)[it]); 
+            _pred->set(target, it);
+            break;
+          case Heap::IN_HEAP:
+            if ((*_edge_order)[it] < (*_heap)[target]) {
+              _heap->decrease(target, (*_edge_order)[it]); 
+              _pred->set(target, it);
+            }
+            break;
+          case Heap::POST_HEAP:
+            break;
           }
         }
-      }
-      (*_level)[source] = -1;      
+        _arborescence->set((*_pred)[source], true);
+      }
     }
 
@@ -552 +806 @@
     mca.arborescenceMap(arborescence);
     mca.run(source);
-    return mca.arborescenceCost();
+    return mca.arborescenceValue();
   }
 
@@ -558 +812 @@
 
 #endif
-
-// Hilbert - Huang