COIN-OR::LEMON - Graph Library

Ignore:
Timestamp:
03/31/06 13:10:44 (13 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:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/min_cost_arborescence.h

    r2017 r2025  
    2727
    2828#include <lemon/list_graph.h>
     29#include <lemon/bin_heap.h>
    2930
    3031namespace lemon {
     
    5758    ///
    5859    /// The type of the map that stores which edges are in the arborescence.
    59     /// It must meet the \ref concept::ReadWriteMap "ReadWriteMap" concept.
    60     /// Initially it will be setted to false on each edge. The algorithm
    61     /// may set each value one time to true and maybe after it to false again.
    62     /// Therefore you cannot use maps like BackInserteBoolMap with this
    63     /// algorithm.   
     60    /// It must meet the \ref concept::WritedMap "WriteMap" concept.
     61    /// Initially it will be setted to false on each edge. After it
     62    /// will set all arborescence edges once.
    6463    typedef typename Graph::template EdgeMap<bool> ArborescenceMap;
    6564
     
    7372    }
    7473
     74    /// \brief The type of the PredMap
     75    ///
     76    /// The type of the PredMap. It is a node map with an edge value type.
     77    typedef typename Graph::template NodeMap<typename Graph::Edge> PredMap;
     78
     79    /// \brief Instantiates a PredMap.
     80    ///
     81    /// This function instantiates a \ref PredMap.
     82    /// \param _graph is the graph, to which we would like to define the
     83    /// PredMap.
     84    static PredMap *createPredMap(const Graph &_graph){
     85      return new PredMap(_graph);
     86    }
     87   
    7588  };
    7689
     
    8699  /// given sources and spans all the nodes which are reachable from the
    87100  /// sources. The time complexity of the algorithm is O(n^2 + e).
     101  ///
     102  /// The algorithm provides also an optimal dual solution to arborescence
     103  /// that way the optimality of the algorithm can be proofed easily.
    88104  ///
    89105  /// \param _Graph The graph type the algorithm runs on. The default value
     
    136152    ///The type of the costs of the edges.
    137153    typedef typename Traits::Value Value;
     154    ///The type of the predecessor map.
     155    typedef typename Traits::PredMap PredMap;
    138156    ///The type of the map that stores which edges are in the arborescence.
    139157    typedef typename Traits::ArborescenceMap ArborescenceMap;
     
    158176    };
    159177
    160     const Graph* graph;
    161     const CostMap* cost;
    162 
    163     ArborescenceMap* _arborescence_map;
    164     bool local_arborescence_map;
    165 
    166     typedef typename Graph::template NodeMap<int> LevelMap;
    167     LevelMap *_level;
     178    const Graph *graph;
     179    const CostMap *cost;
     180
     181    PredMap *_pred;
     182    bool local_pred;
     183
     184    ArborescenceMap *_arborescence;
     185    bool local_arborescence;
     186
     187    typedef typename Graph::template EdgeMap<int> EdgeOrder;
     188    EdgeOrder *_edge_order;
     189   
     190    typedef typename Graph::template NodeMap<int> NodeOrder;
     191    NodeOrder *_node_order;
    168192
    169193    typedef typename Graph::template NodeMap<CostEdge> CostEdgeMap;
     
    180204    std::vector<Node> queue;
    181205
    182     int node_counter;
     206    typedef std::vector<typename Graph::Node> DualNodeList;
     207
     208    DualNodeList _dual_node_list;
     209
     210    struct DualVariable {
     211      int begin, end;
     212      Value value;
     213     
     214      DualVariable(int _begin, int _end, Value _value)
     215        : begin(_begin), end(_end), value(_value) {}
     216
     217    };
     218
     219    typedef std::vector<DualVariable> DualVariables;
     220
     221    DualVariables _dual_variables;
     222
     223    typedef typename Graph::template NodeMap<int> HeapCrossRef;
     224   
     225    HeapCrossRef *_heap_cross_ref;
     226
     227    typedef BinHeap<Node, int, HeapCrossRef> Heap;
     228
     229    Heap *_heap;
    183230
    184231  public:
     
    209256                                   DefArborescenceMapTraits<T> > Create;
    210257    };
     258
     259    template <class T>
     260    struct DefPredMapTraits : public Traits {
     261      typedef T PredMap;
     262      static PredMap *createPredMap(const Graph &)
     263      {
     264        throw UninitializedParameter();
     265      }
     266    };
     267
     268    /// \brief \ref named-templ-param "Named parameter" for
     269    /// setting PredMap type
     270    ///
     271    /// \ref named-templ-param "Named parameter" for setting
     272    /// PredMap type
     273    template <class T>
     274    struct DefPredMap
     275      : public MinCostArborescence<Graph, CostMap, DefPredMapTraits<T> > {
     276      typedef MinCostArborescence<Graph, CostMap, DefPredMapTraits<T> > Create;
     277    };
    211278   
    212279    /// @}
     
    217284    /// \param _cost The cost map used by the algorithm.
    218285    MinCostArborescence(const Graph& _graph, const CostMap& _cost)
    219       : graph(&_graph), cost(&_cost),
    220         _arborescence_map(0), local_arborescence_map(false),
    221         _level(0), _cost_edges(0) {}
     286      : graph(&_graph), cost(&_cost), _pred(0), local_pred(false),
     287        _arborescence(0), local_arborescence(false),
     288        _edge_order(0), _node_order(0), _cost_edges(0),
     289        _heap_cross_ref(0), _heap(0) {}
    222290
    223291    /// \brief Destructor.
     
    231299    /// \return \c (*this)
    232300    MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
    233       _arborescence_map = &m;
     301      if (local_arborescence) {
     302        delete _arborescence;
     303      }
     304      local_arborescence = false;
     305      _arborescence = &m;
     306      return *this;
     307    }
     308
     309    /// \brief Sets the arborescence map.
     310    ///
     311    /// Sets the arborescence map.
     312    /// \return \c (*this)
     313    MinCostArborescence& predMap(PredMap& m) {
     314      if (local_pred) {
     315        delete _pred;
     316      }
     317      local_pred = false;
     318      _pred = &m;
    234319      return *this;
    235320    }
     
    247332    /// Returns a reference to the arborescence map.
    248333    const ArborescenceMap& arborescenceMap() const {
    249       return *_arborescence_map;
     334      return *_arborescence;
    250335    }
    251336
     
    255340    /// \param edge The edge of the graph.
    256341    /// \pre \ref run() must be called before using this function.
    257     bool arborescenceEdge(Edge edge) const {
    258       return (*_arborescence_map)[edge];
     342    bool arborescence(Edge edge) const {
     343      return (*_pred)[graph->target(edge)] == edge;
     344    }
     345
     346    /// \brief Returns a reference to the pred map.
     347    ///
     348    /// Returns a reference to the pred map.
     349    const PredMap& predMap() const {
     350      return *_pred;
     351    }
     352
     353    /// \brief Returns the predecessor edge of the given node.
     354    ///
     355    /// Returns the predecessor edge of the given node.
     356    bool pred(Node node) const {
     357      return (*_pred)[node];
    259358    }
    260359 
     
    262361    ///
    263362    /// Returns the cost of the arborescence.
    264    Value arborescenceCost() const {
     363    Value arborescenceValue() const {
    265364      Value sum = 0;
    266365      for (EdgeIt it(*graph); it != INVALID; ++it) {
    267         if (arborescenceEdge(it)) {
     366        if (arborescence(it)) {
    268367          sum += (*cost)[it];
    269368        }
     
    271370      return sum;
    272371    }
     372
     373    /// \brief Indicates that a node is reachable from the sources.
     374    ///
     375    /// Indicates that a node is reachable from the sources.
     376    bool reached(Node node) const {
     377      return (*_node_order)[node] != -3;
     378    }
     379
     380    /// \brief Indicates that a node is processed.
     381    ///
     382    /// Indicates that a node is processed. The arborescence path exists
     383    /// from the source to the given node.
     384    bool processed(Node node) const {
     385      return (*_node_order)[node] == -1;
     386    }
     387
     388    /// \brief Returns the number of the dual variables in basis.
     389    ///
     390    /// Returns the number of the dual variables in basis.
     391    int dualSize() const {
     392      return _dual_variables.size();
     393    }
     394
     395    /// \brief Returns the value of the dual solution.
     396    ///
     397    /// Returns the value of the dual solution. It should be
     398    /// equal to the arborescence value.
     399    Value dualValue() const {
     400      Value sum = 0;
     401      for (int i = 0; i < (int)_dual_variables.size(); ++i) {
     402        sum += _dual_variables[i].value;
     403      }
     404      return sum;
     405    }
     406
     407    /// \brief Returns the number of the nodes in the dual variable.
     408    ///
     409    /// Returns the number of the nodes in the dual variable.
     410    int dualSize(int k) const {
     411      return _dual_variables[k].end - _dual_variables[k].begin;
     412    }
     413
     414    /// \brief Returns the value of the dual variable.
     415    ///
     416    /// Returns the the value of the dual variable.
     417    const Value& dualValue(int k) const {
     418      return _dual_variables[k].value;
     419    }
     420
     421    /// \brief Lemon iterator for get a dual variable.
     422    ///
     423    /// Lemon iterator for get a dual variable. This class provides
     424    /// a common style lemon iterator which gives back a subset of
     425    /// the nodes.
     426    class DualIt {
     427    public:
     428
     429      /// \brief Constructor.
     430      ///
     431      /// Constructor for get the nodeset of the variable.
     432      DualIt(const MinCostArborescence& algorithm, int variable)
     433        : _algorithm(&algorithm), _variable(variable)
     434      {
     435        _index = _algorithm->_dual_variables[_variable].begin;
     436      }
     437
     438      /// \brief Invalid constructor.
     439      ///
     440      /// Invalid constructor.
     441      DualIt(Invalid) : _algorithm(0) {}
     442
     443      /// \brief Conversion to node.
     444      ///
     445      /// Conversion to node.
     446      operator Node() const {
     447        return _algorithm ? _algorithm->_dual_node_list[_index] : INVALID;
     448      }
     449
     450      /// \brief Increment operator.
     451      ///
     452      /// Increment operator.
     453      DualIt& operator++() {
     454        ++_index;
     455        if (_algorithm->_dual_variables[_variable].end == _index) {
     456          _algorithm = 0;
     457        }
     458        return *this;
     459      }
     460
     461      bool operator==(const DualIt& it) const {
     462        return (Node)(*this) == (Node)it;
     463      }
     464      bool operator!=(const DualIt& it) const {
     465        return (Node)(*this) != (Node)it;
     466      }
     467      bool operator<(const DualIt& it) const {
     468        return (Node)(*this) < (Node)it;
     469      }
     470     
     471    private:
     472      const MinCostArborescence* _algorithm;
     473      int _variable;
     474      int _index;
     475    };
    273476
    274477    /// @}
     
    280483    /// first you must call \ref init(), then you can add several
    281484    /// source nodes with \ref addSource().
    282     /// Finally \ref start() will perform the actual path
     485    /// Finally \ref start() will perform the arborescence
    283486    /// computation.
    284487
     
    291494    void init() {
    292495      initStructures();
     496      _heap->clear();
    293497      for (NodeIt it(*graph); it != INVALID; ++it) {
    294498        (*_cost_edges)[it].edge = INVALID;
    295         (*_level)[it] = -3;
     499        _node_order->set(it, -3);
     500        _heap_cross_ref->set(it, Heap::PRE_HEAP);
    296501      }
    297502      for (EdgeIt it(*graph); it != INVALID; ++it) {
    298         _arborescence_map->set(it, false);
    299       }
     503        _arborescence->set(it, false);
     504        _edge_order->set(it, -1);
     505      }
     506      _dual_node_list.clear();
     507      _dual_variables.clear();
    300508    }
    301509
     
    310518        nodes.pop_back();
    311519        for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {
    312           if ((*_level)[graph->target(it)] == -3) {
    313             (*_level)[graph->target(it)] = -2;
    314             nodes.push_back(graph->target(it));
    315             queue.push_back(graph->target(it));
     520          Node target = graph->target(it);
     521          if ((*_node_order)[target] == -3) {
     522            (*_node_order)[target] = -2;
     523            nodes.push_back(target);
     524            queue.push_back(target);
    316525          }
    317526        }
    318527      }
    319       (*_level)[source] = -1;
     528      (*_node_order)[source] = -1;
    320529    }
    321530
     
    328537    /// \warning The queue must not be empty!
    329538    Node processNextNode() {
    330       node_counter = 0;
    331539      Node node = queue.back();
    332540      queue.pop_back();
    333       if ((*_level)[node] == -2) {
     541      if ((*_node_order)[node] == -2) {
    334542        Edge edge = prepare(node);
    335         while ((*_level)[graph->source(edge)] != -1) {
    336           if ((*_level)[graph->source(edge)] >= 0) {
    337             edge = contract(bottom((*_level)[graph->source(edge)]));
     543        Node source = graph->source(edge);
     544        while ((*_node_order)[source] != -1) {
     545          if ((*_node_order)[source] >= 0) {
     546            edge = contract(source);
    338547          } else {
    339             edge = prepare(graph->source(edge));
     548            edge = prepare(source);
    340549          }
     550          source = graph->source(edge);
    341551        }
    342         finalize(graph->target(edge));
     552        finalize(edge);
    343553        level_stack.clear();       
    344554      }
     
    401611
    402612    void initStructures() {
    403       if (!_arborescence_map) {
    404         local_arborescence_map = true;
    405         _arborescence_map = Traits::createArborescenceMap(*graph);
    406       }
    407       if (!_level) {
    408         _level = new LevelMap(*graph);
     613      if (!_pred) {
     614        local_pred = true;
     615        _pred = Traits::createPredMap(*graph);
     616      }
     617      if (!_arborescence) {
     618        local_arborescence = true;
     619        _arborescence = Traits::createArborescenceMap(*graph);
     620      }
     621      if (!_edge_order) {
     622        _edge_order = new EdgeOrder(*graph);
     623      }
     624      if (!_node_order) {
     625        _node_order = new NodeOrder(*graph);
    409626      }
    410627      if (!_cost_edges) {
    411628        _cost_edges = new CostEdgeMap(*graph);
    412629      }
     630      if (!_heap_cross_ref) {
     631        _heap_cross_ref = new HeapCrossRef(*graph, -1);
     632      }
     633      if (!_heap) {
     634        _heap = new Heap(*_heap_cross_ref);
     635      }
    413636    }
    414637
    415638    void destroyStructures() {
    416       if (_level) {
    417         delete _level;
     639      if (local_arborescence) {
     640        delete _arborescence;
     641      }
     642      if (local_pred) {
     643        delete _pred;
     644      }
     645      if (!_edge_order) {
     646        delete _edge_order;
     647      }
     648      if (_node_order) {
     649        delete _node_order;
    418650      }
    419651      if (!_cost_edges) {
    420652        delete _cost_edges;
    421653      }
    422       if (local_arborescence_map) {
    423         delete _arborescence_map;
     654      if (!_heap) {
     655        delete _heap;
     656      }
     657      if (!_heap_cross_ref) {
     658        delete _heap_cross_ref;
    424659      }
    425660    }
     
    427662    Edge prepare(Node node) {
    428663      std::vector<Node> nodes;
    429       (*_level)[node] = node_counter;
     664      (*_node_order)[node] = _dual_node_list.size();
     665      StackLevel level;
     666      level.node_level = _dual_node_list.size();
     667      _dual_node_list.push_back(node);
    430668      for (InEdgeIt it(*graph, node); it != INVALID; ++it) {
    431669        Edge edge = it;
     670        Node source = graph->source(edge);
    432671        Value value = (*cost)[it];
    433         if (graph->source(edge) == node ||
    434             (*_level)[graph->source(edge)] == -3) continue;
    435         if ((*_cost_edges)[graph->source(edge)].edge == INVALID) {
    436           (*_cost_edges)[graph->source(edge)].edge = edge;
    437           (*_cost_edges)[graph->source(edge)].value = value;
    438           nodes.push_back(graph->source(edge));
     672        if (source == node || (*_node_order)[source] == -3) continue;
     673        if ((*_cost_edges)[source].edge == INVALID) {
     674          (*_cost_edges)[source].edge = edge;
     675          (*_cost_edges)[source].value = value;
     676          nodes.push_back(source);
    439677        } else {
    440           if ((*_cost_edges)[graph->source(edge)].value > value) {
    441             (*_cost_edges)[graph->source(edge)].edge = edge;
    442             (*_cost_edges)[graph->source(edge)].value = value;
     678          if ((*_cost_edges)[source].value > value) {
     679            (*_cost_edges)[source].edge = edge;
     680            (*_cost_edges)[source].value = value;
    443681          }
    444682        }     
     
    450688        }
    451689      }
    452       StackLevel level;
    453       level.node_level = node_counter;
     690      _edge_order->set(minimum.edge, _dual_variables.size());
     691      DualVariable var(_dual_node_list.size() - 1,
     692                       _dual_node_list.size(), minimum.value);
     693      _dual_variables.push_back(var);
    454694      for (int i = 0; i < (int)nodes.size(); ++i) {
    455695        (*_cost_edges)[nodes[i]].value -= minimum.value;
     
    458698      }
    459699      level_stack.push_back(level);
    460       ++node_counter;
    461       _arborescence_map->set(minimum.edge, true);
    462700      return minimum.edge;
    463701    }
    464702 
    465     Edge contract(int node_bottom) {
     703    Edge contract(Node node) {
     704      int node_bottom = bottom(node);
    466705      std::vector<Node> nodes;
    467706      while (!level_stack.empty() &&
     
    469708        for (int i = 0; i < (int)level_stack.back().edges.size(); ++i) {
    470709          Edge edge = level_stack.back().edges[i].edge;
     710          Node source = graph->source(edge);
    471711          Value value = level_stack.back().edges[i].value;
    472           if ((*_level)[graph->source(edge)] >= node_bottom) continue;
    473           if ((*_cost_edges)[graph->source(edge)].edge == INVALID) {
    474             (*_cost_edges)[graph->source(edge)].edge = edge;
    475             (*_cost_edges)[graph->source(edge)].value = value;
    476             nodes.push_back(graph->source(edge));
     712          if ((*_node_order)[source] >= node_bottom) continue;
     713          if ((*_cost_edges)[source].edge == INVALID) {
     714            (*_cost_edges)[source].edge = edge;
     715            (*_cost_edges)[source].value = value;
     716            nodes.push_back(source);
    477717          } else {
    478             if ((*_cost_edges)[graph->source(edge)].value > value) {
    479               (*_cost_edges)[graph->source(edge)].edge = edge;
    480               (*_cost_edges)[graph->source(edge)].value = value;
     718            if ((*_cost_edges)[source].value > value) {
     719              (*_cost_edges)[source].edge = edge;
     720              (*_cost_edges)[source].value = value;
    481721            }
    482722          }
     
    490730        }
    491731      }
     732      _edge_order->set(minimum.edge, _dual_variables.size());
     733      DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
     734      _dual_variables.push_back(var);
    492735      StackLevel level;
    493736      level.node_level = node_bottom;
     
    498741      }
    499742      level_stack.push_back(level);
    500       _arborescence_map->set(minimum.edge, true);
    501743      return minimum.edge;
    502744    }
    503745
    504     int bottom(int level) {
     746    int bottom(Node node) {
    505747      int k = level_stack.size() - 1;
    506       while (level_stack[k].node_level > level) {
     748      while (level_stack[k].node_level > (*_node_order)[node]) {
    507749        --k;
    508750      }
     
    510752    }
    511753
    512     void finalize(Node source) {
    513       std::vector<Node> nodes;
    514       nodes.push_back(source);
    515       while (!nodes.empty()) {
    516         Node node = nodes.back();
    517         nodes.pop_back();
    518         for (OutEdgeIt it(*graph, node); it != INVALID; ++it) {
    519           if ((*_level)[graph->target(it)] >= 0 && (*_arborescence_map)[it]) {
    520             (*_level)[graph->target(it)] = -1;
    521             nodes.push_back(graph->target(it));
    522           } else {
    523             _arborescence_map->set(it, false);
     754    void finalize(Edge edge) {
     755      Node node = graph->target(edge);
     756      _heap->push(node, (*_edge_order)[edge]);
     757      _pred->set(node, edge);
     758      while (!_heap->empty()) {
     759        Node source = _heap->top();
     760        _heap->pop();
     761        _node_order->set(source, -1);
     762        for (OutEdgeIt it(*graph, source); it != INVALID; ++it) {
     763          if ((*_edge_order)[it] < 0) continue;
     764          Node target = graph->target(it);
     765          switch(_heap->state(target)) {
     766          case Heap::PRE_HEAP:
     767            _heap->push(target, (*_edge_order)[it]);
     768            _pred->set(target, it);
     769            break;
     770          case Heap::IN_HEAP:
     771            if ((*_edge_order)[it] < (*_heap)[target]) {
     772              _heap->decrease(target, (*_edge_order)[it]);
     773              _pred->set(target, it);
     774            }
     775            break;
     776          case Heap::POST_HEAP:
     777            break;
    524778          }
    525779        }
    526       }
    527       (*_level)[source] = -1;     
     780        _arborescence->set((*_pred)[source], true);
     781      }
    528782    }
    529783
     
    552806    mca.arborescenceMap(arborescence);
    553807    mca.run(source);
    554     return mca.arborescenceCost();
     808    return mca.arborescenceValue();
    555809  }
    556810
     
    558812
    559813#endif
    560 
    561 // Hilbert - Huang
Note: See TracChangeset for help on using the changeset viewer.