gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Minor doc improvements
0 5 0
default
5 files changed with 8 insertions and 7 deletions:
↑ Collapse diff ↑
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2010
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_CAPACITY_SCALING_H
20 20
#define LEMON_CAPACITY_SCALING_H
21 21

	
22 22
/// \ingroup min_cost_flow_algs
23 23
///
24 24
/// \file
25 25
/// \brief Capacity Scaling algorithm for finding a minimum cost flow.
26 26

	
27 27
#include <vector>
28 28
#include <limits>
29 29
#include <lemon/core.h>
30 30
#include <lemon/bin_heap.h>
31 31

	
32 32
namespace lemon {
33 33

	
34 34
  /// \brief Default traits class of CapacityScaling algorithm.
35 35
  ///
36 36
  /// Default traits class of CapacityScaling algorithm.
37 37
  /// \tparam GR Digraph type.
38 38
  /// \tparam V The number type used for flow amounts, capacity bounds
39 39
  /// and supply values. By default it is \c int.
40 40
  /// \tparam C The number type used for costs and potentials.
41 41
  /// By default it is the same as \c V.
42 42
  template <typename GR, typename V = int, typename C = V>
43 43
  struct CapacityScalingDefaultTraits
44 44
  {
45 45
    /// The type of the digraph
46 46
    typedef GR Digraph;
47 47
    /// The type of the flow amounts, capacity bounds and supply values
48 48
    typedef V Value;
49 49
    /// The type of the arc costs
50 50
    typedef C Cost;
51 51

	
52 52
    /// \brief The type of the heap used for internal Dijkstra computations.
53 53
    ///
54 54
    /// The type of the heap used for internal Dijkstra computations.
55 55
    /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
56 56
    /// its priority type must be \c Cost and its cross reference type
57 57
    /// must be \ref RangeMap "RangeMap<int>".
58 58
    typedef BinHeap<Cost, RangeMap<int> > Heap;
59 59
  };
60 60

	
61 61
  /// \addtogroup min_cost_flow_algs
62 62
  /// @{
63 63

	
64 64
  /// \brief Implementation of the Capacity Scaling algorithm for
65 65
  /// finding a \ref min_cost_flow "minimum cost flow".
66 66
  ///
67 67
  /// \ref CapacityScaling implements the capacity scaling version
68 68
  /// of the successive shortest path algorithm for finding a
69 69
  /// \ref min_cost_flow "minimum cost flow" \ref amo93networkflows,
70 70
  /// \ref edmondskarp72theoretical. It is an efficient dual
71 71
  /// solution method.
72 72
  ///
73 73
  /// Most of the parameters of the problem (except for the digraph)
74 74
  /// can be given using separate functions, and the algorithm can be
75 75
  /// executed using the \ref run() function. If some parameters are not
76 76
  /// specified, then default values will be used.
77 77
  ///
78 78
  /// \tparam GR The digraph type the algorithm runs on.
79 79
  /// \tparam V The number type used for flow amounts, capacity bounds
80 80
  /// and supply values in the algorithm. By default, it is \c int.
81 81
  /// \tparam C The number type used for costs and potentials in the
82 82
  /// algorithm. By default, it is the same as \c V.
83 83
  /// \tparam TR The traits class that defines various types used by the
84 84
  /// algorithm. By default, it is \ref CapacityScalingDefaultTraits
85 85
  /// "CapacityScalingDefaultTraits<GR, V, C>".
86 86
  /// In most cases, this parameter should not be set directly,
87 87
  /// consider to use the named template parameters instead.
88 88
  ///
89
  /// \warning Both number types must be signed and all input data must
89
  /// \warning Both \c V and \c C must be signed number types.
90
  /// \warning All input data (capacities, supply values, and costs) must
90 91
  /// be integer.
91 92
  /// \warning This algorithm does not support negative costs for such
92 93
  /// arcs that have infinite upper bound.
93 94
#ifdef DOXYGEN
94 95
  template <typename GR, typename V, typename C, typename TR>
95 96
#else
96 97
  template < typename GR, typename V = int, typename C = V,
97 98
             typename TR = CapacityScalingDefaultTraits<GR, V, C> >
98 99
#endif
99 100
  class CapacityScaling
100 101
  {
101 102
  public:
102 103

	
103 104
    /// The type of the digraph
104 105
    typedef typename TR::Digraph Digraph;
105 106
    /// The type of the flow amounts, capacity bounds and supply values
106 107
    typedef typename TR::Value Value;
107 108
    /// The type of the arc costs
108 109
    typedef typename TR::Cost Cost;
109 110

	
110 111
    /// The type of the heap used for internal Dijkstra computations
111 112
    typedef typename TR::Heap Heap;
112 113

	
113 114
    /// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm
114 115
    typedef TR Traits;
115 116

	
116 117
  public:
117 118

	
118 119
    /// \brief Problem type constants for the \c run() function.
119 120
    ///
120 121
    /// Enum type containing the problem type constants that can be
121 122
    /// returned by the \ref run() function of the algorithm.
122 123
    enum ProblemType {
123 124
      /// The problem has no feasible solution (flow).
124 125
      INFEASIBLE,
125 126
      /// The problem has optimal solution (i.e. it is feasible and
126 127
      /// bounded), and the algorithm has found optimal flow and node
127 128
      /// potentials (primal and dual solutions).
128 129
      OPTIMAL,
129 130
      /// The digraph contains an arc of negative cost and infinite
130 131
      /// upper bound. It means that the objective function is unbounded
131 132
      /// on that arc, however, note that it could actually be bounded
132 133
      /// over the feasible flows, but this algroithm cannot handle
133 134
      /// these cases.
134 135
      UNBOUNDED
135 136
    };
136 137

	
137 138
  private:
138 139

	
139 140
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
140 141

	
141 142
    typedef std::vector<int> IntVector;
142 143
    typedef std::vector<Value> ValueVector;
143 144
    typedef std::vector<Cost> CostVector;
144 145
    typedef std::vector<char> BoolVector;
145 146
    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
146 147

	
147 148
  private:
148 149

	
149 150
    // Data related to the underlying digraph
150 151
    const GR &_graph;
151 152
    int _node_num;
152 153
    int _arc_num;
153 154
    int _res_arc_num;
154 155
    int _root;
155 156

	
156 157
    // Parameters of the problem
157 158
    bool _have_lower;
158 159
    Value _sum_supply;
159 160

	
160 161
    // Data structures for storing the digraph
161 162
    IntNodeMap _node_id;
162 163
    IntArcMap _arc_idf;
163 164
    IntArcMap _arc_idb;
164 165
    IntVector _first_out;
165 166
    BoolVector _forward;
166 167
    IntVector _source;
167 168
    IntVector _target;
168 169
    IntVector _reverse;
169 170

	
170 171
    // Node and arc data
171 172
    ValueVector _lower;
172 173
    ValueVector _upper;
173 174
    CostVector _cost;
174 175
    ValueVector _supply;
175 176

	
176 177
    ValueVector _res_cap;
177 178
    CostVector _pi;
178 179
    ValueVector _excess;
179 180
    IntVector _excess_nodes;
180 181
    IntVector _deficit_nodes;
181 182

	
182 183
    Value _delta;
183 184
    int _factor;
184 185
    IntVector _pred;
185 186

	
186 187
  public:
187 188

	
188 189
    /// \brief Constant for infinite upper bounds (capacities).
189 190
    ///
190 191
    /// Constant for infinite upper bounds (capacities).
191 192
    /// It is \c std::numeric_limits<Value>::infinity() if available,
192 193
    /// \c std::numeric_limits<Value>::max() otherwise.
193 194
    const Value INF;
194 195

	
195 196
  private:
196 197

	
197 198
    // Special implementation of the Dijkstra algorithm for finding
198 199
    // shortest paths in the residual network of the digraph with
199 200
    // respect to the reduced arc costs and modifying the node
200 201
    // potentials according to the found distance labels.
201 202
    class ResidualDijkstra
202 203
    {
203 204
    private:
204 205

	
205 206
      int _node_num;
206 207
      bool _geq;
207 208
      const IntVector &_first_out;
208 209
      const IntVector &_target;
209 210
      const CostVector &_cost;
210 211
      const ValueVector &_res_cap;
211 212
      const ValueVector &_excess;
212 213
      CostVector &_pi;
213 214
      IntVector &_pred;
214 215

	
215 216
      IntVector _proc_nodes;
216 217
      CostVector _dist;
217 218

	
218 219
    public:
219 220

	
220 221
      ResidualDijkstra(CapacityScaling& cs) :
221 222
        _node_num(cs._node_num), _geq(cs._sum_supply < 0),
222 223
        _first_out(cs._first_out), _target(cs._target), _cost(cs._cost),
223 224
        _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),
224 225
        _pred(cs._pred), _dist(cs._node_num)
225 226
      {}
226 227

	
227 228
      int run(int s, Value delta = 1) {
228 229
        RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP);
229 230
        Heap heap(heap_cross_ref);
230 231
        heap.push(s, 0);
231 232
        _pred[s] = -1;
232 233
        _proc_nodes.clear();
233 234

	
234 235
        // Process nodes
235 236
        while (!heap.empty() && _excess[heap.top()] > -delta) {
236 237
          int u = heap.top(), v;
237 238
          Cost d = heap.prio() + _pi[u], dn;
238 239
          _dist[u] = heap.prio();
239 240
          _proc_nodes.push_back(u);
240 241
          heap.pop();
241 242

	
242 243
          // Traverse outgoing residual arcs
243 244
          int last_out = _geq ? _first_out[u+1] : _first_out[u+1] - 1;
244 245
          for (int a = _first_out[u]; a != last_out; ++a) {
245 246
            if (_res_cap[a] < delta) continue;
246 247
            v = _target[a];
247 248
            switch (heap.state(v)) {
248 249
              case Heap::PRE_HEAP:
249 250
                heap.push(v, d + _cost[a] - _pi[v]);
250 251
                _pred[v] = a;
251 252
                break;
252 253
              case Heap::IN_HEAP:
253 254
                dn = d + _cost[a] - _pi[v];
254 255
                if (dn < heap[v]) {
255 256
                  heap.decrease(v, dn);
256 257
                  _pred[v] = a;
257 258
                }
258 259
                break;
259 260
              case Heap::POST_HEAP:
260 261
                break;
261 262
            }
262 263
          }
263 264
        }
264 265
        if (heap.empty()) return -1;
265 266

	
266 267
        // Update potentials of processed nodes
267 268
        int t = heap.top();
268 269
        Cost dt = heap.prio();
269 270
        for (int i = 0; i < int(_proc_nodes.size()); ++i) {
270 271
          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
271 272
        }
272 273

	
273 274
        return t;
274 275
      }
275 276

	
276 277
    }; //class ResidualDijkstra
277 278

	
278 279
  public:
