Google OR-Tools v9.14
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/log/log.h"
177#include "absl/strings/string_view.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;
414
415 // Initialize a MinCostFlow instance on the given graph. The graph does not
416 // need to be fully built yet, but its capacity reservation is used to
417 // initialize the memory of this class.
418 explicit GenericMinCostFlow(const Graph* graph);
419
420#ifndef SWIG
421 // This type is neither copyable nor movable.
424#endif
425
426 // Returns the graph associated to the current object.
427 const Graph* graph() const { return graph_; }
428
429 // Returns the status of last call to Solve(). NOT_SOLVED is returned if
430 // Solve() has never been called or if the problem has been modified in such a
431 // way that the previous solution becomes invalid.
432 Status status() const { return status_; }
433
434 // Sets the supply corresponding to node. A demand is modeled as a negative
435 // supply.
436 void SetNodeSupply(NodeIndex node, FlowQuantity supply);
437
438 // Sets the unit cost for the given arc.
439 void SetArcUnitCost(ArcIndex arc, ArcScaledCostType unit_cost);
440
441 // Sets the capacity for the given arc.
442 void SetArcCapacity(ArcIndex arc, ArcFlowType new_capacity);
443
444 // Solves the problem, returning true if a min-cost flow could be found.
445 bool Solve();
446
447 // Returns the cost of the minimum-cost flow found by the algorithm. This
448 // works in O(num_arcs). This will only work if the last Solve() call was
449 // successful and returned true, otherwise it will return 0. Note that the
450 // computation might overflow, in which case we will cap the cost at
451 // std::numeric_limits<CostValue>::max().
453
454 // Returns the flow on the given arc using the equations given in the
455 // comment on residual_arc_capacity_.
456 FlowQuantity Flow(ArcIndex arc) const;
457
458 // Returns the capacity of the given arc.
459 //
460 // Warning: If the capacity were close to ArcFlowType::max() we might have
461 // adjusted them in order to avoid overflow.
462 FlowQuantity Capacity(ArcIndex arc) const;
463
464 // Returns the unscaled cost for the given arc.
465 CostValue UnitCost(ArcIndex arc) const;
466
467 // Returns the supply at a given node.
468 // Demands are modelled as negative supplies.
469 FlowQuantity Supply(NodeIndex node) const;
470
471 // Whether to check the feasibility of the problem with a max-flow, prior to
472 // solving it. This uses about twice as much memory, but detects infeasible
473 // problems (where the flow can't be satisfied) and makes Solve() return
474 // INFEASIBLE. If you disable this check, you will spare memory but you must
475 // make sure that your problem is feasible, otherwise the code can loop
476 // forever.
477 void SetCheckFeasibility(bool value) { check_feasibility_ = value; }
478
479 // Algorithm options.
480 void SetUseUpdatePrices(bool value) { use_price_update_ = value; }
481 void SetPriceScaling(bool value) { scale_prices_ = value; }
482
483 private:
484 // Checks for feasibility, i.e., that all the supplies and demands can be
485 // matched without exceeding bottlenecks in the network.
486 // Note that CheckFeasibility is called by Solve() when SetCheckFeasibility()
487 // is set to true, which is the default.
488 bool CheckFeasibility();
489
490 // Returns true if the given arc is admissible i.e. if its residual capacity
491 // is strictly positive, and its reduced cost strictly negative, i.e., pushing
492 // more flow into it will result in a reduction of the total cost.
493 bool IsAdmissible(ArcIndex arc) const;
494 bool FastIsAdmissible(ArcIndex arc, CostValue tail_potential) const;
495
496 // Returns true if node is active, i.e., if its supply is positive.
497 bool IsActive(NodeIndex node) const;
498
499 // Returns the reduced cost for a given arc.
500 CostValue ReducedCost(ArcIndex arc) const;
501 CostValue FastReducedCost(ArcIndex arc, CostValue tail_potential) const;
502
503 // Returns the first incident arc of a given node.
504 ArcIndex GetFirstOutgoingOrOppositeIncomingArc(NodeIndex node) const;
505
506 // Checks the consistency of the input, i.e., whether the sum of the supplies
507 // for all nodes is equal to zero. To be used in a DCHECK.
508 bool CheckInputConsistency();
509
510 // Checks whether the result is valid, i.e. whether for each arc,
511 // residual_arc_capacity_[arc] == 0 || ReducedCost(arc) >= -epsilon_.
512 // (A solution is epsilon-optimal if ReducedCost(arc) >= -epsilon.)
513 // To be used in a DCHECK.
514 bool CheckResult() const;
515
516 // Checks the relabel precondition (to be used in a DCHECK):
517 // - The node must be active, or have a 0 excess (relaxation for the Push
518 // Look-Ahead heuristic).
519 // - The node must have no admissible arcs.
520 bool CheckRelabelPrecondition(NodeIndex node) const;
521
522 // Returns context concatenated with information about a given arc
523 // in a human-friendly way.
524 std::string DebugString(absl::string_view context, ArcIndex arc) const;
525
526 // Resets the first_admissible_arc_ array to the first incident arc of each
527 // node.
528 void ResetFirstAdmissibleArcs();
529
530 // Scales the costs, by multiplying them by (graph_->num_nodes() + 1).
531 // Returns false on integer overflow.
532 bool ScaleCosts();
533
534 // Unscales the costs, by dividing them by (graph_->num_nodes() + 1).
535 void UnscaleCosts();
536
537 // Optimizes the cost by dividing epsilon_ by alpha_ and calling Refine().
538 // Returns false on integer overflow.
539 bool Optimize();
540
541 // Saturates the admissible arcs, i.e., push as much flow as possible.
542 void SaturateAdmissibleArcs();
543
544 // Pushes flow on a given arc, i.e., consumes flow on
545 // residual_arc_capacity_[arc], and consumes -flow on
546 // residual_arc_capacity_[Opposite(arc)]. Updates node_excess_ at the tail
547 // and head of the arc accordingly.
548 void PushFlow(FlowQuantity flow, ArcIndex arc);
549 void FastPushFlow(FlowQuantity flow, ArcIndex arc, NodeIndex tail);
550
551 // Initializes the stack active_nodes_.
552 void InitializeActiveNodeStack();
553
554 // Price update heuristics as described in A.V. Goldberg, "An Efficient
555 // Implementation of a Scaling Minimum-Cost Flow Algorithm," Journal of
556 // Algorithms, (1997) 22:1-29
557 // http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.31.258
558 // Returns false on overflow or infeasibility.
559 bool UpdatePrices();
560
561 // Performs an epsilon-optimization step by saturating admissible arcs
562 // and discharging the active nodes.
563 // Returns false on overflow or infeasibility.
564 bool Refine();
565
566 // Discharges an active node by saturating its admissible adjacent arcs,
567 // if any, and by relabelling it when it becomes inactive.
568 // Returns false on overflow or infeasibility.
569 bool Discharge(NodeIndex node);
570
571 // Part of the Push LookAhead heuristic.
572 // Returns true iff the node has admissible arc.
573 bool NodeHasAdmissibleArc(NodeIndex node);
574
575 // Relabels node, i.e., decreases its potential while keeping the
576 // epsilon-optimality of the pseudo flow. See CheckRelabelPrecondition() for
577 // details on the preconditions.
578 // Returns false on overflow or infeasibility.
579 bool Relabel(NodeIndex node);
580
581 // Handy member functions to make the code more compact.
582 NodeIndex Head(ArcIndex arc) const { return graph_->Head(arc); }
583 NodeIndex Tail(ArcIndex arc) const { return graph_->Tail(arc); }
584 ArcIndex Opposite(ArcIndex arc) const;
585 bool IsArcDirect(ArcIndex arc) const;
586 bool IsArcValid(ArcIndex arc) const;
587
588 // Pointer to the graph passed as argument.
589 const Graph* graph_;
590
591 // An array representing the supply (if > 0) or the demand (if < 0)
592 // for each node in graph_.
593 std::unique_ptr<FlowQuantity[]> node_excess_;
594
595 // An array representing the potential (or price function) for
596 // each node in graph_.
597 std::unique_ptr<CostValue[]> node_potential_;
598
599 // An array representing the residual_capacity for each arc in graph_.
600 // Residual capacities enable one to represent the capacity and flow for all
601 // arcs in the graph in the following manner.
602 // For all arcs, residual_arc_capacity_[arc] = capacity[arc] - flow[arc]
603 // Moreover, for reverse arcs, capacity[arc] = 0 by definition.
604 // Also flow[Opposite(arc)] = -flow[arc] by definition.
605 // Therefore:
606 // - for a direct arc:
607 // flow[arc] = 0 - flow[Opposite(arc)]
608 // = capacity[Opposite(arc)] - flow[Opposite(arc)]
609 // = residual_arc_capacity_[Opposite(arc)]
610 // - for a reverse arc:
611 // flow[arc] = -residual_arc_capacity_[arc]
612 // Using these facts enables one to only maintain residual_arc_capacity_,
613 // instead of both capacity and flow, for each direct and indirect arc. This
614 // reduces the amount of memory for this information by a factor 2.
615 // Note that the sum of the largest capacity of an arc in the graph and of
616 // the total flow in the graph mustn't exceed the largest 64 bit integer
617 // to avoid errors. CheckInputConsistency() verifies this constraint.
618 ZVector<ArcFlowType> residual_arc_capacity_;
619
620 // An array representing the first admissible arc for each node in graph_.
621 std::unique_ptr<ArcIndex[]> first_admissible_arc_;
622
623 // A stack used for managing active nodes in the algorithm.
624 // Note that the papers cited above recommend the use of a queue, but
625 // benchmarking so far has not proved it is better.
626 std::stack<NodeIndex> active_nodes_;
627
628 // epsilon_ is the tolerance for optimality.
629 CostValue epsilon_ = 0;
630
631 // alpha_ is the factor by which epsilon_ is divided at each iteration of
632 // Refine().
633 const int64_t alpha_;
634
635 // cost_scaling_factor_ is the scaling factor for cost.
636 CostValue cost_scaling_factor_ = 1;
637
638 // Our node potentials should be >= this at all time.
639 // The first time a potential cross this, BAD_COST_RANGE status is returned.
640 CostValue overflow_threshold_;
641
642 // An array representing the scaled unit cost for each arc in graph_.
643 ZVector<ArcScaledCostType> scaled_arc_unit_cost_;
644
645 // The status of the problem.
646 Status status_ = NOT_SOLVED;
647
648 // An array containing the initial excesses (i.e. the supplies) for each
649 // node. This is used to create the max-flow-based feasibility checker.
650 std::unique_ptr<FlowQuantity[]> initial_node_excess_;
651
652 // Statistics about this class.
653 StatsGroup stats_;
654
655 // Number of Relabel() since last UpdatePrices().
656 int num_relabels_since_last_price_update_ = 0;
657
658 // A Boolean which is true when feasibility has been checked.
659 bool feasibility_checked_ = false;
660
661 // Whether to use the UpdatePrices() heuristic.
662 bool use_price_update_ = false;
663
664 // Whether to check the problem feasibility with a max-flow.
665 bool check_feasibility_ = false;
666
667 // Whether to scale prices, see SimpleMinCostFlow::SetPriceScaling().
668 bool scale_prices_ = true;
669};
670
671#if !SWIG
672
673// Note: SWIG does not seem to understand explicit template specialization and
674// instantiation declarations.
675
678extern template class GenericMinCostFlow<
679 ::util::ReverseArcStaticGraph<uint16_t, int32_t>>;
680extern template class GenericMinCostFlow<
681 ::util::ReverseArcListGraph<int64_t, int64_t>, int64_t, int64_t>;
682extern template class GenericMinCostFlow<
683 ::util::ReverseArcStaticGraph<uint16_t, int32_t>,
684 /*ArcFlowType=*/int16_t,
685 /*ArcScaledCostType=*/int32_t>;
686
687// TODO(b/385094969): Remove this alias after 2025-07-01 to give or-tools users
688// a grace period.
690 template <typename = void>
692 LOG(FATAL) << "MinCostFlow is deprecated. Use `SimpleMinCostFlow` or "
693 "`GenericMinCostFlow` with a specific graph type instead.";
694 }
695};
696
697#endif // SWIG
698
699} // namespace operations_research
700#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::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.
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:1718
NodeIndexType Head(ArcIndexType arc) const
Definition graph.h:1710
In SWIG mode, we don't want anything besides these top-level includes.
BlossomGraph::CostValue CostValue
BlossomGraph::NodeIndex NodeIndex
util::ReverseArcStaticGraph Graph
Type of graph to use.