lemon/min_cost_arborescence.h
author Peter Kovacs <kpeter@inf.elte.hu>
Wed, 06 May 2009 14:46:05 +0200
changeset 649 76cbcb3e9bbb
parent 584 33c6b6e755cd
child 713 4ac30454f1c1
permissions -rw-r--r--
Add a test file for the connectivity functions (#285)
The Euler tools have a separate test file.
deba@490
     1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
deba@490
     2
 *
deba@490
     3
 * This file is a part of LEMON, a generic C++ optimization library.
deba@490
     4
 *
deba@490
     5
 * Copyright (C) 2003-2008
deba@490
     6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@490
     7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@490
     8
 *
deba@490
     9
 * Permission to use, modify and distribute this software is granted
deba@490
    10
 * provided that this copyright notice appears in all copies. For
deba@490
    11
 * precise terms see the accompanying LICENSE file.
deba@490
    12
 *
deba@490
    13
 * This software is provided "AS IS" with no warranty of any kind,
deba@490
    14
 * express or implied, and with no claim as to its suitability for any
deba@490
    15
 * purpose.
deba@490
    16
 *
deba@490
    17
 */
deba@490
    18
deba@490
    19
#ifndef LEMON_MIN_COST_ARBORESCENCE_H
deba@490
    20
#define LEMON_MIN_COST_ARBORESCENCE_H
deba@490
    21
deba@490
    22
///\ingroup spantree
deba@490
    23
///\file
deba@490
    24
///\brief Minimum Cost Arborescence algorithm.
deba@490
    25
deba@490
    26
#include <vector>
deba@490
    27
deba@490
    28
#include <lemon/list_graph.h>
deba@490
    29
#include <lemon/bin_heap.h>
deba@490
    30
#include <lemon/assert.h>
deba@490
    31
deba@490
    32
namespace lemon {
deba@490
    33
deba@490
    34
deba@490
    35
  /// \brief Default traits class for MinCostArborescence class.
deba@490
    36
  ///
deba@490
    37
  /// Default traits class for MinCostArborescence class.
kpeter@559
    38
  /// \param GR Digraph type.
kpeter@625
    39
  /// \param CM Type of the cost map.
kpeter@559
    40
  template <class GR, class CM>
deba@490
    41
  struct MinCostArborescenceDefaultTraits{
deba@490
    42
deba@490
    43
    /// \brief The digraph type the algorithm runs on.
kpeter@559
    44
    typedef GR Digraph;
deba@490
    45
deba@490
    46
    /// \brief The type of the map that stores the arc costs.
deba@490
    47
    ///
deba@490
    48
    /// The type of the map that stores the arc costs.
kpeter@625
    49
    /// It must conform to the \ref concepts::ReadMap "ReadMap" concept.
kpeter@559
    50
    typedef CM CostMap;
deba@490
    51
deba@490
    52
    /// \brief The value type of the costs.
deba@490
    53
    ///
deba@490
    54
    /// The value type of the costs.
deba@490
    55
    typedef typename CostMap::Value Value;
deba@490
    56
deba@490
    57
    /// \brief The type of the map that stores which arcs are in the
deba@490
    58
    /// arborescence.
deba@490
    59
    ///
deba@490
    60
    /// The type of the map that stores which arcs are in the
kpeter@625
    61
    /// arborescence.  It must conform to the \ref concepts::WriteMap
kpeter@625
    62
    /// "WriteMap" concept, and its value type must be \c bool
kpeter@625
    63
    /// (or convertible). Initially it will be set to \c false on each
kpeter@625
    64
    /// arc, then it will be set on each arborescence arc once.
deba@490
    65
    typedef typename Digraph::template ArcMap<bool> ArborescenceMap;
deba@490
    66
kpeter@559
    67
    /// \brief Instantiates a \c ArborescenceMap.
deba@490
    68
    ///
kpeter@559
    69
    /// This function instantiates a \c ArborescenceMap.
kpeter@625
    70
    /// \param digraph The digraph to which we would like to calculate
kpeter@625
    71
    /// the \c ArborescenceMap.
deba@490
    72
    static ArborescenceMap *createArborescenceMap(const Digraph &digraph){
deba@490
    73
      return new ArborescenceMap(digraph);
deba@490
    74
    }
deba@490
    75
kpeter@559
    76
    /// \brief The type of the \c PredMap
deba@490
    77
    ///
kpeter@625
    78
    /// The type of the \c PredMap. It must confrom to the
kpeter@625
    79
    /// \ref concepts::WriteMap "WriteMap" concept, and its value type
kpeter@625
    80
    /// must be the \c Arc type of the digraph.
deba@490
    81
    typedef typename Digraph::template NodeMap<typename Digraph::Arc> PredMap;
deba@490
    82
kpeter@559
    83
    /// \brief Instantiates a \c PredMap.
deba@490
    84
    ///
kpeter@559
    85
    /// This function instantiates a \c PredMap.
kpeter@559
    86
    /// \param digraph The digraph to which we would like to define the
kpeter@559
    87
    /// \c PredMap.
deba@490
    88
    static PredMap *createPredMap(const Digraph &digraph){
deba@490
    89
      return new PredMap(digraph);
deba@490
    90
    }
deba@490
    91
deba@490
    92
  };
deba@490
    93
deba@490
    94
  /// \ingroup spantree
deba@490
    95
  ///
kpeter@584
    96
  /// \brief Minimum Cost Arborescence algorithm class.
deba@490
    97
  ///
kpeter@625
    98
  /// This class provides an efficient implementation of the
kpeter@584
    99
  /// Minimum Cost Arborescence algorithm. The arborescence is a tree
deba@490
   100
  /// which is directed from a given source node of the digraph. One or
kpeter@625
   101
  /// more sources should be given to the algorithm and it will calculate
kpeter@625
   102
  /// the minimum cost subgraph that is the union of arborescences with the
deba@490
   103
  /// given sources and spans all the nodes which are reachable from the
kpeter@559
   104
  /// sources. The time complexity of the algorithm is O(n<sup>2</sup>+e).
deba@490
   105
  ///
kpeter@625
   106
  /// The algorithm also provides an optimal dual solution, therefore
deba@490
   107
  /// the optimality of the solution can be checked.
deba@490
   108
  ///
kpeter@625
   109
  /// \param GR The digraph type the algorithm runs on.
kpeter@625
   110
  /// \param CM A read-only arc map storing the costs of the
deba@490
   111
  /// arcs. It is read once for each arc, so the map may involve in
kpeter@625
   112
  /// relatively time consuming process to compute the arc costs if
deba@490
   113
  /// it is necessary. The default map type is \ref
deba@490
   114
  /// concepts::Digraph::ArcMap "Digraph::ArcMap<int>".
kpeter@559
   115
  /// \param TR Traits class to set various data types used
deba@490
   116
  /// by the algorithm. The default traits class is
deba@490
   117
  /// \ref MinCostArborescenceDefaultTraits
kpeter@625
   118
  /// "MinCostArborescenceDefaultTraits<GR, CM>".
deba@490
   119
#ifndef DOXYGEN
kpeter@625
   120
  template <typename GR,
kpeter@559
   121
            typename CM = typename GR::template ArcMap<int>,
kpeter@559
   122
            typename TR =
kpeter@559
   123
              MinCostArborescenceDefaultTraits<GR, CM> >
deba@490
   124
#else
kpeter@559
   125
  template <typename GR, typename CM, typedef TR>
deba@490
   126
#endif
deba@490
   127
  class MinCostArborescence {
deba@490
   128
  public:
deba@490
   129
kpeter@625
   130
    /// \brief The \ref MinCostArborescenceDefaultTraits "traits class" 
kpeter@625
   131
    /// of the algorithm. 
kpeter@559
   132
    typedef TR Traits;
deba@490
   133
    /// The type of the underlying digraph.
deba@490
   134
    typedef typename Traits::Digraph Digraph;
deba@490
   135
    /// The type of the map that stores the arc costs.
deba@490
   136
    typedef typename Traits::CostMap CostMap;
deba@490
   137
    ///The type of the costs of the arcs.
deba@490
   138
    typedef typename Traits::Value Value;
deba@490
   139
    ///The type of the predecessor map.
deba@490
   140
    typedef typename Traits::PredMap PredMap;
deba@490
   141
    ///The type of the map that stores which arcs are in the arborescence.
deba@490
   142
    typedef typename Traits::ArborescenceMap ArborescenceMap;
deba@490
   143
deba@490
   144
    typedef MinCostArborescence Create;
deba@490
   145
deba@490
   146
  private:
deba@490
   147
deba@490
   148
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
deba@490
   149
deba@490
   150
    struct CostArc {
deba@490
   151
deba@490
   152
      Arc arc;
deba@490
   153
      Value value;
deba@490
   154
deba@490
   155
      CostArc() {}
deba@490
   156
      CostArc(Arc _arc, Value _value) : arc(_arc), value(_value) {}
deba@490
   157
deba@490
   158
    };
deba@490
   159
deba@490
   160
    const Digraph *_digraph;
deba@490
   161
    const CostMap *_cost;
deba@490
   162
deba@490
   163
    PredMap *_pred;
deba@490
   164
    bool local_pred;
deba@490
   165
deba@490
   166
    ArborescenceMap *_arborescence;
deba@490
   167
    bool local_arborescence;
deba@490
   168
deba@490
   169
    typedef typename Digraph::template ArcMap<int> ArcOrder;
deba@490
   170
    ArcOrder *_arc_order;
deba@490
   171
deba@490
   172
    typedef typename Digraph::template NodeMap<int> NodeOrder;
deba@490
   173
    NodeOrder *_node_order;
deba@490
   174
deba@490
   175
    typedef typename Digraph::template NodeMap<CostArc> CostArcMap;
deba@490
   176
    CostArcMap *_cost_arcs;
deba@490
   177
deba@490
   178
    struct StackLevel {
deba@490
   179
deba@490
   180
      std::vector<CostArc> arcs;
deba@490
   181
      int node_level;
deba@490
   182
deba@490
   183
    };
deba@490
   184
deba@490
   185
    std::vector<StackLevel> level_stack;
deba@490
   186
    std::vector<Node> queue;
deba@490
   187
deba@490
   188
    typedef std::vector<typename Digraph::Node> DualNodeList;
deba@490
   189
deba@490
   190
    DualNodeList _dual_node_list;
deba@490
   191
deba@490
   192
    struct DualVariable {
deba@490
   193
      int begin, end;
deba@490
   194
      Value value;
deba@490
   195
deba@490
   196
      DualVariable(int _begin, int _end, Value _value)
deba@490
   197
        : begin(_begin), end(_end), value(_value) {}
deba@490
   198
deba@490
   199
    };
deba@490
   200
deba@490
   201
    typedef std::vector<DualVariable> DualVariables;
deba@490
   202
deba@490
   203
    DualVariables _dual_variables;
deba@490
   204
deba@490
   205
    typedef typename Digraph::template NodeMap<int> HeapCrossRef;
deba@490
   206
deba@490
   207
    HeapCrossRef *_heap_cross_ref;
deba@490
   208
deba@490
   209
    typedef BinHeap<int, HeapCrossRef> Heap;
deba@490
   210
deba@490
   211
    Heap *_heap;
deba@490
   212
deba@490
   213
  protected:
deba@490
   214
deba@490
   215
    MinCostArborescence() {}
deba@490
   216
deba@490
   217
  private:
deba@490
   218
deba@490
   219
    void createStructures() {
deba@490
   220
      if (!_pred) {
deba@490
   221
        local_pred = true;
deba@490
   222
        _pred = Traits::createPredMap(*_digraph);
deba@490
   223
      }
deba@490
   224
      if (!_arborescence) {
deba@490
   225
        local_arborescence = true;
deba@490
   226
        _arborescence = Traits::createArborescenceMap(*_digraph);
deba@490
   227
      }
deba@490
   228
      if (!_arc_order) {
deba@490
   229
        _arc_order = new ArcOrder(*_digraph);
deba@490
   230
      }
deba@490
   231
      if (!_node_order) {
deba@490
   232
        _node_order = new NodeOrder(*_digraph);
deba@490
   233
      }
deba@490
   234
      if (!_cost_arcs) {
deba@490
   235
        _cost_arcs = new CostArcMap(*_digraph);
deba@490
   236
      }
deba@490
   237
      if (!_heap_cross_ref) {
deba@490
   238
        _heap_cross_ref = new HeapCrossRef(*_digraph, -1);
deba@490
   239
      }
deba@490
   240
      if (!_heap) {
deba@490
   241
        _heap = new Heap(*_heap_cross_ref);
deba@490
   242
      }
deba@490
   243
    }
deba@490
   244
deba@490
   245
    void destroyStructures() {
deba@490
   246
      if (local_arborescence) {
deba@490
   247
        delete _arborescence;
deba@490
   248
      }
deba@490
   249
      if (local_pred) {
deba@490
   250
        delete _pred;
deba@490
   251
      }
deba@490
   252
      if (_arc_order) {
deba@490
   253
        delete _arc_order;
deba@490
   254
      }
deba@490
   255
      if (_node_order) {
deba@490
   256
        delete _node_order;
deba@490
   257
      }
deba@490
   258
      if (_cost_arcs) {
deba@490
   259
        delete _cost_arcs;
deba@490
   260
      }
deba@490
   261
      if (_heap) {
deba@490
   262
        delete _heap;
deba@490
   263
      }
deba@490
   264
      if (_heap_cross_ref) {
deba@490
   265
        delete _heap_cross_ref;
deba@490
   266
      }
deba@490
   267
    }
deba@490
   268
deba@490
   269
    Arc prepare(Node node) {
deba@490
   270
      std::vector<Node> nodes;
deba@490
   271
      (*_node_order)[node] = _dual_node_list.size();
deba@490
   272
      StackLevel level;
deba@490
   273
      level.node_level = _dual_node_list.size();
deba@490
   274
      _dual_node_list.push_back(node);
deba@490
   275
      for (InArcIt it(*_digraph, node); it != INVALID; ++it) {
deba@490
   276
        Arc arc = it;
deba@490
   277
        Node source = _digraph->source(arc);
deba@490
   278
        Value value = (*_cost)[it];
deba@490
   279
        if (source == node || (*_node_order)[source] == -3) continue;
deba@490
   280
        if ((*_cost_arcs)[source].arc == INVALID) {
deba@490
   281
          (*_cost_arcs)[source].arc = arc;
deba@490
   282
          (*_cost_arcs)[source].value = value;
deba@490
   283
          nodes.push_back(source);
deba@490
   284
        } else {
deba@490
   285
          if ((*_cost_arcs)[source].value > value) {
deba@490
   286
            (*_cost_arcs)[source].arc = arc;
deba@490
   287
            (*_cost_arcs)[source].value = value;
deba@490
   288
          }
deba@490
   289
        }
deba@490
   290
      }
deba@490
   291
      CostArc minimum = (*_cost_arcs)[nodes[0]];
deba@490
   292
      for (int i = 1; i < int(nodes.size()); ++i) {
deba@490
   293
        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
deba@490
   294
          minimum = (*_cost_arcs)[nodes[i]];
deba@490
   295
        }
deba@490
   296
      }
kpeter@581
   297
      (*_arc_order)[minimum.arc] = _dual_variables.size();
deba@490
   298
      DualVariable var(_dual_node_list.size() - 1,
deba@490
   299
                       _dual_node_list.size(), minimum.value);
deba@490
   300
      _dual_variables.push_back(var);
deba@490
   301
      for (int i = 0; i < int(nodes.size()); ++i) {
deba@490
   302
        (*_cost_arcs)[nodes[i]].value -= minimum.value;
deba@490
   303
        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
deba@490
   304
        (*_cost_arcs)[nodes[i]].arc = INVALID;
deba@490
   305
      }
deba@490
   306
      level_stack.push_back(level);
deba@490
   307
      return minimum.arc;
deba@490
   308
    }
deba@490
   309
deba@490
   310
    Arc contract(Node node) {
deba@490
   311
      int node_bottom = bottom(node);
deba@490
   312
      std::vector<Node> nodes;
deba@490
   313
      while (!level_stack.empty() &&
deba@490
   314
             level_stack.back().node_level >= node_bottom) {
deba@490
   315
        for (int i = 0; i < int(level_stack.back().arcs.size()); ++i) {
deba@490
   316
          Arc arc = level_stack.back().arcs[i].arc;
deba@490
   317
          Node source = _digraph->source(arc);
deba@490
   318
          Value value = level_stack.back().arcs[i].value;
deba@490
   319
          if ((*_node_order)[source] >= node_bottom) continue;
deba@490
   320
          if ((*_cost_arcs)[source].arc == INVALID) {
deba@490
   321
            (*_cost_arcs)[source].arc = arc;
deba@490
   322
            (*_cost_arcs)[source].value = value;
deba@490
   323
            nodes.push_back(source);
deba@490
   324
          } else {
deba@490
   325
            if ((*_cost_arcs)[source].value > value) {
deba@490
   326
              (*_cost_arcs)[source].arc = arc;
deba@490
   327
              (*_cost_arcs)[source].value = value;
deba@490
   328
            }
deba@490
   329
          }
deba@490
   330
        }
deba@490
   331
        level_stack.pop_back();
deba@490
   332
      }
deba@490
   333
      CostArc minimum = (*_cost_arcs)[nodes[0]];
deba@490
   334
      for (int i = 1; i < int(nodes.size()); ++i) {
deba@490
   335
        if ((*_cost_arcs)[nodes[i]].value < minimum.value) {
deba@490
   336
          minimum = (*_cost_arcs)[nodes[i]];
deba@490
   337
        }
deba@490
   338
      }
kpeter@581
   339
      (*_arc_order)[minimum.arc] = _dual_variables.size();
deba@490
   340
      DualVariable var(node_bottom, _dual_node_list.size(), minimum.value);
deba@490
   341
      _dual_variables.push_back(var);
deba@490
   342
      StackLevel level;
deba@490
   343
      level.node_level = node_bottom;
deba@490
   344
      for (int i = 0; i < int(nodes.size()); ++i) {
deba@490
   345
        (*_cost_arcs)[nodes[i]].value -= minimum.value;
deba@490
   346
        level.arcs.push_back((*_cost_arcs)[nodes[i]]);
deba@490
   347
        (*_cost_arcs)[nodes[i]].arc = INVALID;
deba@490
   348
      }
deba@490
   349
      level_stack.push_back(level);
deba@490
   350
      return minimum.arc;
deba@490
   351
    }
deba@490
   352
deba@490
   353
    int bottom(Node node) {
deba@490
   354
      int k = level_stack.size() - 1;
deba@490
   355
      while (level_stack[k].node_level > (*_node_order)[node]) {
deba@490
   356
        --k;
deba@490
   357
      }
deba@490
   358
      return level_stack[k].node_level;
deba@490
   359
    }
deba@490
   360
deba@490
   361
    void finalize(Arc arc) {
deba@490
   362
      Node node = _digraph->target(arc);
deba@490
   363
      _heap->push(node, (*_arc_order)[arc]);
deba@490
   364
      _pred->set(node, arc);
deba@490
   365
      while (!_heap->empty()) {
deba@490
   366
        Node source = _heap->top();
deba@490
   367
        _heap->pop();
kpeter@581
   368
        (*_node_order)[source] = -1;
deba@490
   369
        for (OutArcIt it(*_digraph, source); it != INVALID; ++it) {
deba@490
   370
          if ((*_arc_order)[it] < 0) continue;
deba@490
   371
          Node target = _digraph->target(it);
deba@490
   372
          switch(_heap->state(target)) {
deba@490
   373
          case Heap::PRE_HEAP:
deba@490
   374
            _heap->push(target, (*_arc_order)[it]);
deba@490
   375
            _pred->set(target, it);
deba@490
   376
            break;
deba@490
   377
          case Heap::IN_HEAP:
deba@490
   378
            if ((*_arc_order)[it] < (*_heap)[target]) {
deba@490
   379
              _heap->decrease(target, (*_arc_order)[it]);
deba@490
   380
              _pred->set(target, it);
deba@490
   381
            }
deba@490
   382
            break;
deba@490
   383
          case Heap::POST_HEAP:
deba@490
   384
            break;
deba@490
   385
          }
deba@490
   386
        }
deba@490
   387
        _arborescence->set((*_pred)[source], true);
deba@490
   388
      }
deba@490
   389
    }
deba@490
   390
deba@490
   391
deba@490
   392
  public:
deba@490
   393
kpeter@584
   394
    /// \name Named Template Parameters
deba@490
   395
deba@490
   396
    /// @{
deba@490
   397
deba@490
   398
    template <class T>
kpeter@625
   399
    struct SetArborescenceMapTraits : public Traits {
deba@490
   400
      typedef T ArborescenceMap;
deba@490
   401
      static ArborescenceMap *createArborescenceMap(const Digraph &)
deba@490
   402
      {
deba@490
   403
        LEMON_ASSERT(false, "ArborescenceMap is not initialized");
deba@490
   404
        return 0; // ignore warnings
deba@490
   405
      }
deba@490
   406
    };
deba@490
   407
deba@490
   408
    /// \brief \ref named-templ-param "Named parameter" for
kpeter@625
   409
    /// setting \c ArborescenceMap type
deba@490
   410
    ///
deba@490
   411
    /// \ref named-templ-param "Named parameter" for setting
kpeter@625
   412
    /// \c ArborescenceMap type.
kpeter@625
   413
    /// It must conform to the \ref concepts::WriteMap "WriteMap" concept,
kpeter@625
   414
    /// and its value type must be \c bool (or convertible).
kpeter@625
   415
    /// Initially it will be set to \c false on each arc,
kpeter@625
   416
    /// then it will be set on each arborescence arc once.
deba@490
   417
    template <class T>
kpeter@625
   418
    struct SetArborescenceMap
deba@490
   419
      : public MinCostArborescence<Digraph, CostMap,
kpeter@625
   420
                                   SetArborescenceMapTraits<T> > {
deba@490
   421
    };
deba@490
   422
deba@490
   423
    template <class T>
kpeter@625
   424
    struct SetPredMapTraits : public Traits {
deba@490
   425
      typedef T PredMap;
deba@490
   426
      static PredMap *createPredMap(const Digraph &)
deba@490
   427
      {
deba@490
   428
        LEMON_ASSERT(false, "PredMap is not initialized");
kpeter@625
   429
        return 0; // ignore warnings
deba@490
   430
      }
deba@490
   431
    };
deba@490
   432
deba@490
   433
    /// \brief \ref named-templ-param "Named parameter" for
kpeter@625
   434
    /// setting \c PredMap type
deba@490
   435
    ///
deba@490
   436
    /// \ref named-templ-param "Named parameter" for setting
kpeter@625
   437
    /// \c PredMap type.
kpeter@625
   438
    /// It must meet the \ref concepts::WriteMap "WriteMap" concept, 
kpeter@625
   439
    /// and its value type must be the \c Arc type of the digraph.
deba@490
   440
    template <class T>
kpeter@625
   441
    struct SetPredMap
kpeter@625
   442
      : public MinCostArborescence<Digraph, CostMap, SetPredMapTraits<T> > {
deba@490
   443
    };
deba@490
   444
deba@490
   445
    /// @}
deba@490
   446
deba@490
   447
    /// \brief Constructor.
deba@490
   448
    ///
kpeter@559
   449
    /// \param digraph The digraph the algorithm will run on.
kpeter@559
   450
    /// \param cost The cost map used by the algorithm.
deba@490
   451
    MinCostArborescence(const Digraph& digraph, const CostMap& cost)
deba@490
   452
      : _digraph(&digraph), _cost(&cost), _pred(0), local_pred(false),
deba@490
   453
        _arborescence(0), local_arborescence(false),
deba@490
   454
        _arc_order(0), _node_order(0), _cost_arcs(0),
deba@490
   455
        _heap_cross_ref(0), _heap(0) {}
deba@490
   456
deba@490
   457
    /// \brief Destructor.
deba@490
   458
    ~MinCostArborescence() {
deba@490
   459
      destroyStructures();
deba@490
   460
    }
deba@490
   461
deba@490
   462
    /// \brief Sets the arborescence map.
deba@490
   463
    ///
deba@490
   464
    /// Sets the arborescence map.
kpeter@559
   465
    /// \return <tt>(*this)</tt>
deba@490
   466
    MinCostArborescence& arborescenceMap(ArborescenceMap& m) {
deba@490
   467
      if (local_arborescence) {
deba@490
   468
        delete _arborescence;
deba@490
   469
      }
deba@490
   470
      local_arborescence = false;
deba@490
   471
      _arborescence = &m;
deba@490
   472
      return *this;
deba@490
   473
    }
deba@490
   474
kpeter@625
   475
    /// \brief Sets the predecessor map.
deba@490
   476
    ///
kpeter@625
   477
    /// Sets the predecessor map.
kpeter@559
   478
    /// \return <tt>(*this)</tt>
deba@490
   479
    MinCostArborescence& predMap(PredMap& m) {
deba@490
   480
      if (local_pred) {
deba@490
   481
        delete _pred;
deba@490
   482
      }
deba@490
   483
      local_pred = false;
deba@490
   484
      _pred = &m;
deba@490
   485
      return *this;
deba@490
   486
    }
deba@490
   487
kpeter@584
   488
    /// \name Execution Control
deba@490
   489
    /// The simplest way to execute the algorithm is to use
deba@490
   490
    /// one of the member functions called \c run(...). \n
deba@490
   491
    /// If you need more control on the execution,
deba@490
   492
    /// first you must call \ref init(), then you can add several
deba@490
   493
    /// source nodes with \ref addSource().
deba@490
   494
    /// Finally \ref start() will perform the arborescence
deba@490
   495
    /// computation.
deba@490
   496
deba@490
   497
    ///@{
deba@490
   498
deba@490
   499
    /// \brief Initializes the internal data structures.
deba@490
   500
    ///
deba@490
   501
    /// Initializes the internal data structures.
deba@490
   502
    ///
deba@490
   503
    void init() {
deba@490
   504
      createStructures();
deba@490
   505
      _heap->clear();
deba@490
   506
      for (NodeIt it(*_digraph); it != INVALID; ++it) {
deba@490
   507
        (*_cost_arcs)[it].arc = INVALID;
kpeter@581
   508
        (*_node_order)[it] = -3;
kpeter@581
   509
        (*_heap_cross_ref)[it] = Heap::PRE_HEAP;
deba@490
   510
        _pred->set(it, INVALID);
deba@490
   511
      }
deba@490
   512
      for (ArcIt it(*_digraph); it != INVALID; ++it) {
deba@490
   513
        _arborescence->set(it, false);
kpeter@581
   514
        (*_arc_order)[it] = -1;
deba@490
   515
      }
deba@490
   516
      _dual_node_list.clear();
deba@490
   517
      _dual_variables.clear();
deba@490
   518
    }
deba@490
   519
deba@490
   520
    /// \brief Adds a new source node.
deba@490
   521
    ///
deba@490
   522
    /// Adds a new source node to the algorithm.
deba@490
   523
    void addSource(Node source) {
deba@490
   524
      std::vector<Node> nodes;
deba@490
   525
      nodes.push_back(source);
deba@490
   526
      while (!nodes.empty()) {
deba@490
   527
        Node node = nodes.back();
deba@490
   528
        nodes.pop_back();
deba@490
   529
        for (OutArcIt it(*_digraph, node); it != INVALID; ++it) {
deba@490
   530
          Node target = _digraph->target(it);
deba@490
   531
          if ((*_node_order)[target] == -3) {
deba@490
   532
            (*_node_order)[target] = -2;
deba@490
   533
            nodes.push_back(target);
deba@490
   534
            queue.push_back(target);
deba@490
   535
          }
deba@490
   536
        }
deba@490
   537
      }
deba@490
   538
      (*_node_order)[source] = -1;
deba@490
   539
    }
deba@490
   540
deba@490
   541
    /// \brief Processes the next node in the priority queue.
deba@490
   542
    ///
deba@490
   543
    /// Processes the next node in the priority queue.
deba@490
   544
    ///
deba@490
   545
    /// \return The processed node.
deba@490
   546
    ///
kpeter@625
   547
    /// \warning The queue must not be empty.
deba@490
   548
    Node processNextNode() {
deba@490
   549
      Node node = queue.back();
deba@490
   550
      queue.pop_back();
deba@490
   551
      if ((*_node_order)[node] == -2) {
deba@490
   552
        Arc arc = prepare(node);
deba@490
   553
        Node source = _digraph->source(arc);
deba@490
   554
        while ((*_node_order)[source] != -1) {
deba@490
   555
          if ((*_node_order)[source] >= 0) {
deba@490
   556
            arc = contract(source);
deba@490
   557
          } else {
deba@490
   558
            arc = prepare(source);
deba@490
   559
          }
deba@490
   560
          source = _digraph->source(arc);
deba@490
   561
        }
deba@490
   562
        finalize(arc);
deba@490
   563
        level_stack.clear();
deba@490
   564
      }
deba@490
   565
      return node;
deba@490
   566
    }
deba@490
   567
deba@490
   568
    /// \brief Returns the number of the nodes to be processed.
deba@490
   569
    ///
kpeter@625
   570
    /// Returns the number of the nodes to be processed in the priority
kpeter@625
   571
    /// queue.
deba@490
   572
    int queueSize() const {
deba@490
   573
      return queue.size();
deba@490
   574
    }
deba@490
   575
deba@490
   576
    /// \brief Returns \c false if there are nodes to be processed.
deba@490
   577
    ///
deba@490
   578
    /// Returns \c false if there are nodes to be processed.
deba@490
   579
    bool emptyQueue() const {
deba@490
   580
      return queue.empty();
deba@490
   581
    }
deba@490
   582
deba@490
   583
    /// \brief Executes the algorithm.
deba@490
   584
    ///
deba@490
   585
    /// Executes the algorithm.
deba@490
   586
    ///
deba@490
   587
    /// \pre init() must be called and at least one node should be added
deba@490
   588
    /// with addSource() before using this function.
deba@490
   589
    ///
deba@490
   590
    ///\note mca.start() is just a shortcut of the following code.
deba@490
   591
    ///\code
deba@490
   592
    ///while (!mca.emptyQueue()) {
deba@490
   593
    ///  mca.processNextNode();
deba@490
   594
    ///}
deba@490
   595
    ///\endcode
deba@490
   596
    void start() {
deba@490
   597
      while (!emptyQueue()) {
deba@490
   598
        processNextNode();
deba@490
   599
      }
deba@490
   600
    }
deba@490
   601
deba@490
   602
    /// \brief Runs %MinCostArborescence algorithm from node \c s.
deba@490
   603
    ///
deba@490
   604
    /// This method runs the %MinCostArborescence algorithm from
deba@490
   605
    /// a root node \c s.
deba@490
   606
    ///
deba@490
   607
    /// \note mca.run(s) is just a shortcut of the following code.
deba@490
   608
    /// \code
deba@490
   609
    /// mca.init();
deba@490
   610
    /// mca.addSource(s);
deba@490
   611
    /// mca.start();
deba@490
   612
    /// \endcode
kpeter@625
   613
    void run(Node s) {
deba@490
   614
      init();
kpeter@625
   615
      addSource(s);
deba@490
   616
      start();
deba@490
   617
    }
deba@490
   618
deba@490
   619
    ///@}
deba@490
   620
kpeter@625
   621
    /// \name Query Functions
kpeter@625
   622
    /// The result of the %MinCostArborescence algorithm can be obtained
kpeter@625
   623
    /// using these functions.\n
kpeter@625
   624
    /// Either run() or start() must be called before using them.
kpeter@625
   625
kpeter@625
   626
    /// @{
kpeter@625
   627
kpeter@625
   628
    /// \brief Returns the cost of the arborescence.
kpeter@625
   629
    ///
kpeter@625
   630
    /// Returns the cost of the arborescence.
kpeter@625
   631
    Value arborescenceCost() const {
kpeter@625
   632
      Value sum = 0;
kpeter@625
   633
      for (ArcIt it(*_digraph); it != INVALID; ++it) {
kpeter@625
   634
        if (arborescence(it)) {
kpeter@625
   635
          sum += (*_cost)[it];
kpeter@625
   636
        }
kpeter@625
   637
      }
kpeter@625
   638
      return sum;
kpeter@625
   639
    }
kpeter@625
   640
kpeter@625
   641
    /// \brief Returns \c true if the arc is in the arborescence.
kpeter@625
   642
    ///
kpeter@625
   643
    /// Returns \c true if the given arc is in the arborescence.
kpeter@625
   644
    /// \param arc An arc of the digraph.
kpeter@625
   645
    /// \pre \ref run() must be called before using this function.
kpeter@625
   646
    bool arborescence(Arc arc) const {
kpeter@625
   647
      return (*_pred)[_digraph->target(arc)] == arc;
kpeter@625
   648
    }
kpeter@625
   649
kpeter@625
   650
    /// \brief Returns a const reference to the arborescence map.
kpeter@625
   651
    ///
kpeter@625
   652
    /// Returns a const reference to the arborescence map.
kpeter@625
   653
    /// \pre \ref run() must be called before using this function.
kpeter@625
   654
    const ArborescenceMap& arborescenceMap() const {
kpeter@625
   655
      return *_arborescence;
kpeter@625
   656
    }
kpeter@625
   657
kpeter@625
   658
    /// \brief Returns the predecessor arc of the given node.
kpeter@625
   659
    ///
kpeter@625
   660
    /// Returns the predecessor arc of the given node.
kpeter@625
   661
    /// \pre \ref run() must be called before using this function.
kpeter@625
   662
    Arc pred(Node node) const {
kpeter@625
   663
      return (*_pred)[node];
kpeter@625
   664
    }
kpeter@625
   665
kpeter@625
   666
    /// \brief Returns a const reference to the pred map.
kpeter@625
   667
    ///
kpeter@625
   668
    /// Returns a const reference to the pred map.
kpeter@625
   669
    /// \pre \ref run() must be called before using this function.
kpeter@625
   670
    const PredMap& predMap() const {
kpeter@625
   671
      return *_pred;
kpeter@625
   672
    }
kpeter@625
   673
kpeter@625
   674
    /// \brief Indicates that a node is reachable from the sources.
kpeter@625
   675
    ///
kpeter@625
   676
    /// Indicates that a node is reachable from the sources.
kpeter@625
   677
    bool reached(Node node) const {
kpeter@625
   678
      return (*_node_order)[node] != -3;
kpeter@625
   679
    }
kpeter@625
   680
kpeter@625
   681
    /// \brief Indicates that a node is processed.
kpeter@625
   682
    ///
kpeter@625
   683
    /// Indicates that a node is processed. The arborescence path exists
kpeter@625
   684
    /// from the source to the given node.
kpeter@625
   685
    bool processed(Node node) const {
kpeter@625
   686
      return (*_node_order)[node] == -1;
kpeter@625
   687
    }
kpeter@625
   688
kpeter@625
   689
    /// \brief Returns the number of the dual variables in basis.
kpeter@625
   690
    ///
kpeter@625
   691
    /// Returns the number of the dual variables in basis.
kpeter@625
   692
    int dualNum() const {
kpeter@625
   693
      return _dual_variables.size();
kpeter@625
   694
    }
kpeter@625
   695
kpeter@625
   696
    /// \brief Returns the value of the dual solution.
kpeter@625
   697
    ///
kpeter@625
   698
    /// Returns the value of the dual solution. It should be
kpeter@625
   699
    /// equal to the arborescence value.
kpeter@625
   700
    Value dualValue() const {
kpeter@625
   701
      Value sum = 0;
kpeter@625
   702
      for (int i = 0; i < int(_dual_variables.size()); ++i) {
kpeter@625
   703
        sum += _dual_variables[i].value;
kpeter@625
   704
      }
kpeter@625
   705
      return sum;
kpeter@625
   706
    }
kpeter@625
   707
kpeter@625
   708
    /// \brief Returns the number of the nodes in the dual variable.
kpeter@625
   709
    ///
kpeter@625
   710
    /// Returns the number of the nodes in the dual variable.
kpeter@625
   711
    int dualSize(int k) const {
kpeter@625
   712
      return _dual_variables[k].end - _dual_variables[k].begin;
kpeter@625
   713
    }
kpeter@625
   714
kpeter@625
   715
    /// \brief Returns the value of the dual variable.
kpeter@625
   716
    ///
kpeter@625
   717
    /// Returns the the value of the dual variable.
kpeter@625
   718
    Value dualValue(int k) const {
kpeter@625
   719
      return _dual_variables[k].value;
kpeter@625
   720
    }
kpeter@625
   721
kpeter@625
   722
    /// \brief LEMON iterator for getting a dual variable.
kpeter@625
   723
    ///
kpeter@625
   724
    /// This class provides a common style LEMON iterator for getting a
kpeter@625
   725
    /// dual variable of \ref MinCostArborescence algorithm.
kpeter@625
   726
    /// It iterates over a subset of the nodes.
kpeter@625
   727
    class DualIt {
kpeter@625
   728
    public:
kpeter@625
   729
kpeter@625
   730
      /// \brief Constructor.
kpeter@625
   731
      ///
kpeter@625
   732
      /// Constructor for getting the nodeset of the dual variable
kpeter@625
   733
      /// of \ref MinCostArborescence algorithm.
kpeter@625
   734
      DualIt(const MinCostArborescence& algorithm, int variable)
kpeter@625
   735
        : _algorithm(&algorithm)
kpeter@625
   736
      {
kpeter@625
   737
        _index = _algorithm->_dual_variables[variable].begin;
kpeter@625
   738
        _last = _algorithm->_dual_variables[variable].end;
kpeter@625
   739
      }
kpeter@625
   740
kpeter@625
   741
      /// \brief Conversion to \c Node.
kpeter@625
   742
      ///
kpeter@625
   743
      /// Conversion to \c Node.
kpeter@625
   744
      operator Node() const {
kpeter@625
   745
        return _algorithm->_dual_node_list[_index];
kpeter@625
   746
      }
kpeter@625
   747
kpeter@625
   748
      /// \brief Increment operator.
kpeter@625
   749
      ///
kpeter@625
   750
      /// Increment operator.
kpeter@625
   751
      DualIt& operator++() {
kpeter@625
   752
        ++_index;
kpeter@625
   753
        return *this;
kpeter@625
   754
      }
kpeter@625
   755
kpeter@625
   756
      /// \brief Validity checking
kpeter@625
   757
      ///
kpeter@625
   758
      /// Checks whether the iterator is invalid.
kpeter@625
   759
      bool operator==(Invalid) const {
kpeter@625
   760
        return _index == _last;
kpeter@625
   761
      }
kpeter@625
   762
kpeter@625
   763
      /// \brief Validity checking
kpeter@625
   764
      ///
kpeter@625
   765
      /// Checks whether the iterator is valid.
kpeter@625
   766
      bool operator!=(Invalid) const {
kpeter@625
   767
        return _index != _last;
kpeter@625
   768
      }
kpeter@625
   769
kpeter@625
   770
    private:
kpeter@625
   771
      const MinCostArborescence* _algorithm;
kpeter@625
   772
      int _index, _last;
kpeter@625
   773
    };
kpeter@625
   774
kpeter@625
   775
    /// @}
kpeter@625
   776
deba@490
   777
  };
deba@490
   778
deba@490
   779
  /// \ingroup spantree
deba@490
   780
  ///
deba@490
   781
  /// \brief Function type interface for MinCostArborescence algorithm.
deba@490
   782
  ///
deba@490
   783
  /// Function type interface for MinCostArborescence algorithm.
kpeter@625
   784
  /// \param digraph The digraph the algorithm runs on.
kpeter@625
   785
  /// \param cost An arc map storing the costs.
kpeter@625
   786
  /// \param source The source node of the arborescence.
kpeter@625
   787
  /// \retval arborescence An arc map with \c bool (or convertible) value
kpeter@625
   788
  /// type that stores the arborescence.
kpeter@625
   789
  /// \return The total cost of the arborescence.
deba@490
   790
  ///
deba@490
   791
  /// \sa MinCostArborescence
deba@490
   792
  template <typename Digraph, typename CostMap, typename ArborescenceMap>
deba@490
   793
  typename CostMap::Value minCostArborescence(const Digraph& digraph,
deba@490
   794
                                              const CostMap& cost,
deba@490
   795
                                              typename Digraph::Node source,
deba@490
   796
                                              ArborescenceMap& arborescence) {
deba@490
   797
    typename MinCostArborescence<Digraph, CostMap>
kpeter@625
   798
      ::template SetArborescenceMap<ArborescenceMap>
deba@490
   799
      ::Create mca(digraph, cost);
deba@490
   800
    mca.arborescenceMap(arborescence);
deba@490
   801
    mca.run(source);
kpeter@625
   802
    return mca.arborescenceCost();
deba@490
   803
  }
deba@490
   804
deba@490
   805
}
deba@490
   806
deba@490
   807
#endif