279 280

	
280 281
    /// \name Named Template Parameters
281 282
    /// @{
282 283

	
283 284
    template <typename T>
284 285
    struct SetHeapTraits : public Traits {
285 286
      typedef T Heap;
286 287
    };
287 288

	
288 289
    /// \brief \ref named-templ-param "Named parameter" for setting
289 290
    /// \c Heap type.
290 291
    ///
291 292
    /// \ref named-templ-param "Named parameter" for setting \c Heap
292 293
    /// type, which is used for internal Dijkstra computations.
293 294
    /// It must conform to the \ref lemon::concepts::Heap "Heap" concept,
294 295
    /// its priority type must be \c Cost and its cross reference type
295 296
    /// must be \ref RangeMap "RangeMap<int>".
296 297
    template <typename T>
297 298
    struct SetHeap
298 299
      : public CapacityScaling<GR, V, C, SetHeapTraits<T> > {
299 300
      typedef  CapacityScaling<GR, V, C, SetHeapTraits<T> > Create;
300 301
    };
301 302

	
302 303
    /// @}
303 304

	
304 305
  protected:
305 306

	
306 307
    CapacityScaling() {}
307 308

	
308 309
  public:
309 310

	
310 311
    /// \brief Constructor.
311 312
    ///
312 313
    /// The constructor of the class.
313 314
    ///
314 315
    /// \param graph The digraph the algorithm runs on.
315 316
    CapacityScaling(const GR& graph) :
316 317
      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
317 318
      INF(std::numeric_limits<Value>::has_infinity ?
318 319
          std::numeric_limits<Value>::infinity() :
319 320
          std::numeric_limits<Value>::max())
320 321
    {
321 322
      // Check the number types
322 323
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
323 324
        "The flow type of CapacityScaling must be signed");
324 325
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
325 326
        "The cost type of CapacityScaling must be signed");
326 327

	
327 328
      // Reset data structures
328 329
      reset();
329 330
    }
330 331

	
331 332
    /// \name Parameters
332 333
    /// The parameters of the algorithm can be specified using these
333 334
    /// functions.
334 335

	
335 336
    /// @{
336 337

	
337 338
    /// \brief Set the lower bounds on the arcs.
338 339
    ///
339 340
    /// This function sets the lower bounds on the arcs.
340 341
    /// If it is not used before calling \ref run(), the lower bounds
341 342
    /// will be set to zero on all arcs.
342 343
    ///
343 344
    /// \param map An arc map storing the lower bounds.
344 345
    /// Its \c Value type must be convertible to the \c Value type
345 346
    /// of the algorithm.
346 347
    ///
347 348
    /// \return <tt>(*this)</tt>
348 349
    template <typename LowerMap>
349 350
    CapacityScaling& lowerMap(const LowerMap& map) {
350 351
      _have_lower = true;
351 352
      for (ArcIt a(_graph); a != INVALID; ++a) {
352 353
        _lower[_arc_idf[a]] = map[a];
353 354
        _lower[_arc_idb[a]] = map[a];
354 355
      }
355 356
      return *this;
356 357
    }
357 358

	
358 359
    /// \brief Set the upper bounds (capacities) on the arcs.
359 360
    ///
360 361
    /// This function sets the upper bounds (capacities) on the arcs.
361 362
    /// If it is not used before calling \ref run(), the upper bounds
362 363
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
363 364
    /// unbounded from above).
364 365
    ///
365 366
    /// \param map An arc map storing the upper bounds.
366 367
    /// Its \c Value type must be convertible to the \c Value type
367 368
    /// of the algorithm.
368 369
    ///
369 370
    /// \return <tt>(*this)</tt>
370 371
    template<typename UpperMap>
371 372
    CapacityScaling& upperMap(const UpperMap& map) {
372 373
      for (ArcIt a(_graph); a != INVALID; ++a) {
373 374
        _upper[_arc_idf[a]] = map[a];
374 375
      }
375 376
      return *this;
376 377
    }
377 378

	
378 379
    /// \brief Set the costs of the arcs.
379 380
    ///
380 381
    /// This function sets the costs of the arcs.
381 382
    /// If it is not used before calling \ref run(), the costs
382 383
    /// will be set to \c 1 on all arcs.
383 384
    ///
384 385
    /// \param map An arc map storing the costs.
385 386
    /// Its \c Value type must be convertible to the \c Cost type
386 387
    /// of the algorithm.
387 388
    ///
388 389
    /// \return <tt>(*this)</tt>
389 390
    template<typename CostMap>
390 391
    CapacityScaling& costMap(const CostMap& map) {
391 392
      for (ArcIt a(_graph); a != INVALID; ++a) {
392 393
        _cost[_arc_idf[a]] =  map[a];
393 394
        _cost[_arc_idb[a]] = -map[a];
394 395
      }
395 396
      return *this;
396 397
    }
397 398

	
398 399
    /// \brief Set the supply values of the nodes.
399 400
    ///
400 401
    /// This function sets the supply values of the nodes.
401 402
    /// If neither this function nor \ref stSupply() is used before
402 403
    /// calling \ref run(), the supply of each node will be set to zero.
403 404
    ///
404 405
    /// \param map A node map storing the supply values.
405 406
    /// Its \c Value type must be convertible to the \c Value type
406 407
    /// of the algorithm.
407 408
    ///
408 409
    /// \return <tt>(*this)</tt>
409 410
    template<typename SupplyMap>
410 411
    CapacityScaling& supplyMap(const SupplyMap& map) {
411 412
      for (NodeIt n(_graph); n != INVALID; ++n) {
412 413
        _supply[_node_id[n]] = map[n];
413 414
      }
414 415
      return *this;
415 416
    }
416 417

	
417 418
    /// \brief Set single source and target nodes and a supply value.
418 419
    ///
419 420
    /// This function sets a single source node and a single target node
420 421
    /// and the required flow value.
421 422
    /// If neither this function nor \ref supplyMap() is used before
422 423
    /// calling \ref run(), the supply of each node will be set to zero.
423 424
    ///
424 425
    /// Using this function has the same effect as using \ref supplyMap()
425 426
    /// with such a map in which \c k is assigned to \c s, \c -k is
426 427
    /// assigned to \c t and all other nodes have zero supply value.
427 428
    ///
428 429
    /// \param s The source node.
429 430
    /// \param t The target node.
430 431
    /// \param k The required amount of flow from node \c s to node \c t
431 432
    /// (i.e. the supply of \c s and the demand of \c t).
432 433
    ///
433 434
    /// \return <tt>(*this)</tt>
434 435
    CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
435 436
      for (int i = 0; i != _node_num; ++i) {
436 437
        _supply[i] = 0;
437 438
      }
438 439
      _supply[_node_id[s]] =  k;
439 440
      _supply[_node_id[t]] = -k;
440 441
      return *this;
441 442
    }
442 443

	
443 444
    /// @}
444 445

	
445 446
    /// \name Execution control
446 447
    /// The algorithm can be executed using \ref run().
447 448

	
448 449
    /// @{
449 450

	
450 451
    /// \brief Run the algorithm.
451 452
    ///
452 453
    /// This function runs the algorithm.
453 454
    /// The paramters can be specified using functions \ref lowerMap(),
454 455
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
455 456
    /// For example,
456 457
    /// \code
457 458
    ///   CapacityScaling<ListDigraph> cs(graph);
458 459
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
459 460
    ///     .supplyMap(sup).run();
460 461
    /// \endcode
461 462
    ///
462 463
    /// This function can be called more than once. All the given parameters
463 464
    /// are kept for the next call, unless \ref resetParams() or \ref reset()
464 465
    /// is used, thus only the modified parameters have to be set again.
465 466
    /// If the underlying digraph was also modified after the construction
466 467
    /// of the class (or the last \ref reset() call), then the \ref reset()
467 468
    /// function must be called.
468 469
    ///
469 470
    /// \param factor The capacity scaling factor. It must be larger than
470 471
    /// one to use scaling. If it is less or equal to one, then scaling
471 472
    /// will be disabled.
472 473
    ///
