rv::som:: added our tests

rogueviz/som/analyzer.cpp

@@ -0,0 +1,229 @@
+// analyze the SOM test results
+// Copyright (C) 2011-2022 Tehora and Zeno Rogue, see 'hyper.cpp' for details
+
+#include "kohonen.h"
+
+namespace rogueviz { 
+
+namespace kohonen_test {
+extern string cg();
+extern int data_scale, embed_scale;
+extern int subdata_value;
+}
+
+double UNKNOWN_RESULT = -3;
+
+using namespace kohonen_test;
+
+namespace analyzer {
+
+using measures::manidata;
+using measures::MCOUNT;
+
+struct maniset {
+  vector<string> names;
+  map<string, manidata> mdata;
+  };
+
+maniset origs, embs;
+
+void load_maniset(maniset& m, int scale) {
+  FILE *f = fopen((som_test_dir + "edgelists-" + its(scale) + ".txt").c_str(), "r");
+  while(true) {
+    char buf[200];
+    if(fscanf(f, "%s", buf) <= 0) break;
+    if(buf[0] != '#')
+      m.names.push_back(buf);
+    auto& md = m.mdata[buf];
+    int N, M;
+    fscanf(f, "%d%d", &N, &M);
+    println(hlog, "reading ", buf, " of size ", N, " and ", M, " edges");
+    md.size = N;
+    auto& ed = md.edges;
+    for(int i=0; i<M; i++) {
+      int a, b;
+      fscanf(f, "%d%d", &a, &b);
+      ed.emplace_back(a, b);
+      }
+    md.distances = measures::build_distance_matrix(N, md.edges);
+    }
+  fclose(f);
+  }
+
+vector<pair<int, int> > vor_edges;
+
+void analyze_test() {
+
+  println(hlog, "loading original manifolds");
+  load_maniset(origs, data_scale);
+  println(hlog, "loading embedding manifolds");
+  load_maniset(embs, embed_scale);
+
+  map<string, map<int, array<double, MCOUNT>>> results;
+
+  string tablefile = som_test_dir + "table" + cg() + ".csv";
+  
+  if(1) {
+    
+    fhstream res(tablefile, "rt");
+    if(!res.f) { println(hlog, "table ", tablefile, " missing"); }
+    else {
+      println(hlog, "reading the old table ", tablefile);
+      string s = scanline(res);
+      while(true) {
+        s = scanline(res);
+        if(s == "") break;
+        int i = 0;
+        while(s[i] != ';') i++;
+        i++;
+        while(s[i] != ';') i++;
+        i++;
+        int id;
+        sscanf(s.c_str()+i, "%d", &id);
+        auto& res = results[s.substr(0, i-1)][id];
+        for(auto& d: res) d = UNKNOWN_RESULT;
+        sscanf(s.c_str()+i, "%d;%lf;%lf;%lf;%lf;%lf;%lf;%lf;%lf;%lf;%lf;%lf;%lf", &id, &res[0], &res[1], &res[2], &res[3], &res[4], &res[5], &res[6], &res[7], &res[8], &res[9], &res[10], &res[11]);
+        }
+      }
+    }
+  
+  for(string orig: origs.names) {
+    
+    string dir = som_test_dir + "pairs" + cg();
+  
+    for(string emb: embs.names) {
+      string fname = dir + "/" + orig + "-" + emb + ".txt";
+      if(!file_exists(fname)) continue;
+      FILE *f = fopen(fname.c_str(), "rt");
+
+      string vorname = dir + "/" + orig + "-" + emb + ".vor";
+      if(!file_exists(vorname)) continue;
+      FILE *vor = fopen(vorname.c_str(), "rt");
+
+      int No = origs.mdata[orig].size;
+      auto& edo = origs.mdata[orig].edges;
+      vector<pair<int, int>> edo_recreated;
+      
+      int new_results = 0;
+
+      int current_row = -1;
+
+      auto& res1 = results[orig + ";" + emb];
+      for(int it=0; it<100; it++) {
+
+        bool know_result = res1.count(it);
+        auto& res = res1[it];
+        if(!know_result)
+          for(int kk=0; kk<MCOUNT; kk++)
+            res[kk] = UNKNOWN_RESULT;
+
+        vector<int> mapp(No, -2);
+        for(int k=0; k<MCOUNT; k++) {        
+          bool recompute = false;
+          ld energy = 0;
+          if(res[k] == UNKNOWN_RESULT || recompute) {
+
+            while(current_row < it) {
+              current_row++;
+              for(int i=0; i<No; i++) fscanf(f, "%d", &mapp[i]);
+              int V = 0; fscanf(vor, "%d", &V);
+              vor_edges.resize(V);
+              for(int i=0; i<V; i++) fscanf(vor, "%d%d", &vor_edges[i].first, &vor_edges[i].second);
+              if(mapp.back() == -2) goto next_pair;
+              if(current_row == it)
+                edo_recreated = measures::recreate_topology(mapp, edo);
+              }
+
+            new_results++;
+            
+            energy = evaluate_measure(embs.mdata[emb], origs.mdata[orig], mapp, vor_edges, edo_recreated, k);
+            
+            if(recompute && res[k] != UNKNOWN_RESULT) {
+              if(abs(res[k] - energy) > 1e-5) {
+                println(hlog, "ERROR in ", orig, "->", emb, " in ", cg(), " index ", it, " : was ", res[k], " is ", energy);
+                if(subdata_value) res[k] = energy;
+                }
+              }
+            else {
+              res[k] = energy;
+              }
+            }
+          }
+        }
+      next_pair:
+      fclose(f);
+      
+      if(new_results) println(hlog, "computed ", new_results, " new results in ", cg(), " for ", orig, " -> ", emb);
+      }
+    }
+
+  fhstream o(som_test_dir + "summary" + cg() + ".html", "wt");
+  println(o, "<html><body>\n");
+  
+  if(1) {
+    println(o, "<table>\n");
+    
+    int manis = 0;
+    
+    for(string orig: origs.names) {
+
+    print(o, "<tr><th>", orig, "</th>");
+    for(string emb: embs.names) print(o, "<th>", emb, "</th>");
+    print(o, "</tr>\n");
+    
+    for(int k: {0, 2, 4, 7, 8, 10, 11}) {
+      print(o, "<tr><td>", measures::catnames[k], "</td>");
+      
+      for(string emb: embs.names) {
+        ld tenergy = 0;
+        int att = 0;
+        auto& res1 = results[orig + ";" + emb];
+
+        for(int it=0; it<100; it++) {
+          if(!res1.count(it)) continue;
+          auto& res = res1[it];
+          if(res[k] == UNKNOWN_RESULT) continue;
+          att++; tenergy += res[k];
+          }
+    
+        if(orig == emb)
+          print(o, "<td bgcolor=\"#C0FFFF\">");
+        else if(manis & 1)
+          print(o, "<td bgcolor=\"#FFFFC0\">");
+        else
+          print(o, "<td>");
+        if(att == 100)
+          print(o, tenergy * 1. / att);
+        else if(att)
+          print(o, tenergy * 1. / att, "?");
+        else
+          print(o, "???");
+        print(o, "</td>");
+        }
+      print(o, "</tr>\n");
+      }
+      manis++;
+      }
+    print(o, "</table>\n");
+    }
+  print(o, "</body></html>\n");
+  
+  fhstream res(som_test_dir + "table" + cg() + ".csv", "wt");
+  print(res, "orig;emb;index");
+  for(int i=0; i<MCOUNT; i++) print(res, ";", measures::catnames[i]); println(res);
+  
+  for(auto& p: results) {
+    for(auto& p2: p.second) {
+      auto &p3 = p2.second;
+      print(res, p.first, ";", p2.first);
+      for(int i=0; i<MCOUNT; i++) print(res, ";", p3[i]); println(res);
+      }
+    }  
+  }
+
+auto khook = arg::add3("-analyze-test", analyze_test);
+
+}}
+
+
+ 

