gravatar
kpeter (Peter Kovacs)
kpeter@inf.elte.hu
Handle graph changes in the MCF algorithms (#327) The reset() functions are renamed to resetParams() and the new reset() functions handle the graph chnages, as well.
0 5 0
default
5 files changed with 423 insertions and 311 deletions:
↑ Collapse diff ↑
Show white space 384 line context
... ...
@@ -125,656 +125,685 @@
125 125
      /// upper bound. It means that the objective function is unbounded
126 126
      /// on that arc, however, note that it could actually be bounded
127 127
      /// over the feasible flows, but this algroithm cannot handle
128 128
      /// these cases.
129 129
      UNBOUNDED
130 130
    };
131 131
  
132 132
  private:
133 133

	
134 134
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
135 135

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

	
141 141
  private:
142 142

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

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

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

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

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

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

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

	
189 189
  private:
190 190

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

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

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

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

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

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

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

	
267 267
        return t;
268 268
      }
269 269

	
270 270
    }; //class ResidualDijkstra
271 271

	
272 272
  public:
273 273

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

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

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

	
296 296
    /// @}
297 297

	
298 298
  public:
299 299

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

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

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

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

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

	
383 321
    /// \name Parameters
384 322
    /// The parameters of the algorithm can be specified using these
385 323
    /// functions.
386 324

	
387 325
    /// @{
388 326

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

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

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

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

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

	
497 435
    /// \name Execution control
498 436
    /// The algorithm can be executed using \ref run().
499 437

	
500 438
    /// @{
501 439

	
502 440
    /// \brief Run the algorithm.
503 441
    ///
504 442
    /// This function runs the algorithm.
505 443
    /// The paramters can be specified using functions \ref lowerMap(),
506 444
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
507 445
    /// For example,
508 446
    /// \code
509 447
    ///   CapacityScaling<ListDigraph> cs(graph);
510 448
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
511 449
    ///     .supplyMap(sup).run();
512 450
    /// \endcode
513 451
    ///
514
    /// This function can be called more than once. All the parameters
515
    /// that have been given are kept for the next call, unless
516
    /// \ref reset() is called, thus only the modified parameters
517
    /// have to be set again. See \ref reset() for examples.
518
    /// However, the underlying digraph must not be modified after this
519
    /// class have been constructed, since it copies and extends the graph.
452
    /// This function can be called more than once. All the given parameters
453
    /// are kept for the next call, unless \ref resetParams() or \ref reset()
454
    /// is used, thus only the modified parameters have to be set again.
455
    /// If the underlying digraph was also modified after the construction
456
    /// of the class (or the last \ref reset() call), then the \ref reset()
457
    /// function must be called.
520 458
    ///
521 459
    /// \param factor The capacity scaling factor. It must be larger than
522 460
    /// one to use scaling. If it is less or equal to one, then scaling
523 461
    /// will be disabled.
524 462
    ///
525 463
    /// \return \c INFEASIBLE if no feasible flow exists,
526 464
    /// \n \c OPTIMAL if the problem has optimal solution
527 465
    /// (i.e. it is feasible and bounded), and the algorithm has found
528 466
    /// optimal flow and node potentials (primal and dual solutions),
529 467
    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
530 468
    /// and infinite upper bound. It means that the objective function
531 469
    /// is unbounded on that arc, however, note that it could actually be
532 470
    /// bounded over the feasible flows, but this algroithm cannot handle
533 471
    /// these cases.
534 472
    ///
535 473
    /// \see ProblemType
474
    /// \see resetParams(), reset()
536 475
    ProblemType run(int factor = 4) {
537 476
      _factor = factor;
538 477
      ProblemType pt = init();
539 478
      if (pt != OPTIMAL) return pt;
540 479
      return start();
541 480
    }
542 481

	
543 482
    /// \brief Reset all the parameters that have been given before.
544 483
    ///
545 484
    /// This function resets all the paramaters that have been given
546 485
    /// before using functions \ref lowerMap(), \ref upperMap(),
547 486
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
548 487
    ///
549
    /// It is useful for multiple run() calls. If this function is not
550
    /// used, all the parameters given before are kept for the next
551
    /// \ref run() call.
552
    /// However, the underlying digraph must not be modified after this
553
    /// class have been constructed, since it copies and extends the graph.
488
    /// It is useful for multiple \ref run() calls. Basically, all the given
489
    /// parameters are kept for the next \ref run() call, unless
490
    /// \ref resetParams() or \ref reset() is used.
491
    /// If the underlying digraph was also modified after the construction
492
    /// of the class or the last \ref reset() call, then the \ref reset()
493
    /// function must be used, otherwise \ref resetParams() is sufficient.
554 494
    ///
555 495
    /// For example,
556 496
    /// \code
557 497
    ///   CapacityScaling<ListDigraph> cs(graph);
558 498
    ///
559 499
    ///   // First run
560 500
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
561 501
    ///     .supplyMap(sup).run();
562 502
    ///
563
    ///   // Run again with modified cost map (reset() is not called,
503
    ///   // Run again with modified cost map (resetParams() is not called,
564 504
    ///   // so only the cost map have to be set again)
565 505
    ///   cost[e] += 100;
566 506
    ///   cs.costMap(cost).run();
567 507
    ///
568
    ///   // Run again from scratch using reset()
508
    ///   // Run again from scratch using resetParams()
569 509
    ///   // (the lower bounds will be set to zero on all arcs)
570
    ///   cs.reset();
510
    ///   cs.resetParams();
571 511
    ///   cs.upperMap(capacity).costMap(cost)
572 512
    ///     .supplyMap(sup).run();
573 513
    /// \endcode
574 514
    ///
575 515
    /// \return <tt>(*this)</tt>
576
    CapacityScaling& reset() {
516
    ///
517
    /// \see reset(), run()
518
    CapacityScaling& resetParams() {
577 519
      for (int i = 0; i != _node_num; ++i) {
578 520
        _supply[i] = 0;
579 521
      }
580 522
      for (int j = 0; j != _res_arc_num; ++j) {
581 523
        _lower[j] = 0;
582 524
        _upper[j] = INF;
583 525
        _cost[j] = _forward[j] ? 1 : -1;
584 526
      }
585 527
      _have_lower = false;
586 528
      return *this;
587 529
    }
588 530

	
531
    /// \brief Reset the internal data structures and all the parameters
532
    /// that have been given before.
533
    ///
534
    /// This function resets the internal data structures and all the
535
    /// paramaters that have been given before using functions \ref lowerMap(),
536
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
537
    ///
538
    /// It is useful for multiple \ref run() calls. Basically, all the given
539
    /// parameters are kept for the next \ref run() call, unless
540
    /// \ref resetParams() or \ref reset() is used.
541
    /// If the underlying digraph was also modified after the construction
542
    /// of the class or the last \ref reset() call, then the \ref reset()
543
    /// function must be used, otherwise \ref resetParams() is sufficient.
544
    ///
545
    /// See \ref resetParams() for examples.
546
    ///
547
    /// \return <tt>(*this)</tt>
548
    ///
549
    /// \see resetParams(), run()
550
    CapacityScaling& reset() {
551
      // Resize vectors
552
      _node_num = countNodes(_graph);
553
      _arc_num = countArcs(_graph);
554
      _res_arc_num = 2 * (_arc_num + _node_num);
555
      _root = _node_num;
556
      ++_node_num;
557

	
558
      _first_out.resize(_node_num + 1);
559
      _forward.resize(_res_arc_num);
560
      _source.resize(_res_arc_num);
561
      _target.resize(_res_arc_num);
562
      _reverse.resize(_res_arc_num);
563

	
564
      _lower.resize(_res_arc_num);
565
      _upper.resize(_res_arc_num);
566
      _cost.resize(_res_arc_num);
567
      _supply.resize(_node_num);
568
      
569
      _res_cap.resize(_res_arc_num);
570
      _pi.resize(_node_num);
571
      _excess.resize(_node_num);
572
      _pred.resize(_node_num);
573

	
574
      // Copy the graph
575
      int i = 0, j = 0, k = 2 * _arc_num + _node_num - 1;
576
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
577
        _node_id[n] = i;
578
      }
579
      i = 0;
580
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
581
        _first_out[i] = j;
582
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
583
          _arc_idf[a] = j;
584
          _forward[j] = true;
585
          _source[j] = i;
586
          _target[j] = _node_id[_graph.runningNode(a)];
587
        }
588
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
589
          _arc_idb[a] = j;
590
          _forward[j] = false;
591
          _source[j] = i;
592
          _target[j] = _node_id[_graph.runningNode(a)];
593
        }
594
        _forward[j] = false;
595
        _source[j] = i;
596
        _target[j] = _root;
597
        _reverse[j] = k;
598
        _forward[k] = true;
599
        _source[k] = _root;
600
        _target[k] = i;
601
        _reverse[k] = j;
602
        ++j; ++k;
603
      }
604
      _first_out[i] = j;
605
      _first_out[_node_num] = k;
606
      for (ArcIt a(_graph); a != INVALID; ++a) {
607
        int fi = _arc_idf[a];
608
        int bi = _arc_idb[a];
609
        _reverse[fi] = bi;
610
        _reverse[bi] = fi;
611
      }
612
      
613
      // Reset parameters
614
      resetParams();
615
      return *this;
616
    }
617

	
589 618
    /// @}
590 619

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

	
596 625
    /// @{
597 626

	
598 627
    /// \brief Return the total cost of the found flow.
599 628
    ///
600 629
    /// This function returns the total cost of the found flow.
601 630
    /// Its complexity is O(e).
602 631
    ///
603 632
    /// \note The return type of the function can be specified as a
604 633
    /// template parameter. For example,
605 634
    /// \code
606 635
    ///   cs.totalCost<double>();
607 636
    /// \endcode
608 637
    /// It is useful if the total cost cannot be stored in the \c Cost
609 638
    /// type of the algorithm, which is the default return type of the
610 639
    /// function.
611 640
    ///
612 641
    /// \pre \ref run() must be called before using this function.
613 642
    template <typename Number>
614 643
    Number totalCost() const {
615 644
      Number c = 0;
616 645
      for (ArcIt a(_graph); a != INVALID; ++a) {
617 646
        int i = _arc_idb[a];
618 647
        c += static_cast<Number>(_res_cap[i]) *
619 648
             (-static_cast<Number>(_cost[i]));
620 649
      }
621 650
      return c;
622 651
    }
623 652

	
624 653
#ifndef DOXYGEN
625 654
    Cost totalCost() const {
626 655
      return totalCost<Cost>();
627 656
    }
628 657
#endif
629 658

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

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

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

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

	
678 707
    /// @}
679 708

	
680 709
  private:
681 710

	
682 711
    // Initialize the algorithm
683 712
    ProblemType init() {
684 713
      if (_node_num <= 1) return INFEASIBLE;
685 714

	
686 715
      // Check the sum of supply values
687 716
      _sum_supply = 0;
688 717
      for (int i = 0; i != _root; ++i) {
689 718
        _sum_supply += _supply[i];
690 719
      }
691 720
      if (_sum_supply > 0) return INFEASIBLE;
692 721
      
693 722
      // Initialize vectors
694 723
      for (int i = 0; i != _root; ++i) {
695 724
        _pi[i] = 0;
696 725
        _excess[i] = _supply[i];
697 726
      }
698 727

	
699 728
      // Remove non-zero lower bounds
700 729
      const Value MAX = std::numeric_limits<Value>::max();
701 730
      int last_out;
702 731
      if (_have_lower) {
703 732
        for (int i = 0; i != _root; ++i) {
704 733
          last_out = _first_out[i+1];
705 734
          for (int j = _first_out[i]; j != last_out; ++j) {
706 735
            if (_forward[j]) {
707 736
              Value c = _lower[j];
708 737
              if (c >= 0) {
709 738
                _res_cap[j] = _upper[j] < MAX ? _upper[j] - c : INF;
710 739
              } else {
711 740
                _res_cap[j] = _upper[j] < MAX + c ? _upper[j] - c : INF;
712 741
              }
713 742
              _excess[i] -= c;
714 743
              _excess[_target[j]] += c;
715 744
            } else {
716 745
              _res_cap[j] = 0;
717 746
            }
718 747
          }
719 748
        }
720 749
      } else {
721 750
        for (int j = 0; j != _res_arc_num; ++j) {
722 751
          _res_cap[j] = _forward[j] ? _upper[j] : 0;
723 752
        }
724 753
      }
725 754

	
726 755
      // Handle negative costs
727 756
      for (int i = 0; i != _root; ++i) {
728 757
        last_out = _first_out[i+1] - 1;
729 758
        for (int j = _first_out[i]; j != last_out; ++j) {
730 759
          Value rc = _res_cap[j];
731 760
          if (_cost[j] < 0 && rc > 0) {
732 761
            if (rc >= MAX) return UNBOUNDED;
733 762
            _excess[i] -= rc;
734 763
            _excess[_target[j]] += rc;
735 764
            _res_cap[j] = 0;
736 765
            _res_cap[_reverse[j]] += rc;
737 766
          }
738 767
        }
739 768
      }
740 769
      
741 770
      // Handle GEQ supply type
742 771
      if (_sum_supply < 0) {
743 772
        _pi[_root] = 0;
744 773
        _excess[_root] = -_sum_supply;
745 774
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
746 775
          int ra = _reverse[a];
747 776
          _res_cap[a] = -_sum_supply + 1;
748 777
          _res_cap[ra] = 0;
749 778
          _cost[a] = 0;
750 779
          _cost[ra] = 0;
751 780
        }
752 781
      } else {
753 782
        _pi[_root] = 0;
754 783
        _excess[_root] = 0;
755 784
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
756 785
          int ra = _reverse[a];
757 786
          _res_cap[a] = 1;
758 787
          _res_cap[ra] = 0;
759 788
          _cost[a] = 0;
760 789
          _cost[ra] = 0;
761 790
        }
