Google OR-Tools: ortools/graph/hamiltonian

Google OR-Tools: ortools/graph/hamiltonian_path.h Source File

519

520 private:

521

522 template <typename T,

523 bool = true >

524

525

526 struct SaturatedArithmetic {

527 static T Add(T a, T b) { return a + b; }

528 static T Sub(T a, T b) { return a - b; }

529 };

530 template <bool Dummy>

531 struct SaturatedArithmetic<int64_t, Dummy> {

532 static int64_t Add(int64_t a, int64_t b) { return CapAdd(a, b); }

533 static int64_t Sub(int64_t a, int64_t b) { return CapSub(a, b); }

534 };

535

536 template <bool Dummy>

537 struct SaturatedArithmetic<int32_t, Dummy> {

538 static int32_t Add(int32_t a, int32_t b) {

539 const int64_t a64 = a;

540 const int64_t b64 = b;

541 const int64_t min_int32 = std::numeric_limits<int32_t>::min();

542 const int64_t max_int32 = std::numeric_limits<int32_t>::max();

543 return static_cast<int32_t>(

544 std::max(min_int32, std::min(max_int32, a64 + b64)));

545 }

546 static int32_t Sub(int32_t a, int32_t b) {

547 const int64_t a64 = a;

548 const int64_t b64 = b;

549 const int64_t min_int32 = std::numeric_limits<int32_t>::min();

550 const int64_t max_int32 = std::numeric_limits<int32_t>::max();

551 return static_cast<int32_t>(

552 std::max(min_int32, std::min(max_int32, a64 - b64)));

553 }

554 };

555

556 template <typename T>

557 using Saturated = SaturatedArithmetic<T>;

558

559

560 CostType Cost(int i, int j) { return cost_(i, j); }

561

562

563 void Solve();

564

565

566 std::vector<int> ComputePath(CostType cost, NodeSet set, int end);

567

568

569 bool PathIsValid(absl::Span<const int> path, CostType cost);

570

571

572 MatrixOrFunction<CostType, CostFunction, true> cost_;

573

574

575 int num_nodes_;

576

577

578 CostType tsp_cost_;

579

580

581 std::vector<CostType> hamiltonian_costs_;

582

583 bool robust_;

584 bool triangle_inequality_ok_;

585 bool robustness_checked_;

586 bool triangle_inequality_checked_;

587 bool solved_;

588 std::vector<int> tsp_path_;

589

590

591

592 std::vector<std::vector<int>> hamiltonian_paths_;

593

594

595

596

597 int best_hamiltonian_path_end_node_;

598

599 LatticeMemoryManager<NodeSet, CostType> mem_;

600};

601

602

603template <typename CostType, typename CostFunction>

605 int num_nodes, CostFunction cost) {

607 std::move(cost));

608}

609

610template <typename CostType, typename CostFunction>

614

615template <typename CostType, typename CostFunction>

617 int num_nodes, CostFunction cost)

618 : cost_(std::move(cost)),

619 num_nodes_(num_nodes),

620 tsp_cost_(0),

621 hamiltonian_costs_(0),

622 robust_(true),

623 triangle_inequality_ok_(true),

624 robustness_checked_(false),

625 triangle_inequality_checked_(false),

626 solved_(false) {

628 CHECK(cost_.Check());

629}

630

631template <typename CostType, typename CostFunction>

633 CostFunction cost) {

635}

636

637template <typename CostType, typename CostFunction>

639 int num_nodes, CostFunction cost) {

640 robustness_checked_ = false;

641 triangle_inequality_checked_ = false;

642 solved_ = false;

643 cost_.Reset(cost);

644 num_nodes_ = num_nodes;

646 CHECK(cost_.Check());

647}

648

649template <typename CostType, typename CostFunction>