rogueviz/som/kohonen.h

@@ -20,6 +20,8 @@ static const flagtype KS_NEURONS_INI = 16;
 
 extern flagtype state;
 
+extern map<int, int> sample_vdata_id;
+
 extern int tmax;
 extern int samples;
 extern int t, tmax, lpct, cells;
@@ -33,6 +35,7 @@ extern int dispersion_count;
 extern double learning_factor, distmul;
 extern double ttpower;
 extern int min_group, max_group, columns;
+extern bool dispersion_long;
 
 struct neuron {
   kohvec net;
@@ -58,6 +61,8 @@ extern vector<int> samples_to_show;
 extern vector<neuron> net;
 extern vector<string> colnames;
 
+extern vector<int> sample_sequence;
+
 void initialize_neurons();
 void initialize_neurons_initial();
 void initialize_dispersion();
@@ -80,6 +85,8 @@ double vnorm(kohvec& a, kohvec& b);
 
 namespace embeddings {
 
+ld gaussian_random();
+
 using kohonen::columns;
 using kohonen::kohvec;
 using kohonen::alloc;
@@ -90,14 +97,59 @@ void mark_signposts(bool full, const vector<cell*>& ac);
 void mark_signposts_subg(int a, int b, const vector<cell*>& ac);
 void generate_rug(int i, bool close);
 void init_landscape(int dimensions);
+void init_landscape_det(const vector<cell*>& ac);
 
 extern map<cell*, hyperpoint> rug_coordinates;
 extern void get_coordinates(kohvec& v, cell *c, cell *c0);
 
 extern vector<cell*> signposts;
 
+extern map<cellwalker, kohvec> delta_at;
+
 }
 
+namespace voronoi {
+
+struct manifold {
+
+  int N; /* the number of tiles */
+  int T; /* the number of triangles */
+
+  /* triangles between three adjacent tiles. triangles[i] are the tiles of triangle i */
+  vector<array<int, 3>> triangles;
+
+  /* triangles_of_tile[t] is ids of triangles which contain t */
+  vector<vector<int>> triangles_of_tile;
+
+  map<pair<int, int>, vector<int>> triangles_of_edge;
+
+  void generate_data();
+  };
+
+inline string debug_str;
+manifold build_manifold(const vector<cell*>& cells);
+vector<pair<int, int> > compute_voronoi_edges(manifold& m);
+}
+
+namespace measures {
+
+struct manidata {
+  int size;
+  vector<vector<int> > distances;
+  vector<pair<int, int> > edges;
+  };
+
+static const int MCOUNT = 12;
+
+extern vector<string> catnames;
+
+vector<pair<int, int>> recreate_topology(const vector<int>& mapp, const vector<pair<int, int> >& edges);
+vector<vector<int>> build_distance_matrix(int N, const vector<pair<int,int>>& vedges);
+ld evaluate_measure(manidata& emb, manidata& orig, vector<int>& mapp, vector<pair<int, int> >& vor_edges, vector<pair<int, int>>& edo_recreated, int k);
+
+}
+
+static const string som_test_dir = "results/";
 }
 
 

