Google OR-Tools v9.12
a fast and portable software suite for combinatorial optimization
Loading...
Searching...
No Matches
min_cost_flow.h
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
14// An implementation of a cost-scaling push-relabel algorithm for
15// the min-cost flow problem.
16//
17// In the following, we consider a graph G = (V,E) where V denotes the set
18// of nodes (vertices) in the graph, E denotes the set of arcs (edges).
19// n = |V| denotes the number of nodes in the graph, and m = |E| denotes the
20// number of arcs in the graph.
21//
22// With each arc (v,w) is associated a nonnegative capacity u(v,w)
23// (where 'u' stands for "upper bound") and a unit cost c(v,w). With
24// each node v is associated a quantity named supply(v), which
25// represents a supply of fluid (if >0) or a demand (if <0).
26// Furthermore, no fluid is created in the graph so
27// sum_{v in V} supply(v) = 0.
28//
29// A flow is a function from E to R such that:
30// a) f(v,w) <= u(v,w) for all (v,w) in E (capacity constraint).
31// b) f(v,w) = -f(w,v) for all (v,w) in E (flow antisymmetry constraint).
32// c) sum on v f(v,w) + supply(w) = 0 (flow conservation).
33//
34// The cost of a flow is sum on (v,w) in E ( f(v,w) * c(v,w) ) [Note:
35// It can be confusing to beginners that the cost is actually double
36// the amount that it might seem at first because of flow
37// antisymmetry.]
38//
39// The problem to solve: find a flow of minimum cost such that all the
40// fluid flows from the supply nodes to the demand nodes.
41//
42// The principles behind this algorithm are the following:
43// 1/ handle pseudo-flows instead of flows and refine pseudo-flows until an
44// epsilon-optimal minimum-cost flow is obtained,
45// 2/ deal with epsilon-optimal pseudo-flows.
46//
47// 1/ A pseudo-flow is like a flow, except that a node's outflow minus
48// its inflow can be different from its supply. If it is the case at a
49// given node v, it is said that there is an excess (or deficit) at
50// node v. A deficit is denoted by a negative excess and inflow =
51// outflow + excess.
52// (Look at ortools/graph/max_flow.h to see that the definition
53// of preflow is more restrictive than the one for pseudo-flow in that a preflow
54// only allows non-negative excesses, i.e., no deficit.)
55// More formally, a pseudo-flow is a function f such that:
56// a) f(v,w) <= u(v,w) for all (v,w) in E (capacity constraint).
57// b) f(v,w) = -f(w,v) for all (v,w) in E (flow antisymmetry constraint).
58//
59// For each v in E, we also define the excess at node v, the algebraic sum of
60// all the incoming preflows at this node, added together with the supply at v.
61// excess(v) = sum on u f(u,v) + supply(v)
62//
63// The goal of the algorithm is to obtain excess(v) = 0 for all v in V, while
64// consuming capacity on some arcs, at the lowest possible cost.
65//
66// 2/ Internally to the algorithm and its analysis (but invisibly to
67// the client), each node has an associated "price" (or potential), in
68// addition to its excess. It is formally a function from E to R (the
69// set of real numbers.). For a given price function p, the reduced
70// cost of an arc (v,w) is:
71// c_p(v,w) = c(v,w) + p(v) - p(w)
72// (c(v,w) is the cost of arc (v,w).) For those familiar with linear
73// programming, the price function can be viewed as a set of dual
74// variables.
75//
76// For a constant epsilon >= 0, a pseudo-flow f is said to be epsilon-optimal
77// with respect to a price function p if for every residual arc (v,w) in E,
78// c_p(v,w) >= -epsilon.
79//
80// A flow f is optimal if and only if there exists a price function p such that
81// no arc is admissible with respect to f and p.
82//
83// If the arc costs are integers, and epsilon < 1/n, any epsilon-optimal flow
84// is optimal. The integer cost case is handled by multiplying all the arc costs
85// and the initial value of epsilon by (n+1). When epsilon reaches 1, and
86// the solution is epsilon-optimal, it means: for all residual arc (v,w) in E,
87// (n+1) * c_p(v,w) >= -1, thus c_p(v,w) >= -1/(n+1) > -1/n, and the
88// solution is optimal.
89//
90// A node v is said to be *active* if excess(v) > 0.
91// In this case the following operations can be applied to it:
92// - if there are *admissible* incident arcs, i.e. arcs which are not saturated,
93// and whose reduced costs are negative, a PushFlow operation can
94// be applied. It consists in sending as much flow as both the excess at the
95// node and the capacity of the arc permit.
96// - if there are no admissible arcs, the active node considered is relabeled,
97// This is implemented in Discharge, which itself calls PushFlow and Relabel.
98//
99// Discharge itself is called by Refine. Refine first saturates all the
100// admissible arcs, then builds a stack of active nodes. It then applies
101// Discharge for each active node, possibly adding new ones in the process,
102// until no nodes are active. In that case an epsilon-optimal flow is obtained.
103//
104// Optimize iteratively calls Refine, while epsilon > 1, and divides epsilon by
105// alpha (set by default to 5) before each iteration.
106//
107// The algorithm starts with epsilon = C, where C is the maximum absolute value
108// of the arc costs. In the integer case which we are dealing with, since all
109// costs are multiplied by (n+1), the initial value of epsilon is (n+1)*C.
110// The algorithm terminates when epsilon = 1, and Refine() has been called.
111// In this case, a minimum-cost flow is obtained.
112//
113// The complexity of the algorithm is O(n^2*m*log(n*C)) where C is the value of
114// the largest arc cost in the graph.
115//
116// IMPORTANT:
117// The algorithm is not able to detect the infeasibility of a problem (i.e.,
118// when a bottleneck in the network prohibits sending all the supplies.)
119// Worse, it could in some cases loop forever. This is why feasibility checking
120// is enabled by default (FLAGS_min_cost_flow_check_feasibility=true.)
121// Feasibility checking is implemented using a max-flow, which has a much lower
122// complexity. The impact on performance is negligible, while the risk of being
123// caught in an endless loop is removed. Note that using the feasibility checker
124// roughly doubles the memory consumption.
125//
126// The starting reference for this class of algorithms is:
127// A.V. Goldberg and R.E. Tarjan, "Finding Minimum-Cost Circulations by
128// Successive Approximation." Mathematics of Operations Research, Vol. 15,
129// 1990:430-466.
130// http://portal.acm.org/citation.cfm?id=92225
131//
132// Implementation issues are tackled in:
133// A.V. Goldberg, "An Efficient Implementation of a Scaling Minimum-Cost Flow
134// Algorithm," Journal of Algorithms, (1997) 22:1-29
135// http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.31.258
136//
137// A.V. Goldberg and M. Kharitonov, "On Implementing Scaling Push-Relabel
138// Algorithms for the Minimum-Cost Flow Problem", Network flows and matching:
139// First DIMACS implementation challenge, DIMACS Series in Discrete Mathematics
140// and Theoretical Computer Science, (1993) 12:157-198.
141// ftp://dimacs.rutgers.edu/pub/netflow/submit/papers/Goldberg-mincost/scalmin.ps
142// and in:
143// U. Bunnagel, B. Korte, and J. Vygen. “Efficient implementation of the
144// Goldberg-Tarjan minimum-cost flow algorithm.” Optimization Methods and
145// Software (1998) vol. 10, no. 2:157-174.
146// http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.84.9897
147//
148// We have tried as much as possible in this implementation to keep the
149// notations and namings of the papers cited above, except for 'demand' or
150// 'balance' which have been replaced by 'supply', with the according sign
151// changes to better accommodate with the API of the rest of our tools. A demand
152// is denoted by a negative supply.
153//
154// TODO(user): See whether the following can bring any improvements on real-life
155// problems.
156// R.K. Ahuja, A.V. Goldberg, J.B. Orlin, and R.E. Tarjan, "Finding minimum-cost
157// flows by double scaling," Mathematical Programming, (1992) 53:243-266.
158// http://www.springerlink.com/index/gu7404218u6kt166.pdf
159//
160// An interesting general reference on network flows is:
161// R. K. Ahuja, T. L. Magnanti, J. B. Orlin, "Network Flows: Theory, Algorithms,
162// and Applications," Prentice Hall, 1993, ISBN: 978-0136175490,
163// http://www.amazon.com/dp/013617549X
164//
165// Keywords: Push-relabel, min-cost flow, network, graph, Goldberg, Tarjan,
166// Dinic, Dinitz.
167
168#ifndef OR_TOOLS_GRAPH_MIN_COST_FLOW_H_
169#define OR_TOOLS_GRAPH_MIN_COST_FLOW_H_
170
171#include <cstdint>
172#include <stack>
173#include <string>
174#include <vector>
175
176#include "absl/strings/string_view.h"
177#include "ortools/base/logging.h"
178#include "ortools/graph/graph.h"
179#include "ortools/util/stats.h"
180#include "ortools/util/zvector.h"
181
182namespace operations_research {
183
184// Forward declaration.
185template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
187
188// Different statuses for a solved problem.
189// We use a base class to share it between our different interfaces.
191 public:
192 enum Status {
199
200 // This is returned when an integer overflow occurred during the algorithm
201 // execution.
202 //
203 // Some details on how to deal with this:
204 // - Since we scale cost, each arc unit cost times (num_nodes + 1) should
205 // not overflow. We detect that at the beginning of the Solve().
206 // - This is however not sufficient as the node potential depends on the
207 // minimum cost of sending 1 unit of flow through the residual graph. If
208 // the maximum arc unit cost is smaller than kint64max / (2 * n ^ 2) then
209 // one shouldn't run into any overflow. But in pratice this bound is quite
210 // loose. So it is possible to try with higher cost, and the algorithm
211 // will detect if an overflow actually happen and return BAD_COST_RANGE,
212 // so we can retry with smaller cost.
213 //
214 // And two remarks:
215 // - Note that the complexity of the algorithm depends on the maximum cost,
216 // so it is usually a good idea to use unit cost that are as small as
217 // possible.
218 // - Even if there is no overflow, note that the total cost can easily not
219 // fit on an int64_t since it is the product of the unit cost times the
220 // actual amount of flow sent. This is easy to detect since the final
221 // optimal cost will be set to kint64max. It is also relatively easy to
222 // deal with since we will still have the proper flow on each arc. It is
223 // thus possible to recompute the total cost in double or using
224 // absl::int128 if the need arise.
226
227 // This is returned in our initial check if the arc capacity are too large.
228 // For each node these quantity should not overflow the FlowQuantity type
229 // which is int64_t by default.
230 // - Its initial excess (+supply or -demand) + sum incoming arc capacity
231 // - Its initial excess (+supply or -demand) - sum outgoing arc capacity.
232 //
233 // Note that these are upper bounds that guarantee that no overflow will
234 // take place during the algorithm execution. It is possible to go over that
235 // and still encounter no overflow, but since we cannot guarantee this, we
236 // rather check early and return a proper status.
237 //
238 // Note that we might cap the capacity of the arcs if we can detect it
239 // is too large in order to avoid returning this status for simple case,
240 // like if a client used int64_t max for all arcs out of the source.
241 //
242 // TODO(user): Not sure this is a good idea, probably better to make sure
243 // client use reasonable capacities. Also we should template by FlowQuantity
244 // and allow use of absl::int128 so we never have issue if the input is
245 // int64_t.
247 };
248};
249
250// A simple and efficient min-cost flow interface. This is as fast as
251// GenericMinCostFlow<ReverseArcStaticGraph>, which is the fastest, but is uses
252// more memory in order to hide the somewhat involved construction of the
253// static graph.
254//
255// TODO(user): If the need arises, extend this interface to support warm start
256// and incrementality between solves. Note that this is already supported by the
257// GenericMinCostFlow<> interface.
259 public:
260 typedef int32_t NodeIndex;
261 typedef int32_t ArcIndex;
262 typedef int64_t FlowQuantity;
263 typedef int64_t CostValue;
264
265 // By default, the constructor takes no size. New node indices are created
266 // lazily by AddArcWithCapacityAndUnitCost() or SetNodeSupply() such that the
267 // set of valid nodes will always be [0, NumNodes()).
268 //
269 // You may pre-reserve the internal data structures with a given expected
270 // number of nodes and arcs, to potentially gain performance.
271 explicit SimpleMinCostFlow(NodeIndex reserve_num_nodes = 0,
272 ArcIndex reserve_num_arcs = 0);
273
274#ifndef SWIG
275 // This type is neither copyable nor movable.
278#endif
279
280 // Adds a directed arc from tail to head to the underlying graph with
281 // a given capacity and cost per unit of flow.
282 // * Node indices and the capacity must be non-negative (>= 0).
283 // * The unit cost can take any integer value (even negative).
284 // * Self-looping and duplicate arcs are supported.
285 // * After the method finishes, NumArcs() == the returned ArcIndex + 1.
287 FlowQuantity capacity,
288 CostValue unit_cost);
289
290 // Modifies the capacity of the given arc. The arc index must be non-negative
291 // (>= 0); it must be an index returned by a previous call to
292 // AddArcWithCapacityAndUnitCost().
293 void SetArcCapacity(ArcIndex arc, FlowQuantity capacity);
294
295 // Sets the supply of the given node. The node index must be non-negative (>=
296 // 0). Nodes implicitly created will have a default supply set to 0. A demand
297 // is modeled as a negative supply.
298 void SetNodeSupply(NodeIndex node, FlowQuantity supply);
299
300 // Solves the problem, and returns the problem status. This function
301 // requires that the sum of all node supply minus node demand is zero and
302 // that the graph has enough capacity to send all supplies and serve all
303 // demands. Otherwise, it will return INFEASIBLE.
305 return SolveWithPossibleAdjustment(SupplyAdjustment::DONT_ADJUST);
306 }
307
308 // Same as Solve(), but does not have the restriction that the supply
309 // must match the demand or that the graph has enough capacity to serve
310 // all the demand or use all the supply. This will compute a maximum-flow
311 // with minimum cost. The value of the maximum-flow will be given by
312 // MaximumFlow().
314 return SolveWithPossibleAdjustment(SupplyAdjustment::ADJUST);
315 }
316
317 // Returns the cost of the minimum-cost flow found by the algorithm when
318 // the returned Status is OPTIMAL.
319 CostValue OptimalCost() const;
320
321 // Returns the total flow of the minimum-cost flow found by the algorithm
322 // when the returned Status is OPTIMAL.
324
325 // Returns the flow on arc, this only make sense for a successful Solve().
326 //
327 // Note: It is possible that there is more than one optimal solution. The
328 // algorithm is deterministic so it will always return the same solution for
329 // a given problem. However, there is no guarantee of this from one code
330 // version to the next (but the code does not change often).
331 FlowQuantity Flow(ArcIndex arc) const;
332
333 // Accessors for the user given data. The implementation will crash if "arc"
334 // is not in [0, NumArcs()) or "node" is not in [0, NumNodes()).
335 NodeIndex NumNodes() const;
336 ArcIndex NumArcs() const;
337 NodeIndex Tail(ArcIndex arc) const;
338 NodeIndex Head(ArcIndex arc) const;
339 FlowQuantity Capacity(ArcIndex arc) const;
340 FlowQuantity Supply(NodeIndex node) const;
341 CostValue UnitCost(ArcIndex arc) const;
342
343 // Advanced usage. The default is true.
344 //
345 // Without cost scaling, the algorithm will return a 1-optimal solution to the
346 // given problem. The solution will be an optimal solution to a perturbed
347 // problem where some of the arc unit costs are changed by at most 1.
348 //
349 // If the cost are initially integer and we scale them by (num_nodes + 1),
350 // then we can show that such 1-optimal solution is actually optimal. This
351 // is what happen by default or when SetPriceScaling(true) is called.
352 //
353 // However, if your cost were originally double, you don't really care to
354 // solve optimally a problem where the weights are approximated in the first
355 // place. It is better to multiply your double by a scaling_factor (prefer a
356 // power of 2) so that the maximum rounded arc unit cost is under kint64max /
357 // (num_nodes + 1) to prevent any overflow. You can then solve without any
358 // cost scaling. The final result will be the optimal to a problem were the
359 // unit cost of some arc as been changed by at most 1 / scaling_factor.
360 void SetPriceScaling(bool value) { scale_prices_ = value; }
361
362 private:
363 typedef ::util::ReverseArcStaticGraph<NodeIndex, ArcIndex> Graph;
364 enum SupplyAdjustment { ADJUST, DONT_ADJUST };
365
366 // Applies the permutation in arc_permutation_ to the given arc index.
367 ArcIndex PermutedArc(ArcIndex arc);
368 // Solves the problem, potentially applying supply and demand adjustment,
369 // and returns the problem status.
370 Status SolveWithPossibleAdjustment(SupplyAdjustment adjustment);
371 void ResizeNodeVectors(NodeIndex node);
372
373 std::vector<NodeIndex> arc_tail_;
374 std::vector<NodeIndex> arc_head_;
375 std::vector<FlowQuantity> arc_capacity_;
376 std::vector<FlowQuantity> node_supply_;
377 std::vector<CostValue> arc_cost_;
378 std::vector<ArcIndex> arc_permutation_;
379 std::vector<FlowQuantity> arc_flow_;
380 CostValue optimal_cost_;
381 FlowQuantity maximum_flow_;
382
383 bool scale_prices_ = true;
384};
385
386// Generic MinCostFlow that works with all the graphs handling reverse arcs from
387// graph.h, see the end of min_cost_flow.cc for the exact types this class is
388// compiled for.
389//
390// One can greatly decrease memory usage by using appropriately small integer
391// types:
392// - For the Graph<> types, i.e. NodeIndexType and ArcIndexType, see graph.h.
393// - ArcFlowType is used for the *per-arc* flow quantity. It must be signed, and
394// large enough to hold the maximum arc capacity and its negation.
395// - ArcScaledCostType is used for a per-arc scaled cost. It must be signed
396// and large enough to hold the maximum unit cost of an arc times
397// (num_nodes + 1).
398//
399// Note that the latter two are different than FlowQuantity and CostValue, which
400// are used for global, aggregated values and may need to be larger.
401//
402// Also uses the Arc*Type where there is no risk of overflow in more places.
403template <typename Graph, typename ArcFlowType = int64_t,
404 typename ArcScaledCostType = int64_t>
406 public:
407 typedef typename Graph::NodeIndex NodeIndex;
408 typedef typename Graph::ArcIndex ArcIndex;
409 typedef int64_t CostValue;
410 typedef int64_t FlowQuantity;
416
417 // Initialize a MinCostFlow instance on the given graph. The graph does not
418 // need to be fully built yet, but its capacity reservation is used to
419 // initialize the memory of this class.
420 explicit GenericMinCostFlow(const Graph* graph);
421
422#ifndef SWIG
423 // This type is neither copyable nor movable.
426#endif
427
428 // Returns the graph associated to the current object.
429 const Graph* graph() const { return graph_; }
430
431 // Returns the status of last call to Solve(). NOT_SOLVED is returned if
432 // Solve() has never been called or if the problem has been modified in such a
433 // way that the previous solution becomes invalid.
434 Status status() const { return status_; }
435
436 // Sets the supply corresponding to node. A demand is modeled as a negative
437 // supply.
438 void SetNodeSupply(NodeIndex node, FlowQuantity supply);
439
440 // Sets the unit cost for the given arc.
441 void SetArcUnitCost(ArcIndex arc, ArcScaledCostType unit_cost);
442
443 // Sets the capacity for the given arc.
444 void SetArcCapacity(ArcIndex arc, ArcFlowType new_capacity);
445
446 // Solves the problem, returning true if a min-cost flow could be found.
447 bool Solve();
448
449 // Returns the cost of the minimum-cost flow found by the algorithm. This
450 // works in O(num_arcs). This will only work if the last Solve() call was
451 // successful and returned true, otherwise it will return 0. Note that the
452 // computation might overflow, in which case we will cap the cost at
453 // std::numeric_limits<CostValue>::max().
455
456 // Returns the flow on the given arc using the equations given in the
457 // comment on residual_arc_capacity_.
458 FlowQuantity Flow(ArcIndex arc) const;
459
460 // Returns the capacity of the given arc.
461 //
462 // Warning: If the capacity were close to ArcFlowType::max() we might have
463 // adjusted them in order to avoid overflow.
464 FlowQuantity Capacity(ArcIndex arc) const;
465
466 // Returns the unscaled cost for the given arc.
467 CostValue UnitCost(ArcIndex arc) const;
468
469 // Returns the supply at a given node.
470 // Demands are modelled as negative supplies.
471 FlowQuantity Supply(NodeIndex node) const;
472
473 // Whether to check the feasibility of the problem with a max-flow, prior to
474 // solving it. This uses about twice as much memory, but detects infeasible
475 // problems (where the flow can't be satisfied) and makes Solve() return
476 // INFEASIBLE. If you disable this check, you will spare memory but you must
477 // make sure that your problem is feasible, otherwise the code can loop
478 // forever.
479 void SetCheckFeasibility(bool value) { check_feasibility_ = value; }
480
481 // Algorithm options.
482 void SetUseUpdatePrices(bool value) { use_price_update_ = value; }
483 void SetPriceScaling(bool value) { scale_prices_ = value; }
484
485 private:
486 // Checks for feasibility, i.e., that all the supplies and demands can be
487 // matched without exceeding bottlenecks in the network.
488 // Note that CheckFeasibility is called by Solve() when SetCheckFeasibility()
489 // is set to true, which is the default.
490 bool CheckFeasibility();
491
492 // Returns true if the given arc is admissible i.e. if its residual capacity
493 // is strictly positive, and its reduced cost strictly negative, i.e., pushing
494 // more flow into it will result in a reduction of the total cost.
495 bool IsAdmissible(ArcIndex arc) const;
496 bool FastIsAdmissible(ArcIndex arc, CostValue tail_potential) const;
497
498 // Returns true if node is active, i.e., if its supply is positive.
499 bool IsActive(NodeIndex node) const;
500
501 // Returns the reduced cost for a given arc.
502 CostValue ReducedCost(ArcIndex arc) const;
503 CostValue FastReducedCost(ArcIndex arc, CostValue tail_potential) const;
504
505 // Returns the first incident arc of a given node.
506 ArcIndex GetFirstOutgoingOrOppositeIncomingArc(NodeIndex node) const;
507
508 // Checks the consistency of the input, i.e., whether the sum of the supplies
509 // for all nodes is equal to zero. To be used in a DCHECK.
510 bool CheckInputConsistency();
511
512 // Checks whether the result is valid, i.e. whether for each arc,
513 // residual_arc_capacity_[arc] == 0 || ReducedCost(arc) >= -epsilon_.
514 // (A solution is epsilon-optimal if ReducedCost(arc) >= -epsilon.)
515 // To be used in a DCHECK.
516 bool CheckResult() const;
517
518 // Checks the relabel precondition (to be used in a DCHECK):
519 // - The node must be active, or have a 0 excess (relaxation for the Push
520 // Look-Ahead heuristic).
521 // - The node must have no admissible arcs.
522 bool CheckRelabelPrecondition(NodeIndex node) const;
523
524 // Returns context concatenated with information about a given arc
525 // in a human-friendly way.
526 std::string DebugString(absl::string_view context, ArcIndex arc) const;
527
528 // Resets the first_admissible_arc_ array to the first incident arc of each
529 // node.
530 void ResetFirstAdmissibleArcs();
531
532 // Scales the costs, by multiplying them by (graph_->num_nodes() + 1).
533 // Returns false on integer overflow.
534 bool ScaleCosts();
535
536 // Unscales the costs, by dividing them by (graph_->num_nodes() + 1).
537 void UnscaleCosts();
538
539 // Optimizes the cost by dividing epsilon_ by alpha_ and calling Refine().
540 // Returns false on integer overflow.
541 bool Optimize();
542
543 // Saturates the admissible arcs, i.e., push as much flow as possible.
544 void SaturateAdmissibleArcs();
545
546 // Pushes flow on a given arc, i.e., consumes flow on
547 // residual_arc_capacity_[arc], and consumes -flow on
548 // residual_arc_capacity_[Opposite(arc)]. Updates node_excess_ at the tail
549 // and head of the arc accordingly.
550 void PushFlow(FlowQuantity flow, ArcIndex arc);
551 void FastPushFlow(FlowQuantity flow, ArcIndex arc, NodeIndex tail);
552
553 // Initializes the stack active_nodes_.
554 void InitializeActiveNodeStack();
555
556 // Price update heuristics as described in A.V. Goldberg, "An Efficient
557 // Implementation of a Scaling Minimum-Cost Flow Algorithm," Journal of
558 // Algorithms, (1997) 22:1-29
559 // http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.31.258
560 // Returns false on overflow or infeasibility.
561 bool UpdatePrices();
562
563 // Performs an epsilon-optimization step by saturating admissible arcs
564 // and discharging the active nodes.
565 // Returns false on overflow or infeasibility.
566 bool Refine();
567
568 // Discharges an active node by saturating its admissible adjacent arcs,
569 // if any, and by relabelling it when it becomes inactive.
570 // Returns false on overflow or infeasibility.
571 bool Discharge(NodeIndex node);
572
573 // Part of the Push LookAhead heuristic.
574 // Returns true iff the node has admissible arc.
575 bool NodeHasAdmissibleArc(NodeIndex node);
576
577 // Relabels node, i.e., decreases its potential while keeping the
578 // epsilon-optimality of the pseudo flow. See CheckRelabelPrecondition() for
579 // details on the preconditions.
580 // Returns false on overflow or infeasibility.
581 bool Relabel(NodeIndex node);
582
583 // Handy member functions to make the code more compact.
584 NodeIndex Head(ArcIndex arc) const { return graph_->Head(arc); }
585 NodeIndex Tail(ArcIndex arc) const { return graph_->Tail(arc); }
586 ArcIndex Opposite(ArcIndex arc) const;
587 bool IsArcDirect(ArcIndex arc) const;
588 bool IsArcValid(ArcIndex arc) const;
589
590 // Pointer to the graph passed as argument.
591 const Graph* graph_;
592
593 // An array representing the supply (if > 0) or the demand (if < 0)
594 // for each node in graph_.
595 std::vector<FlowQuantity> node_excess_;
596
597 // An array representing the potential (or price function) for
598 // each node in graph_.
599 std::vector<CostValue> node_potential_;
600
601 // An array representing the residual_capacity for each arc in graph_.
602 // Residual capacities enable one to represent the capacity and flow for all
603 // arcs in the graph in the following manner.
604 // For all arcs, residual_arc_capacity_[arc] = capacity[arc] - flow[arc]
605 // Moreover, for reverse arcs, capacity[arc] = 0 by definition.
606 // Also flow[Opposite(arc)] = -flow[arc] by definition.
607 // Therefore:
608 // - for a direct arc:
609 // flow[arc] = 0 - flow[Opposite(arc)]
610 // = capacity[Opposite(arc)] - flow[Opposite(arc)]
611 // = residual_arc_capacity_[Opposite(arc)]
612 // - for a reverse arc:
613 // flow[arc] = -residual_arc_capacity_[arc]
614 // Using these facts enables one to only maintain residual_arc_capacity_,
615 // instead of both capacity and flow, for each direct and indirect arc. This
616 // reduces the amount of memory for this information by a factor 2.
617 // Note that the sum of the largest capacity of an arc in the graph and of
618 // the total flow in the graph mustn't exceed the largest 64 bit integer
619 // to avoid errors. CheckInputConsistency() verifies this constraint.
620 ZVector<ArcFlowType> residual_arc_capacity_;
621
622 // An array representing the first admissible arc for each node in graph_.
623 std::vector<ArcIndex> first_admissible_arc_;
624
625 // A stack used for managing active nodes in the algorithm.
626 // Note that the papers cited above recommend the use of a queue, but
627 // benchmarking so far has not proved it is better.
628 std::stack<NodeIndex> active_nodes_;
629
630 // epsilon_ is the tolerance for optimality.
631 CostValue epsilon_ = 0;
632
633 // alpha_ is the factor by which epsilon_ is divided at each iteration of
634 // Refine().
635 const int64_t alpha_;
636
637 // cost_scaling_factor_ is the scaling factor for cost.
638 CostValue cost_scaling_factor_ = 1;
639
640 // Our node potentials should be >= this at all time.
641 // The first time a potential cross this, BAD_COST_RANGE status is returned.
642 CostValue overflow_threshold_;
643
644 // An array representing the scaled unit cost for each arc in graph_.
645 ZVector<ArcScaledCostType> scaled_arc_unit_cost_;
646
647 // The status of the problem.
648 Status status_ = NOT_SOLVED;
649
650 // An array containing the initial excesses (i.e. the supplies) for each
651 // node. This is used to create the max-flow-based feasibility checker.
652 std::vector<FlowQuantity> initial_node_excess_;
653
654 // Statistics about this class.
655 StatsGroup stats_;
656
657 // Number of Relabel() since last UpdatePrices().
658 int num_relabels_since_last_price_update_ = 0;
659
660 // A Boolean which is true when feasibility has been checked.
661 bool feasibility_checked_ = false;
662
663 // Whether to use the UpdatePrices() heuristic.
664 bool use_price_update_ = false;
665
666 // Whether to check the problem feasibility with a max-flow.
667 bool check_feasibility_ = false;
668
669 // Whether to scale prices, see SimpleMinCostFlow::SetPriceScaling().
670 bool scale_prices_ = true;
671};
672
673#if !SWIG
674
675// Note: SWIG does not seem to understand explicit template specialization and
676// instantiation declarations.
677
681extern template class GenericMinCostFlow<
682 ::util::ReverseArcStaticGraph<uint16_t, int32_t>>;
683extern template class GenericMinCostFlow<
684 ::util::ReverseArcListGraph<int64_t, int64_t>, int64_t, int64_t>;
685extern template class GenericMinCostFlow<
686 ::util::ReverseArcStaticGraph<uint16_t, int32_t>,
687 /*ArcFlowType=*/int16_t,
688 /*ArcScaledCostType=*/int32_t>;
689
690// TODO(b/385094969): Remove this alias after 2025-07-01 to give or-tools users
691// a grace period.
693 template <typename = void>
695 LOG(FATAL) << "MinCostFlow is deprecated. Use `SimpleMinCostFlow` or "
696 "`GenericMinCostFlow` with a specific graph type instead.";
697 }
698};
699
700#endif // SWIG
701
702} // namespace operations_research
703#endif // OR_TOOLS_GRAPH_MIN_COST_FLOW_H_
FlowQuantity Flow(ArcIndex arc) const
void SetArcUnitCost(ArcIndex arc, ArcScaledCostType unit_cost)
Sets the unit cost for the given arc.
Graph::IncomingArcIterator IncomingArcIterator
Graph::OutgoingOrOppositeIncomingArcIterator OutgoingOrOppositeIncomingArcIterator
GenericMinCostFlow & operator=(const GenericMinCostFlow &)=delete
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.
Graph::OutgoingArcIterator OutgoingArcIterator
GenericMinCostFlow(const GenericMinCostFlow &)=delete
This type is neither copyable nor movable.
void SetNodeSupply(NodeIndex node, FlowQuantity supply)
void SetUseUpdatePrices(bool value)
Algorithm options.
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_.
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)
SimpleMinCostFlow(const SimpleMinCostFlow &)=delete
This type is neither copyable nor movable.
NodeIndex Head(ArcIndex arc) const
SimpleMinCostFlow(NodeIndex reserve_num_nodes=0, ArcIndex reserve_num_arcs=0)
SimpleMinCostFlow & operator=(const SimpleMinCostFlow &)=delete
FlowQuantity Supply(NodeIndex node) const
NodeIndexType Tail(ArcIndexType arc) const
Definition graph.h:1782
NodeIndexType Head(ArcIndexType arc) const
Definition graph.h:1774
In SWIG mode, we don't want anything besides these top-level includes.
util::ReverseArcStaticGraph Graph
Type of graph to use.