gravatar
alpar (Alpar Juttner)
alpar@cs.elte.hu
Merge
0 1 0
merge default
1 file changed with 71 insertions and 64 deletions:
↑ Collapse diff ↑
Ignore white space 96 line context
... ...
@@ -317,415 +317,422 @@
317 317
          }
318 318
        }
319 319
        return min < 0;
320 320
      }
321 321

	
322 322
    }; //class BestEligiblePivotRule
323 323

	
324 324

	
325 325
    // Implementation of the Block Search pivot rule
326 326
    class BlockSearchPivotRule
327 327
    {
328 328
    private:
329 329

	
330 330
      // References to the NetworkSimplex class
331 331
      const IntVector  &_source;
332 332
      const IntVector  &_target;
333 333
      const CostVector &_cost;
334 334
      const IntVector  &_state;
335 335
      const CostVector &_pi;
336 336
      int &_in_arc;
337 337
      int _search_arc_num;
338 338

	
339 339
      // Pivot rule data
340 340
      int _block_size;
341 341
      int _next_arc;
342 342

	
343 343
    public:
344 344

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

	
356 356
        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
357 357
                                    std::sqrt(double(_search_arc_num))),
358 358
                                MIN_BLOCK_SIZE );
359 359
      }
360 360

	
361 361
      // Find next entering arc
362 362
      bool findEnteringArc() {
363 363
        Cost c, min = 0;
364 364
        int cnt = _block_size;
365
        int e, min_arc = _next_arc;
365
        int e;
366 366
        for (e = _next_arc; e < _search_arc_num; ++e) {
367 367
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
368 368
          if (c < min) {
369 369
            min = c;
370
            min_arc = e;
370
            _in_arc = e;
371 371
          }
372 372
          if (--cnt == 0) {
373
            if (min < 0) break;
373
            if (min < 0) goto search_end;
374 374
            cnt = _block_size;
375 375
          }
376 376
        }
377
        if (min == 0 || cnt > 0) {
378
          for (e = 0; e < _next_arc; ++e) {
379
            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
380
            if (c < min) {
381
              min = c;
382
              min_arc = e;
383
            }
384
            if (--cnt == 0) {
385
              if (min < 0) break;
386
              cnt = _block_size;
387
            }
377
        for (e = 0; e < _next_arc; ++e) {
378
          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
379
          if (c < min) {
380
            min = c;
381
            _in_arc = e;
382
          }
383
          if (--cnt == 0) {
384
            if (min < 0) goto search_end;
385
            cnt = _block_size;
388 386
          }
389 387
        }
390 388
        if (min >= 0) return false;
391
        _in_arc = min_arc;
389

	
390
      search_end:
392 391
        _next_arc = e;
393 392
        return true;
394 393
      }
395 394

	
396 395
    }; //class BlockSearchPivotRule
397 396

	
398 397

	
399 398
    // Implementation of the Candidate List pivot rule
400 399
    class CandidateListPivotRule
401 400
    {
402 401
    private:
403 402

	
404 403
      // References to the NetworkSimplex class
405 404
      const IntVector  &_source;
406 405
      const IntVector  &_target;
407 406
      const CostVector &_cost;
408 407
      const IntVector  &_state;
409 408
      const CostVector &_pi;
410 409
      int &_in_arc;
411 410
      int _search_arc_num;
412 411

	
413 412
      // Pivot rule data
414 413
      IntVector _candidates;
415 414
      int _list_length, _minor_limit;
416 415
      int _curr_length, _minor_count;
417 416
      int _next_arc;
418 417

	
419 418
    public:
420 419

	
421 420
      /// Constructor
422 421
      CandidateListPivotRule(NetworkSimplex &ns) :
423 422
        _source(ns._source), _target(ns._target),
424 423
        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
425 424
        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
426 425
        _next_arc(0)
427 426
      {
428 427
        // The main parameters of the pivot rule
429
        const double LIST_LENGTH_FACTOR = 1.0;
428
        const double LIST_LENGTH_FACTOR = 0.25;
430 429
        const int MIN_LIST_LENGTH = 10;
431 430
        const double MINOR_LIMIT_FACTOR = 0.1;
432 431
        const int MIN_MINOR_LIMIT = 3;
433 432

	
434 433
        _list_length = std::max( int(LIST_LENGTH_FACTOR *
435 434
                                     std::sqrt(double(_search_arc_num))),
436 435
                                 MIN_LIST_LENGTH );
437 436
        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
438 437
                                 MIN_MINOR_LIMIT );
439 438
        _curr_length = _minor_count = 0;
440 439
        _candidates.resize(_list_length);
441 440
      }
442 441

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

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

	
503 498
    }; //class CandidateListPivotRule
504 499

	
505 500

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

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

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

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

	
538 533
      SortFunc _sort_func;
539 534

	
540 535
    public:
541 536

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

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

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

	
577 572
        // Extend the list
578 573
        int cnt = _block_size;
579
        int last_arc = 0;
580 574
        int limit = _head_length;
581 575

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

	
613 604
        // Make heap of the candidate list (approximating a partial sort)
614 605
        make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
615 606
                   _sort_func );
616 607

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

	
625 617
    }; //class AlteringListPivotRule
626 618

	
627 619
  public:
628 620

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

	
652 648
      _source.resize(max_arc_num);
653 649
      _target.resize(max_arc_num);
654 650

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

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

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

	
690 697
    /// \name Parameters
691 698
    /// The parameters of the algorithm can be specified using these
692 699
    /// functions.
693 700

	
694 701
    /// @{
695 702

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

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