always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

graph cleanup

Browse Source
This commit is contained in:
Laurent Perron
2025-01-14 18:27:09 +01:00
parent 59a9fa179e
commit 9beb040c68
10 changed files with 3 additions and 1469 deletions
Download Patch File Download Diff File Expand all files Collapse all files

29
ortools/graph/BUILD.bazel
View File

@@ -528,35 +528,6 @@ cc_test(
],
)
# Min Cost Flow with floating point flow.
cc_library(
name = "fp_min_cost_flow",
srcs = ["fp_min_cost_flow.cc"],
hdrs = ["fp_min_cost_flow.h"],
deps = [
":min_cost_flow",
"//ortools/util:fp_roundtrip_conv",
"//ortools/util:saturated_arithmetic",
"@com_google_absl//absl/algorithm:container",
"@com_google_absl//absl/log",
"@com_google_absl//absl/log:check",
],
)
# cc_test(
# name = "fp_min_cost_flow_test",
# srcs = ["fp_min_cost_flow_test.cc"],
# deps = [
# ":fp_min_cost_flow",
# ":min_cost_flow",
# "//absl/base:no_destructor",
# "//ortools/base:gmock_main",
# "//ortools/util:fp_roundtrip_conv",
# "@com_google_absl//absl/base:log_severity",
# "@com_google_absl//absl/strings",
# ],
# )
# Flow-problem solver
cc_binary(
name = "solve_flow_model",

1
ortools/graph/CMakeLists.txt
View File

@@ -43,7 +43,6 @@ target_link_libraries(${NAME} PRIVATE
if(BUILD_TESTING)
file(GLOB _TEST_SRCS "*_test.cc")
list(FILTER _TEST_SRCS EXCLUDE REGEX "max_flow_test.cc")
list(FILTER _TEST_SRCS EXCLUDE REGEX "fp_min_cost_flow_test.cc")
foreach(_FULL_FILE_NAME IN LISTS _TEST_SRCS)
get_filename_component(_NAME ${_FULL_FILE_NAME} NAME_WE)
get_filename_component(_FILE_NAME ${_FULL_FILE_NAME} NAME)

373
ortools/graph/fp_min_cost_flow.cc
View File

@@ -1,373 +0,0 @@
// Copyright 2010-2025 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ortools/graph/fp_min_cost_flow.h"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <limits>
#include <optional>
#include <ostream>
#include <type_traits>
#include <vector>
#include "absl/algorithm/container.h"
#include "absl/log/check.h"
#include "ortools/base/logging.h"
#include "ortools/graph/min_cost_flow.h"
#include "ortools/util/fp_roundtrip_conv.h"
#include "ortools/util/saturated_arithmetic.h"
namespace operations_research {
namespace {
constexpr SimpleMinCostFlow::FlowQuantity kMaxFlowQuantity =
std::numeric_limits<SimpleMinCostFlow::FlowQuantity>::max();
constexpr SimpleFloatingPointMinCostFlow::FPFlowQuantity kMaxFPFlow =
std::numeric_limits<SimpleFloatingPointMinCostFlow::FPFlowQuantity>::max();
using ::operations_research::internal::ScaleFlow;
// Returns the scaling value computed from log2_scale.
double Scale(const int log2_scale) { return std::ldexp(1.0, log2_scale); }
// Returns the inverse of the scaling value computed from log2_scale.
double InvScale(const int log2_scale) { return std::ldexp(1.0, -log2_scale); }
// Returns true if the max in-flow or the max out-flow are >= kMaxFlowQuantity.
bool AreInOrOutFlowsOverflowing(const SimpleMinCostFlow& min_cost_flow) {
const SimpleMinCostFlow::NodeIndex num_nodes = min_cost_flow.NumNodes();
const SimpleMinCostFlow::ArcIndex num_arcs = min_cost_flow.NumArcs();
if (num_nodes == 0) {
return false;
}
std::vector<SimpleMinCostFlow::FlowQuantity> max_node_in_flow(num_nodes);
std::vector<SimpleMinCostFlow::FlowQuantity> max_node_out_flow(num_nodes);
for (SimpleMinCostFlow::NodeIndex node = 0; node < num_nodes; ++node) {
SimpleMinCostFlow::FlowQuantity supply = min_cost_flow.Supply(node);
if (supply < 0) {
max_node_in_flow[node] = -supply;
} else {
max_node_out_flow[node] = supply;
}
}
for (SimpleMinCostFlow::ArcIndex arc = 0; arc < num_arcs; ++arc) {
const SimpleMinCostFlow::NodeIndex head = min_cost_flow.Head(arc);
const SimpleMinCostFlow::NodeIndex tail = min_cost_flow.Tail(arc);
const SimpleMinCostFlow::FlowQuantity capacity =
min_cost_flow.Capacity(arc);
static_assert(std::is_same_v<SimpleMinCostFlow::FlowQuantity, int64_t>,
"CapAdd() only exists for int64_t");
max_node_in_flow[head] = CapAdd(max_node_in_flow[head], capacity);
max_node_out_flow[tail] = CapAdd(max_node_out_flow[tail], capacity);
}
// It is safe to assume absl::c_max_element() returns a valid iterator as we
// exit this function immediately when `num_nodes == 0`.
const SimpleMinCostFlow::FlowQuantity max_nodes_in_or_out_flow =
std::max(*absl::c_max_element(max_node_in_flow),
*absl::c_max_element(max_node_out_flow));
return max_nodes_in_or_out_flow == kMaxFlowQuantity;
}
} // namespace
SimpleFloatingPointMinCostFlow::SimpleFloatingPointMinCostFlow(
const NodeIndex reserve_num_nodes, const ArcIndex reserve_num_arcs)
: integer_flow_(/*reserve_num_nodes=*/reserve_num_nodes,
/*reserve_num_arcs=*/reserve_num_arcs) {
if (reserve_num_nodes > 0) {
node_supply_.reserve(reserve_num_nodes);
}
if (reserve_num_arcs > 0) {
arc_capacity_.reserve(reserve_num_arcs);
arc_flow_.reserve(reserve_num_arcs);
}
}
SimpleFloatingPointMinCostFlow::ArcIndex
SimpleFloatingPointMinCostFlow::AddArcWithCapacityAndUnitCost(
const NodeIndex tail, const NodeIndex head, const FPFlowQuantity capacity,
const CostValue unit_cost) {
// Add an arc in the integer flow with a temporary capacity of 0. We will
// update it when SolveMaxFlowWithMinCost() is called.
const ArcIndex arc = integer_flow_.AddArcWithCapacityAndUnitCost(
/*tail=*/tail, /*head=*/head,
/*capacity=*/0, /*unit_cost=*/unit_cost);
CHECK_EQ(arc, arc_capacity_.size());
arc_capacity_.push_back(capacity);
arc_flow_.push_back(0.0);
// AddArcWithCapacityAndUnitCost() may have added new nodes based on tail and
// head; we need to take them into account.
const NodeIndex new_num_nodes = integer_flow_.NumNodes();
if (new_num_nodes > node_supply_.size()) {
node_supply_.resize(new_num_nodes);
}
return arc;
}
void SimpleFloatingPointMinCostFlow::SetNodeSupply(
const NodeIndex node, const FPFlowQuantity supply) {
// Set a capacity on the integer flow.
integer_flow_.SetNodeSupply(node, 0);
if (node >= node_supply_.size()) {
node_supply_.resize(node + 1);
}
node_supply_[node] = supply;
}
SimpleFloatingPointMinCostFlow::Status
SimpleFloatingPointMinCostFlow::SolveMaxFlowWithMinCost() {
if (!ScaleSupplyAndCapacity()) {
// Resets the previously computed flow.
absl::c_fill(arc_flow_, 0.0);
return Status::BAD_CAPACITY_RANGE;
}
const Status solve_status = integer_flow_.SolveMaxFlowWithMinCost();
UpdateFlowFromIntegerFlow(solve_status);
return solve_status;
}
std::optional<SimpleFloatingPointMinCostFlow::FPFlowQuantity>
SimpleFloatingPointMinCostFlow::ComputeMaxInOrOutFlow() {
const NodeIndex num_nodes = integer_flow_.NumNodes();
const ArcIndex num_arcs = integer_flow_.NumArcs();
CHECK_EQ(num_nodes, node_supply_.size());
CHECK_EQ(num_arcs, arc_capacity_.size());
if (num_nodes == 0) {
return 0.0;
}
// Compute the max in-flow and max out-flow for each node.
std::vector<FPFlowQuantity> max_node_in_flow(num_nodes, 0.0);
std::vector<FPFlowQuantity> max_node_out_flow(num_nodes, 0.0);
for (NodeIndex node = 0; node < num_nodes; ++node) {
const FPFlowQuantity node_supply = node_supply_[node];
if (!std::isfinite(node_supply)) {
LOG(ERROR) << "Node " << node << " supply is not finite: " << node_supply;
return std::nullopt;
}
if (node_supply < 0.0) {
// Negative supply is demand, thus an input.
max_node_in_flow[node] = -node_supply;
} else {
max_node_out_flow[node] = node_supply;
}
}
for (ArcIndex arc = 0; arc < num_arcs; ++arc) {
const FPFlowQuantity arc_capacity = arc_capacity_[arc];
if (!std::isfinite(arc_capacity)) {
LOG(ERROR) << "Arc " << arc
<< " capacity is not finite: " << arc_capacity;
return std::nullopt;
}
// Negative capacities are considered zero.
if (arc_capacity <= 0.0) {
continue;
}
max_node_in_flow[integer_flow_.Head(arc)] += arc_capacity_[arc];
max_node_out_flow[integer_flow_.Tail(arc)] += arc_capacity_[arc];
}
// We already exited when num_nodes == 0 so we know that returned iterators
// are valid.
return std::max(*absl::c_max_element(max_node_in_flow),
*absl::c_max_element(max_node_out_flow));
}
bool SimpleFloatingPointMinCostFlow::ScaleSupplyAndCapacity() {
const NodeIndex num_nodes = integer_flow_.NumNodes();
const ArcIndex num_arcs = integer_flow_.NumArcs();
CHECK_EQ(num_nodes, node_supply_.size());
CHECK_EQ(num_arcs, arc_capacity_.size());
// Compute the scaling factor for flows.
//
// We choose to use the largest scaling that would not produce an integer
// overflow when solving integer_flow_.
//
// We could use a smaller scaling though, as long as it would not lose any
// non-zero bits. This would be a bit more complex though. On top of that
// always using the largest value may help finding overflow bugs even with
// simple test data.
//
// We want to make sure that the resulting integer flow does not produce a
// BAD_CAPACITY_RANGE. To do so we must make sure that for each node:
// * the sum of incoming arcs capacities + the max(0, node_supply),
// * and the sum of outgoing arcs capacities + the max(0, -node_supply)
// are less (<) than the max value of FlowQuantity.
//
// We thus compute these maximum values using floating point computations and
// use them to compute a scaling factor. Since floating point computations are
// rounded the end result may not be correct and the integer sum may still
// overflow. When that is the case we simply divide the scale by 2 and
// retry. We could do better by changing the rounding mode of floating point
// computation to FE_UPWARD instead of FE_TONEAREST to get an upper-bound on
// the sums but this requires the compiler to properly handle changing
// rounding modes, which is not reliable.
num_tested_scales_ = 0; // Always reset.
if (num_nodes == 0) {
// If we don't have nodes, we don't have arcs either. Thus we have nothing
// to scale.
log2_scale_ = 0;
return true;
}
const std::optional<FPFlowQuantity> max_nodes_in_or_out_flow =
ComputeMaxInOrOutFlow();
if (!max_nodes_in_or_out_flow.has_value()) {
// A LOG(ERROR) already occurred in ComputeMaxInOrOutFlow().
return false;
}
if (!std::isfinite(max_nodes_in_or_out_flow.value())) {
// Note that here we could scale down the floating point values to make the
// sum not overflow. But in practice the caller should avoid this situation.
LOG(ERROR) << "The computed max node in or out flow is not finite: "
<< max_nodes_in_or_out_flow.value();
return false;
}
// Compute the initial scale based on the max in or out flow over all nodes.
// We want that:
//
// static_cast<SimpleMinCostFlow::FlowQuantity>(
// std::round(Scale(log2_scale_) * max_nodes_in_or_out_flow))
// < kMaxFlowQuantity
//
// To make scaling and unscaling more precise, we will only use power-of-two
// values for scale_. Thus we compute p such that:
//
// 2^p <= kMaxFlowQuantity / max_nodes_in_or_out_flow
//
// `f = std::frexp(x, &p)` returns the floating point number of the form `f *
// 2^e` such that:
//
// 2^(e-1) ≤ x < 2^e
//
// Thus we can use it to compute a good starting candidate for the scaling
// factor. Note though that since the division and the computation of
// max_nodes_in_or_out_flow have been subject to rounding, this starting value
// may still lead to overflow of the scaled values. We will thus have a loop
// to try to lower p until we find a value of the scale that works.
if (max_nodes_in_or_out_flow == 0.0) {
// When we have no flow on any node, we can simply scale with 2^0 = 1.
log2_scale_ = 0;
} else {
// Note that if max_nodes_in_or_out_flow is a very small number (< 2^-960)
// then the following division can overflow. If this is the case we simply
// replace the scale by the max possible value.
const double scale_upper_bound =
std::min(std::numeric_limits<double>::max(),
static_cast<double>(kMaxFlowQuantity) /
max_nodes_in_or_out_flow.value());
const double f = std::frexp(scale_upper_bound, &log2_scale_);
// The result of std::frexp() is such that:
// 2^(p-1) <= scale_upper_bound < 2^p
// we thus need to decrement the resulting value to get the lower bound.
//
// When f == 0.5, this means that actually scale_upper_bound == 2^(p-1). In
// that case we know that using this value as the scale will lead to an
// overflow so we simply decrement it by 2 instead.
log2_scale_ -= f == 0.5 ? 2 : 1;
}
// Iterate on values of p until we find one that does not overflow in
// integers. We use saturated_arithmetic.h and detect the issue when we
// overflow FlowQuantity.
//
// Note that we don't expect to do more than two loops usually. If we reach
// cases where we need more loops, then maybe we should a binary search
// approach here. Note that in any case we do a maximum of ~2000 loops (from
// the highest representable power-of-two to the smallest one).
const int initial_log2_scale = log2_scale_;
// We stop when the scale inverse is not representable anymore in a `double`
// (which occurs when we reach denormal number, a.k.a. very close to zero).
for (/*empty*/; std::isfinite(InvScale(log2_scale_)); --log2_scale_) {
const double scale = Scale(log2_scale_);
++num_tested_scales_;
// Sets the integer supply quantities and capacities.
for (NodeIndex node = 0; node < num_nodes; ++node) {
integer_flow_.SetNodeSupply(node, ScaleFlow(node_supply_[node], scale));
}
for (ArcIndex arc = 0; arc < num_arcs; ++arc) {
// Note that here it is safe to use std::max() as we have tested above
// that capacities are not NaN.
integer_flow_.SetArcCapacity(
arc, ScaleFlow(std::max(0.0, arc_capacity_[arc]), scale));
}
// Test the loop end condition.
if (!AreInOrOutFlowsOverflowing(integer_flow_)) {
return true;
}
VLOG(1) << "scale = " << RoundTripDoubleFormat(scale) << " (i.e. 2^"
<< log2_scale_
<< ") lead to an integer overflow; decrementing log2_scale_ and "
"trying again";
}
// It may not be possible to reach this code. If we ever do, we still want
// to treat this as an error.
LOG(ERROR) << "Failed to compute a positive scale that works; started with "
"log2_scale_ = "
<< initial_log2_scale
<< " and stopped at log2_scale_ = " << log2_scale_
<< " with scale_ = " << RoundTripDoubleFormat(Scale(log2_scale_))
<< " 1.0/scale_ = "
<< RoundTripDoubleFormat(InvScale(log2_scale_));
return false;
}
void SimpleFloatingPointMinCostFlow::UpdateFlowFromIntegerFlow(
const Status solve_status) {
switch (solve_status) {
case Status::OPTIMAL:
case Status::FEASIBLE: {
CHECK_EQ(integer_flow_.NumArcs(), arc_flow_.size());
// ScaleSupplyAndCapacity() only selects log2_scale_ values for which
// InvScale() is finite.
const double inv_scale = InvScale(log2_scale_);
for (ArcIndex arc = 0; arc < integer_flow_.NumArcs(); ++arc) {
arc_flow_[arc] = inv_scale * integer_flow_.Flow(arc);
}
break;
}
default:
// SimpleMinCostFlow::arc_flow_ is usually not set in errors cases so we
// simply reset the flow.
absl::c_fill(arc_flow_, 0.0);
break;
}
}
SimpleFloatingPointMinCostFlow::SolveStats
SimpleFloatingPointMinCostFlow::LastSolveStats() const {
return {
.scale = Scale(log2_scale_),
.num_tested_scales = num_tested_scales_,
};
}
std::ostream& operator<<(
std::ostream& out,
const SimpleFloatingPointMinCostFlow::SolveStats& stats) {
out << "{ scale: " << RoundTripDoubleFormat(stats.scale)
<< ", num_tested_scales: " << stats.num_tested_scales << " }";
return out;
}
} // namespace operations_research

226
ortools/graph/fp_min_cost_flow.h
View File

@@ -1,226 +0,0 @@
// Copyright 2010-2025 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef OR_TOOLS_GRAPH_FP_MIN_COST_FLOW_H_
#define OR_TOOLS_GRAPH_FP_MIN_COST_FLOW_H_
#include <cmath>
#include <limits>
#include <optional>
#include <ostream>
#include <vector>
#include "absl/log/check.h"
#include "ortools/graph/min_cost_flow.h"
namespace operations_research {
// An implementation of an approximate min-cost-max-flow algorithm supporting
// floating point numbers for flow capacities.
//
// It uses the integer algorithm of SimpleMinCostFlow by scaling and rounding
// floating point supply quantities and capacities to make them fit on
// integers. This can be seen as using fixed-point arithmetic.
//
// Note that this only supports min-cost-max-flow and not min-cost-flow. The
// reason is that with floating point numbers it is harder to define that all
// demand and supply are met. This would require some tolerances and testing
// those tolerances would require solving the max-flow anyway.
class SimpleFloatingPointMinCostFlow {
public:
using NodeIndex = SimpleMinCostFlow::NodeIndex;
using ArcIndex = SimpleMinCostFlow::ArcIndex;
using CostValue = SimpleMinCostFlow::CostValue;
using FPFlowQuantity = double;
using Status = SimpleMinCostFlow::Status;
// Statistics associated with a call to SolveMaxFlowWithMinCost().
//
// Returned by LastSolveStats().
struct SolveStats {
// The scaling factor used to convert the FPFlowQuantity to a
// SimpleMinCostFlow::FlowQuantity.
double scale = 1.0;
// The number of values tested for the scaling factor.
//
// Internally SolveMaxFlowWithMinCost() will compute a first scaling factor
// with floating point arithmetic. Due to the approximate nature of this
// computation, it may still be too high. If the resulting integer numbers
// lead to integer overflows, a new lower scaling factor is computing.
int num_tested_scales = 0;
// Prints all the stats values on a single line.
friend std::ostream& operator<<(std::ostream& out, const SolveStats& stats);
};
explicit SimpleFloatingPointMinCostFlow(NodeIndex reserve_num_nodes = 0,
ArcIndex reserve_num_arcs = 0);
// This type is neither copyable nor movable.
SimpleFloatingPointMinCostFlow(const SimpleFloatingPointMinCostFlow&) =
delete;
SimpleFloatingPointMinCostFlow& operator=(
const SimpleFloatingPointMinCostFlow&) = delete;
// Adds a directed arc from tail to head to the underlying graph with
// a given capacity and cost per unit of flow.
// * Node indices must be non-negative (>= 0).
// * The capacity must be finite. When not, SolveMaxFlowWithMinCost() returns
// BAD_CAPACITY_RANGE. Negative values are OK and will be considered zero
// (which is useful when computed values are close to zero but negative).
// * The unit cost can take any integer value (even negative).
// * Self-looping and duplicate arcs are supported.
// * After the method finishes, NumArcs() == the returned ArcIndex + 1.
ArcIndex AddArcWithCapacityAndUnitCost(NodeIndex tail, NodeIndex head,
FPFlowQuantity capacity,
CostValue unit_cost);
// Sets the supply of the given node. The node index must be non-negative (>=
// 0). Nodes implicitly created will have a default supply set to 0. A demand
// is modeled as a negative supply.
//
// The supply quantity must be finite. When not, SolveMaxFlowWithMinCost()
// returns BAD_CAPACITY_RANGE.
void SetNodeSupply(NodeIndex node, FPFlowQuantity supply);
// Computes a maximum-flow with minimum cost.
//
// This function returns the Status of the underlying
// SimpleMinCostFlow::SolveMaxFlowWithMinCost().
//
// It also returns BAD_CAPACITY_RANGE:
// * either when the arc capacities or node supply quantities are NaN or
// infinite,
// * or when the computed in-flow or out-flow of a node results in an infinite
// value,
// * or if it failed to find a scale factor to make capacities and supply
// quantities fit in integers.
// In case of failure, a LOG(ERROR) should contain the explanation.
Status SolveMaxFlowWithMinCost();
// Returns the flow on arc, this only make sense for a successful
// SolveMaxFlowWithMinCost().
//
// If called before SolveMaxFlowWithMinCost(), it will return 0.0.
//
// Note: It is possible that there is more than one optimal solution. The
// algorithm is deterministic so it will always return the same solution for
// a given problem. However, there is no guarantee of this from one code
// version to the next (but the code does not change often).
FPFlowQuantity Flow(ArcIndex arc) const { return arc_flow_[arc]; }
// Returns the statistics of the last call to SolveMaxFlowWithMinCost().
SolveStats LastSolveStats() const;
// Accessors for the user given data. The implementation will crash if "arc"
// is not in [0, NumArcs()) or "node" is not in [0, NumNodes()).
NodeIndex NumNodes() const { return integer_flow_.NumNodes(); }
ArcIndex NumArcs() const { return integer_flow_.NumArcs(); }
NodeIndex Tail(ArcIndex arc) const { return integer_flow_.Tail(arc); }
NodeIndex Head(ArcIndex arc) const { return integer_flow_.Head(arc); }
FPFlowQuantity Capacity(ArcIndex arc) const { return arc_capacity_[arc]; }
FPFlowQuantity Supply(NodeIndex node) const { return node_supply_[node]; }
CostValue UnitCost(ArcIndex arc) const { return integer_flow_.UnitCost(arc); }
private:
// Returns the max value all in-flows or out-flows of each nodes. If some
// nodes or arcs have non-finite supply or capacities, returns std::nullopt
// after the rationale has been printed to LOG(ERROR).
//
// Precisely, returns max(max(in_flow(n) ∀ n), max(out_flow(n) ∀ n)). Returns
// 0.0 when there are no nodes.
std::optional<FPFlowQuantity> ComputeMaxInOrOutFlow();
// Sets the integer supply and capacity values on integer_flow_ from
// node_supply_ and arc_capacity_ floating point values after computing the
// log2_scale_. It will also updates num_tested_scales_.
//
// The integer quantities are computed by:
//
// static_cast<FlowQuantity>(
// std::round(Scale(log2_scale_) * fp_flow_quantity));
//
// Returns true in case of success, false otherwise (after the rational was
// printed to LOG(ERROR)).
bool ScaleSupplyAndCapacity();
// Updates arc_flow_ using the values of integer_flow_ and the status of the
// solve.
void UpdateFlowFromIntegerFlow(Status solve_status);
// The underlying integer min-cost-max-flow solver.
SimpleMinCostFlow integer_flow_;
// The log2 of the scale applied to FPFlowQuantity values to get the integer
// FlowQuantity ones in integer_flow_. When ScaleSupplyAndCapacity() succeeds
// this will contain a value for which both Scale() and InvScale() are finite
// and non-zero.
//
// See Scale() and InvScale().
int log2_scale_ = 0;
// The number of values of scale_ tested during the clal to
// ScaleSupplyAndCapacity().
int num_tested_scales_ = 0;
// Invariant on the following vectors: their size matches
// integer_flow_.NumNodes() or integer_flow_.NumArcs().
std::vector<FPFlowQuantity> arc_capacity_;
std::vector<FPFlowQuantity> node_supply_;
std::vector<FPFlowQuantity> arc_flow_;
};
// This namespace contains internal functions only there for tests.
namespace internal {
// Scale the input floating point flow to the integer flow, clamping the
// result to [-kMaxFlowQuantity, kMaxFlowQuantity].
//
// The scale and fp_flow must be finite (CHECKed).
//
// This function is internal. It is only exposed since by construction the input
// fp_flow and scale of ScaleFlow() should never trigger the code paths that
// lead to overflows. But we still want to make sure that if they are ever
// triggered they work as expected.
inline SimpleMinCostFlow::FlowQuantity ScaleFlow(
const SimpleFloatingPointMinCostFlow::FPFlowQuantity fp_flow,
const double scale) {
CHECK(std::isfinite(scale)) << scale;
CHECK(std::isfinite(fp_flow)) << fp_flow;
const double rounded_scaled_flow = std::round(scale * fp_flow);
// Note that we compare to kMaxFlowQuantity with >= and not > as:
// * the comparison will convert kMaxFlowQuantity to `double` first (see
// https://en.cppreference.com/w/cpp/language/usual_arithmetic_conversions),
// * and kMaxFlowQuantity (2^63 - 1) is not exactly representable in a
// `double` (that only has 53 bits of mantissa),
// * thus it will be rounded to the nearest `double`, i.e. 2^64,
// * so if we were to compare with > only, we would not reject the `double`
// 2^64 that can't fit in an int64_t; comparing with >= will reject this
// number and only accept `double` that are less or equal to the predecessor
// of 2^64, which all fit in an int64_t.
constexpr SimpleMinCostFlow::FlowQuantity kMaxFlowQuantity =
std::numeric_limits<SimpleMinCostFlow::FlowQuantity>::max();
if (rounded_scaled_flow >= kMaxFlowQuantity) {
return kMaxFlowQuantity;
}
if (rounded_scaled_flow <= -kMaxFlowQuantity) {
return -kMaxFlowQuantity;
}
return static_cast<SimpleMinCostFlow::FlowQuantity>(rounded_scaled_flow);
}
} // namespace internal
} // namespace operations_research
#endif // OR_TOOLS_GRAPH_FP_MIN_COST_FLOW_H_

604
ortools/graph/fp_min_cost_flow_test.cc
View File

@@ -1,604 +0,0 @@
// Copyright 2010-2025 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ortools/graph/fp_min_cost_flow.h"
#include <cmath>
#include <limits>
#include <ostream>
#include <sstream>
#include <string>
#include <vector>
#include "absl/base/log_severity.h"
#include "absl/base/no_destructor.h"
#include "absl/strings/escaping.h"
#include "gtest/gtest.h"
#include "ortools/base/gmock.h"
#include "ortools/graph/min_cost_flow.h"
#include "ortools/util/fp_roundtrip_conv.h"
#include "testing/base/public/mock-log.h"
namespace operations_research {
namespace {
using ::testing::AnyNumber;
using ::testing::DoubleNear;
using ::testing::HasSubstr;
using ::testing::kDoNotCaptureLogsYet;
using ::testing::Mock;
using ::testing::ScopedMockLog;
using Status = SimpleFloatingPointMinCostFlow::Status;
using FPFlowQuantity = SimpleFloatingPointMinCostFlow::FPFlowQuantity;
using CostValue = SimpleFloatingPointMinCostFlow::CostValue;
using NodeIndex = SimpleFloatingPointMinCostFlow::NodeIndex;
using ArcIndex = SimpleFloatingPointMinCostFlow::ArcIndex;
using ::operations_research::internal::ScaleFlow;
constexpr double kInf = std::numeric_limits<double>::infinity();
constexpr double kNaN = std::numeric_limits<double>::quiet_NaN();
constexpr SimpleFloatingPointMinCostFlow::CostValue kMaxCostValue =
std::numeric_limits<SimpleFloatingPointMinCostFlow::CostValue>::max();
const SimpleFloatingPointMinCostFlow::FPFlowQuantity kMaxFPFlow =
std::numeric_limits<SimpleFloatingPointMinCostFlow::FPFlowQuantity>::max();
constexpr SimpleFloatingPointMinCostFlow::FPFlowQuantity kMinFPFlow =
std::numeric_limits<
SimpleFloatingPointMinCostFlow::FPFlowQuantity>::denorm_min();
// Returns the path of the source file that emits LOG(ERROR) in case of failure.
//
// It is meant to be used for log mocking assertions. It needs to be a `const
// std::string&` and can't be an absl::string_view. Thus this function returns a
// reference on a static.
const std::string& SourceFilePath() {
static const absl::NoDestructor<std::string> kModule(
"ortools/graph/fp_min_cost_flow.cc");
return *kModule;
}
TEST(SimpleFloatingPointMinCostFlowTest, Empty) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(), Status::OPTIMAL);
EXPECT_EQ(fp_min_cost_flow.NumNodes(), 0);
EXPECT_EQ(fp_min_cost_flow.NumArcs(), 0);
}
TEST(SimpleFloatingPointMinCostFlowTest, SetNodeSupplyCreateNodes) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
// Setting the supply of node 3 should create four nodes, even if the supply
// is 0.0.
fp_min_cost_flow.SetNodeSupply(/*node=*/3, /*supply=*/0.0);
ASSERT_EQ(fp_min_cost_flow.NumNodes(), 4);
EXPECT_EQ(fp_min_cost_flow.Supply(0), 0.0);
EXPECT_EQ(fp_min_cost_flow.Supply(1), 0.0);
EXPECT_EQ(fp_min_cost_flow.Supply(2), 0.0);
EXPECT_EQ(fp_min_cost_flow.Supply(3), 0.0);
// Setting arbitrary positive and negative values.
fp_min_cost_flow.SetNodeSupply(/*node=*/1, /*supply=*/13.75);
fp_min_cost_flow.SetNodeSupply(/*node=*/0, /*supply=*/-28.125);
ASSERT_EQ(fp_min_cost_flow.NumNodes(), 4);
EXPECT_EQ(fp_min_cost_flow.Supply(0), -28.125);
EXPECT_EQ(fp_min_cost_flow.Supply(1), 13.75);
EXPECT_EQ(fp_min_cost_flow.Supply(2), 0.0);
EXPECT_EQ(fp_min_cost_flow.Supply(3), 0.0);
}
TEST(SimpleFloatingPointMinCostFlowTest, AddArcCreateNodes) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
// Adding an arc between nodes 1 and 3 should create four nodes, with zero
// supply.
const ArcIndex arc0 = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/3, /*head=*/1, /*capacity=*/0.0, /*unit_cost=*/0);
const ArcIndex arc1 = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/0, /*head=*/1, /*capacity=*/15.25, /*unit_cost=*/3);
ASSERT_EQ(fp_min_cost_flow.NumNodes(), 4);
EXPECT_EQ(fp_min_cost_flow.Supply(0), 0.0);
EXPECT_EQ(fp_min_cost_flow.Supply(1), 0.0);
EXPECT_EQ(fp_min_cost_flow.Supply(2), 0.0);
EXPECT_EQ(fp_min_cost_flow.Supply(3), 0.0);
ASSERT_EQ(fp_min_cost_flow.NumArcs(), 2);
EXPECT_EQ(fp_min_cost_flow.Tail(arc0), 3);
EXPECT_EQ(fp_min_cost_flow.Head(arc0), 1);
EXPECT_EQ(fp_min_cost_flow.Capacity(arc0), 0.0);
EXPECT_EQ(fp_min_cost_flow.UnitCost(arc0), 0.0);
EXPECT_EQ(fp_min_cost_flow.Tail(arc1), 0);
EXPECT_EQ(fp_min_cost_flow.Head(arc1), 1);
EXPECT_EQ(fp_min_cost_flow.Capacity(arc1), 15.25);
EXPECT_EQ(fp_min_cost_flow.UnitCost(arc1), 3);
}
// Value-parameterized test of invalid node supply values.
class BadNodeSupplyTest : public testing::TestWithParam<double> {};
TEST_P(BadNodeSupplyTest, FailingSolve) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
fp_min_cost_flow.SetNodeSupply(/*node=*/0, /*supply=*/GetParam());
ScopedMockLog log(kDoNotCaptureLogsYet);
EXPECT_CALL(log, Log).Times(AnyNumber()); // Ignore unexpected logs.
EXPECT_CALL(log, Log(base_logging::ERROR, SourceFilePath(),
HasSubstr("supply is not finite")));
log.StartCapturingLogs();
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(),
SimpleFloatingPointMinCostFlow::Status::BAD_CAPACITY_RANGE);
log.StopCapturingLogs();
Mock::VerifyAndClearExpectations(&log);
}
// Tests that failures properly resets computed arc flows before returning.
TEST_P(BadNodeSupplyTest, FailureResetsArcFlows) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
// First make a successful solve.
const ArcIndex arc = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/0, /*head=*/1, /*capacity=*/5.25, /*unit_cost=*/0);
fp_min_cost_flow.SetNodeSupply(0, 10.0);
fp_min_cost_flow.SetNodeSupply(1, -10.0);
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(),
SimpleFloatingPointMinCostFlow::Status::OPTIMAL);
EXPECT_EQ(fp_min_cost_flow.Flow(arc), 5.25);
// Then set the invalid supply, the next solve should fail and resets the arc
// flow.
fp_min_cost_flow.SetNodeSupply(/*node=*/0, /*supply=*/GetParam());
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(),
SimpleFloatingPointMinCostFlow::Status::BAD_CAPACITY_RANGE);
EXPECT_EQ(fp_min_cost_flow.Flow(arc), 0.0);
}
INSTANTIATE_TEST_SUITE_P(All, BadNodeSupplyTest,
testing::Values(kInf, -kInf, kNaN));
// Value-parameterized test of invalid arc capacity values.
class BadArcCapacityTest : public testing::TestWithParam<double> {};
TEST_P(BadArcCapacityTest, FailingSolve) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/0, /*head=*/1, /*capacity=*/GetParam(), /*unit_cost=*/0);
ScopedMockLog log(kDoNotCaptureLogsYet);
EXPECT_CALL(log, Log).Times(AnyNumber()); // Ignore unexpected logs.
EXPECT_CALL(log, Log(base_logging::ERROR, SourceFilePath(),
HasSubstr("capacity is not finite")));
log.StartCapturingLogs();
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(),
SimpleFloatingPointMinCostFlow::Status::BAD_CAPACITY_RANGE);
log.StopCapturingLogs();
Mock::VerifyAndClearExpectations(&log);
}
// Tests that failures properly resets computed arc flows before returning.
TEST_P(BadArcCapacityTest, FailureResetsArcFlows) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
// First make a successful solve.
const ArcIndex arc = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/0, /*head=*/1, /*capacity=*/5.25, /*unit_cost=*/0);
fp_min_cost_flow.SetNodeSupply(0, 10.0);
fp_min_cost_flow.SetNodeSupply(1, -10.0);
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(),
SimpleFloatingPointMinCostFlow::Status::OPTIMAL);
EXPECT_EQ(fp_min_cost_flow.Flow(arc), 5.25);
// Then add an arc with the invalid capacity, the next solve should fail and
// resets the arc flow.
const ArcIndex bad_arc = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/0, /*head=*/1, /*capacity=*/GetParam(), /*unit_cost=*/0);
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(),
SimpleFloatingPointMinCostFlow::Status::BAD_CAPACITY_RANGE);
EXPECT_EQ(fp_min_cost_flow.Flow(arc), 0.0);
EXPECT_EQ(fp_min_cost_flow.Flow(bad_arc), 0.0);
}
INSTANTIATE_TEST_SUITE_P(All, BadArcCapacityTest,
testing::Values(kInf, -kInf, kNaN));
// Parameters for the TwoNodesOneArcTest value-parameterized test.
struct TwoNodesOneArcTestParams {
// Name of the test; this will be used as a suffix to the test name and thus
// should be snake_case with no spaces.
std::string name;
FPFlowQuantity head_node_supply = 0.0;
FPFlowQuantity tail_node_supply = 0.0;
FPFlowQuantity arc_capacity = 0.0;
CostValue arc_unit_cost = 0;
// The expected computed flow on the arc; it is tested using
// expected_flow_tolerance absolute tolerance.
FPFlowQuantity expected_flow = 0.0;
FPFlowQuantity expected_flow_tolerance = 0.0;
friend std::ostream& operator<<(std::ostream& out,
const TwoNodesOneArcTestParams& params) {
out << "{ name: '" << absl::CEscape(params.name) << "', head_node_supply: "
<< RoundTripDoubleFormat(params.head_node_supply)
<< ", tail_node_supply: "
<< RoundTripDoubleFormat(params.tail_node_supply)
<< ", arc_capacity: " << RoundTripDoubleFormat(params.arc_capacity)
<< ", arc_unit_cost: " << params.arc_unit_cost
<< ", expected_flow: " << RoundTripDoubleFormat(params.expected_flow)
<< ", expected_flow_tolerance: "
<< RoundTripDoubleFormat(params.expected_flow_tolerance) << " }";
return out;
}
};
// Value-parameterized tests with a graph containing a single arc and two nodes.
class TwoNodesOneArcTest
: public testing::TestWithParam<TwoNodesOneArcTestParams> {};
TEST_P(TwoNodesOneArcTest, SolveMaxFlowWithMinCost) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
const NodeIndex kTailNode = 0;
const NodeIndex kHeadNode = 1;
const ArcIndex arc = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/kTailNode, /*head=*/kHeadNode,
/*capacity=*/GetParam().arc_capacity,
/*unit_cost=*/GetParam().arc_unit_cost);
fp_min_cost_flow.SetNodeSupply(kTailNode, GetParam().tail_node_supply);
fp_min_cost_flow.SetNodeSupply(kHeadNode, GetParam().head_node_supply);
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(), Status::OPTIMAL);
EXPECT_THAT(
fp_min_cost_flow.Flow(arc),
DoubleNear(GetParam().expected_flow, GetParam().expected_flow_tolerance));
}
INSTANTIATE_TEST_SUITE_P(
All, TwoNodesOneArcTest,
testing::Values(
// All parameters at 0. Expect no flow.
TwoNodesOneArcTestParams{.name = "empty"},
// A flow of 1.0 with small integer capacities and supply.
TwoNodesOneArcTestParams{
.name = "small_integers",
.head_node_supply = -1.0,
.tail_node_supply = 2.0,
.arc_capacity = 3.0,
.arc_unit_cost = 0,
.expected_flow = 1.0,
.expected_flow_tolerance = 0.0,
},
// Some big but not huge capacity and supply numbers.
TwoNodesOneArcTestParams{
.name = "big_numbers",
.head_node_supply = -1.0e18,
.tail_node_supply = 2.0e16,
.arc_capacity = 5.0e16,
.arc_unit_cost = 0,
.expected_flow = 2.0e16,
.expected_flow_tolerance = 0.0,
},
// Big capacity and small supply. We use -1.1 so that we have a mantissa
// with a lot of non-zeros (in base-2 -1.1 is -1.0001100110011...).
//
// With a value of 5.3e10, the largest scaling factor so that:
//
// std::round(2^p * 5.3e10) < (2^63 - 1)
//
// is 2**27. With this value, the 1.1 theoretical flow after scaling and
// unscaling will be:
//
// std::round(1.1 * 2.0^27) * 2.0^-27 = 1.1000000014901161
//
// the error is thus:
//
// 1.4901160305669237e-09
//
// We simply use 1.5e-9 as tolerance.
TwoNodesOneArcTestParams{
.name = "big_capacity_small_supply",
.head_node_supply = -1.1,
.tail_node_supply = 2.5,
.arc_capacity = 5.3e10,
.arc_unit_cost = 0,
.expected_flow = 1.1,
.expected_flow_tolerance = 1.5e-9, // See comment above.
},
// Small capacity and big supply. See previous test for rationale.
TwoNodesOneArcTestParams{
.name = "small_capacity_big_supply",
.head_node_supply = -5.3e10,
.tail_node_supply = 4.2e8,
.arc_capacity = 1.1,
.arc_unit_cost = 0,
.expected_flow = 1.1,
.expected_flow_tolerance = 1.0e-5, // See comment above.
},
// tiny capacity (1.5) and huge supply (2^100). We expect scaling to
// result in a 0 integer capacity as the scale factor that makes 2^100
// fit into an int64_t will make the small capacity close to zero.
TwoNodesOneArcTestParams{
.name = "tiny_capacity_huge_supply",
.head_node_supply = -std::pow(2.0, 100.0),
.tail_node_supply = std::pow(2.0, 80.0),
.arc_capacity = 1.5,
.arc_unit_cost = 0,
.expected_flow = 0,
.expected_flow_tolerance = 0.0,
},
// Huge capacity (2^100) and tiny supply values (1.5). See previous test
// for rationale.
TwoNodesOneArcTestParams{
.name = "huge_capacity_tiny_supply",
.head_node_supply = -1.5,
.tail_node_supply = 1.5,
.arc_capacity = std::pow(2.0, 100.0),
.arc_unit_cost = 0,
.expected_flow = 0,
.expected_flow_tolerance = 0.0,
},
// A capacity flow should be considered 0.
TwoNodesOneArcTestParams{
.name = "negative_capacity",
.head_node_supply = -1.0,
.tail_node_supply = 2.0,
.arc_capacity = -3.0,
.arc_unit_cost = 0,
.expected_flow = 0.0,
.expected_flow_tolerance = 0.0,
},
// Test scaling when the numbers are the smallest possible. The computed
// flow is likely 0. We use 4 * kMinFPFlow as a reasonable tolerance.
TwoNodesOneArcTestParams{
.name = "min_denormal_flow",
.head_node_supply = -kMinFPFlow,
.tail_node_supply = kMinFPFlow,
.arc_capacity = kMinFPFlow,
.arc_unit_cost = 0,
.expected_flow = kMinFPFlow,
.expected_flow_tolerance = 4 * kMinFPFlow,
},
// Test when supply are the smallest possible and capacity is zero (so
// that for each node we only have kMinFPFlow in or out flow).
TwoNodesOneArcTestParams{
.name = "min_denormal_supply_no_flow",
.head_node_supply = -kMinFPFlow,
.tail_node_supply = kMinFPFlow,
.arc_capacity = 0.0,
.arc_unit_cost = 0,
.expected_flow = 0.0,
.expected_flow_tolerance = 0.0,
}, ),
[](const testing::TestParamInfo<TwoNodesOneArcTestParams>& info) {
return info.param.name;
});
TEST(SimpleFloatingPointMinCostFlowTest, IntegerSolverFailing) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
// To generate a failure, we use a unit_cost that is too high.
const ArcIndex arc = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/0, /*head=*/1,
/*capacity=*/1.0,
/*unit_cost=*/kMaxCostValue);
fp_min_cost_flow.SetNodeSupply(0, 1.0);
fp_min_cost_flow.SetNodeSupply(1, -1.0);
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(), Status::BAD_COST_RANGE);
EXPECT_EQ(fp_min_cost_flow.Flow(arc), 0.0);
}
TEST(SimpleFloatingPointMinCostFlowTest, IntegerSolverFailingResetsFlow) {
// First make a successful solve.
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
const ArcIndex arc = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/0, /*head=*/1,
/*capacity=*/1.0,
/*unit_cost=*/0);
fp_min_cost_flow.SetNodeSupply(0, 1.0);
fp_min_cost_flow.SetNodeSupply(1, -1.0);
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(), Status::OPTIMAL);
EXPECT_EQ(fp_min_cost_flow.Flow(arc), 1.0);
// Then generate a failure using a too high cost value.
const ArcIndex failing_arc = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/0, /*head=*/1, /*capacity=*/1.0, /*unit_cost=*/kMaxCostValue);
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(), Status::BAD_COST_RANGE);
EXPECT_EQ(fp_min_cost_flow.Flow(arc), 0.0);
EXPECT_EQ(fp_min_cost_flow.Flow(failing_arc), 0.0);
}
TEST(SimpleFloatingPointMinCostFlowTest, FloatingPointSumOverflow) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
// We use two arcs pointing the same node with floating point capacities that
// overflows the `double` exponent range. This will result in an overflow when
// computing the maximum intake of the node which should result in an error.
//
// Note that we could update SimpleFloatingPointMinCostFlow to deal with that
// by scaling the numbers down eventually.
fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/0, /*head=*/1,
/*capacity=*/kMaxFPFlow,
/*unit_cost=*/0);
fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/0, /*head=*/1,
/*capacity=*/kMaxFPFlow,
/*unit_cost=*/0);
ScopedMockLog log(kDoNotCaptureLogsYet);
EXPECT_CALL(log, Log).Times(AnyNumber()); // Ignore unexpected logs.
EXPECT_CALL(
log, Log(base_logging::ERROR, SourceFilePath(), HasSubstr("not finite")));
log.StartCapturingLogs();
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(),
Status::BAD_CAPACITY_RANGE);
log.StopCapturingLogs();
Mock::VerifyAndClearExpectations(&log);
}
TEST(SimpleFloatingPointMinCostFlowTest, FirstScaleFailed) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
// We build a corner case where the evaluation of the in/out-flow of nodes in
// double would lead to a starting scale value that is too high.
//
// The chosen `supply` value is such that adding arc capacities one by one
// (each arc has 1.0 capacity) does not change it. Thus the computed in-flow
// and out-flow of both nodes is `supply`.
//
// On top of that computing the initial scale will lead to a scale of 2, which
// would be OK as 2 * supply is a valid integer supply. The issue is that each
// arc integer capacity is 2 and we have 2048 of them. And now:
//
// 2 * supply + 2 * 2048 > 2^63 - 1
//
// We thus will have to loop and try using a scale of 1.0 which will work.
//
// Note that this test we relies on the implementation adding nodes supply
// first, then arc capacities.
const FPFlowQuantity supply = std::nextafter(std::ldexp(1.0, 62), -kInf);
fp_min_cost_flow.SetNodeSupply(/*node=*/0, supply);
fp_min_cost_flow.SetNodeSupply(/*node=*/1, -supply);
std::vector<ArcIndex> arcs;
for (int i = 0; i < 2048; ++i) {
arcs.push_back(fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/0, /*head=*/1,
/*capacity=*/1.0, /*unit_cost=*/0));
}
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(), Status::OPTIMAL);
EXPECT_EQ(fp_min_cost_flow.LastSolveStats().num_tested_scales, 2);
for (const ArcIndex arc : arcs) {
// Use ASSERT_EQ() here to prevent having too many errors. All arcs are
// equivalent so if there is an error it is likely to happen for each arc
// anyway.
ASSERT_EQ(fp_min_cost_flow.Flow(arc), 1.0) << "arc: " << arc;
}
}
TEST(SimpleFloatingPointMinCostFlowTest, FirstScaleSucceeds) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
// We build a corner case where the computed initial scale is good enough and
// we test that we get the expected result. This happens when the max
// in/out-flow is a power of two.
fp_min_cost_flow.SetNodeSupply(/*node=*/0, 2.0);
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(), Status::OPTIMAL);
EXPECT_EQ(fp_min_cost_flow.LastSolveStats().num_tested_scales, 1);
}
TEST(SimpleFloatingPointMinCostFlowTest, Assignment) {
SimpleFloatingPointMinCostFlow fp_min_cost_flow;
// Tests that costs are taken into account by doing a simple assignment:
//
// cap:5.1 cost:2
// n0:4 ---------------> n3:-5
// ^
// / cap:4
// / cost:1
// n1:3 -------------+
// \ cap:5
// \ cost:2
// v
// n2:5 ---------------> n4:-4
// cap:3 cost:3
//
// Here we expect:
// * n3 to match its demand of 5 from n0 and n1, using n1 first as it has the
// lowest cost, then n0.
// * n4 to match its demand of 4 from n1 and n2, using n1 first as it has the
// lower cost, then n2.
// * n3 to consume from n1 only 2 and not 3 as doing so would starve n4 that
// would not be able to match its demand.
fp_min_cost_flow.SetNodeSupply(/*node=*/0, 4.0);
fp_min_cost_flow.SetNodeSupply(/*node=*/1, 3.0);
fp_min_cost_flow.SetNodeSupply(/*node=*/2, 4.0);
fp_min_cost_flow.SetNodeSupply(/*node=*/3, -5.0);
fp_min_cost_flow.SetNodeSupply(/*node=*/4, -4.0);
const ArcIndex a03 = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/0, /*head=*/3, /*capacity=*/5.1, /*unit_cost=*/2);
const ArcIndex a13 = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/1, /*head=*/3, /*capacity=*/4.0, /*unit_cost=*/1);
const ArcIndex a14 = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/1, /*head=*/4, /*capacity=*/5.0, /*unit_cost=*/2);
const ArcIndex a24 = fp_min_cost_flow.AddArcWithCapacityAndUnitCost(
/*tail=*/2, /*head=*/4, /*capacity=*/3.0, /*unit_cost=*/3);
EXPECT_EQ(fp_min_cost_flow.SolveMaxFlowWithMinCost(), Status::OPTIMAL);
EXPECT_EQ(fp_min_cost_flow.Flow(a03), 3);
EXPECT_EQ(fp_min_cost_flow.Flow(a13), 2);
EXPECT_EQ(fp_min_cost_flow.Flow(a14), 1);
EXPECT_EQ(fp_min_cost_flow.Flow(a24), 3);
}
TEST(ScaleFlowTest, AllValues) {
constexpr SimpleMinCostFlow::FlowQuantity kMaxFlowQuantity =
std::numeric_limits<SimpleMinCostFlow::FlowQuantity>::max();
EXPECT_EQ(ScaleFlow(/*fp_flow=*/kMaxFPFlow,
/*scale=*/std::numeric_limits<double>::max()),
kMaxFlowQuantity);
EXPECT_EQ(ScaleFlow(/*fp_flow=*/-kMaxFPFlow,
/*scale=*/std::numeric_limits<double>::max()),
-kMaxFlowQuantity);
EXPECT_EQ(ScaleFlow(/*fp_flow=*/kMaxFPFlow, /*scale=*/1.0), kMaxFlowQuantity);
EXPECT_EQ(ScaleFlow(/*fp_flow=*/-kMaxFPFlow, /*scale=*/1.0),
-kMaxFlowQuantity);
EXPECT_EQ(ScaleFlow(/*fp_flow=*/2.0 * kMaxFlowQuantity, /*scale=*/1.0),
kMaxFlowQuantity);
EXPECT_EQ(ScaleFlow(/*fp_flow=*/-2.0 * kMaxFlowQuantity, /*scale=*/1.0),
-kMaxFlowQuantity);
EXPECT_EQ(ScaleFlow(/*fp_flow=*/kMaxFlowQuantity, /*scale=*/1.0),
kMaxFlowQuantity);
EXPECT_EQ(ScaleFlow(/*fp_flow=*/-kMaxFlowQuantity, /*scale=*/1.0),
-kMaxFlowQuantity);
// The integer that is the `double` just before the rounded value of
// kMaxFlowQuantity. This rounded value does not fit in an integer but its
// predecessor will.
const FPFlowQuantity kMaxFlowQuantityPredecessor =
std::nextafter(static_cast<FPFlowQuantity>(kMaxFlowQuantity), 0.0);
const SimpleMinCostFlow::FlowQuantity kMaxFlowQuantityPredecessorAsInt =
static_cast<SimpleMinCostFlow::FlowQuantity>(kMaxFlowQuantityPredecessor);
ASSERT_LT(kMaxFlowQuantityPredecessorAsInt, kMaxFlowQuantity);
EXPECT_EQ(ScaleFlow(/*fp_flow=*/kMaxFlowQuantityPredecessor, /*scale=*/1.0),
kMaxFlowQuantityPredecessor);
EXPECT_EQ(ScaleFlow(/*fp_flow=*/-kMaxFlowQuantityPredecessor, /*scale=*/1.0),
-kMaxFlowQuantityPredecessor);
}
TEST(ScaleFlowDeathTest, NaNScale) {
EXPECT_DEATH(ScaleFlow(/*fp_flow=*/1.0, /*scale=*/kNaN), "scale");
}
TEST(ScaleFlowDeathTest, InfScale) {
EXPECT_DEATH(ScaleFlow(/*fp_flow=*/1.0, /*scale=*/kInf), "scale");
}
TEST(ScaleFlowDeathTest, NaNFlow) {
EXPECT_DEATH(ScaleFlow(/*fp_flow=*/kNaN, /*scale=*/1.0), "fp_flow");
}
TEST(ScaleFlowDeathTest, InfFlow) {
EXPECT_DEATH(ScaleFlow(/*fp_flow=*/kInf, /*scale=*/1.0), "fp_flow");
}
TEST(SolveStatsTest, Stream) {
std::ostringstream oss;
oss << SimpleFloatingPointMinCostFlow::SolveStats{
.scale = 0.1,
.num_tested_scales = 4,
};
EXPECT_EQ(oss.str(), "{ scale: 0.1, num_tested_scales: 4 }");
}
} // namespace
} // namespace operations_research