473 474
    /// \return \c INFEASIBLE if no feasible flow exists,
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2010
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_COST_SCALING_H
20 20
#define LEMON_COST_SCALING_H
21 21

	
22 22
/// \ingroup min_cost_flow_algs
23 23
/// \file
24 24
/// \brief Cost scaling algorithm for finding a minimum cost flow.
25 25

	
26 26
#include <vector>
27 27
#include <deque>
28 28
#include <limits>
29 29

	
30 30
#include <lemon/core.h>
31 31
#include <lemon/maps.h>
32 32
#include <lemon/math.h>
33 33
#include <lemon/static_graph.h>
34 34
#include <lemon/circulation.h>
35 35
#include <lemon/bellman_ford.h>
36 36

	
37 37
namespace lemon {
38 38

	
39 39
  /// \brief Default traits class of CostScaling algorithm.
40 40
  ///
41 41
  /// Default traits class of CostScaling algorithm.
42 42
  /// \tparam GR Digraph type.
43 43
  /// \tparam V The number type used for flow amounts, capacity bounds
44 44
  /// and supply values. By default it is \c int.
45 45
  /// \tparam C The number type used for costs and potentials.
46 46
  /// By default it is the same as \c V.
47 47
#ifdef DOXYGEN
48 48
  template <typename GR, typename V = int, typename C = V>
49 49
#else
50 50
  template < typename GR, typename V = int, typename C = V,
51 51
             bool integer = std::numeric_limits<C>::is_integer >
52 52
#endif
53 53
  struct CostScalingDefaultTraits
54 54
  {
55 55
    /// The type of the digraph
56 56
    typedef GR Digraph;
57 57
    /// The type of the flow amounts, capacity bounds and supply values
58 58
    typedef V Value;
59 59
    /// The type of the arc costs
60 60
    typedef C Cost;
61 61

	
62 62
    /// \brief The large cost type used for internal computations
63 63
    ///
64 64
    /// The large cost type used for internal computations.
65 65
    /// It is \c long \c long if the \c Cost type is integer,
66 66
    /// otherwise it is \c double.
67 67
    /// \c Cost must be convertible to \c LargeCost.
68 68
    typedef double LargeCost;
69 69
  };
70 70

	
71 71
  // Default traits class for integer cost types
72 72
  template <typename GR, typename V, typename C>
73 73
  struct CostScalingDefaultTraits<GR, V, C, true>
74 74
  {
75 75
    typedef GR Digraph;
76 76
    typedef V Value;
77 77
    typedef C Cost;
78 78
#ifdef LEMON_HAVE_LONG_LONG
79 79
    typedef long long LargeCost;
80 80
#else
81 81
    typedef long LargeCost;
82 82
#endif
83 83
  };
84 84

	
85 85

	
86 86
  /// \addtogroup min_cost_flow_algs
87 87
  /// @{
88 88

	
89 89
  /// \brief Implementation of the Cost Scaling algorithm for
90 90
  /// finding a \ref min_cost_flow "minimum cost flow".
91 91
  ///
92 92
  /// \ref CostScaling implements a cost scaling algorithm that performs
93 93
  /// push/augment and relabel operations for finding a \ref min_cost_flow
94 94
  /// "minimum cost flow" \ref amo93networkflows, \ref goldberg90approximation,
95 95
  /// \ref goldberg97efficient, \ref bunnagel98efficient.
96 96
  /// It is a highly efficient primal-dual solution method, which
97 97
  /// can be viewed as the generalization of the \ref Preflow
98 98
  /// "preflow push-relabel" algorithm for the maximum flow problem.
99 99
  ///
100 100
  /// Most of the parameters of the problem (except for the digraph)
101 101
  /// can be given using separate functions, and the algorithm can be
102 102
  /// executed using the \ref run() function. If some parameters are not
103 103
  /// specified, then default values will be used.
104 104
  ///
105 105
  /// \tparam GR The digraph type the algorithm runs on.
106 106
  /// \tparam V The number type used for flow amounts, capacity bounds
107 107
  /// and supply values in the algorithm. By default, it is \c int.
108 108
  /// \tparam C The number type used for costs and potentials in the
109 109
  /// algorithm. By default, it is the same as \c V.
110 110
  /// \tparam TR The traits class that defines various types used by the
111 111
  /// algorithm. By default, it is \ref CostScalingDefaultTraits
112 112
  /// "CostScalingDefaultTraits<GR, V, C>".
113 113
  /// In most cases, this parameter should not be set directly,
114 114
  /// consider to use the named template parameters instead.
115 115
  ///
116
  /// \warning Both number types must be signed and all input data must
116
  /// \warning Both \c V and \c C must be signed number types.
117
  /// \warning All input data (capacities, supply values, and costs) must
117 118
  /// be integer.
118 119
  /// \warning This algorithm does not support negative costs for such
119 120
  /// arcs that have infinite upper bound.
120 121
  ///
121 122
  /// \note %CostScaling provides three different internal methods,
122 123
  /// from which the most efficient one is used by default.
123 124
  /// For more information, see \ref Method.
124 125
#ifdef DOXYGEN
125 126
  template <typename GR, typename V, typename C, typename TR>
126 127
#else
127 128
  template < typename GR, typename V = int, typename C = V,
128 129
             typename TR = CostScalingDefaultTraits<GR, V, C> >
129 130
#endif
130 131
  class CostScaling
131 132
  {
132 133
  public:
133 134

	
134 135
    /// The type of the digraph
135 136
    typedef typename TR::Digraph Digraph;
136 137
    /// The type of the flow amounts, capacity bounds and supply values
137 138
    typedef typename TR::Value Value;
138 139
    /// The type of the arc costs
139 140
    typedef typename TR::Cost Cost;
140 141

	
141 142
    /// \brief The large cost type
142 143
    ///
143 144
    /// The large cost type used for internal computations.
144 145
    /// By default, it is \c long \c long if the \c Cost type is integer,
145 146
    /// otherwise it is \c double.
146 147
    typedef typename TR::LargeCost LargeCost;
147 148

	
148 149
    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
149 150
    typedef TR Traits;
150 151

	
151 152
  public:
152 153

	
153 154
    /// \brief Problem type constants for the \c run() function.
154 155
    ///
155 156
    /// Enum type containing the problem type constants that can be
156 157
    /// returned by the \ref run() function of the algorithm.
157 158
    enum ProblemType {
158 159
      /// The problem has no feasible solution (flow).
159 160
      INFEASIBLE,
160 161
      /// The problem has optimal solution (i.e. it is feasible and
161 162
      /// bounded), and the algorithm has found optimal flow and node
162 163
      /// potentials (primal and dual solutions).
163 164
      OPTIMAL,
164 165
      /// The digraph contains an arc of negative cost and infinite
165 166
      /// upper bound. It means that the objective function is unbounded
166 167
      /// on that arc, however, note that it could actually be bounded
167 168
      /// over the feasible flows, but this algroithm cannot handle
168 169
      /// these cases.
169 170
      UNBOUNDED
170 171
    };
171 172

	
172 173
    /// \brief Constants for selecting the internal method.
173 174
    ///
174 175
    /// Enum type containing constants for selecting the internal method
175 176
    /// for the \ref run() function.
176 177
    ///
177 178
    /// \ref CostScaling provides three internal methods that differ mainly
178 179
    /// in their base operations, which are used in conjunction with the
179 180
    /// relabel operation.
180 181
    /// By default, the so called \ref PARTIAL_AUGMENT
181 182
    /// "Partial Augment-Relabel" method is used, which proved to be
182 183
    /// the most efficient and the most robust on various test inputs.
183 184
    /// However, the other methods can be selected using the \ref run()
184 185
    /// function with the proper parameter.
185 186
    enum Method {
186 187
      /// Local push operations are used, i.e. flow is moved only on one
187 188
      /// admissible arc at once.
188 189
      PUSH,
189 190
      /// Augment operations are used, i.e. flow is moved on admissible
190 191
      /// paths from a node with excess to a node with deficit.
191 192
      AUGMENT,
192 193
      /// Partial augment operations are used, i.e. flow is moved on
193 194
      /// admissible paths started from a node with excess, but the
194 195
      /// lengths of these paths are limited. This method can be viewed
195 196
      /// as a combined version of the previous two operations.
196 197
      PARTIAL_AUGMENT
197 198
    };
198 199

	
199 200
  private:
200 201

	
201 202
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
202 203

	
203 204
    typedef std::vector<int> IntVector;
204 205
    typedef std::vector<Value> ValueVector;
205 206
    typedef std::vector<Cost> CostVector;
206 207
    typedef std::vector<LargeCost> LargeCostVector;
207 208
    typedef std::vector<char> BoolVector;
208 209
    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
209 210

	
210 211
  private:
211 212

	
212 213
    template <typename KT, typename VT>
213 214
    class StaticVectorMap {
214 215
    public:
215 216
      typedef KT Key;
216 217
      typedef VT Value;
217 218

	
218 219
      StaticVectorMap(std::vector<Value>& v) : _v(v) {}
219 220

	
220 221
      const Value& operator[](const Key& key) const {
221 222
        return _v[StaticDigraph::id(key)];
222 223
      }
223 224

	
224 225
      Value& operator[](const Key& key) {
225 226
        return _v[StaticDigraph::id(key)];
226 227
      }
227 228

	
228 229
      void set(const Key& key, const Value& val) {
229 230
        _v[StaticDigraph::id(key)] = val;
230 231
      }
231 232

	
232 233
    private:
233 234
      std::vector<Value>& _v;
234 235
    };
235 236

	
236 237
    typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
237 238
    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
238 239

	
239 240
  private:
240 241

	
241 242
    // Data related to the underlying digraph
242 243
    const GR &_graph;
243 244
    int _node_num;
244 245
    int _arc_num;
245 246
    int _res_node_num;
246 247
    int _res_arc_num;
247 248
    int _root;
248 249

	
249 250
    // Parameters of the problem
250 251
    bool _have_lower;
251 252
    Value _sum_supply;
252 253
    int _sup_node_num;
253 254

	
254 255
    // Data structures for storing the digraph
255 256
    IntNodeMap _node_id;
256 257
    IntArcMap _arc_idf;
257 258
    IntArcMap _arc_idb;
258 259
    IntVector _first_out;
259 260
    BoolVector _forward;
260 261
    IntVector _source;
261 262
    IntVector _target;
262 263
    IntVector _reverse;
263 264

	
264 265
    // Node and arc data
265 266
    ValueVector _lower;
266 267
    ValueVector _upper;
267 268
    CostVector _scost;
268 269
    ValueVector _supply;
269 270

	
270 271
    ValueVector _res_cap;
271 272
    LargeCostVector _cost;
272 273
    LargeCostVector _pi;
273 274
    ValueVector _excess;
274 275
    IntVector _next_out;
275 276
    std::deque<int> _active_nodes;
276 277

	
277 278
    // Data for scaling
278 279
    LargeCost _epsilon;
279 280
    int _alpha;
280 281

	
281 282
    IntVector _buckets;
282 283
    IntVector _bucket_next;
283 284
    IntVector _bucket_prev;
284 285
    IntVector _rank;
285 286
    int _max_rank;
286 287

	
287 288
    // Data for a StaticDigraph structure
288 289
    typedef std::pair<int, int> IntPair;
289 290
    StaticDigraph _sgr;
290 291
    std::vector<IntPair> _arc_vec;
291 292
    std::vector<LargeCost> _cost_vec;
292 293
    LargeCostArcMap _cost_map;
293 294
    LargeCostNodeMap _pi_map;
294 295

	
295 296
  public:
296 297

	
297 298
    /// \brief Constant for infinite upper bounds (capacities).
298 299
    ///
299 300
    /// Constant for infinite upper bounds (capacities).
300 301
    /// It is \c std::numeric_limits<Value>::infinity() if available,
301 302
    /// \c std::numeric_limits<Value>::max() otherwise.
302 303
    const Value INF;
303 304

	
304 305
  public:
305 306

	
306 307
    /// \name Named Template Parameters
307 308
    /// @{
308 309

	
309 310
    template <typename T>
310 311
    struct SetLargeCostTraits : public Traits {
311 312
      typedef T LargeCost;
312 313
    };
313 314

	
314 315
    /// \brief \ref named-templ-param "Named parameter" for setting
315 316
    /// \c LargeCost type.
316 317
    ///
317 318
    /// \ref named-templ-param "Named parameter" for setting \c LargeCost
318 319
    /// type, which is used for internal computations in the algorithm.
319 320
    /// \c Cost must be convertible to \c LargeCost.
320 321
    template <typename T>
321 322
    struct SetLargeCost
322 323
      : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
323 324
      typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
324 325
    };
325 326

	
326 327
    /// @}
327 328

	
328 329
  protected:
329 330

	
330 331
    CostScaling() {}
331 332

	
332 333
  public:
333 334

	
334 335
    /// \brief Constructor.
335 336
    ///
336 337
    /// The constructor of the class.
337 338
    ///
338 339
    /// \param graph The digraph the algorithm runs on.
339 340
    CostScaling(const GR& graph) :
340 341
      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
341 342
      _cost_map(_cost_vec), _pi_map(_pi),
342 343
      INF(std::numeric_limits<Value>::has_infinity ?
343 344
          std::numeric_limits<Value>::infinity() :
344 345
          std::numeric_limits<Value>::max())
345 346
    {
346 347
      // Check the number types
347 348
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
348 349
        "The flow type of CostScaling must be signed");
349 350
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
350 351
        "The cost type of CostScaling must be signed");
351 352

	
352 353
      // Reset data structures
353 354
      reset();
354 355
    }