650void HamiltonianPathSolver<CostType, CostFunction>::Solve() {

651 if (solved_) return;

652 if (num_nodes_ == 0) {

653 tsp_cost_ = 0;

654 tsp_path_ = {0};

655 hamiltonian_paths_.resize(1);

656 hamiltonian_costs_.resize(1);

657 best_hamiltonian_path_end_node_ = 0;

658 hamiltonian_costs_[0] = 0;

659 hamiltonian_paths_[0] = {0};

660 return;

661 }

662 mem_.Init(num_nodes_);

663

664

665 for (int dest = 0; dest < num_nodes_; ++dest) {

666 DCHECK_EQ(dest, mem_.BaseOffset(1, NodeSet::Singleton(dest)));

667 mem_.SetValueAtOffset(dest, Cost(0, dest));

668 }

669

670

671

672 for (int card = 2; card <= num_nodes_; ++card) {

673

674 for (NodeSet set :

676

677

678 const uint64_t set_offset = mem_.BaseOffset(card, set);

679

680

681

682 uint64_t subset_offset =

683 mem_.BaseOffset(card - 1, set.RemoveSmallestElement());

684 int prev_dest = set.SmallestElement();

685 int dest_rank = 0;

686 for (int dest : set) {

687 CostType min_cost = std::numeric_limits<CostType>::max();

688 const NodeSet subset = set.RemoveElement(dest);

689

690

691

692 subset_offset += mem_.OffsetDelta(card - 1, prev_dest, dest, dest_rank);

693 int src_rank = 0;

694 for (int src : subset) {

695 min_cost = std::min(

696 min_cost, Saturated<CostType>::Add(

697 Cost(src, dest),

698 mem_.ValueAtOffset(subset_offset + src_rank)));

699 ++src_rank;

700 }

701 prev_dest = dest;

702 mem_.SetValueAtOffset(set_offset + dest_rank, min_cost);

703 ++dest_rank;

704 }

705 }

706 }

707

708 const NodeSet full_set = NodeSet::FullSet(num_nodes_);

709

710

711

712 tsp_cost_ = mem_.Value(full_set, 0);

713 tsp_path_ = ComputePath(tsp_cost_, full_set, 0);

714

715 hamiltonian_paths_.resize(num_nodes_);

716 hamiltonian_costs_.resize(num_nodes_);

717

718

719

720 CostType min_hamiltonian_cost = std::numeric_limits<CostType>::max();

721 const NodeSet hamiltonian_set = full_set.RemoveElement(0);

722 for (int end_node : hamiltonian_set) {

723 const CostType cost = mem_.Value(hamiltonian_set, end_node);

724 hamiltonian_costs_[end_node] = cost;

725 if (cost <= min_hamiltonian_cost) {

726 min_hamiltonian_cost = cost;

727 best_hamiltonian_path_end_node_ = end_node;

728 }

729 DCHECK_LE(tsp_cost_, Saturated<CostType>::Add(cost, Cost(end_node, 0)));

730

731 hamiltonian_paths_[end_node] =

732 ComputePath(hamiltonian_costs_[end_node], hamiltonian_set, end_node);

733 }

734

735 solved_ = true;

736}

737

738template <typename CostType, typename CostFunction>

739std::vector<int> HamiltonianPathSolver<CostType, CostFunction>::ComputePath(

740 CostType cost, NodeSet set, int end_node) {

741 DCHECK(set.Contains(end_node));

742 const int path_size = set.Cardinality() + 1;

743 std::vector<int> path(path_size, 0);

744 NodeSet subset = set.RemoveElement(end_node);

745 path[path_size - 1] = end_node;

746 int dest = end_node;

747 CostType current_cost = cost;

748 for (int rank = path_size - 2; rank >= 0; --rank) {

749 for (int src : subset) {

750 const CostType partial_cost = mem_.Value(subset, src);

751 const CostType incumbent_cost =

752 Saturated<CostType>::Add(partial_cost, Cost(src, dest));

753

754

755 if (std::abs(Saturated<CostType>::Sub(current_cost, incumbent_cost)) <=

756 std::numeric_limits<CostType>::epsilon() * current_cost) {

757 subset = subset.RemoveElement(src);

758 current_cost = partial_cost;

759 path[rank] = src;

760 dest = src;

761 break;

762 }

763 }

764 }

765 DCHECK_EQ(0, subset.value());

766 DCHECK(PathIsValid(path, cost));

767 return path;

768}