762 791
      }
763 792

	
764 793
      // Initialize delta value
765 794
      if (_factor > 1) {
766 795
        // With scaling
767 796
        Value max_sup = 0, max_dem = 0;
768 797
        for (int i = 0; i != _node_num; ++i) {
769 798
          Value ex = _excess[i];
770 799
          if ( ex > max_sup) max_sup =  ex;
771 800
          if (-ex > max_dem) max_dem = -ex;
772 801
        }
773 802
        Value max_cap = 0;
774 803
        for (int j = 0; j != _res_arc_num; ++j) {
775 804
          if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
776 805
        }
777 806
        max_sup = std::min(std::min(max_sup, max_dem), max_cap);
778 807
        for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
779 808
      } else {
780 809
        // Without scaling
Show white space 384 line context
... ...
@@ -144,668 +144,690 @@
144 144
    /// The \ref CostScalingDefaultTraits "traits class" of the algorithm
145 145
    typedef TR Traits;
146 146

	
147 147
  public:
148 148

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

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

	
195 195
  private:
196 196

	
197 197
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
198 198

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

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

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

	
227 227
    private:
228 228
      std::vector<Value>& _v;
229 229
    };
230 230

	
231 231
    typedef StaticVectorMap<StaticDigraph::Node, LargeCost> LargeCostNodeMap;
232 232
    typedef StaticVectorMap<StaticDigraph::Arc, LargeCost> LargeCostArcMap;
233 233

	
234 234
  private:
235 235

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

	
244 244
    // Parameters of the problem
245 245
    bool _have_lower;
246 246
    Value _sum_supply;
247 247

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

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

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

	
271 271
    // Data for scaling
272 272
    LargeCost _epsilon;
273 273
    int _alpha;
274 274

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

	
292 292
  public:
293 293

	
294 294
    /// \name Named Template Parameters
295 295
    /// @{
296 296

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

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

	
314 314
    /// @}
315 315

	
316 316
  public:
317 317

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

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

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

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

	
360
      _arc_vec.reserve(_res_arc_num);
361
      _cost_vec.reserve(_res_arc_num);
362

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

	
406 340
    /// \name Parameters
407 341
    /// The parameters of the algorithm can be specified using these
408 342
    /// functions.
409 343

	
410 344
    /// @{
411 345

	
412 346
    /// \brief Set the lower bounds on the arcs.
413 347
    ///
414 348
    /// This function sets the lower bounds on the arcs.
415 349
    /// If it is not used before calling \ref run(), the lower bounds
416 350
    /// will be set to zero on all arcs.
417 351
    ///
418 352
    /// \param map An arc map storing the lower bounds.
419 353
    /// Its \c Value type must be convertible to the \c Value type
420 354
    /// of the algorithm.
421 355
    ///
422 356
    /// \return <tt>(*this)</tt>
423 357
    template <typename LowerMap>
424 358
    CostScaling& lowerMap(const LowerMap& map) {
425 359
      _have_lower = true;
426 360
      for (ArcIt a(_graph); a != INVALID; ++a) {
427 361
        _lower[_arc_idf[a]] = map[a];
428 362
        _lower[_arc_idb[a]] = map[a];
429 363
      }
430 364
      return *this;
431 365
    }
432 366

	
433 367
    /// \brief Set the upper bounds (capacities) on the arcs.
434 368
    ///
435 369
    /// This function sets the upper bounds (capacities) on the arcs.
436 370
    /// If it is not used before calling \ref run(), the upper bounds
437 371
    /// will be set to \ref INF on all arcs (i.e. the flow value will be
438 372
    /// unbounded from above).
439 373
    ///
440 374
    /// \param map An arc map storing the upper bounds.
441 375
    /// Its \c Value type must be convertible to the \c Value type
442 376
    /// of the algorithm.
443 377
    ///
444 378
    /// \return <tt>(*this)</tt>
445 379
    template<typename UpperMap>
446 380
    CostScaling& upperMap(const UpperMap& map) {
447 381
      for (ArcIt a(_graph); a != INVALID; ++a) {
448 382
        _upper[_arc_idf[a]] = map[a];
449 383
      }
450 384
      return *this;
451 385
    }
452 386

	
453 387
    /// \brief Set the costs of the arcs.
454 388
    ///
455 389
    /// This function sets the costs of the arcs.
456 390
    /// If it is not used before calling \ref run(), the costs
457 391
    /// will be set to \c 1 on all arcs.
458 392
    ///
459 393
    /// \param map An arc map storing the costs.
460 394
    /// Its \c Value type must be convertible to the \c Cost type
461 395
    /// of the algorithm.
462 396
    ///
463 397
    /// \return <tt>(*this)</tt>
464 398
    template<typename CostMap>
465 399
    CostScaling& costMap(const CostMap& map) {
466 400
      for (ArcIt a(_graph); a != INVALID; ++a) {
467 401
        _scost[_arc_idf[a]] =  map[a];
468 402
        _scost[_arc_idb[a]] = -map[a];
469 403
      }
470 404
      return *this;
471 405
    }
472 406

	
473 407
    /// \brief Set the supply values of the nodes.
474 408
    ///
475 409
    /// This function sets the supply values of the nodes.
476 410
    /// If neither this function nor \ref stSupply() is used before
477 411
    /// calling \ref run(), the supply of each node will be set to zero.
478 412
    ///
479 413
    /// \param map A node map storing the supply values.
480 414
    /// Its \c Value type must be convertible to the \c Value type
481 415
    /// of the algorithm.
482 416
    ///
483 417
    /// \return <tt>(*this)</tt>
484 418
    template<typename SupplyMap>
485 419
    CostScaling& supplyMap(const SupplyMap& map) {
486 420
      for (NodeIt n(_graph); n != INVALID; ++n) {
487 421
        _supply[_node_id[n]] = map[n];
488 422
      }
489 423
      return *this;
490 424
    }
491 425

	
492 426
    /// \brief Set single source and target nodes and a supply value.
493 427
    ///
494 428
    /// This function sets a single source node and a single target node
495 429
    /// and the required flow value.
496 430
    /// If neither this function nor \ref supplyMap() is used before
497 431
    /// calling \ref run(), the supply of each node will be set to zero.
498 432
    ///
499 433
    /// Using this function has the same effect as using \ref supplyMap()
500 434
    /// with such a map in which \c k is assigned to \c s, \c -k is
501 435
    /// assigned to \c t and all other nodes have zero supply value.
502 436
    ///
503 437
    /// \param s The source node.
504 438
    /// \param t The target node.
505 439
    /// \param k The required amount of flow from node \c s to node \c t
506 440
    /// (i.e. the supply of \c s and the demand of \c t).
507 441
    ///
508 442
    /// \return <tt>(*this)</tt>
509 443
    CostScaling& stSupply(const Node& s, const Node& t, Value k) {
510 444
      for (int i = 0; i != _res_node_num; ++i) {
511 445
        _supply[i] = 0;
512 446
      }
513 447
      _supply[_node_id[s]] =  k;
514 448
      _supply[_node_id[t]] = -k;
515 449
      return *this;
516 450
    }
517 451
    
518 452
    /// @}
519 453

	
520 454
    /// \name Execution control
521 455
    /// The algorithm can be executed using \ref run().
522 456

	
523 457
    /// @{
524 458

	
525 459
    /// \brief Run the algorithm.
526 460
    ///
527 461
    /// This function runs the algorithm.
528 462
    /// The paramters can be specified using functions \ref lowerMap(),
529 463
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
530 464
    /// For example,
531 465
    /// \code
532 466
    ///   CostScaling<ListDigraph> cs(graph);
533 467
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
534 468
    ///     .supplyMap(sup).run();
535 469
    /// \endcode
536 470
    ///
537
    /// This function can be called more than once. All the parameters
538
    /// that have been given are kept for the next call, unless
539
    /// \ref reset() is called, thus only the modified parameters
540
    /// have to be set again. See \ref reset() for examples.
541
    /// However, the underlying digraph must not be modified after this
542
    /// class have been constructed, since it copies and extends the graph.
471
    /// This function can be called more than once. All the given parameters
472
    /// are kept for the next call, unless \ref resetParams() or \ref reset()
473
    /// is used, thus only the modified parameters have to be set again.
474
    /// If the underlying digraph was also modified after the construction
475
    /// of the class (or the last \ref reset() call), then the \ref reset()
476
    /// function must be called.
543 477
    ///
544 478
    /// \param method The internal method that will be used in the
545 479
    /// algorithm. For more information, see \ref Method.
546 480
    /// \param factor The cost scaling factor. It must be larger than one.
547 481
    ///
548 482
    /// \return \c INFEASIBLE if no feasible flow exists,
549 483
    /// \n \c OPTIMAL if the problem has optimal solution
550 484
    /// (i.e. it is feasible and bounded), and the algorithm has found
551 485
    /// optimal flow and node potentials (primal and dual solutions),
552 486
    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
553 487
    /// and infinite upper bound. It means that the objective function
554 488
    /// is unbounded on that arc, however, note that it could actually be
555 489
    /// bounded over the feasible flows, but this algroithm cannot handle
556 490
    /// these cases.
557 491
    ///
558 492
    /// \see ProblemType, Method
493
    /// \see resetParams(), reset()
559 494
    ProblemType run(Method method = PARTIAL_AUGMENT, int factor = 8) {
560 495
      _alpha = factor;
561 496
      ProblemType pt = init();
562 497
      if (pt != OPTIMAL) return pt;
563 498
      start(method);
564 499
      return OPTIMAL;
565 500
    }
566 501

	
567 502
    /// \brief Reset all the parameters that have been given before.
568 503
    ///
569 504
    /// This function resets all the paramaters that have been given
570 505
    /// before using functions \ref lowerMap(), \ref upperMap(),
571 506
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
572 507
    ///
573
    /// It is useful for multiple run() calls. If this function is not
574
    /// used, all the parameters given before are kept for the next
575
    /// \ref run() call.
576
    /// However, the underlying digraph must not be modified after this
577
    /// class have been constructed, since it copies and extends the graph.
508
    /// It is useful for multiple \ref run() calls. Basically, all the given
509
    /// parameters are kept for the next \ref run() call, unless
510
    /// \ref resetParams() or \ref reset() is used.
511
    /// If the underlying digraph was also modified after the construction
512
    /// of the class or the last \ref reset() call, then the \ref reset()
513
    /// function must be used, otherwise \ref resetParams() is sufficient.
578 514
    ///
579 515
    /// For example,
580 516
    /// \code
581 517
    ///   CostScaling<ListDigraph> cs(graph);
582 518
    ///
583 519
    ///   // First run
584 520
    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)
585 521
    ///     .supplyMap(sup).run();
586 522
    ///
587
    ///   // Run again with modified cost map (reset() is not called,
523
    ///   // Run again with modified cost map (resetParams() is not called,
588 524
    ///   // so only the cost map have to be set again)
589 525
    ///   cost[e] += 100;
590 526
    ///   cs.costMap(cost).run();
591 527
    ///
592
    ///   // Run again from scratch using reset()
528
    ///   // Run again from scratch using resetParams()
593 529
    ///   // (the lower bounds will be set to zero on all arcs)
594
    ///   cs.reset();
530
    ///   cs.resetParams();
595 531
    ///   cs.upperMap(capacity).costMap(cost)
596 532
    ///     .supplyMap(sup).run();
597 533
    /// \endcode
598 534
    ///
599 535
    /// \return <tt>(*this)</tt>
