Skip to content

Instantly share code, notes, and snippets.

@hakomo

hakomo/PointsOnTheCircle.cpp

Last active February 1, 2018 03:32
Show Gist options
  • Save hakomo/23c07ab9b823324a5ef6e6581a7fa6b3 to your computer and use it in GitHub Desktop.
Save hakomo/23c07ab9b823324a5ef6e6581a7fa6b3 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
PointsOnTheCircle.cpp
#include <algorithm>
#include <cmath>
#include <functional>
#include <iostream>
#include <vector>
using namespace std;
#define U(v) cerr << #v << ": " << (v) << endl
double getTime() {
unsigned long long a, d;
__asm__ volatile ("rdtsc" : "=a" (a), "=d" (d));
return (d << 32 | a) / 2.8e9;
}
int rando(int b) {
static unsigned y = 2463534242;
y ^= (y ^= (y ^= y << 13) >> 17) << 5;
return y % b;
}
int N;
double sint[400];
double diffs[200];
double transitionTh[256];
vector<pair<int, int>> edges;
vector<unsigned char> valueToValues[200];
double valueToScore[200];
int valueToIndex[200];
int indexToValue[200];
int bestIndexToValue[200];
int rangei[100];
int ranget[100][200];
void init() {
for (int i = 0; i < N; ++i)
valueToIndex[indexToValue[i]] = i;
fill_n(valueToScore, N, 0);
for (int v = 0; v < N; ++v)
for (unsigned char w : valueToValues[v])
valueToScore[v] += sint[valueToIndex[w] - valueToIndex[v] + 200];
random_shuffle(edges.begin(), edges.end(), rando);
for (int d = 1; d < N / 2; ++d)
random_shuffle(ranget[d], ranget[d] + d * 2, rando);
}
double eval() {
double score = 0;
for (int i = 0; i < N; ++i)
for (int j : valueToValues[i])
score += sint[valueToIndex[i] - valueToIndex[j] + 200];
return score / 2;
}
struct PointsOnTheCircle {
vector<int> permute(const vector<int>& matrix) const {
double startTime = getTime();
N = (int)sqrt(matrix.size());
{
double pi = atan(1) * 4;
for (int i = 1; i < N; ++i)
sint[200 + i] = sint[200 - i] = sin(pi / N * i);
for (int i = 1; i < N; ++i)
diffs[i] = sin(pi / N * (i + 1)) - sin(pi / N * i);
}
for (int d = 1; d < N / 2; ++d) {
int i = 0;
for (int e = 1; e <= d; ++e) {
ranget[d][i++] = e;
ranget[d][i++] = -e;
}
}
for (int y = 0; y < N; ++y) {
for (int x = 0; x < N; ++x) {
if (matrix[y * N + x]) {
valueToValues[y].emplace_back(x);
edges.emplace_back(x, y);
}
}
}
double score;
double bestScore = 1e9;
const int MultiStartNumber = (int)(25 * pow(40 / max(1.0, edges.size() / (double)N), 0.25));
for (int _ = MultiStartNumber; _--; ) {
if (_) {
for (int i = 0; i < N; ++i)
indexToValue[i] = i;
random_shuffle(indexToValue, indexToValue + N, rando);
if (edges.empty())
return { indexToValue, indexToValue + N };
} else {
copy_n(bestIndexToValue, N, indexToValue);
}
init();
score = eval();
if (bestScore > score) {
bestScore = score;
copy_n(indexToValue, N, bestIndexToValue);
}
int ei = -1;
for (int iter = -1; ; ) {
if ((++iter & 0x1fff) == 0) {
double p = 1 - (getTime() - startTime - 8.85 / (MultiStartNumber - 1) * (MultiStartNumber - _ - 1)) / (_ ? 8.85 / (MultiStartNumber - 1) : 1);
if (p < 0) break;
const double InitialHeat = _ ? 0.4 : 0.1;
double heat = -p * InitialHeat;
for (int i = 0; i < 256; ++i)
transitionTh[i] = log((i + 1) / 257.0) * heat;
random_shuffle(transitionTh, transitionTh + 256, rando);
}
if (++ei == edges.size()) ei = 0;
int v1 = edges[ei].first;
int v2 = edges[ei].second;
int i1 = valueToIndex[v1];
int i2 = valueToIndex[v2];
int d = abs(i1 - i2);
d = min(d, N - d) - 1;
if (d == 0) continue;
if (++rangei[d] == d << 1) rangei[d] = 0;
i2 += ranget[d][rangei[d]];
if (i2 < 0) {
i2 += N;
} else if (i2 >= N) {
i2 -= N;
}
v2 = indexToValue[i2];
auto sin1 = sint - i1 + 200;
auto sin2 = sint - i2 + 200;
double d1, d2;
d1 = d2 = matrix[v1 * N + v2] ? sin1[i2] : 0;
for (unsigned char v : valueToValues[v1])
d1 += sin2[valueToIndex[v]];
for (unsigned char v : valueToValues[v2])
d2 += sin1[valueToIndex[v]];
double diff = d1 + d2 - valueToScore[v1] - valueToScore[v2];
if (diff < transitionTh[iter & 255]) {
score += diff;
valueToIndex[v1] = i2;
valueToIndex[v2] = i1;
indexToValue[i1] = v2;
indexToValue[i2] = v1;
for (unsigned char v : valueToValues[v2])
valueToScore[v] += sin1[valueToIndex[v]] - sin2[valueToIndex[v]];
for (unsigned char v : valueToValues[v1])
valueToScore[v] += sin2[valueToIndex[v]] - sin1[valueToIndex[v]];
valueToScore[v1] = d1;
valueToScore[v2] = d2;
if (bestScore > score) {
bestScore = score;
copy_n(indexToValue, N, bestIndexToValue);
}
}
}
}
copy_n(bestIndexToValue, N, indexToValue);
for (int i = 0; i < N; ++i)
valueToIndex[indexToValue[i]] = i;
score = bestScore;
while (true) {
bool updated = false;
for (int i = 0; i < N; ++i) {
pair<double, int> best{ 0, i };
double diff = 0;
for (int j = i + 1; j < N; ++j) {
for (unsigned char v : valueToValues[indexToValue[j - 1]]) {
if (valueToIndex[v] < j) {
diff += diffs[j - valueToIndex[v] - 1];
} else {
diff -= diffs[valueToIndex[v] - j];
}
}
for (unsigned char v : valueToValues[indexToValue[j]]) {
if (valueToIndex[v] > j - 1) {
diff += diffs[valueToIndex[v] - j];
} else {
diff -= diffs[j - valueToIndex[v] - 1];
}
}
best = min(best, { diff, j });
++valueToIndex[indexToValue[j - 1]];
--valueToIndex[indexToValue[j]];
swap(indexToValue[j - 1], indexToValue[j]);
}
for (int j = N - 1; j > i; --j) {
++valueToIndex[indexToValue[j - 1]];
--valueToIndex[indexToValue[j]];
swap(indexToValue[j - 1], indexToValue[j]);
}
diff = 0;
for (int j = i; j; --j) {
for (unsigned char v : valueToValues[indexToValue[j - 1]]) {
if (valueToIndex[v] < j) {
diff += diffs[j - valueToIndex[v] - 1];
} else {
diff -= diffs[valueToIndex[v] - j];
}
}
for (unsigned char v : valueToValues[indexToValue[j]]) {
if (valueToIndex[v] > j - 1) {
diff += diffs[valueToIndex[v] - j];
} else {
diff -= diffs[j - valueToIndex[v] - 1];
}
}
best = min(best, { diff, j - 1 });
++valueToIndex[indexToValue[j - 1]];
--valueToIndex[indexToValue[j]];
swap(indexToValue[j - 1], indexToValue[j]);
}
score += best.first;
for (int j = 1; j <= best.second; ++j) {
++valueToIndex[indexToValue[j - 1]];
--valueToIndex[indexToValue[j]];
swap(indexToValue[j - 1], indexToValue[j]);
}
updated |= best.first < -1e-9;
if (getTime() - startTime > 9.95) break;
}
constexpr int L = 6;
for (int beg = 0; beg <= N - L; ++beg) {
int end = beg + L;
double ps[L][L]{};
for (int i = beg; i < end; ++i) {
for (int j = beg; j < end; ++j) {
for (unsigned char v : valueToValues[indexToValue[i]]) {
if (valueToIndex[v] < beg) {
ps[i - beg][j - beg] += sint[j - valueToIndex[v] + 200];
} else if (valueToIndex[v] >= end) {
ps[i - beg][j - beg] += sint[valueToIndex[v] - j + 200];
}
}
}
}
double best = 1e9;
int currs[L];
int bests[L];
static vector<int> ii(L);
for (int i = 0; i < L; ++i)
ii[i] = i;
int v2i[L];
static vector<pair<int, int>> es;
es.clear();
for (int i = 1; i < L; ++i) {
for (int j = 0; j < i; ++j) {
if (matrix[indexToValue[beg + i] * N + indexToValue[beg + j]])
es.emplace_back(j, i);
}
}
function<void (int, double)> dfs = [&](int d, double s) {
if (d == L) {
for (auto& e : es)
s += sint[v2i[e.first] - v2i[e.second] + 200];
if (best > s) {
best = s;
copy_n(currs, L, bests);
}
} else {
for (int i = 0; i < (int)ii.size(); ++i) {
int iii = ii[i];
double ss = s + ps[iii][d];
v2i[iii] = d;
currs[d] = indexToValue[beg + iii];
ii[i] = ii.back();
ii.pop_back();
dfs(d + 1, ss);
ii.emplace_back(ii[i]);
ii[i] = iii;
}
}
};
dfs(0, 0);
for (int i = 0; i < L; ++i) {
if (indexToValue[beg + i] != bests[i]) {
indexToValue[beg + i] = bests[i];
updated = true;
}
}
for (int i = beg; i < end; ++i)
valueToIndex[indexToValue[i]] = i;
if (getTime() - startTime > 9.95) break;
}
for (int i = 2; i < N; ++i) {
for (int j = 0; j < i - 1; ++j) {
double diff = matrix[indexToValue[i] * N + indexToValue[j]] ? sint[i - j + 200] * 2 : 0;
for (unsigned char v : valueToValues[indexToValue[i]])
diff += sint[valueToIndex[v] - j + 200] - sint[valueToIndex[v] - i + 200];
for (unsigned char v : valueToValues[indexToValue[j]])
diff += sint[valueToIndex[v] - i + 200] - sint[valueToIndex[v] - j + 200];
if (diff < -1e-9) {
score += diff;
swap(valueToIndex[indexToValue[i]], valueToIndex[indexToValue[j]]);
swap(indexToValue[i], indexToValue[j]);
updated = true;
}
}
if (getTime() - startTime > 9.95) break;
}
if (!updated || getTime() - startTime > 9.95) break;
}
return { indexToValue, indexToValue + N };
}
};
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment