lemon/min_cost_arborescence.h
author Peter Kovacs <kpeter@inf.elte.hu>
Wed, 25 Mar 2009 15:58:44 +0100
changeset 605 5232721b3f14
child 559 c5fd2d996909
permissions -rw-r--r--
Rework the interface of NetworkSimplex (#234)

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