355 356

	
356 357
    /// \name Parameters
357 358
    /// The parameters of the algorithm can be specified using these
358 359
    /// functions.
359 360

	
360 361
    /// @{
361 362

	
362 363
    /// \brief Set the lower bounds on the arcs.
363 364
    ///
364 365
    /// This function sets the lower bounds on the arcs.
365 366
    /// If it is not used before calling \ref run(), the lower bounds
366 367
    /// will be set to zero on all arcs.
367 368
    ///
368 369
    /// \param map An arc map storing the lower bounds.
369 370
    /// Its \c Value type must be convertible to the \c Value type
370 371
    /// of the algorithm.
371 372
    ///
372 373
    /// \return <tt>(*this)</tt>
373 374
    template <typename LowerMap>
374 375
    CostScaling& lowerMap(const LowerMap& map) {
375 376
      _have_lower = true;
376 377
      for (ArcIt a(_graph); a != INVALID; ++a) {
377 378
        _lower[_arc_idf[a]] = map[a];
378 379
        _lower[_arc_idb[a]] = map[a];
379 380
      }
380 381
      return *this;
381 382
    }
382 383

	
383 384
    /// \brief Set the upper bounds (capacities) on the arcs.
384 385
    ///
385 386
    /// This function sets the upper bounds (capacities) on the arcs.
386 387
    /// If it is not used before calling \ref run(), the upper bounds
387 388
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
388 389
    /// unbounded from above).
389 390
    ///
390 391
    /// \param map An arc map storing the upper bounds.
391 392
    /// Its \c Value type must be convertible to the \c Value type
392 393
    /// of the algorithm.
393 394
    ///
394 395
    /// \return <tt>(*this)</tt>
395 396
    template<typename UpperMap>
396 397
    CostScaling& upperMap(const UpperMap& map) {
397 398
      for (ArcIt a(_graph); a != INVALID; ++a) {
398 399
        _upper[_arc_idf[a]] = map[a];
399 400
      }
400 401
      return *this;
401 402
    }
402 403

	
403 404
    /// \brief Set the costs of the arcs.
404 405
    ///
405 406
    /// This function sets the costs of the arcs.
406 407
    /// If it is not used before calling \ref run(), the costs
407 408
    /// will be set to \c 1 on all arcs.
408 409
    ///
409 410
    /// \param map An arc map storing the costs.
410 411
    /// Its \c Value type must be convertible to the \c Cost type
411 412
    /// of the algorithm.
412 413
    ///
413 414
    /// \return <tt>(*this)</tt>
414 415
    template<typename CostMap>
415 416
    CostScaling& costMap(const CostMap& map) {
416 417
      for (ArcIt a(_graph); a != INVALID; ++a) {
417 418
        _scost[_arc_idf[a]] =  map[a];
418 419
        _scost[_arc_idb[a]] = -map[a];
419 420
      }
420 421
      return *this;
421 422
    }
422 423

	
423 424
    /// \brief Set the supply values of the nodes.
424 425
    ///
425 426
    /// This function sets the supply values of the nodes.
426 427
    /// If neither this function nor \ref stSupply() is used before
427 428
    /// calling \ref run(), the supply of each node will be set to zero.
428 429
    ///
429 430
    /// \param map A node map storing the supply values.
430 431
    /// Its \c Value type must be convertible to the \c Value type
431 432
    /// of the algorithm.
432 433
    ///
433 434
    /// \return <tt>(*this)</tt>
434 435
    template<typename SupplyMap>
435 436
    CostScaling& supplyMap(const SupplyMap& map) {
436 437
      for (NodeIt n(_graph); n != INVALID; ++n) {
437 438
        _supply[_node_id[n]] = map[n];
438 439
      }
439 440
      return *this;
440 441
    }
441 442

	
442 443
    /// \brief Set single source and target nodes and a supply value.
443 444
    ///
444 445
    /// This function sets a single source node and a single target node
445 446
    /// and the required flow value.
446 447
    /// If neither this function nor \ref supplyMap() is used before
447 448
    /// calling \ref run(), the supply of each node will be set to zero.
448 449
    ///
449 450
    /// Using this function has the same effect as using \ref supplyMap()
450 451
    /// with such a map in which \c k is assigned to \c s, \c -k is
451 452
    /// assigned to \c t and all other nodes have zero supply value.
452 453
    ///
453 454
    /// \param s The source node.
454 455
    /// \param t The target node.
455 456
    /// \param k The required amount of flow from node \c s to node \c t
456 457
    /// (i.e. the supply of \c s and the demand of \c t).
457 458
    ///
458 459
    /// \return <tt>(*this)</tt>
459 460
    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
460 461
      for (int i = 0; i != _res_node_num; ++i) {
461 462
        _supply[i] = 0;
462 463
      }
463 464
      _supply[_node_id[s]] =  k;
464 465
      _supply[_node_id[t]] = -k;
465 466
      return *this;
466 467
    }
467 468

	
468 469
    /// @}
469 470

	
470 471
    /// \name Execution control
471 472
    /// The algorithm can be executed using \ref run().
472 473

	
473 474
    /// @{
474 475

	
475 476
    /// \brief Run the algorithm.
476 477
    ///
477 478
    /// This function runs the algorithm.
478 479
    /// The paramters can be specified using functions \ref lowerMap(),
479 480
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
480 481
    /// For example,
481 482
    /// \code
482 483
    ///   CostScaling<ListDigraph> cs(graph);
483 484
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
484 485
    ///     .supplyMap(sup).run();
485 486
    /// \endcode
486 487
    ///
487 488
    /// This function can be called more than once. All the given parameters
488 489
    /// are kept for the next call, unless \ref resetParams() or \ref reset()
489 490
    /// is used, thus only the modified parameters have to be set again.
490 491
    /// If the underlying digraph was also modified after the construction
491 492
    /// of the class (or the last \ref reset() call), then the \ref reset()
492 493
    /// function must be called.
493 494
    ///
494 495
    /// \param method The internal method that will be used in the
495 496
    /// algorithm. For more information, see \ref Method.
496 497
    /// \param factor The cost scaling factor. It must be larger than one.
497 498
    ///
498 499
    /// \return \c INFEASIBLE if no feasible flow exists,
499 500
    /// \n \c OPTIMAL if the problem has optimal solution
500 501
    /// (i.e. it is feasible and bounded), and the algorithm has found
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2010
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_CYCLE_CANCELING_H
20 20
#define LEMON_CYCLE_CANCELING_H
21 21

	
22 22
/// \ingroup min_cost_flow_algs
23 23
/// \file
24 24
/// \brief Cycle-canceling algorithms for finding a minimum cost flow.
25 25

	
26 26
#include <vector>
27 27
#include <limits>
28 28

	
29 29
#include <lemon/core.h>
30 30
#include <lemon/maps.h>
31 31
#include <lemon/path.h>
32 32
#include <lemon/math.h>
33 33
#include <lemon/static_graph.h>
34 34
#include <lemon/adaptors.h>
35 35
#include <lemon/circulation.h>
36 36
#include <lemon/bellman_ford.h>
37 37
#include <lemon/howard_mmc.h>
38 38

	
39 39
namespace lemon {
40 40

	
41 41
  /// \addtogroup min_cost_flow_algs
42 42
  /// @{
43 43

	
44 44
  /// \brief Implementation of cycle-canceling algorithms for
45 45
  /// finding a \ref min_cost_flow "minimum cost flow".
46 46
  ///
47 47
  /// \ref CycleCanceling implements three different cycle-canceling
48 48
  /// algorithms for finding a \ref min_cost_flow "minimum cost flow"
49 49
  /// \ref amo93networkflows, \ref klein67primal,
50 50
  /// \ref goldberg89cyclecanceling.
51 51
  /// The most efficent one (both theoretically and practically)
52 52
  /// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm,
53 53
  /// thus it is the default method.
54 54
  /// It is strongly polynomial, but in practice, it is typically much
55 55
  /// slower than the scaling algorithms and NetworkSimplex.
56 56
  ///
57 57
  /// Most of the parameters of the problem (except for the digraph)
58 58
  /// can be given using separate functions, and the algorithm can be
59 59
  /// executed using the \ref run() function. If some parameters are not
60 60
  /// specified, then default values will be used.
61 61
  ///
62 62
  /// \tparam GR The digraph type the algorithm runs on.
63 63
  /// \tparam V The number type used for flow amounts, capacity bounds
64 64
  /// and supply values in the algorithm. By default, it is \c int.
65 65
  /// \tparam C The number type used for costs and potentials in the
66 66
  /// algorithm. By default, it is the same as \c V.
67 67
  ///
68
  /// \warning Both number types must be signed and all input data must
68
  /// \warning Both \c V and \c C must be signed number types.
69
  /// \warning All input data (capacities, supply values, and costs) must
69 70
  /// be integer.
70 71
  /// \warning This algorithm does not support negative costs for such
71 72
  /// arcs that have infinite upper bound.
72 73
  ///
73 74
  /// \note For more information about the three available methods,
74 75
  /// see \ref Method.
75 76
#ifdef DOXYGEN
76 77
  template <typename GR, typename V, typename C>
77 78
#else
78 79
  template <typename GR, typename V = int, typename C = V>
79 80
#endif
80 81
  class CycleCanceling
81 82
  {
82 83
  public:
83 84

	
84 85
    /// The type of the digraph
85 86
    typedef GR Digraph;
86 87
    /// The type of the flow amounts, capacity bounds and supply values
87 88
    typedef V Value;
88 89
    /// The type of the arc costs
89 90
    typedef C Cost;
90 91

	
91 92
  public:
92 93

	
93 94
    /// \brief Problem type constants for the \c run() function.
94 95
    ///
95 96
    /// Enum type containing the problem type constants that can be
96 97
    /// returned by the \ref run() function of the algorithm.
97 98
    enum ProblemType {
98 99
      /// The problem has no feasible solution (flow).
99 100
      INFEASIBLE,
100 101
      /// The problem has optimal solution (i.e. it is feasible and
101 102
      /// bounded), and the algorithm has found optimal flow and node
102 103
      /// potentials (primal and dual solutions).
103 104
      OPTIMAL,
104 105
      /// The digraph contains an arc of negative cost and infinite
105 106
      /// upper bound. It means that the objective function is unbounded
106 107
      /// on that arc, however, note that it could actually be bounded
107 108
      /// over the feasible flows, but this algroithm cannot handle
108 109
      /// these cases.
109 110
      UNBOUNDED
110 111
    };
111 112

	
112 113
    /// \brief Constants for selecting the used method.
113 114
    ///
114 115
    /// Enum type containing constants for selecting the used method
115 116
    /// for the \ref run() function.
116 117
    ///
117 118
    /// \ref CycleCanceling provides three different cycle-canceling
118 119
    /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten"
119 120
    /// is used, which proved to be the most efficient and the most robust
120 121
    /// on various test inputs.
121 122
    /// However, the other methods can be selected using the \ref run()
122 123
    /// function with the proper parameter.
123 124
    enum Method {
124 125
      /// A simple cycle-canceling method, which uses the
125 126
      /// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration
126 127
      /// number for detecting negative cycles in the residual network.
127 128
      SIMPLE_CYCLE_CANCELING,
128 129
      /// The "Minimum Mean Cycle-Canceling" algorithm, which is a
129 130
      /// well-known strongly polynomial method
130 131
      /// \ref goldberg89cyclecanceling. It improves along a
131 132
      /// \ref min_mean_cycle "minimum mean cycle" in each iteration.
132 133
      /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)).
133 134
      MINIMUM_MEAN_CYCLE_CANCELING,
134 135
      /// The "Cancel And Tighten" algorithm, which can be viewed as an
135 136
      /// improved version of the previous method
136 137
      /// \ref goldberg89cyclecanceling.
137 138
      /// It is faster both in theory and in practice, its running time
138 139
      /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)).