600
    CostScaling& reset() {
536
    ///
537
    /// \see reset(), run()
538
    CostScaling& resetParams() {
601 539
      for (int i = 0; i != _res_node_num; ++i) {
602 540
        _supply[i] = 0;
603 541
      }
604 542
      int limit = _first_out[_root];
605 543
      for (int j = 0; j != limit; ++j) {
606 544
        _lower[j] = 0;
607 545
        _upper[j] = INF;
608 546
        _scost[j] = _forward[j] ? 1 : -1;
609 547
      }
610 548
      for (int j = limit; j != _res_arc_num; ++j) {
611 549
        _lower[j] = 0;
612 550
        _upper[j] = INF;
613 551
        _scost[j] = 0;
614 552
        _scost[_reverse[j]] = 0;
615 553
      }      
616 554
      _have_lower = false;
617 555
      return *this;
618 556
    }
619 557

	
558
    /// \brief Reset all the parameters that have been given before.
559
    ///
560
    /// This function resets all the paramaters that have been given
561
    /// before using functions \ref lowerMap(), \ref upperMap(),
562
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
563
    ///
564
    /// It is useful for multiple run() calls. If this function is not
565
    /// used, all the parameters given before are kept for the next
566
    /// \ref run() call.
567
    /// However, the underlying digraph must not be modified after this
568
    /// class have been constructed, since it copies and extends the graph.
569
    /// \return <tt>(*this)</tt>
570
    CostScaling& reset() {
571
      // Resize vectors
572
      _node_num = countNodes(_graph);
573
      _arc_num = countArcs(_graph);
574
      _res_node_num = _node_num + 1;
575
      _res_arc_num = 2 * (_arc_num + _node_num);
576
      _root = _node_num;
577

	
578
      _first_out.resize(_res_node_num + 1);
579
      _forward.resize(_res_arc_num);
580
      _source.resize(_res_arc_num);
581
      _target.resize(_res_arc_num);
582
      _reverse.resize(_res_arc_num);
583

	
584
      _lower.resize(_res_arc_num);
585
      _upper.resize(_res_arc_num);
586
      _scost.resize(_res_arc_num);
587
      _supply.resize(_res_node_num);
588
      
589
      _res_cap.resize(_res_arc_num);
590
      _cost.resize(_res_arc_num);
591
      _pi.resize(_res_node_num);
592
      _excess.resize(_res_node_num);
593
      _next_out.resize(_res_node_num);
594

	
595
      _arc_vec.reserve(_res_arc_num);
596
      _cost_vec.reserve(_res_arc_num);
597

	
598
      // Copy the graph
599
      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
600
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
601
        _node_id[n] = i;
602
      }
603
      i = 0;
604
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
605
        _first_out[i] = j;
606
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
607
          _arc_idf[a] = j;
608
          _forward[j] = true;
609
          _source[j] = i;
610
          _target[j] = _node_id[_graph.runningNode(a)];
611
        }
612
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
613
          _arc_idb[a] = j;
614
          _forward[j] = false;
615
          _source[j] = i;
616
          _target[j] = _node_id[_graph.runningNode(a)];
617
        }
618
        _forward[j] = false;
619
        _source[j] = i;
620
        _target[j] = _root;
621
        _reverse[j] = k;
622
        _forward[k] = true;
623
        _source[k] = _root;
624
        _target[k] = i;
625
        _reverse[k] = j;
626
        ++j; ++k;
627
      }
628
      _first_out[i] = j;
629
      _first_out[_res_node_num] = k;
630
      for (ArcIt a(_graph); a != INVALID; ++a) {
631
        int fi = _arc_idf[a];
632
        int bi = _arc_idb[a];
633
        _reverse[fi] = bi;
634
        _reverse[bi] = fi;
635
      }
636
      
637
      // Reset parameters
638
      resetParams();
639
      return *this;
640
    }
641

	
620 642
    /// @}
621 643

	
622 644
    /// \name Query Functions
623 645
    /// The results of the algorithm can be obtained using these
624 646
    /// functions.\n
625 647
    /// The \ref run() function must be called before using them.
626 648

	
627 649
    /// @{
628 650

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

	
655 677
#ifndef DOXYGEN
656 678
    Cost totalCost() const {
657 679
      return totalCost<Cost>();
658 680
    }
659 681
#endif
660 682

	
661 683
    /// \brief Return the flow on the given arc.
662 684
    ///
663 685
    /// This function returns the flow on the given arc.
664 686
    ///
665 687
    /// \pre \ref run() must be called before using this function.
666 688
    Value flow(const Arc& a) const {
667 689
      return _res_cap[_arc_idb[a]];
668 690
    }
669 691

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

	
684 706
    /// \brief Return the potential (dual value) of the given node.
685 707
    ///
686 708
    /// This function returns the potential (dual value) of the
687 709
    /// given node.
688 710
    ///
689 711
    /// \pre \ref run() must be called before using this function.
690 712
    Cost potential(const Node& n) const {
691 713
      return static_cast<Cost>(_pi[_node_id[n]]);
692 714
    }
693 715

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

	
709 731
    /// @}
710 732

	
711 733
  private:
712 734

	
713 735
    // Initialize the algorithm
714 736
    ProblemType init() {
715 737
      if (_res_node_num <= 1) return INFEASIBLE;
716 738

	
717 739
      // Check the sum of supply values
718 740
      _sum_supply = 0;
719 741
      for (int i = 0; i != _root; ++i) {
720 742
        _sum_supply += _supply[i];
721 743
      }
722 744
      if (_sum_supply > 0) return INFEASIBLE;
723 745
      
724 746

	
725 747
      // Initialize vectors
726 748
      for (int i = 0; i != _res_node_num; ++i) {
727 749
        _pi[i] = 0;
728 750
        _excess[i] = _supply[i];
729 751
      }
730 752
      
731 753
      // Remove infinite upper bounds and check negative arcs
732 754
      const Value MAX = std::numeric_limits<Value>::max();
733 755
      int last_out;
734 756
      if (_have_lower) {
735 757
        for (int i = 0; i != _root; ++i) {
736 758
          last_out = _first_out[i+1];
737 759
          for (int j = _first_out[i]; j != last_out; ++j) {
738 760
            if (_forward[j]) {
739 761
              Value c = _scost[j] < 0 ? _upper[j] : _lower[j];
740 762
              if (c >= MAX) return UNBOUNDED;
741 763
              _excess[i] -= c;
742 764
              _excess[_target[j]] += c;
743 765
            }
744 766
          }
745 767
        }
746 768
      } else {
747 769
        for (int i = 0; i != _root; ++i) {
748 770
          last_out = _first_out[i+1];
749 771
          for (int j = _first_out[i]; j != last_out; ++j) {
750 772
            if (_forward[j] && _scost[j] < 0) {
751 773
              Value c = _upper[j];
752 774
              if (c >= MAX) return UNBOUNDED;
753 775
              _excess[i] -= c;
754 776
              _excess[_target[j]] += c;
755 777
            }
756 778
          }
757 779
        }
758 780
      }
759 781
      Value ex, max_cap = 0;
760 782
      for (int i = 0; i != _res_node_num; ++i) {
761 783
        ex = _excess[i];
762 784
        _excess[i] = 0;
763 785
        if (ex < 0) max_cap -= ex;
764 786
      }
765 787
      for (int j = 0; j != _res_arc_num; ++j) {
766 788
        if (_upper[j] >= MAX) _upper[j] = max_cap;
767 789
      }
768 790

	
769 791
      // Initialize the large cost vector and the epsilon parameter
770 792
      _epsilon = 0;
771 793
      LargeCost lc;
772 794
      for (int i = 0; i != _root; ++i) {
773 795
        last_out = _first_out[i+1];
774 796
        for (int j = _first_out[i]; j != last_out; ++j) {
775 797
          lc = static_cast<LargeCost>(_scost[j]) * _res_node_num * _alpha;
776 798
          _cost[j] = lc;
777 799
          if (lc > _epsilon) _epsilon = lc;
778 800
        }
779 801
      }
780 802
      _epsilon /= _alpha;
781 803

	
782 804
      // Initialize maps for Circulation and remove non-zero lower bounds
783 805
      ConstMap<Arc, Value> low(0);
784 806
      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
785 807
      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
786 808
      ValueArcMap cap(_graph), flow(_graph);
787 809
      ValueNodeMap sup(_graph);
788 810
      for (NodeIt n(_graph); n != INVALID; ++n) {
789 811
        sup[n] = _supply[_node_id[n]];
790 812
      }
791 813
      if (_have_lower) {
792 814
        for (ArcIt a(_graph); a != INVALID; ++a) {
793 815
          int j = _arc_idf[a];
794 816
          Value c = _lower[j];
795 817
          cap[a] = _upper[j] - c;
796 818
          sup[_graph.source(a)] -= c;
797 819
          sup[_graph.target(a)] += c;
798 820
        }
799 821
      } else {
800 822
        for (ArcIt a(_graph); a != INVALID; ++a) {
801 823
          cap[a] = _upper[_arc_idf[a]];
802 824
        }
803 825
      }
804 826

	
805 827
      // Find a feasible flow using Circulation
806 828
      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
807 829
        circ(_graph, low, cap, sup);
808 830
      if (!circ.flowMap(flow).run()) return INFEASIBLE;
809 831

	
810 832
      // Set residual capacities and handle GEQ supply type
811 833
      if (_sum_supply < 0) {
Show white space 384 line context
... ...
@@ -61,664 +61,693 @@
61 61
  ///
62 62
  /// \tparam GR The digraph type the algorithm runs on.
63 63
  /// \tparam V The number type used for flow amounts, capacity bounds
64 64
  /// and supply values in the algorithm. By default, it is \c int.
65 65
  /// \tparam C The number type used for costs and potentials in the
66 66
  /// algorithm. By default, it is the same as \c V.
67 67
  ///
68 68
  /// \warning Both number types must be signed and all input data must
69 69
  /// be integer.
70 70
  /// \warning This algorithm does not support negative costs for such
71 71
  /// arcs that have infinite upper bound.
72 72
  ///
73 73
  /// \note For more information about the three available methods,
74 74
  /// see \ref Method.
75 75
#ifdef DOXYGEN
76 76
  template <typename GR, typename V, typename C>
77 77
#else
78 78
  template <typename GR, typename V = int, typename C = V>
79 79
#endif
80 80
  class CycleCanceling
81 81
  {
82 82
  public:
83 83

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

	
91 91
  public:
92 92

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

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

	
142 142
  private:
143 143

	
144 144
    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
145 145
    
146 146
    typedef std::vector<int> IntVector;
147 147
    typedef std::vector<char> CharVector;
148 148
    typedef std::vector<double> DoubleVector;
149 149
    typedef std::vector<Value> ValueVector;
150 150
    typedef std::vector<Cost> CostVector;
151 151

	
152 152
  private:
153 153
  
154 154
    template <typename KT, typename VT>
155 155
    class StaticVectorMap {
156 156
    public:
157 157
      typedef KT Key;
158 158
      typedef VT Value;
159 159
      
160 160
      StaticVectorMap(std::vector<Value>& v) : _v(v) {}
161 161
      
162 162
      const Value& operator[](const Key& key) const {
163 163
        return _v[StaticDigraph::id(key)];
164 164
      }
165 165

	
166 166
      Value& operator[](const Key& key) {
167 167
        return _v[StaticDigraph::id(key)];
168 168
      }
169 169
      
170 170
      void set(const Key& key, const Value& val) {
171 171
        _v[StaticDigraph::id(key)] = val;
172 172
      }
173 173

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

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

	
181 181
  private:
182 182

	
183 183

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

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

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

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

	
212 212
    ValueVector _res_cap;
213 213
    CostVector _pi;
214 214

	
215 215
    // Data for a StaticDigraph structure
216 216
    typedef std::pair<int, int> IntPair;
217 217
    StaticDigraph _sgr;
218 218
    std::vector<IntPair> _arc_vec;
219 219
    std::vector<Cost> _cost_vec;
220 220
    IntVector _id_vec;
221 221
    CostArcMap _cost_map;
222 222
    CostNodeMap _pi_map;
223 223
  
224 224
  public:
225 225
  
226 226
    /// \brief Constant for infinite upper bounds (capacities).
227 227
    ///
228 228
    /// Constant for infinite upper bounds (capacities).
229 229
    /// It is \c std::numeric_limits<Value>::infinity() if available,
230 230
    /// \c std::numeric_limits<Value>::max() otherwise.
231 231
    const Value INF;
232 232

	
233 233
  public:
234 234

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

	
253
      // Resize vectors
254
      _node_num = countNodes(_graph);
255
      _arc_num = countArcs(_graph);
256
      _res_node_num = _node_num + 1;
257
      _res_arc_num = 2 * (_arc_num + _node_num);
258
      _root = _node_num;
259

	
260
      _first_out.resize(_res_node_num + 1);
261
      _forward.resize(_res_arc_num);
262
      _source.resize(_res_arc_num);
263
      _target.resize(_res_arc_num);
264
      _reverse.resize(_res_arc_num);
265

	
266
      _lower.resize(_res_arc_num);
267
      _upper.resize(_res_arc_num);
268
      _cost.resize(_res_arc_num);
269
      _supply.resize(_res_node_num);
270
      
271
      _res_cap.resize(_res_arc_num);
272
      _pi.resize(_res_node_num);
273

	
274
      _arc_vec.reserve(_res_arc_num);
275
      _cost_vec.reserve(_res_arc_num);
276
      _id_vec.reserve(_res_arc_num);
277

	
278
      // Copy the graph
279
      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
280
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
281
        _node_id[n] = i;
282
      }
283
      i = 0;
284
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
285
        _first_out[i] = j;
286
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
287
          _arc_idf[a] = j;
288
          _forward[j] = true;
289
          _source[j] = i;
290
          _target[j] = _node_id[_graph.runningNode(a)];
291
        }
292
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
293
          _arc_idb[a] = j;
294
          _forward[j] = false;
295
          _source[j] = i;
296
          _target[j] = _node_id[_graph.runningNode(a)];
297
        }
298
        _forward[j] = false;
299
        _source[j] = i;
300
        _target[j] = _root;
301
        _reverse[j] = k;
302
        _forward[k] = true;
303
        _source[k] = _root;
304
        _target[k] = i;
305
        _reverse[k] = j;
306
        ++j; ++k;
307
      }
