Google OR-Tools v9.14
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
min_cost_flow.cc
Go to the documentation of this file.
1// Copyright 2010-2025 Google LLC
2// Licensed under the Apache License, Version 2.0 (the "License");
3// you may not use this file except in compliance with the License.
4// You may obtain a copy of the License at
5//
6// http://www.apache.org/licenses/LICENSE-2.0
7//
8// Unless required by applicable law or agreed to in writing, software
9// distributed under the License is distributed on an "AS IS" BASIS,
10// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11// See the License for the specific language governing permissions and
12// limitations under the License.
13
15
16#include <algorithm>
17#include <cmath>
18#include <cstdint>
19#include <cstdlib>
20#include <limits>
21#include <memory>
22#include <string>
23#include <vector>
24
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"
33#include "ortools/graph/graph.h"
35#include "ortools/util/stats.h"
36
37// TODO(user): Remove these flags and expose the parameters in the API.
38// New clients, please do not use these flags!
39ABSL_FLAG(int64_t, min_cost_flow_alpha, 5,
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.");
45ABSL_FLAG(bool, min_cost_flow_check_result, true,
46 "Check that the result is valid.");
47
48namespace operations_research {
49
50template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
52 const Graph* graph)
53 : graph_(graph),
54 alpha_(absl::GetFlag(FLAGS_min_cost_flow_alpha)),
55 stats_("MinCostFlow"),
56 check_feasibility_(absl::GetFlag(FLAGS_min_cost_flow_check_feasibility)) {
57 // This class assumes we have negative reverse arcs.
58 static_assert(Graph::kHasNegativeReverseArcs);
59
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(),
64 first_admissible_arc_.get() + max_num_nodes, Graph::kNilArc);
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);
68 }
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);
75 }
76}
77
78template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
80 NodeIndex node, FlowQuantity supply) {
81 DCHECK(graph_->IsNodeValid(node));
82 node_excess_[node] = supply;
83 initial_node_excess_[node] = supply;
84 status_ = NOT_SOLVED;
85 feasibility_checked_ = false;
86}
87
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]);
94 status_ = NOT_SOLVED;
95 feasibility_checked_ = false;
96}
97
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];
104 const FlowQuantity capacity_delta = new_capacity - Capacity(arc);
105 if (capacity_delta == 0) {
106 return; // Nothing to do.
107 }
108 status_ = NOT_SOLVED;
109 feasibility_checked_ = false;
110 const FlowQuantity new_availability = free_capacity + capacity_delta;
111 if (new_availability >= 0) {
112 // The above condition is true when one of two following holds:
113 // 1/ (capacity_delta > 0), meaning we are increasing the capacity
114 // 2/ (capacity_delta < 0 && free_capacity + capacity_delta >= 0)
115 // meaning we are reducing the capacity, but that the capacity
116 // reduction is not larger than the free capacity.
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]);
121 } else {
122 // We have to reduce the flow on the arc, and update the excesses
123 // accordingly.
124 const FlowQuantity flow = residual_arc_capacity_[Opposite(arc)];
125 const FlowQuantity flow_excess = flow - new_capacity;
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)]);
132 }
133}
134
135template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
136bool GenericMinCostFlow<Graph, ArcFlowType,
137 ArcScaledCostType>::CheckInputConsistency() {
138 const FlowQuantity kMaxFlow = std::numeric_limits<FlowQuantity>::max();
139
140 // First lets make sure supply == demand and the total supply/demand do not
141 // overflow.
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];
146 if (excess > 0) {
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";
151 return false;
152 }
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";
158 return false;
159 }
160 }
161 }
162 if (sum_of_supplies != sum_of_demands) {
163 status_ = UNBALANCED;
164 LOG(ERROR) << "Input consistency error: unbalanced problem";
165 return false;
166 }
167
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);
179 }
180
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) {
185 // Try to fix it.
186 // Some user just use arc with infinite capacity out of the source/sink.
187 // This is why we use CappedCapacity() for a name.
188 if (max_node_excess[node] < std::numeric_limits<ArcFlowType>::max()) {
189 // There is no point having an outgoing arc with more than this.
190 const ArcFlowType upper_bound =
191 std::max<ArcFlowType>(0, max_node_excess[node]);
192
193 // Adjust and recompute min_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]);
200 }
201 if (min_node_excess[node] > -kMaxFlow) continue;
202 }
203 if (min_node_excess[node] > -std::numeric_limits<ArcFlowType>::max()) {
204 // There is no point having an incoming arc with more than this.
205 const ArcFlowType upper_bound =
206 std::max<ArcFlowType>(0, -min_node_excess[node]);
207
208 // Adjust and recompute max_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]);
215 }
216 if (max_node_excess[node] < kMaxFlow) continue;
217 }
218
219 status_ = BAD_CAPACITY_RANGE;
220 LOG(ERROR) << "Maximum in or out flow of node + excess " << node
221 << " overflow the FlowQuantity type (int64_t).";
222 return false;
223 }
224 }
225 return true;
226}
227
228template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
229bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::CheckResult()
230 const {
231 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
232 if (node_excess_[node] != 0) {
233 LOG(DFATAL) << "node_excess_[" << node << "] != 0";
234 return false;
235 }
236 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
237 it.Next()) {
238 const ArcIndex arc = it.Index();
239 bool ok = true;
240 if (residual_arc_capacity_[arc] < 0) {
241 LOG(DFATAL) << "residual_arc_capacity_[" << arc << "] < 0";
242 ok = false;
243 }
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_ << ").";
248 ok = false;
249 }
250 if (!ok) {
251 LOG(DFATAL) << DebugString("CheckResult ", arc);
252 return false;
253 }
254 }
255 }
256 return true;
257}
258
259template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
262 // Note that the classical Relabel precondition assumes IsActive(node), i.e.,
263 // the node_excess_[node] > 0. However, to implement the Push Look-Ahead
264 // heuristic, we can relax this condition as explained in the section 4.3 of
265 // the article "An Efficient Implementation of a Scaling Minimum-Cost Flow
266 // Algorithm", A.V. Goldberg, Journal of Algorithms 22(1), January 1997, pp.
267 // 1-29.
268 DCHECK_GE(node_excess_[node], 0);
269 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
270 it.Next()) {
271 const ArcIndex arc = it.Index();
272 DCHECK(!IsAdmissible(arc)) << DebugString("CheckRelabelPrecondition:", arc);
273 }
274 return true;
275}
276
277template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
278std::string
279GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::DebugString(
280 absl::string_view context, ArcIndex arc) const {
281 const NodeIndex tail = Tail(arc);
282 const NodeIndex head = Head(arc);
283 // Reduced cost is computed directly without calling ReducedCost to avoid
284 // recursive calls between ReducedCost and DebugString in case a DCHECK in
285 // ReducedCost fails.
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]),
299 reduced_cost);
300}
301
302template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
303bool GenericMinCostFlow<Graph, ArcFlowType,
304 ArcScaledCostType>::CheckFeasibility() {
305 SCOPED_TIME_STAT(&stats_);
306 // Create a new graph, which is a copy of graph_, with the following
307 // modifications:
308 // Two nodes are added: a source and a sink.
309 // The source is linked to each supply node (whose supply > 0) by an arc whose
310 // capacity is equal to the supply at the supply node.
311 // The sink is linked to each demand node (whose supply < 0) by an arc whose
312 // capacity is the demand (-supply) at the demand node.
313 // There are no supplies or demands or costs in the graph, as we will run
314 // max-flow.
315 // TODO(user): make it possible to share a graph by MaxFlow and MinCostFlow.
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) {
320 ++num_extra_arcs;
321 }
322 }
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);
329 GenericMaxFlow<CheckerGraph> checker(&checker_graph, source, sink);
330 // Copy graph_ to checker_graph.
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));
336 }
337 FlowQuantity total_demand = 0;
338 FlowQuantity total_supply = 0;
339 // Create the source-to-supply node arcs and the demand-node-to-sink arcs.
340 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
341 const FlowQuantity supply = initial_node_excess_[node];
342 if (supply > 0) {
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;
350 }
351 }
352 if (total_supply != total_demand) {
353 LOG(DFATAL) << "total_supply(" << total_supply << ") != total_demand("
354 << total_demand << ").";
355 return false;
356 }
357 if (!checker.Solve()) {
358 LOG(DFATAL) << "Max flow could not be computed.";
359 return false;
360 }
361 const FlowQuantity optimal_max_flow = checker.GetOptimalFlow();
362 feasibility_checked_ = true;
363 return optimal_max_flow == total_supply;
364}
365
366template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
368 ArcIndex arc) const -> FlowQuantity {
369 if (IsArcDirect(arc)) {
370 return residual_arc_capacity_[Opposite(arc)];
371 } else {
372 return -residual_arc_capacity_[arc];
373 }
374}
375
376// We use the equations given in the comment of residual_arc_capacity_.
377template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
379 ArcIndex arc) const -> FlowQuantity {
380 if (IsArcDirect(arc)) {
381 return residual_arc_capacity_[arc] + residual_arc_capacity_[Opposite(arc)];
382 } else {
383 return 0;
384 }
385}
386
387template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
389 ArcIndex arc) const -> CostValue {
390 DCHECK(IsArcValid(arc));
391 DCHECK_EQ(uint64_t{1}, cost_scaling_factor_);
392 return scaled_arc_unit_cost_[arc];
393}
394
395template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
397 NodeIndex node) const -> FlowQuantity {
398 DCHECK(graph_->IsNodeValid(node));
399 return node_excess_[node];
400}
401
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)]);
406}
407
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;
414}
415
416template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
417bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsActive(
418 NodeIndex node) const {
419 return node_excess_[node] > 0;
420}
421
422template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
423auto GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::ReducedCost(
424 ArcIndex arc) const -> CostValue {
425 return FastReducedCost(arc, node_potential_[Tail(arc)]);
426}
427
428template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
429auto GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::FastReducedCost(
430 ArcIndex arc, CostValue tail_potential) const -> CostValue {
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);
436
437 // Note that there should be no overflow independently of the order of
438 // operations if potentials are in [overflow_threshold_, 0].
439 return scaled_arc_unit_cost_[arc] + tail_potential -
440 node_potential_[Head(arc)];
441}
442
443template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
446 OutgoingOrOppositeIncomingArcIterator arc_it(*graph_, node);
447 return arc_it.Index();
448}
449
450template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
452 if (!CheckInputConsistency()) {
453 return false;
454 }
455 if (check_feasibility_ && !CheckFeasibility()) {
456 status_ = INFEASIBLE;
457 return false;
458 }
459
460 status_ = NOT_SOLVED;
461 std::fill(node_potential_.get(),
462 node_potential_.get() + graph_->node_capacity(), 0);
463
464 ResetFirstAdmissibleArcs();
465 if (!ScaleCosts()) return false;
466 if (!Optimize()) return false;
467 DCHECK_EQ(status_, NOT_SOLVED);
468 status_ = OPTIMAL;
469
470 if (absl::GetFlag(FLAGS_min_cost_flow_check_result) && !CheckResult()) {
471 status_ = BAD_RESULT;
472 UnscaleCosts();
473 return false;
474 }
475 UnscaleCosts();
476
477 IF_STATS_ENABLED(VLOG(1) << stats_.StatString());
478 return true;
479}
480
481template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
483 -> CostValue {
484 if (status_ != OPTIMAL) {
485 return 0;
486 }
487
488 // The total cost of the flow.
489 // We cap the result if its overflow.
490 CostValue total_flow_cost = 0;
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)];
495 const CostValue flow_cost =
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) {
500 return kMaxCost;
501 }
502 }
503 return total_flow_cost;
504}
505
506template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
507void GenericMinCostFlow<Graph, ArcFlowType,
508 ArcScaledCostType>::ResetFirstAdmissibleArcs() {
509 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
510 first_admissible_arc_[node] = GetFirstOutgoingOrOppositeIncomingArc(node);
511 }
512}
513
514template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
515bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::ScaleCosts() {
516 SCOPED_TIME_STAT(&stats_);
517 cost_scaling_factor_ = scale_prices_ ? graph_->num_nodes() + 1 : 1;
518 epsilon_ = 1LL;
519 VLOG(3) << "Number of nodes in the graph = " << graph_->num_nodes();
520 VLOG(3) << "Number of arcs in the graph = " << graph_->num_arcs();
521 const CostValue threshold =
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;
527 return false;
528 }
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);
532 epsilon_ = std::max(epsilon_, MathUtil::Abs(cost));
533 }
534
535 // Since epsilon_ is currently the largest scaled cost, as long as our
536 // node potentials stay above this threshold, we can correctly compute
537 // potential - epsilon or cost + potential_diff.
538 overflow_threshold_ = std::numeric_limits<CostValue>::min() + epsilon_;
539
540 VLOG(3) << "Initial epsilon = " << epsilon_;
541 VLOG(3) << "Cost scaling factor = " << cost_scaling_factor_;
542 return true;
543}
544
545template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
546void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::UnscaleCosts() {
547 SCOPED_TIME_STAT(&stats_);
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);
552 }
553 cost_scaling_factor_ = 1;
554}
555
556template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
557bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Optimize() {
558 const CostValue kEpsilonMin = 1LL;
559 num_relabels_since_last_price_update_ = 0;
560 do {
561 // Avoid epsilon_ == 0.
562 epsilon_ = std::max(epsilon_ / alpha_, kEpsilonMin);
563 VLOG(3) << "Epsilon changed to: " << epsilon_;
564 if (!Refine()) return false;
565 } while (epsilon_ != 1LL);
566 return true;
567}
568
569template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
570void GenericMinCostFlow<Graph, ArcFlowType,
571 ArcScaledCostType>::SaturateAdmissibleArcs() {
572 SCOPED_TIME_STAT(&stats_);
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);
581 }
582 }
583
584 // We just saturated all the admissible arcs, so there are no arcs with a
585 // positive residual capacity that are incident to the current node.
586 // Moreover, during the course of the algorithm, if the residual capacity of
587 // such an arc becomes positive again, then the arc is still not admissible
588 // until we relabel the node (because the reverse arc was admissible for
589 // this to happen). In conclusion, the optimization below is correct.
590 first_admissible_arc_[node] = Graph::kNilArc;
591 }
592}
593
594template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
595void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::PushFlow(
596 FlowQuantity flow, ArcIndex arc) {
597 SCOPED_TIME_STAT(&stats_);
598 FastPushFlow(flow, arc, Tail(arc));
599}
600
601template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
602void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::FastPushFlow(
603 FlowQuantity flow, ArcIndex arc, NodeIndex tail) {
604 SCOPED_TIME_STAT(&stats_);
605 DCHECK_EQ(Tail(arc), tail);
606 DCHECK_GT(residual_arc_capacity_[arc], 0);
607 DCHECK_LE(flow, residual_arc_capacity_[arc]);
608 // Reduce the residual capacity on the arc by flow.
609 residual_arc_capacity_.Set(arc, residual_arc_capacity_[arc] - flow);
610 // Increase the residual capacity on the opposite arc by flow.
611 const ArcIndex opposite = Opposite(arc);
612 residual_arc_capacity_.Set(opposite, residual_arc_capacity_[opposite] + flow);
613 // Update the excesses at the tail and head of the arc.
614 node_excess_[tail] -= flow;
615 node_excess_[Head(arc)] += flow;
616}
617
618template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
619void GenericMinCostFlow<Graph, ArcFlowType,
620 ArcScaledCostType>::InitializeActiveNodeStack() {
621 SCOPED_TIME_STAT(&stats_);
622 DCHECK(active_nodes_.empty());
623 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
624 if (IsActive(node)) {
625 active_nodes_.push(node);
626 }
627 }
628}
629
630template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
631bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::UpdatePrices() {
632 SCOPED_TIME_STAT(&stats_);
633
634 // The algorithm works as follows. Start with a set of nodes S containing all
635 // the nodes with negative excess. Expand the set along reverse admissible
636 // arcs. If at the end, the complement of S contains at least one node with
637 // positive excess, relabel all the nodes in the complement of S by
638 // subtracting epsilon from their current potential. See the paper cited in
639 // the .h file.
640 //
641 // After this relabeling is done, the heuristic is reapplied by extending S as
642 // much as possible, relabeling the complement of S, and so on until there is
643 // no node with positive excess that is not in S. Note that this is not
644 // described in the paper.
645 //
646 // Note(user): The triggering mechanism of this UpdatePrices() is really
647 // important; if it is not done properly it may degrade performance!
648
649 // This represents the set S.
650 const NodeIndex num_nodes = graph_->num_nodes();
651 std::vector<NodeIndex> bfs_queue;
652 std::vector<bool> node_in_queue(num_nodes, false);
653
654 // This is used to update the potential of the nodes not in S.
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;
658
659 // Sum of the positive excesses out of S, used for early exit.
660 FlowQuantity remaining_excess = 0;
661
662 // First consider the nodes which have a negative excess.
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;
667
668 // This uses the fact that the sum of excesses is always 0.
669 remaining_excess -= node_excess_[node];
670 }
671 }
672
673 // All the nodes not yet in the bfs_queue will have their potential changed by
674 // +potential_delta (which becomes more and more negative at each pass). This
675 // update is applied when a node is pushed into the queue and at the end of
676 // the function for the nodes that are still unprocessed.
677 CostValue potential_delta = 0;
678
679 int queue_index = 0;
680 while (remaining_excess > 0) {
681 // Reverse BFS that expands S as much as possible in the reverse admissible
682 // graph. Once S cannot be expanded anymore, perform a relabeling on the
683 // nodes not in S but that can reach it in one arc and try to expand S
684 // again.
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();
689 it.Next()) {
690 const NodeIndex head = Head(it.Index());
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;
697 return false;
698 }
699 if (ReducedCost(opposite_arc) < 0) {
700 DCHECK(IsAdmissible(opposite_arc));
701
702 // TODO(user): Try to steal flow if node_excess_[head] > 0.
703 // An initial experiment didn't show a big speedup though.
704
705 remaining_excess -= node_excess_[head];
706 if (remaining_excess == 0) {
707 node_potential_[head] -= potential_delta;
708 break;
709 }
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);
715 }
716 } else {
717 // The opposite_arc is not admissible but is in the residual graph;
718 // this updates its min_non_admissible_potential.
719 node_potential_[head] -= potential_delta;
720 if (min_non_admissible_potential[head] == kMinCostValue) {
721 nodes_to_process.push_back(head);
722 }
723 min_non_admissible_potential[head] = std::max(
724 min_non_admissible_potential[head],
725 node_potential_[node] - scaled_arc_unit_cost_[opposite_arc]);
726 }
727 }
728 }
729 if (remaining_excess == 0) break;
730 }
731 if (remaining_excess == 0) break;
732
733 // Decrease by as much as possible instead of decreasing by epsilon.
734 // TODO(user): Is it worth the extra loop?
735 CostValue max_potential_diff = kMinCostValue;
736 for (int i = 0; i < nodes_to_process.size(); ++i) {
737 const NodeIndex node = nodes_to_process[i];
738 if (node_in_queue[node]) continue;
739 max_potential_diff =
740 std::max(max_potential_diff,
741 min_non_admissible_potential[node] - node_potential_[node]);
742 if (max_potential_diff == potential_delta) break;
743 }
744 DCHECK_LE(max_potential_diff, potential_delta);
745 potential_delta = max_potential_diff - epsilon_;
746
747 // Loop over nodes_to_process_ and for each node, apply the first of the
748 // rules below that match or leave it in the queue for later iteration:
749 // - Remove it if it is already in the queue.
750 // - If the node is connected to S by an admissible arc after it is
751 // relabeled by +potential_delta, add it to bfs_queue_ and remove it from
752 // nodes_to_process.
753 int index = 0;
754 for (int i = 0; i < nodes_to_process.size(); ++i) {
755 const NodeIndex node = nodes_to_process[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;
762 return false;
763 }
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];
769 continue;
770 }
771
772 // Keep the node for later iteration.
773 nodes_to_process[index] = node;
774 ++index;
775 }
776 nodes_to_process.resize(index);
777 }
778
779 // Update the potentials of the nodes not yet processed.
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;
786 return false;
787 }
788 first_admissible_arc_[node] = GetFirstOutgoingOrOppositeIncomingArc(node);
789 }
790 }
791 return true;
792}
793
794template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
795bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Refine() {
796 SCOPED_TIME_STAT(&stats_);
797 SaturateAdmissibleArcs();
798 InitializeActiveNodeStack();
799
800 const NodeIndex num_nodes = graph_->num_nodes();
801 while (!active_nodes_.empty()) {
802 // TODO(user): Experiment with different factors in front of num_nodes.
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;
807 }
808 }
809 const NodeIndex node = active_nodes_.top();
810 active_nodes_.pop();
811 DCHECK(IsActive(node));
812 if (!Discharge(node)) return false;
813 }
814 return true;
815}
816
817template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
818bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Discharge(
819 NodeIndex node) {
820 SCOPED_TIME_STAT(&stats_);
821 while (true) {
822 // The node is initially active, and we exit as soon as it becomes
823 // inactive.
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;
831
832 // We look ahead to see if this node can accept the flow or will need
833 // to be relabeled later in which case we do it now. Note that this will
834 // just be skipped for self-loop so it is fine.
835 const NodeIndex head = Head(arc);
836 if (node_excess_[head] >= 0 && !NodeHasAdmissibleArc(head)) {
837 // Relabel head and skip if the arc is not admissible anymore.
838 if (!Relabel(head)) return false;
839 if (!FastIsAdmissible(arc, tail_potential)) continue;
840 }
841
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);
849 }
850
851 if (node_excess_[node] == 0) {
852 // arc may still be admissible.
853 first_admissible_arc_[node] = arc;
854 return true;
855 }
856 }
857 if (!Relabel(node)) return false;
858 }
859}
860
861template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
864 SCOPED_TIME_STAT(&stats_);
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;
872 return true;
873 }
874 }
875 return false;
876}
877
878// Note that we only relabel if there is no leaving admissible arc.
879template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
880bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Relabel(
881 NodeIndex node) {
882 SCOPED_TIME_STAT(&stats_);
883 DCHECK(CheckRelabelPrecondition(node));
884 ++num_relabels_since_last_price_update_;
885
886 // By setting node_potential_[node] to the guaranteed_new_potential we are
887 // sure to keep epsilon-optimality of the pseudo-flow. Note that we could
888 // return right away with this value, but we prefer to check that this value
889 // will lead to at least one admissible arc, and if not, to decrease the
890 // potential as much as possible.
891 const CostValue guaranteed_new_potential = node_potential_[node] - epsilon_;
892 if (guaranteed_new_potential < overflow_threshold_) {
893 status_ = BAD_COST_RANGE;
894 return false;
895 }
896
897 // This will be updated to contain the minimum node potential for which
898 // the node has no admissible arc. We know that:
899 // - min_non_admissible_potential <= node_potential_[node]
900 // - We can set the new node potential to min_non_admissible_potential -
901 // epsilon_ and still keep the epsilon-optimality of the pseudo flow.
902 const CostValue kMinCostValue = std::numeric_limits<CostValue>::min();
903 CostValue min_non_admissible_potential = kMinCostValue;
904
905 // The following variables help setting the first_admissible_arc_[node] to a
906 // value different from GetFirstOutgoingOrOppositeIncomingArc(node) which
907 // avoids looking again at some arcs.
908 CostValue previous_min_non_admissible_potential = kMinCostValue;
909 ArcIndex first_arc = Graph::kNilArc;
910
911 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
912 it.Next()) {
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];
917 // We test this condition first as it is more probable that the first arc
918 // with residual capacity will be admissible if we reduce the node
919 // potential by epsilon.
920 if (min_non_admissible_potential_for_arc > guaranteed_new_potential) {
921 // We found an admissible arc for the guaranteed_new_potential. We
922 // stop right now instead of trying to compute the minimum possible
923 // new potential that keeps the epsilon-optimality of the pseudo flow.
924 node_potential_[node] = guaranteed_new_potential;
925 first_admissible_arc_[node] = arc;
926 return true;
927 }
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;
931 first_arc = arc;
932 }
933 }
934 }
935
936 // No residual arc leaves this node!
937 //
938 // TODO(user): This can be dealt with before the aglorithm start so that we
939 // do not need to test it here.
940 if (min_non_admissible_potential == kMinCostValue) {
941 if (node_excess_[node] != 0) {
942 // Note that this infeasibility detection is incomplete.
943 // Only max flow can detect that a min-cost flow problem is infeasible.
944 status_ = INFEASIBLE;
945 return false;
946 } else {
947 // This source saturates all its arcs, we can actually decrease the
948 // potential by as much as we want.
949 //
950 // TODO(user): Set it to a minimum value? but be careful of overflow.
951 node_potential_[node] = guaranteed_new_potential;
952 first_admissible_arc_[node] = Graph::kNilArc;
953 }
954 return true;
955 }
956
957 // We decrease the potential as much as possible, but we do not know the first
958 // admissible arc (most of the time). Keeping the
959 // previous_min_non_admissible_potential makes it faster by a few percent.
960 if (min_non_admissible_potential < overflow_threshold_) {
961 status_ = BAD_COST_RANGE;
962 return false;
963 }
964 const CostValue new_potential = min_non_admissible_potential - epsilon_;
965 if (new_potential < overflow_threshold_) {
966 status_ = BAD_COST_RANGE;
967 return false;
968 }
969 node_potential_[node] = new_potential;
970 if (previous_min_non_admissible_potential <= new_potential) {
971 first_admissible_arc_[node] = first_arc;
972 } else {
973 // We have no indication of what may be the first admissible arc.
974 first_admissible_arc_[node] = GetFirstOutgoingOrOppositeIncomingArc(node);
975 }
976 return true;
977}
978
979template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
980typename Graph::ArcIndex
981GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Opposite(
982 ArcIndex arc) const {
983 return graph_->OppositeArc(arc);
984}
985
986template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
987bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsArcValid(
988 ArcIndex arc) const {
989 return graph_->IsArcValid(arc);
990}
991
992template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
993bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsArcDirect(
994 ArcIndex arc) const {
995 DCHECK(IsArcValid(arc));
996 return arc >= 0;
997}
998
999// Explicit instantiations that can be used by a client.
1000//
1001// TODO(user): Move this code out of a .cc file and include it at the end of
1002// the header so it can work with any graph implementation?
1005template class GenericMinCostFlow<
1006 ::util::ReverseArcStaticGraph<uint16_t, int32_t>>;
1008 int64_t, int64_t>;
1009
1010// A more memory-efficient version for large graphs.
1011template class GenericMinCostFlow<
1012 ::util::ReverseArcStaticGraph<uint16_t, int32_t>,
1013 /*ArcFlowType=*/int16_t,
1014 /*ArcScaledCostType=*/int32_t>;
1015
1017 ArcIndex reserve_num_arcs) {
1018 if (reserve_num_nodes > 0) {
1019 node_supply_.reserve(reserve_num_nodes);
1020 }
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);
1028 }
1029}
1030
1032 ResizeNodeVectors(node);
1033 node_supply_[node] = supply;
1034}
1035
1037 NodeIndex tail, NodeIndex head, FlowQuantity capacity,
1038 CostValue unit_cost) {
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);
1045 return arc;
1046}
1047
1049 arc_capacity_[arc] = capacity;
1050}
1051
1052SimpleMinCostFlow::ArcIndex SimpleMinCostFlow::PermutedArc(ArcIndex arc) {
1053 return arc < arc_permutation_.size() ? arc_permutation_[arc] : arc;
1054}
1055
1056SimpleMinCostFlow::Status SimpleMinCostFlow::SolveWithPossibleAdjustment(
1057 SupplyAdjustment adjustment) {
1058 optimal_cost_ = 0;
1059 maximum_flow_ = 0;
1060 arc_flow_.clear();
1061 const NodeIndex num_nodes = node_supply_.size();
1062 const ArcIndex num_arcs = arc_capacity_.size();
1063 if (num_nodes == 0) return OPTIMAL;
1064
1065 int supply_node_count = 0, demand_node_count = 0;
1066 FlowQuantity total_supply = 0, total_demand = 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];
1074 }
1075 }
1076 if (adjustment == DONT_ADJUST && total_supply != total_demand) {
1077 return UNBALANCED;
1078 }
1079
1080 // Feasibility checking, and possible supply/demand adjustment, is done by:
1081 // 1. Creating a new source and sink node.
1082 // 2. Taking all nodes that have a non-zero supply or demand and
1083 // connecting them to the source or sink respectively. The arc thus
1084 // added has a capacity of the supply or demand.
1085 // 3. Computing the max flow between the new source and sink.
1086 // 4. If adjustment isn't being done, checking that the max flow is equal
1087 // to the total supply/demand (and returning INFEASIBLE if it isn't).
1088 // 5. Running min-cost max-flow on this augmented graph, using the max
1089 // flow computed in step 3 as the supply of the source and demand of
1090 // the sink.
1091 const ArcIndex augmented_num_arcs =
1092 num_arcs + supply_node_count + demand_node_count;
1093 const NodeIndex source = num_nodes;
1094 const NodeIndex sink = num_nodes + 1;
1095 const NodeIndex augmented_num_nodes = num_nodes + 2;
1096
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]);
1100 }
1101
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);
1107 }
1108 }
1109
1110 graph.Build(&arc_permutation_);
1111
1112 {
1113 GenericMaxFlow<Graph> max_flow(&graph, source, sink);
1114 ArcIndex arc;
1115 for (arc = 0; arc < num_arcs; ++arc) {
1116 max_flow.SetArcCapacity(PermutedArc(arc), arc_capacity_[arc]);
1117 }
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]));
1121 ++arc;
1122 }
1123 }
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()) {
1129 return NOT_SOLVED;
1131 LOG(ERROR)
1132 << "Max flow failed but claimed to have an optimal solution";
1133 ABSL_FALLTHROUGH_INTENDED;
1134 default:
1135 return BAD_RESULT;
1136 }
1137 }
1138 maximum_flow_ = max_flow.GetOptimalFlow();
1139 }
1140
1141 if (adjustment == DONT_ADJUST && maximum_flow_ != total_supply) {
1142 return INFEASIBLE;
1143 }
1144
1145 GenericMinCostFlow<Graph> min_cost_flow(&graph);
1146 ArcIndex arc;
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]);
1151 }
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);
1157 ++arc;
1158 }
1159 }
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_);
1164
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));
1170 }
1171 }
1172
1173 return min_cost_flow.status();
1174}
1175
1177 return optimal_cost_;
1178}
1179
1181 return maximum_flow_;
1182}
1183
1185 return arc_flow_[arc];
1186}
1187
1188SimpleMinCostFlow::SimpleMinCostFlow::NodeIndex SimpleMinCostFlow::NumNodes()
1189 const {
1190 return node_supply_.size();
1191}
1192
1194 return arc_tail_.size();
1195}
1196
1198 return arc_tail_[arc];
1199}
1200
1202 return arc_head_[arc];
1203}
1204
1206 ArcIndex arc) const {
1207 return arc_capacity_[arc];
1208}
1209
1211 return arc_cost_[arc];
1212}
1213
1215 NodeIndex node) const {
1216 return node_supply_[node];
1217}
1218
1219void SimpleMinCostFlow::ResizeNodeVectors(NodeIndex node) {
1220 if (node < node_supply_.size()) return;
1221 node_supply_.resize(node + 1);
1222}
1223
1224} // namespace operations_research
FlowQuantity Flow(ArcIndex arc) const
void SetArcUnitCost(ArcIndex arc, ArcScaledCostType unit_cost)
Sets the unit cost for the given arc.
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
FlowQuantity Capacity(ArcIndex arc) const
We use the equations given in the comment of residual_arc_capacity_.
static T Abs(const T x)
Definition mathutil.h:95
NodeIndex Tail(ArcIndex arc) 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
ArcIndex AddArcWithCapacityAndUnitCost(NodeIndex tail, NodeIndex head, FlowQuantity capacity, CostValue unit_cost)
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
Definition graph.h:293
static constexpr bool kHasNegativeReverseArcs
Definition graph.h:215
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)
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)
Definition stats.h:414
#define SCOPED_TIME_STAT(stats)
Definition stats.h:421