OR-Tools  8.2
min_cost_flow.cc
Go to the documentation of this file.
1// Copyright 2010-2018 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 <limits>
19
20#include "absl/strings/str_format.h"
23#include "ortools/graph/graph.h"
26
27// TODO(user): Remove these flags and expose the parameters in the API.
28// New clients, please do not use these flags!
29ABSL_FLAG(int64, min_cost_flow_alpha, 5,
30 "Divide factor for epsilon at each refine step.");
31ABSL_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.");
35ABSL_FLAG(bool, min_cost_flow_check_balance, true,
36 "Check that the sum of supplies is equal to the sum of demands.");
37ABSL_FLAG(bool, min_cost_flow_check_costs, true,
38 "Check that the magnitude of the costs will not exceed the "
39 "precision of the machine when scaled (multiplied) by the number "
40 "of nodes");
41ABSL_FLAG(bool, min_cost_flow_check_result, true,
42 "Check that the result is valid.");
43
44namespace operations_research {
45
46template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
48 const Graph* graph)
49 : graph_(graph),
50 node_excess_(),
51 node_potential_(),
52 residual_arc_capacity_(),
53 first_admissible_arc_(),
54 active_nodes_(),
55 epsilon_(0),
56 alpha_(absl::GetFlag(FLAGS_min_cost_flow_alpha)),
57 cost_scaling_factor_(1),
58 scaled_arc_unit_cost_(),
59 total_flow_cost_(0),
60 status_(NOT_SOLVED),
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)) {
67 const NodeIndex max_num_nodes = Graphs<Graph>::NodeReservation(*graph_);
68 if (max_num_nodes > 0) {
69 node_excess_.Reserve(0, max_num_nodes - 1);
70 node_excess_.SetAll(0);
71 node_potential_.Reserve(0, max_num_nodes - 1);
72 node_potential_.SetAll(0);
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);
79 }
80 const ArcIndex max_num_arcs = Graphs<Graph>::ArcReservation(*graph_);
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);
86 }
87}
88
89template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
91 NodeIndex node, FlowQuantity supply) {
92 DCHECK(graph_->IsNodeValid(node));
93 node_excess_.Set(node, supply);
94 initial_node_excess_.Set(node, supply);
95 status_ = NOT_SOLVED;
96 feasibility_checked_ = false;
97}
98
99template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
101 ArcIndex arc, ArcScaledCostType unit_cost) {
102 DCHECK(IsArcDirect(arc));
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;
107}
108
109template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
111 ArcIndex arc, ArcFlowType new_capacity) {
112 DCHECK_LE(0, new_capacity);
113 DCHECK(IsArcDirect(arc));
114 const FlowQuantity free_capacity = residual_arc_capacity_[arc];
115 const FlowQuantity capacity_delta = new_capacity - Capacity(arc);
116 if (capacity_delta == 0) {
117 return; // Nothing to do.
118 }
119 status_ = NOT_SOLVED;
120 feasibility_checked_ = false;
121 const FlowQuantity new_availability = free_capacity + capacity_delta;
122 if (new_availability >= 0) {
123 // The above condition is true when one of two following holds:
124 // 1/ (capacity_delta > 0), meaning we are increasing the capacity
125 // 2/ (capacity_delta < 0 && free_capacity + capacity_delta >= 0)
126 // meaning we are reducing the capacity, but that the capacity
127 // reduction is not larger than the free capacity.
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]);
132 } else {
133 // We have to reduce the flow on the arc, and update the excesses
134 // accordingly.
135 const FlowQuantity flow = residual_arc_capacity_[Opposite(arc)];
136 const FlowQuantity flow_excess = flow - new_capacity;
137 residual_arc_capacity_.Set(arc, 0);
138 residual_arc_capacity_.Set(Opposite(arc), new_capacity);
139 const NodeIndex tail = Tail(arc);
140 node_excess_.Set(tail, node_excess_[tail] + flow_excess);
141 const NodeIndex head = Head(arc);
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)]);
145 }
146}
147
148template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
150 ArcIndex arc, ArcFlowType new_flow) {
151 DCHECK(IsArcValid(arc));
152 const FlowQuantity capacity = Capacity(arc);
153 DCHECK_GE(capacity, 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;
158}
159
160template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
161bool GenericMinCostFlow<Graph, ArcFlowType,
162 ArcScaledCostType>::CheckInputConsistency() const {
163 FlowQuantity total_supply = 0;
164 uint64 max_capacity = 0; // uint64 because it is positive and will be used
165 // to check against FlowQuantity overflows.
166 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
167 const uint64 capacity = static_cast<uint64>(residual_arc_capacity_[arc]);
168 max_capacity = std::max(capacity, max_capacity);
169 }
170 uint64 total_flow = 0; // uint64 for the same reason as max_capacity.
171 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
172 const FlowQuantity excess = node_excess_[node];
173 total_supply += excess;
174 if (excess > 0) {
175 total_flow += excess;
177 max_capacity + total_flow) {
178 LOG(DFATAL) << "Input consistency error: max capacity + flow exceed "
179 << "precision";
180 return false;
181 }
182 }
183 }
184 if (total_supply != 0) {
185 LOG(DFATAL) << "Input consistency error: unbalanced problem";
186 return false;
187 }
188 return true;
189}
190
191template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
192bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::CheckResult()
193 const {
194 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
195 if (node_excess_[node] != 0) {
196 LOG(DFATAL) << "node_excess_[" << node << "] != 0";
197 return false;
198 }
199 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
200 it.Next()) {
201 const ArcIndex arc = it.Index();
202 bool ok = true;
203 if (residual_arc_capacity_[arc] < 0) {
204 LOG(DFATAL) << "residual_arc_capacity_[" << arc << "] < 0";
205 ok = false;
206 }
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_ << ").";
211 ok = false;
212 }
213 if (!ok) {
214 LOG(DFATAL) << DebugString("CheckResult ", arc);
215 return false;
216 }
217 }
218 }
219 return true;
220}
221
222template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
223bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::CheckCostRange()
224 const {
225 CostValue min_cost_magnitude = std::numeric_limits<CostValue>::max();
226 CostValue max_cost_magnitude = 0;
227 // Traverse the initial arcs of the graph:
228 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
229 const CostValue cost_magnitude = MathUtil::Abs(scaled_arc_unit_cost_[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);
233 }
234 }
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.";
242 return false;
243 }
244#endif
245 return true;
246}
247
248template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
249bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::
250 CheckRelabelPrecondition(NodeIndex node) const {
251 // Note that the classical Relabel precondition assumes IsActive(node), i.e.,
252 // the node_excess_[node] > 0. However, to implement the Push Look-Ahead
253 // heuristic, we can relax this condition as explained in the section 4.3 of
254 // the article "An Efficient Implementation of a Scaling Minimum-Cost Flow
255 // Algorithm", A.V. Goldberg, Journal of Algorithms 22(1), January 1997, pp.
256 // 1-29.
257 DCHECK_GE(node_excess_[node], 0);
258 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
259 it.Next()) {
260 const ArcIndex arc = it.Index();
261 DCHECK(!IsAdmissible(arc)) << DebugString("CheckRelabelPrecondition:", arc);
262 }
263 return true;
264}
265
266template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
267std::string
268GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::DebugString(
269 const std::string& context, ArcIndex arc) const {
270 const NodeIndex tail = Tail(arc);
271 const NodeIndex head = Head(arc);
272 // Reduced cost is computed directly without calling ReducedCost to avoid
273 // recursive calls between ReducedCost and DebugString in case a DCHECK in
274 // ReducedCost fails.
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, ",
284 context, arc, tail, head, Capacity(arc),
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]),
288 reduced_cost);
289}
290
291template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
293 CheckFeasibility(std::vector<NodeIndex>* const infeasible_supply_node,
294 std::vector<NodeIndex>* const infeasible_demand_node) {
295 SCOPED_TIME_STAT(&stats_);
296 // Create a new graph, which is a copy of graph_, with the following
297 // modifications:
298 // Two nodes are added: a source and a sink.
299 // The source is linked to each supply node (whose supply > 0) by an arc whose
300 // capacity is equal to the supply at the supply node.
301 // The sink is linked to each demand node (whose supply < 0) by an arc whose
302 // capacity is the demand (-supply) at the demand node.
303 // There are no supplies or demands or costs in the graph, as we will run
304 // max-flow.
305 // TODO(user): make it possible to share a graph by MaxFlow and MinCostFlow.
306 // For this it is necessary to make StarGraph resizable.
307 feasibility_checked_ = false;
308 ArcIndex num_extra_arcs = 0;
309 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
310 if (initial_node_excess_[node] != 0) {
311 ++num_extra_arcs;
312 }
313 }
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);
320 checker.SetCheckInput(false);
321 checker.SetCheckResult(false);
322 // Copy graph_ to checker_graph.
323 for (ArcIndex arc = 0; arc < graph_->num_arcs(); ++arc) {
324 const ArcIndex new_arc =
325 checker_graph.AddArc(graph_->Tail(arc), graph_->Head(arc));
326 DCHECK_EQ(arc, new_arc);
327 checker.SetArcCapacity(new_arc, Capacity(arc));
328 }
329 FlowQuantity total_demand = 0;
330 FlowQuantity total_supply = 0;
331 // Create the source-to-supply node arcs and the demand-node-to-sink arcs.
332 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
333 const FlowQuantity supply = initial_node_excess_[node];
334 if (supply > 0) {
335 const ArcIndex new_arc = checker_graph.AddArc(source, node);
336 checker.SetArcCapacity(new_arc, supply);
337 total_supply += supply;
338 } else if (supply < 0) {
339 const ArcIndex new_arc = checker_graph.AddArc(node, sink);
340 checker.SetArcCapacity(new_arc, -supply);
341 total_demand -= supply;
342 }
343 }
344 if (total_supply != total_demand) {
345 LOG(DFATAL) << "total_supply(" << total_supply << ") != total_demand("
346 << total_demand << ").";
347 return false;
348 }
349 if (!checker.Solve()) {
350 LOG(DFATAL) << "Max flow could not be computed.";
351 return false;
352 }
353 const FlowQuantity optimal_max_flow = checker.GetOptimalFlow();
354 feasible_node_excess_.SetAll(0);
355 for (StarGraph::OutgoingArcIterator it(checker_graph, source); it.Ok();
356 it.Next()) {
357 const ArcIndex arc = it.Index();
358 const NodeIndex node = checker_graph.Head(arc);
359 const FlowQuantity flow = checker.Flow(arc);
360 feasible_node_excess_.Set(node, flow);
361 if (infeasible_supply_node != nullptr) {
362 infeasible_supply_node->push_back(node);
364 }
365 for (StarGraph::IncomingArcIterator it(checker_graph, sink); it.Ok();
366 it.Next()) {
367 const ArcIndex arc = it.Index();
368 const NodeIndex node = checker_graph.Tail(arc);
369 const FlowQuantity flow = checker.Flow(arc);
370 feasible_node_excess_.Set(node, -flow);
371 if (infeasible_demand_node != nullptr) {
372 infeasible_demand_node->push_back(node);
374 }
375 feasibility_checked_ = true;
376 return optimal_max_flow == total_supply;
377}
378
379template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
381 if (!feasibility_checked_) {
382 return false;
383 }
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);
389 return true;
390}
392template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
394 ArcIndex arc) const {
395 if (IsArcDirect(arc)) {
396 return residual_arc_capacity_[Opposite(arc)];
397 } else {
398 return -residual_arc_capacity_[arc];
399 }
400}
402// We use the equations given in the comment of residual_arc_capacity_.
403template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
406 ArcIndex arc) const {
407 if (IsArcDirect(arc)) {
408 return residual_arc_capacity_[arc] + residual_arc_capacity_[Opposite(arc)];
409 } else {
410 return 0;
411 }
412}
413
414template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
416 ArcIndex arc) const {
417 DCHECK(IsArcValid(arc));
418 DCHECK_EQ(uint64{1}, cost_scaling_factor_);
419 return scaled_arc_unit_cost_[arc];
420}
421
422template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
424 NodeIndex node) const {
425 DCHECK(graph_->IsNodeValid(node));
426 return node_excess_[node];
427}
428
429template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
432 NodeIndex node) const {
433 return initial_node_excess_[node];
434}
435
436template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
439 NodeIndex node) const {
440 return feasible_node_excess_[node];
441}
442
443template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
445 ArcIndex arc) const {
446 return FastIsAdmissible(arc, node_potential_[Tail(arc)]);
447}
448
449template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
450bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::
451 FastIsAdmissible(ArcIndex arc, CostValue tail_potential) const {
452 DCHECK_EQ(node_potential_[Tail(arc)], tail_potential);
453 return residual_arc_capacity_[arc] > 0 &&
454 FastReducedCost(arc, tail_potential) < 0;
455}
456
457template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
458bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsActive(
459 NodeIndex node) const {
460 return node_excess_[node] > 0;
461}
462
463template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
465GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::ReducedCost(
466 ArcIndex arc) const {
467 return FastReducedCost(arc, node_potential_[Tail(arc)]);
468}
469
470template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
472GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::FastReducedCost(
473 ArcIndex arc, CostValue tail_potential) const {
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)];
481}
482
483template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
485GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::
486 GetFirstOutgoingOrOppositeIncomingArc(NodeIndex node) const {
487 OutgoingOrOppositeIncomingArcIterator arc_it(*graph_, node);
488 return arc_it.Index();
489}
490
491template <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;
497 return false;
498 }
499 if (absl::GetFlag(FLAGS_min_cost_flow_check_costs) && !CheckCostRange()) {
500 status_ = BAD_COST_RANGE;
501 return false;
502 }
503 if (check_feasibility_ && !CheckFeasibility(nullptr, nullptr)) {
504 status_ = INFEASIBLE;
505 return false;
506 }
507 node_potential_.SetAll(0);
508 ResetFirstAdmissibleArcs();
509 ScaleCosts();
510 Optimize();
511 if (absl::GetFlag(FLAGS_min_cost_flow_check_result) && !CheckResult()) {
512 status_ = BAD_RESULT;
513 UnscaleCosts();
514 return false;
515 }
516 UnscaleCosts();
517 if (status_ != OPTIMAL) {
518 LOG(DFATAL) << "Status != OPTIMAL";
519 total_flow_cost_ = 0;
520 return false;
521 }
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;
526 }
527 status_ = OPTIMAL;
528 IF_STATS_ENABLED(VLOG(1) << stats_.StatString());
529 return true;
530}
531
532template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
533void GenericMinCostFlow<Graph, ArcFlowType,
534 ArcScaledCostType>::ResetFirstAdmissibleArcs() {
535 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
536 first_admissible_arc_.Set(node,
537 GetFirstOutgoingOrOppositeIncomingArc(node));
538 }
539}
540
541template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
542void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::ScaleCosts() {
543 SCOPED_TIME_STAT(&stats_);
544 cost_scaling_factor_ = graph_->num_nodes() + 1;
545 epsilon_ = 1LL;
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);
552 epsilon_ = std::max(epsilon_, MathUtil::Abs(cost));
553 }
554 VLOG(3) << "Initial epsilon = " << epsilon_;
555 VLOG(3) << "Cost scaling factor = " << cost_scaling_factor_;
556}
557
558template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
559void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::UnscaleCosts() {
560 SCOPED_TIME_STAT(&stats_);
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);
565 }
566 cost_scaling_factor_ = 1;
567}
568
569template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
570void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Optimize() {
571 const CostValue kEpsilonMin = 1LL;
572 num_relabels_since_last_price_update_ = 0;
573 do {
574 // Avoid epsilon_ == 0.
575 epsilon_ = std::max(epsilon_ / alpha_, kEpsilonMin);
576 VLOG(3) << "Epsilon changed to: " << epsilon_;
577 Refine();
578 } while (epsilon_ != 1LL && status_ != INFEASIBLE);
579 if (status_ == NOT_SOLVED) {
580 status_ = OPTIMAL;
581 }
582}
583
584template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
585void GenericMinCostFlow<Graph, ArcFlowType,
586 ArcScaledCostType>::SaturateAdmissibleArcs() {
587 SCOPED_TIME_STAT(&stats_);
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()) {
593 const ArcIndex arc = it.Index();
594 if (FastIsAdmissible(arc, tail_potential)) {
595 FastPushFlow(residual_arc_capacity_[arc], arc, node);
596 }
597 }
598
599 // We just saturated all the admissible arcs, so there are no arcs with a
600 // positive residual capacity that are incident to the current node.
601 // Moreover, during the course of the algorithm, if the residual capacity of
602 // such an arc becomes positive again, then the arc is still not admissible
603 // until we relabel the node (because the reverse arc was admissible for
604 // this to happen). In conclusion, the optimization below is correct.
605 first_admissible_arc_[node] = Graph::kNilArc;
606 }
607}
608
609template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
610void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::PushFlow(
611 FlowQuantity flow, ArcIndex arc) {
612 SCOPED_TIME_STAT(&stats_);
613 FastPushFlow(flow, arc, Tail(arc));
614}
615
616template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
617void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::FastPushFlow(
618 FlowQuantity flow, ArcIndex arc, NodeIndex tail) {
619 SCOPED_TIME_STAT(&stats_);
620 DCHECK_EQ(Tail(arc), tail);
621 DCHECK_GT(residual_arc_capacity_[arc], 0);
622 DCHECK_LE(flow, residual_arc_capacity_[arc]);
623 // Reduce the residual capacity on the arc by flow.
624 residual_arc_capacity_.Set(arc, residual_arc_capacity_[arc] - flow);
625 // Increase the residual capacity on the opposite arc by flow.
626 const ArcIndex opposite = Opposite(arc);
627 residual_arc_capacity_.Set(opposite, residual_arc_capacity_[opposite] + flow);
628 // Update the excesses at the tail and head of the arc.
629 node_excess_.Set(tail, node_excess_[tail] - flow);
630 const NodeIndex head = Head(arc);
631 node_excess_.Set(head, node_excess_[head] + flow);
632}
633
634template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
635void GenericMinCostFlow<Graph, ArcFlowType,
636 ArcScaledCostType>::InitializeActiveNodeStack() {
637 SCOPED_TIME_STAT(&stats_);
638 DCHECK(active_nodes_.empty());
639 for (NodeIndex node = 0; node < graph_->num_nodes(); ++node) {
640 if (IsActive(node)) {
641 active_nodes_.push(node);
642 }
643 }
644}
645
646template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
647void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::UpdatePrices() {
648 SCOPED_TIME_STAT(&stats_);
649
650 // The algorithm works as follows. Start with a set of nodes S containing all
651 // the nodes with negative excess. Expand the set along reverse admissible
652 // arcs. If at the end, the complement of S contains at least one node with
653 // positive excess, relabel all the nodes in the complement of S by
654 // subtracting epsilon from their current potential. See the paper cited in
655 // the .h file.
656 //
657 // After this relabeling is done, the heuristic is reapplied by extending S as
658 // much as possible, relabeling the complement of S, and so on until there is
659 // no node with positive excess that is not in S. Note that this is not
660 // described in the paper.
661 //
662 // Note(user): The triggering mechanism of this UpdatePrices() is really
663 // important; if it is not done properly it may degrade performance!
664
665 // This represents the set S.
666 const NodeIndex num_nodes = graph_->num_nodes();
667 std::vector<NodeIndex> bfs_queue;
668 std::vector<bool> node_in_queue(num_nodes, false);
669
670 // This is used to update the potential of the nodes not in S.
671 const CostValue kMinCostValue = std::numeric_limits<CostValue>::min();
672 std::vector<CostValue> min_non_admissible_potential(num_nodes, kMinCostValue);
673 std::vector<NodeIndex> nodes_to_process;
674
675 // Sum of the positive excesses out of S, used for early exit.
676 FlowQuantity remaining_excess = 0;
677
678 // First consider the nodes which have a negative excess.
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;
683
684 // This uses the fact that the sum of excesses is always 0.
685 remaining_excess -= node_excess_[node];
686 }
687 }
688
689 // All the nodes not yet in the bfs_queue will have their potential changed by
690 // +potential_delta (which becomes more and more negative at each pass). This
691 // update is applied when a node is pushed into the queue and at the end of
692 // the function for the nodes that are still unprocessed.
693 CostValue potential_delta = 0;
694
695 int queue_index = 0;
696 while (remaining_excess > 0) {
697 // Reverse BFS that expands S as much as possible in the reverse admissible
698 // graph. Once S cannot be expanded anymore, perform a relabeling on the
699 // nodes not in S but that can reach it in one arc and try to expand S
700 // again.
701 for (; queue_index < bfs_queue.size(); ++queue_index) {
702 DCHECK_GE(num_nodes, bfs_queue.size());
703 const NodeIndex node = bfs_queue[queue_index];
704 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
705 it.Next()) {
706 const NodeIndex head = Head(it.Index());
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));
713
714 // TODO(user): Try to steal flow if node_excess_[head] > 0.
715 // An initial experiment didn't show a big speedup though.
716
717 remaining_excess -= node_excess_[head];
718 if (remaining_excess == 0) {
719 node_potential_[head] -= potential_delta;
720 break;
721 }
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);
727 }
728 } else {
729 // The opposite_arc is not admissible but is in the residual graph;
730 // this updates its min_non_admissible_potential.
731 node_potential_[head] -= potential_delta;
732 if (min_non_admissible_potential[head] == kMinCostValue) {
733 nodes_to_process.push_back(head);
734 }
735 min_non_admissible_potential[head] = std::max(
736 min_non_admissible_potential[head],
737 node_potential_[node] - scaled_arc_unit_cost_[opposite_arc]);
738 }
739 }
740 }
741 if (remaining_excess == 0) break;
742 }
743 if (remaining_excess == 0) break;
744
745 // Decrease by as much as possible instead of decreasing by epsilon.
746 // TODO(user): Is it worth the extra loop?
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;
751 max_potential_diff =
752 std::max(max_potential_diff,
753 min_non_admissible_potential[node] - node_potential_[node]);
754 if (max_potential_diff == potential_delta) break;
755 }
756 DCHECK_LE(max_potential_diff, potential_delta);
757 potential_delta = max_potential_diff - epsilon_;
758
759 // Loop over nodes_to_process_ and for each node, apply the first of the
760 // rules below that match or leave it in the queue for later iteration:
761 // - Remove it if it is already in the queue.
762 // - If the node is connected to S by an admissible arc after it is
763 // relabeled by +potential_delta, add it to bfs_queue_ and remove it from
764 // nodes_to_process.
765 int index = 0;
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];
777 continue;
778 }
779
780 // Keep the node for later iteration.
781 nodes_to_process[index] = node;
782 ++index;
783 }
784 nodes_to_process.resize(index);
785 }
786
787 // Update the potentials of the nodes not yet processed.
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);
793 }
794 }
795}
796
797template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
798void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Refine() {
799 SCOPED_TIME_STAT(&stats_);
800 SaturateAdmissibleArcs();
801 InitializeActiveNodeStack();
802
803 const NodeIndex num_nodes = graph_->num_nodes();
804 while (status_ != INFEASIBLE && !active_nodes_.empty()) {
805 // TODO(user): Experiment with different factors in front of num_nodes.
806 if (num_relabels_since_last_price_update_ >= num_nodes) {
807 num_relabels_since_last_price_update_ = 0;
808 if (use_price_update_) {
809 UpdatePrices();
810 }
811 }
812 const NodeIndex node = active_nodes_.top();
813 active_nodes_.pop();
814 DCHECK(IsActive(node));
815 Discharge(node);
816 }
817}
818
819template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
820void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Discharge(
821 NodeIndex node) {
822 SCOPED_TIME_STAT(&stats_);
823 do {
824 // The node is initially active, and we exit as soon as it becomes
825 // inactive.
826 DCHECK(IsActive(node));
827 const CostValue tail_potential = node_potential_[node];
828 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node,
829 first_admissible_arc_[node]);
830 it.Ok(); it.Next()) {
831 const ArcIndex arc = it.Index();
832 if (FastIsAdmissible(arc, tail_potential)) {
833 const NodeIndex head = Head(arc);
834 if (!LookAhead(arc, tail_potential, head)) continue;
835 const bool head_active_before_push = IsActive(head);
836 const FlowQuantity delta =
837 std::min(node_excess_[node],
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);
842 }
843 if (node_excess_[node] == 0) {
844 // arc may still be admissible.
845 first_admissible_arc_.Set(node, arc);
846 return;
847 }
848 }
849 }
850 Relabel(node);
851 } while (status_ != INFEASIBLE);
852}
853
854template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
855bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::LookAhead(
856 ArcIndex in_arc, CostValue in_tail_potential, NodeIndex node) {
857 SCOPED_TIME_STAT(&stats_);
858 DCHECK_EQ(Head(in_arc), node);
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()) {
865 const ArcIndex arc = it.Index();
866 if (FastIsAdmissible(arc, tail_potential)) {
867 first_admissible_arc_.Set(node, arc);
868 return true;
869 }
870 }
871
872 // The node we looked ahead has no admissible arc at its current potential.
873 // We relabel it and return true if the original arc is still admissible.
874 Relabel(node);
875 return FastIsAdmissible(in_arc, in_tail_potential);
876}
877
878template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
879void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Relabel(
880 NodeIndex node) {
881 SCOPED_TIME_STAT(&stats_);
882 DCHECK(CheckRelabelPrecondition(node));
883 ++num_relabels_since_last_price_update_;
884
885 // By setting node_potential_[node] to the guaranteed_new_potential we are
886 // sure to keep epsilon-optimality of the pseudo-flow. Note that we could
887 // return right away with this value, but we prefer to check that this value
888 // will lead to at least one admissible arc, and if not, to decrease the
889 // potential as much as possible.
890 const CostValue guaranteed_new_potential = node_potential_[node] - epsilon_;
891
892 // This will be updated to contain the minimum node potential for which
893 // the node has no admissible arc. We know that:
894 // - min_non_admissible_potential <= node_potential_[node]
895 // - We can set the new node potential to min_non_admissible_potential -
896 // epsilon_ and still keep the epsilon-optimality of the pseudo flow.
897 const CostValue kMinCostValue = std::numeric_limits<CostValue>::min();
898 CostValue min_non_admissible_potential = kMinCostValue;
899
900 // The following variables help setting the first_admissible_arc_[node] to a
901 // value different from GetFirstOutgoingOrOppositeIncomingArc(node) which
902 // avoids looking again at some arcs.
903 CostValue previous_min_non_admissible_potential = kMinCostValue;
904 ArcIndex first_arc = Graph::kNilArc;
905
906 for (OutgoingOrOppositeIncomingArcIterator it(*graph_, node); it.Ok();
907 it.Next()) {
908 const ArcIndex arc = it.Index();
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) {
914 // We found an admissible arc for the guaranteed_new_potential. We
915 // stop right now instead of trying to compute the minimum possible
916 // new potential that keeps the epsilon-optimality of the pseudo flow.
917 node_potential_.Set(node, guaranteed_new_potential);
918 first_admissible_arc_.Set(node, arc);
919 return;
920 }
921 previous_min_non_admissible_potential = min_non_admissible_potential;
922 min_non_admissible_potential = min_non_admissible_potential_for_arc;
923 first_arc = arc;
924 }
925 }
926 }
927
928 // No admissible arc leaves this node!
929 if (min_non_admissible_potential == kMinCostValue) {
930 if (node_excess_[node] != 0) {
931 // Note that this infeasibility detection is incomplete.
932 // Only max flow can detect that a min-cost flow problem is infeasible.
933 status_ = INFEASIBLE;
934 LOG(ERROR) << "Infeasible problem.";
935 } else {
936 // This source saturates all its arcs, we can actually decrease the
937 // potential by as much as we want.
938 // TODO(user): Set it to a minimum value, but be careful of overflow.
939 node_potential_.Set(node, guaranteed_new_potential);
940 first_admissible_arc_.Set(node,
941 GetFirstOutgoingOrOppositeIncomingArc(node));
942 }
943 return;
944 }
945
946 // We decrease the potential as much as possible, but we do not know the first
947 // admissible arc (most of the time). Keeping the
948 // previous_min_non_admissible_potential makes it faster by a few percent.
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);
953 } else {
954 // We have no indication of what may be the first admissible arc.
955 first_admissible_arc_.Set(node,
956 GetFirstOutgoingOrOppositeIncomingArc(node));
957 }
958}
959
960template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
961typename Graph::ArcIndex
962GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Opposite(
963 ArcIndex arc) const {
964 return Graphs<Graph>::OppositeArc(*graph_, arc);
965}
966
967template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
968bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsArcValid(
969 ArcIndex arc) const {
970 return Graphs<Graph>::IsArcValid(*graph_, arc);
971}
972
973template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
974bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsArcDirect(
975 ArcIndex arc) const {
976 DCHECK(IsArcValid(arc));
977 return arc >= 0;
978}
979
980// Explicit instantiations that can be used by a client.
981//
982// TODO(user): Move this code out of a .cc file and include it at the end of
983// the header so it can work with any graph implementation?
984template class GenericMinCostFlow<StarGraph>;
985template class GenericMinCostFlow<::util::ReverseArcListGraph<>>;
986template class GenericMinCostFlow<::util::ReverseArcStaticGraph<>>;
987template class GenericMinCostFlow<::util::ReverseArcMixedGraph<>>;
988template class GenericMinCostFlow<::util::ReverseArcStaticGraph<uint16, int32>>;
989
990// A more memory-efficient version for large graphs.
991template class GenericMinCostFlow<::util::ReverseArcStaticGraph<uint16, int32>,
992 /*ArcFlowType=*/int16,
993 /*ArcScaledCostType=*/int32>;
994
996 ArcIndex reserve_num_arcs) {
997 if (reserve_num_nodes > 0) {
998 node_supply_.reserve(reserve_num_nodes);
999 }
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);
1007 }
1008}
1009
1011 ResizeNodeVectors(node);
1012 node_supply_[node] = supply;
1013}
1014
1018 CostValue unit_cost) {
1019 ResizeNodeVectors(std::max(tail, head));
1020 const ArcIndex arc = arc_tail_.size();
1021 arc_tail_.push_back(tail);
1022 arc_head_.push_back(head);
1023 arc_capacity_.push_back(capacity);
1024 arc_cost_.push_back(unit_cost);
1025 return arc;
1026}
1027
1028ArcIndex SimpleMinCostFlow::PermutedArc(ArcIndex arc) {
1029 return arc < arc_permutation_.size() ? arc_permutation_[arc] : arc;
1030}
1031
1032SimpleMinCostFlow::Status SimpleMinCostFlow::SolveWithPossibleAdjustment(
1033 SupplyAdjustment adjustment) {
1034 optimal_cost_ = 0;
1035 maximum_flow_ = 0;
1036 arc_flow_.clear();
1037 const NodeIndex num_nodes = node_supply_.size();
1038 const ArcIndex num_arcs = arc_capacity_.size();
1039 if (num_nodes == 0) return OPTIMAL;
1040
1041 int supply_node_count = 0, demand_node_count = 0;
1042 FlowQuantity total_supply = 0, total_demand = 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];
1050 }
1051 }
1052 if (adjustment == DONT_ADJUST && total_supply != total_demand) {
1053 return UNBALANCED;
1054 }
1055
1056 // Feasibility checking, and possible supply/demand adjustment, is done by:
1057 // 1. Creating a new source and sink node.
1058 // 2. Taking all nodes that have a non-zero supply or demand and
1059 // connecting them to the source or sink respectively. The arc thus
1060 // added has a capacity of the supply or demand.
1061 // 3. Computing the max flow between the new source and sink.
1062 // 4. If adjustment isn't being done, checking that the max flow is equal
1063 // to the total supply/demand (and returning INFEASIBLE if it isn't).
1064 // 5. Running min-cost max-flow on this augmented graph, using the max
1065 // flow computed in step 3 as the supply of the source and demand of
1066 // the sink.
1067 const ArcIndex augmented_num_arcs =
1068 num_arcs + supply_node_count + demand_node_count;
1069 const NodeIndex source = num_nodes;
1070 const NodeIndex sink = num_nodes + 1;
1071 const NodeIndex augmented_num_nodes = num_nodes + 2;
1072
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]);
1076 }
1077
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);
1083 }
1084 }
1085
1086 graph.Build(&arc_permutation_);
1087
1088 {
1089 GenericMaxFlow<Graph> max_flow(&graph, source, sink);
1090 ArcIndex arc;
1091 for (arc = 0; arc < num_arcs; ++arc) {
1092 max_flow.SetArcCapacity(PermutedArc(arc), arc_capacity_[arc]);
1093 }
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]));
1097 ++arc;
1098 }
1099 }
1100 CHECK_EQ(arc, augmented_num_arcs);
1101 if (!max_flow.Solve()) {
1102 LOG(ERROR) << "Max flow could not be computed.";
1103 switch (max_flow.status()) {
1105 return NOT_SOLVED;
1107 LOG(ERROR)
1108 << "Max flow failed but claimed to have an optimal solution";
1109 ABSL_FALLTHROUGH_INTENDED;
1110 default:
1111 return BAD_RESULT;
1112 }
1113 }
1114 maximum_flow_ = max_flow.GetOptimalFlow();
1115 }
1116
1117 if (adjustment == DONT_ADJUST && maximum_flow_ != total_supply) {
1118 return INFEASIBLE;
1119 }
1120
1121 GenericMinCostFlow<Graph> min_cost_flow(&graph);
1122 ArcIndex arc;
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]);
1127 }
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);
1133 ++arc;
1134 }
1135 }
1136 min_cost_flow.SetNodeSupply(source, maximum_flow_);
1137 min_cost_flow.SetNodeSupply(sink, -maximum_flow_);
1138 min_cost_flow.SetCheckFeasibility(false);
1139
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));
1145 }
1146 }
1147 return min_cost_flow.status();
1148}
1149
1150CostValue SimpleMinCostFlow::OptimalCost() const { return optimal_cost_; }
1151
1152FlowQuantity SimpleMinCostFlow::MaximumFlow() const { return maximum_flow_; }
1153
1155 return arc_flow_[arc];
1156}
1157
1158NodeIndex SimpleMinCostFlow::NumNodes() const { return node_supply_.size(); }
1159
1160ArcIndex SimpleMinCostFlow::NumArcs() const { return arc_tail_.size(); }
1161
1162ArcIndex SimpleMinCostFlow::Tail(ArcIndex arc) const { return arc_tail_[arc]; }
1163
1164ArcIndex SimpleMinCostFlow::Head(ArcIndex arc) const { return arc_head_[arc]; }
1165
1167 return arc_capacity_[arc];
1168}
1169
1171 return arc_cost_[arc];
1172}
1173
1175 return node_supply_[node];
1176}
1177
1178void SimpleMinCostFlow::ResizeNodeVectors(NodeIndex node) {
1179 if (node < node_supply_.size()) return;
1180 node_supply_.resize(node + 1);
1181}
1182
1183} // namespace operations_research
int64 min
Definition: alldiff_cst.cc:138
int64 max
Definition: alldiff_cst.cc:139
#define DCHECK_LE(val1, val2)
Definition: base/logging.h:887
#define CHECK_EQ(val1, val2)
Definition: base/logging.h:697
#define DCHECK_GE(val1, val2)
Definition: base/logging.h:889
#define DCHECK_GT(val1, val2)
Definition: base/logging.h:890
#define LOG(severity)
Definition: base/logging.h:420
#define DCHECK(condition)
Definition: base/logging.h:884
#define DCHECK_EQ(val1, val2)
Definition: base/logging.h:885
#define VLOG(verboselevel)
Definition: base/logging.h:978
ArcIndexType AddArc(NodeIndexType tail, NodeIndexType head)
Definition: ebert_graph.h:1001
NodeIndexType Tail(const ArcIndexType arc) const
Definition: ebert_graph.h:1376
FlowQuantity GetOptimalFlow() const
Definition: max_flow.h:361
void SetArcCapacity(ArcIndex arc, FlowQuantity new_capacity)
Definition: max_flow.cc:175
FlowQuantity Flow(ArcIndex arc) const
Definition: max_flow.h:365
CostValue UnitCost(ArcIndex arc) const
FlowQuantity FeasibleSupply(NodeIndex node) const
FlowQuantity Flow(ArcIndex arc) const
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)
static T Abs(const T x)
Definition: mathutil.h:95
CostValue UnitCost(ArcIndex arc) const
ArcIndex AddArcWithCapacityAndUnitCost(NodeIndex tail, NodeIndex head, FlowQuantity capacity, CostValue unit_cost)
FlowQuantity Flow(ArcIndex arc) 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
NodeIndexType Head(const ArcIndexType arc) const
Definition: ebert_graph.h:297
bool Reserve(int64 new_min_index, int64 new_max_index)
Definition: zvector.h:98
GurobiMPCallbackContext * context
short int16
int int32
int64_t int64
uint64_t uint64
const int ERROR
Definition: log_severity.h:32
ABSL_FLAG(int64, min_cost_flow_alpha, 5, "Divide factor for epsilon at each refine step.")
Definition: cleanup.h:22
The vehicle routing library lets one model and solve generic vehicle routing problems ranging from th...
ListGraph Graph
Definition: graph.h:2360
int index
Definition: pack.cc:508
int64 delta
Definition: resource.cc:1684
int64 tail
int64 cost
int64 head
int64 capacity
#define IF_STATS_ENABLED(instructions)
Definition: stats.h:437
#define SCOPED_TIME_STAT(stats)
Definition: stats.h:438
static ArcIndex ArcReservation(const Graph &graph)
Definition: graphs.h:39
static NodeIndex NodeReservation(const Graph &graph)
Definition: graphs.h:36
static bool IsArcValid(const Graph &graph, ArcIndex arc)
Definition: graphs.h:33
static ArcIndex OppositeArc(const Graph &graph, ArcIndex arc)
Definition: graphs.h:30