308
      _first_out[i] = j;
309
      _first_out[_res_node_num] = k;
310
      for (ArcIt a(_graph); a != INVALID; ++a) {
311
        int fi = _arc_idf[a];
312
        int bi = _arc_idb[a];
313
        _reverse[fi] = bi;
314
        _reverse[bi] = fi;
315
      }
316
      
317
      // Reset parameters
253
      // Reset data structures
318 254
      reset();
319 255
    }
320 256

	
321 257
    /// \name Parameters
322 258
    /// The parameters of the algorithm can be specified using these
323 259
    /// functions.
324 260

	
325 261
    /// @{
326 262

	
327 263
    /// \brief Set the lower bounds on the arcs.
328 264
    ///
329 265
    /// This function sets the lower bounds on the arcs.
330 266
    /// If it is not used before calling \ref run(), the lower bounds
331 267
    /// will be set to zero on all arcs.
332 268
    ///
333 269
    /// \param map An arc map storing the lower bounds.
334 270
    /// Its \c Value type must be convertible to the \c Value type
335 271
    /// of the algorithm.
336 272
    ///
337 273
    /// \return <tt>(*this)</tt>
338 274
    template <typename LowerMap>
339 275
    CycleCanceling& lowerMap(const LowerMap& map) {
340 276
      _have_lower = true;
341 277
      for (ArcIt a(_graph); a != INVALID; ++a) {
342 278
        _lower[_arc_idf[a]] = map[a];
343 279
        _lower[_arc_idb[a]] = map[a];
344 280
      }
345 281
      return *this;
346 282
    }
347 283

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

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

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

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

	
435 371
    /// \name Execution control
436 372
    /// The algorithm can be executed using \ref run().
437 373

	
438 374
    /// @{
439 375

	
440 376
    /// \brief Run the algorithm.
441 377
    ///
442 378
    /// This function runs the algorithm.
443 379
    /// The paramters can be specified using functions \ref lowerMap(),
444 380
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
445 381
    /// For example,
446 382
    /// \code
447 383
    ///   CycleCanceling<ListDigraph> cc(graph);
448 384
    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
449 385
    ///     .supplyMap(sup).run();
450 386
    /// \endcode
451 387
    ///
452
    /// This function can be called more than once. All the parameters
453
    /// that have been given are kept for the next call, unless
454
    /// \ref reset() is called, thus only the modified parameters
455
    /// have to be set again. See \ref reset() for examples.
456
    /// However, the underlying digraph must not be modified after this
457
    /// class have been constructed, since it copies and extends the graph.
388
    /// This function can be called more than once. All the given parameters
389
    /// are kept for the next call, unless \ref resetParams() or \ref reset()
390
    /// is used, thus only the modified parameters have to be set again.
391
    /// If the underlying digraph was also modified after the construction
392
    /// of the class (or the last \ref reset() call), then the \ref reset()
393
    /// function must be called.
458 394
    ///
459 395
    /// \param method The cycle-canceling method that will be used.
460 396
    /// For more information, see \ref Method.
461 397
    ///
462 398
    /// \return \c INFEASIBLE if no feasible flow exists,
463 399
    /// \n \c OPTIMAL if the problem has optimal solution
464 400
    /// (i.e. it is feasible and bounded), and the algorithm has found
465 401
    /// optimal flow and node potentials (primal and dual solutions),
466 402
    /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
467 403
    /// and infinite upper bound. It means that the objective function
468 404
    /// is unbounded on that arc, however, note that it could actually be
469 405
    /// bounded over the feasible flows, but this algroithm cannot handle
470 406
    /// these cases.
471 407
    ///
472 408
    /// \see ProblemType, Method
409
    /// \see resetParams(), reset()
473 410
    ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
474 411
      ProblemType pt = init();
475 412
      if (pt != OPTIMAL) return pt;
476 413
      start(method);
477 414
      return OPTIMAL;
478 415
    }
479 416

	
480 417
    /// \brief Reset all the parameters that have been given before.
481 418
    ///
482 419
    /// This function resets all the paramaters that have been given
483 420
    /// before using functions \ref lowerMap(), \ref upperMap(),
484 421
    /// \ref costMap(), \ref supplyMap(), \ref stSupply().
485 422
    ///
486
    /// It is useful for multiple run() calls. If this function is not
487
    /// used, all the parameters given before are kept for the next
488
    /// \ref run() call.
489
    /// However, the underlying digraph must not be modified after this
490
    /// class have been constructed, since it copies and extends the graph.
423
    /// It is useful for multiple \ref run() calls. Basically, all the given
424
    /// parameters are kept for the next \ref run() call, unless
425
    /// \ref resetParams() or \ref reset() is used.
426
    /// If the underlying digraph was also modified after the construction
427
    /// of the class or the last \ref reset() call, then the \ref reset()
428
    /// function must be used, otherwise \ref resetParams() is sufficient.
491 429
    ///
492 430
    /// For example,
493 431
    /// \code
494 432
    ///   CycleCanceling<ListDigraph> cs(graph);
495 433
    ///
496 434
    ///   // First run
497 435
    ///   cc.lowerMap(lower).upperMap(upper).costMap(cost)
498 436
    ///     .supplyMap(sup).run();
499 437
    ///
500
    ///   // Run again with modified cost map (reset() is not called,
438
    ///   // Run again with modified cost map (resetParams() is not called,
501 439
    ///   // so only the cost map have to be set again)
502 440
    ///   cost[e] += 100;
503 441
    ///   cc.costMap(cost).run();
504 442
    ///
505
    ///   // Run again from scratch using reset()
443
    ///   // Run again from scratch using resetParams()
506 444
    ///   // (the lower bounds will be set to zero on all arcs)
507
    ///   cc.reset();
445
    ///   cc.resetParams();
508 446
    ///   cc.upperMap(capacity).costMap(cost)
509 447
    ///     .supplyMap(sup).run();
510 448
    /// \endcode
511 449
    ///
512 450
    /// \return <tt>(*this)</tt>
513
    CycleCanceling& reset() {
451
    ///
452
    /// \see reset(), run()
453
    CycleCanceling& resetParams() {
514 454
      for (int i = 0; i != _res_node_num; ++i) {
515 455
        _supply[i] = 0;
516 456
      }
517 457
      int limit = _first_out[_root];
518 458
      for (int j = 0; j != limit; ++j) {
519 459
        _lower[j] = 0;
520 460
        _upper[j] = INF;
521 461
        _cost[j] = _forward[j] ? 1 : -1;
522 462
      }
523 463
      for (int j = limit; j != _res_arc_num; ++j) {
524 464
        _lower[j] = 0;
525 465
        _upper[j] = INF;
526 466
        _cost[j] = 0;
527 467
        _cost[_reverse[j]] = 0;
528 468
      }      
529 469
      _have_lower = false;
530 470
      return *this;
531 471
    }
532 472

	
473
    /// \brief Reset the internal data structures and all the parameters
474
    /// that have been given before.
475
    ///
476
    /// This function resets the internal data structures and all the
477
    /// paramaters that have been given before using functions \ref lowerMap(),
478
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
479
    ///
480
    /// It is useful for multiple \ref run() calls. Basically, all the given
481
    /// parameters are kept for the next \ref run() call, unless
482
    /// \ref resetParams() or \ref reset() is used.
483
    /// If the underlying digraph was also modified after the construction
484
    /// of the class or the last \ref reset() call, then the \ref reset()
485
    /// function must be used, otherwise \ref resetParams() is sufficient.
486
    ///
487
    /// See \ref resetParams() for examples.
488
    ///
489
    /// \return <tt>(*this)</tt>
490
    ///
491
    /// \see resetParams(), run()
492
    CycleCanceling& reset() {
493
      // Resize vectors
494
      _node_num = countNodes(_graph);
495
      _arc_num = countArcs(_graph);
496
      _res_node_num = _node_num + 1;
497
      _res_arc_num = 2 * (_arc_num + _node_num);
498
      _root = _node_num;
499

	
500
      _first_out.resize(_res_node_num + 1);
501
      _forward.resize(_res_arc_num);
502
      _source.resize(_res_arc_num);
503
      _target.resize(_res_arc_num);
504
      _reverse.resize(_res_arc_num);
505

	
506
      _lower.resize(_res_arc_num);
507
      _upper.resize(_res_arc_num);
508
      _cost.resize(_res_arc_num);
509
      _supply.resize(_res_node_num);
510
      
511
      _res_cap.resize(_res_arc_num);
512
      _pi.resize(_res_node_num);
513

	
514
      _arc_vec.reserve(_res_arc_num);
515
      _cost_vec.reserve(_res_arc_num);
516
      _id_vec.reserve(_res_arc_num);
517

	
518
      // Copy the graph
519
      int i = 0, j = 0, k = 2 * _arc_num + _node_num;
520
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
521
        _node_id[n] = i;
522
      }
523
      i = 0;
524
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
525
        _first_out[i] = j;
526
        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
527
          _arc_idf[a] = j;
528
          _forward[j] = true;
529
          _source[j] = i;
530
          _target[j] = _node_id[_graph.runningNode(a)];
531
        }
532
        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
533
          _arc_idb[a] = j;
534
          _forward[j] = false;
535
          _source[j] = i;
536
          _target[j] = _node_id[_graph.runningNode(a)];
537
        }
538
        _forward[j] = false;
539
        _source[j] = i;
540
        _target[j] = _root;
541
        _reverse[j] = k;
542
        _forward[k] = true;
543
        _source[k] = _root;
544
        _target[k] = i;
545
        _reverse[k] = j;
546
        ++j; ++k;
547
      }
548
      _first_out[i] = j;
549
      _first_out[_res_node_num] = k;
550
      for (ArcIt a(_graph); a != INVALID; ++a) {
551
        int fi = _arc_idf[a];
552
        int bi = _arc_idb[a];
553
        _reverse[fi] = bi;
554
        _reverse[bi] = fi;
555
      }
556
      
557
      // Reset parameters
558
      resetParams();
559
      return *this;
560
    }
561

	
533 562
    /// @}
534 563

	
535 564
    /// \name Query Functions
536 565
    /// The results of the algorithm can be obtained using these
537 566
    /// functions.\n
538 567
    /// The \ref run() function must be called before using them.
539 568

	
540 569
    /// @{
541 570

	
542 571
    /// \brief Return the total cost of the found flow.
543 572
    ///
544 573
    /// This function returns the total cost of the found flow.
545 574
    /// Its complexity is O(e).
546 575
    ///
547 576
    /// \note The return type of the function can be specified as a
548 577
    /// template parameter. For example,
549 578
    /// \code
550 579
    ///   cc.totalCost<double>();
551 580
    /// \endcode
552 581
    /// It is useful if the total cost cannot be stored in the \c Cost
553 582
    /// type of the algorithm, which is the default return type of the
554 583
    /// function.
555 584
    ///
556 585
    /// \pre \ref run() must be called before using this function.
557 586
    template <typename Number>
558 587
    Number totalCost() const {
559 588
      Number c = 0;
560 589
      for (ArcIt a(_graph); a != INVALID; ++a) {
561 590
        int i = _arc_idb[a];
562 591
        c += static_cast<Number>(_res_cap[i]) *
563 592
             (-static_cast<Number>(_cost[i]));
564 593
      }
565 594
      return c;
566 595
    }
567 596

	
568 597
#ifndef DOXYGEN
569 598
    Cost totalCost() const {
570 599
      return totalCost<Cost>();
571 600
    }
572 601
#endif
573 602

	
574 603
    /// \brief Return the flow on the given arc.
575 604
    ///
576 605
    /// This function returns the flow on the given arc.
577 606
    ///
578 607
    /// \pre \ref run() must be called before using this function.
579 608
    Value flow(const Arc& a) const {
580 609
      return _res_cap[_arc_idb[a]];
581 610
    }
582 611

	
583 612
    /// \brief Return the flow map (the primal solution).
584 613
    ///
585 614
    /// This function copies the flow value on each arc into the given
586 615
    /// map. The \c Value type of the algorithm must be convertible to
587 616
    /// the \c Value type of the map.
588 617
    ///
589 618
    /// \pre \ref run() must be called before using this function.
590 619
    template <typename FlowMap>
591 620
    void flowMap(FlowMap &map) const {
592 621
      for (ArcIt a(_graph); a != INVALID; ++a) {
593 622
        map.set(a, _res_cap[_arc_idb[a]]);
594 623
      }
595 624
    }
