20 #include "absl/strings/str_format.h"
30 "Divide factor for epsilon at each refine step.");
31 ABSL_FLAG(
bool, min_cost_flow_check_feasibility,
true,
32 "Check that the graph has enough capacity to send all supplies "
33 "and serve all demands. Also check that the sum of supplies "
34 "is equal to the sum of demands.");
36 "Check that the sum of supplies is equal to the sum of demands.");
38 "Check that the magnitude of the costs will not exceed the "
39 "precision of the machine when scaled (multiplied) by the number "
42 "Check that the result is valid.");
46 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
52 residual_arc_capacity_(),
53 first_admissible_arc_(),
56 alpha_(
absl::GetFlag(FLAGS_min_cost_flow_alpha)),
57 cost_scaling_factor_(1),
58 scaled_arc_unit_cost_(),
61 initial_node_excess_(),
62 feasible_node_excess_(),
63 stats_(
"MinCostFlow"),
64 feasibility_checked_(false),
65 use_price_update_(false),
66 check_feasibility_(
absl::GetFlag(FLAGS_min_cost_flow_check_feasibility)) {
68 if (max_num_nodes > 0) {
69 node_excess_.
Reserve(0, max_num_nodes - 1);
71 node_potential_.
Reserve(0, max_num_nodes - 1);
73 first_admissible_arc_.
Reserve(0, max_num_nodes - 1);
74 first_admissible_arc_.
SetAll(Graph::kNilArc);
75 initial_node_excess_.
Reserve(0, max_num_nodes - 1);
76 initial_node_excess_.
SetAll(0);
77 feasible_node_excess_.
Reserve(0, max_num_nodes - 1);
78 feasible_node_excess_.
SetAll(0);
81 if (max_num_arcs > 0) {
82 residual_arc_capacity_.Reserve(-max_num_arcs, max_num_arcs - 1);
83 residual_arc_capacity_.SetAll(0);
84 scaled_arc_unit_cost_.Reserve(-max_num_arcs, max_num_arcs - 1);
85 scaled_arc_unit_cost_.SetAll(0);
89 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
92 DCHECK(graph_->IsNodeValid(node));
93 node_excess_.Set(node, supply);
94 initial_node_excess_.Set(node, supply);
96 feasibility_checked_ =
false;
99 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
101 ArcIndex arc, ArcScaledCostType unit_cost) {
103 scaled_arc_unit_cost_.Set(arc, unit_cost);
104 scaled_arc_unit_cost_.Set(Opposite(arc), -scaled_arc_unit_cost_[arc]);
105 status_ = NOT_SOLVED;
106 feasibility_checked_ =
false;
109 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
111 ArcIndex arc, ArcFlowType new_capacity) {
114 const FlowQuantity free_capacity = residual_arc_capacity_[arc];
115 const FlowQuantity capacity_delta = new_capacity - Capacity(arc);
116 if (capacity_delta == 0) {
119 status_ = NOT_SOLVED;
120 feasibility_checked_ =
false;
121 const FlowQuantity new_availability = free_capacity + capacity_delta;
122 if (new_availability >= 0) {
128 DCHECK((capacity_delta > 0) ||
129 (capacity_delta < 0 && new_availability >= 0));
130 residual_arc_capacity_.Set(arc, new_availability);
131 DCHECK_LE(0, residual_arc_capacity_[arc]);
135 const FlowQuantity flow = residual_arc_capacity_[Opposite(arc)];
137 residual_arc_capacity_.Set(arc, 0);
138 residual_arc_capacity_.Set(Opposite(arc), new_capacity);
140 node_excess_.Set(
tail, node_excess_[
tail] + flow_excess);
142 node_excess_.Set(
head, node_excess_[
head] - flow_excess);
143 DCHECK_LE(0, residual_arc_capacity_[arc]);
144 DCHECK_LE(0, residual_arc_capacity_[Opposite(arc)]);
148 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
150 ArcIndex arc, ArcFlowType new_flow) {
154 residual_arc_capacity_.Set(Opposite(arc), new_flow);
155 residual_arc_capacity_.Set(arc,
capacity - new_flow);
156 status_ = NOT_SOLVED;
157 feasibility_checked_ =
false;
160 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
162 ArcScaledCostType>::CheckInputConsistency()
const {
166 for (
ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
171 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
173 total_supply += excess;
175 total_flow += excess;
177 max_capacity + total_flow) {
178 LOG(DFATAL) <<
"Input consistency error: max capacity + flow exceed "
184 if (total_supply != 0) {
185 LOG(DFATAL) <<
"Input consistency error: unbalanced problem";
191 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
192 bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::CheckResult()
194 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
195 if (node_excess_[node] != 0) {
196 LOG(DFATAL) <<
"node_excess_[" << node <<
"] != 0";
199 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
203 if (residual_arc_capacity_[arc] < 0) {
204 LOG(DFATAL) <<
"residual_arc_capacity_[" << arc <<
"] < 0";
207 if (residual_arc_capacity_[arc] > 0 && ReducedCost(arc) < -epsilon_) {
208 LOG(DFATAL) <<
"residual_arc_capacity_[" << arc
209 <<
"] > 0 && ReducedCost(" << arc <<
") < " << -epsilon_
210 <<
". (epsilon_ = " << epsilon_ <<
").";
214 LOG(DFATAL) << DebugString(
"CheckResult ", arc);
222 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
223 bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::CheckCostRange()
228 for (
ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
230 max_cost_magnitude =
std::max(max_cost_magnitude, cost_magnitude);
231 if (cost_magnitude != 0.0) {
232 min_cost_magnitude =
std::min(min_cost_magnitude, cost_magnitude);
235 VLOG(3) <<
"Min cost magnitude = " << min_cost_magnitude
236 <<
", Max cost magnitude = " << max_cost_magnitude;
237 #if !defined(_MSC_VER)
239 log(max_cost_magnitude + 1) + log(graph_->num_nodes() + 1)) {
240 LOG(DFATAL) <<
"Maximum cost magnitude " << max_cost_magnitude <<
" is too "
241 <<
"high for the number of nodes. Try changing the data.";
248 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
249 bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::
250 CheckRelabelPrecondition(
NodeIndex node)
const {
258 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
261 DCHECK(!IsAdmissible(arc)) << DebugString(
"CheckRelabelPrecondition:", arc);
266 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
268 GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::DebugString(
275 const CostValue reduced_cost = scaled_arc_unit_cost_[arc] +
276 node_potential_[
tail] - node_potential_[
head];
277 return absl::StrFormat(
278 "%s Arc %d, from %d to %d, "
279 "Capacity = %d, Residual capacity = %d, "
280 "Flow = residual capacity for reverse arc = %d, "
281 "Height(tail) = %d, Height(head) = %d, "
282 "Excess(tail) = %d, Excess(head) = %d, "
283 "Cost = %d, Reduced cost = %d, ",
285 static_cast<FlowQuantity>(residual_arc_capacity_[arc]), Flow(arc),
286 node_potential_[
tail], node_potential_[
head], node_excess_[
tail],
287 node_excess_[
head],
static_cast<CostValue>(scaled_arc_unit_cost_[arc]),
291 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
294 std::vector<NodeIndex>*
const infeasible_demand_node) {
307 feasibility_checked_ =
false;
309 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
310 if (initial_node_excess_[node] != 0) {
314 const NodeIndex num_nodes_in_max_flow = graph_->num_nodes() + 2;
315 const ArcIndex num_arcs_in_max_flow = graph_->num_arcs() + num_extra_arcs;
316 const NodeIndex source = num_nodes_in_max_flow - 2;
317 const NodeIndex sink = num_nodes_in_max_flow - 1;
318 StarGraph checker_graph(num_nodes_in_max_flow, num_arcs_in_max_flow);
319 MaxFlow checker(&checker_graph, source, sink);
323 for (
ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
325 checker_graph.
AddArc(graph_->Tail(arc), graph_->Head(arc));
332 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
337 total_supply += supply;
338 }
else if (supply < 0) {
341 total_demand -= supply;
344 if (total_supply != total_demand) {
345 LOG(DFATAL) <<
"total_supply(" << total_supply <<
") != total_demand("
346 << total_demand <<
").";
349 if (!checker.
Solve()) {
350 LOG(DFATAL) <<
"Max flow could not be computed.";
354 feasible_node_excess_.SetAll(0);
355 for (StarGraph::OutgoingArcIterator it(checker_graph, source); it.Ok();
360 feasible_node_excess_.Set(node, flow);
361 if (infeasible_supply_node !=
nullptr) {
362 infeasible_supply_node->push_back(node);
365 for (StarGraph::IncomingArcIterator it(checker_graph, sink); it.Ok();
370 feasible_node_excess_.Set(node, -flow);
371 if (infeasible_demand_node !=
nullptr) {
372 infeasible_demand_node->push_back(node);
375 feasibility_checked_ =
true;
376 return optimal_max_flow == total_supply;
379 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
381 if (!feasibility_checked_) {
384 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
385 const FlowQuantity excess = feasible_node_excess_[node];
386 node_excess_.Set(node, excess);
387 initial_node_excess_.Set(node, excess);
392 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
395 if (IsArcDirect(arc)) {
396 return residual_arc_capacity_[Opposite(arc)];
398 return -residual_arc_capacity_[arc];
403 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
407 if (IsArcDirect(arc)) {
408 return residual_arc_capacity_[arc] + residual_arc_capacity_[Opposite(arc)];
414 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
419 return scaled_arc_unit_cost_[arc];
422 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
425 DCHECK(graph_->IsNodeValid(node));
426 return node_excess_[node];
429 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
433 return initial_node_excess_[node];
436 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
440 return feasible_node_excess_[node];
443 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
446 return FastIsAdmissible(arc, node_potential_[Tail(arc)]);
449 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
450 bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::
452 DCHECK_EQ(node_potential_[Tail(arc)], tail_potential);
453 return residual_arc_capacity_[arc] > 0 &&
454 FastReducedCost(arc, tail_potential) < 0;
457 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
458 bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsActive(
460 return node_excess_[node] > 0;
463 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
465 GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::ReducedCost(
467 return FastReducedCost(arc, node_potential_[Tail(arc)]);
470 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
472 GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::FastReducedCost(
474 DCHECK_EQ(node_potential_[Tail(arc)], tail_potential);
475 DCHECK(graph_->IsNodeValid(Tail(arc)));
476 DCHECK(graph_->IsNodeValid(Head(arc)));
477 DCHECK_LE(node_potential_[Tail(arc)], 0) << DebugString(
"ReducedCost:", arc);
478 DCHECK_LE(node_potential_[Head(arc)], 0) << DebugString(
"ReducedCost:", arc);
479 return scaled_arc_unit_cost_[arc] + tail_potential -
480 node_potential_[Head(arc)];
483 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
485 GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::
486 GetFirstOutgoingOrOppositeIncomingArc(
NodeIndex node)
const {
487 OutgoingOrOppositeIncomingArcIterator arc_it(*graph_, node);
488 return arc_it.Index();
491 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
493 status_ = NOT_SOLVED;
494 if (absl::GetFlag(FLAGS_min_cost_flow_check_balance) &&
495 !CheckInputConsistency()) {
496 status_ = UNBALANCED;
499 if (absl::GetFlag(FLAGS_min_cost_flow_check_costs) && !CheckCostRange()) {
500 status_ = BAD_COST_RANGE;
503 if (check_feasibility_ && !CheckFeasibility(
nullptr,
nullptr)) {
507 node_potential_.SetAll(0);
508 ResetFirstAdmissibleArcs();
511 if (absl::GetFlag(FLAGS_min_cost_flow_check_result) && !CheckResult()) {
512 status_ = BAD_RESULT;
518 LOG(DFATAL) <<
"Status != OPTIMAL";
519 total_flow_cost_ = 0;
522 total_flow_cost_ = 0;
523 for (
ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
524 const FlowQuantity flow_on_arc = residual_arc_capacity_[Opposite(arc)];
525 total_flow_cost_ += scaled_arc_unit_cost_[arc] * flow_on_arc;
532 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
534 ArcScaledCostType>::ResetFirstAdmissibleArcs() {
535 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
536 first_admissible_arc_.Set(node,
537 GetFirstOutgoingOrOppositeIncomingArc(node));
541 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
542 void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::ScaleCosts() {
544 cost_scaling_factor_ = graph_->num_nodes() + 1;
546 VLOG(3) <<
"Number of nodes in the graph = " << graph_->num_nodes();
547 VLOG(3) <<
"Number of arcs in the graph = " << graph_->num_arcs();
548 for (
ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
549 const CostValue cost = scaled_arc_unit_cost_[arc] * cost_scaling_factor_;
550 scaled_arc_unit_cost_.Set(arc,
cost);
551 scaled_arc_unit_cost_.Set(Opposite(arc), -
cost);
554 VLOG(3) <<
"Initial epsilon = " << epsilon_;
555 VLOG(3) <<
"Cost scaling factor = " << cost_scaling_factor_;
558 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
559 void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::UnscaleCosts() {
561 for (
ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
562 const CostValue cost = scaled_arc_unit_cost_[arc] / cost_scaling_factor_;
563 scaled_arc_unit_cost_.Set(arc,
cost);
564 scaled_arc_unit_cost_.Set(Opposite(arc), -
cost);
566 cost_scaling_factor_ = 1;
569 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
570 void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Optimize() {
572 num_relabels_since_last_price_update_ = 0;
575 epsilon_ =
std::max(epsilon_ / alpha_, kEpsilonMin);
576 VLOG(3) <<
"Epsilon changed to: " << epsilon_;
578 }
while (epsilon_ != 1LL && status_ != INFEASIBLE);
579 if (status_ == NOT_SOLVED) {
584 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
585 void GenericMinCostFlow<
Graph, ArcFlowType,
586 ArcScaledCostType>::SaturateAdmissibleArcs() {
588 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
589 const CostValue tail_potential = node_potential_[node];
590 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node,
591 first_admissible_arc_[node]);
592 it.Ok(); it.Next()) {
594 if (FastIsAdmissible(arc, tail_potential)) {
595 FastPushFlow(residual_arc_capacity_[arc], arc, node);
605 first_admissible_arc_[node] = Graph::kNilArc;
609 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
610 void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::PushFlow(
613 FastPushFlow(flow, arc, Tail(arc));
616 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
617 void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::FastPushFlow(
621 DCHECK_GT(residual_arc_capacity_[arc], 0);
622 DCHECK_LE(flow, residual_arc_capacity_[arc]);
624 residual_arc_capacity_.Set(arc, residual_arc_capacity_[arc] - flow);
626 const ArcIndex opposite = Opposite(arc);
627 residual_arc_capacity_.Set(opposite, residual_arc_capacity_[opposite] + flow);
629 node_excess_.Set(
tail, node_excess_[
tail] - flow);
631 node_excess_.Set(
head, node_excess_[
head] + flow);
634 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
635 void GenericMinCostFlow<
Graph, ArcFlowType,
636 ArcScaledCostType>::InitializeActiveNodeStack() {
638 DCHECK(active_nodes_.empty());
639 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
640 if (IsActive(node)) {
641 active_nodes_.push(node);
646 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
647 void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::UpdatePrices() {
666 const NodeIndex num_nodes = graph_->num_nodes();
667 std::vector<NodeIndex> bfs_queue;
668 std::vector<bool> node_in_queue(num_nodes,
false);
672 std::vector<CostValue> min_non_admissible_potential(num_nodes, kMinCostValue);
673 std::vector<NodeIndex> nodes_to_process;
679 for (
NodeIndex node = 0; node < num_nodes; ++node) {
680 if (node_excess_[node] < 0) {
681 bfs_queue.push_back(node);
682 node_in_queue[node] =
true;
685 remaining_excess -= node_excess_[node];
696 while (remaining_excess > 0) {
701 for (; queue_index < bfs_queue.size(); ++queue_index) {
703 const NodeIndex node = bfs_queue[queue_index];
704 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
707 if (node_in_queue[
head])
continue;
708 const ArcIndex opposite_arc = Opposite(it.Index());
709 if (residual_arc_capacity_[opposite_arc] > 0) {
710 node_potential_[
head] += potential_delta;
711 if (ReducedCost(opposite_arc) < 0) {
712 DCHECK(IsAdmissible(opposite_arc));
717 remaining_excess -= node_excess_[
head];
718 if (remaining_excess == 0) {
719 node_potential_[
head] -= potential_delta;
722 bfs_queue.push_back(
head);
723 node_in_queue[
head] =
true;
724 if (potential_delta < 0) {
725 first_admissible_arc_[
head] =
726 GetFirstOutgoingOrOppositeIncomingArc(
head);
731 node_potential_[
head] -= potential_delta;
732 if (min_non_admissible_potential[
head] == kMinCostValue) {
733 nodes_to_process.push_back(
head);
736 min_non_admissible_potential[
head],
737 node_potential_[node] - scaled_arc_unit_cost_[opposite_arc]);
741 if (remaining_excess == 0)
break;
743 if (remaining_excess == 0)
break;
747 CostValue max_potential_diff = kMinCostValue;
748 for (
int i = 0; i < nodes_to_process.size(); ++i) {
749 const NodeIndex node = nodes_to_process[i];
750 if (node_in_queue[node])
continue;
753 min_non_admissible_potential[node] - node_potential_[node]);
754 if (max_potential_diff == potential_delta)
break;
756 DCHECK_LE(max_potential_diff, potential_delta);
757 potential_delta = max_potential_diff - epsilon_;
766 for (
int i = 0; i < nodes_to_process.size(); ++i) {
767 const NodeIndex node = nodes_to_process[i];
768 if (node_in_queue[node])
continue;
769 if (node_potential_[node] + potential_delta <
770 min_non_admissible_potential[node]) {
771 node_potential_[node] += potential_delta;
772 first_admissible_arc_[node] =
773 GetFirstOutgoingOrOppositeIncomingArc(node);
774 bfs_queue.push_back(node);
775 node_in_queue[node] =
true;
776 remaining_excess -= node_excess_[node];
781 nodes_to_process[
index] = node;
784 nodes_to_process.resize(
index);
788 if (potential_delta == 0)
return;
789 for (
NodeIndex node = 0; node < num_nodes; ++node) {
790 if (!node_in_queue[node]) {
791 node_potential_[node] += potential_delta;
792 first_admissible_arc_[node] = GetFirstOutgoingOrOppositeIncomingArc(node);
797 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
798 void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Refine() {
800 SaturateAdmissibleArcs();
801 InitializeActiveNodeStack();
803 const NodeIndex num_nodes = graph_->num_nodes();
804 while (status_ != INFEASIBLE && !active_nodes_.empty()) {
806 if (num_relabels_since_last_price_update_ >= num_nodes) {
807 num_relabels_since_last_price_update_ = 0;
808 if (use_price_update_) {
812 const NodeIndex node = active_nodes_.top();
819 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
820 void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Discharge(
827 const CostValue tail_potential = node_potential_[node];
828 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node,
829 first_admissible_arc_[node]);
830 it.Ok(); it.Next()) {
832 if (FastIsAdmissible(arc, tail_potential)) {
834 if (!LookAhead(arc, tail_potential,
head))
continue;
835 const bool head_active_before_push = IsActive(
head);
838 static_cast<FlowQuantity>(residual_arc_capacity_[arc]));
839 FastPushFlow(
delta, arc, node);
840 if (IsActive(
head) && !head_active_before_push) {
841 active_nodes_.push(
head);
843 if (node_excess_[node] == 0) {
845 first_admissible_arc_.Set(node, arc);
851 }
while (status_ != INFEASIBLE);
854 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
855 bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::LookAhead(
859 DCHECK_EQ(node_potential_[Tail(in_arc)], in_tail_potential);
860 if (node_excess_[node] < 0)
return true;
861 const CostValue tail_potential = node_potential_[node];
862 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node,
863 first_admissible_arc_[node]);
864 it.Ok(); it.Next()) {
866 if (FastIsAdmissible(arc, tail_potential)) {
867 first_admissible_arc_.Set(node, arc);
875 return FastIsAdmissible(in_arc, in_tail_potential);
878 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
879 void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Relabel(
882 DCHECK(CheckRelabelPrecondition(node));
883 ++num_relabels_since_last_price_update_;
890 const CostValue guaranteed_new_potential = node_potential_[node] - epsilon_;
898 CostValue min_non_admissible_potential = kMinCostValue;
903 CostValue previous_min_non_admissible_potential = kMinCostValue;
904 ArcIndex first_arc = Graph::kNilArc;
906 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
909 if (residual_arc_capacity_[arc] > 0) {
910 const CostValue min_non_admissible_potential_for_arc =
911 node_potential_[Head(arc)] - scaled_arc_unit_cost_[arc];
912 if (min_non_admissible_potential_for_arc > min_non_admissible_potential) {
913 if (min_non_admissible_potential_for_arc > guaranteed_new_potential) {
917 node_potential_.Set(node, guaranteed_new_potential);
918 first_admissible_arc_.Set(node, arc);
921 previous_min_non_admissible_potential = min_non_admissible_potential;
922 min_non_admissible_potential = min_non_admissible_potential_for_arc;
929 if (min_non_admissible_potential == kMinCostValue) {
930 if (node_excess_[node] != 0) {
934 LOG(
ERROR) <<
"Infeasible problem.";
939 node_potential_.Set(node, guaranteed_new_potential);
940 first_admissible_arc_.Set(node,
941 GetFirstOutgoingOrOppositeIncomingArc(node));
949 const CostValue new_potential = min_non_admissible_potential - epsilon_;
950 node_potential_.Set(node, new_potential);
951 if (previous_min_non_admissible_potential <= new_potential) {
952 first_admissible_arc_.Set(node, first_arc);
955 first_admissible_arc_.Set(node,
956 GetFirstOutgoingOrOppositeIncomingArc(node));
960 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
962 GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Opposite(
967 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
968 bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsArcValid(
973 template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
974 bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsArcDirect(
984 template class GenericMinCostFlow<StarGraph>;
985 template class GenericMinCostFlow<::util::ReverseArcListGraph<>>;
986 template class GenericMinCostFlow<::util::ReverseArcStaticGraph<>>;
987 template class GenericMinCostFlow<::util::ReverseArcMixedGraph<>>;
988 template class GenericMinCostFlow<::util::ReverseArcStaticGraph<uint16, int32>>;
991 template class GenericMinCostFlow<::util::ReverseArcStaticGraph<uint16, int32>,
997 if (reserve_num_nodes > 0) {
998 node_supply_.reserve(reserve_num_nodes);
1000 if (reserve_num_arcs > 0) {
1001 arc_tail_.reserve(reserve_num_arcs);
1002 arc_head_.reserve(reserve_num_arcs);
1003 arc_capacity_.reserve(reserve_num_arcs);
1004 arc_cost_.reserve(reserve_num_arcs);
1005 arc_permutation_.reserve(reserve_num_arcs);
1006 arc_flow_.reserve(reserve_num_arcs);
1011 ResizeNodeVectors(node);
1012 node_supply_[node] = supply;
1020 const ArcIndex arc = arc_tail_.size();
1021 arc_tail_.push_back(
tail);
1022 arc_head_.push_back(
head);
1024 arc_cost_.push_back(unit_cost);
1029 return arc < arc_permutation_.size() ? arc_permutation_[arc] : arc;
1033 SupplyAdjustment adjustment) {
1037 const NodeIndex num_nodes = node_supply_.size();
1038 const ArcIndex num_arcs = arc_capacity_.size();
1039 if (num_nodes == 0)
return OPTIMAL;
1041 int supply_node_count = 0, demand_node_count = 0;
1043 for (
NodeIndex node = 0; node < num_nodes; ++node) {
1044 if (node_supply_[node] > 0) {
1045 ++supply_node_count;
1046 total_supply += node_supply_[node];
1047 }
else if (node_supply_[node] < 0) {
1048 ++demand_node_count;
1049 total_demand -= node_supply_[node];
1052 if (adjustment == DONT_ADJUST && total_supply != total_demand) {
1067 const ArcIndex augmented_num_arcs =
1068 num_arcs + supply_node_count + demand_node_count;
1071 const NodeIndex augmented_num_nodes = num_nodes + 2;
1073 Graph graph(augmented_num_nodes, augmented_num_arcs);
1074 for (
ArcIndex arc = 0; arc < num_arcs; ++arc) {
1075 graph.AddArc(arc_tail_[arc], arc_head_[arc]);
1078 for (
NodeIndex node = 0; node < num_nodes; ++node) {
1079 if (node_supply_[node] > 0) {
1080 graph.AddArc(source, node);
1081 }
else if (node_supply_[node] < 0) {
1082 graph.AddArc(node, sink);
1086 graph.Build(&arc_permutation_);
1089 GenericMaxFlow<Graph> max_flow(&graph, source, sink);
1091 for (arc = 0; arc < num_arcs; ++arc) {
1092 max_flow.SetArcCapacity(PermutedArc(arc), arc_capacity_[arc]);
1094 for (
NodeIndex node = 0; node < num_nodes; ++node) {
1095 if (node_supply_[node] != 0) {
1096 max_flow.SetArcCapacity(PermutedArc(arc), std::abs(node_supply_[node]));
1101 if (!max_flow.Solve()) {
1102 LOG(
ERROR) <<
"Max flow could not be computed.";
1103 switch (max_flow.status()) {
1108 <<
"Max flow failed but claimed to have an optimal solution";
1109 ABSL_FALLTHROUGH_INTENDED;
1114 maximum_flow_ = max_flow.GetOptimalFlow();
1117 if (adjustment == DONT_ADJUST && maximum_flow_ != total_supply) {
1121 GenericMinCostFlow<Graph> min_cost_flow(&graph);
1123 for (arc = 0; arc < num_arcs; ++arc) {
1124 ArcIndex permuted_arc = PermutedArc(arc);
1125 min_cost_flow.SetArcUnitCost(permuted_arc, arc_cost_[arc]);
1126 min_cost_flow.SetArcCapacity(permuted_arc, arc_capacity_[arc]);
1128 for (
NodeIndex node = 0; node < num_nodes; ++node) {
1129 if (node_supply_[node] != 0) {
1130 ArcIndex permuted_arc = PermutedArc(arc);
1131 min_cost_flow.SetArcCapacity(permuted_arc, std::abs(node_supply_[node]));
1132 min_cost_flow.SetArcUnitCost(permuted_arc, 0);
1136 min_cost_flow.SetNodeSupply(source, maximum_flow_);
1137 min_cost_flow.SetNodeSupply(sink, -maximum_flow_);
1138 min_cost_flow.SetCheckFeasibility(
false);
1140 arc_flow_.resize(num_arcs);
1141 if (min_cost_flow.Solve()) {
1142 optimal_cost_ = min_cost_flow.GetOptimalCost();
1143 for (arc = 0; arc < num_arcs; ++arc) {
1144 arc_flow_[arc] = min_cost_flow.Flow(PermutedArc(arc));
1147 return min_cost_flow.status();
1155 return arc_flow_[arc];
1167 return arc_capacity_[arc];
1171 return arc_cost_[arc];
1175 return node_supply_[node];
1178 void SimpleMinCostFlow::ResizeNodeVectors(
NodeIndex node) {
1179 if (node < node_supply_.size())
return;
1180 node_supply_.resize(node + 1);
#define DCHECK_LE(val1, val2)
#define CHECK_EQ(val1, val2)
#define DCHECK_GE(val1, val2)
#define DCHECK_GT(val1, val2)
#define DCHECK(condition)
#define DCHECK_EQ(val1, val2)
#define VLOG(verboselevel)
ArcIndexType AddArc(NodeIndexType tail, NodeIndexType head)
NodeIndexType Tail(const ArcIndexType arc) const
FlowQuantity GetOptimalFlow() const
void SetArcCapacity(ArcIndex arc, FlowQuantity new_capacity)
FlowQuantity Flow(ArcIndex arc) const
void SetCheckResult(bool value)
void SetCheckInput(bool value)
CostValue UnitCost(ArcIndex arc) const
FlowQuantity FeasibleSupply(NodeIndex node) const
FlowQuantity Flow(ArcIndex arc) const
GenericMinCostFlow(const Graph *graph)
Graph::NodeIndex NodeIndex
void SetNodeSupply(NodeIndex node, FlowQuantity supply)
FlowQuantity InitialSupply(NodeIndex node) const
void SetArcFlow(ArcIndex arc, ArcFlowType new_flow)
bool CheckFeasibility(std::vector< NodeIndex > *const infeasible_supply_node, std::vector< NodeIndex > *const infeasible_demand_node)
FlowQuantity Capacity(ArcIndex arc) const
FlowQuantity Supply(NodeIndex node) const
void SetArcUnitCost(ArcIndex arc, ArcScaledCostType unit_cost)
void SetArcCapacity(ArcIndex arc, ArcFlowType new_capacity)
CostValue UnitCost(ArcIndex arc) const
ArcIndex AddArcWithCapacityAndUnitCost(NodeIndex tail, NodeIndex head, FlowQuantity capacity, CostValue unit_cost)
FlowQuantity Flow(ArcIndex arc) const
FlowQuantity MaximumFlow() const
NodeIndex NumNodes() const
SimpleMinCostFlow(NodeIndex reserve_num_nodes=0, ArcIndex reserve_num_arcs=0)
void SetNodeSupply(NodeIndex node, FlowQuantity supply)
NodeIndex Tail(ArcIndex arc) const
FlowQuantity Capacity(ArcIndex arc) const
FlowQuantity Supply(NodeIndex node) const
NodeIndex Head(ArcIndex arc) const
CostValue OptimalCost() const
NodeIndexType Head(const ArcIndexType arc) const
bool Reserve(int64 new_min_index, int64 new_max_index)
GurobiMPCallbackContext * context
ABSL_FLAG(int64, min_cost_flow_alpha, 5, "Divide factor for epsilon at each refine step.")
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
#define IF_STATS_ENABLED(instructions)
#define SCOPED_TIME_STAT(stats)
static ArcIndex ArcReservation(const Graph &graph)
static NodeIndex NodeReservation(const Graph &graph)
static bool IsArcValid(const Graph &graph, ArcIndex arc)
static ArcIndex OppositeArc(const Graph &graph, ArcIndex arc)