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

qap_sat solver and qap reader

Browse Source
This commit is contained in:
Laurent Perron
2022-10-10 15:03:10 +02:00
parent 7efb7892f7
commit a25d93fcfa
6 changed files with 447 additions and 100 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
examples/cpp/BUILD.bazel
View File

@@ -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"],

155
examples/cpp/qap_sat.cc Normal file
View File

@@ -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 <limits>
#include <string>
#include <vector>
#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<std::vector<BoolVar>> CreatePlaceVars(int n,
CpModelBuilder& cp_model) {
const BoolVar false_literal = cp_model.FalseVar();
std::vector<std::vector<BoolVar>> 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<std::vector<BoolVar>>& 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<BoolVar> 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<BoolVar> 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<std::vector<BoolVar>>& 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<std::vector<BoolVar>> 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), &parameters))
<< 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;
}

44
ortools/util/BUILD.bazel
View File

@@ -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"],

72
ortools/util/qap_reader.cc Normal file
View File

@@ -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 <limits>
#include <string>
#include <vector>
#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<std::string> 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<int64_t>::max());
qap_problem.distances[j].assign(n, std::numeric_limits<int64_t>::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

66
ortools/util/qap_reader.h Normal file
View File

@@ -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 <string>
#include <vector>
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<std::vector<int64_t>> weights;
// `distances[i][j]` is the distance from location i to location j.
std::vector<std::vector<int64_t>> 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_

196
ortools/util/saturated_arithmetic.h
View File

@@ -14,13 +14,53 @@
#ifndef OR_TOOLS_UTIL_SATURATED_ARITHMETIC_H_
#define OR_TOOLS_UTIL_SATURATED_ARITHMETIC_H_
#include <cstdint>
#include <limits>
#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<int64_t>::min() ||
x == std::numeric_limits<int64_t>::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<int64_t>::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<int64_t>(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<int64_t>::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<int64_t>::min() ||
x == -std::numeric_limits<int64_t>::max();
}
} // namespace operations_research
#endif // OR_TOOLS_UTIL_SATURATED_ARITHMETIC_H_