596 625

	
597 626
    /// \brief Return the potential (dual value) of the given node.
598 627
    ///
599 628
    /// This function returns the potential (dual value) of the
600 629
    /// given node.
601 630
    ///
602 631
    /// \pre \ref run() must be called before using this function.
603 632
    Cost potential(const Node& n) const {
604 633
      return static_cast<Cost>(_pi[_node_id[n]]);
605 634
    }
606 635

	
607 636
    /// \brief Return the potential map (the dual solution).
608 637
    ///
609 638
    /// This function copies the potential (dual value) of each node
610 639
    /// into the given map.
611 640
    /// The \c Cost type of the algorithm must be convertible to the
612 641
    /// \c Value type of the map.
613 642
    ///
614 643
    /// \pre \ref run() must be called before using this function.
615 644
    template <typename PotentialMap>
616 645
    void potentialMap(PotentialMap &map) const {
617 646
      for (NodeIt n(_graph); n != INVALID; ++n) {
618 647
        map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
619 648
      }
620 649
    }
621 650

	
622 651
    /// @}
623 652

	
624 653
  private:
625 654

	
626 655
    // Initialize the algorithm
627 656
    ProblemType init() {
628 657
      if (_res_node_num <= 1) return INFEASIBLE;
629 658

	
630 659
      // Check the sum of supply values
631 660
      _sum_supply = 0;
632 661
      for (int i = 0; i != _root; ++i) {
633 662
        _sum_supply += _supply[i];
634 663
      }
635 664
      if (_sum_supply > 0) return INFEASIBLE;
636 665
      
637 666

	
638 667
      // Initialize vectors
639 668
      for (int i = 0; i != _res_node_num; ++i) {
640 669
        _pi[i] = 0;
641 670
      }
642 671
      ValueVector excess(_supply);
643 672
      
644 673
      // Remove infinite upper bounds and check negative arcs
645 674
      const Value MAX = std::numeric_limits<Value>::max();
646 675
      int last_out;
647 676
      if (_have_lower) {
648 677
        for (int i = 0; i != _root; ++i) {
649 678
          last_out = _first_out[i+1];
650 679
          for (int j = _first_out[i]; j != last_out; ++j) {
651 680
            if (_forward[j]) {
652 681
              Value c = _cost[j] < 0 ? _upper[j] : _lower[j];
653 682
              if (c >= MAX) return UNBOUNDED;
654 683
              excess[i] -= c;
655 684
              excess[_target[j]] += c;
656 685
            }
657 686
          }
658 687
        }
659 688
      } else {
660 689
        for (int i = 0; i != _root; ++i) {
661 690
          last_out = _first_out[i+1];
662 691
          for (int j = _first_out[i]; j != last_out; ++j) {
663 692
            if (_forward[j] && _cost[j] < 0) {
664 693
              Value c = _upper[j];
665 694
              if (c >= MAX) return UNBOUNDED;
666 695
              excess[i] -= c;
667 696
              excess[_target[j]] += c;
668 697
            }
669 698
          }
670 699
        }
671 700
      }
672 701
      Value ex, max_cap = 0;
673 702
      for (int i = 0; i != _res_node_num; ++i) {
674 703
        ex = excess[i];
675 704
        if (ex < 0) max_cap -= ex;
676 705
      }
677 706
      for (int j = 0; j != _res_arc_num; ++j) {
678 707
        if (_upper[j] >= MAX) _upper[j] = max_cap;
679 708
      }
680 709

	
681 710
      // Initialize maps for Circulation and remove non-zero lower bounds
682 711
      ConstMap<Arc, Value> low(0);
683 712
      typedef typename Digraph::template ArcMap<Value> ValueArcMap;
684 713
      typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
685 714
      ValueArcMap cap(_graph), flow(_graph);
686 715
      ValueNodeMap sup(_graph);
687 716
      for (NodeIt n(_graph); n != INVALID; ++n) {
688 717
        sup[n] = _supply[_node_id[n]];
689 718
      }
690 719
      if (_have_lower) {
691 720
        for (ArcIt a(_graph); a != INVALID; ++a) {
692 721
          int j = _arc_idf[a];
693 722
          Value c = _lower[j];
694 723
          cap[a] = _upper[j] - c;
695 724
          sup[_graph.source(a)] -= c;
696 725
          sup[_graph.target(a)] += c;
697 726
        }
698 727
      } else {
699 728
        for (ArcIt a(_graph); a != INVALID; ++a) {
700 729
          cap[a] = _upper[_arc_idf[a]];
701 730
        }
702 731
      }
703 732

	
704 733
      // Find a feasible flow using Circulation
705 734
      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
706 735
        circ(_graph, low, cap, sup);
707 736
      if (!circ.flowMap(flow).run()) return INFEASIBLE;
708 737

	
709 738
      // Set residual capacities and handle GEQ supply type
710 739
      if (_sum_supply < 0) {
711 740
        for (ArcIt a(_graph); a != INVALID; ++a) {
712 741
          Value fa = flow[a];
713 742
          _res_cap[_arc_idf[a]] = cap[a] - fa;
714 743
          _res_cap[_arc_idb[a]] = fa;
715 744
          sup[_graph.source(a)] -= fa;
716 745
          sup[_graph.target(a)] += fa;
717 746
        }
718 747
        for (NodeIt n(_graph); n != INVALID; ++n) {
719 748
          excess[_node_id[n]] = sup[n];
720 749
        }
721 750
        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
722 751
          int u = _target[a];
723 752
          int ra = _reverse[a];
724 753
          _res_cap[a] = -_sum_supply + 1;
Show white space 384 line context
... ...
@@ -5,384 +5,385 @@
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 46
  /// This algorithm is a highly efficient specialized version of the
47 47
  /// linear programming simplex method directly for the minimum cost
48 48
  /// flow problem.
49 49
  ///
50 50
  /// In general, %NetworkSimplex is the fastest implementation available
51 51
  /// in LEMON for this problem.
52 52
  /// Moreover, it supports both directions of the supply/demand inequality
53 53
  /// constraints. For more information, see \ref SupplyType.
54 54
  ///
55 55
  /// Most of the parameters of the problem (except for the digraph)
56 56
  /// can be given using separate functions, and the algorithm can be
57 57
  /// executed using the \ref run() function. If some parameters are not
58 58
  /// specified, then default values will be used.
59 59
  ///
60 60
  /// \tparam GR The digraph type the algorithm runs on.
61 61
  /// \tparam V The number type used for flow amounts, capacity bounds
62 62
  /// and supply values in the algorithm. By default, it is \c int.
63 63
  /// \tparam C The number type used for costs and potentials in the
64 64
  /// algorithm. By default, it is the same as \c V.
65 65
  ///
66 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 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
    bool _arc_mixing;
197 198

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

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

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

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

	
236 237
  private:
237 238

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

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

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

	
255 256
    public:
256 257

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

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

	
287 288
    }; //class FirstEligiblePivotRule
288 289

	
289 290

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

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

	
304 305
    public:
305 306

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

	
313 314
      // Find next entering arc
314 315
      bool findEnteringArc() {
315 316
        Cost c, min = 0;
316 317
        for (int e = 0; e < _search_arc_num; ++e) {
317 318
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
318 319
          if (c < min) {
319 320
            min = c;
320 321
            _in_arc = e;
321 322
          }
322 323
        }
323 324
        return min < 0;
324 325
      }
325 326

	
326 327
    }; //class BestEligiblePivotRule
327 328

	
328 329

	
329 330
    // Implementation of the Block Search pivot rule
330 331
    class BlockSearchPivotRule
331 332
    {
332 333
    private:
333 334

	
334 335
      // References to the NetworkSimplex class
335 336
      const IntVector  &_source;
336 337
      const IntVector  &_target;
337 338
      const CostVector &_cost;
338 339
      const CharVector &_state;
339 340
      const CostVector &_pi;
340 341
      int &_in_arc;
341 342
      int _search_arc_num;
342 343

	
343 344
      // Pivot rule data
344 345
      int _block_size;
345 346
      int _next_arc;
346 347

	
347 348
    public:
348 349

	
349 350
      // Constructor
350 351
      BlockSearchPivotRule(NetworkSimplex &ns) :
351 352
        _source(ns._source), _target(ns._target),
352 353
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
353 354
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
354 355
        _next_arc(0)
355 356
      {
356 357
        // The main parameters of the pivot rule
357 358
        const double BLOCK_SIZE_FACTOR = 0.5;
358 359
        const int MIN_BLOCK_SIZE = 10;
359 360

	
360 361
        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
361 362
                                    std::sqrt(double(_search_arc_num))),
362 363
                                MIN_BLOCK_SIZE );
363 364
      }
364 365

	
365 366
      // Find next entering arc
366 367
      bool findEnteringArc() {
367 368
        Cost c, min = 0;
368 369
        int cnt = _block_size;
369 370
        int e;
370 371
        for (e = _next_arc; e < _search_arc_num; ++e) {
371 372
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
372 373
          if (c < min) {
373 374
            min = c;
374 375
            _in_arc = e;
375 376
          }
376 377
          if (--cnt == 0) {
377 378
            if (min < 0) goto search_end;
378 379
            cnt = _block_size;
379 380
          }
380 381
        }
381 382
        for (e = 0; e < _next_arc; ++e) {
382 383
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
383 384
          if (c < min) {
384 385
            min = c;
385 386
            _in_arc = e;
386 387
          }
387 388
          if (--cnt == 0) {
388 389
            if (min < 0) goto search_end;
... ...
@@ -444,664 +445,695 @@
444 445
      }
445 446

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

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

	
502 503
    }; //class CandidateListPivotRule
503 504

	
504 505

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

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

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

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

	
537 538
      SortFunc _sort_func;
538 539

	
539 540
    public:
540 541

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

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

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

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

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

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

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

	
621 622
    }; //class AlteringListPivotRule
622 623

	
623 624
  public:
624 625

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

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

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

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

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

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

	
705 656
    /// @{
706 657

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

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

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

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

	
785 736
    /// \brief Set single source and target nodes and a supply value.
786 737
    ///
787 738
    /// This function sets a single source node and a single target node
788 739
    /// and the required flow value.
789 740
    /// If neither this function nor \ref supplyMap() is used before
790 741
    /// calling \ref run(), the supply of each node will be set to zero.
791 742
    ///
792 743
    /// Using this function has the same effect as using \ref supplyMap()
793 744
    /// with such a map in which \c k is assigned to \c s, \c -k is
794 745
    /// assigned to \c t and all other nodes have zero supply value.
795 746
    ///
796 747
    /// \param s The source node.
797 748
    /// \param t The target node.
798 749
    /// \param k The required amount of flow from node \c s to node \c t
799 750
    /// (i.e. the supply of \c s and the demand of \c t).
800 751
    ///
801 752
    /// \return <tt>(*this)</tt>
802 753
    NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) {
803 754
      for (int i = 0; i != _node_num; ++i) {
804 755
        _supply[i] = 0;
805 756
      }
806 757
      _supply[_node_id[s]] =  k;
807 758
      _supply[_node_id[t]] = -k;
808 759
      return *this;
809 760
    }
810 761
    
811 762
    /// \brief Set the type of the supply constraints.
812 763
    ///
813 764
    /// This function sets the type of the supply/demand constraints.
814 765
    /// If it is not used before calling \ref run(), the \ref GEQ supply
815 766
    /// type will be used.
816 767
    ///
817 768
    /// For more information, see \ref SupplyType.
818 769
    ///
819 770
    /// \return <tt>(*this)</tt>
820 771
    NetworkSimplex& supplyType(SupplyType supply_type) {
821 772
      _stype = supply_type;
822 773
      return *this;
823 774
    }
824 775

	
825 776
    /// @}
826 777

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

	
830 781
    /// @{
831 782

	
832 783
    /// \brief Run the algorithm.
833 784
    ///
834 785
    /// This function runs the algorithm.
835 786
    /// The paramters can be specified using functions \ref lowerMap(),
836 787
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 
837 788
    /// \ref supplyType().
838 789
    /// For example,
839 790
    /// \code
840 791
    ///   NetworkSimplex<ListDigraph> ns(graph);
841 792
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
842 793
    ///     .supplyMap(sup).run();
843 794
    /// \endcode
844 795
    ///
845
    /// This function can be called more than once. All the parameters
846
    /// that have been given are kept for the next call, unless
847
    /// \ref reset() is called, thus only the modified parameters
848
    /// have to be set again. See \ref reset() for examples.
849
    /// However, the underlying digraph must not be modified after this
850
    /// class have been constructed, since it copies and extends the graph.
796
    /// This function can be called more than once. All the given parameters
797
    /// are kept for the next call, unless \ref resetParams() or \ref reset()
798
    /// is used, thus only the modified parameters have to be set again.
799
    /// If the underlying digraph was also modified after the construction
800
    /// of the class (or the last \ref reset() call), then the \ref reset()
801
    /// function must be called.
851 802
    ///
852 803
    /// \param pivot_rule The pivot rule that will be used during the
