1
0
mirror of https://github.com/zenorogue/hyperrogue.git synced 2024-12-11 03:50:27 +00:00
hyperrogue/rogueviz/som/voronoi.cpp
2022-04-21 11:56:01 +02:00

259 lines
7.1 KiB
C++

// Voronoi used to measure the quality of the embedding (Villman's measure)
// Copyright (C) 2011-2022 Tehora and Zeno Rogue, see 'hyper.cpp' for details
#include "kohonen.h"
namespace rogueviz {
namespace voronoi {
void manifold::generate_data() {
T = isize(triangles);
triangles_of_tile.resize(N);
for(int i=0; i<T; i++) {
for(auto v: triangles[i])
triangles_of_tile[v].push_back(i);
for(int j=0; j<3; j++) {
int e0 = triangles[i][j%3];
int e1 = triangles[i][(j+1)%3];
if(e1<e0) swap(e0, e1);
auto p = make_pair(e0, e1);
triangles_of_edge[p].push_back(i);
}
}
}
manifold build_manifold(const vector<cell*>& cells) {
map<cell*, int> neuron_id;
int N = isize(cells);
for(int i=0; i<N; i++)
neuron_id[cells[i]] = i;
set<vector<int> > faces_seen;
for(auto c: cells) {
for(int i=0; i<c->type; i++) {
cellwalker cw1(c, i);
cellwalker cw2 = cw1;
vector<int> tlist;
do {
cw2 += wstep;
cw2++;
auto p = at_or_null(neuron_id, cw2.at);
if(!p) goto next;
tlist.push_back(*p);
}
while(cw2 != cw1);
if(1) {
int minv = 0;
for(int i=0; i<isize(tlist); i++)
if(tlist[i] < tlist[minv])
minv = i;
vector<int> tlist_sorted;
for(int i=minv; i<isize(tlist); i++)
tlist_sorted.push_back(tlist[i]);
for(int i=0; i<minv; i++)
tlist_sorted.push_back(tlist[i]);
if(tlist_sorted[1] > tlist_sorted.back())
reverse(tlist_sorted.begin()+1, tlist_sorted.end());
faces_seen.insert(tlist_sorted);
}
next: ;
}
}
manifold m;
m.N = N;
for(const auto& v: faces_seen)
for(int i=2; i<isize(v); i++)
m.triangles.emplace_back(make_array(v[0], v[i-1], v[i]));
m.generate_data();
return m;
}
vector<pair<int, int> > compute_voronoi_edges(manifold& m) {
using kohonen::net;
using kohonen::vnorm;
using kohonen::vshift;
using kohonen::data;
using kohonen::kohvec;
using kohonen::samples;
vector<int> best_tile;
/* for every neuron, compute its best tile */
int N = isize(net);
for(int i=0; i<N; i++) {
ld bestval = HUGE_VAL, best_id = -1;
for(int j=0; j<samples; j++) {
ld val = vnorm(net[i].net, data[j].val);
if(val < bestval) bestval = val, best_id = j;
}
best_tile.push_back(best_id);
}
constexpr int SUBD = 8;
using neuron_triangle_pair = pair<int, int>;
set< neuron_triangle_pair > visited;
queue<neuron_triangle_pair> q;
vector<kohvec> projected(N);
auto visit = [&] (neuron_triangle_pair p) {
if(visited.count(p)) return;
visited.insert(p);
q.push(p);
};
kohvec at;
kohonen::alloc(at);
auto project = [&] (kohvec& at, const array<int, 3>& tri, int i, int j) {
int k = SUBD-i-j;
for(auto& x: at) x = 0;
vshift(at, data[tri[0]].val, i * 1. / SUBD);
vshift(at, data[tri[1]].val, j * 1. / SUBD);
vshift(at, data[tri[2]].val, k * 1. / SUBD);
};
set<kohvec> already_picked;
map<int, string> which_best;
/* project all the net[ni].net on the manifold */
for(int ni=0; ni<N; ni++) {
kohvec best;
int best_tri;
ld best_dist = HUGE_VAL;
reaction_t better = [] {};
set<int> triangles_to_visit;
queue<int> triangles_queue;
auto visit1 = [&] (int tri) {
if(triangles_to_visit.count(tri)) return;
triangles_to_visit.insert(tri);
triangles_queue.push(tri);
};
for(int tr: m.triangles_of_tile[best_tile[ni]])
visit1(tr);
auto& bes = which_best[ni];
while(!triangles_queue.empty()) {
int tri = triangles_queue.front();
triangles_queue.pop();
for(int i=0; i<=SUBD; i++)
for(int j=0; j<=SUBD-i; j++) {
project(at, m.triangles[tri], i, j);
ld dist = vnorm(at, net[ni].net);
if(dist < best_dist && !already_picked.count(at)) {
best_dist = dist, best = at, best_tri = tri;
bes = lalign(0, tie(tri, i, j));
better = [&tri, i, j, &m, &visit1] () {
auto flip_edge = [&] (int t1, int t2) {
if(t2 < t1) swap(t1, t2);
for(auto tri1: m.triangles_of_edge[{t1, t2}])
visit1(tri1);
};
auto& tria = m.triangles[tri];
if(i == 0) flip_edge(tria[1], tria[2]);
if(j == 0) flip_edge(tria[0], tria[2]);
if(i+j == SUBD) flip_edge(tria[0], tria[1]);
};
}
}
better();
}
projected[ni] = best;
already_picked.insert(best);
visit({ni, best_tri});
}
struct triangle_data {
double dist[SUBD+1][SUBD+1];
int which[SUBD+1][SUBD+1];
triangle_data() {
for(int i=0; i<=SUBD; i++)
for(int j=0; j<=SUBD; j++)
dist[i][j] = HUGE_VAL, which[i][j] = -1;
}
};
vector<triangle_data> tdatas(m.T);
while(!q.empty()) {
auto ntp = q.front();
q.pop();
auto ni = ntp.first;
auto& tri = m.triangles[ntp.second];
auto& td = tdatas[ntp.second];
for(int i=0; i<=SUBD; i++)
for(int j=0; j<=SUBD-i; j++) {
project(at, tri, i, j);
ld dist = vnorm(at, projected[ni]);
auto& odist = td.dist[i][j];
bool tie = abs(dist - odist) < 1e-6;
if(tie ? ni < td.which[i][j] : dist < odist) {
td.dist[i][j] = dist,
td.which[i][j] = ni;
auto flip_edge = [&] (int t1, int t2) {
if(t2 < t1) swap(t1, t2);
for(auto tr: m.triangles_of_edge[{t1, t2}]) {
visit({ni, tr});
}
};
if(i == 0) flip_edge(tri[1], tri[2]);
if(j == 0) flip_edge(tri[0], tri[2]);
if(i+j == SUBD) flip_edge(tri[0], tri[1]);
}
}
}
set<pair<int, int> > voronoi_edges;
auto add_edge = [&] (int i, int j) {
if(i>j) swap(i, j);
if(i==j) return;
voronoi_edges.insert({i, j});
};
for(auto& td: tdatas) {
for(int i=0; i<=SUBD; i++)
for(int j=0; j<=SUBD-i; j++) {
if(i>0) add_edge(td.which[i][j], td.which[i-1][j]);
if(j>0) add_edge(td.which[i][j], td.which[i][j-1]);
if(j>0) add_edge(td.which[i][j], td.which[i+1][j-1]);
}
}
if(1) {
vector<int> degs(N, 0);
for(auto e: voronoi_edges) degs[e.first]++, degs[e.second]++;
for(int v=0; v<N; v++) if(degs[v] == 0) {
fhstream vorerr("voronoi-error.txt", "at");
println(vorerr, "error: degree 0 vertex ", v, " in ", debug_str);
println(vorerr, "best is: ", which_best[v]);
int id = 0;
for(auto& td: tdatas) {
for(int i=0; i<=SUBD; i++)
for(int j=0; j<=SUBD-i; j++)
if(td.which[i][j] == v) println(vorerr, "Found at ", tie(id, i, j));
id++;
}
}
}
vector<pair<int, int> > result;
for(auto ve: voronoi_edges) result.push_back(ve);
return result;
}
}
}