rogueviz/som/measures.cpp

@@ -0,0 +1,239 @@
+// test measures
+// Copyright (C) 2011-2022 Tehora and Zeno Rogue, see 'hyper.cpp' for details
+
+#include "kohonen.h"
+
+namespace rogueviz {
+namespace measures {
+
+double kendall(const vector<pair<int, int>>& allp) {
+  int maxo = 0, maxe = 0;
+  for(const auto& a: allp) maxo = max(maxo, a.first), maxe = max(maxe, a.second);
+  maxo++; maxe++;
+  
+  if(maxo >= 256 || maxe >= 256) {
+    println(hlog, "more than 256!\n");
+    exit(1);
+    }
+  int cnt[256][256];
+  for(int a=0; a<maxo; a++)
+  for(int b=0; b<maxe; b++)
+    cnt[a][b] = 0;
+  for(const auto& a: allp) cnt[a.first][a.second]++;
+
+  // int i1 = 0, i2 = 0;
+  int K = isize(allp);
+  // allp.emplace_back(maxo, maxe);
+  vector<int> counts(maxe, 0);
+  vector<int> totals(maxe);
+  double tau = 0;
+  for(int i=0; i<maxo; i++) {
+    totals[0] = 0;
+    for(int ii=1; ii<maxe; ii++)
+      totals[0] -= counts[ii];
+    for(int ii=1; ii<maxe; ii++)
+      totals[ii] = totals[ii-1] + counts[ii] + counts[ii-1];
+    
+    for(int b=0; b<maxe; b++) {
+      tau += totals[b] * 1. * cnt[i][b];
+      counts[b] += cnt[i][b];
+      }
+    }
+  ld par = (K * (K-1.) / 2);
+  return tau / par;
+  }
+
+vector<pair<int, int>> recreate_topology(const vector<int>& mapp, const vector<pair<int, int> >& edges) {
+
+  auto cmapp = mapp;
+  for(int i=0; i<isize(cmapp); i++) if(cmapp[i] >= 0) cmapp[i] = i;
+  
+  while(true) {
+    vector<pair<int, int> > changes;
+    for(auto e: edges) {
+      if(cmapp[e.first] == -1 && cmapp[e.second] >= 0) changes.emplace_back(e.first, cmapp[e.second]);
+      if(cmapp[e.second] == -1 && cmapp[e.first] >= 0) changes.emplace_back(e.second, cmapp[e.first]);
+      }
+    if(changes.empty()) break;
+    for(auto ch: changes) cmapp[ch.first] = ch.second;
+    }
+
+  set<pair<int, int>> subedges;
+  for(auto e: edges) {
+      auto a = cmapp[e.first], b = cmapp[e.second];
+      if(a==b) continue;
+      if(a<b) swap(a, b);
+      subedges.emplace(a, b);
+      }
+
+  vector<pair<int, int>> result;
+  for(auto sube: subedges)
+    result.emplace_back(sube.first, sube.second);
+
+  return result;
+  }
+
+vector<vector<int>> build_distance_matrix(int N, const vector<pair<int,int>>& vedges) {
+  vector<vector<int>> res;
+  res.resize(N);
+  for(auto& r: res) r.resize(N, 1000);
+  vector<vector<int>> edges_of(N);
+  for(auto [a, b]: vedges) edges_of[a].push_back(b), edges_of[b].push_back(a);  
+   
+  for(int i=0; i<N; i++) {
+    auto& d = res[i];
+    queue<int> bfs;
+    auto visit = [&] (int id, int dist) {
+      if(d[id] <= dist) return;
+      d[id] = dist;
+      bfs.push(id);
+      };
+    visit(i, 0);
+    while(!bfs.empty()) {
+      auto j = bfs.front();
+      for(auto e: edges_of[j]) visit(e, d[j]+1);
+      bfs.pop();
+      }
+    }
+  
+  return res;
+  }
+
+vector<string> catnames = {"energy", "rips", "kendall", "empty", "emptyx", "subedges", "emptyxx", "villman1", "villman2", "vedges", "maxvill1", "maxvill2"};
+
+ld evaluate_measure(manidata& emb, manidata& orig, vector<int>& mapp, vector<pair<int, int> >& vor_edges, vector<pair<int, int>>& edo_recreated, int k) {
+
+  int No = orig.size;
+  int Ne = emb.size;
+  auto& edo = orig.edges;
+  auto& ede = emb.edges;
+  auto& dise = emb.distances;
+  auto& diso = orig.distances;
+  ld energy = 0;
+
+  if(k == 2) {
+    vector<pair<int, int> > allp;
+    allp.reserve(No * No);
+    for(int i=0; i<No; i++) if(mapp[i] != -1)
+    for(int j=0; j<i; j++) if(mapp[j] != -1)
+      allp.emplace_back(diso[i][j], dise[mapp[i]][mapp[j]]);
+    energy = kendall(allp);
+    }
+  else if(k == 3) {
+    vector<bool> empty(Ne, true);
+    for(int i=0; i<No; i++) 
+      if(mapp[i] != -1) empty[mapp[i]] = false;
+    for(auto b: empty) if(b) energy++;
+    }
+  else if(k == 4 || k == 6) {
+    vector<vector<int>> adj(Ne);
+    vector<vector<int>> on(Ne);
+    for(int i=0; i<No; i++) if(mapp[i] != -1)
+      on[mapp[i]].push_back(i);
+    for(auto [a,b]: ede) adj[a].push_back(b), adj[b].push_back(a);
+    for(auto p: edo_recreated)
+      diso[p.first][p.second] -= 100,
+      diso[p.second][p.first] -= 100;
+    for(int i=0; i<Ne; i++) {
+      bool empty = on[i].empty();
+      if(empty) {
+        for(auto i1: adj[i])
+        for(auto i2: adj[i])
+        for(auto oi1: on[i1])
+        for(auto oi2: on[i2])
+        if(dise[i1][i2] > 1 && diso[oi1][oi2] <= 1)
+          empty = false;
+        }
+      if(empty && k == 6) {
+        for(auto i1: adj[i])
+        for(auto i2: adj[i])
+        for(auto i11: adj[i1]) {
+          for(auto i21: adj[i2])
+          for(auto oi1: on[i11])
+          for(auto oi2: on[i21])
+          if(dise[i11][i21] == dise[i11][i] + dise[i21][i] && diso[oi1][oi2] <= 1)
+            empty = false;
+
+          for(auto oi1: on[i11])
+          for(auto oi2: on[i2])
+          if(dise[i11][i2] == dise[i11][i] + dise[i2][i] && diso[oi1][oi2] <= 1)
+            empty = false;
+          }
+        }
+      if(empty) energy++;
+      }
+    for(auto p: edo_recreated)
+      diso[p.first][p.second] += 100,
+      diso[p.second][p.first] += 100;
+    }
+  else if(k == 5) energy = isize(edo_recreated);
+  else if(k == 1) {
+    for(auto [a,b]: edo_recreated) {
+      if(mapp[a] == -1 || mapp[b] == -1) continue;
+      int d = dise[mapp[a]][mapp[b]];
+      if(k == 9)
+        energy += d * d;
+      else if(k == 1)
+        energy += (d < 2 ? 0 : d == 2 ? 1 : 100);
+      }                
+    }
+  else if(k == 0) {
+    for(auto [a,b]: edo) {
+      if(mapp[a] == -1 || mapp[b] == -1) continue;
+      int d = dise[mapp[a]][mapp[b]];
+      energy += d * d;
+      }                
+    }
+  else if(k == 7) {
+    for(auto [a,b]: vor_edges) {
+      int d = dise[a][b] - 1;
+      energy += d*d;
+      }
+    }
+  else if(k == 8) {
+    auto disv = build_distance_matrix(Ne, vor_edges);
+    int bugs = 0;
+    for(auto [a,b]: ede) {
+      int d = disv[a][b] - 1;
+      if(d >= 900) bugs++;
+      energy += d*d;
+      }
+    if(bugs) println(hlog, "bugs=", bugs);
+    }
+  else if(k == 10) {
+    for(auto [a,b]: vor_edges)
+      energy = max<ld>(energy, dise[a][b] - 1);
+    }
+  else if(k == 11) {
+    auto disv = build_distance_matrix(Ne, vor_edges);
+    for(auto [a,b]: ede) 
+      energy = max<ld>(energy, disv[a][b] - 1);
+    }
+  else if(k == 9)
+    energy = isize(vor_edges);  
+  return energy;
+  }
+
+/*
+void test_kendall() {
+  for(string orig: origs.names) {
+    int No = origs.mdata[orig].size;
+    vector<pair<int, int> > allp;
+    auto& diso = origs.mdata[orig].distances;
+    
+    for(int i=0; i<No; i++)
+    for(int j=0; j<i; j++)
+      allp.emplace_back(diso[i][j], diso[i][j]);
+      
+    allp.emplace_back(0, 0);
+    allp.emplace_back(1, 1);
+    allp.emplace_back(2, 2);
+    allp.emplace_back(3, 3);
+
+    println(hlog, orig, " : ", measures::kendall(allp));
+    }  
+  }
+*/
+
+}
+}

rogueviz/som/voronoi.cpp

@@ -0,0 +1,258 @@
+// 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;  
+  }
+
+}
+}