853 804
    /// algorithm. For more information, see \ref PivotRule.
854 805
    ///
855 806
    /// \return \c INFEASIBLE if no feasible flow exists,
856 807
    /// \n \c OPTIMAL if the problem has optimal solution
857 808
    /// (i.e. it is feasible and bounded), and the algorithm has found
858 809
    /// optimal flow and node potentials (primal and dual solutions),
859 810
    /// \n \c UNBOUNDED if the objective function of the problem is
860 811
    /// unbounded, i.e. there is a directed cycle having negative total
861 812
    /// cost and infinite upper bound.
862 813
    ///
863 814
    /// \see ProblemType, PivotRule
815
    /// \see resetParams(), reset()
864 816
    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {
865 817
      if (!init()) return INFEASIBLE;
866 818
      return start(pivot_rule);
867 819
    }
868 820

	
869 821
    /// \brief Reset all the parameters that have been given before.
870 822
    ///
871 823
    /// This function resets all the paramaters that have been given
872 824
    /// before using functions \ref lowerMap(), \ref upperMap(),
873 825
    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType().
874 826
    ///
875
    /// It is useful for multiple run() calls. If this function is not
876
    /// used, all the parameters given before are kept for the next
877
    /// \ref run() call.
878
    /// However, the underlying digraph must not be modified after this
879
    /// class have been constructed, since it copies and extends the graph.
827
    /// It is useful for multiple \ref run() calls. Basically, all the given
828
    /// parameters are kept for the next \ref run() call, unless
829
    /// \ref resetParams() or \ref reset() is used.
830
    /// If the underlying digraph was also modified after the construction
831
    /// of the class or the last \ref reset() call, then the \ref reset()
832
    /// function must be used, otherwise \ref resetParams() is sufficient.
880 833
    ///
881 834
    /// For example,
882 835
    /// \code
883 836
    ///   NetworkSimplex<ListDigraph> ns(graph);
884 837
    ///
885 838
    ///   // First run
886 839
    ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
887 840
    ///     .supplyMap(sup).run();
888 841
    ///
889
    ///   // Run again with modified cost map (reset() is not called,
842
    ///   // Run again with modified cost map (resetParams() is not called,
890 843
    ///   // so only the cost map have to be set again)
891 844
    ///   cost[e] += 100;
892 845
    ///   ns.costMap(cost).run();
893 846
    ///
894
    ///   // Run again from scratch using reset()
847
    ///   // Run again from scratch using resetParams()
895 848
    ///   // (the lower bounds will be set to zero on all arcs)
896
    ///   ns.reset();
849
    ///   ns.resetParams();
897 850
    ///   ns.upperMap(capacity).costMap(cost)
898 851
    ///     .supplyMap(sup).run();
899 852
    /// \endcode
900 853
    ///
901 854
    /// \return <tt>(*this)</tt>
902
    NetworkSimplex& reset() {
855
    ///
856
    /// \see reset(), run()
857
    NetworkSimplex& resetParams() {
903 858
      for (int i = 0; i != _node_num; ++i) {
904 859
        _supply[i] = 0;
905 860
      }
906 861
      for (int i = 0; i != _arc_num; ++i) {
907 862
        _lower[i] = 0;
908 863
        _upper[i] = INF;
909 864
        _cost[i] = 1;
910 865
      }
911 866
      _have_lower = false;
912 867
      _stype = GEQ;
913 868
      return *this;
914 869
    }
915 870

	
871
    /// \brief Reset the internal data structures and all the parameters
872
    /// that have been given before.
873
    ///
874
    /// This function resets the internal data structures and all the
875
    /// paramaters that have been given before using functions \ref lowerMap(),
876
    /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
877
    /// \ref supplyType().
878
    ///
879
    /// It is useful for multiple \ref run() calls. Basically, all the given
880
    /// parameters are kept for the next \ref run() call, unless
881
    /// \ref resetParams() or \ref reset() is used.
882
    /// If the underlying digraph was also modified after the construction
883
    /// of the class or the last \ref reset() call, then the \ref reset()
884
    /// function must be used, otherwise \ref resetParams() is sufficient.
885
    ///
886
    /// See \ref resetParams() for examples.
887
    ///
888
    /// \return <tt>(*this)</tt>
889
    ///
890
    /// \see resetParams(), run()
891
    NetworkSimplex& reset() {
892
      // Resize vectors
893
      _node_num = countNodes(_graph);
894
      _arc_num = countArcs(_graph);
895
      int all_node_num = _node_num + 1;
896
      int max_arc_num = _arc_num + 2 * _node_num;
897

	
898
      _source.resize(max_arc_num);
899
      _target.resize(max_arc_num);
900

	
901
      _lower.resize(_arc_num);
902
      _upper.resize(_arc_num);
903
      _cap.resize(max_arc_num);
904
      _cost.resize(max_arc_num);
905
      _supply.resize(all_node_num);
906
      _flow.resize(max_arc_num);
907
      _pi.resize(all_node_num);
908

	
909
      _parent.resize(all_node_num);
910
      _pred.resize(all_node_num);
911
      _forward.resize(all_node_num);
912
      _thread.resize(all_node_num);
913
      _rev_thread.resize(all_node_num);
914
      _succ_num.resize(all_node_num);
915
      _last_succ.resize(all_node_num);
916
      _state.resize(max_arc_num);
917

	
918
      // Copy the graph
919
      int i = 0;
920
      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
921
        _node_id[n] = i;
922
      }
923
      if (_arc_mixing) {
924
        // Store the arcs in a mixed order
925
        int k = std::max(int(std::sqrt(double(_arc_num))), 10);
926
        int i = 0, j = 0;
927
        for (ArcIt a(_graph); a != INVALID; ++a) {
928
          _arc_id[a] = i;
929
          _source[i] = _node_id[_graph.source(a)];
930
          _target[i] = _node_id[_graph.target(a)];
931
          if ((i += k) >= _arc_num) i = ++j;
932
        }
933
      } else {
934
        // Store the arcs in the original order
935
        int i = 0;
936
        for (ArcIt a(_graph); a != INVALID; ++a, ++i) {
937
          _arc_id[a] = i;
938
          _source[i] = _node_id[_graph.source(a)];
939
          _target[i] = _node_id[_graph.target(a)];
940
        }
941
      }
942
      
943
      // Reset parameters
944
      resetParams();
945
      return *this;
946
    }
947
    
916 948
    /// @}
917 949

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

	
923 955
    /// @{
924 956

	
925 957
    /// \brief Return the total cost of the found flow.
926 958
    ///
927 959
    /// This function returns the total cost of the found flow.
928 960
    /// Its complexity is O(e).
929 961
    ///
930 962
    /// \note The return type of the function can be specified as a
931 963
    /// template parameter. For example,
932 964
    /// \code
933 965
    ///   ns.totalCost<double>();
934 966
    /// \endcode
935 967
    /// It is useful if the total cost cannot be stored in the \c Cost
936 968
    /// type of the algorithm, which is the default return type of the
937 969
    /// function.
938 970
    ///
939 971
    /// \pre \ref run() must be called before using this function.
940 972
    template <typename Number>
941 973
    Number totalCost() const {
942 974
      Number c = 0;
943 975
      for (ArcIt a(_graph); a != INVALID; ++a) {
944 976
        int i = _arc_id[a];
945 977
        c += Number(_flow[i]) * Number(_cost[i]);
946 978
      }
947 979
      return c;
948 980
    }
949 981

	
950 982
#ifndef DOXYGEN
951 983
    Cost totalCost() const {
952 984
      return totalCost<Cost>();
953 985
    }
954 986
#endif
955 987

	
956 988
    /// \brief Return the flow on the given arc.
957 989
    ///
958 990
    /// This function returns the flow on the given arc.
959 991
    ///
960 992
    /// \pre \ref run() must be called before using this function.
961 993
    Value flow(const Arc& a) const {
962 994
      return _flow[_arc_id[a]];
963 995
    }
964 996

	
965 997
    /// \brief Return the flow map (the primal solution).
966 998
    ///
967 999
    /// This function copies the flow value on each arc into the given
968 1000
    /// map. The \c Value type of the algorithm must be convertible to
969 1001
    /// the \c Value type of the map.
970 1002
    ///
971 1003
    /// \pre \ref run() must be called before using this function.
972 1004
    template <typename FlowMap>
973 1005
    void flowMap(FlowMap &map) const {
974 1006
      for (ArcIt a(_graph); a != INVALID; ++a) {
975 1007
        map.set(a, _flow[_arc_id[a]]);
976 1008
      }
977 1009
    }
978 1010

	
979 1011
    /// \brief Return the potential (dual value) of the given node.
980 1012
    ///
981 1013
    /// This function returns the potential (dual value) of the
982 1014
    /// given node.
983 1015
    ///
984 1016
    /// \pre \ref run() must be called before using this function.
985 1017
    Cost potential(const Node& n) const {
986 1018
      return _pi[_node_id[n]];
987 1019
    }
988 1020

	
989 1021
    /// \brief Return the potential map (the dual solution).
990 1022
    ///
991 1023
    /// This function copies the potential (dual value) of each node
992 1024
    /// into the given map.
993 1025
    /// The \c Cost type of the algorithm must be convertible to the
994 1026
    /// \c Value type of the map.
995 1027
    ///
996 1028
    /// \pre \ref run() must be called before using this function.
997 1029
    template <typename PotentialMap>
998 1030
    void potentialMap(PotentialMap &map) const {
999 1031
      for (NodeIt n(_graph); n != INVALID; ++n) {
1000 1032
        map.set(n, _pi[_node_id[n]]);
1001 1033
      }
1002 1034
    }
1003 1035

	
1004 1036
    /// @}
1005 1037

	
1006 1038
  private:
1007 1039

	
1008 1040
    // Initialize internal data structures
1009 1041
    bool init() {
1010 1042
      if (_node_num == 0) return false;
1011 1043

	
1012 1044
      // Check the sum of supply values
1013 1045
      _sum_supply = 0;
1014 1046
      for (int i = 0; i != _node_num; ++i) {
1015 1047
        _sum_supply += _supply[i];
1016 1048
      }
1017 1049
      if ( !((_stype == GEQ && _sum_supply <= 0) ||
1018 1050
             (_stype == LEQ && _sum_supply >= 0)) ) return false;
1019 1051

	
1020 1052
      // Remove non-zero lower bounds
1021 1053
      if (_have_lower) {
1022 1054
        for (int i = 0; i != _arc_num; ++i) {
1023 1055
          Value c = _lower[i];
1024 1056
          if (c >= 0) {
1025 1057
            _cap[i] = _upper[i] < MAX ? _upper[i] - c : INF;
1026 1058
          } else {
1027 1059
            _cap[i] = _upper[i] < MAX + c ? _upper[i] - c : INF;
1028 1060
          }
1029 1061
          _supply[_source[i]] -= c;
1030 1062
          _supply[_target[i]] += c;
1031 1063
        }
1032 1064
      } else {
1033 1065
        for (int i = 0; i != _arc_num; ++i) {
1034 1066
          _cap[i] = _upper[i];
1035 1067
        }
1036 1068
      }
1037 1069

	
1038 1070
      // Initialize artifical cost
1039 1071
      Cost ART_COST;
1040 1072
      if (std::numeric_limits<Cost>::is_exact) {
1041 1073
        ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
1042 1074
      } else {
1043 1075
        ART_COST = std::numeric_limits<Cost>::min();
1044 1076
        for (int i = 0; i != _arc_num; ++i) {
1045 1077
          if (_cost[i] > ART_COST) ART_COST = _cost[i];
1046 1078
        }
1047 1079
        ART_COST = (ART_COST + 1) * _node_num;
1048 1080
      }
1049 1081

	
1050 1082
      // Initialize arc maps
1051 1083
      for (int i = 0; i != _arc_num; ++i) {
1052 1084
        _flow[i] = 0;
1053 1085
        _state[i] = STATE_LOWER;
1054 1086
      }
1055 1087
      
1056 1088
      // Set data for the artificial root node
1057 1089
      _root = _node_num;
1058 1090
      _parent[_root] = -1;
1059 1091
      _pred[_root] = -1;
1060 1092
      _thread[_root] = 0;
1061 1093
      _rev_thread[0] = _root;
1062 1094
      _succ_num[_root] = _node_num + 1;
1063 1095
      _last_succ[_root] = _root - 1;
1064 1096
      _supply[_root] = -_sum_supply;
1065 1097
      _pi[_root] = 0;
1066 1098

	
1067 1099
      // Add artificial arcs and initialize the spanning tree data structure
1068 1100
      if (_sum_supply == 0) {
1069 1101
        // EQ supply constraints
1070 1102
        _search_arc_num = _arc_num;
1071 1103
        _all_arc_num = _arc_num + _node_num;
1072 1104
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1073 1105
          _parent[u] = _root;
1074 1106
          _pred[u] = e;
1075 1107
          _thread[u] = u + 1;
1076 1108
          _rev_thread[u + 1] = u;
1077 1109
          _succ_num[u] = 1;
1078 1110
          _last_succ[u] = u;
1079 1111
          _cap[e] = INF;
1080 1112
          _state[e] = STATE_TREE;
1081 1113
          if (_supply[u] >= 0) {
1082 1114
            _forward[u] = true;
1083 1115
            _pi[u] = 0;
1084 1116
            _source[e] = u;
1085 1117
            _target[e] = _root;
1086 1118
            _flow[e] = _supply[u];
1087 1119
            _cost[e] = 0;
1088 1120
          } else {
1089 1121
            _forward[u] = false;
1090 1122
            _pi[u] = ART_COST;
1091 1123
            _source[e] = _root;
1092 1124
            _target[e] = u;
1093 1125
            _flow[e] = -_supply[u];
1094 1126
            _cost[e] = ART_COST;
1095 1127
          }
1096 1128
        }
1097 1129
      }
