diff --git a/examples/cpp/BUILD.bazel b/examples/cpp/BUILD.bazel index 96783238d4..a9478fd08e 100644 --- a/examples/cpp/BUILD.bazel +++ b/examples/cpp/BUILD.bazel @@ -434,6 +434,20 @@ cc_binary( ], ) + +cc_binary( + name = "qap_sat", + srcs = ["qap_sat.cc"], + deps = [ + "//ortools/base", + "//ortools/sat:cp_model", + "//ortools/sat:model", + "//ortools/sat:sat_parameters_cc_proto", + "//ortools/util:qap_reader", + "@com_google_absl//absl/strings", + ], +) + cc_binary( name = "random_tsp", srcs = ["random_tsp.cc"], diff --git a/examples/cpp/qap_sat.cc b/examples/cpp/qap_sat.cc new file mode 100644 index 0000000000..a70541d13c --- /dev/null +++ b/examples/cpp/qap_sat.cc @@ -0,0 +1,155 @@ +// Copyright 2010-2022 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 +#include +#include + +#include "absl/flags/flag.h" +#include "absl/strings/str_cat.h" +#include "google/protobuf/text_format.h" +#include "ortools/base/commandlineflags.h" +#include "ortools/base/init_google.h" +#include "ortools/base/logging.h" +#include "ortools/sat/cp_model.h" +#include "ortools/sat/sat_parameters.pb.h" +#include "ortools/util/qap_reader.h" + +ABSL_FLAG(std::string, input, "", "Input file name containing a QAP instance."); +ABSL_FLAG(std::string, params, "", "Specific params to use with CP-SAT."); + +namespace operations_research { +namespace sat { +namespace { + +// place_vars[f][l] contains the binary variable that decides whether to put +// factory f in location l. +std::vector> CreatePlaceVars(int n, + CpModelBuilder& cp_model) { + const BoolVar false_literal = cp_model.FalseVar(); + std::vector> place_vars; + place_vars.resize(n); + for (int f = 0; f < n; ++f) { + place_vars[f].assign(n, false_literal); + } + + for (int f = 0; f < n; ++f) { + for (int l = 0; l < n; ++l) { + place_vars[f][l] = + cp_model.NewBoolVar().WithName(absl::StrCat("place_", f, "_", l)); + } + } + return place_vars; +} + +void CreatePlacementConstraints( + const std::vector>& place_vars, + CpModelBuilder& cp_model) { + const int n = place_vars.size(); + + // Place each factory exactly once. + for (int f = 0; f < n; ++f) { + std::vector tmp; + for (int l = 0; l < n; ++l) { + tmp.push_back(place_vars[f][l]); + } + cp_model.AddExactlyOne(tmp); + } + + // Occupy each location exactly once. + for (int l = 0; l < n; ++l) { + std::vector tmp; + for (int f = 0; f < n; ++f) { + tmp.push_back(place_vars[f][l]); + } + cp_model.AddExactlyOne(tmp); + } +} + +void CreateObjective(const QapProblem& qap, + const std::vector>& place_vars, + CpModelBuilder& cp_model) { + const int n = place_vars.size(); + + LinearExpr objective; + for (int f1 = 0; f1 < n; ++f1) { + for (int f2 = 0; f2 < n; ++f2) { + if (f1 == f2) continue; + if (qap.weights[f1][f2] == 0) continue; + for (int l1 = 0; l1 < n; ++l1) { + for (int l2 = 0; l2 < n; ++l2) { + if (l1 == l2) continue; + if (qap.distances[l1][l2] == 0) continue; + + BoolVar product = cp_model.NewBoolVar(); + cp_model.AddMultiplicationEquality(product, place_vars[f1][l1], + place_vars[f2][l2]); + + objective += product * qap.weights[f1][f2] * qap.distances[l1][l2]; + } + } + } + } + + for (int f = 0; f < n; ++f) { + for (int l = 0; l < n; ++l) { + objective += place_vars[f][l] * qap.weights[f][f] * qap.distances[l][l]; + } + } + + cp_model.Minimize(objective); +} + +} // namespace + +void SolveQap() { + LOG(INFO) << "Reading QAP problem from '" << absl::GetFlag(FLAGS_input) + << "'."; + QapProblem qap = ReadQapProblemOrDie(absl::GetFlag(FLAGS_input)); + const int n = qap.weights.size(); + + CpModelBuilder cp_model; + std::vector> place_vars = CreatePlaceVars(n, cp_model); + CreatePlacementConstraints(place_vars, cp_model); + CreateObjective(qap, place_vars, cp_model); + + // Setup parameters. + SatParameters parameters; + parameters.set_log_search_progress(true); + // Parse the --params flag. + if (!absl::GetFlag(FLAGS_params).empty()) { + CHECK(google::protobuf::TextFormat::MergeFromString( + absl::GetFlag(FLAGS_params), ¶meters)) + << absl::GetFlag(FLAGS_params); + } + + // Solve. + const CpSolverResponse response = + SolveWithParameters(cp_model.Build(), parameters); +} + +} // namespace sat +} // namespace operations_research + +static const char kUsageStr[] = + "Solves quadratic assigment problems with CP-SAT. " + "The input file should have the format described in the QAPLIB (see " + "http://anjos.mgi.polymtl.ca/qaplib/)."; + +int main(int argc, char** argv) { + InitGoogle(kUsageStr, &argc, &argv, /*remove_flags=*/true); + CHECK(!absl::GetFlag(FLAGS_input).empty()) << "--input is required"; + + operations_research::sat::SolveQap(); + return EXIT_SUCCESS; +} diff --git a/ortools/util/BUILD.bazel b/ortools/util/BUILD.bazel index ddc26dd82a..9900da37ea 100644 --- a/ortools/util/BUILD.bazel +++ b/ortools/util/BUILD.bazel @@ -412,40 +412,6 @@ cc_library( ], ) -#cc_library( -# name = "tsplib_parser", -# srcs = ["tsplib_parser.cc"], -# hdrs = ["tsplib_parser.h"], -# visibility = ["//visibility:public"], -# deps = [ -# "@com_google_absl//absl/strings", -# ":filelineiter", -# "//ortools/base", -# "//ortools/base:file", -# "//ortools/base:file:path", -# "//ortools/base:map_util", -# "//ortools/base:mathutil", -# "/@com_google_absl//absl/strings", -# "//file/zipfile", -# "@com_google_re2//:re2", -# ], -#) - -#cc_library( -# name = "pdtsp_parser", -# srcs = ["pdtsp_parser.cc"], -# hdrs = ["pdtsp_parser.h"], -# visibility = ["//visibility:public"], -# deps = [ -# "@com_google_absl//absl/strings", -# ":filelineiter", -# "//ortools/base", -# "//ortools/base:file", -# "//ortools/base:file:path", -# "//ortools/base:mathutil", -# ], -#) - #cc_library( # name = "bp_parser", # srcs = ["bp_parser.cc"], @@ -459,6 +425,16 @@ cc_library( # ], #) +cc_library( + name = "qap_reader", + srcs = ["qap_reader.cc"], + hdrs = ["qap_reader.h"], + deps = [ + "//ortools/util:filelineiter", + "@com_google_absl//absl/strings", + ], +) + cc_library( name = "sort", hdrs = ["sort.h"], diff --git a/ortools/util/qap_reader.cc b/ortools/util/qap_reader.cc new file mode 100644 index 0000000000..3a6bf2d03f --- /dev/null +++ b/ortools/util/qap_reader.cc @@ -0,0 +1,72 @@ +// Copyright 2010-2022 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/util/qap_reader.h" + +#include +#include +#include + +#include "absl/strings/numbers.h" +#include "absl/strings/str_split.h" +#include "ortools/util/filelineiter.h" + +namespace operations_research { + +// TODO(user): Unit test cases when the function dies, or return +// (and test) a status instead. +QapProblem ReadQapProblemOrDie(const std::string& filepath) { + QapProblem qap_problem; + + int n = 0; + int k = 0; + for (const std::string& line : FileLines(filepath)) { + const std::vector tokens = + absl::StrSplit(line, ' ', absl::SkipEmpty()); + if (tokens.empty()) continue; + if (k == 0) { + CHECK_EQ(tokens.size(), 1); + CHECK(absl::SimpleAtoi(tokens[0], &n)); + qap_problem.weights.resize(n); + qap_problem.distances.resize(n); + for (int j = 0; j < n; ++j) { + qap_problem.weights[j].assign(n, std::numeric_limits::max()); + qap_problem.distances[j].assign(n, std::numeric_limits::max()); + } + ++k; + } else if (k <= n * n) { + for (const std::string& token : tokens) { + const int i = (k - 1) / n; + const int j = (k - 1) % n; + int64_t v = 0.0; + CHECK(absl::SimpleAtoi(token, &v)); + qap_problem.weights[i][j] = v; + ++k; + } + } else if (k <= 2 * n * n) { + for (const std::string& token : tokens) { + const int i = (k - n * n - 1) / n; + const int j = (k - n * n - 1) % n; + int64_t v = 0.0; + CHECK(absl::SimpleAtoi(token, &v)); + qap_problem.distances[i][j] = v; + ++k; + } + } else { + CHECK(false) << "File contains more than 1 + 2 * N^2 entries."; + } + } + return qap_problem; +} + +} // namespace operations_research diff --git a/ortools/util/qap_reader.h b/ortools/util/qap_reader.h new file mode 100644 index 0000000000..ea6449d762 --- /dev/null +++ b/ortools/util/qap_reader.h @@ -0,0 +1,66 @@ +// Copyright 2010-2022 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_UTIL_QAP_READER_H_ +#define OR_TOOLS_UTIL_QAP_READER_H_ + +#include +#include + +namespace operations_research { + +// Quadratic assignment problem (QAP) is a classical combinatorial optimization +// problem. See https://en.wikipedia.org/wiki/Quadratic_assignment_problem. +// In short, there are a set of n facilities and a set of n locations. +// For each pair of locations, a `distance` is specified and for each pair of +// facilities a `weight` is specified (e.g., the amount of supplies transported +// between the two facilities). The problem is to assign all facilities to +// different locations with the goal of minimizing the sum of the distances +// multiplied by the corresponding flows. +struct QapProblem { + // `weights[i][j]` is the amount of flow from facility i to facility j. + std::vector> weights; + + // `distances[i][j]` is the distance from location i to location j. + std::vector> distances; + + // For testing. + bool operator==(const QapProblem& q) const { + return weights == q.weights && distances == q.distances; + } +}; + +// Reads a QAP problem from file in a format used in QAPLIB. See +// anjos.mgi.polymtl.ca/qaplib/ for more context. The format is "n W D", where +// n is the number of factories/locations, and W and D are weight and +// distance matrices, respectively. Both W and D are square matrices of size +// N x N. Each entry of the matrices must parse as double (we CHECK if it +// does). Multiple spaces, or '\n' are disregarded. For example: +// +// 2 +// +// 0 2 +// 3 0 +// +// 0 2 +// 1 0 +// +// defines a problem with two factories and two locations. There are 2 units of +// supply flowing from factory 0 to factory 1, and 3 units of supply flowing +// from factory 1 to 0. The distance from location 0 to location 1 is equal +// to 2, and the distance from location 1 to 0 is equal to 1. +QapProblem ReadQapProblemOrDie(const std::string& filepath); + +} // namespace operations_research + +#endif // OR_TOOLS_UTIL_QAP_READER_H_ diff --git a/ortools/util/saturated_arithmetic.h b/ortools/util/saturated_arithmetic.h index c2ec1c1498..986d902ded 100644 --- a/ortools/util/saturated_arithmetic.h +++ b/ortools/util/saturated_arithmetic.h @@ -14,13 +14,53 @@ #ifndef OR_TOOLS_UTIL_SATURATED_ARITHMETIC_H_ #define OR_TOOLS_UTIL_SATURATED_ARITHMETIC_H_ +#include #include #include "absl/base/casts.h" #include "ortools/base/integral_types.h" #include "ortools/util/bitset.h" +// This file contains implementations for saturated addition, subtraction and +// multiplication. +// Currently, there are three versions of the code. +// The code using the built-ins provided since GCC 5.0 now compiles really well +// with clang/LLVM on both x86_64 and ARM. It should therefore be the standard, +// but not for multiplication, see below. +// +// For example, on ARM, only 4 instructions are needed, two of which are +// additions that can be executed in parallel. +// +// On x86_64, we're keeping the code with inline assembly for GCC as GCC does +// manage to compile the code with built-ins properly. +// On x86_64, the product of two 64-bit registers is a 128-bit integer +// stored in two 64-bit registers. It's the carry flag that is set when the +// result exceeds 64 bits, not the overflow flag. Since the built-in uses the +// overflow flag, we have to resort on the assembly-based version of the code. +// +// Sadly, MSVC does not support the built-ins nor does it support inline +// assembly. We have to rely on the generic, C++ only version of the code which +// is much slower. +// +// TODO(user): make this implementation the default everywhere. +// TODO(user): investigate the code generated by MSVC. + namespace operations_research { + +// Checks if x is equal to the min or the max value of an int64_t. +inline bool AtMinOrMaxInt64(int64_t x) { + return x == std::numeric_limits::min() || + x == std::numeric_limits::max(); +} + +// Note(user): -kint64min != kint64max, but kint64max == ~kint64min. +inline int64_t CapOpp(int64_t v) { return v == kint64min ? ~v : -v; } + +inline int64_t CapAbs(int64_t v) { + return v == kint64min ? std::numeric_limits::max() + : (v < 0 ? -v : v); +} + // ---------- Overflow utility functions ---------- // Implement two's complement addition and subtraction on int64s. @@ -103,14 +143,12 @@ inline int64_t CapWithSignOf(int64_t x) { return TwosComplementAddition(kint64max, static_cast(x < 0)); } -inline int64_t CapAddGeneric(int64_t x, int64_t y) { - const int64_t result = TwosComplementAddition(x, y); - return AddHadOverflow(x, y, result) ? CapWithSignOf(x) : result; -} - -#if defined(__GNUC__) && defined(__x86_64__) -// TODO(user): port this to other architectures. -inline int64_t CapAddFast(int64_t x, int64_t y) { +// The following implementations are here for GCC because it does not +// compiled the built-ins correctly. The code is either too long without +// branches or contains jumps. These implementations are probably optimal +// on x86_64. +#if defined(__GNUC__) && !defined(__clang_) && defined(__x86_64__) +inline int64_t CapAddAsm(int64_t x, int64_t y) { const int64_t cap = CapWithSignOf(x); int64_t result = x; // clang-format off @@ -123,26 +161,8 @@ inline int64_t CapAddFast(int64_t x, int64_t y) { // clang-format on return result; } -#endif -inline int64_t CapAdd(int64_t x, int64_t y) { -#if defined(__GNUC__) && defined(__x86_64__) - return CapAddFast(x, y); -#else - return CapAddGeneric(x, y); -#endif -} - -inline void CapAddTo(int64_t x, int64_t* y) { *y = CapAdd(*y, x); } - -inline int64_t CapSubGeneric(int64_t x, int64_t y) { - const int64_t result = TwosComplementSubtraction(x, y); - return SubHadOverflow(x, y, result) ? CapWithSignOf(x) : result; -} - -#if defined(__GNUC__) && defined(__x86_64__) -// TODO(user): port this to other architectures. -inline int64_t CapSubFast(int64_t x, int64_t y) { +inline int64_t CapSubAsm(int64_t x, int64_t y) { const int64_t cap = CapWithSignOf(x); int64_t result = x; // clang-format off @@ -155,22 +175,67 @@ inline int64_t CapSubFast(int64_t x, int64_t y) { // clang-format on return result; } + +// Note that on x86_64, we have to use this code because it's the carry flag +// that is set when the product of two 64-bit integers does not fit in 64-bit. +inline int64_t CapProdAsm(int64_t x, int64_t y) { + // cap = kint64max if x and y have the same sign, cap = kint64min + // otherwise. + const int64_t cap = CapWithSignOf(x ^ y); + int64_t result = x; + // Here, we use the fact that imul of two signed 64-integers returns a 128-bit + // result -- we care about the lower 64 bits. More importantly, imul also sets + // the carry flag if 64 bits were not enough. + // We therefore use cmovc to return cap if the carry was set. + // clang-format off + asm volatile( // 'volatile': ask compiler optimizer "keep as is". + "\n\t" "imulq %[y],%[result]" + "\n\t" "cmovcq %[cap],%[result]" // Conditional move if carry. + : [result] "=r"(result) // Output + : "[result]" (result), [y] "r"(y), [cap] "r"(cap) // Input. + : "cc" /* Clobbered registers */); + // clang-format on + return result; +} #endif -inline int64_t CapSub(int64_t x, int64_t y) { -#if defined(__GNUC__) && defined(__x86_64__) - return CapSubFast(x, y); -#else - return CapSubGeneric(x, y); -#endif +// Simple implementations which use the built-ins provided by both GCC and +// clang. clang compiles to very good code for both x86_64 and ARM. This is the +// preferred implementation in general. +#if defined(__clang__) +inline int64_t CapAddBuiltIn(int64_t x, int64_t y) { + const int64_t cap = CapWithSignOf(x); + int64_t result; + const bool overflowed = __builtin_add_overflow(x, y, &result); + return overflowed ? cap : result; } -// Note(user): -kint64min != kint64max, but kint64max == ~kint64min. -inline int64_t CapOpp(int64_t v) { return v == kint64min ? ~v : -v; } +inline int64_t CapSubBuiltIn(int64_t x, int64_t y) { + const int64_t cap = CapWithSignOf(x); + int64_t result; + const bool overflowed = __builtin_sub_overflow(x, y, &result); + return overflowed ? cap : result; +} -inline int64_t CapAbs(int64_t v) { - return v == kint64min ? std::numeric_limits::max() - : (v < 0 ? -v : v); +// As said above, this is useless on x86_64. +inline int64_t CapProdBuiltIn(int64_t x, int64_t y) { + const int64_t cap = CapWithSignOf(x ^ y); + int64_t result; + const bool overflowed = __builtin_mul_overflow(x, y, &result); + return overflowed ? cap : result; +} +#endif + +// Generic implementations. They are very good for addition and subtraction, +// less so for multiplication. +inline int64_t CapAddGeneric(int64_t x, int64_t y) { + const int64_t result = TwosComplementAddition(x, y); + return AddHadOverflow(x, y, result) ? CapWithSignOf(x) : result; +} + +inline int64_t CapSubGeneric(int64_t x, int64_t y) { + const int64_t result = TwosComplementSubtraction(x, y); + return SubHadOverflow(x, y, result) ? CapWithSignOf(x) : result; } namespace cap_prod_util { @@ -216,43 +281,42 @@ inline int64_t CapProdGeneric(int64_t x, int64_t y) { return cap < 0 ? -abs_result : abs_result; } -#if defined(__GNUC__) && defined(__x86_64__) -// TODO(user): port this to other architectures. -inline int64_t CapProdFast(int64_t x, int64_t y) { - // cap = kint64max if x and y have the same sign, cap = kint64min - // otherwise. - const int64_t cap = CapWithSignOf(x ^ y); - int64_t result = x; - // Here, we use the fact that imul of two signed 64-integers returns a 128-bit - // result -- we care about the lower 64 bits. More importantly, imul also sets - // the carry flag if 64 bits were not enough. - // We therefore use cmovc to return cap if the carry was set. - // clang-format off - asm volatile( // 'volatile': ask compiler optimizer "keep as is". - "\n\t" "imulq %[y],%[result]" - "\n\t" "cmovcq %[cap],%[result]" // Conditional move if carry. - : [result] "=r"(result) // Output - : "[result]" (result), [y] "r"(y), [cap] "r"(cap) // Input. - : "cc" /* Clobbered registers */); - // clang-format on - return result; -} +inline int64_t CapAdd(int64_t x, int64_t y) { +#if defined(__GNUC__) && !defined(__clang__) && defined(__x86_64__) + return CapAddAsm(x, y); +#elif defined(__clang__) + return CapAddBuiltIn(x, y); +#else + return CapAddGeneric(x, y); #endif +} + +inline void CapAddTo(int64_t x, int64_t* y) { *y = CapAdd(*y, x); } + +inline int64_t CapSub(int64_t x, int64_t y) { +#if defined(__GNUC__) && !defined(__clang__) && defined(__x86_64__) + return CapSubAsm(x, y); +#elif defined(__clang__) + return CapSubBuiltIn(x, y); +#else + return CapSubGeneric(x, y); +#endif +} inline int64_t CapProd(int64_t x, int64_t y) { #if defined(__GNUC__) && defined(__x86_64__) - return CapProdFast(x, y); + // On x86_64, the product of two 64-bit registeres is a 128-bit integer, + // stored in two 64-bit registers. It's the carry flag that is set when the + // result exceeds 64 bits, not the overflow flag. We therefore have to resort + // to the assembly-based version of the code. + return CapProdAsm(x, y); +#elif defined(__clang__) + return CapProdBuiltIn(x, y); #else return CapProdGeneric(x, y); #endif } -// Checks is x is equal to the min or the max value of an int64_t. -inline bool AtMinOrMaxInt64(int64_t x) { - return x == std::numeric_limits::min() || - x == -std::numeric_limits::max(); -} - } // namespace operations_research #endif // OR_TOOLS_UTIL_SATURATED_ARITHMETIC_H_