diff --git a/ortools/graph/BUILD.bazel b/ortools/graph/BUILD.bazel index 0d55839620..35525277bf 100644 --- a/ortools/graph/BUILD.bazel +++ b/ortools/graph/BUILD.bazel @@ -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", diff --git a/ortools/graph/CMakeLists.txt b/ortools/graph/CMakeLists.txt index 153fbaf09e..0a54e293e7 100644 --- a/ortools/graph/CMakeLists.txt +++ b/ortools/graph/CMakeLists.txt @@ -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) diff --git a/ortools/graph/fp_min_cost_flow.cc b/ortools/graph/fp_min_cost_flow.cc deleted file mode 100644 index 19579ba9f1..0000000000 --- a/ortools/graph/fp_min_cost_flow.cc +++ /dev/null @@ -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 -#include -#include -#include -#include -#include -#include -#include - -#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::max(); -constexpr SimpleFloatingPointMinCostFlow::FPFlowQuantity kMaxFPFlow = - std::numeric_limits::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 max_node_in_flow(num_nodes); - std::vector 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, - "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::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 max_node_in_flow(num_nodes, 0.0); - std::vector 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 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( - // 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::max(), - static_cast(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 diff --git a/ortools/graph/fp_min_cost_flow.h b/ortools/graph/fp_min_cost_flow.h deleted file mode 100644 index 24e61fc4a8..0000000000 --- a/ortools/graph/fp_min_cost_flow.h +++ /dev/null @@ -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 -#include -#include -#include -#include - -#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 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( - // 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 arc_capacity_; - std::vector node_supply_; - std::vector 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::max(); - if (rounded_scaled_flow >= kMaxFlowQuantity) { - return kMaxFlowQuantity; - } - if (rounded_scaled_flow <= -kMaxFlowQuantity) { - return -kMaxFlowQuantity; - } - return static_cast(rounded_scaled_flow); -} - -} // namespace internal -} // namespace operations_research - -#endif // OR_TOOLS_GRAPH_FP_MIN_COST_FLOW_H_ diff --git a/ortools/graph/fp_min_cost_flow_test.cc b/ortools/graph/fp_min_cost_flow_test.cc deleted file mode 100644 index 76ae33f4b0..0000000000 --- a/ortools/graph/fp_min_cost_flow_test.cc +++ /dev/null @@ -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 -#include -#include -#include -#include -#include - -#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::infinity(); -constexpr double kNaN = std::numeric_limits::quiet_NaN(); -constexpr SimpleFloatingPointMinCostFlow::CostValue kMaxCostValue = - std::numeric_limits::max(); -const SimpleFloatingPointMinCostFlow::FPFlowQuantity kMaxFPFlow = - std::numeric_limits::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 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 {}; - -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 {}; - -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 {}; - -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& 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 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::max(); - EXPECT_EQ(ScaleFlow(/*fp_flow=*/kMaxFPFlow, - /*scale=*/std::numeric_limits::max()), - kMaxFlowQuantity); - EXPECT_EQ(ScaleFlow(/*fp_flow=*/-kMaxFPFlow, - /*scale=*/std::numeric_limits::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(kMaxFlowQuantity), 0.0); - const SimpleMinCostFlow::FlowQuantity kMaxFlowQuantityPredecessorAsInt = - static_cast(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 diff --git a/ortools/graph/graph.h b/ortools/graph/graph.h index e9840249e1..7f3a43fe77 100644 --- a/ortools/graph/graph.h +++ b/ortools/graph/graph.h @@ -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 - 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* v); diff --git a/ortools/graph/linear_assignment.h b/ortools/graph/linear_assignment.h index 110d48c4ef..562f0a8069 100644 --- a/ortools/graph/linear_assignment.h +++ b/ortools/graph/linear_assignment.h @@ -948,7 +948,7 @@ LinearSumAssignment::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( new ActiveNodeStack()) diff --git a/ortools/graph/min_cost_flow.h b/ortools/graph/min_cost_flow.h index edce7a179c..90964c7893 100644 --- a/ortools/graph/min_cost_flow.h +++ b/ortools/graph/min_cost_flow.h @@ -174,6 +174,7 @@ #include #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 MinCostFlow() { LOG(FATAL) << "MinCostFlow is deprecated. Use `SimpleMinCostFlow` or " - "`GenericMinCostFlow` with a specific graph type instead."); + "`GenericMinCostFlow` with a specific graph type instead."; } }; diff --git a/ortools/graph/python/fp_min_cost_flow.cc b/ortools/graph/python/fp_min_cost_flow.cc deleted file mode 100644 index abe43168fd..0000000000 --- a/ortools/graph/python/fp_min_cost_flow.cc +++ /dev/null @@ -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 - -#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_ 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()); - 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_ 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_(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(); -} diff --git a/ortools/graph/python/fp_min_cost_flow_test.py b/ortools/graph/python/fp_min_cost_flow_test.py deleted file mode 100644 index 88b39826b1..0000000000 --- a/ortools/graph/python/fp_min_cost_flow_test.py +++ /dev/null @@ -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()