1098 1130
      else if (_sum_supply > 0) {
1099 1131
        // LEQ supply constraints
1100 1132
        _search_arc_num = _arc_num + _node_num;
1101 1133
        int f = _arc_num + _node_num;
1102 1134
        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
1103 1135
          _parent[u] = _root;
1104 1136
          _thread[u] = u + 1;
1105 1137
          _rev_thread[u + 1] = u;
1106 1138
          _succ_num[u] = 1;
1107 1139
          _last_succ[u] = u;
Show 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
#include <iostream>
20 20
#include <fstream>
21 21
#include <limits>
22 22

	
23 23
#include <lemon/list_graph.h>
24 24
#include <lemon/lgf_reader.h>
25 25

	
26 26
#include <lemon/network_simplex.h>
27 27
#include <lemon/capacity_scaling.h>
28 28
#include <lemon/cost_scaling.h>
29 29
#include <lemon/cycle_canceling.h>
30 30

	
31 31
#include <lemon/concepts/digraph.h>
32 32
#include <lemon/concepts/heap.h>
33 33
#include <lemon/concept_check.h>
34 34

	
35 35
#include "test_tools.h"
36 36

	
37 37
using namespace lemon;
38 38

	
39 39
// Test networks
40 40
char test_lgf[] =
41 41
  "@nodes\n"
42 42
  "label  sup1 sup2 sup3 sup4 sup5 sup6\n"
43 43
  "    1    20   27    0   30   20   30\n"
44 44
  "    2    -4    0    0    0   -8   -3\n"
45 45
  "    3     0    0    0    0    0    0\n"
46 46
  "    4     0    0    0    0    0    0\n"
47 47
  "    5     9    0    0    0    6   11\n"
48 48
  "    6    -6    0    0    0   -5   -6\n"
49 49
  "    7     0    0    0    0    0    0\n"
50 50
  "    8     0    0    0    0    0    3\n"
51 51
  "    9     3    0    0    0    0    0\n"
52 52
  "   10    -2    0    0    0   -7   -2\n"
53 53
  "   11     0    0    0    0  -10    0\n"
54 54
  "   12   -20  -27    0  -30  -30  -20\n"
55 55
  "\n"                
56 56
  "@arcs\n"
57 57
  "       cost  cap low1 low2 low3\n"
58 58
  " 1  2    70   11    0    8    8\n"
59 59
  " 1  3   150    3    0    1    0\n"
60 60
  " 1  4    80   15    0    2    2\n"
61 61
  " 2  8    80   12    0    0    0\n"
62 62
  " 3  5   140    5    0    3    1\n"
63 63
  " 4  6    60   10    0    1    0\n"
64 64
  " 4  7    80    2    0    0    0\n"
65 65
  " 4  8   110    3    0    0    0\n"
66 66
  " 5  7    60   14    0    0    0\n"
67 67
  " 5 11   120   12    0    0    0\n"
68 68
  " 6  3     0    3    0    0    0\n"
69 69
  " 6  9   140    4    0    0    0\n"
70 70
  " 6 10    90    8    0    0    0\n"
71 71
  " 7  1    30    5    0    0   -5\n"
72 72
  " 8 12    60   16    0    4    3\n"
73 73
  " 9 12    50    6    0    0    0\n"
74 74
  "10 12    70   13    0    5    2\n"
75 75
  "10  2   100    7    0    0    0\n"
76 76
  "10  7    60   10    0    0   -3\n"
77 77
  "11 10    20   14    0    6  -20\n"
78 78
  "12 11    30   10    0    0  -10\n"
79 79
  "\n"
80 80
  "@attributes\n"
81 81
  "source 1\n"
82 82
  "target 12\n";
83 83

	
84 84
char test_neg1_lgf[] =
85 85
  "@nodes\n"
86 86
  "label   sup\n"
87 87
  "    1   100\n"
88 88
  "    2     0\n"
89 89
  "    3     0\n"
90 90
  "    4  -100\n"
91 91
  "    5     0\n"
92 92
  "    6     0\n"
93 93
  "    7     0\n"
94 94
  "@arcs\n"
95 95
  "      cost   low1   low2\n"
96 96
  "1 2    100      0      0\n"
97 97
  "1 3     30      0      0\n"
98 98
  "2 4     20      0      0\n"
99 99
  "3 4     80      0      0\n"
100 100
  "3 2     50      0      0\n"
101 101
  "5 3     10      0      0\n"
102 102
  "5 6     80      0   1000\n"
103 103
  "6 7     30      0  -1000\n"
104 104
  "7 5   -120      0      0\n";
105 105
  
106 106
char test_neg2_lgf[] =
107 107
  "@nodes\n"
108 108
  "label   sup\n"
109 109
  "    1   100\n"
110 110
  "    2  -300\n"
111 111
  "@arcs\n"
112 112
  "      cost\n"
113 113
  "1 2     -1\n";
114 114

	
115 115

	
116 116
// Test data
117 117
typedef ListDigraph Digraph;
118 118
DIGRAPH_TYPEDEFS(ListDigraph);
119 119

	
120 120
Digraph gr;
121 121
Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
122 122
Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
123 123
ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
124 124
Node v, w;
125 125

	
126 126
Digraph neg1_gr;
127 127
Digraph::ArcMap<int> neg1_c(neg1_gr), neg1_l1(neg1_gr), neg1_l2(neg1_gr);
128 128
ConstMap<Arc, int> neg1_u1(std::numeric_limits<int>::max()), neg1_u2(5000);
129 129
Digraph::NodeMap<int> neg1_s(neg1_gr);
130 130

	
131 131
Digraph neg2_gr;
132 132
Digraph::ArcMap<int> neg2_c(neg2_gr);
133 133
ConstMap<Arc, int> neg2_l(0), neg2_u(1000);
134 134
Digraph::NodeMap<int> neg2_s(neg2_gr);
135 135

	
136 136

	
137 137
enum SupplyType {
138 138
  EQ,
139 139
  GEQ,
140 140
  LEQ
141 141
};
142 142

	
143 143

	
144 144
// Check the interface of an MCF algorithm
145 145
template <typename GR, typename Value, typename Cost>
146 146
class McfClassConcept
147 147
{
148 148
public:
149 149

	
150 150
  template <typename MCF>
151 151
  struct Constraints {
152 152
    void constraints() {
153 153
      checkConcept<concepts::Digraph, GR>();
154 154
      
155 155
      const Constraints& me = *this;
156 156

	
157 157
      MCF mcf(me.g);
158 158
      const MCF& const_mcf = mcf;
159 159

	
160
      b = mcf.reset()
160
      b = mcf.reset().resetParams()
161 161
             .lowerMap(me.lower)
162 162
             .upperMap(me.upper)
163 163
             .costMap(me.cost)
164 164
             .supplyMap(me.sup)
165 165
             .stSupply(me.n, me.n, me.k)
166 166
             .run();
167 167

	
168 168
      c = const_mcf.totalCost();
169 169
      x = const_mcf.template totalCost<double>();
170 170
      v = const_mcf.flow(me.a);
171 171
      c = const_mcf.potential(me.n);
172 172
      const_mcf.flowMap(fm);
173 173
      const_mcf.potentialMap(pm);
174 174
    }
175 175

	
176 176
    typedef typename GR::Node Node;
177 177
    typedef typename GR::Arc Arc;
178 178
    typedef concepts::ReadMap<Node, Value> NM;
179 179
    typedef concepts::ReadMap<Arc, Value> VAM;
180 180
    typedef concepts::ReadMap<Arc, Cost> CAM;
181 181
    typedef concepts::WriteMap<Arc, Value> FlowMap;
182 182
    typedef concepts::WriteMap<Node, Cost> PotMap;
183 183
  
184 184
    GR g;
185 185
    VAM lower;
186 186
    VAM upper;
187 187
    CAM cost;
188 188
    NM sup;
189 189
    Node n;
190 190
    Arc a;
191 191
    Value k;
192 192

	
193 193
    FlowMap fm;
194 194
    PotMap pm;
195 195
    bool b;
196 196
    double x;
197 197
    typename MCF::Value v;
198 198
    typename MCF::Cost c;
199 199
  };
200 200

	
201 201
};
202 202

	
203 203

	
204 204
// Check the feasibility of the given flow (primal soluiton)
205 205
template < typename GR, typename LM, typename UM,
206 206
           typename SM, typename FM >
207 207
bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
208 208
                const SM& supply, const FM& flow,
209 209
                SupplyType type = EQ )
210 210
{
211 211
  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
212 212

	
213 213
  for (ArcIt e(gr); e != INVALID; ++e) {
214 214
    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
215 215
  }
216 216

	
217 217
  for (NodeIt n(gr); n != INVALID; ++n) {
218 218
    typename SM::Value sum = 0;
219 219
    for (OutArcIt e(gr, n); e != INVALID; ++e)
220 220
      sum += flow[e];
221 221
    for (InArcIt e(gr, n); e != INVALID; ++e)
222 222
      sum -= flow[e];
223 223
    bool b = (type ==  EQ && sum == supply[n]) ||
224 224
             (type == GEQ && sum >= supply[n]) ||
225 225
             (type == LEQ && sum <= supply[n]);
226 226
    if (!b) return false;
227 227
  }
228 228

	
229 229
  return true;
230 230
}
231 231

	
232 232
// Check the feasibility of the given potentials (dual soluiton)
233 233
// using the "Complementary Slackness" optimality condition
234 234
template < typename GR, typename LM, typename UM,
235 235
           typename CM, typename SM, typename FM, typename PM >
236 236
bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
237 237
                     const CM& cost, const SM& supply, const FM& flow, 
238 238
                     const PM& pi, SupplyType type )
239 239
{
240 240
  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
241 241

	
242 242
  bool opt = true;
243 243
  for (ArcIt e(gr); opt && e != INVALID; ++e) {
244 244
    typename CM::Value red_cost =
245 245
      cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
246 246
    opt = red_cost == 0 ||
247 247
          (red_cost > 0 && flow[e] == lower[e]) ||
248 248
          (red_cost < 0 && flow[e] == upper[e]);
249 249
  }
250 250
  
251 251
  for (NodeIt n(gr); opt && n != INVALID; ++n) {
252 252
    typename SM::Value sum = 0;
253 253
    for (OutArcIt e(gr, n); e != INVALID; ++e)
254 254
      sum += flow[e];
255 255
    for (InArcIt e(gr, n); e != INVALID; ++e)
256 256
      sum -= flow[e];
257 257
    if (type != LEQ) {
258 258
      opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0);
259 259
    } else {
260 260
      opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0);
261 261
    }
262 262
  }
263 263
  
264 264
  return opt;
265 265
}
266 266

	
267 267
// Check whether the dual cost is equal to the primal cost
268 268
template < typename GR, typename LM, typename UM,
269 269
           typename CM, typename SM, typename PM >
270 270
bool checkDualCost( const GR& gr, const LM& lower, const UM& upper,
271 271
                    const CM& cost, const SM& supply, const PM& pi,
272 272
                    typename CM::Value total )
273 273
{
274 274
  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
275 275

	
276 276
  typename CM::Value dual_cost = 0;
277 277
  SM red_supply(gr);
278 278
  for (NodeIt n(gr); n != INVALID; ++n) {
279 279
    red_supply[n] = supply[n];
280 280
  }
281 281
  for (ArcIt a(gr); a != INVALID; ++a) {
282 282
    if (lower[a] != 0) {
283 283
      dual_cost += lower[a] * cost[a];
284 284
      red_supply[gr.source(a)] -= lower[a];
285 285
      red_supply[gr.target(a)] += lower[a];
286 286
    }
287 287
  }
288 288
  
289 289
  for (NodeIt n(gr); n != INVALID; ++n) {
290 290
    dual_cost -= red_supply[n] * pi[n];
291 291
  }
292 292
  for (ArcIt a(gr); a != INVALID; ++a) {
293 293
    typename CM::Value red_cost =
294 294
      cost[a] + pi[gr.source(a)] - pi[gr.target(a)];
295 295
    dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0);
296 296
  }
297 297
  
298 298
  return dual_cost == total;
299 299
}
300 300

	
301 301
// Run a minimum cost flow algorithm and check the results
302 302
template < typename MCF, typename GR,
303 303
           typename LM, typename UM,