139 140
      CANCEL_AND_TIGHTEN
140 141
    };
141 142

	
142 143
  private:
143 144

	
144 145
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
145 146

	
146 147
    typedef std::vector<int> IntVector;
147 148
    typedef std::vector<double> DoubleVector;
148 149
    typedef std::vector<Value> ValueVector;
149 150
    typedef std::vector<Cost> CostVector;
150 151
    typedef std::vector<char> BoolVector;
151 152
    // Note: vector<char> is used instead of vector<bool> for efficiency reasons
152 153

	
153 154
  private:
154 155

	
155 156
    template <typename KT, typename VT>
156 157
    class StaticVectorMap {
157 158
    public:
158 159
      typedef KT Key;
159 160
      typedef VT Value;
160 161

	
161 162
      StaticVectorMap(std::vector<Value>& v) : _v(v) {}
162 163

	
163 164
      const Value& operator[](const Key& key) const {
164 165
        return _v[StaticDigraph::id(key)];
165 166
      }
166 167

	
167 168
      Value& operator[](const Key& key) {
168 169
        return _v[StaticDigraph::id(key)];
169 170
      }
170 171

	
171 172
      void set(const Key& key, const Value& val) {
172 173
        _v[StaticDigraph::id(key)] = val;
173 174
      }
174 175

	
175 176
    private:
176 177
      std::vector<Value>& _v;
177 178
    };
178 179

	
179 180
    typedef StaticVectorMap<StaticDigraph::Node, Cost> CostNodeMap;
180 181
    typedef StaticVectorMap<StaticDigraph::Arc, Cost> CostArcMap;
181 182

	
182 183
  private:
183 184

	
184 185

	
185 186
    // Data related to the underlying digraph
186 187
    const GR &_graph;
187 188
    int _node_num;
188 189
    int _arc_num;
189 190
    int _res_node_num;
190 191
    int _res_arc_num;
191 192
    int _root;
192 193

	
193 194
    // Parameters of the problem
194 195
    bool _have_lower;
195 196
    Value _sum_supply;
196 197

	
197 198
    // Data structures for storing the digraph
198 199
    IntNodeMap _node_id;
199 200
    IntArcMap _arc_idf;
200 201
    IntArcMap _arc_idb;
201 202
    IntVector _first_out;
202 203
    BoolVector _forward;
203 204
    IntVector _source;
204 205
    IntVector _target;
205 206
    IntVector _reverse;
206 207

	
207 208
    // Node and arc data
208 209
    ValueVector _lower;
209 210
    ValueVector _upper;
210 211
    CostVector _cost;
211 212
    ValueVector _supply;
212 213

	
213 214
    ValueVector _res_cap;
214 215
    CostVector _pi;
215 216

	
216 217
    // Data for a StaticDigraph structure
217 218
    typedef std::pair<int, int> IntPair;
218 219
    StaticDigraph _sgr;
219 220
    std::vector<IntPair> _arc_vec;
220 221
    std::vector<Cost> _cost_vec;
221 222
    IntVector _id_vec;
222 223
    CostArcMap _cost_map;
223 224
    CostNodeMap _pi_map;
224 225

	
225 226
  public:
226 227

	
227 228
    /// \brief Constant for infinite upper bounds (capacities).
228 229
    ///
229 230
    /// Constant for infinite upper bounds (capacities).
230 231
    /// It is \c std::numeric_limits<Value>::infinity() if available,
231 232
    /// \c std::numeric_limits<Value>::max() otherwise.
232 233
    const Value INF;
233 234

	
234 235
  public:
235 236

	
236 237
    /// \brief Constructor.
237 238
    ///
238 239
    /// The constructor of the class.
239 240
    ///
240 241
    /// \param graph The digraph the algorithm runs on.
241 242
    CycleCanceling(const GR& graph) :
242 243
      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
243 244
      _cost_map(_cost_vec), _pi_map(_pi),
244 245
      INF(std::numeric_limits<Value>::has_infinity ?
245 246
          std::numeric_limits<Value>::infinity() :
246 247
          std::numeric_limits<Value>::max())
247 248
    {
248 249
      // Check the number types
249 250
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
250 251
        "The flow type of CycleCanceling must be signed");
251 252
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
252 253
        "The cost type of CycleCanceling must be signed");
253 254

	
254 255
      // Reset data structures
255 256
      reset();
256 257
    }
257 258

	
258 259
    /// \name Parameters
259 260
    /// The parameters of the algorithm can be specified using these
260 261
    /// functions.
261 262

	
262 263
    /// @{
263 264

	
264 265
    /// \brief Set the lower bounds on the arcs.
265 266
    ///
266 267
    /// This function sets the lower bounds on the arcs.
267 268
    /// If it is not used before calling \ref run(), the lower bounds
268 269
    /// will be set to zero on all arcs.
269 270
    ///
270 271
    /// \param map An arc map storing the lower bounds.
271 272
    /// Its \c Value type must be convertible to the \c Value type
272 273
    /// of the algorithm.
273 274
    ///
274 275
    /// \return <tt>(*this)</tt>
275 276
    template <typename LowerMap>
276 277
    CycleCanceling& lowerMap(const LowerMap& map) {
277 278
      _have_lower = true;
278 279
      for (ArcIt a(_graph); a != INVALID; ++a) {
279 280
        _lower[_arc_idf[a]] = map[a];
280 281
        _lower[_arc_idb[a]] = map[a];
281 282
      }
282 283
      return *this;
283 284
    }
284 285

	
285 286
    /// \brief Set the upper bounds (capacities) on the arcs.
286 287
    ///
287 288
    /// This function sets the upper bounds (capacities) on the arcs.
288 289
    /// If it is not used before calling \ref run(), the upper bounds
289 290
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
290 291
    /// unbounded from above).
291 292
    ///
292 293
    /// \param map An arc map storing the upper bounds.
293 294
    /// Its \c Value type must be convertible to the \c Value type
294 295
    /// of the algorithm.
295 296
    ///
296 297
    /// \return <tt>(*this)</tt>
297 298
    template<typename UpperMap>
298 299
    CycleCanceling& upperMap(const UpperMap& map) {
299 300
      for (ArcIt a(_graph); a != INVALID; ++a) {
300 301
        _upper[_arc_idf[a]] = map[a];
301 302
      }
302 303
      return *this;
303 304
    }
304 305

	
305 306
    /// \brief Set the costs of the arcs.
306 307
    ///
307 308
    /// This function sets the costs of the arcs.
308 309
    /// If it is not used before calling \ref run(), the costs
309 310
    /// will be set to \c 1 on all arcs.
310 311
    ///
311 312
    /// \param map An arc map storing the costs.
312 313
    /// Its \c Value type must be convertible to the \c Cost type
313 314
    /// of the algorithm.
314 315
    ///
315 316
    /// \return <tt>(*this)</tt>
316 317
    template<typename CostMap>
317 318
    CycleCanceling& costMap(const CostMap& map) {
318 319
      for (ArcIt a(_graph); a != INVALID; ++a) {
319 320
        _cost[_arc_idf[a]] =  map[a];
320 321
        _cost[_arc_idb[a]] = -map[a];
321 322
      }
322 323
      return *this;
323 324
    }
324 325

	
325 326
    /// \brief Set the supply values of the nodes.
326 327
    ///
327 328
    /// This function sets the supply values of the nodes.
328 329
    /// If neither this function nor \ref stSupply() is used before
329 330
    /// calling \ref run(), the supply of each node will be set to zero.
330 331
    ///
331 332
    /// \param map A node map storing the supply values.
332 333
    /// Its \c Value type must be convertible to the \c Value type
333 334
    /// of the algorithm.
334 335
    ///
335 336
    /// \return <tt>(*this)</tt>
336 337
    template<typename SupplyMap>
337 338
    CycleCanceling& supplyMap(const SupplyMap& map) {
338 339
      for (NodeIt n(_graph); n != INVALID; ++n) {
339 340
        _supply[_node_id[n]] = map[n];
340 341
      }
341 342
      return *this;
342 343
    }
343 344

	
344 345
    /// \brief Set single source and target nodes and a supply value.
345 346
    ///
346 347
    /// This function sets a single source node and a single target node
347 348
    /// and the required flow value.
348 349
    /// If neither this function nor \ref supplyMap() is used before
349 350
    /// calling \ref run(), the supply of each node will be set to zero.
350 351
    ///
351 352
    /// Using this function has the same effect as using \ref supplyMap()
352 353
    /// with such a map in which \c k is assigned to \c s, \c -k is
353 354
    /// assigned to \c t and all other nodes have zero supply value.
354 355
    ///
355 356
    /// \param s The source node.
356 357
    /// \param t The target node.
357 358
    /// \param k The required amount of flow from node \c s to node \c t
358 359
    /// (i.e. the supply of \c s and the demand of \c t).
359 360
    ///
360 361
    /// \return <tt>(*this)</tt>
361 362
    CycleCanceling& stSupply(const Node& s, const Node& t, Value k) {
362 363
      for (int i = 0; i != _res_node_num; ++i) {
363 364
        _supply[i] = 0;
364 365
      }
365 366
      _supply[_node_id[s]] =  k;
366 367
      _supply[_node_id[t]] = -k;
367 368
      return *this;
368 369
    }
369 370

	
370 371
    /// @}
371 372

	
372 373
    /// \name Execution control
373 374
    /// The algorithm can be executed using \ref run().
374 375

	
375 376
    /// @{
376 377

	
377 378
    /// \brief Run the algorithm.
378 379
    ///
379 380
    /// This function runs the algorithm.
380 381
    /// The paramters can be specified using functions \ref lowerMap(),
381 382
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
382 383
    /// For example,
383 384
    /// \code
384 385
    ///   CycleCanceling<ListDigraph> cc(graph);
385 386
    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
386 387
    ///     .supplyMap(sup).run();
