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 96 line context
... ...
@@ -269,159 +269,97 @@
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;
... ...
@@ -466,171 +404,262 @@
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]];
Show white space 96 line context
... ...
@@ -288,163 +288,97 @@
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;
... ...
@@ -489,179 +423,267 @@
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]];
Show white space 96 line context
... ...
@@ -205,161 +205,97 @@
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;
... ...
@@ -404,177 +340,270 @@
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]];
Show white space 96 line context
... ...
@@ -149,96 +149,97 @@
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;
... ...
@@ -588,158 +589,108 @@
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
    }
... ...
@@ -797,167 +748,248 @@
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
    }
Show white space 96 line context
... ...
@@ -112,97 +112,97 @@
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,
... ...
@@ -301,131 +301,131 @@
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)
0 comments (0 inline)