// 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; } } }