25#include "absl/base/attributes.h"
26#include "absl/flags/flag.h"
27#include "absl/log/check.h"
28#include "absl/log/log.h"
29#include "absl/strings/str_format.h"
30#include "absl/strings/string_view.h"
40 "Divide factor for epsilon at each refine step.");
41ABSL_FLAG(
bool, min_cost_flow_check_feasibility,
true,
42 "Check that the graph has enough capacity to send all supplies "
43 "and serve all demands. Also check that the sum of supplies "
44 "is equal to the sum of demands.");
46 "Check that the result is valid.");
50template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
54 alpha_(
absl::GetFlag(FLAGS_min_cost_flow_alpha)),
55 stats_(
"MinCostFlow"),
56 check_feasibility_(
absl::GetFlag(FLAGS_min_cost_flow_check_feasibility)) {
60 const NodeIndex max_num_nodes = graph_->node_capacity();
61 if (max_num_nodes > 0) {
62 first_admissible_arc_ = std::make_unique<ArcIndex[]>(max_num_nodes);
63 std::fill(first_admissible_arc_.get(),
65 node_potential_ = std::make_unique<CostValue[]>(max_num_nodes);
66 node_excess_ = std::make_unique<FlowQuantity[]>(max_num_nodes);
67 initial_node_excess_ = std::make_unique<FlowQuantity[]>(max_num_nodes);
69 const ArcIndex max_num_arcs = graph_->arc_capacity();
70 if (max_num_arcs > 0) {
71 residual_arc_capacity_.Reserve(-max_num_arcs, max_num_arcs - 1);
72 residual_arc_capacity_.SetAll(0);
73 scaled_arc_unit_cost_.Reserve(-max_num_arcs, max_num_arcs - 1);
74 scaled_arc_unit_cost_.SetAll(0);
78template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
81 DCHECK(graph_->IsNodeValid(node));
82 node_excess_[node] = supply;
83 initial_node_excess_[node] = supply;
85 feasibility_checked_ =
false;
88template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
90 ArcIndex arc, ArcScaledCostType unit_cost) {
91 DCHECK(IsArcDirect(arc));
92 scaled_arc_unit_cost_.Set(arc, unit_cost);
93 scaled_arc_unit_cost_.Set(Opposite(arc), -scaled_arc_unit_cost_[arc]);
95 feasibility_checked_ =
false;
98template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
100 ArcIndex arc, ArcFlowType new_capacity) {
101 DCHECK_LE(0, new_capacity);
102 DCHECK(IsArcDirect(arc));
103 const FlowQuantity free_capacity = residual_arc_capacity_[arc];
105 if (capacity_delta == 0) {
109 feasibility_checked_ =
false;
110 const FlowQuantity new_availability = free_capacity + capacity_delta;
111 if (new_availability >= 0) {
117 DCHECK((capacity_delta > 0) ||
118 (capacity_delta < 0 && new_availability >= 0));
119 residual_arc_capacity_.Set(arc, new_availability);
120 DCHECK_LE(0, residual_arc_capacity_[arc]);
124 const FlowQuantity flow = residual_arc_capacity_[Opposite(arc)];
126 residual_arc_capacity_.Set(arc, 0);
127 residual_arc_capacity_.Set(Opposite(arc), new_capacity);
128 node_excess_[Tail(arc)] += flow_excess;
129 node_excess_[Head(arc)] -= flow_excess;
130 DCHECK_LE(0, residual_arc_capacity_[arc]);
131 DCHECK_LE(0, residual_arc_capacity_[Opposite(arc)]);
135template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
137 ArcScaledCostType>::CheckInputConsistency() {
138 const FlowQuantity kMaxFlow = std::numeric_limits<FlowQuantity>::max();
142 FlowQuantity sum_of_supplies = 0;
143 FlowQuantity sum_of_demands = 0;
144 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
145 const FlowQuantity excess = node_excess_[node];
147 sum_of_supplies =
CapAdd(sum_of_supplies, excess);
148 if (sum_of_supplies >= kMaxFlow) {
149 status_ = BAD_CAPACITY_RANGE;
150 LOG(ERROR) <<
"Input consistency error: sum of supplies overflow";
153 }
else if (excess < 0) {
154 sum_of_demands =
CapAdd(sum_of_demands, -excess);
155 if (sum_of_demands >= kMaxFlow) {
156 status_ = BAD_CAPACITY_RANGE;
157 LOG(ERROR) <<
"Input consistency error: sum of demands overflow";
162 if (sum_of_supplies != sum_of_demands) {
163 status_ = UNBALANCED;
164 LOG(ERROR) <<
"Input consistency error: unbalanced problem";
168 std::vector<FlowQuantity> max_node_excess(
169 node_excess_.get(), node_excess_.get() + graph_->num_nodes());
170 std::vector<FlowQuantity> min_node_excess = max_node_excess;
171 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
172 const FlowQuantity capacity = residual_arc_capacity_[arc];
173 CHECK_GE(capacity, 0);
174 if (capacity == 0)
continue;
175 const int tail = graph_->Tail(arc);
176 const int head = graph_->Head(arc);
177 min_node_excess[tail] =
CapSub(min_node_excess[tail], capacity);
178 max_node_excess[head] =
CapAdd(max_node_excess[head], capacity);
181 const int num_nodes = graph_->num_nodes();
182 for (
NodeIndex node = 0; node < num_nodes; ++node) {
183 if (max_node_excess[node] >= kMaxFlow ||
184 min_node_excess[node] <= -kMaxFlow) {
188 if (max_node_excess[node] < std::numeric_limits<ArcFlowType>::max()) {
190 const ArcFlowType upper_bound =
191 std::max<ArcFlowType>(0, max_node_excess[node]);
194 min_node_excess[node] = node_excess_[node];
195 for (
const ArcIndex arc : graph_->OutgoingArcs(node)) {
196 residual_arc_capacity_[arc] =
197 std::min(residual_arc_capacity_[arc], upper_bound);
198 min_node_excess[node] =
199 CapSub(min_node_excess[node], residual_arc_capacity_[arc]);
201 if (min_node_excess[node] > -kMaxFlow)
continue;
203 if (min_node_excess[node] > -std::numeric_limits<ArcFlowType>::max()) {
205 const ArcFlowType upper_bound =
206 std::max<ArcFlowType>(0, -min_node_excess[node]);
209 max_node_excess[node] = node_excess_[node];
210 for (
const ArcIndex arc : graph_->IncomingArcs(node)) {
211 residual_arc_capacity_[arc] =
212 std::min(residual_arc_capacity_[arc], upper_bound);
213 max_node_excess[node] =
214 CapAdd(max_node_excess[node], residual_arc_capacity_[arc]);
216 if (max_node_excess[node] < kMaxFlow)
continue;
219 status_ = BAD_CAPACITY_RANGE;
220 LOG(ERROR) <<
"Maximum in or out flow of node + excess " << node
221 <<
" overflow the FlowQuantity type (int64_t).";
228template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
229bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::CheckResult()
231 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
232 if (node_excess_[node] != 0) {
233 LOG(DFATAL) <<
"node_excess_[" << node <<
"] != 0";
236 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
238 const ArcIndex arc = it.Index();
240 if (residual_arc_capacity_[arc] < 0) {
241 LOG(DFATAL) <<
"residual_arc_capacity_[" << arc <<
"] < 0";
244 if (residual_arc_capacity_[arc] > 0 && ReducedCost(arc) < -epsilon_) {
245 LOG(DFATAL) <<
"residual_arc_capacity_[" << arc
246 <<
"] > 0 && ReducedCost(" << arc <<
") < " << -epsilon_
247 <<
". (epsilon_ = " << epsilon_ <<
").";
251 LOG(DFATAL) << DebugString(
"CheckResult ", arc);
259template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
268 DCHECK_GE(node_excess_[node], 0);
269 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
271 const ArcIndex arc = it.Index();
272 DCHECK(!IsAdmissible(arc)) << DebugString(
"CheckRelabelPrecondition:", arc);
277template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
279GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::DebugString(
280 absl::string_view context, ArcIndex arc)
const {
286 const CostValue reduced_cost = scaled_arc_unit_cost_[arc] +
287 node_potential_[tail] - node_potential_[head];
288 return absl::StrFormat(
289 "%s Arc %d, from %d to %d, "
290 "Capacity = %d, Residual capacity = %d, "
291 "Flow = residual capacity for reverse arc = %d, "
292 "Height(tail) = %d, Height(head) = %d, "
293 "Excess(tail) = %d, Excess(head) = %d, "
294 "Cost = %d, Reduced cost = %d, ",
295 context, arc, tail, head, Capacity(arc),
296 static_cast<FlowQuantity
>(residual_arc_capacity_[arc]), Flow(arc),
297 node_potential_[tail], node_potential_[head], node_excess_[tail],
298 node_excess_[head],
static_cast<CostValue>(scaled_arc_unit_cost_[arc]),
302template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
304 ArcScaledCostType>::CheckFeasibility() {
316 feasibility_checked_ =
false;
317 ArcIndex num_extra_arcs = 0;
318 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
319 if (initial_node_excess_[node] != 0) {
323 const NodeIndex num_nodes_in_max_flow = graph_->num_nodes() + 2;
324 const ArcIndex num_arcs_in_max_flow = graph_->num_arcs() + num_extra_arcs;
325 const NodeIndex source = num_nodes_in_max_flow - 2;
326 const NodeIndex sink = num_nodes_in_max_flow - 1;
327 using CheckerGraph = ::util::ReverseArcListGraph<>;
328 CheckerGraph checker_graph(num_nodes_in_max_flow, num_arcs_in_max_flow);
331 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
332 const ArcIndex new_arc =
333 checker_graph.AddArc(graph_->Tail(arc), graph_->Head(arc));
334 DCHECK_EQ(arc, new_arc);
335 checker.SetArcCapacity(new_arc, Capacity(arc));
337 FlowQuantity total_demand = 0;
338 FlowQuantity total_supply = 0;
340 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
341 const FlowQuantity supply = initial_node_excess_[node];
343 const ArcIndex new_arc = checker_graph.AddArc(source, node);
344 checker.SetArcCapacity(new_arc, supply);
345 total_supply += supply;
346 }
else if (supply < 0) {
347 const ArcIndex new_arc = checker_graph.AddArc(node, sink);
348 checker.SetArcCapacity(new_arc, -supply);
349 total_demand -= supply;
352 if (total_supply != total_demand) {
353 LOG(DFATAL) <<
"total_supply(" << total_supply <<
") != total_demand("
354 << total_demand <<
").";
357 if (!checker.Solve()) {
358 LOG(DFATAL) <<
"Max flow could not be computed.";
361 const FlowQuantity optimal_max_flow = checker.GetOptimalFlow();
362 feasibility_checked_ =
true;
363 return optimal_max_flow == total_supply;
366template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
369 if (IsArcDirect(arc)) {
370 return residual_arc_capacity_[Opposite(arc)];
372 return -residual_arc_capacity_[arc];
377template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
380 if (IsArcDirect(arc)) {
381 return residual_arc_capacity_[arc] + residual_arc_capacity_[Opposite(arc)];
387template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
390 DCHECK(IsArcValid(arc));
391 DCHECK_EQ(uint64_t{1}, cost_scaling_factor_);
392 return scaled_arc_unit_cost_[arc];
395template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
398 DCHECK(graph_->IsNodeValid(node));
399 return node_excess_[node];
402template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
403bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsAdmissible(
404 ArcIndex arc)
const {
405 return FastIsAdmissible(arc, node_potential_[Tail(arc)]);
408template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
409bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::
410 FastIsAdmissible(ArcIndex arc,
CostValue tail_potential)
const {
411 DCHECK_EQ(node_potential_[Tail(arc)], tail_potential);
412 return residual_arc_capacity_[arc] > 0 &&
413 FastReducedCost(arc, tail_potential) < 0;
416template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
417bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsActive(
419 return node_excess_[node] > 0;
422template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
423auto GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::ReducedCost(
425 return FastReducedCost(arc, node_potential_[Tail(arc)]);
428template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
429auto GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::FastReducedCost(
431 DCHECK_EQ(node_potential_[Tail(arc)], tail_potential);
432 DCHECK(graph_->IsNodeValid(Tail(arc)));
433 DCHECK(graph_->IsNodeValid(Head(arc)));
434 DCHECK_LE(node_potential_[Tail(arc)], 0) << DebugString(
"ReducedCost:", arc);
435 DCHECK_LE(node_potential_[Head(arc)], 0) << DebugString(
"ReducedCost:", arc);
439 return scaled_arc_unit_cost_[arc] + tail_potential -
440 node_potential_[Head(arc)];
443template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
446 OutgoingOrOppositeIncomingArcIterator arc_it(*graph_, node);
447 return arc_it.Index();
450template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
452 if (!CheckInputConsistency()) {
455 if (check_feasibility_ && !CheckFeasibility()) {
461 std::fill(node_potential_.get(),
462 node_potential_.get() + graph_->node_capacity(), 0);
464 ResetFirstAdmissibleArcs();
465 if (!ScaleCosts())
return false;
466 if (!Optimize())
return false;
470 if (absl::GetFlag(FLAGS_min_cost_flow_check_result) && !CheckResult()) {
481template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
491 const CostValue kMaxCost = std::numeric_limits<CostValue>::max();
492 const CostValue kMinCost = std::numeric_limits<CostValue>::min();
493 for (
ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
494 const CostValue flow_on_arc = residual_arc_capacity_[Opposite(arc)];
496 CapProd(scaled_arc_unit_cost_[arc], flow_on_arc);
497 if (flow_cost == kMaxCost || flow_cost == kMinCost)
return kMaxCost;
498 total_flow_cost =
CapAdd(flow_cost, total_flow_cost);
499 if (total_flow_cost == kMaxCost || total_flow_cost == kMinCost) {
503 return total_flow_cost;
506template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
508 ArcScaledCostType>::ResetFirstAdmissibleArcs() {
509 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
510 first_admissible_arc_[node] = GetFirstOutgoingOrOppositeIncomingArc(node);
514template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
515bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::ScaleCosts() {
517 cost_scaling_factor_ = scale_prices_ ? graph_->num_nodes() + 1 : 1;
519 VLOG(3) <<
"Number of nodes in the graph = " << graph_->num_nodes();
520 VLOG(3) <<
"Number of arcs in the graph = " << graph_->num_arcs();
522 std::numeric_limits<CostValue>::max() / (2 * cost_scaling_factor_);
523 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
524 if (scaled_arc_unit_cost_[arc] > threshold ||
525 scaled_arc_unit_cost_[arc] < -threshold) {
526 status_ = BAD_COST_RANGE;
529 const CostValue cost = scaled_arc_unit_cost_[arc] * cost_scaling_factor_;
530 scaled_arc_unit_cost_.Set(arc, cost);
531 scaled_arc_unit_cost_.Set(Opposite(arc), -cost);
538 overflow_threshold_ = std::numeric_limits<CostValue>::min() + epsilon_;
540 VLOG(3) <<
"Initial epsilon = " << epsilon_;
541 VLOG(3) <<
"Cost scaling factor = " << cost_scaling_factor_;
545template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
546void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::UnscaleCosts() {
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);
553 cost_scaling_factor_ = 1;
556template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
557bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Optimize() {
559 num_relabels_since_last_price_update_ = 0;
562 epsilon_ = std::max(epsilon_ / alpha_, kEpsilonMin);
563 VLOG(3) <<
"Epsilon changed to: " << epsilon_;
564 if (!Refine())
return false;
565 }
while (epsilon_ != 1LL);
569template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
571 ArcScaledCostType>::SaturateAdmissibleArcs() {
573 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
574 const CostValue tail_potential = node_potential_[node];
575 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node,
576 first_admissible_arc_[node]);
577 it.Ok(); it.Next()) {
578 const ArcIndex arc = it.Index();
579 if (FastIsAdmissible(arc, tail_potential)) {
580 FastPushFlow(residual_arc_capacity_[arc], arc, node);
594template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
595void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::PushFlow(
596 FlowQuantity flow, ArcIndex arc) {
598 FastPushFlow(flow, arc, Tail(arc));
601template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
602void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::FastPushFlow(
603 FlowQuantity flow, ArcIndex arc,
NodeIndex tail) {
605 DCHECK_EQ(Tail(arc), tail);
606 DCHECK_GT(residual_arc_capacity_[arc], 0);
607 DCHECK_LE(flow, residual_arc_capacity_[arc]);
609 residual_arc_capacity_.Set(arc, residual_arc_capacity_[arc] - flow);
611 const ArcIndex opposite = Opposite(arc);
612 residual_arc_capacity_.Set(opposite, residual_arc_capacity_[opposite] + flow);
614 node_excess_[tail] -= flow;
615 node_excess_[Head(arc)] += flow;
618template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
620 ArcScaledCostType>::InitializeActiveNodeStack() {
622 DCHECK(active_nodes_.empty());
623 for (
NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
624 if (IsActive(node)) {
625 active_nodes_.push(node);
630template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
631bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::UpdatePrices() {
650 const NodeIndex num_nodes = graph_->num_nodes();
651 std::vector<NodeIndex> bfs_queue;
652 std::vector<bool> node_in_queue(num_nodes,
false);
655 const CostValue kMinCostValue = std::numeric_limits<CostValue>::min();
656 std::vector<CostValue> min_non_admissible_potential(num_nodes, kMinCostValue);
657 std::vector<NodeIndex> nodes_to_process;
660 FlowQuantity remaining_excess = 0;
663 for (
NodeIndex node = 0; node < num_nodes; ++node) {
664 if (node_excess_[node] < 0) {
665 bfs_queue.push_back(node);
666 node_in_queue[node] =
true;
669 remaining_excess -= node_excess_[node];
680 while (remaining_excess > 0) {
685 for (; queue_index < bfs_queue.size(); ++queue_index) {
686 DCHECK_GE(num_nodes, bfs_queue.size());
687 const NodeIndex node = bfs_queue[queue_index];
688 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
691 if (node_in_queue[head])
continue;
692 const ArcIndex opposite_arc = Opposite(it.Index());
693 if (residual_arc_capacity_[opposite_arc] > 0) {
694 node_potential_[head] += potential_delta;
695 if (node_potential_[head] < overflow_threshold_) {
696 status_ = BAD_COST_RANGE;
699 if (ReducedCost(opposite_arc) < 0) {
700 DCHECK(IsAdmissible(opposite_arc));
705 remaining_excess -= node_excess_[head];
706 if (remaining_excess == 0) {
707 node_potential_[head] -= potential_delta;
710 bfs_queue.push_back(head);
711 node_in_queue[head] =
true;
712 if (potential_delta < 0) {
713 first_admissible_arc_[head] =
714 GetFirstOutgoingOrOppositeIncomingArc(head);
719 node_potential_[head] -= potential_delta;
720 if (min_non_admissible_potential[head] == kMinCostValue) {
721 nodes_to_process.push_back(head);
723 min_non_admissible_potential[head] = std::max(
724 min_non_admissible_potential[head],
725 node_potential_[node] - scaled_arc_unit_cost_[opposite_arc]);
729 if (remaining_excess == 0)
break;
731 if (remaining_excess == 0)
break;
735 CostValue max_potential_diff = kMinCostValue;
736 for (
int i = 0;
i < nodes_to_process.size(); ++
i) {
738 if (node_in_queue[node])
continue;
740 std::max(max_potential_diff,
741 min_non_admissible_potential[node] - node_potential_[node]);
742 if (max_potential_diff == potential_delta)
break;
744 DCHECK_LE(max_potential_diff, potential_delta);
745 potential_delta = max_potential_diff - epsilon_;
754 for (
int i = 0;
i < nodes_to_process.size(); ++
i) {
756 if (node_in_queue[node])
continue;
757 if (node_potential_[node] + potential_delta <
758 min_non_admissible_potential[node]) {
759 node_potential_[node] += potential_delta;
760 if (node_potential_[node] < overflow_threshold_) {
761 status_ = BAD_COST_RANGE;
764 first_admissible_arc_[node] =
765 GetFirstOutgoingOrOppositeIncomingArc(node);
766 bfs_queue.push_back(node);
767 node_in_queue[node] =
true;
768 remaining_excess -= node_excess_[node];
773 nodes_to_process[index] = node;
776 nodes_to_process.resize(index);
780 if (potential_delta == 0)
return true;
781 for (
NodeIndex node = 0; node < num_nodes; ++node) {
782 if (!node_in_queue[node]) {
783 node_potential_[node] += potential_delta;
784 if (node_potential_[node] < overflow_threshold_) {
785 status_ = BAD_COST_RANGE;
788 first_admissible_arc_[node] = GetFirstOutgoingOrOppositeIncomingArc(node);
794template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
795bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Refine() {
797 SaturateAdmissibleArcs();
798 InitializeActiveNodeStack();
800 const NodeIndex num_nodes = graph_->num_nodes();
801 while (!active_nodes_.empty()) {
803 if (num_relabels_since_last_price_update_ >= num_nodes) {
804 num_relabels_since_last_price_update_ = 0;
805 if (use_price_update_) {
806 if (!UpdatePrices())
return false;
809 const NodeIndex node = active_nodes_.top();
811 DCHECK(IsActive(node));
812 if (!Discharge(node))
return false;
817template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
818bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Discharge(
824 DCHECK(IsActive(node));
825 const CostValue tail_potential = node_potential_[node];
826 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node,
827 first_admissible_arc_[node]);
828 it.Ok(); it.Next()) {
829 const ArcIndex arc = it.Index();
830 if (!FastIsAdmissible(arc, tail_potential))
continue;
836 if (node_excess_[head] >= 0 && !NodeHasAdmissibleArc(head)) {
838 if (!Relabel(head))
return false;
839 if (!FastIsAdmissible(arc, tail_potential))
continue;
842 const FlowQuantity delta =
843 std::min(node_excess_[node],
844 static_cast<FlowQuantity
>(residual_arc_capacity_[arc]));
845 const bool head_active_before_push = IsActive(head);
846 FastPushFlow(delta, arc, node);
847 if (IsActive(head) && !head_active_before_push) {
848 active_nodes_.push(head);
851 if (node_excess_[node] == 0) {
853 first_admissible_arc_[node] = arc;
857 if (!Relabel(node))
return false;
861template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
865 const CostValue tail_potential = node_potential_[node];
866 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node,
867 first_admissible_arc_[node]);
868 it.Ok(); it.Next()) {
869 const ArcIndex arc = it.Index();
870 if (FastIsAdmissible(arc, tail_potential)) {
871 first_admissible_arc_[node] = arc;
879template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
880bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Relabel(
883 DCHECK(CheckRelabelPrecondition(node));
884 ++num_relabels_since_last_price_update_;
891 const CostValue guaranteed_new_potential = node_potential_[node] - epsilon_;
892 if (guaranteed_new_potential < overflow_threshold_) {
893 status_ = BAD_COST_RANGE;
902 const CostValue kMinCostValue = std::numeric_limits<CostValue>::min();
903 CostValue min_non_admissible_potential = kMinCostValue;
908 CostValue previous_min_non_admissible_potential = kMinCostValue;
911 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
913 const ArcIndex arc = it.Index();
914 if (residual_arc_capacity_[arc] > 0) {
915 const CostValue min_non_admissible_potential_for_arc =
916 node_potential_[Head(arc)] - scaled_arc_unit_cost_[arc];
920 if (min_non_admissible_potential_for_arc > guaranteed_new_potential) {
924 node_potential_[node] = guaranteed_new_potential;
925 first_admissible_arc_[node] = arc;
928 if (min_non_admissible_potential_for_arc > min_non_admissible_potential) {
929 previous_min_non_admissible_potential = min_non_admissible_potential;
930 min_non_admissible_potential = min_non_admissible_potential_for_arc;
940 if (min_non_admissible_potential == kMinCostValue) {
941 if (node_excess_[node] != 0) {
951 node_potential_[node] = guaranteed_new_potential;
960 if (min_non_admissible_potential < overflow_threshold_) {
961 status_ = BAD_COST_RANGE;
964 const CostValue new_potential = min_non_admissible_potential - epsilon_;
965 if (new_potential < overflow_threshold_) {
966 status_ = BAD_COST_RANGE;
969 node_potential_[node] = new_potential;
970 if (previous_min_non_admissible_potential <= new_potential) {
971 first_admissible_arc_[node] = first_arc;
974 first_admissible_arc_[node] = GetFirstOutgoingOrOppositeIncomingArc(node);
979template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
981GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Opposite(
982 ArcIndex arc)
const {
983 return graph_->OppositeArc(arc);
986template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
987bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsArcValid(
988 ArcIndex arc)
const {
989 return graph_->IsArcValid(arc);
992template <
typename Graph,
typename ArcFlowType,
typename ArcScaledCostType>
993bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsArcDirect(
994 ArcIndex arc)
const {
995 DCHECK(IsArcValid(arc));
1006 ::util::ReverseArcStaticGraph<uint16_t, int32_t>>;
1012 ::util::ReverseArcStaticGraph<uint16_t, int32_t>,
1018 if (reserve_num_nodes > 0) {
1019 node_supply_.reserve(reserve_num_nodes);
1021 if (reserve_num_arcs > 0) {
1022 arc_tail_.reserve(reserve_num_arcs);
1023 arc_head_.reserve(reserve_num_arcs);
1024 arc_capacity_.reserve(reserve_num_arcs);
1025 arc_cost_.reserve(reserve_num_arcs);
1026 arc_permutation_.reserve(reserve_num_arcs);
1027 arc_flow_.reserve(reserve_num_arcs);
1032 ResizeNodeVectors(node);
1033 node_supply_[node] = supply;
1039 ResizeNodeVectors(std::max(tail, head));
1040 const ArcIndex arc = arc_tail_.size();
1041 arc_tail_.push_back(tail);
1042 arc_head_.push_back(head);
1043 arc_capacity_.push_back(capacity);
1044 arc_cost_.push_back(unit_cost);
1049 arc_capacity_[arc] = capacity;
1053 return arc < arc_permutation_.size() ? arc_permutation_[arc] : arc;
1057 SupplyAdjustment adjustment) {
1061 const NodeIndex num_nodes = node_supply_.size();
1062 const ArcIndex num_arcs = arc_capacity_.size();
1063 if (num_nodes == 0)
return OPTIMAL;
1065 int supply_node_count = 0, demand_node_count = 0;
1067 for (
NodeIndex node = 0; node < num_nodes; ++node) {
1068 if (node_supply_[node] > 0) {
1069 ++supply_node_count;
1070 total_supply += node_supply_[node];
1071 }
else if (node_supply_[node] < 0) {
1072 ++demand_node_count;
1073 total_demand -= node_supply_[node];
1076 if (adjustment == DONT_ADJUST && total_supply != total_demand) {
1091 const ArcIndex augmented_num_arcs =
1092 num_arcs + supply_node_count + demand_node_count;
1095 const NodeIndex augmented_num_nodes = num_nodes + 2;
1097 Graph graph(augmented_num_nodes, augmented_num_arcs);
1098 for (
ArcIndex arc = 0; arc < num_arcs; ++arc) {
1099 graph.AddArc(arc_tail_[arc], arc_head_[arc]);
1102 for (
NodeIndex node = 0; node < num_nodes; ++node) {
1103 if (node_supply_[node] > 0) {
1104 graph.AddArc(source, node);
1105 }
else if (node_supply_[node] < 0) {
1106 graph.AddArc(node, sink);
1110 graph.Build(&arc_permutation_);
1113 GenericMaxFlow<Graph> max_flow(&graph, source, sink);
1115 for (arc = 0; arc < num_arcs; ++arc) {
1116 max_flow.SetArcCapacity(PermutedArc(arc), arc_capacity_[arc]);
1118 for (
NodeIndex node = 0; node < num_nodes; ++node) {
1119 if (node_supply_[node] != 0) {
1120 max_flow.SetArcCapacity(PermutedArc(arc), std::abs(node_supply_[node]));
1124 CHECK_EQ(arc, augmented_num_arcs);
1125 if (!max_flow.Solve()) {
1126 LOG(ERROR) <<
"Max flow could not be computed.";
1127 switch (max_flow.status()) {
1132 <<
"Max flow failed but claimed to have an optimal solution";
1133 ABSL_FALLTHROUGH_INTENDED;
1138 maximum_flow_ = max_flow.GetOptimalFlow();
1141 if (adjustment == DONT_ADJUST && maximum_flow_ != total_supply) {
1145 GenericMinCostFlow<Graph> min_cost_flow(&graph);
1147 for (arc = 0; arc < num_arcs; ++arc) {
1148 ArcIndex permuted_arc = PermutedArc(arc);
1149 min_cost_flow.SetArcUnitCost(permuted_arc, arc_cost_[arc]);
1150 min_cost_flow.SetArcCapacity(permuted_arc, arc_capacity_[arc]);
1152 for (
NodeIndex node = 0; node < num_nodes; ++node) {
1153 if (node_supply_[node] != 0) {
1154 ArcIndex permuted_arc = PermutedArc(arc);
1155 min_cost_flow.SetArcCapacity(permuted_arc, std::abs(node_supply_[node]));
1156 min_cost_flow.SetArcUnitCost(permuted_arc, 0);
1160 min_cost_flow.SetNodeSupply(source, maximum_flow_);
1161 min_cost_flow.SetNodeSupply(sink, -maximum_flow_);
1162 min_cost_flow.SetCheckFeasibility(
false);
1163 min_cost_flow.SetPriceScaling(scale_prices_);
1165 arc_flow_.resize(num_arcs);
1166 if (min_cost_flow.Solve()) {
1167 optimal_cost_ = min_cost_flow.GetOptimalCost();
1168 for (arc = 0; arc < num_arcs; ++arc) {
1169 arc_flow_[arc] = min_cost_flow.Flow(PermutedArc(arc));
1173 return min_cost_flow.status();
1177 return optimal_cost_;
1181 return maximum_flow_;
1185 return arc_flow_[arc];
1190 return node_supply_.size();
1194 return arc_tail_.size();
1198 return arc_tail_[arc];
1202 return arc_head_[arc];
1207 return arc_capacity_[arc];
1211 return arc_cost_[arc];
1216 return node_supply_[node];
1219void SimpleMinCostFlow::ResizeNodeVectors(
NodeIndex node) {
1220 if (node < node_supply_.size())
return;
1221 node_supply_.resize(node + 1);
FlowQuantity Flow(ArcIndex arc) const
void SetArcUnitCost(ArcIndex arc, ArcScaledCostType unit_cost)
Sets the unit cost for the given arc.
GenericMinCostFlow(const Graph *graph)
bool Solve()
Solves the problem, returning true if a min-cost flow could be found.
void SetArcCapacity(ArcIndex arc, ArcFlowType new_capacity)
Sets the capacity for the given arc.
void SetNodeSupply(NodeIndex node, FlowQuantity supply)
const Graph * graph() const
Returns the graph associated to the current object.
CostValue UnitCost(ArcIndex arc) const
Returns the unscaled cost for the given arc.
FlowQuantity Supply(NodeIndex node) const
Graph::NodeIndex NodeIndex
FlowQuantity Capacity(ArcIndex arc) const
We use the equations given in the comment of residual_arc_capacity_.
CostValue GetOptimalCost()
NodeIndex Tail(ArcIndex arc) const
CostValue OptimalCost() const
void SetNodeSupply(NodeIndex node, FlowQuantity supply)
void SetArcCapacity(ArcIndex arc, FlowQuantity capacity)
FlowQuantity Flow(ArcIndex arc) const
FlowQuantity Capacity(ArcIndex arc) const
CostValue UnitCost(ArcIndex arc) const
NodeIndex NumNodes() const
ArcIndex AddArcWithCapacityAndUnitCost(NodeIndex tail, NodeIndex head, FlowQuantity capacity, CostValue unit_cost)
FlowQuantity MaximumFlow() const
NodeIndex Head(ArcIndex arc) const
SimpleMinCostFlow(NodeIndex reserve_num_nodes=0, ArcIndex reserve_num_arcs=0)
FlowQuantity Supply(NodeIndex node) const
static constexpr int32_t kNilArc
static constexpr bool kHasNegativeReverseArcs
ABSL_FLAG(int64_t, min_cost_flow_alpha, 5, "Divide factor for epsilon at each refine step.")
In SWIG mode, we don't want anything besides these top-level includes.
int64_t CapAdd(int64_t x, int64_t y)
@ INFEASIBLE
A solution could not be found.
int64_t CapSub(int64_t x, int64_t y)
BlossomGraph::CostValue CostValue
BlossomGraph::NodeIndex NodeIndex
int64_t CapProd(int64_t x, int64_t y)
util::ReverseArcStaticGraph Graph
Type of graph to use.
#define IF_STATS_ENABLED(instructions)
#define SCOPED_TIME_STAT(stats)