Page MenuHome
Differential D5594 Diff 18116 extern/quadriflow/src/localsat.cpp

Changeset View

Standalone View

View Options

extern/quadriflow/src/localsat.cpp

  • This file was added.
#ifdef NDEBUG
#undef NDEBUG
#endif
#include "localsat.hpp"
#include "config.hpp"
#include "dedge.hpp"
#include "field-math.hpp"
#include <Eigen/Core>
#include <deque>
#include <memory>
#include <utility>
#include <vector>
namespace qflow {
const int max_depth = 0;
using namespace Eigen;
SolverStatus RunCNF(const std::string &fin_name, int n_variable, int timeout,
const std::vector<std::vector<int>> &sat_clause, std::vector<int> &value) {
int n_sat_variable = 3 * n_variable;
auto fout_name = fin_name + ".result.txt";
FILE *fout = fopen(fin_name.c_str(), "w");
fprintf(fout, "p cnf %d %d\n", n_sat_variable, (int)sat_clause.size());
for (auto &c : sat_clause) {
for (auto e : c) fprintf(fout, "%d ", e);
fputs("0\n", fout);
}
fclose(fout);
char cmd[100];
snprintf(cmd, 99, "rm %s > /dev/null 2>&1", fout_name.c_str());
system(cmd);
snprintf(cmd, 99, "timeout %d minisat %s %s > /dev/null 2>&1", timeout, fin_name.c_str(),
fout_name.c_str());
int exit_code = system(cmd);
FILE *fin = fopen(fout_name.c_str(), "r");
char buf[16] = {0};
fscanf(fin, "%15s", buf);
lprintf(" MiniSAT:");
if (strcmp(buf, "SAT") != 0) {
fclose(fin);
if (exit_code == 124) {
lprintf(" Timeout! ");
return SolverStatus::Timeout;
}
lprintf(" Unsatisfiable! ");
return SolverStatus::Unsat;
};
lprintf(" Satisfiable! ");
for (int i = 0; i < n_variable; ++i) {
int sign[3];
fscanf(fin, "%d %d %d", sign + 0, sign + 1, sign + 2);
int nvalue = -2;
for (int j = 0; j < 3; ++j) {
assert(abs(sign[j]) == 3 * i + j + 1);
if ((sign[j] > 0) == (value[i] != j - 1)) {
assert(nvalue == -2);
nvalue = j - 1;
}
}
value[i] = nvalue;
}
fclose(fin);
return SolverStatus::Sat;
}
SolverStatus SolveSatProblem(int n_variable, std::vector<int> &value,
const std::vector<bool> flexible, // NOQA
const std::vector<Vector3i> &variable_eq,
const std::vector<Vector3i> &constant_eq,
const std::vector<Vector4i> &variable_ge,
const std::vector<Vector2i> &constant_ge,
int timeout) {
for (int v : value) assert(-1 <= v && v <= +1);
auto VAR = [&](int i, int v) {
int index = 1 + 3 * i + v + 1;
// We initialize the SAT problem by setting all the variable to false.
// This is because minisat by default will try false first.
if (v == value[i]) index = -index;
return index;
};
int n_flexible = 0;
std::vector<std::vector<int>> sat_clause;
std::vector<bool> sat_ishard;
auto add_clause = [&](const std::vector<int> &clause, bool hard) {
sat_clause.push_back(clause);
sat_ishard.push_back(hard);
};
for (int i = 0; i < n_variable; ++i) {
add_clause({-VAR(i, -1), -VAR(i, 0)}, true);
add_clause({-VAR(i, +1), -VAR(i, 0)}, true);
add_clause({-VAR(i, -1), -VAR(i, +1)}, true);
add_clause({VAR(i, -1), VAR(i, 0), VAR(i, +1)}, true);
if (!flexible[i]) {
add_clause({VAR(i, value[i])}, true);
} else {
++n_flexible;
}
}
for (int i = 0; i < (int)variable_eq.size(); ++i) {
auto &var = variable_eq[i];
auto &cst = constant_eq[i];
for (int v0 = -1; v0 <= 1; ++v0)
for (int v1 = -1; v1 <= 1; ++v1)
for (int v2 = -1; v2 <= 1; ++v2)
if (cst[0] * v0 + cst[1] * v1 + cst[2] * v2 != 0) {
add_clause({-VAR(var[0], v0), -VAR(var[1], v1), -VAR(var[2], v2)}, true);
}
}
for (int i = 0; i < (int)variable_ge.size(); ++i) {
auto &var = variable_ge[i];
auto &cst = constant_ge[i];
for (int v0 = -1; v0 <= 1; ++v0)
for (int v1 = -1; v1 <= 1; ++v1)
for (int v2 = -1; v2 <= 1; ++v2)
for (int v3 = -1; v3 <= 1; ++v3)
if (cst[0] * v0 * v1 - cst[1] * v2 * v3 < 0) {
add_clause({-VAR(var[0], v0), -VAR(var[1], v1), -VAR(var[2], v2),
-VAR(var[3], v3)},
false);
}
}
int nflip_before = 0, nflip_after = 0;
for (int i = 0; i < (int)variable_ge.size(); ++i) {
auto &var = variable_ge[i];
auto &cst = constant_ge[i];
if (value[var[0]] * value[var[1]] * cst[0] - value[var[2]] * value[var[3]] * cst[1] < 0)
nflip_before++;
}
lprintf(" [SAT] nvar: %6d nflip: %3d ", n_flexible * 2, nflip_before);
auto rcnf = RunCNF("test.out", n_variable, timeout, sat_clause, value);
for (int i = 0; i < (int)variable_eq.size(); ++i) {
auto &var = variable_eq[i];
auto &cst = constant_eq[i];
assert(cst[0] * value[var[0]] + cst[1] * value[var[1]] + cst[2] * value[var[2]] == 0);
}
for (int i = 0; i < (int)variable_ge.size(); ++i) {
auto &var = variable_ge[i];
auto &cst = constant_ge[i];
int area = value[var[0]] * value[var[1]] * cst[0] - value[var[2]] * value[var[3]] * cst[1];
if (area < 0) ++nflip_after;
}
lprintf("nflip: %3d\n", nflip_after);
return rcnf;
}
void ExportLocalSat(std::vector<Vector2i> &edge_diff, const std::vector<Vector3i> &face_edgeIds,
const std::vector<Vector3i> &face_edgeOrients, const MatrixXi &F,
const VectorXi &V2E, const VectorXi &E2E) {
int flip_count = 0;
int flip_count1 = 0;
std::vector<int> value(2 * edge_diff.size());
for (int i = 0; i < (int)edge_diff.size(); ++i) {
value[2 * i + 0] = edge_diff[i][0];
value[2 * i + 1] = edge_diff[i][1];
}
std::deque<std::pair<int, int>> Q;
std::vector<bool> mark_vertex(V2E.size(), false);
assert(F.cols() == (int)face_edgeIds.size());
std::vector<Vector3i> variable_eq(face_edgeIds.size() * 2);
std::vector<Vector3i> constant_eq(face_edgeIds.size() * 2);
std::vector<Vector4i> variable_ge(face_edgeIds.size());
std::vector<Vector2i> constant_ge(face_edgeIds.size());
VectorXd face_area(F.cols());
for (int i = 0; i < (int)face_edgeIds.size(); ++i) {
Vector2i diff[3];
Vector2i var[3];
Vector2i cst[3];
for (int j = 0; j < 3; ++j) {
int edgeid = face_edgeIds[i][j];
diff[j] = rshift90(edge_diff[edgeid], face_edgeOrients[i][j]);
var[j] = rshift90(Vector2i(edgeid * 2 + 1, edgeid * 2 + 2), face_edgeOrients[i][j]);
cst[j] = var[j].array().sign();
var[j] = var[j].array().abs() - 1;
}
assert(diff[0] + diff[1] + diff[2] == Vector2i::Zero());
variable_eq[2 * i + 0] = Vector3i(var[0][0], var[1][0], var[2][0]);
constant_eq[2 * i + 0] = Vector3i(cst[0][0], cst[1][0], cst[2][0]);
variable_eq[2 * i + 1] = Vector3i(var[0][1], var[1][1], var[2][1]);
constant_eq[2 * i + 1] = Vector3i(cst[0][1], cst[1][1], cst[2][1]);
face_area[i] = diff[0][0] * diff[1][1] - diff[0][1] * diff[1][0];
if (face_area[i] < 0) {
printf("[SAT] Face %d's area < 0\n", i);
for (int j = 0; j < 3; ++j) {
int v = F(j, i);
if (mark_vertex[v]) continue;
Q.push_back(std::make_pair(v, 0));
mark_vertex[v] = true;
}
flip_count += 1;
}
variable_ge[i] = Vector4i(var[0][0], var[1][1], var[0][1], var[1][0]);
constant_ge[i] = Vector2i(cst[0][0] * cst[1][1], cst[0][1] * cst[1][0]);
}
for (int i = 0; i < (int)variable_eq.size(); ++i) {
auto &var = variable_eq[i];
auto &cst = constant_eq[i];
assert((0 <= var.array()).all());
assert((var.array() < value.size()).all());
assert(cst[0] * value[var[0]] + cst[1] * value[var[1]] + cst[2] * value[var[2]] == 0);
}
for (int i = 0; i < (int)variable_ge.size(); ++i) {
auto &var = variable_ge[i];
auto &cst = constant_ge[i];
assert((0 <= variable_ge[i].array()).all());
assert((variable_ge[i].array() < value.size()).all());
if (value[var[0]] * value[var[1]] * cst[0] - value[var[2]] * value[var[3]] * cst[1] < 0) {
assert(face_area[i] < 0);
flip_count1++;
}
}
assert(flip_count == flip_count1);
// BFS
printf("[SAT] Start BFS: Q.size() = %d\n", (int)Q.size());
int mark_count = Q.size();
while (!Q.empty()) {
int vertex = Q.front().first;
int depth = Q.front().second;
Q.pop_front();
mark_count++;
int e0 = V2E(vertex);
for (int e = e0;;) {
int v = F((e + 1) % 3, e / 3);
if (!mark_vertex[v]) {
int undirected_edge_id = face_edgeIds[e / 3][e % 3];
int undirected_edge_length = edge_diff[undirected_edge_id].array().abs().sum() > 0;
int ndepth = depth + undirected_edge_length;
if (ndepth <= max_depth) {
if (undirected_edge_length == 0)
Q.push_front(std::make_pair(v, ndepth));
else
Q.push_back(std::make_pair(v, ndepth));
mark_vertex[v] = true;
}
}
e = dedge_next_3(E2E(e));
if (e == e0) break;
}
}
printf("[SAT] Mark %d vertices out of %d\n", mark_count, (int)V2E.size());
std::vector<bool> flexible(value.size(), false);
for (int i = 0; i < (int)face_edgeIds.size(); ++i) {
for (int j = 0; j < 3; ++j) {
int edgeid = face_edgeIds[i][j];
if (mark_vertex[F(j, i)] || mark_vertex[F((j + 1) % 3, i)]) {
flexible[edgeid * 2 + 0] = true;
flexible[edgeid * 2 + 1] = true;
} else {
assert(face_area[i] >= 0);
}
}
}
SolveSatProblem(value.size(), value, flexible, variable_eq, constant_eq, variable_ge,
constant_ge);
for (int i = 0; i < edge_diff.size(); ++i) {
edge_diff[i][0] = value[2 * i + 0];
edge_diff[i][1] = value[2 * i + 1];
}
}
} // namespace qflow