diff options
Diffstat (limited to 'extern/quadriflow/src/localsat.cpp')
| -rw-r--r-- | extern/quadriflow/src/localsat.cpp | 295 |
1 files changed, 295 insertions, 0 deletions
diff --git a/extern/quadriflow/src/localsat.cpp b/extern/quadriflow/src/localsat.cpp new file mode 100644 index 00000000000..1ab665175da --- /dev/null +++ b/extern/quadriflow/src/localsat.cpp @@ -0,0 +1,295 @@ +#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 |