304 304
           typename CM, typename SM,
305 305
           typename PT >
306 306
void checkMcf( const MCF& mcf, PT mcf_result,
307 307
               const GR& gr, const LM& lower, const UM& upper,
308 308
               const CM& cost, const SM& supply,
309 309
               PT result, bool optimal, typename CM::Value total,
310 310
               const std::string &test_id = "",
311 311
               SupplyType type = EQ )
312 312
{
313 313
  check(mcf_result == result, "Wrong result " + test_id);
314 314
  if (optimal) {
315 315
    typename GR::template ArcMap<typename SM::Value> flow(gr);
316 316
    typename GR::template NodeMap<typename CM::Value> pi(gr);
317 317
    mcf.flowMap(flow);
318 318
    mcf.potentialMap(pi);
319 319
    check(checkFlow(gr, lower, upper, supply, flow, type),
320 320
          "The flow is not feasible " + test_id);
321 321
    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
322 322
    check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type),
323 323
          "Wrong potentials " + test_id);
324 324
    check(checkDualCost(gr, lower, upper, cost, supply, pi, total),
325 325
          "Wrong dual cost " + test_id);
326 326
  }
327 327
}
328 328

	
329 329
template < typename MCF, typename Param >
330 330
void runMcfGeqTests( Param param,
331 331
                     const std::string &test_str = "",
332 332
                     bool full_neg_cost_support = false )
333 333
{
334 334
  MCF mcf1(gr), mcf2(neg1_gr), mcf3(neg2_gr);
335 335
  
336 336
  // Basic tests
337 337
  mcf1.upperMap(u).costMap(c).supplyMap(s1);
338 338
  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s1,
339 339
           mcf1.OPTIMAL, true,     5240, test_str + "-1");
340 340
  mcf1.stSupply(v, w, 27);
341 341
  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s2,
342 342
           mcf1.OPTIMAL, true,     7620, test_str + "-2");
343 343
  mcf1.lowerMap(l2).supplyMap(s1);
344 344
  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s1,
345 345
           mcf1.OPTIMAL, true,     5970, test_str + "-3");
346 346
  mcf1.stSupply(v, w, 27);
347 347
  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s2,
348 348
           mcf1.OPTIMAL, true,     8010, test_str + "-4");
349
  mcf1.reset().supplyMap(s1);
349
  mcf1.resetParams().supplyMap(s1);
350 350
  checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s1,
351 351
           mcf1.OPTIMAL, true,       74, test_str + "-5");
352 352
  mcf1.lowerMap(l2).stSupply(v, w, 27);
353 353
  checkMcf(mcf1, mcf1.run(param), gr, l2, cu, cc, s2,
354 354
           mcf1.OPTIMAL, true,       94, test_str + "-6");
355 355
  mcf1.reset();
356 356
  checkMcf(mcf1, mcf1.run(param), gr, l1, cu, cc, s3,
357 357
           mcf1.OPTIMAL, true,        0, test_str + "-7");
358 358
  mcf1.lowerMap(l2).upperMap(u);
359 359
  checkMcf(mcf1, mcf1.run(param), gr, l2, u, cc, s3,
360 360
           mcf1.INFEASIBLE, false,    0, test_str + "-8");
361 361
  mcf1.lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
362 362
  checkMcf(mcf1, mcf1.run(param), gr, l3, u, c, s4,
363 363
           mcf1.OPTIMAL, true,     6360, test_str + "-9");
364 364

	
365 365
  // Tests for the GEQ form
366
  mcf1.reset().upperMap(u).costMap(c).supplyMap(s5);
366
  mcf1.resetParams().upperMap(u).costMap(c).supplyMap(s5);
367 367
  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s5,
368 368
           mcf1.OPTIMAL, true,     3530, test_str + "-10", GEQ);
369 369
  mcf1.lowerMap(l2);
370 370
  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,
371 371
           mcf1.OPTIMAL, true,     4540, test_str + "-11", GEQ);
372 372
  mcf1.supplyMap(s6);
373 373
  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,
374 374
           mcf1.INFEASIBLE, false,    0, test_str + "-12", GEQ);
375 375

	
376 376
  // Tests with negative costs
377 377
  mcf2.lowerMap(neg1_l1).costMap(neg1_c).supplyMap(neg1_s);
378 378
  checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u1, neg1_c, neg1_s,
379 379
           mcf2.UNBOUNDED, false,     0, test_str + "-13");
380 380
  mcf2.upperMap(neg1_u2);
381 381
  checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l1, neg1_u2, neg1_c, neg1_s,
382 382
           mcf2.OPTIMAL, true,   -40000, test_str + "-14");
383
  mcf2.reset().lowerMap(neg1_l2).costMap(neg1_c).supplyMap(neg1_s);
383
  mcf2.resetParams().lowerMap(neg1_l2).costMap(neg1_c).supplyMap(neg1_s);
384 384
  checkMcf(mcf2, mcf2.run(param), neg1_gr, neg1_l2, neg1_u1, neg1_c, neg1_s,
385 385
           mcf2.UNBOUNDED, false,     0, test_str + "-15");
386 386

	
387 387
  mcf3.costMap(neg2_c).supplyMap(neg2_s);
388 388
  if (full_neg_cost_support) {
389 389
    checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
390 390
             mcf3.OPTIMAL, true,   -300, test_str + "-16", GEQ);
391 391
  } else {
392 392
    checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
393 393
             mcf3.UNBOUNDED, false,   0, test_str + "-17", GEQ);
394 394
  }
395 395
  mcf3.upperMap(neg2_u);
396 396
  checkMcf(mcf3, mcf3.run(param), neg2_gr, neg2_l, neg2_u, neg2_c, neg2_s,
397 397
           mcf3.OPTIMAL, true,     -300, test_str + "-18", GEQ);
398 398
}
399 399

	
400 400
template < typename MCF, typename Param >
401 401
void runMcfLeqTests( Param param,
402 402
                     const std::string &test_str = "" )
403 403
{
404 404
  // Tests for the LEQ form
405 405
  MCF mcf1(gr);
406 406
  mcf1.supplyType(mcf1.LEQ);
407 407
  mcf1.upperMap(u).costMap(c).supplyMap(s6);
408 408
  checkMcf(mcf1, mcf1.run(param), gr, l1, u, c, s6,
409 409
           mcf1.OPTIMAL, true,   5080, test_str + "-19", LEQ);
410 410
  mcf1.lowerMap(l2);
411 411
  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s6,
412 412
           mcf1.OPTIMAL, true,   5930, test_str + "-20", LEQ);
413 413
  mcf1.supplyMap(s5);
414 414
  checkMcf(mcf1, mcf1.run(param), gr, l2, u, c, s5,
415 415
           mcf1.INFEASIBLE, false,  0, test_str + "-21", LEQ);
416 416
}
417 417

	
418 418

	
419 419
int main()
420 420
{
421 421
  // Read the test networks
422 422
  std::istringstream input(test_lgf);
423 423
  DigraphReader<Digraph>(gr, input)
424 424
    .arcMap("cost", c)
425 425
    .arcMap("cap", u)
426 426
    .arcMap("low1", l1)
427 427
    .arcMap("low2", l2)
428 428
    .arcMap("low3", l3)
429 429
    .nodeMap("sup1", s1)
430 430
    .nodeMap("sup2", s2)
431 431
    .nodeMap("sup3", s3)
432 432
    .nodeMap("sup4", s4)
433 433
    .nodeMap("sup5", s5)
434 434
    .nodeMap("sup6", s6)
435 435
    .node("source", v)
436 436
    .node("target", w)
437 437
    .run();
438 438
  
439 439
  std::istringstream neg_inp1(test_neg1_lgf);
440 440
  DigraphReader<Digraph>(neg1_gr, neg_inp1)
441 441
    .arcMap("cost", neg1_c)
442 442
    .arcMap("low1", neg1_l1)
443 443
    .arcMap("low2", neg1_l2)
444 444
    .nodeMap("sup", neg1_s)
445 445
    .run();
446 446
  
447 447
  std::istringstream neg_inp2(test_neg2_lgf);
448 448
  DigraphReader<Digraph>(neg2_gr, neg_inp2)
449 449
    .arcMap("cost", neg2_c)
450 450
    .nodeMap("sup", neg2_s)
451 451
    .run();
452 452
  
453 453
  // Check the interface of NetworkSimplex
454 454
  {
455 455
    typedef concepts::Digraph GR;
456 456
    checkConcept< McfClassConcept<GR, int, int>,
457 457
                  NetworkSimplex<GR> >();
458 458
    checkConcept< McfClassConcept<GR, double, double>,
459 459
                  NetworkSimplex<GR, double> >();
460 460
    checkConcept< McfClassConcept<GR, int, double>,
461 461
                  NetworkSimplex<GR, int, double> >();
462 462
  }
463 463

	
464 464
  // Check the interface of CapacityScaling
465 465
  {
466 466
    typedef concepts::Digraph GR;
467 467
    checkConcept< McfClassConcept<GR, int, int>,
468 468
                  CapacityScaling<GR> >();
469 469
    checkConcept< McfClassConcept<GR, double, double>,
470 470
                  CapacityScaling<GR, double> >();
471 471
    checkConcept< McfClassConcept<GR, int, double>,
472 472
                  CapacityScaling<GR, int, double> >();
473 473
    typedef CapacityScaling<GR>::
474 474
      SetHeap<concepts::Heap<int, RangeMap<int> > >::Create CAS;
475 475
    checkConcept< McfClassConcept<GR, int, int>, CAS >();
476 476
  }
477 477

	
478 478
  // Check the interface of CostScaling
479 479
  {
480 480
    typedef concepts::Digraph GR;
481 481
    checkConcept< McfClassConcept<GR, int, int>,
482 482
                  CostScaling<GR> >();
483 483
    checkConcept< McfClassConcept<GR, double, double>,
484 484
                  CostScaling<GR, double> >();
485 485
    checkConcept< McfClassConcept<GR, int, double>,
486 486
                  CostScaling<GR, int, double> >();
487 487
    typedef CostScaling<GR>::
488 488
      SetLargeCost<double>::Create COS;
489 489
    checkConcept< McfClassConcept<GR, int, int>, COS >();
490 490
  }
491 491

	
492 492
  // Check the interface of CycleCanceling
493 493
  {
494 494
    typedef concepts::Digraph GR;
495 495
    checkConcept< McfClassConcept<GR, int, int>,
496 496
                  CycleCanceling<GR> >();
497 497
    checkConcept< McfClassConcept<GR, double, double>,
498 498
                  CycleCanceling<GR, double> >();
499 499
    checkConcept< McfClassConcept<GR, int, double>,
500 500
                  CycleCanceling<GR, int, double> >();
501 501
  }
502 502

	
503 503
  // Test NetworkSimplex
504 504
  { 
505 505
    typedef NetworkSimplex<Digraph> MCF;
506 506
    runMcfGeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE", true);
507 507
    runMcfLeqTests<MCF>(MCF::FIRST_ELIGIBLE, "NS-FE");
508 508
    runMcfGeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE", true);
509 509
    runMcfLeqTests<MCF>(MCF::BEST_ELIGIBLE,  "NS-BE");
510 510
    runMcfGeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS", true);
511 511
    runMcfLeqTests<MCF>(MCF::BLOCK_SEARCH,   "NS-BS");
512 512
    runMcfGeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL", true);
513 513
    runMcfLeqTests<MCF>(MCF::CANDIDATE_LIST, "NS-CL");
514 514
    runMcfGeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL", true);
515 515
    runMcfLeqTests<MCF>(MCF::ALTERING_LIST,  "NS-AL");
516 516
  }
517 517
  
518 518
  // Test CapacityScaling
519 519
  {
520 520
    typedef CapacityScaling<Digraph> MCF;
521 521
    runMcfGeqTests<MCF>(0, "SSP");
522 522
    runMcfGeqTests<MCF>(2, "CAS");
523 523
  }
524 524

	
525 525
  // Test CostScaling
526 526
  {
527 527
    typedef CostScaling<Digraph> MCF;
528 528
    runMcfGeqTests<MCF>(MCF::PUSH, "COS-PR");
529 529
    runMcfGeqTests<MCF>(MCF::AUGMENT, "COS-AR");
530 530
    runMcfGeqTests<MCF>(MCF::PARTIAL_AUGMENT, "COS-PAR");
531 531
  }
532 532

	
533 533
  // Test CycleCanceling
534 534
  {
535 535
    typedef CycleCanceling<Digraph> MCF;
536 536
    runMcfGeqTests<MCF>(MCF::SIMPLE_CYCLE_CANCELING, "SCC");
537 537
    runMcfGeqTests<MCF>(MCF::MINIMUM_MEAN_CYCLE_CANCELING, "MMCC");
538 538
    runMcfGeqTests<MCF>(MCF::CANCEL_AND_TIGHTEN, "CAT");
539 539
  }
540 540

	
541 541
  return 0;
542 542
}
0 comments (0 inline)