8
ortools/graph/graph.h
View File

@@ -272,14 +272,6 @@ class BaseGraph {
static const NodeIndexType kNilNode;
static const ArcIndexType kNilArc;
// TODO(user): remove the public functions below. They are just here during
// the transition from the old ebert_graph api to this new graph api.
template <typename A, typename B>
void GroupForwardArcsByFunctor(const A& a, B* b) {
LOG(FATAL) << "Not supported";
}
ArcIndexType max_end_arc_index() const { return arc_capacity_; }
protected:
// Functions commented when defined because they are implementation details.
void ComputeCumulativeSum(std::vector<ArcIndexType>* v);

2
ortools/graph/linear_assignment.h
View File

@@ -948,7 +948,7 @@ LinearSumAssignment<GraphType, CostValue>::LinearSumAssignment(
price_(num_left_nodes, 2 * num_left_nodes - 1),
matched_arc_(num_left_nodes, 0),
matched_node_(num_left_nodes, 2 * num_left_nodes - 1),
scaled_arc_cost_(graph.max_end_arc_index(), 0),
scaled_arc_cost_(graph.arc_capacity(), 0),
active_nodes_(absl::GetFlag(FLAGS_assignment_stack_order)
? static_cast<ActiveNodeContainerInterface*>(
new ActiveNodeStack())

3
ortools/graph/min_cost_flow.h
View File

@@ -174,6 +174,7 @@
#include <vector>
#include "absl/strings/string_view.h"
#include "ortools/base/logging.h"
#include "ortools/graph/graph.h"
#include "ortools/util/stats.h"
#include "ortools/util/zvector.h"
@@ -692,7 +693,7 @@ struct MinCostFlow : public MinCostFlowBase {
template <typename = void>
MinCostFlow() {
LOG(FATAL) << "MinCostFlow is deprecated. Use `SimpleMinCostFlow` or "
"`GenericMinCostFlow` with a specific graph type instead.");
"`GenericMinCostFlow` with a specific graph type instead.";
}
};

91
ortools/graph/python/fp_min_cost_flow.cc
View File

@@ -1,91 +0,0 @@
// Copyright 2010-2025 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ortools/graph/fp_min_cost_flow.h"
#include <sstream>
#include "ortools/graph/min_cost_flow.h"
#include "pybind11/numpy.h"
#include "pybind11/pybind11.h"
#include "pybind11/stl.h"
using ::operations_research::MinCostFlowBase;
using ::operations_research::SimpleFloatingPointMinCostFlow;
using ::pybind11::arg;
PYBIND11_MODULE(fp_min_cost_flow, m) {
pybind11::class_<SimpleFloatingPointMinCostFlow> smcf(
m, "SimpleFloatingPointMinCostFlow");
smcf.def(pybind11::init<>());
smcf.def("add_arc_with_capacity_and_unit_cost",
&SimpleFloatingPointMinCostFlow::AddArcWithCapacityAndUnitCost,
arg("tail"), arg("head"), arg("capacity"), arg("unit_cost"));
smcf.def("add_arcs_with_capacity_and_unit_cost",
pybind11::vectorize(
&SimpleFloatingPointMinCostFlow::AddArcWithCapacityAndUnitCost),
arg("tails"), arg("heads"), arg("capacities"), arg("unit_costs"));
smcf.def("set_node_supply", &SimpleFloatingPointMinCostFlow::SetNodeSupply,
arg("node"), arg("supply"));
smcf.def("set_nodes_supplies",
pybind11::vectorize(&SimpleFloatingPointMinCostFlow::SetNodeSupply),
arg("nodes"), arg("supplies"));
smcf.def("num_nodes", &SimpleFloatingPointMinCostFlow::NumNodes);
smcf.def("num_arcs", &SimpleFloatingPointMinCostFlow::NumArcs);
smcf.def("tail", &SimpleFloatingPointMinCostFlow::Tail, arg("arc"));
smcf.def("head", &SimpleFloatingPointMinCostFlow::Head, arg("arc"));
smcf.def("capacity", &SimpleFloatingPointMinCostFlow::Capacity, arg("arc"));
smcf.def("supply", &SimpleFloatingPointMinCostFlow::Supply, arg("node"));
smcf.def("unit_cost", &SimpleFloatingPointMinCostFlow::UnitCost, arg("arc"));
smcf.def("solve_max_flow_with_min_cost",
&SimpleFloatingPointMinCostFlow::SolveMaxFlowWithMinCost,
// Release the Python GIL during the solve. This function does not
// callback Python code so that is OK.
pybind11::call_guard<pybind11::gil_scoped_release>());
smcf.def("flow", &SimpleFloatingPointMinCostFlow::Flow, arg("arc"));
smcf.def("flows", pybind11::vectorize(&SimpleFloatingPointMinCostFlow::Flow),
arg("arcs"));
smcf.def("last_solve_stats", &SimpleFloatingPointMinCostFlow::LastSolveStats);
pybind11::class_<SimpleFloatingPointMinCostFlow::SolveStats> solve_stats(
smcf, "SolveStats");
solve_stats.def(pybind11::init<>());
solve_stats.def_readwrite("scale",
&SimpleFloatingPointMinCostFlow::SolveStats::scale);
solve_stats.def_readwrite(
"num_tested_scales",
&SimpleFloatingPointMinCostFlow::SolveStats::num_tested_scales);
solve_stats.def("__repr__",
[](const SimpleFloatingPointMinCostFlow::SolveStats& stats) {
std::ostringstream oss;
oss << "SolveStats" << stats;
return oss.str();
});
solve_stats.def("__str__",
[](const SimpleFloatingPointMinCostFlow::SolveStats& stats) {
std::ostringstream oss;
oss << stats;
return oss.str();
});
pybind11::enum_<SimpleFloatingPointMinCostFlow::Status>(smcf, "Status")
.value("BAD_COST_RANGE", MinCostFlowBase::Status::BAD_COST_RANGE)
.value("BAD_CAPACITY_RANGE", MinCostFlowBase::Status::BAD_CAPACITY_RANGE)
.value("BAD_RESULT", MinCostFlowBase::Status::BAD_RESULT)
.value("FEASIBLE", MinCostFlowBase::Status::FEASIBLE)
.value("INFEASIBLE", MinCostFlowBase::Status::INFEASIBLE)
.value("NOT_SOLVED", MinCostFlowBase::Status::NOT_SOLVED)
.value("OPTIMAL", MinCostFlowBase::Status::OPTIMAL)
.value("UNBALANCED", MinCostFlowBase::Status::UNBALANCED)
.export_values();
}

135
ortools/graph/python/fp_min_cost_flow_test.py
View File

@@ -1,135 +0,0 @@
#!/usr/bin/env python3
# Copyright 2010-2025 Google LLC
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""Test of the fp_min_cost_flow Python wrapper."""
import math
import numpy as np
from absl.testing import absltest
from ortools.graph.python import fp_min_cost_flow
class FpMinCostFlowTest(absltest.TestCase):
def test_success(self) -> None:
"""Test a successful solve.
We use the same test as scenario as the
SimpleFloatingPointMinCostFlowTest.Assignment test in
fp_min_cost_flow_test.cc.
"""
solver = fp_min_cost_flow.SimpleFloatingPointMinCostFlow()
solver.set_node_supply(node=0, supply=4.0)
solver.set_node_supply(node=1, supply=3.0)
solver.set_node_supply(node=2, supply=4.0)
solver.set_node_supply(node=3, supply=-5.0)
solver.set_node_supply(node=4, supply=-4.0)
a03 = solver.add_arc_with_capacity_and_unit_cost(
tail=0, head=3, capacity=5.1, unit_cost=2
)
a13 = solver.add_arc_with_capacity_and_unit_cost(
tail=1, head=3, capacity=4.0, unit_cost=1
)
a14 = solver.add_arc_with_capacity_and_unit_cost(
tail=1, head=4, capacity=5.0, unit_cost=2
)
a24 = solver.add_arc_with_capacity_and_unit_cost(
tail=2, head=4, capacity=3.0, unit_cost=3
)
status = solver.solve_max_flow_with_min_cost()
self.assertEqual(
status, fp_min_cost_flow.SimpleFloatingPointMinCostFlow.Status.OPTIMAL
)
self.assertEqual(solver.flow(a03), 3)
self.assertEqual(solver.flow(a13), 2)
self.assertEqual(solver.flow(a14), 1)
self.assertEqual(solver.flow(a24), 3)
def test_success_with_numpy(self) -> None:
"""Test a successful solve using numpy inputs/outputs.
Same test as test_success() but with numpy vector functions.
"""
solver = fp_min_cost_flow.SimpleFloatingPointMinCostFlow()
solver.set_nodes_supplies(
nodes=np.arange(5), supplies=np.array([4.0, 3.0, 4.0, -5.0, -4.0])
)
arcs = solver.add_arcs_with_capacity_and_unit_cost(
tails=np.array([0, 1, 1, 2]),
heads=np.array([3, 3, 4, 4]),
capacities=np.array([5.1, 4.0, 5.0, 3.0]),
unit_costs=np.array([2, 1, 2, 3]),
)
status = solver.solve_max_flow_with_min_cost()
self.assertEqual(
status, fp_min_cost_flow.SimpleFloatingPointMinCostFlow.Status.OPTIMAL
)
np.testing.assert_array_equal(
solver.flows(arcs), np.array([3.0, 2.0, 1.0, 3.0])
)
def test_failure(self) -> None:
"""Test a failing solve.
We use an infinite supply on a node, which is not supported.
"""
solver = fp_min_cost_flow.SimpleFloatingPointMinCostFlow()
solver.set_node_supply(node=0, supply=math.inf)
status = solver.solve_max_flow_with_min_cost()
self.assertEqual(
status,
fp_min_cost_flow.SimpleFloatingPointMinCostFlow.Status.BAD_CAPACITY_RANGE,
)
def test_solve_stats(self) -> None:
solve_stats = fp_min_cost_flow.SimpleFloatingPointMinCostFlow.SolveStats()
solve_stats.scale = 0.1
solve_stats.num_tested_scales = 3
self.assertEqual(
repr(solve_stats), "SolveStats{ scale: 0.1, num_tested_scales: 3 }"
)
self.assertEqual(str(solve_stats), "{ scale: 0.1, num_tested_scales: 3 }")
def test_last_solve_stats(self) -> None:
"""Test getting the last solve statistics."""
solver = fp_min_cost_flow.SimpleFloatingPointMinCostFlow()
solver.set_node_supply(node=0, supply=2.0)
status = solver.solve_max_flow_with_min_cost()
self.assertEqual(
status, fp_min_cost_flow.SimpleFloatingPointMinCostFlow.Status.OPTIMAL
)
solve_stats = solver.last_solve_stats()
self.assertIsInstance(
solve_stats, fp_min_cost_flow.SimpleFloatingPointMinCostFlow.SolveStats
)
# 2.0**61 is the last power of two such as:
#
# (2.0 * 2.0**61) < 2**63 - 1
#
# This is the largest scale we expect to use.
self.assertEqual(solve_stats.scale, 2.0**61)
# The code should have directly tested the correct largest scale.
self.assertEqual(solve_stats.num_tested_scales, 1)
if __name__ == "__main__":
absltest.main()