387 388
    /// \endcode
388 389
    ///
389 390
    /// This function can be called more than once. All the given parameters
390 391
    /// are kept for the next call, unless \ref resetParams() or \ref reset()
391 392
    /// is used, thus only the modified parameters have to be set again.
392 393
    /// If the underlying digraph was also modified after the construction
393 394
    /// of the class (or the last \ref reset() call), then the \ref reset()
394 395
    /// function must be called.
395 396
    ///
396 397
    /// \param method The cycle-canceling method that will be used.
397 398
    /// For more information, see \ref Method.
398 399
    ///
399 400
    /// \return \c INFEASIBLE if no feasible flow exists,
400 401
    /// \n \c OPTIMAL if the problem has optimal solution
401 402
    /// (i.e. it is feasible and bounded), and the algorithm has found
402 403
    /// optimal flow and node potentials (primal and dual solutions),
403 404
    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
404 405
    /// and infinite upper bound. It means that the objective function
405 406
    /// is unbounded on that arc, however, note that it could actually be
406 407
    /// bounded over the feasible flows, but this algroithm cannot handle
407 408
    /// these cases.
408 409
    ///
409 410
    /// \see ProblemType, Method
410 411
    /// \see resetParams(), reset()
411 412
    ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
412 413
      ProblemType pt = init();
413 414
      if (pt != OPTIMAL) return pt;
414 415
      start(method);
415 416
      return OPTIMAL;
416 417
    }
417 418

	
418 419
    /// \brief Reset all the parameters that have been given before.
419 420
    ///
420 421
    /// This function resets all the paramaters that have been given
421 422
    /// before using functions \ref lowerMap(), \ref upperMap(),
422 423
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
423 424
    ///
424 425
    /// It is useful for multiple \ref run() calls. Basically, all the given
425 426
    /// parameters are kept for the next \ref run() call, unless
426 427
    /// \ref resetParams() or \ref reset() is used.
427 428
    /// If the underlying digraph was also modified after the construction
428 429
    /// of the class or the last \ref reset() call, then the \ref reset()
429 430
    /// function must be used, otherwise \ref resetParams() is sufficient.
430 431
    ///
431 432
    /// For example,
432 433
    /// \code
433 434
    ///   CycleCanceling<ListDigraph> cs(graph);
434 435
    ///
435 436
    ///   // First run
436 437
    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
437 438
    ///     .supplyMap(sup).run();
438 439
    ///
439 440
    ///   // Run again with modified cost map (resetParams() is not called,
440 441
    ///   // so only the cost map have to be set again)
441 442
    ///   cost[e] += 100;
442 443
    ///   cc.costMap(cost).run();
443 444
    ///
444 445
    ///   // Run again from scratch using resetParams()
445 446
    ///   // (the lower bounds will be set to zero on all arcs)
446 447
    ///   cc.resetParams();
447 448
    ///   cc.upperMap(capacity).costMap(cost)
448 449
    ///     .supplyMap(sup).run();
449 450
    /// \endcode
450 451
    ///
451 452
    /// \return <tt>(*this)</tt>
452 453
    ///
