gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Small doc improvements + unifications in MCF classes (#180)
0 3 0
default
3 files changed with 32 insertions and 32 deletions:
↑ Collapse diff ↑
Ignore white space 384 line context
1 1
/* -*- C++ -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library
4 4
 *
5 5
 * Copyright (C) 2003-2008
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
  /// \tparam V The value type used for flow amounts, capacity bounds
38
  /// \tparam V The number type used for flow amounts, capacity bounds
39 39
  /// and supply values. By default it is \c int.
40
  /// \tparam C The value type used for costs and potentials.
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". It is an efficient dual
70 70
  /// solution method.
71 71
  ///
72 72
  /// Most of the parameters of the problem (except for the digraph)
73 73
  /// can be given using separate functions, and the algorithm can be
74 74
  /// executed using the \ref run() function. If some parameters are not
75 75
  /// specified, then default values will be used.
76 76
  ///
77 77
  /// \tparam GR The digraph type the algorithm runs on.
78
  /// \tparam V The value type used for flow amounts, capacity bounds
78
  /// \tparam V The number type used for flow amounts, capacity bounds
79 79
  /// and supply values in the algorithm. By default it is \c int.
80
  /// \tparam C The value type used for costs and potentials in the
80
  /// \tparam C The number type used for costs and potentials in the
81 81
  /// algorithm. By default it is the same as \c V.
82 82
  ///
83
  /// \warning Both value types must be signed and all input data must
83
  /// \warning Both number types must be signed and all input data must
84 84
  /// be integer.
85 85
  /// \warning This algorithm does not support negative costs for such
86 86
  /// arcs that have infinite upper bound.
87 87
#ifdef DOXYGEN
88 88
  template <typename GR, typename V, typename C, typename TR>
89 89
#else
90 90
  template < typename GR, typename V = int, typename C = V,
91 91
             typename TR = CapacityScalingDefaultTraits<GR, V, C> >
92 92
#endif
93 93
  class CapacityScaling
94 94
  {
95 95
  public:
96 96

	
97 97
    /// The type of the digraph
98 98
    typedef typename TR::Digraph Digraph;
99 99
    /// The type of the flow amounts, capacity bounds and supply values
100 100
    typedef typename TR::Value Value;
101 101
    /// The type of the arc costs
102 102
    typedef typename TR::Cost Cost;
103 103

	
104 104
    /// The type of the heap used for internal Dijkstra computations
105 105
    typedef typename TR::Heap Heap;
106 106

	
107 107
    /// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm
108 108
    typedef TR Traits;
109 109

	
110 110
  public:
111 111

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

	
133 133
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
134 134

	
135 135
    typedef std::vector<int> IntVector;
136 136
    typedef std::vector<char> BoolVector;
137 137
    typedef std::vector<Value> ValueVector;
138 138
    typedef std::vector<Cost> CostVector;
139 139

	
140 140
  private:
141 141

	
142 142
    // Data related to the underlying digraph
143 143
    const GR &_graph;
144 144
    int _node_num;
145 145
    int _arc_num;
146 146
    int _res_arc_num;
147 147
    int _root;
148 148

	
149 149
    // Parameters of the problem
150 150
    bool _have_lower;
151 151
    Value _sum_supply;
152 152

	
153 153
    // Data structures for storing the digraph
154 154
    IntNodeMap _node_id;
155 155
    IntArcMap _arc_idf;
156 156
    IntArcMap _arc_idb;
157 157
    IntVector _first_out;
158 158
    BoolVector _forward;
159 159
    IntVector _source;
160 160
    IntVector _target;
161 161
    IntVector _reverse;
162 162

	
163 163
    // Node and arc data
164 164
    ValueVector _lower;
165 165
    ValueVector _upper;
166 166
    CostVector _cost;
167 167
    ValueVector _supply;
168 168

	
169 169
    ValueVector _res_cap;
170 170
    CostVector _pi;
171 171
    ValueVector _excess;
172 172
    IntVector _excess_nodes;
173 173
    IntVector _deficit_nodes;
174 174

	
175 175
    Value _delta;
176 176
    int _factor;
177 177
    IntVector _pred;
178 178

	
179 179
  public:
180 180
  
181 181
    /// \brief Constant for infinite upper bounds (capacities).
182 182
    ///
183 183
    /// Constant for infinite upper bounds (capacities).
184 184
    /// It is \c std::numeric_limits<Value>::infinity() if available,
185 185
    /// \c std::numeric_limits<Value>::max() otherwise.
186 186
    const Value INF;
187 187

	
188 188
  private:
189 189

	
190 190
    // Special implementation of the Dijkstra algorithm for finding
191 191
    // shortest paths in the residual network of the digraph with
192 192
    // respect to the reduced arc costs and modifying the node
193 193
    // potentials according to the found distance labels.
194 194
    class ResidualDijkstra
195 195
    {
196 196
    private:
197 197

	
198 198
      int _node_num;
199 199
      bool _geq;
200 200
      const IntVector &_first_out;
201 201
      const IntVector &_target;
202 202
      const CostVector &_cost;
203 203
      const ValueVector &_res_cap;
204 204
      const ValueVector &_excess;
205 205
      CostVector &_pi;
206 206
      IntVector &_pred;
207 207
      
208 208
      IntVector _proc_nodes;
209 209
      CostVector _dist;
210 210
      
211 211
    public:
212 212

	
213 213
      ResidualDijkstra(CapacityScaling& cs) :
214 214
        _node_num(cs._node_num), _geq(cs._sum_supply < 0),
215 215
        _first_out(cs._first_out), _target(cs._target), _cost(cs._cost),
216 216
        _res_cap(cs._res_cap), _excess(cs._excess), _pi(cs._pi),
217 217
        _pred(cs._pred), _dist(cs._node_num)
218 218
      {}
219 219

	
220 220
      int run(int s, Value delta = 1) {
221 221
        RangeMap<int> heap_cross_ref(_node_num, Heap::PRE_HEAP);
222 222
        Heap heap(heap_cross_ref);
223 223
        heap.push(s, 0);
224 224
        _pred[s] = -1;
225 225
        _proc_nodes.clear();
226 226

	
227 227
        // Process nodes
228 228
        while (!heap.empty() && _excess[heap.top()] > -delta) {
229 229
          int u = heap.top(), v;
230 230
          Cost d = heap.prio() + _pi[u], dn;
231 231
          _dist[u] = heap.prio();
232 232
          _proc_nodes.push_back(u);
233 233
          heap.pop();
234 234

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

	
259 259
        // Update potentials of processed nodes
260 260
        int t = heap.top();
261 261
        Cost dt = heap.prio();
262 262
        for (int i = 0; i < int(_proc_nodes.size()); ++i) {
263 263
          _pi[_proc_nodes[i]] += _dist[_proc_nodes[i]] - dt;
264 264
        }
265 265

	
266 266
        return t;
267 267
      }
268 268

	
269 269
    }; //class ResidualDijkstra
270 270

	
271 271
  public:
272 272

	
273 273
    /// \name Named Template Parameters
274 274
    /// @{
275 275

	
276 276
    template <typename T>
277 277
    struct SetHeapTraits : public Traits {
278 278
      typedef T Heap;
279 279
    };
280 280

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

	
295 295
    /// @}
296 296

	
297 297
  public:
298 298

	
299 299
    /// \brief Constructor.
300 300
    ///
301 301
    /// The constructor of the class.
302 302
    ///
303 303
    /// \param graph The digraph the algorithm runs on.
304 304
    CapacityScaling(const GR& graph) :
305 305
      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
306 306
      INF(std::numeric_limits<Value>::has_infinity ?
307 307
          std::numeric_limits<Value>::infinity() :
308 308
          std::numeric_limits<Value>::max())
309 309
    {
310
      // Check the value types
310
      // Check the number types
311 311
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
312 312
        "The flow type of CapacityScaling must be signed");
313 313
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
314 314
        "The cost type of CapacityScaling must be signed");
315 315

	
316 316
      // Resize vectors
317 317
      _node_num = countNodes(_graph);
318 318
      _arc_num = countArcs(_graph);
319 319
      _res_arc_num = 2 * (_arc_num + _node_num);
320 320
      _root = _node_num;
321 321
      ++_node_num;
322 322

	
323 323
      _first_out.resize(_node_num + 1);
324 324
      _forward.resize(_res_arc_num);
325 325
      _source.resize(_res_arc_num);
326 326
      _target.resize(_res_arc_num);
327 327
      _reverse.resize(_res_arc_num);
328 328

	
329 329
      _lower.resize(_res_arc_num);
330 330
      _upper.resize(_res_arc_num);
331 331
      _cost.resize(_res_arc_num);
332 332
      _supply.resize(_node_num);
333 333
      
334 334
      _res_cap.resize(_res_arc_num);
335 335
      _pi.resize(_node_num);
336 336
      _excess.resize(_node_num);
337 337
      _pred.resize(_node_num);
338 338

	
339 339
      // Copy the graph
340 340
      int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
341 341
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
342 342
        _node_id[n] = i;
343 343
      }
344 344
      i = 0;
345 345
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
346 346
        _first_out[i] = j;
347 347
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
348 348
          _arc_idf[a] = j;
349 349
          _forward[j] = true;
350 350
          _source[j] = i;
351 351
          _target[j] = _node_id[_graph.runningNode(a)];
352 352
        }
353 353
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
354 354
          _arc_idb[a] = j;
355 355
          _forward[j] = false;
356 356
          _source[j] = i;
357 357
          _target[j] = _node_id[_graph.runningNode(a)];
358 358
        }
359 359
        _forward[j] = false;
360 360
        _source[j] = i;
361 361
        _target[j] = _root;
362 362
        _reverse[j] = k;
363 363
        _forward[k] = true;
364 364
        _source[k] = _root;
365 365
        _target[k] = i;
366 366
        _reverse[k] = j;
367 367
        ++j; ++k;
368 368
      }
369 369
      _first_out[i] = j;
370 370
      _first_out[_node_num] = k;
371 371
      for (ArcIt a(_graph); a != INVALID; ++a) {
372 372
        int fi = _arc_idf[a];
373 373
        int bi = _arc_idb[a];
374 374
        _reverse[fi] = bi;
375 375
        _reverse[bi] = fi;
376 376
      }
377 377
      
378 378
      // Reset parameters
379 379
      reset();
380 380
    }
381 381

	
382 382
    /// \name Parameters
383 383
    /// The parameters of the algorithm can be specified using these
384 384
    /// functions.
385 385

	
386 386
    /// @{
387 387

	
388 388
    /// \brief Set the lower bounds on the arcs.
389 389
    ///
390 390
    /// This function sets the lower bounds on the arcs.
391 391
    /// If it is not used before calling \ref run(), the lower bounds
392 392
    /// will be set to zero on all arcs.
393 393
    ///
394 394
    /// \param map An arc map storing the lower bounds.
395 395
    /// Its \c Value type must be convertible to the \c Value type
396 396
    /// of the algorithm.
397 397
    ///
398 398
    /// \return <tt>(*this)</tt>
399 399
    template <typename LowerMap>
400 400
    CapacityScaling& lowerMap(const LowerMap& map) {
401 401
      _have_lower = true;
402 402
      for (ArcIt a(_graph); a != INVALID; ++a) {
403 403
        _lower[_arc_idf[a]] = map[a];
404 404
        _lower[_arc_idb[a]] = map[a];
405 405
      }
406 406
      return *this;
407 407
    }
408 408

	
409 409
    /// \brief Set the upper bounds (capacities) on the arcs.
410 410
    ///
411 411
    /// This function sets the upper bounds (capacities) on the arcs.
412 412
    /// If it is not used before calling \ref run(), the upper bounds
413 413
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
414
    /// unbounded from above on each arc).
414
    /// unbounded from above).
415 415
    ///
416 416
    /// \param map An arc map storing the upper bounds.
417 417
    /// Its \c Value type must be convertible to the \c Value type
418 418
    /// of the algorithm.
419 419
    ///
420 420
    /// \return <tt>(*this)</tt>
421 421
    template<typename UpperMap>
422 422
    CapacityScaling& upperMap(const UpperMap& map) {
423 423
      for (ArcIt a(_graph); a != INVALID; ++a) {
424 424
        _upper[_arc_idf[a]] = map[a];
425 425
      }
426 426
      return *this;
427 427
    }
428 428

	
429 429
    /// \brief Set the costs of the arcs.
430 430
    ///
431 431
    /// This function sets the costs of the arcs.
432 432
    /// If it is not used before calling \ref run(), the costs
433 433
    /// will be set to \c 1 on all arcs.
434 434
    ///
435 435
    /// \param map An arc map storing the costs.
436 436
    /// Its \c Value type must be convertible to the \c Cost type
437 437
    /// of the algorithm.
438 438
    ///
439 439
    /// \return <tt>(*this)</tt>
440 440
    template<typename CostMap>
441 441
    CapacityScaling& costMap(const CostMap& map) {
442 442
      for (ArcIt a(_graph); a != INVALID; ++a) {
443 443
        _cost[_arc_idf[a]] =  map[a];
444 444
        _cost[_arc_idb[a]] = -map[a];
445 445
      }
446 446
      return *this;
447 447
    }
448 448

	
449 449
    /// \brief Set the supply values of the nodes.
450 450
    ///
451 451
    /// This function sets the supply values of the nodes.
452 452
    /// If neither this function nor \ref stSupply() is used before
453 453
    /// calling \ref run(), the supply of each node will be set to zero.
454 454
    ///
455 455
    /// \param map A node map storing the supply values.
456 456
    /// Its \c Value type must be convertible to the \c Value type
457 457
    /// of the algorithm.
458 458
    ///
459 459
    /// \return <tt>(*this)</tt>
460 460
    template<typename SupplyMap>
461 461
    CapacityScaling& supplyMap(const SupplyMap& map) {
462 462
      for (NodeIt n(_graph); n != INVALID; ++n) {
463 463
        _supply[_node_id[n]] = map[n];
464 464
      }
465 465
      return *this;
466 466
    }
467 467

	
468 468
    /// \brief Set single source and target nodes and a supply value.
469 469
    ///
470 470
    /// This function sets a single source node and a single target node
471 471
    /// and the required flow value.
472 472
    /// If neither this function nor \ref supplyMap() is used before
473 473
    /// calling \ref run(), the supply of each node will be set to zero.
474 474
    ///
475 475
    /// Using this function has the same effect as using \ref supplyMap()
476 476
    /// with such a map in which \c k is assigned to \c s, \c -k is
477 477
    /// assigned to \c t and all other nodes have zero supply value.
478 478
    ///
479 479
    /// \param s The source node.
480 480
    /// \param t The target node.
481 481
    /// \param k The required amount of flow from node \c s to node \c t
482 482
    /// (i.e. the supply of \c s and the demand of \c t).
483 483
    ///
484 484
    /// \return <tt>(*this)</tt>
485 485
    CapacityScaling& stSupply(const Node& s, const Node& t, Value k) {
486 486
      for (int i = 0; i != _node_num; ++i) {
487 487
        _supply[i] = 0;
488 488
      }
489 489
      _supply[_node_id[s]] =  k;
490 490
      _supply[_node_id[t]] = -k;
491 491
      return *this;
492 492
    }
493 493
    
494 494
    /// @}
495 495

	
496 496
    /// \name Execution control
497 497
    /// The algorithm can be executed using \ref run().
498 498

	
499 499
    /// @{
500 500

	
501 501
    /// \brief Run the algorithm.
502 502
    ///
503 503
    /// This function runs the algorithm.
504 504
    /// The paramters can be specified using functions \ref lowerMap(),
505 505
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
506 506
    /// For example,
507 507
    /// \code
508 508
    ///   CapacityScaling<ListDigraph> cs(graph);
509 509
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
510 510
    ///     .supplyMap(sup).run();
511 511
    /// \endcode
512 512
    ///
513 513
    /// This function can be called more than once. All the parameters
514 514
    /// that have been given are kept for the next call, unless
515 515
    /// \ref reset() is called, thus only the modified parameters
516 516
    /// have to be set again. See \ref reset() for examples.
517
    /// However the underlying digraph must not be modified after this
517
    /// However, the underlying digraph must not be modified after this
518 518
    /// class have been constructed, since it copies and extends the graph.
519 519
    ///
520 520
    /// \param factor The capacity scaling factor. It must be larger than
521 521
    /// one to use scaling. If it is less or equal to one, then scaling
522 522
    /// will be disabled.
523 523
    ///
524 524
    /// \return \c INFEASIBLE if no feasible flow exists,
525 525
    /// \n \c OPTIMAL if the problem has optimal solution
526 526
    /// (i.e. it is feasible and bounded), and the algorithm has found
527 527
    /// optimal flow and node potentials (primal and dual solutions),
528 528
    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
529 529
    /// and infinite upper bound. It means that the objective function
530
    /// is unbounded on that arc, however note that it could actually be
530
    /// is unbounded on that arc, however, note that it could actually be
531 531
    /// bounded over the feasible flows, but this algroithm cannot handle
532 532
    /// these cases.
533 533
    ///
534 534
    /// \see ProblemType
535 535
    ProblemType run(int factor = 4) {
536 536
      _factor = factor;
537 537
      ProblemType pt = init();
538 538
      if (pt != OPTIMAL) return pt;
539 539
      return start();
540 540
    }
541 541

	
542 542
    /// \brief Reset all the parameters that have been given before.
543 543
    ///
544 544
    /// This function resets all the paramaters that have been given
545 545
    /// before using functions \ref lowerMap(), \ref upperMap(),
546 546
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
547 547
    ///
548 548
    /// It is useful for multiple run() calls. If this function is not
549 549
    /// used, all the parameters given before are kept for the next
550 550
    /// \ref run() call.
551 551
    /// However, the underlying digraph must not be modified after this
552 552
    /// class have been constructed, since it copies and extends the graph.
553 553
    ///
554 554
    /// For example,
555 555
    /// \code
556 556
    ///   CapacityScaling<ListDigraph> cs(graph);
557 557
    ///
558 558
    ///   // First run
559 559
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
560 560
    ///     .supplyMap(sup).run();
561 561
    ///
562 562
    ///   // Run again with modified cost map (reset() is not called,
563 563
    ///   // so only the cost map have to be set again)
564 564
    ///   cost[e] += 100;
565 565
    ///   cs.costMap(cost).run();
566 566
    ///
567 567
    ///   // Run again from scratch using reset()
568 568
    ///   // (the lower bounds will be set to zero on all arcs)
569 569
    ///   cs.reset();
570 570
    ///   cs.upperMap(capacity).costMap(cost)
571 571
    ///     .supplyMap(sup).run();
572 572
    /// \endcode
573 573
    ///
574 574
    /// \return <tt>(*this)</tt>
575 575
    CapacityScaling& reset() {
576 576
      for (int i = 0; i != _node_num; ++i) {
577 577
        _supply[i] = 0;
578 578
      }
579 579
      for (int j = 0; j != _res_arc_num; ++j) {
580 580
        _lower[j] = 0;
581 581
        _upper[j] = INF;
582 582
        _cost[j] = _forward[j] ? 1 : -1;
583 583
      }
584 584
      _have_lower = false;
585 585
      return *this;
586 586
    }
587 587

	
588 588
    /// @}
589 589

	
590 590
    /// \name Query Functions
591 591
    /// The results of the algorithm can be obtained using these
592 592
    /// functions.\n
593 593
    /// The \ref run() function must be called before using them.
594 594

	
595 595
    /// @{
596 596

	
597 597
    /// \brief Return the total cost of the found flow.
598 598
    ///
599 599
    /// This function returns the total cost of the found flow.
600 600
    /// Its complexity is O(e).
601 601
    ///
602 602
    /// \note The return type of the function can be specified as a
603 603
    /// template parameter. For example,
604 604
    /// \code
605 605
    ///   cs.totalCost<double>();
606 606
    /// \endcode
607 607
    /// It is useful if the total cost cannot be stored in the \c Cost
608 608
    /// type of the algorithm, which is the default return type of the
609 609
    /// function.
610 610
    ///
611 611
    /// \pre \ref run() must be called before using this function.
612 612
    template <typename Number>
613 613
    Number totalCost() const {
614 614
      Number c = 0;
615 615
      for (ArcIt a(_graph); a != INVALID; ++a) {
616 616
        int i = _arc_idb[a];
617 617
        c += static_cast<Number>(_res_cap[i]) *
618 618
             (-static_cast<Number>(_cost[i]));
619 619
      }
620 620
      return c;
621 621
    }
622 622

	
623 623
#ifndef DOXYGEN
624 624
    Cost totalCost() const {
625 625
      return totalCost<Cost>();
626 626
    }
627 627
#endif
628 628

	
629 629
    /// \brief Return the flow on the given arc.
630 630
    ///
631 631
    /// This function returns the flow on the given arc.
632 632
    ///
633 633
    /// \pre \ref run() must be called before using this function.
634 634
    Value flow(const Arc& a) const {
635 635
      return _res_cap[_arc_idb[a]];
636 636
    }
637 637

	
638 638
    /// \brief Return the flow map (the primal solution).
639 639
    ///
640 640
    /// This function copies the flow value on each arc into the given
641 641
    /// map. The \c Value type of the algorithm must be convertible to
642 642
    /// the \c Value type of the map.
643 643
    ///
644 644
    /// \pre \ref run() must be called before using this function.
645 645
    template <typename FlowMap>
646 646
    void flowMap(FlowMap &map) const {
647 647
      for (ArcIt a(_graph); a != INVALID; ++a) {
648 648
        map.set(a, _res_cap[_arc_idb[a]]);
649 649
      }
650 650
    }
651 651

	
652 652
    /// \brief Return the potential (dual value) of the given node.
653 653
    ///
654 654
    /// This function returns the potential (dual value) of the
655 655
    /// given node.
656 656
    ///
657 657
    /// \pre \ref run() must be called before using this function.
658 658
    Cost potential(const Node& n) const {
659 659
      return _pi[_node_id[n]];
660 660
    }
661 661

	
662 662
    /// \brief Return the potential map (the dual solution).
663 663
    ///
664 664
    /// This function copies the potential (dual value) of each node
665 665
    /// into the given map.
666 666
    /// The \c Cost type of the algorithm must be convertible to the
667 667
    /// \c Value type of the map.
668 668
    ///
669 669
    /// \pre \ref run() must be called before using this function.
670 670
    template <typename PotentialMap>
671 671
    void potentialMap(PotentialMap &map) const {
672 672
      for (NodeIt n(_graph); n != INVALID; ++n) {
673 673
        map.set(n, _pi[_node_id[n]]);
674 674
      }
675 675
    }
676 676

	
677 677
    /// @}
678 678

	
679 679
  private:
680 680

	
681 681
    // Initialize the algorithm
682 682
    ProblemType init() {
683 683
      if (_node_num == 0) return INFEASIBLE;
684 684

	
685 685
      // Check the sum of supply values
686 686
      _sum_supply = 0;
687 687
      for (int i = 0; i != _root; ++i) {
688 688
        _sum_supply += _supply[i];
689 689
      }
690 690
      if (_sum_supply > 0) return INFEASIBLE;
691 691
      
692 692
      // Initialize vectors
693 693
      for (int i = 0; i != _root; ++i) {
694 694
        _pi[i] = 0;
695 695
        _excess[i] = _supply[i];
696 696
      }
697 697

	
698 698
      // Remove non-zero lower bounds
699 699
      const Value MAX = std::numeric_limits<Value>::max();
700 700
      int last_out;
701 701
      if (_have_lower) {
702 702
        for (int i = 0; i != _root; ++i) {
703 703
          last_out = _first_out[i+1];
704 704
          for (int j = _first_out[i]; j != last_out; ++j) {
705 705
            if (_forward[j]) {
706 706
              Value c = _lower[j];
707 707
              if (c >= 0) {
708 708
                _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
709 709
              } else {
710 710
                _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
711 711
              }
712 712
              _excess[i] -= c;
713 713
              _excess[_target[j]] += c;
714 714
            } else {
715 715
              _res_cap[j] = 0;
716 716
            }
717 717
          }
718 718
        }
719 719
      } else {
720 720
        for (int j = 0; j != _res_arc_num; ++j) {
721 721
          _res_cap[j] = _forward[j] ? _upper[j] : 0;
722 722
        }
Ignore white space 384 line context
1 1
/* -*- C++ -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library
4 4
 *
5 5
 * Copyright (C) 2003-2008
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
  /// \tparam V The value type used for flow amounts, capacity bounds
43
  /// \tparam V The number type used for flow amounts, capacity bounds
44 44
  /// and supply values. By default it is \c int.
45
  /// \tparam C The value type used for costs and potentials.
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 minimum cost
94 94
  /// flow. It is an efficient primal-dual solution method, which
95 95
  /// can be viewed as the generalization of the \ref Preflow
96 96
  /// "preflow push-relabel" algorithm for the maximum flow problem.
97 97
  ///
98 98
  /// Most of the parameters of the problem (except for the digraph)
99 99
  /// can be given using separate functions, and the algorithm can be
100 100
  /// executed using the \ref run() function. If some parameters are not
101 101
  /// specified, then default values will be used.
102 102
  ///
103 103
  /// \tparam GR The digraph type the algorithm runs on.
104
  /// \tparam V The value type used for flow amounts, capacity bounds
104
  /// \tparam V The number type used for flow amounts, capacity bounds
105 105
  /// and supply values in the algorithm. By default it is \c int.
106
  /// \tparam C The value type used for costs and potentials in the
106
  /// \tparam C The number type used for costs and potentials in the
107 107
  /// algorithm. By default it is the same as \c V.
108 108
  ///
109
  /// \warning Both value types must be signed and all input data must
109
  /// \warning Both number types must be signed and all input data must
110 110
  /// be integer.
111 111
  /// \warning This algorithm does not support negative costs for such
112 112
  /// arcs that have infinite upper bound.
113 113
  ///
114 114
  /// \note %CostScaling provides three different internal methods,
115 115
  /// from which the most efficient one is used by default.
116 116
  /// For more information, see \ref Method.
117 117
#ifdef DOXYGEN
118 118
  template <typename GR, typename V, typename C, typename TR>
119 119
#else
120 120
  template < typename GR, typename V = int, typename C = V,
121 121
             typename TR = CostScalingDefaultTraits<GR, V, C> >
122 122
#endif
123 123
  class CostScaling
124 124
  {
125 125
  public:
126 126

	
127 127
    /// The type of the digraph
128 128
    typedef typename TR::Digraph Digraph;
129 129
    /// The type of the flow amounts, capacity bounds and supply values
130 130
    typedef typename TR::Value Value;
131 131
    /// The type of the arc costs
132 132
    typedef typename TR::Cost Cost;
133 133

	
134 134
    /// \brief The large cost type
135 135
    ///
136 136
    /// The large cost type used for internal computations.
137 137
    /// Using the \ref CostScalingDefaultTraits "default traits class",
138 138
    /// it is \c long \c long if the \c Cost type is integer,
139 139
    /// otherwise it is \c double.
140 140
    typedef typename TR::LargeCost LargeCost;
141 141

	
142 142
    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
143 143
    typedef TR Traits;
144 144

	
145 145
  public:
146 146

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

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

	
193 193
  private:
194 194

	
195 195
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
196 196

	
197 197
    typedef std::vector<int> IntVector;
198 198
    typedef std::vector<char> BoolVector;
199 199
    typedef std::vector<Value> ValueVector;
200 200
    typedef std::vector<Cost> CostVector;
201 201
    typedef std::vector<LargeCost> LargeCostVector;
202 202

	
203 203
  private:
204 204
  
205 205
    template <typename KT, typename VT>
206 206
    class VectorMap {
207 207
    public:
208 208
      typedef KT Key;
209 209
      typedef VT Value;
210 210
      
211 211
      VectorMap(std::vector<Value>& v) : _v(v) {}
212 212
      
213 213
      const Value& operator[](const Key& key) const {
214 214
        return _v[StaticDigraph::id(key)];
215 215
      }
216 216

	
217 217
      Value& operator[](const Key& key) {
218 218
        return _v[StaticDigraph::id(key)];
219 219
      }
220 220
      
221 221
      void set(const Key& key, const Value& val) {
222 222
        _v[StaticDigraph::id(key)] = val;
223 223
      }
224 224

	
225 225
    private:
226 226
      std::vector<Value>& _v;
227 227
    };
228 228

	
229 229
    typedef VectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
230 230
    typedef VectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
231 231

	
232 232
  private:
233 233

	
234 234
    // Data related to the underlying digraph
235 235
    const GR &_graph;
236 236
    int _node_num;
237 237
    int _arc_num;
238 238
    int _res_node_num;
239 239
    int _res_arc_num;
240 240
    int _root;
241 241

	
242 242
    // Parameters of the problem
243 243
    bool _have_lower;
244 244
    Value _sum_supply;
245 245

	
246 246
    // Data structures for storing the digraph
247 247
    IntNodeMap _node_id;
248 248
    IntArcMap _arc_idf;
249 249
    IntArcMap _arc_idb;
250 250
    IntVector _first_out;
251 251
    BoolVector _forward;
252 252
    IntVector _source;
253 253
    IntVector _target;
254 254
    IntVector _reverse;
255 255

	
256 256
    // Node and arc data
257 257
    ValueVector _lower;
258 258
    ValueVector _upper;
259 259
    CostVector _scost;
260 260
    ValueVector _supply;
261 261

	
262 262
    ValueVector _res_cap;
263 263
    LargeCostVector _cost;
264 264
    LargeCostVector _pi;
265 265
    ValueVector _excess;
266 266
    IntVector _next_out;
267 267
    std::deque<int> _active_nodes;
268 268

	
269 269
    // Data for scaling
270 270
    LargeCost _epsilon;
271 271
    int _alpha;
272 272

	
273 273
    // Data for a StaticDigraph structure
274 274
    typedef std::pair<int, int> IntPair;
275 275
    StaticDigraph _sgr;
276 276
    std::vector<IntPair> _arc_vec;
277 277
    std::vector<LargeCost> _cost_vec;
278 278
    LargeCostArcMap _cost_map;
279 279
    LargeCostNodeMap _pi_map;
280 280
  
281 281
  public:
282 282
  
283 283
    /// \brief Constant for infinite upper bounds (capacities).
284 284
    ///
285 285
    /// Constant for infinite upper bounds (capacities).
286 286
    /// It is \c std::numeric_limits<Value>::infinity() if available,
287 287
    /// \c std::numeric_limits<Value>::max() otherwise.
288 288
    const Value INF;
289 289

	
290 290
  public:
291 291

	
292 292
    /// \name Named Template Parameters
293 293
    /// @{
294 294

	
295 295
    template <typename T>
296 296
    struct SetLargeCostTraits : public Traits {
297 297
      typedef T LargeCost;
298 298
    };
299 299

	
300 300
    /// \brief \ref named-templ-param "Named parameter" for setting
301 301
    /// \c LargeCost type.
302 302
    ///
303 303
    /// \ref named-templ-param "Named parameter" for setting \c LargeCost
304 304
    /// type, which is used for internal computations in the algorithm.
305 305
    /// \c Cost must be convertible to \c LargeCost.
306 306
    template <typename T>
307 307
    struct SetLargeCost
308 308
      : public CostScaling<GR, V, C, SetLargeCostTraits<T> > {
309 309
      typedef  CostScaling<GR, V, C, SetLargeCostTraits<T> > Create;
310 310
    };
311 311

	
312 312
    /// @}
313 313

	
314 314
  public:
315 315

	
316 316
    /// \brief Constructor.
317 317
    ///
318 318
    /// The constructor of the class.
319 319
    ///
320 320
    /// \param graph The digraph the algorithm runs on.
321 321
    CostScaling(const GR& graph) :
322 322
      _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
323 323
      _cost_map(_cost_vec), _pi_map(_pi),
324 324
      INF(std::numeric_limits<Value>::has_infinity ?
325 325
          std::numeric_limits<Value>::infinity() :
326 326
          std::numeric_limits<Value>::max())
327 327
    {
328
      // Check the value types
328
      // Check the number types
329 329
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
330 330
        "The flow type of CostScaling must be signed");
331 331
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
332 332
        "The cost type of CostScaling must be signed");
333 333

	
334 334
      // Resize vectors
335 335
      _node_num = countNodes(_graph);
336 336
      _arc_num = countArcs(_graph);
337 337
      _res_node_num = _node_num + 1;
338 338
      _res_arc_num = 2 * (_arc_num + _node_num);
339 339
      _root = _node_num;
340 340

	
341 341
      _first_out.resize(_res_node_num + 1);
342 342
      _forward.resize(_res_arc_num);
343 343
      _source.resize(_res_arc_num);
344 344
      _target.resize(_res_arc_num);
345 345
      _reverse.resize(_res_arc_num);
346 346

	
347 347
      _lower.resize(_res_arc_num);
348 348
      _upper.resize(_res_arc_num);
349 349
      _scost.resize(_res_arc_num);
350 350
      _supply.resize(_res_node_num);
351 351
      
352 352
      _res_cap.resize(_res_arc_num);
353 353
      _cost.resize(_res_arc_num);
354 354
      _pi.resize(_res_node_num);
355 355
      _excess.resize(_res_node_num);
356 356
      _next_out.resize(_res_node_num);
357 357

	
358 358
      _arc_vec.reserve(_res_arc_num);
359 359
      _cost_vec.reserve(_res_arc_num);
360 360

	
361 361
      // Copy the graph
362 362
      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
363 363
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
364 364
        _node_id[n] = i;
365 365
      }
366 366
      i = 0;
367 367
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
368 368
        _first_out[i] = j;
369 369
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
370 370
          _arc_idf[a] = j;
371 371
          _forward[j] = true;
372 372
          _source[j] = i;
373 373
          _target[j] = _node_id[_graph.runningNode(a)];
374 374
        }
375 375
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
376 376
          _arc_idb[a] = j;
377 377
          _forward[j] = false;
378 378
          _source[j] = i;
379 379
          _target[j] = _node_id[_graph.runningNode(a)];
380 380
        }
381 381
        _forward[j] = false;
382 382
        _source[j] = i;
383 383
        _target[j] = _root;
384 384
        _reverse[j] = k;
385 385
        _forward[k] = true;
386 386
        _source[k] = _root;
387 387
        _target[k] = i;
388 388
        _reverse[k] = j;
389 389
        ++j; ++k;
390 390
      }
391 391
      _first_out[i] = j;
392 392
      _first_out[_res_node_num] = k;
393 393
      for (ArcIt a(_graph); a != INVALID; ++a) {
394 394
        int fi = _arc_idf[a];
395 395
        int bi = _arc_idb[a];
396 396
        _reverse[fi] = bi;
397 397
        _reverse[bi] = fi;
398 398
      }
399 399
      
400 400
      // Reset parameters
401 401
      reset();
402 402
    }
403 403

	
404 404
    /// \name Parameters
405 405
    /// The parameters of the algorithm can be specified using these
406 406
    /// functions.
407 407

	
408 408
    /// @{
409 409

	
410 410
    /// \brief Set the lower bounds on the arcs.
411 411
    ///
412 412
    /// This function sets the lower bounds on the arcs.
413 413
    /// If it is not used before calling \ref run(), the lower bounds
414 414
    /// will be set to zero on all arcs.
415 415
    ///
416 416
    /// \param map An arc map storing the lower bounds.
417 417
    /// Its \c Value type must be convertible to the \c Value type
418 418
    /// of the algorithm.
419 419
    ///
420 420
    /// \return <tt>(*this)</tt>
421 421
    template <typename LowerMap>
422 422
    CostScaling& lowerMap(const LowerMap& map) {
423 423
      _have_lower = true;
424 424
      for (ArcIt a(_graph); a != INVALID; ++a) {
425 425
        _lower[_arc_idf[a]] = map[a];
426 426
        _lower[_arc_idb[a]] = map[a];
427 427
      }
428 428
      return *this;
429 429
    }
430 430

	
431 431
    /// \brief Set the upper bounds (capacities) on the arcs.
432 432
    ///
433 433
    /// This function sets the upper bounds (capacities) on the arcs.
434 434
    /// If it is not used before calling \ref run(), the upper bounds
435 435
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
436
    /// unbounded from above on each arc).
436
    /// unbounded from above).
437 437
    ///
438 438
    /// \param map An arc map storing the upper bounds.
439 439
    /// Its \c Value type must be convertible to the \c Value type
440 440
    /// of the algorithm.
441 441
    ///
442 442
    /// \return <tt>(*this)</tt>
443 443
    template<typename UpperMap>
444 444
    CostScaling& upperMap(const UpperMap& map) {
445 445
      for (ArcIt a(_graph); a != INVALID; ++a) {
446 446
        _upper[_arc_idf[a]] = map[a];
447 447
      }
448 448
      return *this;
449 449
    }
450 450

	
451 451
    /// \brief Set the costs of the arcs.
452 452
    ///
453 453
    /// This function sets the costs of the arcs.
454 454
    /// If it is not used before calling \ref run(), the costs
455 455
    /// will be set to \c 1 on all arcs.
456 456
    ///
457 457
    /// \param map An arc map storing the costs.
458 458
    /// Its \c Value type must be convertible to the \c Cost type
459 459
    /// of the algorithm.
460 460
    ///
461 461
    /// \return <tt>(*this)</tt>
462 462
    template<typename CostMap>
463 463
    CostScaling& costMap(const CostMap& map) {
464 464
      for (ArcIt a(_graph); a != INVALID; ++a) {
465 465
        _scost[_arc_idf[a]] =  map[a];
466 466
        _scost[_arc_idb[a]] = -map[a];
467 467
      }
468 468
      return *this;
469 469
    }
470 470

	
471 471
    /// \brief Set the supply values of the nodes.
472 472
    ///
473 473
    /// This function sets the supply values of the nodes.
474 474
    /// If neither this function nor \ref stSupply() is used before
475 475
    /// calling \ref run(), the supply of each node will be set to zero.
476 476
    ///
477 477
    /// \param map A node map storing the supply values.
478 478
    /// Its \c Value type must be convertible to the \c Value type
479 479
    /// of the algorithm.
480 480
    ///
481 481
    /// \return <tt>(*this)</tt>
482 482
    template<typename SupplyMap>
483 483
    CostScaling& supplyMap(const SupplyMap& map) {
484 484
      for (NodeIt n(_graph); n != INVALID; ++n) {
485 485
        _supply[_node_id[n]] = map[n];
486 486
      }
487 487
      return *this;
488 488
    }
489 489

	
490 490
    /// \brief Set single source and target nodes and a supply value.
491 491
    ///
492 492
    /// This function sets a single source node and a single target node
493 493
    /// and the required flow value.
494 494
    /// If neither this function nor \ref supplyMap() is used before
495 495
    /// calling \ref run(), the supply of each node will be set to zero.
496 496
    ///
497 497
    /// Using this function has the same effect as using \ref supplyMap()
498 498
    /// with such a map in which \c k is assigned to \c s, \c -k is
499 499
    /// assigned to \c t and all other nodes have zero supply value.
500 500
    ///
501 501
    /// \param s The source node.
502 502
    /// \param t The target node.
503 503
    /// \param k The required amount of flow from node \c s to node \c t
504 504
    /// (i.e. the supply of \c s and the demand of \c t).
505 505
    ///
506 506
    /// \return <tt>(*this)</tt>
507 507
    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
508 508
      for (int i = 0; i != _res_node_num; ++i) {
509 509
        _supply[i] = 0;
510 510
      }
511 511
      _supply[_node_id[s]] =  k;
512 512
      _supply[_node_id[t]] = -k;
513 513
      return *this;
514 514
    }
515 515
    
516 516
    /// @}
517 517

	
518 518
    /// \name Execution control
519 519
    /// The algorithm can be executed using \ref run().
520 520

	
521 521
    /// @{
522 522

	
523 523
    /// \brief Run the algorithm.
524 524
    ///
525 525
    /// This function runs the algorithm.
526 526
    /// The paramters can be specified using functions \ref lowerMap(),
527 527
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
528 528
    /// For example,
529 529
    /// \code
530 530
    ///   CostScaling<ListDigraph> cs(graph);
531 531
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
532 532
    ///     .supplyMap(sup).run();
533 533
    /// \endcode
534 534
    ///
535 535
    /// This function can be called more than once. All the parameters
536 536
    /// that have been given are kept for the next call, unless
537 537
    /// \ref reset() is called, thus only the modified parameters
538 538
    /// have to be set again. See \ref reset() for examples.
539 539
    /// However, the underlying digraph must not be modified after this
540 540
    /// class have been constructed, since it copies and extends the graph.
541 541
    ///
542 542
    /// \param method The internal method that will be used in the
543 543
    /// algorithm. For more information, see \ref Method.
544 544
    /// \param factor The cost scaling factor. It must be larger than one.
545 545
    ///
546 546
    /// \return \c INFEASIBLE if no feasible flow exists,
547 547
    /// \n \c OPTIMAL if the problem has optimal solution
548 548
    /// (i.e. it is feasible and bounded), and the algorithm has found
549 549
    /// optimal flow and node potentials (primal and dual solutions),
550 550
    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
551 551
    /// and infinite upper bound. It means that the objective function
552
    /// is unbounded on that arc, however note that it could actually be
552
    /// is unbounded on that arc, however, note that it could actually be
553 553
    /// bounded over the feasible flows, but this algroithm cannot handle
554 554
    /// these cases.
555 555
    ///
556 556
    /// \see ProblemType, Method
557 557
    ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
558 558
      _alpha = factor;
559 559
      ProblemType pt = init();
560 560
      if (pt != OPTIMAL) return pt;
561 561
      start(method);
562 562
      return OPTIMAL;
563 563
    }
564 564

	
565 565
    /// \brief Reset all the parameters that have been given before.
566 566
    ///
567 567
    /// This function resets all the paramaters that have been given
568 568
    /// before using functions \ref lowerMap(), \ref upperMap(),
569 569
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
570 570
    ///
571 571
    /// It is useful for multiple run() calls. If this function is not
572 572
    /// used, all the parameters given before are kept for the next
573 573
    /// \ref run() call.
574
    /// However the underlying digraph must not be modified after this
574
    /// However, the underlying digraph must not be modified after this
575 575
    /// class have been constructed, since it copies and extends the graph.
576 576
    ///
577 577
    /// For example,
578 578
    /// \code
579 579
    ///   CostScaling<ListDigraph> cs(graph);
580 580
    ///
581 581
    ///   // First run
582 582
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
583 583
    ///     .supplyMap(sup).run();
584 584
    ///
585 585
    ///   // Run again with modified cost map (reset() is not called,
586 586
    ///   // so only the cost map have to be set again)
587 587
    ///   cost[e] += 100;
588 588
    ///   cs.costMap(cost).run();
589 589
    ///
590 590
    ///   // Run again from scratch using reset()
591 591
    ///   // (the lower bounds will be set to zero on all arcs)
592 592
    ///   cs.reset();
593 593
    ///   cs.upperMap(capacity).costMap(cost)
594 594
    ///     .supplyMap(sup).run();
595 595
    /// \endcode
596 596
    ///
597 597
    /// \return <tt>(*this)</tt>
598 598
    CostScaling& reset() {
599 599
      for (int i = 0; i != _res_node_num; ++i) {
600 600
        _supply[i] = 0;
601 601
      }
602 602
      int limit = _first_out[_root];
603 603
      for (int j = 0; j != limit; ++j) {
604 604
        _lower[j] = 0;
605 605
        _upper[j] = INF;
606 606
        _scost[j] = _forward[j] ? 1 : -1;
607 607
      }
608 608
      for (int j = limit; j != _res_arc_num; ++j) {
609 609
        _lower[j] = 0;
610 610
        _upper[j] = INF;
611 611
        _scost[j] = 0;
612 612
        _scost[_reverse[j]] = 0;
613 613
      }      
614 614
      _have_lower = false;
615 615
      return *this;
616 616
    }
617 617

	
618 618
    /// @}
619 619

	
620 620
    /// \name Query Functions
621 621
    /// The results of the algorithm can be obtained using these
622 622
    /// functions.\n
623 623
    /// The \ref run() function must be called before using them.
624 624

	
625 625
    /// @{
626 626

	
627 627
    /// \brief Return the total cost of the found flow.
628 628
    ///
629 629
    /// This function returns the total cost of the found flow.
630 630
    /// Its complexity is O(e).
631 631
    ///
632 632
    /// \note The return type of the function can be specified as a
633 633
    /// template parameter. For example,
634 634
    /// \code
635 635
    ///   cs.totalCost<double>();
636 636
    /// \endcode
637 637
    /// It is useful if the total cost cannot be stored in the \c Cost
638 638
    /// type of the algorithm, which is the default return type of the
639 639
    /// function.
640 640
    ///
641 641
    /// \pre \ref run() must be called before using this function.
642 642
    template <typename Number>
643 643
    Number totalCost() const {
644 644
      Number c = 0;
645 645
      for (ArcIt a(_graph); a != INVALID; ++a) {
646 646
        int i = _arc_idb[a];
647 647
        c += static_cast<Number>(_res_cap[i]) *
648 648
             (-static_cast<Number>(_scost[i]));
649 649
      }
650 650
      return c;
651 651
    }
652 652

	
653 653
#ifndef DOXYGEN
654 654
    Cost totalCost() const {
655 655
      return totalCost<Cost>();
656 656
    }
657 657
#endif
658 658

	
659 659
    /// \brief Return the flow on the given arc.
660 660
    ///
661 661
    /// This function returns the flow on the given arc.
662 662
    ///
663 663
    /// \pre \ref run() must be called before using this function.
664 664
    Value flow(const Arc& a) const {
665 665
      return _res_cap[_arc_idb[a]];
666 666
    }
667 667

	
668 668
    /// \brief Return the flow map (the primal solution).
669 669
    ///
670 670
    /// This function copies the flow value on each arc into the given
671 671
    /// map. The \c Value type of the algorithm must be convertible to
672 672
    /// the \c Value type of the map.
673 673
    ///
674 674
    /// \pre \ref run() must be called before using this function.
675 675
    template <typename FlowMap>
676 676
    void flowMap(FlowMap &map) const {
677 677
      for (ArcIt a(_graph); a != INVALID; ++a) {
678 678
        map.set(a, _res_cap[_arc_idb[a]]);
679 679
      }
680 680
    }
681 681

	
682 682
    /// \brief Return the potential (dual value) of the given node.
683 683
    ///
684 684
    /// This function returns the potential (dual value) of the
685 685
    /// given node.
686 686
    ///
687 687
    /// \pre \ref run() must be called before using this function.
688 688
    Cost potential(const Node& n) const {
689 689
      return static_cast<Cost>(_pi[_node_id[n]]);
690 690
    }
691 691

	
692 692
    /// \brief Return the potential map (the dual solution).
693 693
    ///
694 694
    /// This function copies the potential (dual value) of each node
695 695
    /// into the given map.
696 696
    /// The \c Cost type of the algorithm must be convertible to the
697 697
    /// \c Value type of the map.
698 698
    ///
699 699
    /// \pre \ref run() must be called before using this function.
700 700
    template <typename PotentialMap>
701 701
    void potentialMap(PotentialMap &map) const {
702 702
      for (NodeIt n(_graph); n != INVALID; ++n) {
703 703
        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
704 704
      }
705 705
    }
706 706

	
707 707
    /// @}
708 708

	
709 709
  private:
710 710

	
711 711
    // Initialize the algorithm
712 712
    ProblemType init() {
713 713
      if (_res_node_num == 0) return INFEASIBLE;
714 714

	
715 715
      // Check the sum of supply values
716 716
      _sum_supply = 0;
717 717
      for (int i = 0; i != _root; ++i) {
718 718
        _sum_supply += _supply[i];
719 719
      }
720 720
      if (_sum_supply > 0) return INFEASIBLE;
721 721
      
722 722

	
723 723
      // Initialize vectors
724 724
      for (int i = 0; i != _res_node_num; ++i) {
725 725
        _pi[i] = 0;
726 726
        _excess[i] = _supply[i];
727 727
      }
728 728
      
729 729
      // Remove infinite upper bounds and check negative arcs
730 730
      const Value MAX = std::numeric_limits<Value>::max();
731 731
      int last_out;
732 732
      if (_have_lower) {
733 733
        for (int i = 0; i != _root; ++i) {
734 734
          last_out = _first_out[i+1];
735 735
          for (int j = _first_out[i]; j != last_out; ++j) {
736 736
            if (_forward[j]) {
737 737
              Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
738 738
              if (c >= MAX) return UNBOUNDED;
739 739
              _excess[i] -= c;
740 740
              _excess[_target[j]] += c;
741 741
            }
742 742
          }
743 743
        }
744 744
      } else {
745 745
        for (int i = 0; i != _root; ++i) {
746 746
          last_out = _first_out[i+1];
747 747
          for (int j = _first_out[i]; j != last_out; ++j) {
748 748
            if (_forward[j] && _scost[j] < 0) {
749 749
              Value c = _upper[j];
750 750
              if (c >= MAX) return UNBOUNDED;
751 751
              _excess[i] -= c;
752 752
              _excess[_target[j]] += c;
753 753
            }
754 754
          }
755 755
        }
756 756
      }
757 757
      Value ex, max_cap = 0;
758 758
      for (int i = 0; i != _res_node_num; ++i) {
759 759
        ex = _excess[i];
760 760
        _excess[i] = 0;
761 761
        if (ex < 0) max_cap -= ex;
762 762
      }
763 763
      for (int j = 0; j != _res_arc_num; ++j) {
764 764
        if (_upper[j] >= MAX) _upper[j] = max_cap;
765 765
      }
766 766

	
Ignore white space 384 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_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
  /// This algorithm is a specialized version of the linear programming
47
  /// simplex method directly for the minimum cost flow problem.
48
  /// It is one of the most efficient solution methods.
46
  /// This algorithm is a highly efficient specialized version of the
47
  /// linear programming simplex method directly for the minimum cost
48
  /// flow problem.
49 49
  ///
50
  /// In general this class is the fastest implementation available
51
  /// in LEMON for the minimum cost flow problem.
52
  /// Moreover it supports both directions of the supply/demand inequality
50
  /// In general, %NetworkSimplex is the fastest implementation available
51
  /// in LEMON for this problem.
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
  /// \tparam V The value type used for flow amounts, capacity bounds
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
  /// \tparam C The value type used for costs and potentials in the
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 value types must be signed and all input data must
66
  /// \warning Both number types must be signed and all input data must
67 67
  /// be integer.
68 68
  ///
69 69
  /// \note %NetworkSimplex provides five different pivot rule
70 70
  /// implementations, from which the most efficient one is used
71 71
  /// by default. For more information, see \ref PivotRule.
72 72
  template <typename GR, typename V = int, typename C = V>
73 73
  class NetworkSimplex
74 74
  {
75 75
  public:
76 76

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

	
82 82
  public:
83 83

	
84 84
    /// \brief Problem type constants for the \c run() function.
85 85
    ///
86 86
    /// Enum type containing the problem type constants that can be
87 87
    /// returned by the \ref run() function of the algorithm.
88 88
    enum ProblemType {
89 89
      /// The problem has no feasible solution (flow).
90 90
      INFEASIBLE,
91 91
      /// The problem has optimal solution (i.e. it is feasible and
92 92
      /// bounded), and the algorithm has found optimal flow and node
93 93
      /// potentials (primal and dual solutions).
94 94
      OPTIMAL,
95 95
      /// The objective function of the problem is unbounded, i.e.
96 96
      /// there is a directed cycle having negative total cost and
97 97
      /// infinite upper bound.
98 98
      UNBOUNDED
99 99
    };
100 100
    
101 101
    /// \brief Constants for selecting the type of the supply constraints.
102 102
    ///
103 103
    /// Enum type containing constants for selecting the supply type,
104 104
    /// i.e. the direction of the inequalities in the supply/demand
105 105
    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
106 106
    ///
107 107
    /// The default supply type is \c GEQ, the \c LEQ type can be
108 108
    /// selected using \ref supplyType().
109 109
    /// The equality form is a special case of both supply types.
110 110
    enum SupplyType {
111 111
      /// This option means that there are <em>"greater or equal"</em>
112 112
      /// supply/demand constraints in the definition of the problem.
113 113
      GEQ,
114 114
      /// This option means that there are <em>"less or equal"</em>
115 115
      /// supply/demand constraints in the definition of the problem.
116 116
      LEQ
117 117
    };
118 118
    
119 119
    /// \brief Constants for selecting the pivot rule.
120 120
    ///
121 121
    /// Enum type containing constants for selecting the pivot rule for
122 122
    /// the \ref run() function.
123 123
    ///
124 124
    /// \ref NetworkSimplex provides five different pivot rule
125 125
    /// implementations that significantly affect the running time
126 126
    /// of the algorithm.
127 127
    /// By default, \ref BLOCK_SEARCH "Block Search" is used, which
128 128
    /// proved to be the most efficient and the most robust on various
129
    /// test inputs according to our benchmark tests.
129
    /// test inputs.
130 130
    /// However, another pivot rule can be selected using the \ref run()
131 131
    /// function with the proper parameter.
132 132
    enum PivotRule {
133 133

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

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

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

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

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

	
164 164
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
165 165

	
166 166
    typedef std::vector<int> IntVector;
167 167
    typedef std::vector<char> CharVector;
168 168
    typedef std::vector<Value> ValueVector;
169 169
    typedef std::vector<Cost> CostVector;
170 170

	
171 171
    // State constants for arcs
172 172
    enum ArcStateEnum {
173 173
      STATE_UPPER = -1,
174 174
      STATE_TREE  =  0,
175 175
      STATE_LOWER =  1
176 176
    };
177 177

	
178 178
  private:
179 179

	
180 180
    // Data related to the underlying digraph
181 181
    const GR &_graph;
182 182
    int _node_num;
183 183
    int _arc_num;
184 184
    int _all_arc_num;
185 185
    int _search_arc_num;
186 186

	
187 187
    // Parameters of the problem
188 188
    bool _have_lower;
189 189
    SupplyType _stype;
190 190
    Value _sum_supply;
191 191

	
192 192
    // Data structures for storing the digraph
193 193
    IntNodeMap _node_id;
194 194
    IntArcMap _arc_id;
195 195
    IntVector _source;
196 196
    IntVector _target;
197 197

	
198 198
    // Node and arc data
199 199
    ValueVector _lower;
200 200
    ValueVector _upper;
201 201
    ValueVector _cap;
202 202
    CostVector _cost;
203 203
    ValueVector _supply;
204 204
    ValueVector _flow;
205 205
    CostVector _pi;
206 206

	
207 207
    // Data for storing the spanning tree structure
208 208
    IntVector _parent;
209 209
    IntVector _pred;
210 210
    IntVector _thread;
211 211
    IntVector _rev_thread;
212 212
    IntVector _succ_num;
213 213
    IntVector _last_succ;
214 214
    IntVector _dirty_revs;
215 215
    CharVector _forward;
216 216
    CharVector _state;
217 217
    int _root;
218 218

	
219 219
    // Temporary data used in the current pivot iteration
220 220
    int in_arc, join, u_in, v_in, u_out, v_out;
221 221
    int first, second, right, last;
222 222
    int stem, par_stem, new_stem;
223 223
    Value delta;
224 224
    
225 225
    const Value MAX;
226 226

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

	
236 236
  private:
237 237

	
238 238
    // Implementation of the First Eligible pivot rule
239 239
    class FirstEligiblePivotRule
240 240
    {
241 241
    private:
242 242

	
243 243
      // References to the NetworkSimplex class
244 244
      const IntVector  &_source;
245 245
      const IntVector  &_target;
246 246
      const CostVector &_cost;
247 247
      const CharVector &_state;
248 248
      const CostVector &_pi;
249 249
      int &_in_arc;
250 250
      int _search_arc_num;
251 251

	
252 252
      // Pivot rule data
253 253
      int _next_arc;
254 254

	
255 255
    public:
256 256

	
257 257
      // Constructor
258 258
      FirstEligiblePivotRule(NetworkSimplex &ns) :
259 259
        _source(ns._source), _target(ns._target),
260 260
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
261 261
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
262 262
        _next_arc(0)
263 263
      {}
264 264

	
265 265
      // Find next entering arc
266 266
      bool findEnteringArc() {
267 267
        Cost c;
268 268
        for (int e = _next_arc; e < _search_arc_num; ++e) {
269 269
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
270 270
          if (c < 0) {
271 271
            _in_arc = e;
272 272
            _next_arc = e + 1;
273 273
            return true;
274 274
          }
275 275
        }
276 276
        for (int e = 0; e < _next_arc; ++e) {
277 277
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
278 278
          if (c < 0) {
279 279
            _in_arc = e;
280 280
            _next_arc = e + 1;
281 281
            return true;
282 282
          }
283 283
        }
284 284
        return false;
285 285
      }
286 286

	
287 287
    }; //class FirstEligiblePivotRule
288 288

	
289 289

	
290 290
    // Implementation of the Best Eligible pivot rule
291 291
    class BestEligiblePivotRule
292 292
    {
293 293
    private:
294 294

	
295 295
      // References to the NetworkSimplex class
296 296
      const IntVector  &_source;
297 297
      const IntVector  &_target;
298 298
      const CostVector &_cost;
299 299
      const CharVector &_state;
300 300
      const CostVector &_pi;
301 301
      int &_in_arc;
302 302
      int _search_arc_num;
303 303

	
304 304
    public:
305 305

	
306 306
      // Constructor
307 307
      BestEligiblePivotRule(NetworkSimplex &ns) :
308 308
        _source(ns._source), _target(ns._target),
309 309
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
310 310
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
311 311
      {}
312 312

	
313 313
      // Find next entering arc
314 314
      bool findEnteringArc() {
315 315
        Cost c, min = 0;
316 316
        for (int e = 0; e < _search_arc_num; ++e) {
317 317
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
318 318
          if (c < min) {
319 319
            min = c;
320 320
            _in_arc = e;
321 321
          }
... ...
@@ -448,477 +448,477 @@
448 448
        Cost min, c;
449 449
        int e;
450 450
        if (_curr_length > 0 && _minor_count < _minor_limit) {
451 451
          // Minor iteration: select the best eligible arc from the
452 452
          // current candidate list
453 453
          ++_minor_count;
454 454
          min = 0;
455 455
          for (int i = 0; i < _curr_length; ++i) {
456 456
            e = _candidates[i];
457 457
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
458 458
            if (c < min) {
459 459
              min = c;
460 460
              _in_arc = e;
461 461
            }
462 462
            else if (c >= 0) {
463 463
              _candidates[i--] = _candidates[--_curr_length];
464 464
            }
465 465
          }
466 466
          if (min < 0) return true;
467 467
        }
468 468

	
469 469
        // Major iteration: build a new candidate list
470 470
        min = 0;
471 471
        _curr_length = 0;
472 472
        for (e = _next_arc; e < _search_arc_num; ++e) {
473 473
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
474 474
          if (c < 0) {
475 475
            _candidates[_curr_length++] = e;
476 476
            if (c < min) {
477 477
              min = c;
478 478
              _in_arc = e;
479 479
            }
480 480
            if (_curr_length == _list_length) goto search_end;
481 481
          }
482 482
        }
483 483
        for (e = 0; e < _next_arc; ++e) {
484 484
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
485 485
          if (c < 0) {
486 486
            _candidates[_curr_length++] = e;
487 487
            if (c < min) {
488 488
              min = c;
489 489
              _in_arc = e;
490 490
            }
491 491
            if (_curr_length == _list_length) goto search_end;
492 492
          }
493 493
        }
494 494
        if (_curr_length == 0) return false;
495 495
      
496 496
      search_end:        
497 497
        _minor_count = 1;
498 498
        _next_arc = e;
499 499
        return true;
500 500
      }
501 501

	
502 502
    }; //class CandidateListPivotRule
503 503

	
504 504

	
505 505
    // Implementation of the Altering Candidate List pivot rule
506 506
    class AlteringListPivotRule
507 507
    {
508 508
    private:
509 509

	
510 510
      // References to the NetworkSimplex class
511 511
      const IntVector  &_source;
512 512
      const IntVector  &_target;
513 513
      const CostVector &_cost;
514 514
      const CharVector &_state;
515 515
      const CostVector &_pi;
516 516
      int &_in_arc;
517 517
      int _search_arc_num;
518 518

	
519 519
      // Pivot rule data
520 520
      int _block_size, _head_length, _curr_length;
521 521
      int _next_arc;
522 522
      IntVector _candidates;
523 523
      CostVector _cand_cost;
524 524

	
525 525
      // Functor class to compare arcs during sort of the candidate list
526 526
      class SortFunc
527 527
      {
528 528
      private:
529 529
        const CostVector &_map;
530 530
      public:
531 531
        SortFunc(const CostVector &map) : _map(map) {}
532 532
        bool operator()(int left, int right) {
533 533
          return _map[left] > _map[right];
534 534
        }
535 535
      };
536 536

	
537 537
      SortFunc _sort_func;
538 538

	
539 539
    public:
540 540

	
541 541
      // Constructor
542 542
      AlteringListPivotRule(NetworkSimplex &ns) :
543 543
        _source(ns._source), _target(ns._target),
544 544
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
545 545
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
546 546
        _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
547 547
      {
548 548
        // The main parameters of the pivot rule
549 549
        const double BLOCK_SIZE_FACTOR = 1.0;
550 550
        const int MIN_BLOCK_SIZE = 10;
551 551
        const double HEAD_LENGTH_FACTOR = 0.1;
552 552
        const int MIN_HEAD_LENGTH = 3;
553 553

	
554 554
        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
555 555
                                    std::sqrt(double(_search_arc_num))),
556 556
                                MIN_BLOCK_SIZE );
557 557
        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
558 558
                                 MIN_HEAD_LENGTH );
559 559
        _candidates.resize(_head_length + _block_size);
560 560
        _curr_length = 0;
561 561
      }
562 562

	
563 563
      // Find next entering arc
564 564
      bool findEnteringArc() {
565 565
        // Check the current candidate list
566 566
        int e;
567 567
        for (int i = 0; i < _curr_length; ++i) {
568 568
          e = _candidates[i];
569 569
          _cand_cost[e] = _state[e] *
570 570
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
571 571
          if (_cand_cost[e] >= 0) {
572 572
            _candidates[i--] = _candidates[--_curr_length];
573 573
          }
574 574
        }
575 575

	
576 576
        // Extend the list
577 577
        int cnt = _block_size;
578 578
        int limit = _head_length;
579 579

	
580 580
        for (e = _next_arc; e < _search_arc_num; ++e) {
581 581
          _cand_cost[e] = _state[e] *
582 582
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
583 583
          if (_cand_cost[e] < 0) {
584 584
            _candidates[_curr_length++] = e;
585 585
          }
586 586
          if (--cnt == 0) {
587 587
            if (_curr_length > limit) goto search_end;
588 588
            limit = 0;
589 589
            cnt = _block_size;
590 590
          }
591 591
        }
592 592
        for (e = 0; e < _next_arc; ++e) {
593 593
          _cand_cost[e] = _state[e] *
594 594
            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
595 595
          if (_cand_cost[e] < 0) {
596 596
            _candidates[_curr_length++] = e;
597 597
          }
598 598
          if (--cnt == 0) {
599 599
            if (_curr_length > limit) goto search_end;
600 600
            limit = 0;
601 601
            cnt = _block_size;
602 602
          }
603 603
        }
604 604
        if (_curr_length == 0) return false;
605 605
        
606 606
      search_end:
607 607

	
608 608
        // Make heap of the candidate list (approximating a partial sort)
609 609
        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
610 610
                   _sort_func );
611 611

	
612 612
        // Pop the first element of the heap
613 613
        _in_arc = _candidates[0];
614 614
        _next_arc = e;
615 615
        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
616 616
                  _sort_func );
617 617
        _curr_length = std::min(_head_length, _curr_length - 1);
618 618
        return true;
619 619
      }
620 620

	
621 621
    }; //class AlteringListPivotRule
622 622

	
623 623
  public:
624 624

	
625 625
    /// \brief Constructor.
626 626
    ///
627 627
    /// The constructor of the class.
628 628
    ///
629 629
    /// \param graph The digraph the algorithm runs on.
630 630
    /// \param arc_mixing Indicate if the arcs have to be stored in a
631 631
    /// mixed order in the internal data structure. 
632 632
    /// In special cases, it could lead to better overall performance,
633 633
    /// but it is usually slower. Therefore it is disabled by default.
634 634
    NetworkSimplex(const GR& graph, bool arc_mixing = false) :
635 635
      _graph(graph), _node_id(graph), _arc_id(graph),
636 636
      MAX(std::numeric_limits<Value>::max()),
637 637
      INF(std::numeric_limits<Value>::has_infinity ?
638 638
          std::numeric_limits<Value>::infinity() : MAX)
639 639
    {
640
      // Check the value types
640
      // Check the number types
641 641
      LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
642 642
        "The flow type of NetworkSimplex must be signed");
643 643
      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
644 644
        "The cost type of NetworkSimplex must be signed");
645 645
        
646 646
      // Resize vectors
647 647
      _node_num = countNodes(_graph);
648 648
      _arc_num = countArcs(_graph);
649 649
      int all_node_num = _node_num + 1;
650 650
      int max_arc_num = _arc_num + 2 * _node_num;
651 651

	
652 652
      _source.resize(max_arc_num);
653 653
      _target.resize(max_arc_num);
654 654

	
655 655
      _lower.resize(_arc_num);
656 656
      _upper.resize(_arc_num);
657 657
      _cap.resize(max_arc_num);
658 658
      _cost.resize(max_arc_num);
659 659
      _supply.resize(all_node_num);
660 660
      _flow.resize(max_arc_num);
661 661
      _pi.resize(all_node_num);
662 662

	
663 663
      _parent.resize(all_node_num);
664 664
      _pred.resize(all_node_num);
665 665
      _forward.resize(all_node_num);
666 666
      _thread.resize(all_node_num);
667 667
      _rev_thread.resize(all_node_num);
668 668
      _succ_num.resize(all_node_num);
669 669
      _last_succ.resize(all_node_num);
670 670
      _state.resize(max_arc_num);
671 671

	
672 672
      // Copy the graph
673 673
      int i = 0;
674 674
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
675 675
        _node_id[n] = i;
676 676
      }
677 677
      if (arc_mixing) {
678 678
        // Store the arcs in a mixed order
679 679
        int k = std::max(int(std::sqrt(double(_arc_num))), 10);
680 680
        int i = 0, j = 0;
681 681
        for (ArcIt a(_graph); a != INVALID; ++a) {
682 682
          _arc_id[a] = i;
683 683
          _source[i] = _node_id[_graph.source(a)];
684 684
          _target[i] = _node_id[_graph.target(a)];
685 685
          if ((i += k) >= _arc_num) i = ++j;
686 686
        }
687 687
      } else {
688 688
        // Store the arcs in the original order
689 689
        int i = 0;
690 690
        for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
691 691
          _arc_id[a] = i;
692 692
          _source[i] = _node_id[_graph.source(a)];
693 693
          _target[i] = _node_id[_graph.target(a)];
694 694
        }
695 695
      }
696 696
      
697 697
      // Reset parameters
698 698
      reset();
699 699
    }
700 700

	
701 701
    /// \name Parameters
702 702
    /// The parameters of the algorithm can be specified using these
703 703
    /// functions.
704 704

	
705 705
    /// @{
706 706

	
707 707
    /// \brief Set the lower bounds on the arcs.
708 708
    ///
709 709
    /// This function sets the lower bounds on the arcs.
710 710
    /// If it is not used before calling \ref run(), the lower bounds
711 711
    /// will be set to zero on all arcs.
712 712
    ///
713 713
    /// \param map An arc map storing the lower bounds.
714 714
    /// Its \c Value type must be convertible to the \c Value type
715 715
    /// of the algorithm.
716 716
    ///
717 717
    /// \return <tt>(*this)</tt>
718 718
    template <typename LowerMap>
719 719
    NetworkSimplex& lowerMap(const LowerMap& map) {
720 720
      _have_lower = true;
721 721
      for (ArcIt a(_graph); a != INVALID; ++a) {
722 722
        _lower[_arc_id[a]] = map[a];
723 723
      }
724 724
      return *this;
725 725
    }
726 726

	
727 727
    /// \brief Set the upper bounds (capacities) on the arcs.
728 728
    ///
729 729
    /// This function sets the upper bounds (capacities) on the arcs.
730 730
    /// If it is not used before calling \ref run(), the upper bounds
731 731
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
732
    /// unbounded from above on each arc).
732
    /// unbounded from above).
733 733
    ///
734 734
    /// \param map An arc map storing the upper bounds.
735 735
    /// Its \c Value type must be convertible to the \c Value type
736 736
    /// of the algorithm.
737 737
    ///
738 738
    /// \return <tt>(*this)</tt>
739 739
    template<typename UpperMap>
740 740
    NetworkSimplex& upperMap(const UpperMap& map) {
741 741
      for (ArcIt a(_graph); a != INVALID; ++a) {
742 742
        _upper[_arc_id[a]] = map[a];
743 743
      }
744 744
      return *this;
745 745
    }
746 746

	
747 747
    /// \brief Set the costs of the arcs.
748 748
    ///
749 749
    /// This function sets the costs of the arcs.
750 750
    /// If it is not used before calling \ref run(), the costs
751 751
    /// will be set to \c 1 on all arcs.
752 752
    ///
753 753
    /// \param map An arc map storing the costs.
754 754
    /// Its \c Value type must be convertible to the \c Cost type
755 755
    /// of the algorithm.
756 756
    ///
757 757
    /// \return <tt>(*this)</tt>
758 758
    template<typename CostMap>
759 759
    NetworkSimplex& costMap(const CostMap& map) {
760 760
      for (ArcIt a(_graph); a != INVALID; ++a) {
761 761
        _cost[_arc_id[a]] = map[a];
762 762
      }
763 763
      return *this;
764 764
    }
765 765

	
766 766
    /// \brief Set the supply values of the nodes.
767 767
    ///
768 768
    /// This function sets the supply values of the nodes.
769 769
    /// If neither this function nor \ref stSupply() is used before
770 770
    /// calling \ref run(), the supply of each node will be set to zero.
771 771
    ///
772 772
    /// \param map A node map storing the supply values.
773 773
    /// Its \c Value type must be convertible to the \c Value type
774 774
    /// of the algorithm.
775 775
    ///
776 776
    /// \return <tt>(*this)</tt>
777 777
    template<typename SupplyMap>
778 778
    NetworkSimplex& supplyMap(const SupplyMap& map) {
779 779
      for (NodeIt n(_graph); n != INVALID; ++n) {
780 780
        _supply[_node_id[n]] = map[n];
781 781
      }
782 782
      return *this;
783 783
    }
784 784

	
785 785
    /// \brief Set single source and target nodes and a supply value.
786 786
    ///
787 787
    /// This function sets a single source node and a single target node
788 788
    /// and the required flow value.
789 789
    /// If neither this function nor \ref supplyMap() is used before
790 790
    /// calling \ref run(), the supply of each node will be set to zero.
791 791
    ///
792 792
    /// Using this function has the same effect as using \ref supplyMap()
793 793
    /// with such a map in which \c k is assigned to \c s, \c -k is
794 794
    /// assigned to \c t and all other nodes have zero supply value.
795 795
    ///
796 796
    /// \param s The source node.
797 797
    /// \param t The target node.
798 798
    /// \param k The required amount of flow from node \c s to node \c t
799 799
    /// (i.e. the supply of \c s and the demand of \c t).
800 800
    ///
801 801
    /// \return <tt>(*this)</tt>
802 802
    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
803 803
      for (int i = 0; i != _node_num; ++i) {
804 804
        _supply[i] = 0;
805 805
      }
806 806
      _supply[_node_id[s]] =  k;
807 807
      _supply[_node_id[t]] = -k;
808 808
      return *this;
809 809
    }
810 810
    
811 811
    /// \brief Set the type of the supply constraints.
812 812
    ///
813 813
    /// This function sets the type of the supply/demand constraints.
814 814
    /// If it is not used before calling \ref run(), the \ref GEQ supply
815 815
    /// type will be used.
816 816
    ///
817 817
    /// For more information, see \ref SupplyType.
818 818
    ///
819 819
    /// \return <tt>(*this)</tt>
820 820
    NetworkSimplex& supplyType(SupplyType supply_type) {
821 821
      _stype = supply_type;
822 822
      return *this;
823 823
    }
824 824

	
825 825
    /// @}
826 826

	
827 827
    /// \name Execution Control
828 828
    /// The algorithm can be executed using \ref run().
829 829

	
830 830
    /// @{
831 831

	
832 832
    /// \brief Run the algorithm.
833 833
    ///
834 834
    /// This function runs the algorithm.
835 835
    /// The paramters can be specified using functions \ref lowerMap(),
836 836
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 
837 837
    /// \ref supplyType().
838 838
    /// For example,
839 839
    /// \code
840 840
    ///   NetworkSimplex<ListDigraph> ns(graph);
841 841
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
842 842
    ///     .supplyMap(sup).run();
843 843
    /// \endcode
844 844
    ///
845 845
    /// This function can be called more than once. All the parameters
846 846
    /// that have been given are kept for the next call, unless
847 847
    /// \ref reset() is called, thus only the modified parameters
848 848
    /// have to be set again. See \ref reset() for examples.
849 849
    /// However, the underlying digraph must not be modified after this
850 850
    /// class have been constructed, since it copies and extends the graph.
851 851
    ///
852 852
    /// \param pivot_rule The pivot rule that will be used during the
853 853
    /// algorithm. For more information, see \ref PivotRule.
854 854
    ///
855 855
    /// \return \c INFEASIBLE if no feasible flow exists,
856 856
    /// \n \c OPTIMAL if the problem has optimal solution
857 857
    /// (i.e. it is feasible and bounded), and the algorithm has found
858 858
    /// optimal flow and node potentials (primal and dual solutions),
859 859
    /// \n \c UNBOUNDED if the objective function of the problem is
860 860
    /// unbounded, i.e. there is a directed cycle having negative total
861 861
    /// cost and infinite upper bound.
862 862
    ///
863 863
    /// \see ProblemType, PivotRule
864 864
    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
865 865
      if (!init()) return INFEASIBLE;
866 866
      return start(pivot_rule);
867 867
    }
868 868

	
869 869
    /// \brief Reset all the parameters that have been given before.
870 870
    ///
871 871
    /// This function resets all the paramaters that have been given
872 872
    /// before using functions \ref lowerMap(), \ref upperMap(),
873 873
    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
874 874
    ///
875 875
    /// It is useful for multiple run() calls. If this function is not
876 876
    /// used, all the parameters given before are kept for the next
877 877
    /// \ref run() call.
878 878
    /// However, the underlying digraph must not be modified after this
879 879
    /// class have been constructed, since it copies and extends the graph.
880 880
    ///
881 881
    /// For example,
882 882
    /// \code
883 883
    ///   NetworkSimplex<ListDigraph> ns(graph);
884 884
    ///
885 885
    ///   // First run
886 886
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
887 887
    ///     .supplyMap(sup).run();
888 888
    ///
889 889
    ///   // Run again with modified cost map (reset() is not called,
890 890
    ///   // so only the cost map have to be set again)
891 891
    ///   cost[e] += 100;
892 892
    ///   ns.costMap(cost).run();
893 893
    ///
894 894
    ///   // Run again from scratch using reset()
895 895
    ///   // (the lower bounds will be set to zero on all arcs)
896 896
    ///   ns.reset();
897 897
    ///   ns.upperMap(capacity).costMap(cost)
898 898
    ///     .supplyMap(sup).run();
899 899
    /// \endcode
900 900
    ///
901 901
    /// \return <tt>(*this)</tt>
902 902
    NetworkSimplex& reset() {
903 903
      for (int i = 0; i != _node_num; ++i) {
904 904
        _supply[i] = 0;
905 905
      }
906 906
      for (int i = 0; i != _arc_num; ++i) {
907 907
        _lower[i] = 0;
908 908
        _upper[i] = INF;
909 909
        _cost[i] = 1;
910 910
      }
911 911
      _have_lower = false;
912 912
      _stype = GEQ;
913 913
      return *this;
914 914
    }
915 915

	
916 916
    /// @}
917 917

	
918 918
    /// \name Query Functions
919 919
    /// The results of the algorithm can be obtained using these
920 920
    /// functions.\n
921 921
    /// The \ref run() function must be called before using them.
922 922

	
923 923
    /// @{
924 924

	
0 comments (0 inline)