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"
22 #include "ortools/base/mathutil.h"
23 #include "ortools/graph/graph.h"
24 #include "ortools/graph/graphs.h"
25 #include "ortools/graph/max_flow.h"
26 
27 // TODO(user): Remove these flags and expose the parameters in the API.
28 // New clients, please do not use these flags!
29 ABSL_FLAG(int64, min_cost_flow_alpha, 5,
30  "Divide factor for epsilon at each refine step.");
31 ABSL_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.");
35 ABSL_FLAG(bool, min_cost_flow_check_balance, true,
36  "Check that the sum of supplies is equal to the sum of demands.");
37 ABSL_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");
41 ABSL_FLAG(bool, min_cost_flow_check_result, true,
42  "Check that the result is valid.");
43 
44 namespace operations_research {
45 
46 template <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 
89 template <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 
99 template <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 
109 template <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 
148 template <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 
160 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
161 bool 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 
191 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
192 bool 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 
222 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
223 bool 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 
248 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
249 bool 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 
266 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
267 std::string
268 GenericMinCostFlow<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 
291 template <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);
363  }
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);
373  }
374  }
375  feasibility_checked_ = true;
376  return optimal_max_flow == total_supply;
377 }
378 
379 template <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);
388  }
389  return true;
390 }
391 
392 template <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 }
401 
402 // We use the equations given in the comment of residual_arc_capacity_.
403 template <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 
414 template <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 
422 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
424  NodeIndex node) const {
425  DCHECK(graph_->IsNodeValid(node));
426  return node_excess_[node];
427 }
428 
429 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
432  NodeIndex node) const {
433  return initial_node_excess_[node];
434 }
435 
436 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
439  NodeIndex node) const {
440  return feasible_node_excess_[node];
441 }
442 
443 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
445  ArcIndex arc) const {
446  return FastIsAdmissible(arc, node_potential_[Tail(arc)]);
447 }
448 
449 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
450 bool 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 
457 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
458 bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsActive(
459  NodeIndex node) const {
460  return node_excess_[node] > 0;
461 }
462 
463 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
464 CostValue
465 GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::ReducedCost(
466  ArcIndex arc) const {
467  return FastReducedCost(arc, node_potential_[Tail(arc)]);
468 }
469 
470 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
471 CostValue
472 GenericMinCostFlow<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 
483 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
485 GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::
486  GetFirstOutgoingOrOppositeIncomingArc(NodeIndex node) const {
487  OutgoingOrOppositeIncomingArcIterator arc_it(*graph_, node);
488  return arc_it.Index();
489 }
490 
491 template <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 
532 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
533 void 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 
541 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
542 void 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 
558 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
559 void 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 
569 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
570 void 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 
584 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
585 void 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 
609 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
610 void GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::PushFlow(
611  FlowQuantity flow, ArcIndex arc) {
612  SCOPED_TIME_STAT(&stats_);
613  FastPushFlow(flow, arc, Tail(arc));
614 }
615 
616 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
617 void 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 
634 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
635 void 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 
646 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
647 void 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 
797 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
798 void 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 
819 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
820 void 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 
854 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
855 bool 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 
878 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
879 void 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 
960 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
961 typename Graph::ArcIndex
962 GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::Opposite(
963  ArcIndex arc) const {
964  return Graphs<Graph>::OppositeArc(*graph_, arc);
965 }
966 
967 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
968 bool GenericMinCostFlow<Graph, ArcFlowType, ArcScaledCostType>::IsArcValid(
969  ArcIndex arc) const {
970  return Graphs<Graph>::IsArcValid(*graph_, arc);
971 }
972 
973 template <typename Graph, typename ArcFlowType, typename ArcScaledCostType>
974 bool 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?
984 template class GenericMinCostFlow<StarGraph>;
985 template class GenericMinCostFlow<::util::ReverseArcListGraph<>>;
986 template class GenericMinCostFlow<::util::ReverseArcStaticGraph<>>;
987 template class GenericMinCostFlow<::util::ReverseArcMixedGraph<>>;
988 template class GenericMinCostFlow<::util::ReverseArcStaticGraph<uint16, int32>>;
989 
990 // A more memory-efficient version for large graphs.
991 template 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 
1016  NodeIndex head,
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 
1028 ArcIndex SimpleMinCostFlow::PermutedArc(ArcIndex arc) {
1029  return arc < arc_permutation_.size() ? arc_permutation_[arc] : arc;
1030 }
1031 
1032 SimpleMinCostFlow::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 
1150 CostValue SimpleMinCostFlow::OptimalCost() const { return optimal_cost_; }
1151 
1152 FlowQuantity SimpleMinCostFlow::MaximumFlow() const { return maximum_flow_; }
1153 
1155  return arc_flow_[arc];
1156 }
1157 
1158 NodeIndex SimpleMinCostFlow::NumNodes() const { return node_supply_.size(); }
1159 
1160 ArcIndex SimpleMinCostFlow::NumArcs() const { return arc_tail_.size(); }
1161 
1162 ArcIndex SimpleMinCostFlow::Tail(ArcIndex arc) const { return arc_tail_[arc]; }
1163 
1164 ArcIndex 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 
1178 void 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