Ignore white space 6 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_KRUSKAL_H
20 20
#define LEMON_KRUSKAL_H
21 21

	
22 22
#include <algorithm>
23 23
#include <vector>
24 24
#include <lemon/unionfind.h>
25 25
#include <lemon/maps.h>
26 26

	
27 27
#include <lemon/core.h>
28 28
#include <lemon/bits/traits.h>
29 29

	
30 30
///\ingroup spantree
31 31
///\file
32 32
///\brief Kruskal's algorithm to compute a minimum cost spanning tree
33
///
34
///Kruskal's algorithm to compute a minimum cost spanning tree.
35
///
36 33

	
37 34
namespace lemon {
38 35

	
39 36
  namespace _kruskal_bits {
40 37

	
41 38
    // Kruskal for directed graphs.
42 39

	
43 40
    template <typename Digraph, typename In, typename Out>
44 41
    typename disable_if<lemon::UndirectedTagIndicator<Digraph>,
45 42
                       typename In::value_type::second_type >::type
46 43
    kruskal(const Digraph& digraph, const In& in, Out& out,dummy<0> = 0) {
47 44
      typedef typename In::value_type::second_type Value;
48 45
      typedef typename Digraph::template NodeMap<int> IndexMap;
49 46
      typedef typename Digraph::Node Node;
50 47

	
51 48
      IndexMap index(digraph);
52 49
      UnionFind<IndexMap> uf(index);
53 50
      for (typename Digraph::NodeIt it(digraph); it != INVALID; ++it) {
54 51
        uf.insert(it);
55 52
      }
56 53

	
57 54
      Value tree_value = 0;
58 55
      for (typename In::const_iterator it = in.begin(); it != in.end(); ++it) {
59 56
        if (uf.join(digraph.target(it->first),digraph.source(it->first))) {
60 57
          out.set(it->first, true);
61 58
          tree_value += it->second;
62 59
        }
63 60
        else {
64 61
          out.set(it->first, false);
65 62
        }
66 63
      }
67 64
      return tree_value;
68 65
    }
69 66

	
70 67
    // Kruskal for undirected graphs.
71 68

	
72 69
    template <typename Graph, typename In, typename Out>
73 70
    typename enable_if<lemon::UndirectedTagIndicator<Graph>,
74 71
                       typename In::value_type::second_type >::type
75 72
    kruskal(const Graph& graph, const In& in, Out& out,dummy<1> = 1) {
76 73
      typedef typename In::value_type::second_type Value;
77 74
      typedef typename Graph::template NodeMap<int> IndexMap;
78 75
      typedef typename Graph::Node Node;
79 76

	
80 77
      IndexMap index(graph);
81 78
      UnionFind<IndexMap> uf(index);
82 79
      for (typename Graph::NodeIt it(graph); it != INVALID; ++it) {
83 80
        uf.insert(it);
84 81
      }
85 82

	
86 83
      Value tree_value = 0;
87 84
      for (typename In::const_iterator it = in.begin(); it != in.end(); ++it) {
88 85
        if (uf.join(graph.u(it->first),graph.v(it->first))) {
89 86
          out.set(it->first, true);
90 87
          tree_value += it->second;
91 88
        }
92 89
        else {
93 90
          out.set(it->first, false);
94 91
        }
95 92
      }
96 93
      return tree_value;
97 94
    }
98 95

	
99 96

	
100 97
    template <typename Sequence>
101 98
    struct PairComp {
102 99
      typedef typename Sequence::value_type Value;
103 100
      bool operator()(const Value& left, const Value& right) {
104 101
        return left.second < right.second;
105 102
      }
106 103
    };
107 104

	
108 105
    template <typename In, typename Enable = void>
109 106
    struct SequenceInputIndicator {
110 107
      static const bool value = false;
111 108
    };
112 109

	
113 110
    template <typename In>
114 111
    struct SequenceInputIndicator<In,
115 112
      typename exists<typename In::value_type::first_type>::type> {
116 113
      static const bool value = true;
117 114
    };
118 115

	
119 116
    template <typename In, typename Enable = void>
120 117
    struct MapInputIndicator {
121 118
      static const bool value = false;
122 119
    };
123 120

	
124 121
    template <typename In>
125 122
    struct MapInputIndicator<In,
126 123
      typename exists<typename In::Value>::type> {
127 124
      static const bool value = true;
128 125
    };
129 126

	
130 127
    template <typename In, typename Enable = void>
131 128
    struct SequenceOutputIndicator {
132 129
      static const bool value = false;
133 130
    };
134 131

	
135 132
    template <typename Out>
136 133
    struct SequenceOutputIndicator<Out,
137 134
      typename exists<typename Out::value_type>::type> {
138 135
      static const bool value = true;
139 136
    };
140 137

	
141 138
    template <typename Out, typename Enable = void>
142 139
    struct MapOutputIndicator {
143 140
      static const bool value = false;
144 141
    };
145 142

	
146 143
    template <typename Out>
147 144
    struct MapOutputIndicator<Out,
148 145
      typename exists<typename Out::Value>::type> {
149 146
      static const bool value = true;
150 147
    };
151 148

	
152 149
    template <typename In, typename InEnable = void>
153 150
    struct KruskalValueSelector {};
154 151

	
155 152
    template <typename In>
156 153
    struct KruskalValueSelector<In,
157 154
      typename enable_if<SequenceInputIndicator<In>, void>::type>
158 155
    {
159 156
      typedef typename In::value_type::second_type Value;
160 157
    };
161 158

	
162 159
    template <typename In>
163 160
    struct KruskalValueSelector<In,
164 161
      typename enable_if<MapInputIndicator<In>, void>::type>
165 162
    {
166 163
      typedef typename In::Value Value;
167 164
    };
168 165

	
169 166
    template <typename Graph, typename In, typename Out,
170 167
              typename InEnable = void>
171 168
    struct KruskalInputSelector {};
172 169

	
173 170
    template <typename Graph, typename In, typename Out,
174 171
              typename InEnable = void>
175 172
    struct KruskalOutputSelector {};
176 173

	
177 174
    template <typename Graph, typename In, typename Out>
178 175
    struct KruskalInputSelector<Graph, In, Out,
179 176
      typename enable_if<SequenceInputIndicator<In>, void>::type >
180 177
    {
181 178
      typedef typename In::value_type::second_type Value;
182 179

	
183 180
      static Value kruskal(const Graph& graph, const In& in, Out& out) {
184 181
        return KruskalOutputSelector<Graph, In, Out>::
185 182
          kruskal(graph, in, out);
186 183
      }
187 184

	
188 185
    };
189 186

	
190 187
    template <typename Graph, typename In, typename Out>
191 188
    struct KruskalInputSelector<Graph, In, Out,
192 189
      typename enable_if<MapInputIndicator<In>, void>::type >
193 190
    {
194 191
      typedef typename In::Value Value;
195 192
      static Value kruskal(const Graph& graph, const In& in, Out& out) {
196 193
        typedef typename In::Key MapArc;
197 194
        typedef typename In::Value Value;
198 195
        typedef typename ItemSetTraits<Graph, MapArc>::ItemIt MapArcIt;
199 196
        typedef std::vector<std::pair<MapArc, Value> > Sequence;
200 197
        Sequence seq;
201 198

	
202 199
        for (MapArcIt it(graph); it != INVALID; ++it) {
203 200
          seq.push_back(std::make_pair(it, in[it]));
204 201
        }
205 202

	
206 203
        std::sort(seq.begin(), seq.end(), PairComp<Sequence>());
207 204
        return KruskalOutputSelector<Graph, Sequence, Out>::
208 205
          kruskal(graph, seq, out);
209 206
      }
210 207
    };
211 208

	
212 209
    template <typename T>
213 210
    struct RemoveConst {
214 211
      typedef T type;
215 212
    };
216 213

	
217 214
    template <typename T>
218 215
    struct RemoveConst<const T> {
219 216
      typedef T type;
220 217
    };
221 218

	
222 219
    template <typename Graph, typename In, typename Out>
223 220
    struct KruskalOutputSelector<Graph, In, Out,
224 221
      typename enable_if<SequenceOutputIndicator<Out>, void>::type >
225 222
    {
226 223
      typedef typename In::value_type::second_type Value;
227 224

	
228 225
      static Value kruskal(const Graph& graph, const In& in, Out& out) {
229 226
        typedef LoggerBoolMap<typename RemoveConst<Out>::type> Map;
230 227
        Map map(out);
231 228
        return _kruskal_bits::kruskal(graph, in, map);
232 229
      }
233 230

	
234 231
    };
235 232

	
236 233
    template <typename Graph, typename In, typename Out>
237 234
    struct KruskalOutputSelector<Graph, In, Out,
238 235
      typename enable_if<MapOutputIndicator<Out>, void>::type >
239 236
    {
240 237
      typedef typename In::value_type::second_type Value;
241 238

	
242 239
      static Value kruskal(const Graph& graph, const In& in, Out& out) {
243 240
        return _kruskal_bits::kruskal(graph, in, out);
244 241
      }
245 242
    };
246 243

	
247 244
  }
248 245

	
249 246
  /// \ingroup spantree
250 247
  ///
251 248
  /// \brief Kruskal's algorithm for finding a minimum cost spanning tree of
252 249
  /// a graph.
253 250
  ///
254 251
  /// This function runs Kruskal's algorithm to find a minimum cost
255 252
  /// spanning tree of a graph.
256 253
  /// Due to some C++ hacking, it accepts various input and output types.
257 254
  ///
258 255
  /// \param g The graph the algorithm runs on.
259 256
  /// It can be either \ref concepts::Digraph "directed" or
260 257
  /// \ref concepts::Graph "undirected".
261 258
  /// If the graph is directed, the algorithm consider it to be
262 259
  /// undirected by disregarding the direction of the arcs.
263 260
  ///
264 261
  /// \param in This object is used to describe the arc/edge costs.
265 262
  /// It can be one of the following choices.
266 263
  /// - An STL compatible 'Forward Container' with
267 264
  /// <tt>std::pair<GR::Arc,C></tt> or
268 265
  /// <tt>std::pair<GR::Edge,C></tt> as its <tt>value_type</tt>, where
269 266
  /// \c C is the type of the costs. The pairs indicates the arcs/edges
270 267
  /// along with the assigned cost. <em>They must be in a
271 268
  /// cost-ascending order.</em>
272 269
  /// - Any readable arc/edge map. The values of the map indicate the
273 270
  /// arc/edge costs.
274 271
  ///
275 272
  /// \retval out Here we also have a choice.
276 273
  /// - It can be a writable arc/edge map with \c bool value type. After
277 274
  /// running the algorithm it will contain the found minimum cost spanning
278 275
  /// tree: the value of an arc/edge will be set to \c true if it belongs
279 276
  /// to the tree, otherwise it will be set to \c false. The value of
280 277
  /// each arc/edge will be set exactly once.
281 278
  /// - It can also be an iteraror of an STL Container with
282 279
  /// <tt>GR::Arc</tt> or <tt>GR::Edge</tt> as its
283 280
  /// <tt>value_type</tt>.  The algorithm copies the elements of the
284 281
  /// found tree into this sequence.  For example, if we know that the
285 282
  /// spanning tree of the graph \c g has say 53 arcs, then we can
286 283
  /// put its arcs into an STL vector \c tree with a code like this.
287 284
  ///\code
288 285
  /// std::vector<Arc> tree(53);
289 286
  /// kruskal(g,cost,tree.begin());
290 287
  ///\endcode
291 288
  /// Or if we don't know in advance the size of the tree, we can
292 289
  /// write this.
293 290
  ///\code
294 291
  /// std::vector<Arc> tree;
295 292
  /// kruskal(g,cost,std::back_inserter(tree));
296 293
  ///\endcode
297 294
  ///
298 295
  /// \return The total cost of the found spanning tree.
299 296
  ///
300 297
  /// \note If the input graph is not (weakly) connected, a spanning
301 298
  /// forest is calculated instead of a spanning tree.
302 299

	
303 300
#ifdef DOXYGEN
304 301
  template <typename Graph, typename In, typename Out>
305 302
  Value kruskal(const Graph& g, const In& in, Out& out)
306 303
#else
307 304
  template <class Graph, class In, class Out>
308 305
  inline typename _kruskal_bits::KruskalValueSelector<In>::Value
309 306
  kruskal(const Graph& graph, const In& in, Out& out)
310 307
#endif
311 308
  {
312 309
    return _kruskal_bits::KruskalInputSelector<Graph, In, Out>::
313 310
      kruskal(graph, in, out);
314 311
  }
315 312

	
316 313

	
317 314
  template <class Graph, class In, class Out>
318 315
  inline typename _kruskal_bits::KruskalValueSelector<In>::Value
319 316
  kruskal(const Graph& graph, const In& in, const Out& out)
320 317
  {
321 318
    return _kruskal_bits::KruskalInputSelector<Graph, In, const Out>::
322 319
      kruskal(graph, in, out);
323 320
  }
324 321

	
325 322
} //namespace lemon
326 323

	
327 324
#endif //LEMON_KRUSKAL_H
Ignore white space 768 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2010
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_NETWORK_SIMPLEX_H
20 20
#define LEMON_NETWORK_SIMPLEX_H
21 21

	
22 22
/// \ingroup min_cost_flow_algs
23 23
///
24 24
/// \file
25 25
/// \brief Network Simplex algorithm for finding a minimum cost flow.
26 26

	
27 27
#include <vector>
28 28
#include <limits>
29 29
#include <algorithm>
30 30

	
31 31
#include <lemon/core.h>
32 32
#include <lemon/math.h>
33 33

	
34 34
namespace lemon {
35 35

	
36 36
  /// \addtogroup min_cost_flow_algs
37 37
  /// @{
38 38

	
39 39
  /// \brief Implementation of the primal Network Simplex algorithm
40 40
  /// for finding a \ref min_cost_flow "minimum cost flow".
41 41
  ///
42 42
  /// \ref NetworkSimplex implements the primal Network Simplex algorithm
43 43
  /// for finding a \ref min_cost_flow "minimum cost flow"
44 44
  /// \ref amo93networkflows, \ref dantzig63linearprog,
45 45
  /// \ref kellyoneill91netsimplex.
46 46
  /// This algorithm is a highly efficient specialized version of the
47 47
  /// linear programming simplex method directly for the minimum cost
48 48
  /// flow problem.
49 49
  ///
50 50
  /// In general, %NetworkSimplex is the fastest implementation available
51 51
  /// in LEMON for this problem.
52 52
  /// Moreover, it supports both directions of the supply/demand inequality
53 53
  /// constraints. For more information, see \ref SupplyType.
54 54
  ///
55 55
  /// Most of the parameters of the problem (except for the digraph)
56 56
  /// can be given using separate functions, and the algorithm can be
57 57
  /// executed using the \ref run() function. If some parameters are not
58 58
  /// specified, then default values will be used.
59 59
  ///
60 60
  /// \tparam GR The digraph type the algorithm runs on.
61 61
  /// \tparam V The number type used for flow amounts, capacity bounds
62 62
  /// and supply values in the algorithm. By default, it is \c int.
63 63
  /// \tparam C The number type used for costs and potentials in the
64 64
  /// algorithm. By default, it is the same as \c V.
65 65
  ///
66
  /// \warning Both number types must be signed and all input data must
66
  /// \warning Both \c V and \c C must be signed number types.
67
  /// \warning All input data (capacities, supply values, and costs) must
67 68
  /// be integer.
68 69
  ///
69 70
  /// \note %NetworkSimplex provides five different pivot rule
70 71
  /// implementations, from which the most efficient one is used
71 72
  /// by default. For more information, see \ref PivotRule.
72 73
  template <typename GR, typename V = int, typename C = V>
73 74
  class NetworkSimplex
74 75
  {
75 76
  public:
76 77

	
77 78
    /// The type of the flow amounts, capacity bounds and supply values
78 79
    typedef V Value;
79 80
    /// The type of the arc costs
80 81
    typedef C Cost;
81 82

	
82 83
  public:
83 84

	
84 85
    /// \brief Problem type constants for the \c run() function.
85 86
    ///
86 87
    /// Enum type containing the problem type constants that can be
87 88
    /// returned by the \ref run() function of the algorithm.
88 89
    enum ProblemType {
89 90
      /// The problem has no feasible solution (flow).
90 91
      INFEASIBLE,
91 92
      /// The problem has optimal solution (i.e. it is feasible and
92 93
      /// bounded), and the algorithm has found optimal flow and node
93 94
      /// potentials (primal and dual solutions).
94 95
      OPTIMAL,
95 96
      /// The objective function of the problem is unbounded, i.e.
96 97
      /// there is a directed cycle having negative total cost and
97 98
      /// infinite upper bound.
98 99
      UNBOUNDED
99 100
    };
100 101

	
101 102
    /// \brief Constants for selecting the type of the supply constraints.
102 103
    ///
103 104
    /// Enum type containing constants for selecting the supply type,
104 105
    /// i.e. the direction of the inequalities in the supply/demand
105 106
    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
106 107
    ///
107 108
    /// The default supply type is \c GEQ, the \c LEQ type can be
108 109
    /// selected using \ref supplyType().
109 110
    /// The equality form is a special case of both supply types.
110 111
    enum SupplyType {
111 112
      /// This option means that there are <em>"greater or equal"</em>
112 113
      /// supply/demand constraints in the definition of the problem.
113 114
      GEQ,
114 115
      /// This option means that there are <em>"less or equal"</em>
115 116
      /// supply/demand constraints in the definition of the problem.
116 117
      LEQ
117 118
    };
118 119

	
119 120
    /// \brief Constants for selecting the pivot rule.
120 121
    ///
121 122
    /// Enum type containing constants for selecting the pivot rule for
122 123
    /// the \ref run() function.
123 124
    ///
124 125
    /// \ref NetworkSimplex provides five different pivot rule
125 126
    /// implementations that significantly affect the running time
126 127
    /// of the algorithm.
127 128
    /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
128 129
    /// proved to be the most efficient and the most robust on various
129 130
    /// test inputs.
130 131
    /// However, another pivot rule can be selected using the \ref run()
131 132
    /// function with the proper parameter.
132 133
    enum PivotRule {
133 134

	
134 135
      /// The \e First \e Eligible pivot rule.
135 136
      /// The next eligible arc is selected in a wraparound fashion
136 137
      /// in every iteration.
137 138
      FIRST_ELIGIBLE,
138 139

	
139 140
      /// The \e Best \e Eligible pivot rule.
140 141
      /// The best eligible arc is selected in every iteration.
141 142
      BEST_ELIGIBLE,
142 143

	
143 144
      /// The \e Block \e Search pivot rule.
144 145
      /// A specified number of arcs are examined in every iteration
145 146
      /// in a wraparound fashion and the best eligible arc is selected
146 147
      /// from this block.
147 148
      BLOCK_SEARCH,
148 149

	
149 150
      /// The \e Candidate \e List pivot rule.
150 151
      /// In a major iteration a candidate list is built from eligible arcs
151 152
      /// in a wraparound fashion and in the following minor iterations
152 153
      /// the best eligible arc is selected from this list.
153 154
      CANDIDATE_LIST,
154 155

	
155 156
      /// The \e Altering \e Candidate \e List pivot rule.
156 157
      /// It is a modified version of the Candidate List method.
157 158
      /// It keeps only the several best eligible arcs from the former
158 159
      /// candidate list and extends this list in every iteration.
159 160
      ALTERING_LIST
160 161
    };
161 162

	
162 163
  private:
163 164

	
164 165
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
165 166

	
166 167
    typedef std::vector<int> IntVector;
167 168
    typedef std::vector<Value> ValueVector;
168 169
    typedef std::vector<Cost> CostVector;
169 170
    typedef std::vector<signed char> CharVector;
170 171
    // Note: vector<signed char> is used instead of vector<ArcState> and 
171 172
    // vector<ArcDirection> for efficiency reasons
172 173

	
173 174
    // State constants for arcs
174 175
    enum ArcState {
175 176
      STATE_UPPER = -1,
176 177
      STATE_TREE  =  0,
177 178
      STATE_LOWER =  1
178 179
    };
179 180

	
180 181
    // Direction constants for tree arcs
181 182
    enum ArcDirection {
182 183
      DIR_DOWN = -1,
183 184
      DIR_UP   =  1
184 185
    };
185 186

	
186 187
  private:
187 188

	
188 189
    // Data related to the underlying digraph
189 190
    const GR &_graph;
190 191
    int _node_num;
191 192
    int _arc_num;
192 193
    int _all_arc_num;
193 194
    int _search_arc_num;
194 195

	
195 196
    // Parameters of the problem
196 197
    bool _have_lower;
197 198
    SupplyType _stype;
198 199
    Value _sum_supply;
199 200

	
200 201
    // Data structures for storing the digraph
201 202
    IntNodeMap _node_id;
202 203
    IntArcMap _arc_id;
203 204
    IntVector _source;
204 205
    IntVector _target;
205 206
    bool _arc_mixing;
206 207

	
207 208
    // Node and arc data
208 209
    ValueVector _lower;
209 210
    ValueVector _upper;
210 211
    ValueVector _cap;
211 212
    CostVector _cost;
212 213
    ValueVector _supply;
213 214
    ValueVector _flow;
214 215
    CostVector _pi;
215 216

	
216 217
    // Data for storing the spanning tree structure
217 218
    IntVector _parent;
218 219
    IntVector _pred;
219 220
    IntVector _thread;
220 221
    IntVector _rev_thread;
221 222
    IntVector _succ_num;
222 223
    IntVector _last_succ;
223 224
    CharVector _pred_dir;
224 225
    CharVector _state;
225 226
    IntVector _dirty_revs;
226 227
    int _root;
227 228

	
228 229
    // Temporary data used in the current pivot iteration
229 230
    int in_arc, join, u_in, v_in, u_out, v_out;
230 231
    Value delta;
231 232

	
232 233
    const Value MAX;
233 234

	
234 235
  public:
235 236

	
236 237
    /// \brief Constant for infinite upper bounds (capacities).
237 238
    ///
238 239
    /// Constant for infinite upper bounds (capacities).
239 240
    /// It is \c std::numeric_limits<Value>::infinity() if available,
240 241
    /// \c std::numeric_limits<Value>::max() otherwise.
241 242
    const Value INF;
242 243

	
243 244
  private:
244 245

	
245 246
    // Implementation of the First Eligible pivot rule
246 247
    class FirstEligiblePivotRule
247 248
    {
248 249
    private:
249 250

	
250 251
      // References to the NetworkSimplex class
251 252
      const IntVector  &_source;
252 253
      const IntVector  &_target;
253 254
      const CostVector &_cost;
254 255
      const CharVector &_state;
255 256
      const CostVector &_pi;
256 257
      int &_in_arc;
257 258
      int _search_arc_num;
258 259

	
259 260
      // Pivot rule data
260 261
      int _next_arc;
261 262

	
262 263
    public:
263 264

	
264 265
      // Constructor
265 266
      FirstEligiblePivotRule(NetworkSimplex &ns) :
266 267
        _source(ns._source), _target(ns._target),
267 268
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
268 269
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
269 270
        _next_arc(0)
270 271
      {}
271 272

	
272 273
      // Find next entering arc
273 274
      bool findEnteringArc() {
274 275
        Cost c;
275 276
        for (int e = _next_arc; e != _search_arc_num; ++e) {
276 277
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
277 278
          if (c < 0) {
278 279
            _in_arc = e;
279 280
            _next_arc = e + 1;
280 281
            return true;
281 282
          }
282 283
        }
283 284
        for (int e = 0; e != _next_arc; ++e) {
284 285
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
285 286
          if (c < 0) {
286 287
            _in_arc = e;
287 288
            _next_arc = e + 1;
288 289
            return true;
289 290
          }
290 291
        }
291 292
        return false;
292 293
      }
293 294

	
294 295
    }; //class FirstEligiblePivotRule
295 296

	
296 297

	
297 298
    // Implementation of the Best Eligible pivot rule
298 299
    class BestEligiblePivotRule
299 300
    {
300 301
    private:
301 302

	
302 303
      // References to the NetworkSimplex class
303 304
      const IntVector  &_source;
304 305
      const IntVector  &_target;
305 306
      const CostVector &_cost;
306 307
      const CharVector &_state;
307 308
      const CostVector &_pi;
308 309
      int &_in_arc;
309 310
      int _search_arc_num;
310 311

	
311 312
    public:
312 313

	
313 314
      // Constructor
314 315
      BestEligiblePivotRule(NetworkSimplex &ns) :
315 316
        _source(ns._source), _target(ns._target),
316 317
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
317 318
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
318 319
      {}
319 320

	
320 321
      // Find next entering arc
321 322
      bool findEnteringArc() {
322 323
        Cost c, min = 0;
323 324
        for (int e = 0; e != _search_arc_num; ++e) {
324 325
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
325 326
          if (c < min) {
326 327
            min = c;
327 328
            _in_arc = e;
328 329
          }
329 330
        }
330 331
        return min < 0;
331 332
      }
332 333

	
333 334
    }; //class BestEligiblePivotRule
334 335

	
335 336

	
336 337
    // Implementation of the Block Search pivot rule
337 338
    class BlockSearchPivotRule
338 339
    {
339 340
    private:
340 341

	
341 342
      // References to the NetworkSimplex class
342 343
      const IntVector  &_source;
343 344
      const IntVector  &_target;
344 345
      const CostVector &_cost;
345 346
      const CharVector &_state;
346 347
      const CostVector &_pi;
347 348
      int &_in_arc;
348 349
      int _search_arc_num;
349 350

	
350 351
      // Pivot rule data
351 352
      int _block_size;
352 353
      int _next_arc;
353 354

	
354 355
    public:
355 356

	
356 357
      // Constructor
357 358
      BlockSearchPivotRule(NetworkSimplex &ns) :
358 359
        _source(ns._source), _target(ns._target),
359 360
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
360 361
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
361 362
        _next_arc(0)
362 363
      {
363 364
        // The main parameters of the pivot rule
364 365
        const double BLOCK_SIZE_FACTOR = 1.0;
365 366
        const int MIN_BLOCK_SIZE = 10;
366 367

	
367 368
        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
368 369
                                    std::sqrt(double(_search_arc_num))),
369 370
                                MIN_BLOCK_SIZE );
370 371
      }
371 372

	
372 373
      // Find next entering arc
373 374
      bool findEnteringArc() {
374 375
        Cost c, min = 0;
375 376
        int cnt = _block_size;
376 377
        int e;
377 378
        for (e = _next_arc; e != _search_arc_num; ++e) {
378 379
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
379 380
          if (c < min) {
380 381
            min = c;
381 382
            _in_arc = e;
382 383
          }
383 384
          if (--cnt == 0) {
384 385
            if (min < 0) goto search_end;
385 386
            cnt = _block_size;
386 387
          }
387 388
        }
388 389
        for (e = 0; e != _next_arc; ++e) {
389 390
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
390 391
          if (c < min) {
391 392
            min = c;
392 393
            _in_arc = e;
393 394
          }
394 395
          if (--cnt == 0) {
395 396
            if (min < 0) goto search_end;
396 397
            cnt = _block_size;
397 398
          }
398 399
        }
399 400
        if (min >= 0) return false;
400 401

	
401 402
      search_end:
402 403
        _next_arc = e;
403 404
        return true;
404 405
      }
405 406

	
406 407
    }; //class BlockSearchPivotRule