769

770template <typename CostType, typename CostFunction>

771bool HamiltonianPathSolver<CostType, CostFunction>::PathIsValid(

772 absl::Span<const int> path, CostType cost) {

773 NodeSet coverage(0);

774 for (int node : path) {

775 coverage = coverage.AddElement(node);

776 }

777 DCHECK_EQ(NodeSet::FullSet(num_nodes_).value(), coverage.value());

778 CostType check_cost = 0;

779 for (int i = 0; i < path.size() - 1; ++i) {

780 check_cost =

781 Saturated<CostType>::Add(check_cost, Cost(path[i], path[i + 1]));

782 }

783 DCHECK_LE(std::abs(Saturated<CostType>::Sub(cost, check_cost)),

784 std::numeric_limits<CostType>::epsilon() * cost)

785 << "cost = " << cost << " check_cost = " << check_cost;

786 return true;

787}

788

789template <typename CostType, typename CostFunction>

791 if (std::numeric_limits<CostType>::is_integer) return true;

792 if (robustness_checked_) return robust_;

793 CostType min_cost = std::numeric_limits<CostType>::max();

794 CostType max_cost = std::numeric_limits<CostType>::min();

795

796

797 for (int i = 0; i < num_nodes_; ++i) {

798 for (int j = 0; j < num_nodes_; ++j) {

799 if (i == j) continue;

800 min_cost = std::min(min_cost, Cost(i, j));

801 max_cost = std::max(max_cost, Cost(i, j));

802 }

803 }

804

805

806 robust_ =

807 min_cost >= 0 && min_cost > num_nodes_ * max_cost *

808 std::numeric_limits<CostType>::epsilon();

809 robustness_checked_ = true;

810 return robust_;

811}

812

813template <typename CostType, typename CostFunction>

816 if (triangle_inequality_checked_) return triangle_inequality_ok_;

817 triangle_inequality_ok_ = true;

818 triangle_inequality_checked_ = true;

819 for (int k = 0; k < num_nodes_; ++k) {

820 for (int i = 0; i < num_nodes_; ++i) {

821 for (int j = 0; j < num_nodes_; ++j) {

822 const CostType detour_cost =

823 Saturated<CostType>::Add(Cost(i, k), Cost(k, j));

824 if (detour_cost < Cost(i, j)) {

825 triangle_inequality_ok_ = false;

826 return triangle_inequality_ok_;

827 }

828 }

829 }

830 }

831 return triangle_inequality_ok_;

832}

837 Solve();

838 return best_hamiltonian_path_end_node_;

839}

843 int end_node) {

844 Solve();

845 return hamiltonian_costs_[end_node];

846}

850 int end_node) {

851 Solve();

852 return hamiltonian_paths_[end_node];

853}

857 std::vector<PathNodeIndex>* path) {

858 *path = HamiltonianPath(best_hamiltonian_path_end_node_);

859}

864 Solve();

865 return tsp_cost_;

866}

871 Solve();

872 return tsp_path_;

873}

877 std::vector<PathNodeIndex>* path) {

879}

883

884

885

886

887

888

889

890

891

892

893

894

895

896

897 public:

900

903

904

906

907

908

909

910 private:

911

912 CostType Cost(int i, int j) { return cost_(i, j); }

913

914

915 void Solve(int end_node);

916

917

918 CostType ComputeFutureLowerBound(NodeSet current_set, int last_visited);

919

920

922

923

924 int num_nodes_;

925

926

927 CostType tsp_cost_;

928

929

930 bool solved_;

931

932

934};

943 int num_nodes, CostFunction cost)

944 : cost_(std::move(cost)),

945 num_nodes_(num_nodes),

946 tsp_cost_(0),

947 solved_(false) {}

1012 int end_node) {

1013 Solve(end_node);

1014 return tsp_cost_;

1015}