407 408

	
408 409

	
409 410
    // Implementation of the Candidate List pivot rule
410 411
    class CandidateListPivotRule
411 412
    {
412 413
    private:
413 414

	
414 415
      // References to the NetworkSimplex class
415 416
      const IntVector  &_source;
416 417
      const IntVector  &_target;
417 418
      const CostVector &_cost;
418 419
      const CharVector &_state;
419 420
      const CostVector &_pi;
420 421
      int &_in_arc;
421 422
      int _search_arc_num;
422 423

	
423 424
      // Pivot rule data
424 425
      IntVector _candidates;
425 426
      int _list_length, _minor_limit;
426 427
      int _curr_length, _minor_count;
427 428
      int _next_arc;
428 429

	
429 430
    public:
430 431

	
431 432
      /// Constructor
432 433
      CandidateListPivotRule(NetworkSimplex &ns) :
433 434
        _source(ns._source), _target(ns._target),
434 435
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
435 436
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
436 437
        _next_arc(0)
437 438
      {
438 439
        // The main parameters of the pivot rule
439 440
        const double LIST_LENGTH_FACTOR = 0.25;
440 441
        const int MIN_LIST_LENGTH = 10;
441 442
        const double MINOR_LIMIT_FACTOR = 0.1;
442 443
        const int MIN_MINOR_LIMIT = 3;
443 444

	
444 445
        _list_length = std::max( int(LIST_LENGTH_FACTOR *
445 446
                                     std::sqrt(double(_search_arc_num))),
446 447
                                 MIN_LIST_LENGTH );
447 448
        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
448 449
                                 MIN_MINOR_LIMIT );
449 450
        _curr_length = _minor_count = 0;
450 451
        _candidates.resize(_list_length);
0 comments (0 inline)