diff --git a/rogueviz/dhrg/betweenness.cpp b/rogueviz/dhrg/betweenness.cpp
new file mode 100644
index 00000000..251f7427
--- /dev/null
+++ b/rogueviz/dhrg/betweenness.cpp
@@ -0,0 +1,406 @@
+#include <iostream>
+#include <unordered_map>
+
+// pseudo-betweenness
+
+namespace dhrg {
+
+int out(int x, string cmt="") {
+  println(hlog, "result = ", x, " ", cmt);
+  return x;
+  }
+
+ld outNZ(ld x, string cmt="") { 
+  return x;
+  }
+
+int segmentlen(segment *s) {
+  auto l = s->left;
+  auto r = s->right;
+  int i = 1;
+  while(l != r) l = l->grightsibling(), i++;
+  return i;
+  }
+
+int get0(qtybox& b) { return b.minv == 0 ? b[0] : 0; }
+
+int tallybox_total(qtybox& box) {
+  int tt = 0;
+  for(int i=box.minv; i<box.maxv; i++)
+    tt += box[i];
+  return tt;
+  }
+
+string segdesc(segment *s) {
+  return format("(%s-%s: len=%d, qty=%d/%d)", get_path(s->left).c_str(), get_path(s->right).c_str(), segmentlen(s), get0(s->qty), tallybox_total(s->qty));
+  }
+
+
+int quickdist(segment *p1, segment *p2) {
+  int d = 0;
+  int d1 = p1->left->lev;
+  int d2 = p2->left->lev;
+  
+  while(d1>d2) { p1 = p1->parent; d1--; d++; }
+  while(d2>d1) { p2 = p2->parent; d2--; d++; }
+  
+  return segmentdist(p1, p2, d);
+  }
+
+pair<ll, ll> betweenness3(mycell *c, int setid = 0) {
+
+  vector<mycell*> neighbors;
+  neighbors = allchildren(c);
+  if(c->lev) {
+    neighbors.push_back(c->grightsibling());
+    neighbors.push_back(c->rightparent);
+    if(c->leftparent != c->rightparent)
+      neighbors.push_back(c->leftparent);
+    neighbors.push_back(c->gleftsibling());
+    }
+      
+  build_ack(c, setid);
+  for(auto c1: neighbors) build_ack(c1, setid);
+  
+  std::unordered_map<segment*, int> info;
+  
+  for(segment* p: acknowledged) {
+
+    int tot = tallybox_total(p->qty);
+
+    info[p] += tot;
+    
+    segment *p2 = p->parent;
+    while(p2) {
+      if(p2->seen != -1) {
+        info[p2] -= tot;
+        break;
+        }
+      p2=p2->parent;
+      }
+    }
+
+  for(segment* p: acknowledged) {
+    p->seen = -1;
+    }
+  
+  acknowledged.clear();
+
+  vector<pair<segment*, int> > info_v;
+  for(auto p: info) if(p.second) info_v.emplace_back(p);
+  
+  int NN = isize(neighbors);
+  int NN2 = 2 * NN;
+  
+  vector<ll> bydir(NN2, 0);
+  
+  int bestdist;  
+  vector<int> bestdir;
+
+  vector<segment*> neighbor_segs;
+  for(auto n: neighbors) neighbor_segs.push_back(getsegment(n, n, setid, true));
+  
+  for(auto p: info_v) {
+    if(p.first->left == c && p.first->right == c) continue;
+    bestdist = 1000;
+    bestdir.clear();
+    for(int k=0; k<NN; k++) {
+      int dist = quickdist(neighbor_segs[k], p.first);
+      if(dist < bestdist) bestdist = dist, bestdir.clear();
+      if(dist == bestdist) bestdir.push_back(k);
+      }    
+    if(isize(bestdir) == 1)
+      bydir[bestdir[0] * 2] += p.second;
+    else if(isize(bestdir) > 2) {
+      println(hlog, "many best dirs\n");
+      throw hr_exception("many best dirs");
+      }
+    else if(bestdir[0] + 1 == bestdir[1])
+      bydir[bestdir[0] + bestdir[1]] += p.second;
+    else if(bestdir[0] == 0 && bestdir[1] == NN-1)
+      bydir.back() += p.second;
+    else {
+      println(hlog, "non-adjacent best dirs\n");
+      for(int k=0; k<NN; k++) print(hlog, quickdist(neighbor_segs[k], p.first), " "); printf("\n");
+      println(hlog, "bestdir = ", bestdir);
+      throw hr_exception("non-adjacent best dirs");
+      }
+    }
+  
+  ll result = 0;
+  
+  int here = getsegment(c, c, setid, true)->qty[0];
+  
+  ll elsewhere = 0;
+  for(int a=0; a<NN2; a++) elsewhere += bydir[a];
+  
+  result += (here + elsewhere) * (here + elsewhere);
+
+  // for(int a=0; a<NN2; a++) printf("%3d ", bydir[a]); printf("HERE %d EW=%d N=%d\n", here, elsewhere, N);
+  
+  for(int a=0; a<NN2; a+=2) {
+    result -= bydir[a] * bydir[a];
+    result -= 2 * bydir[a] * bydir[(a+1) % NN2];
+    result -= 2 * bydir[a] * bydir[(a+2) % NN2];
+    result -= 2 * bydir[a] * bydir[(a+3) % NN2];
+    result -= 2 * bydir[a] * bydir[(a+NN2-1) % NN2];
+    result -= 2 * bydir[a] * bydir[(a+NN2-3) % NN2];
+    
+    result -= bydir[a+1] * bydir[(a+1) % NN2];
+    result -= 2 * bydir[a+1] * bydir[(a+3) % NN2];
+    result -= 2 * bydir[a+1] * bydir[(a+5) % NN2];
+    }
+    
+  ll result2 = here * here * NN + 2 * here * elsewhere * NN;
+  
+  for(int a=0; a<NN2; a++) if(bydir[a])
+  for(int b=0; b<a; b++) if(bydir[b])
+    result2 += 2 * bydir[a] * bydir[b] * min(a-b, NN2+b-a);
+  
+  return {result, result2};  
+  }
+
+bool dependent(segment *p1, segment *p2) {
+
+  if(is_insegment(p1->left, p2) || is_insegment(p1->right, p2) || is_insegment(p2->left, p1) || is_insegment(p2->right, p1))
+    return true;
+
+  mycell *mright;
+  
+  mright = p1->right;
+  for(int u=0; u<=cgi.expansion->sibling_limit; u++) {
+    mright->build();
+    if(mright == p2->left) return true;
+    mright = mright->grightsibling();
+    }
+
+  mright = p2->right;
+  for(int u=0; u<=cgi.expansion->sibling_limit; u++) {
+    mright->build();
+    if(mright == p1->left) return true;
+    mright = mright->grightsibling();
+    }
+  
+  return false;
+  }
+
+template<class T> void children(segment *s, const T& f) {
+  s=s->firstchild;
+  while(s) {
+    f(s);
+    s = s->nextchild;
+    }
+  }
+
+typedef long double betweenness_type;
+
+betweenness_type ack(int d01, int d02, int d12) {
+  return pow(cgi.expansion->get_growth(), -(d01 + d02 - d12));
+  // if(d01 + d02 == d12+2) return 1;
+  // return 0;
+  }
+
+int sd(segment *s1, segment *s2) {
+  /*
+  for(auto l = s1->left;; l = l->grightsibling()) {
+    if(is_insegment(l, s2)) return 0;
+    if(l == s1->right) break;
+    }
+  
+  for(auto l = s2->left;; l = l->grightsibling()) {
+    if(is_insegment(l, s1)) return 0;
+    if(l == s2->right) break;
+    }
+  */
+  
+  return quickdist(s1, s2);
+  }
+
+// di0: distance from Hidden to 0
+// di0: distance from Hidden to 1
+betweenness_type brec_fix_other(int di0, const vector<segment*>& lst, int pos, int di1, segment *s1) {
+  // indent_measure im("brec-fo " + its(pos) + " in s1: " + segdesc(s1) + " di0= " + its(di0) + " di1= " + its(di1) );
+  segment *s0 = lst[pos];
+  if(!dependent(s0, s1) || pos == 0)
+    return outNZ(tallybox_total(s1->qty) * ack(di0+pos, sd(s0, s1)+pos, di1));
+  
+  betweenness_type total = 0;
+  if(get0(s1->qty))
+    total += get0(s1->qty) * ack(di0+pos, sd(s0,s1)+pos, di1);
+
+  for(segment *c1 = s1->firstchild; c1; c1 = c1->nextchild)
+    total += brec_fix_other(di0+1, lst, pos-1, di1+1, c1);
+  
+  return outNZ(total);
+  }
+
+map<tuple<int, segment*, segment* >, betweenness_type> memo[BOXSIZE];
+
+betweenness_type brec_fix_main(int d1, segment *s1, int d2, segment *s2);
+
+betweenness_type brec_fix_main_actual(int d1, segment *s1, int d2, segment *s2) {
+  betweenness_type total = 0;
+  if(get0(s1->qty))
+    total += get0(s1->qty) * tallybox_total(s2->qty) * ack(d1, d2, sd(s1, s2));
+  if(get0(s2->qty))
+    total += get0(s2->qty) * tallybox_total(s1->qty) * ack(d1, d2, sd(s1, s2));
+  if(get0(s1->qty) && get0(s2->qty))
+    total -= get0(s1->qty) * get0(s2->qty) * ack(d1, d2, sd(s1, s2));
+  
+  for(segment *c1 = s1->firstchild; c1; c1 = c1->nextchild)
+  for(segment *c2 = s2->firstchild; c2; c2 = c2->nextchild)
+    total += brec_fix_main(d1+1, c1, d2+1, c2);
+  
+  return total;
+  }
+
+betweenness_type brec_fix_main(int d1, segment *s1, int d2, segment *s2) {
+  // indent_measure im("brec-main " + its(d1) + " in s1: " + segdesc(s1) + " " + its(d2) + " in s2: " + segdesc(s2));
+  if(!dependent(s1, s2))
+    return outNZ(tallybox_total(s1->qty) * tallybox_total(s2->qty) * ack(d1, d2, sd(s1, s2)));
+
+  // if(s1->left->lev >= 10) return brec_fix_main_actual(d1, s1, d2, s2);
+
+  auto& mem = memo[s1->left->lev][make_tuple(d1+d2, s1, s2)];
+  if(mem) return mem-1;
+
+  auto total = brec_fix_main_actual(d1, s1, d2, s2);
+  mem = total + 1;
+  return total;
+  }
+
+betweenness_type brec(const vector<segment*>& lst, int pos, segment *s1, segment *s2) {
+  segment *s0 = lst[pos];
+  // indent_measure im("brec " + its(pos) + " in s1: " + segdesc(s1) + " in s2: " + segdesc(s2) + " s0 = " + segdesc(s0) );
+  bool id01 = !dependent(s0, s1);
+  bool id02 = !dependent(s0, s2);
+  if(id01 && id02) {
+    int di01 = sd(s0, s1);
+    int di02 = sd(s0, s2);
+    return outNZ(brec_fix_main(di01+pos, s1, di02+pos, s2));
+    }
+  else {
+    bool id12 = (id01 || id02) && !dependent(s1, s2);
+    if(id12 && id02) {
+      //01 not
+      int di12 = sd(s1, s2);
+      int di02 = sd(s0, s2);
+      return outNZ(tallybox_total(s2->qty) * brec_fix_other(di02, lst, pos, di12, s1));
+      }
+    else if(id12 && id01) {
+      int di12 = sd(s1, s2);
+      int di01 = sd(s0, s1);
+      return outNZ(tallybox_total(s1->qty) * brec_fix_other(di01, lst, pos, di12, s2));
+      }
+    else if(pos == 0)
+      return outNZ(brec_fix_main(sd(s0, s1), s1, sd(s0, s2), s2));
+    else {
+      betweenness_type total = 0;
+      if(get0(s1->qty))
+        total += get0(s1->qty) * brec_fix_other(sd(s0, s1), lst, pos, sd(s1, s2), s2);
+      if(get0(s2->qty))
+        total += get0(s2->qty) * brec_fix_other(sd(s0, s2), lst, pos, sd(s1, s2), s1);
+      if(get0(s1->qty) && get0(s2->qty))
+        total -= get0(s1->qty) * get0(s2->qty) * ack(pos + sd(s0, s1), pos + sd(s0, s2), sd(s1, s2));
+      for(segment *c1 = s1->firstchild; c1; c1 = c1->nextchild)
+      for(segment *c2 = s2->firstchild; c2; c2 = c2->nextchild)
+        total += brec(lst, pos-1, c1, c2);
+      return outNZ(total);
+      }
+    }
+  }
+
+betweenness_type betweenness4(mycell *c, int setid = 0) {
+  segment *s = getsegment(c, c, setid, true);
+  vector<segment*> lst;
+  while(s) { lst.push_back(s); s = s->parent; }
+  segment *rs = getsegment(mroot, mroot, setid, true);
+  return brec(lst, isize(lst)-1, rs, rs);
+  }
+
+bool neq(betweenness_type a, betweenness_type b) {
+  return abs(a-b) > 1e-6;
+  }
+
+void compute_betweenness(bool verify) {
+  progressbar pb(N, "compute_betweenness");
+  int errorcount = 0, errorcount2 = 0;
+  for(int i=0; i<N; i++) {
+    long long total_cache = 0, x_total_cache = 0;
+    int clean_level = -1;
+    for(auto& mm: memo) {
+      ll s = mm.size();
+      if(!s) continue;
+      total_cache += s;
+      if(total_cache >= 8000000000ll / 64) { mm.clear(); }
+      else { x_total_cache += s; clean_level++; }
+      }
+    if(total_cache != x_total_cache) printf("cleanup from %lld to %lld cache items, at level %d\n", total_cache, x_total_cache, clean_level);
+    // if(i != 1) continue;
+    mycell *c1 = vertices[i];
+    // add_to_set(c1, -1, 0);
+    auto b = betweenness3(c1);
+    // add_to_set(c1, 1, 0);
+    auto b4 = betweenness4(c1);
+    print(hlog, format("B;%10Ld;%10Ld;%20.10Lf;%3d;%-40s", b.first, b.second, b4, vertices[i]->lev, rogueviz::vdata[i].name.c_str()));
+    if(verify) {
+      /*
+      betweenness_type a = b.first;
+      betweenness_type b = 0;
+      for(int j=0; j<N; j++)
+      for(int k=0; k<N; k++)
+        b += ack(quickdist(vertices[j], c1), quickdist(vertices[k], c1), quickdist(vertices[j], vertices[k]));
+        if(quickdist(vertices[j], vertices[k]) == quickdist(vertices[j], c1) + quickdist(vertices[k], c1)) {
+          // iprintf("  %-40s %-40s %d/%d/%d\n", rogueviz::vdata[j].name.c_str(), rogueviz::vdata[k].name.c_str(), quickdist(vertices[j], vertices[k]), quickdist(vertices[j], c1), quickdist(vertices[k], c1));
+          b++;
+          }
+      if(neq(a, b)) iprintf(" ERROR: %f\x1b[K\n", double(b));
+      else iprintf(" CORRECT %f\x1b[K\n", double(b));
+      if(b-a>0) errorcount += (b-a);
+      if(b-a<0) errorcount2 += (a-b);
+      if(neq(a, b)) exit(1);
+      */
+      }
+    else printf("\n");
+    pb++;
+    }
+  if(verify) println(hlog, format("errorcount = %d/%d\n", errorcount, errorcount2));
+  }
+
+void build(mycell *c, int lev, string s) {
+  int id = 0;
+  if(lev) for(mycell *c1: allchildren(c)) { build(c1, lev-1, s + std::to_string(id)); id++; }
+  printf("* %s\n", s.c_str());
+  vertices.push_back(c);
+  rogueviz::vdata.emplace_back();  
+  rogueviz::vdata.back().name = s;
+  /*
+  vertices.push_back(c);
+  rogueviz::vdata.emplace_back();  
+  rogueviz::vdata.back().name = s;
+  */
+  }
+
+void build_all(int d) {
+  build(mroot, d, "");
+  N = isize(vertices);
+  counttallies();
+  }
+
+void load_test() {
+  string s;
+  while(getline(std::cin, s)) {
+    mycell *mc = mroot;
+    if(s[0] == '#') continue;
+    for(char c: s) if(c >= '0' && c <= '9') mc = allchildren(mc) [c - '0'];
+    vertices.push_back(mc);
+    rogueviz::vdata.emplace_back();
+    rogueviz::vdata.back().name = "PATH:" + s;
+    }
+  // build(mroot, 5, "");
+  N = isize(vertices);
+  counttallies();
+  // add_to_set(vertices[0], -1, 0);
+  }
+
+}
diff --git a/rogueviz/dhrg/dhrg-128.cpp b/rogueviz/dhrg/dhrg-128.cpp
new file mode 100644
index 00000000..50bdaf6c
--- /dev/null
+++ b/rogueviz/dhrg/dhrg-128.cpp
@@ -0,0 +1,2 @@
+#define BOXSIZE 128
+#include "dhrg.cpp"
diff --git a/rogueviz/dhrg/dhrg-64.cpp b/rogueviz/dhrg/dhrg-64.cpp
new file mode 100644
index 00000000..f1fb5fc1
--- /dev/null
+++ b/rogueviz/dhrg/dhrg-64.cpp
@@ -0,0 +1,2 @@
+#define BOXSIZE 64
+#include "dhrg.cpp"
diff --git a/rogueviz/dhrg/dhrg-utils.cpp b/rogueviz/dhrg/dhrg-utils.cpp
new file mode 100644
index 00000000..3b996de2
--- /dev/null
+++ b/rogueviz/dhrg/dhrg-utils.cpp
@@ -0,0 +1,33 @@
+// general utilities
+#include <sys/time.h>
+
+namespace dhrg {
+
+typedef long long ll;
+
+struct progressbar : indenter_finish {
+  string name;
+  static const int PBSIZE = 64;
+  int step = -1, total, drawat = 0, count = -1;
+  
+  void operator ++ (int) { 
+    step++;
+    while(step >= drawat && total) {
+      count++;
+      drawat = (total * (count+1)) / PBSIZE;
+      fprintf(stderr, "%s [", get_stamp().c_str());
+      for(int k=0; k<count; k++) fprintf(stderr, "#");
+      for(int k=count; k<64; k++) fprintf(stderr, "-");
+      fprintf(stderr, "] %s\x1b[K\r", name.c_str());
+      fflush(stderr);
+      }
+    }
+  
+  ~progressbar() {
+    fprintf(stderr, "\x1b[K");
+    }
+
+  progressbar(int t, string n) : indenter_finish(n) { hlog.indentation -= 2; println(hlog, name); hlog.indentation += 2; total = t; (*this)++; }
+  };
+
+}
diff --git a/rogueviz/dhrg/dhrg.cpp b/rogueviz/dhrg/dhrg.cpp
new file mode 100644
index 00000000..65651608
--- /dev/null
+++ b/rogueviz/dhrg/dhrg.cpp
@@ -0,0 +1,276 @@
+// see this paper: https://arxiv.org/abs/2109.11772
+
+#define DHRGVER "7.1"
+#include "../rogueviz.h"
+
+#define LONG_BRACKETS
+
+namespace rogueviz { extern string fname; }
+
+namespace dhrg {
+
+using namespace hr;
+
+#ifndef BOXSIZE
+static const int BOXSIZE = 32;
+#endif
+static const int MAXDIST = (2*BOXSIZE);
+static const int SETS = 4;
+
+struct segment;
+
+cell *croot() { return currentmap->gamestart(); }
+
+int M;
+vector<struct mycell*> vertices;
+
+vector<ld> disttable0, disttable1;
+
+void memoryInfo();
+
+void cellcoords();
+void origcoords();
+void build_disttable();
+
+void dhrg_init();
+bool dhrg_animate(int sym, int uni);
+}
+
+#include "readgraph.cpp"
+#include "dhrg-utils.cpp"
+#include "regular.cpp"
+#include "gridmapping.cpp"
+#include "mycell.cpp"
+#include "segment.cpp"
+#include "dynamic.cpp"
+#include "loglik.cpp"
+#include "paths.cpp"
+#include "embedder.cpp"
+#include "tests.cpp"
+#include "betweenness.cpp"
+#include "groundtruth.cpp"
+#include "routing.cpp"
+
+namespace dhrg {
+
+void memoryInfo() {
+  string s = "";
+  
+  ll totalmemory = 0;
+
+  for(int it=0; it<2; it++) {
+    auto pct = [&] (auto x, auto y) { 
+      if(it == 0) totalmemory += x*y;
+      else s += " " + its(x) + "x" + its(y) + "B=" + its((x*y*100)/totalmemory) + "%";
+      };
+    pct(cellcount, sizeof(cell));
+    pct(heptacount, sizeof(heptagon));
+    pct(mycellcount, sizeof(mycell));
+    pct(segmentcount, sizeof(segment));
+    }
+  println(hlog, "Memory info: ", s, " (", int(totalmemory/1048576), " MB)");
+  fflush(stdout);
+  }
+
+void debugtally() {
+  print(hlog, "T"); for(int i=0; i<MAXDIST; i++) if(tally[i]) print(hlog, format(" %4d/%4lld", edgetally[i], tally[i]));
+  println(hlog, " .. %", loglikopt()); 
+  }
+
+void place_rogueviz_vertices() {
+  progressbar pb(N, "place_rogueviz_vertices");
+  // transmatrix In = inverse(ggmatrix(currentmap->gamestart()));
+  using namespace rogueviz;
+/*  for(int i=0; i<N; i++) vdata[i].m->base = currentmap->gamestart();
+  for(int i=0; i<N; i++) vdata[i].m->at = In * shmup::ggmatrix(vertices[i]->ascell()); */
+  for(int i=0; i<N; i++) vdata[i].m->base = vertices[i]->ascell();
+  for(int i=0; i<N; i++) vdata[i].m->at = Id;
+  fixedges(); 
+  for(int i=N; i<isize(vdata); i++)
+    if(vdata[i].m)
+      vdata[i].m->store();
+  }
+
+int destroy;
+
+void clear() {
+  coords.clear();  
+  // destroytallies();
+  //  todo: correct cleanup!
+#ifdef BUILD_ON_HR
+  mymap.clear();
+#else
+  delete mroot;
+  mroot = NULL;
+#endif
+  }
+  
+// actual octiles of normal are roughly: -1.15 -.674 -.318 0 .318 .674 1.15
+
+void dhrg_init() {
+  if(!mroot) {
+    println(hlog, "DHRG version " DHRGVER "\n");
+    rogueviz::init(0);
+    rogueviz::rv_hook(hooks_handleKey, 100, dhrg_animate);
+    regular_info();
+    generate_root();
+    }
+  }
+
+int next_timestamp;
+int ts_rogueviz;
+int ts_rbase;
+int ts_coords;
+int ts_vertices;
+bool stored;
+
+int dhrgArgs() {
+  using namespace arg;
+           
+  if(argis("-dhrg")) {
+    PHASE(3); shift(); dhrg_init(); read_graph_full(args());
+    next_timestamp++;
+    ts_rogueviz = next_timestamp;
+    ts_vertices = next_timestamp;
+    }
+    
+  else if(argis("-graph")) {
+    PHASE(3); shift(); dhrg_init(); read_graph(args(), false, false, false);
+    next_timestamp++;
+    ts_rogueviz = next_timestamp;
+    // stored = true;
+    }
+
+  else if(argis("-graphv")) {
+    PHASE(3); shift(); dhrg_init(); read_graph(args(), true, true, true);
+    next_timestamp++;
+    ts_rogueviz = next_timestamp;
+    ts_rbase = next_timestamp;
+    stored = true;
+    }
+  
+  else if(argis("-analyze_grid")) {
+    PHASE(3); shift(); dhrg_init(); do_analyze_grid(argi());
+    }
+
+  else if(argis("-analyze_dists")) {
+    PHASE(3); dhrg_init(); shift(); do_analyze_dists(argi());
+    }
+
+  else if(argis("-test_paths")) {
+    PHASE(3); dhrg_init(); shift(); test_paths(argi());
+    }
+
+  else if(argis("-contll")) {
+    if(ts_rogueviz >= ts_rbase && ts_rogueviz >= ts_vertices && ts_rogueviz > ts_coords) {
+      origcoords();
+      build_disttable_approx();
+      ts_coords = ts_rogueviz;
+      }
+    else if(ts_rbase >= ts_vertices && ts_rbase > ts_coords) {
+      ts_coords = ts_rbase;
+      rvcoords();
+      build_disttable_approx();
+      }
+    else if(ts_vertices > ts_coords) {
+      ts_coords = ts_vertices;
+      cellcoords();
+      build_disttable_approx();
+      }
+    logistic cont;
+    cont.setRT(graph_R, graph_T);
+    fast_loglik_cont(cont, loglik_cont_approx, "lcont", 1, 1e-6);
+    // return loglik_cont();
+    }
+
+  else if(argis("-dhrgview")) {
+    if(ts_vertices > ts_rbase) {
+      place_rogueviz_vertices();
+      ts_rbase = ts_vertices;
+      }
+    if(!stored) rogueviz::storeall(), stored = true;
+    else shmup::fixStorage();
+    mainloop(); quitmainloop = false;
+    }
+
+  else if(argis("-iterate")) {
+  
+    if(ts_rbase > ts_vertices) {
+      PHASE(3); dhrg_init(); graph_from_rv(); 
+      ts_vertices = ts_rbase;
+      }
+    if(!ts_vertices) {
+      printf("Error: read vertices with -dhrg or -graph\n");
+      exit(1);
+      }
+    shift();
+    embedder_loop(argi());
+    next_timestamp++;
+    ts_vertices = next_timestamp;
+    }
+
+  else if(argis("-dorestart")) {
+    dorestart = true;
+    }
+
+  else if(argis("-dontrestart")) {
+    dorestart = false;
+    }
+
+  else if(argis("-lctype")) {
+    shift(); lc_type = args()[0];
+    }
+
+  else if(argis("-loadtest")) {
+    dhrg_init(); load_test();
+    }
+
+  else if(argis("-buildtest")) {
+    shift(); dhrg_init(); build_all(argi());
+    }
+
+  else if(argis("-pbv")) {
+    compute_betweenness(true);
+    }
+
+  else if(argis("-pb")) {
+    compute_betweenness(false);
+    }
+
+  else if(argis("-gtt")) {
+    shift(); ground_truth_test(args());
+    }
+
+  else if(argis("-routing")) {
+    shift(); routing_test(args());
+    }
+
+  else if(argis("-esaveas")) {
+    shift(); save_embedding(args());
+    }
+
+  else if(argis("-esave")) {
+    save_embedding(rogueviz::fname + "-dhrg.txt");
+    }
+
+  else if(argis("-eload")) {
+    PHASE(3); shift(); dhrg_init(); load_embedded(args());
+    next_timestamp++;
+    ts_rogueviz = next_timestamp;
+    ts_vertices = next_timestamp;
+    }
+
+  else return 1;
+
+  return 0;
+  }
+
+auto hook = 
+    addHook(hooks_args, 50, dhrgArgs)
+  + addHook(hooks_clearmemory, 200, clear);
+
+#if CAP_SDL
+#include "visualize.cpp"
+#endif
+
+} 
diff --git a/rogueviz/dhrg/dynamic.cpp b/rogueviz/dhrg/dynamic.cpp
new file mode 100644
index 00000000..25998004
--- /dev/null
+++ b/rogueviz/dhrg/dynamic.cpp
@@ -0,0 +1,146 @@
+// computing pairs of vertices in each distance using dynamic programming (as described in the paper)
+
+namespace dhrg {
+
+bool segmentValid(mycell *cl, mycell *cr) {
+  if(cl == cr) return true;
+  mycell *c1 = cl; 
+  int d = 0;
+  for(; d<7 && c1 != cr; d++) {
+    c1->build(); c1 = c1->grightsibling();
+    }
+  if(d == 7) return false;
+  cr->build();
+  cl = cl->grightsibling()->gleftchild();
+  cr = cr->gleftchild();
+  return segmentValid(cl, cr);
+  }
+
+vector<segment*> all_child_segments(segment *s) {
+  vector<segment*> res;
+  for(auto m1: allchildren(s->left, -1))
+  for(auto m2: allchildren(s->right, +1))
+    if(segmentValid(m1, m2))
+      res.push_back(getsegment(m1,m2,0,true));
+  return res;
+  }
+
+// returns 0 if not in segment, 1-based index if in segment
+int insegment(mycell *mc, segment *s1) {
+  mycell *l = s1->left;
+  int i = 1;
+  while(true) {
+    if(l == mc) return i;
+    if(l == s1->right) return 0;
+    l->build(); l = l->grightsibling();
+    i++;
+    }
+  }
+
+int segmentcode(segment *s) {
+  mycell *l = s->left;
+  int code = 0;
+  while(true) {
+    code += l->gettype();
+    if(l == s->right) return code;
+    l->build(); l = l->grightsibling();
+    code *= 8;
+    }
+  }
+
+int compute_descendants(segment *s, int d) {
+
+  auto memokey = make_tuple(segmentcode(s), d);
+  static map<decltype(memokey), int> mem;
+  if(mem.count(memokey)) return mem[memokey];
+
+  if(d == 0) return s->left == s->right ? 1 : 0;
+  int total = 0;
+  for(auto s1: all_child_segments(s))
+    total += compute_descendants(s1, d-1);
+  return mem[memokey] = total;
+  }
+
+// returns 0 if segments are not crossing, positive number if crossing
+// for segments with equal codes, equal numbers = the same way of crossing
+int segmentcross(segment *s1, segment *s2) {
+  int i1;
+  i1 = insegment(s1->left, s2);
+  if(i1) return 4*i1+1;
+  i1 = insegment(s1->right, s2);
+  if(i1) return 4*i1+2;
+  i1 = insegment(s2->left, s1);
+  if(i1) return 4*i1+3;
+  i1 = insegment(s2->right, s1);
+  if(i1) return 4*i1+4;
+  return 0;
+  }
+
+set<int> allsegments;
+
+int compute_in_dist(segment *s1, segment *s2, int d1, int d2, int dex) {
+  // if(d1 + d2 + 4 < dex) return 0;
+  int d = -segmentcross(s1,s2);
+  if(!d) d = segmentdist(s1, s2, 0);
+
+  if(d > 2 || d1 == 0 || d2 == 0) {
+    if(d < 0) d = 0;
+    return d1+d2+d == dex ? compute_descendants(s1,d1) * compute_descendants(s2,d2) : 0;
+    }
+  else {
+    mycell *ss1 = s1->right, *ss2 = s2->right;
+    if(d > 0) {
+      int side = 0;
+      while(true) {
+        ss1->build(); ss1 = ss1->grightsibling(); if(ss1 == s2->left) { side=1; break; }
+        ss2->build(); ss2 = ss2->grightsibling(); if(ss2 == s1->left) { side=2; break; }
+        }
+      d += 100 * side;
+      }
+    allsegments.insert(segmentcode(s1));
+    allsegments.insert(segmentcode(s2));
+    auto memokey = make_tuple(segmentcode(s1), segmentcode(s2), d, d1, d2, dex);
+    static map<decltype(memokey), int> mem;
+    if(mem.count(memokey)) return mem[memokey];
+    int total = 0;
+    for(auto s3: all_child_segments(s1))
+    for(auto s4: all_child_segments(s2))
+      total += compute_in_dist(s3, s4, d1-1, d2-1, dex);
+    if(0) if(mem.count(memokey) && mem[memokey] != total) {
+      printf("%d vs %d :: %x %x d=%d %d,%d,%d\n", mem[memokey], total, segmentcode(s1), segmentcode(s2), d, d1, d2, dex);
+      return mem[memokey];
+      }
+    return mem[memokey] = total;
+    }
+  }
+
+void do_analyze_dists(int rad) {
+
+  println(hlog, "do_analyze_dists (", rad, ")");
+  indenter_finish indf;
+  auto m = mroot;
+  
+  auto seg = getsegment(m, m, 0, true);
+
+  // compute the correct answer, but not if this requires creating more than 1500 cells
+  celllister cl(croot(), rad, 1500, NULL);
+  vector<int> correct(2*rad+4, 0);
+  
+  for(cell *c1: cl.lst) if(celldist(c1) == rad)
+  for(cell *c2: cl.lst) if(celldist(c2) == rad)
+    correct[celldistance(c1,c2)]++;
+    
+  int total = 0;
+  for(int a=0; a<2*rad+4; a++) {
+    int cd = compute_in_dist(seg, seg, rad, rad, a);
+    printf("%2d: %d/%d\n", a, cd, correct[a]);
+    total += cd;
+    }
+  
+  printf("total = %d (%d)\n", total, cgi.expansion->get_descendants(5).approx_int() * cgi.expansion->get_descendants(5).approx_int());
+  printf("all segments = %d\n", isize(allsegments));
+  
+  // printf("descendants = %d (%d)\n", compute_descendants(seg, 5), int(.1+expansion.get_descendants(5).approx_int()));
+  }
+
+}
diff --git a/rogueviz/dhrg/embedder.cpp b/rogueviz/dhrg/embedder.cpp
new file mode 100644
index 00000000..47b0b155
--- /dev/null
+++ b/rogueviz/dhrg/embedder.cpp
@@ -0,0 +1,340 @@
+// local search
+
+namespace dhrg {
+
+int newmoves = 0;
+int iterations = 0;
+int distlimit;
+
+void dispnewmoves() {
+  if(newmoves == 0) printf(".");
+  else if(newmoves < 10) printf("%c", '0'+newmoves);
+  else printf("%c", 'a' + newmoves/10);
+  newmoves = 0;
+  }
+
+vector<char> tomove;
+
+bool smartmove = false;
+bool dorestart = false;
+
+int lastmoves;
+
+int movearound() {
+  indenter_finish im("movearound");
+  int total = 0;
+  if(smartmove) for(bool b: tomove) if(b) total++;
+  if(total == 0) {
+    tomove.resize(0), tomove.resize(N, true);
+    }
+  int moves = 0;
+  ld llo = loglik_chosen();
+//  int im = 0;
+  
+  {progressbar pb(N, "tomove: " + its(total) + " (last: " + its(lastmoves) + ")");
+  for(int i=0; i<N; i++) {
+    if(!tomove[i]) { pb++; continue; }
+    tomove[i] = false;
+    // if(i && i % 100 == 0) dispnewmoves();
+    mycell *mc = vertices[i];
+    tallyedgesof(i, -1, mc);
+    add_to_set(mc, -1, 0);
+    add_to_tally(mc, -1, 0);
+//    im += size(rogueviz::vdata[i].edges);
+    mycell *mc2[FULL_EDGE];
+    ld llo2[FULL_EDGE];
+    int bestd = -1;
+    ld bestllo = llo;
+    
+    auto nei = allneighbors(mc);
+    
+    for(int d=0; d<isize(nei); d++) {
+      mc2[d] = nei[d];
+      if(mc2[d]->lev >= distlimit) continue;
+      tallyedgesof(i, 1, mc2[d]);
+      add_to_tally(mc2[d], 1, 0);
+      if(lc_type == 'C')
+        add_to_set(mc2[d], 1, 0);
+      llo2[d] = loglik_chosen();
+      if(lc_type == 'C')
+        add_to_set(mc2[d], -1, 0);
+      if(llo2[d] > bestllo) bestd = d, bestllo = llo2[d];
+      
+      add_to_tally(mc2[d], -1, 0);
+      tallyedgesof(i, -1, mc2[d]);
+      }
+    if(bestd >= 0) {
+      moves++; newmoves++;
+      vertices[i] = mc = mc2[bestd];
+      llo = llo2[bestd];
+      tomove[i] = true;
+      for(auto p: rogueviz::vdata[i].edges) {
+        int j = p.second->i ^ p.second->j ^ i;
+        tomove[j] = true;
+        }
+      }
+    tallyedgesof(i, 1, mc);
+    add_to_tally(mc, 1, 0);
+    add_to_set(mc, 1, 0);
+    pb++;
+    }}
+  // dispnewmoves();
+  println(hlog, " moves = ", moves);
+  return lastmoves = moves;
+  }
+
+int move_restart() {
+  indenter_finish im("move_restart");
+  ld llo = loglik_chosen();
+  array<array<int, 128>, 2> distances_map = {0};
+  int moves = 0;
+//  int im = 0;
+  
+  {progressbar pb(N, "move_restart");
+  for(int i=0; i<N; i++) {
+   // if(i && i % 100 == 0) dispnewmoves();
+    mycell *mc = vertices[i];
+    tallyedgesof(i, -1, mc);
+    add_to_set(mc, -1, 0);
+    add_to_tally(mc, -1, 0);
+    
+    auto getllo = [&] (mycell *m) {
+      tallyedgesof(i, 1, m);
+      add_to_tally(m, 1, 0);
+      if(lc_type == 'C')
+        add_to_set(m, 1, 0);
+      ld res = loglik_chosen();
+      if(lc_type == 'C')
+        add_to_set(m, -1, 0);
+      add_to_tally(m, -1, 0);
+      tallyedgesof(i, -1, m);
+      return res;
+      };
+
+    mycell* whereto = mroot;
+    ld bestllo = getllo(whereto);
+    bool changed = true;
+    while(changed) {
+      changed = false;
+      mycell *mc2 = whereto;
+      auto nei2 = allneighbors(whereto);
+      for(mycell *mcn: nei2) {
+        ld newllo = getllo(mcn);
+        if(newllo > bestllo)
+          bestllo = newllo, mc2 = mcn, changed = true;
+        }
+      if(changed) whereto = mc2;
+      }
+    bool better = bestllo > llo;
+    distances_map[better][quickdist(whereto, mc)]++;
+    if(better) {
+      llo = bestllo;
+      vertices[i] = mc = whereto;
+      moves++;
+      }
+    tallyedgesof(i, 1, mc);
+    add_to_tally(mc, 1, 0);
+    add_to_set(mc, 1, 0);
+    pb++;
+    }}
+  // dispnewmoves();
+  println(hlog, " moves = ", moves);
+  print(hlog, " stats:");
+  for(int a=0; a<2; a++) for(int b=0; b<128; b++) {
+    int d = distances_map[a][b];
+    if(d) print(hlog, format(" %d/%d:%d", a,b, d));
+    }
+  println(hlog, "\n");
+  return lastmoves = moves;
+  }
+
+// 7: 12694350 3847975716
+// 6: 39472959 11969080911
+
+void verifycs() {
+  long long edgecs = 0, totalcs = 0;
+  for(int u=0; u<MAXDIST; u++) 
+    edgecs += edgetally[u] * u*u,
+    totalcs += tally[u] * u*u;
+  
+  print(hlog, "edgecs=", format("%lld", edgecs), " totalcs=", format("%lld", totalcs));
+  }
+
+void preparegraph() {
+  indenter_finish im("preparegraph");
+  using namespace rogueviz;
+  M = 0;
+  vertices.resize(N);
+
+  if(1) {
+    progressbar pb(N, "Finding all");
+    for(int i=0; i<N; i++) {
+      M += isize(vdata[i].edges); 
+      pb++;
+      }
+    M /= 2;
+    }
+
+  memoryInfo();
+  
+  cont_logistic.setRT(graph_R, graph_T);
+
+  counttallies();
+  
+  memoryInfo();
+  
+  verifycs();
+
+  if(1) {
+    indenter_finish im("optimizing parameters"); 
+    ld factor = 1/log(cgi.expansion->get_growth());
+    current_logistic.setRT(factor * graph_R, factor * graph_T);
+    saved_logistic = current_logistic; 
+
+    // for(int u=0; u<MAXDIST; u++) iprintf("%d/%Ld\n", edgetally[u], tally[u]);
+
+    fix_logistic_parameters(current_logistic, loglik_logistic, "logistic", 1e-6);
+    writestats();
+    fflush(stdout);
+    }
+  
+  iterations = 0;
+
+  distlimit = 0;
+  for(int j=0; j<BOXSIZE; j++) 
+    if(getsegment(mroot,mroot,0)->qty[j])
+      distlimit = min(j+2, BOXSIZE-1);
+
+  println(hlog, "Using distlimit = ", distlimit);
+  }
+
+void read_graph_full(const string& fname) {
+  using namespace rogueviz;
+  
+  memoryInfo();
+
+  if(true) {
+    indenter_finish im("Read graph");
+  
+    // N = isize(vdata);
+  
+    read_graph(fname, false, false, false);
+    vertices.resize(N);
+    progressbar pb(N, "Translating to cells");
+
+    for(int i=0; i<N; i++) {
+#if BUILD_ON_HR    
+      cell *c = currentmap->gamestart();
+      hyperpoint T0 = vdata[i].m->at * C0;
+      virtualRebase2(vdata[i].m->base, T0, true); 
+      vertices[i] = find_mycell(vdata[i].m->base);
+#else
+      vertices[i] = find_mycell_by_path(computePath(vdata[i].m->at));
+#endif
+      vdata[i].m->at = Id;
+      pb++;
+      // printf("%s\n", computePath(vdata[i].m->base).c_str());
+      }    
+    }
+
+  recycle_compute_map();
+  preparegraph();  
+  }
+
+void graph_from_rv() {
+  using namespace rogueviz;
+  
+  memoryInfo();
+  
+  vertices.resize(N);
+  progressbar pb(N, "converting RogueViz to DHRG");
+
+  for(int i=0; i<N; i++) {
+#if BUILD_ON_HR    
+    vertices[i] = find_mycell(vdata[i].m->base);
+#else
+    auto path1 = computePath(vdata[i].m->base);
+    vertices[i] = find_mycell_by_path(path1);
+#endif
+    vdata[i].m->at = Id;
+    pb++;
+    }    
+
+  preparegraph();
+  }
+
+bool iteration() {
+  iterations++;
+  indenter_finish im("Iteration #" + its(iterations));
+  int m = movearound();
+  if(!m && dorestart) m = move_restart();
+  if(!m) return false;
+  fix_logistic_parameters(current_logistic, loglik_logistic, "logistic", 1e-6);
+  writestats();
+  fflush(stdout);
+  return true;
+  }
+
+void embedder_loop(int max) {
+  indenter_finish im("embedder_loop");
+  ld lastopt = loglik_chosen();
+  
+  while(max-- && iteration()) {
+    ld curopt = loglik_chosen();
+    println(hlog, "current = %", curopt);
+    if(curopt <= lastopt) {
+      println(hlog, "Not enough improvement -- breaking"); 
+      break;
+      }
+    lastopt = curopt;
+    }
+  }
+
+void save_embedding(const string s) {
+  FILE *f = fopen(s.c_str(), "wt");
+  for(int i=0; i<N; i++) {
+    string p = get_path(vertices[i]);
+    if(p == "") p = "X";
+    fprintf(f, "%s %s\n", rogueviz::vdata[i].name.c_str(), p.c_str());
+    }
+  fclose(f);
+  }
+
+void load_embedded(const string s) {
+  if(true) {
+    read_graph(s, false, false, false);
+    indenter_finish im("Read graph");
+    }
+
+  string t = rogueviz::fname + "-dhrg.txt";
+      
+  if(true) {
+    progressbar pb(N, "reading embedding");
+    vertices.resize(N, NULL);
+    
+    map<string, int> ids;
+    for(int i=0; i<N; i++) ids[rogueviz::vdata[i].name] = i;
+
+    FILE *f = fopen(t.c_str(), "rt");
+    while(true) {
+      char who[500], where[500];
+      who[0] = 0;
+      fscanf(f, "%s%s", who, where);
+      if(who[0] == 0) break;
+      if(!ids.count(who)) printf("unknown vertex: %s\n", who);
+      string wh = where;
+      if(wh == "X") wh = "";
+      vertices[ids[who]] = find_mycell_by_path(wh);
+      pb++;
+      }
+    fclose(f);
+    
+    for(int i=0; i<N; i++) if(vertices[i] == NULL) {
+      printf("unmapped: %s\n", rogueviz::vdata[i].name.c_str());
+      exit(1);
+      }
+    }
+  preparegraph();
+  }
+
+}
diff --git a/rogueviz/dhrg/gridmapping.cpp b/rogueviz/dhrg/gridmapping.cpp
new file mode 100644
index 00000000..b4407ea4
--- /dev/null
+++ b/rogueviz/dhrg/gridmapping.cpp
@@ -0,0 +1,82 @@
+// test the conjecture comparing triangulation distances and hyperbolic distances
+
+namespace dhrg {
+
+ld stats[32][3], wstats[32][3];
+  
+hyperpoint celltopoint(cell *c) {
+  return tC0(calc_relative_matrix(c, croot(), C0));
+  }
+
+void do_analyze_grid(int maxv) {
+  cell *root = croot();
+  celllister cl(root, 32, maxv, NULL);
+  // if this works too slow, use a smaller number
+  // (you can also use a larger number if you have time of course)
+
+  // int rot = 0;
+  
+  vector<ld> distances[128];
+  
+  for(cell *c: cl.lst) {
+    hyperpoint h = celltopoint(c);
+    ld dd = hdist0(h);
+    int d = celldist(c);
+    stats[d][0] ++;
+    stats[d][1] += dd;
+    stats[d][2] += dd*dd;
+    distances[d].push_back(dd);
+
+    if(d>0) {
+      ld alpha[2];
+      int qalpha = 0;
+      forCellCM(c2, c) if(celldist(c2) == d) {
+        hyperpoint h1 = celltopoint(c2);
+        alpha[qalpha++] = atan2(h1[0], h1[1]);
+        }
+      if(qalpha != 2) printf("Error: qalpha = %d\n", qalpha);
+      ld df = alpha[0] - alpha[1];
+      if(df<0) df = -df;
+      while(df > 2*M_PI) df -= 2*M_PI;
+      while(df > M_PI) df = 2*M_PI - df;
+      df /= 4*M_PI;
+      wstats[d][0] += df;
+      if(d==2) printf("df == %" PLDF " dd = %" PLDF "\n", df, dd);
+      wstats[d][1] += df*dd;
+      wstats[d][2] += df*dd*dd;
+      }
+    }
+  
+  println(hlog, "log(gamma) = ", log(cgi.expansion->get_growth()));
+  
+  ld lE, dif, lwE;
+  for(int d=0; d<32; d++) if(stats[d][0]) {
+    int q = stats[d][0];
+    if(q != cgi.expansion->get_descendants(d).approx_int()) continue; 
+    ld E = stats[d][1] / q;
+    ld E2 = stats[d][2] / q;
+    ld Vr = E2 - E * E;
+    if(Vr < 0) Vr = 0;
+    dif = E- lE; lE = E;
+    ld Vd = d > 1 ? Vr/(d-1) : 0;
+    
+    ld wE = wstats[d][1];
+    ld wE2 = wstats[d][2];
+    ld wVr = wE2 - wE * wE;
+
+    print(hlog, format("d=%2d: q = %8d E = %12.8" PLDF " dif = %12.8" PLDF " Vr = %12.8" PLDF " Vr/(d-1)=%12.8" PLDF,
+      d, q, E, dif, Vr, Vd));
+
+    if(0) print(hlog, format(" | <%" PLDF "> ex = %12.8" PLDF " d.ex = %12.8" PLDF " Vr = %12.8" PLDF, wstats[d][0], wE, wE - lwE, wVr));
+    
+    ld Sigma = sqrt(Vr);
+    sort(distances[d].begin(), distances[d].end());
+    if(Sigma) for(int u=1; u<8; u++)
+      print(hlog, format(" %8.5" PLDF, (distances[d][u * isize(distances[d]) / 8] - E) / Sigma));
+    
+    println(hlog);
+    lwE = wE;
+    }
+  }
+
+}
diff --git a/rogueviz/dhrg/groundtruth.cpp b/rogueviz/dhrg/groundtruth.cpp
new file mode 100644
index 00000000..01a7e40e
--- /dev/null
+++ b/rogueviz/dhrg/groundtruth.cpp
@@ -0,0 +1,121 @@
+
+//	$(VARIANT) -nogui -dhrg embedded-graphs/facebook_combined_result -contll -iterate 99 -contll -esave > $@
+
+namespace dhrg {
+
+bool is(char *where, const char *what) {
+  while(*where && *what && *where == *what) where++, what++;
+  return !*what;
+  }
+
+void ground_truth_test(string s) {
+
+  logistic cont;
+
+  vector<ld> logliks;
+  vector<string> reps;
+
+  auto report = [&] (string s, ld val) {
+    logliks.push_back(val);
+    reps.push_back(s);
+    println(hlog, "REPORT ", s, " = ", val);
+    };
+
+  if(1) {
+    FILE *f = fopen(("embout/" + s).c_str(), "rb");
+    char buf[999999];
+    int siz = fread(buf, 1, 999999, f);
+    
+    for(int i=0; i<siz; i++) if(is(buf+i, "Embedded Log-likelihood: "))
+      report("bfkl", atof(buf+i+25));
+
+    fclose(f);
+    }
+if(1) {
+  dhrg_init(); read_graph_full("data/simg-" + s);
+
+
+  origcoords();
+  build_disttable();
+
+  cont.setRT(graph_R, graph_T);
+
+  report("gz", loglik_cont(cont));
+  fix_logistic_parameters(cont, loglik_cont, "lcont", 1e-3);
+  report("grc", loglik_cont(cont));
+
+  cellcoords();
+  build_disttable();
+
+  report("gco", loglikopt());
+  report("gcm", loglikopt_mono());
+  report("gcrt", loglik_logistic());
+
+  cont.setRT(graph_R, graph_T);
+  report("gcz", loglik_cont(cont));
+  fix_logistic_parameters(cont, loglik_cont, "lcont", 1e-3);
+  report("gcrc", loglik_cont(cont));
+
+  embedder_loop(20);
+
+  report("geo", loglikopt());
+  report("gem", loglikopt_mono());
+  report("gert", loglik_logistic());
+
+  cellcoords();
+  build_disttable();
+
+  cont.setRT(graph_R, graph_T);
+  fix_logistic_parameters(cont, loglik_cont, "lcont", 1e-3);  
+  report("gerc", loglik_cont(cont));
+
+  delete mroot;
+  mroot = NULL;
+  segmentcount = 0;
+  for(int i=0; i<MAXDIST; i++) tally[i] = 0;
+  for(int i=0; i<MAXDIST; i++) edgetally[i] = 0;
+  vertices.clear();
+  rogueviz::close();
+}
+
+  dhrg_init(); read_graph_full("data/sime-" + s);
+  origcoords();
+  build_disttable();
+  cont.setRT(graph_R, graph_T);
+  report("ez", loglik_cont(cont));
+  fix_logistic_parameters(cont, loglik_cont, "lcont", 1e-3);
+  report("erc", loglik_cont(cont));
+
+  report("eco", loglikopt());
+  report("ecm", loglikopt_mono());
+  report("ecrt", loglik_logistic());
+
+  cellcoords();
+  build_disttable();
+
+  cont.setRT(graph_R, graph_T);
+  report("ecz", loglik_cont(cont));
+  fix_logistic_parameters(cont, loglik_cont, "lcont", 1e-3);
+  report("ecrc", loglik_cont(cont));
+
+  long long a = SDL_GetTicks();
+  embedder_loop(20);
+  long long v = SDL_GetTicks();
+  long long tim = v - a;
+
+  report("eeo", loglikopt());
+  report("eem", loglikopt_mono());
+  report("eert", loglik_logistic());
+
+  cellcoords();
+  build_disttable();
+
+  cont.setRT(graph_R, graph_T);
+  fix_logistic_parameters(cont, loglik_cont, "lcont", 1e-3);  
+  report("eerc", loglik_cont(cont));
+
+  print(hlog, "HDR;", separated(";", reps), ";TIME;N;GROWTH\n");
+  print(hlog, "RES'", separated(";", logliks), ";", int(tim), ";", N, ";", cgi.expansion->get_growth());
+  }
+
+}
diff --git a/rogueviz/dhrg/loglik.cpp b/rogueviz/dhrg/loglik.cpp
new file mode 100644
index 00000000..3b3a1156
--- /dev/null
+++ b/rogueviz/dhrg/loglik.cpp
@@ -0,0 +1,434 @@
+// log-likelihood computation
+
+#include <thread>
+#define USE_THREADS
+int threads = 32;
+
+namespace dhrg {
+
+ld llcont_approx_prec = 10000;
+
+// tally edges of the given vertex at the given index
+
+int edgetally[MAXDIST];
+
+void tallyedgesof(int i, int delta, mycell *mc) {
+  using namespace rogueviz;
+  for(auto p: vdata[i].edges) {
+    int j = p.second->i ^ p.second->j ^ i;
+    if(j==i) printf("LOOP!\n");
+    edgetally[quickdist(mc, vertices[j], 0)] += delta;
+    }
+  }
+
+// --- count all edge tallies
+
+void counttallies() {
+  using namespace rogueviz;
+  
+  {
+  progressbar pb(N, "Tallying pairs");
+  for(int i=0; i<N; i++) {
+    mycell* mc = vertices[i];
+    add_to_tally(mc, 1, 0);
+    add_to_set(mc, 1, 0); pb++;
+    if(i % ((N/10)+1) == 0) {
+      memoryInfo();
+      }
+    }
+  }
+
+  {
+  progressbar pb(M, "Tallying edges");
+  
+  for(int u=0; u<MAXDIST; u++) edgetally[u] = 0;
+  
+  for(int i=0; i<N; i++) for(auto p: vdata[i].edges) {
+    int j = p.second->i ^ p.second->j ^ i;
+    if(j < i) { edgetally[quickdist(vertices[i], vertices[j], 0)]++; pb++; }
+    }
+  }
+
+  }
+
+void destroytallies() {
+  progressbar pb(N, "Destroying tallies");
+  for(int i=0; i<N; i++) add_to_set(vertices[i], -1, 0), pb++;
+  for(int i=0; i<MAXDIST; i++)
+    tally[i] = edgetally[i] = 0;
+  }
+
+// log likelihood utilities
+//--------------------------
+
+// MLE of the binomial distribution (a successes, b failures)
+
+ld bestll(ld a, ld b) {
+  if(a == 0 || b == 0) return 0;
+  return a * log(a/(a+b)) + b * log(b/(a+b));
+  }
+
+// a successes, ab total
+
+ld bestll2(ld a, ld ab) { return bestll(a, ab-a); }
+
+// various methods of loglikelihood computation
+
+struct logistic {
+  ld R, T;
+  ld yes(ld d) { return 1/(1 + exp((d-R) / 2 / T)); }
+  ld no(ld d) { return 1/(1 + exp(-(d-R) / 2 / T)); }
+  ld lyes(ld d) { return log(yes(d)); }
+  ld lno(ld d) { return log(no(d)); }
+  void setRT(ld _R, ld _T) { R = _R; T = _T; }
+  };
+
+template<class T> void fix_logistic_parameters(logistic& l, const T& f, const char *name, ld eps) {
+  indenter_finish im("fix_logistic_parameters");
+  ld cur = f(l);
+  println(hlog, format("%s = %20.10" PLDF " (R=%10.5" PLDF " T=%" PLDF ")", name, cur, l.R, l.T));
+
+  for(ld step=1; step>eps; step /= 2) {
+    
+    while(true) { l.R += step; ld t = f(l); if(t <= cur) break; cur = t; }
+    l.R -= step;
+    
+    while(true) { l.R -= step; ld t = f(l); if(t <= cur) break; cur = t; }
+    l.R += step;
+
+    while(true) { l.T += step; ld t = f(l); if(t <= cur) break; cur = t; }
+    l.T -= step;
+    
+    while(true) { l.T -= step; ld t = f(l); if(t <= cur) break; cur = t; }
+    l.T += step;
+
+    println(hlog, format("%s = %20.10" PLDF " (R=%10.5" PLDF " T=%10.5" PLDF ")", name, cur, l.R, l.T));
+    fflush(stdout);
+    }
+  }
+
+logistic current_logistic;
+logistic saved_logistic;
+logistic cont_logistic;
+
+// --- continuous logistic loglikelihood
+// --------------------------------------
+
+vector<hyperpoint> vertexcoords;
+
+// compute vertexcoords from the original embedding data
+void origcoords() {
+  indenter_finish im("origcoords");
+  using namespace rogueviz;
+  vertexcoords.resize(N);
+  for(int i=0; i<N; i++) 
+    vertexcoords[i] = spin(coords[i].second * 2 * M_PI / 360) * xpush(coords[i].first) * C0;
+  }
+
+// compute vertexcoords from the RogueViz representation
+void rvcoords() {
+  indenter_finish im("rvcoords");
+  using namespace rogueviz;
+  vertexcoords.resize(N);
+  for(int i=0; i<N; i++) 
+    vertexcoords[i] = calc_relative_matrix(rogueviz::vdata[i].m->base, currentmap->gamestart(), C0) * rogueviz::vdata[i].m->at * C0;
+  }
+
+// compute vertexcoords from vertices
+void cellcoords() {
+  indenter_finish im("cellcoords");
+  vertexcoords.resize(N);
+  for(int i=0; i<N; i++) {
+    vertexcoords[i] = celltopoint(vertices[i]->ascell()); // calc_relative_matrix(vertices[i]->ascell(), currentmap->gamestart(), C0) * C0;
+    if(isnan(vertexcoords[i][0])) println(hlog, "got NAN for i=", i, " in lev= ", vertices[i]->lev);
+    }
+  }
+
+// needs cellcoords/rvcoords/origcoords
+void build_disttable() {
+  indenter_finish im("build_disttable");
+  int tab[N];
+  for(int i=0; i<N; i++) tab[i] = N;
+  disttable0.clear();
+  disttable1.clear();
+
+  using namespace rogueviz;
+  
+  for(int i=0; i<N; i++) {
+    for(auto p: vdata[i].edges) {
+      int j = p.second->i ^ p.second->j ^ i;
+      if(j<i) tab[j] = i;
+      }
+    for(int j=0; j<i; j++) {
+      ld dist = hdist(vertexcoords[i], vertexcoords[j]);
+      if(dist < 0) continue;
+      (tab[j] == i ? disttable1:disttable0).push_back(dist);
+      }
+    }
+  
+  sort(disttable0.begin(), disttable0.end());
+  sort(disttable1.begin(), disttable1.end());
+  }
+
+// needs build_disttable
+ld loglik_cont(logistic& l = current_logistic) {
+    
+  ld llh = 1;
+  for(auto p: disttable1) llh += l.lyes(p);
+  for(auto p: disttable0) {
+    ld lp = l.lno(p);
+    llh += lp;
+    if(lp > -1e-10) break;
+    }
+  
+  return llh;
+  }
+
+// --- placement loglikelihood
+
+ld loglik_placement() {
+  mycell *root = mroot;
+  ld placement_loglik = 0;  
+  auto seg = getsegment(root,root,0);
+  for(int j=0; j<BOXSIZE; j++) {
+    int qj = seg->qty[j];
+    if(!qj) continue;
+    placement_loglik += qj * (log(qj*1./N) - cgi.expansion->get_descendants(j).log_approx());
+    }
+  return placement_loglik;
+  }
+
+// --- logistic loglikelihood
+
+ld loglik_logistic(logistic& l = current_logistic) {
+  ld result = 0;
+  for(int u=0; u<MAXDIST; u++) if(edgetally[u] && tally[u]-edgetally[u]) {
+    result += edgetally[u] * l.lyes(u) +
+      (tally[u]-edgetally[u]) * l.lno(u);
+    }
+  return result;
+  }
+
+// --- optimal loglikelihood
+
+ld loglikopt() {
+  ld result = 0;
+  for(int u=0; u<MAXDIST; u++) result += bestll2(edgetally[u], tally[u]);
+  return result;
+  }
+
+// --- optimal monotonic loglikelihood
+
+ld loglikopt_mono() {
+  vector<pair<ld, ld> > pairs;
+  ld result = 0;
+  for(int u=0; u<MAXDIST; u++) {
+    auto p = make_pair<ld,ld> (edgetally[u], tally[u]);
+    if(p.second == 0) continue;
+    while(isize(pairs)) {
+      auto pb = pairs.back();
+      if(p.first / p.second > pb.first / pb.second)
+        p.first += pb.first, p.second += pb.second, pairs.pop_back();      
+      else break;
+      }
+    pairs.push_back(p);
+    }
+  for(auto p: pairs)      
+    result += bestll2(p.first, p.second);
+  return result;
+  }
+
+// --- compute loglikelihood according to current method
+
+char lc_type = 'R';
+
+ld loglik_chosen() {
+  switch(lc_type) {
+    case 'O':
+      return loglikopt();
+    case 'R':
+      return loglik_logistic();
+    case 'M':
+      return loglikopt_mono();
+    case 'C':
+      return loglikopt() + loglik_placement();
+    case 'D':
+      return loglikopt_mono() + loglik_placement();
+    default:
+      return loglikopt();
+    }
+  }
+
+// 1e-3 (cont), 1e-6 (normal)
+
+// statistics
+
+void writestats() {
+  indenter_finish im("writestats");
+  memoryInfo();
+  println(hlog, "Vertices by distance (N = ", N, "):");
+  mycell *root = mroot;
+  for(int j=0; j<BOXSIZE; j++) {
+    int qj = getsegment(root,root,0)->qty[j];
+    if(!qj) continue;
+    print(hlog, " ", j, ":", qj);
+    }
+  println(hlog);
+  ld placement_loglik = loglik_placement();
+  
+  for(int u=0; u<MAXDIST; u++) if(tally[u]) {
+    println(hlog, format("* %4d: %8d / %12Ld = %lf %.10" PLDF " %.10" PLDF, 
+      u, edgetally[u], tally[u], double(edgetally[u]) / tally[u], 
+      saved_logistic.yes(u), current_logistic.yes(u)));
+    }
+    
+  println(hlog, "log likelihood\n");
+  ld loglik_chaos = bestll2(M, N*(N-1)/2);
+  ld loglik_opt = loglikopt();
+  ld loglik_mono = loglikopt_mono();
+  ld loglik_rt = loglik_logistic();
+  println(hlog, "  placement     = ", placement_loglik);
+  println(hlog, "  chaos         = ", loglik_chaos);
+  println(hlog, "  optimal any   = ", loglik_opt);
+  println(hlog, "  optimal mono  = ", loglik_mono);
+  println(hlog, "  estimated R/T = ", loglik_logistic(saved_logistic), " (R=", saved_logistic.R, " T=", saved_logistic.T);
+  println(hlog, "  optimal   R/T = ", loglik_rt, " (R=", current_logistic.R, " T=", current_logistic.T);
+  
+  println(hlog, "Compression ratio = %", (placement_loglik+loglik_opt)/loglik_chaos);
+  }
+
+template<class T> auto parallelize(long long N, T action) -> decltype(action(0,0)) {
+#ifndef USE_THREADS
+  return action(0,N);
+#else
+  if(threads == 1) return action(0,N);
+  std::vector<std::thread> v;
+  typedef decltype(action(0,0)) Res;
+  std::vector<Res> results(threads);
+  for(int k=0; k<threads; k++)
+    v.emplace_back([&,k] () { 
+      results[k] = action(N*k/threads, N*(k+1)/threads); 
+      });
+  for(std::thread& t:v) t.join();
+  Res res = 0;
+  for(Res r: results) res += r;
+  return res;
+#endif
+  }
+
+vector<array<ll, 2>> disttable_approx;
+
+ld precise_hdist(hyperpoint vi, hyperpoint vj) {
+  ld da = acosh(vi[2]);
+  ld db = acosh(vj[2]);
+  
+  ld phia = atan2(vi[0], vi[1]);
+  ld phib = atan2(vj[0], vj[1]);
+
+  ld co = sinh(da) * sinh(db) * (1 - cos(phia-phib));
+  // - (vi[0]*vj[0] + vi[1]*vj[1]);
+  
+  ld v = cosh(da - db) + co;
+  if(v < 1) return 0;
+  
+  return acosh(v);
+  }
+
+void build_disttable_approx() {
+  indenter_finish im("build_disttable_approx");
+
+  array<ll, 2> zero = {0, 0};
+
+  using namespace rogueviz;
+
+  std::vector<vector<array<ll, 2>>> results(threads);
+  std::vector<std::thread> v;
+  for(int k=0; k<threads; k++)
+    v.emplace_back([&,k] () {
+      auto& dt = results[k];
+      int tab[N];
+      for(int i=0; i<N; i++) tab[i] = N;
+      auto p = k ? nullptr : new progressbar(N/threads, "build_disttable_approx");
+      for(int i=k; i<N; i+=threads) {
+        if(p) (*p)++;
+        for(auto p: vdata[i].edges) {
+          int j = p.second->i ^ p.second->j ^ i;
+          if(j<i) tab[j] = i;
+          }
+        for(int j=0; j<i; j++) {
+          ld dist = precise_hdist(vertexcoords[i], vertexcoords[j]);
+          if(dist < 0) continue;
+          int dista = dist * llcont_approx_prec;
+          if(isize(dt) < dista+1)
+            dt.resize(dista+1, zero);
+          dt[dista][(tab[j] == i) ? 1 : 0]++;
+          }
+        }
+      if(p) delete p;
+      });
+  for(std::thread& t:v) t.join();
+  
+  int mx = 0;
+  for(auto& r: results) mx = max(mx, isize(r));
+  disttable_approx.clear();
+  disttable_approx.resize(mx, zero);
+  
+  for(auto& r: results)
+    for(int i=0; i<isize(r); i++) 
+      for(int j=0; j<2; j++)
+        disttable_approx[i][j] += r[i][j];
+  }
+
+ld loglik_cont_approx(logistic& l = current_logistic) {
+    
+  ld llh = 0;
+  int N = isize(disttable_approx);
+  for(int i=0; i<N; i++) {
+    if(disttable_approx[i][0])
+      llh += l.lno((i+.5)/llcont_approx_prec) * disttable_approx[i][0];
+    if(disttable_approx[i][1])
+      llh += l.lyes((i+.5)/llcont_approx_prec) * disttable_approx[i][1];
+    }
+  
+  return llh;
+  }
+
+template<class T> void fast_loglik_cont(logistic& l, const T& f, const char *name, ld start, ld eps) {
+
+  indenter_finish im("fix_logistic_parameters");
+  ld cur = f(l);
+  println(hlog, format("%s = %20.10" PLDF " (R=%10.5" PLDF " T=%" PLDF ")\n", name, cur, l.R, l.T));
+
+  map<pair<double, double>, double> memo;
+  auto ff = [&] () {
+    if(memo.count(make_pair(l.R, l.T)))
+      return memo[make_pair(l.R, l.T)];
+    return memo[make_pair(l.R, l.T)] = f(l);
+    };
+
+  for(ld step=start; step>eps; step /= 2) {
+  
+    loop:
+    bool changed = false;        
+    
+    while(true) { l.R += step; ld t = ff(); if(t <= cur) break; cur = t; changed = true; }
+    l.R -= step;
+    
+    while(true) { l.R -= step; ld t = ff(); if(t <= cur) break; cur = t; changed = true; }
+    l.R += step;
+
+    while(true) { l.T += step; ld t = ff(); if(t <= cur) break; cur = t; changed = true; }
+    l.T -= step;
+    
+    while(true) { l.T -= step; ld t = ff(); if(t <= cur) break; cur = t; changed = true; }
+    l.T += step;
+    
+    if(changed) goto loop;
+
+    println(hlog, format("%s = %20.10" PLDF " (R=%10.5" PLDF " T=%10.5" PLDF ")\n", name, cur, l.R, l.T));
+    fflush(stdout);
+    }
+  }
+
+
+}
diff --git a/rogueviz/dhrg/mycell.cpp b/rogueviz/dhrg/mycell.cpp
new file mode 100644
index 00000000..4438276f
--- /dev/null
+++ b/rogueviz/dhrg/mycell.cpp
@@ -0,0 +1,284 @@
+// mycell -- information about the given vertex of a triangulation
+// cell is the relevant struct from HyperRogue; we do not use cell directly to conserve memory
+
+namespace dhrg {
+
+int mycellcount;
+
+struct segmentlist {
+  segment *s;
+  segmentlist *next;
+  };
+
+struct mycell {
+#ifdef BUILD_ON_HR
+  cell *c;
+#else
+  int type;
+#endif
+  int lev;
+  mycell *leftparent, *rightparent;
+  mycell *leftsibling, *rightsibling;
+  mycell *leftchild;
+#ifdef LONG_BRACKETS
+  segmentlist* bracketing;
+#endif
+
+#ifdef BUILD_ON_HR
+  mycell(cell *_c) : c(_c) { 
+    leftparent = rightparent = 
+    leftsibling = rightsibling = 
+    NULL;
+#ifdef LONG_BRACKETS
+    bracketing = NULL;
+#endif
+    mycellcount++;
+    for(int i=0; i<SETS; i++)
+      byleft[i] = byright[i] = NULL;
+    }
+#else
+  mycell() { 
+    leftsibling = rightsibling = leftchild = NULL;
+    mycellcount++;
+    for(int i=0; i<SETS; i++)
+      byleft[i] = byright[i] = NULL;
+#ifdef LONG_BRACKETS
+    bracketing = NULL;
+#endif
+    }
+#endif
+  segment *byleft[SETS];
+  segment *byright[SETS];
+  ~mycell() { mycellcount--; }
+#ifdef BUILD_ON_HR
+  void build();
+  mycell *grightsibling() { return rightsibling; }
+  mycell *gleftsibling() { return leftsibling; }
+  mycell *gleftchild() { return leftchild; }
+  void gchildren() {}
+  cell *ascell() { return c; }
+  int gettype() { return celltype(c); }
+#else
+  void build() {}
+  mycell *grightsibling();
+  mycell *gleftsibling();
+  mycell *gleftchild();
+  void gchildren();
+  cell *ascell();
+  int gettype() { return type; }
+#endif
+  };
+
+#ifdef BUILD_ON_HR
+map<cell*, mycell*> mymap;
+
+mycell *find_mycell(cell *c) {
+  mycell*& mc = mymap[c];
+  if(mc == NULL) mc = new mycell(c);
+  return mc;
+  }
+
+void mycell::build() {
+  const auto m = this;
+  if(m->leftsibling) return; // already computed
+  cell *c2[MAX_EDGE+1];
+  int dist[MAX_EDGE+1];
+  int t = m->c->type;
+  int d = celldist(m->c);
+  m->lev = d;
+  if(d == 0) {
+    m->leftsibling = m->rightsibling = m;
+    m->leftchild = find_mycell(createMov(m->c,0));
+    forCellCM(c2, croot()) find_mycell(c2)->build();
+    }
+  else {    
+    for(int i=0; i<t; i++) c2[i] = createMov(m->c, i);
+    for(int i=0; i<t; i++) dist[i] = celldist(c2[i]);
+    dist[t] = dist[0]; c2[t] = c2[0];
+    for(int i=0; i<t; i++) {
+      if(dist[i] < d && dist[i+1] == d) {
+        m->leftparent = find_mycell(c2[i]);
+        m->leftsibling = find_mycell(c2[i+1]);
+        m->leftchild = find_mycell(c2[(i+2)%t]);
+        }
+      if(dist[i] == d && dist[i+1] < d) {
+        m->rightparent = find_mycell(c2[i+1]);
+        m->rightsibling = find_mycell(c2[i]);
+        }
+      }
+    }
+  }
+#endif
+
+mycell *mroot;
+
+void generate_root() {
+#if BUILD_ON_HR
+  mroot = find_mycell(croot());
+#else
+  int origtype = cgi.expansion->rootid;
+  
+  mroot = new mycell();
+  mroot->lev = 0;
+  mroot->type = origtype;
+  mroot->leftsibling = mroot->rightsibling = mroot;
+  mroot->leftparent = mroot->rightparent = NULL;
+  
+  mycell *child;
+  bool first = true;
+  for(int c: cgi.expansion->children[origtype]) {
+    if(first) {   
+      first = false;
+      mroot->leftchild = child = new mycell();
+      }
+    else {
+      child->rightsibling = new mycell();
+      child->rightsibling->leftsibling = child;
+      child = child->rightsibling;
+      }
+    
+    child->leftparent = child->rightparent = mroot;
+    child->type = c;
+    child->lev = 1;
+    }
+  child->rightsibling = mroot->leftchild;
+  mroot->leftchild->leftsibling = child;
+#endif
+  }
+
+#ifndef BUILD_ON_HR
+
+/* mycell *find_mycell(cell *c) {
+  printf("find_mycell not implemented\n");
+  exit(1);
+  } */
+
+mycell* mycell::gleftsibling() {
+  if(leftsibling) return leftsibling;
+  leftparent->gchildren();
+  if(!leftsibling) {
+    printf("error: no left sibling found\n");
+    exit(1);
+    }
+  return leftsibling;
+  }
+
+mycell* mycell::grightsibling() {
+  if(rightsibling) return rightsibling;
+  rightparent->gchildren();
+  if(!rightsibling) {
+    printf("error: no right sibling found\n");
+    }
+  return rightsibling;
+  }
+
+mycell* mycell::gleftchild() {  
+  if(leftchild) return leftchild;
+  leftchild = new mycell();
+  leftchild->leftparent = gleftsibling();
+  leftchild->rightparent = this;
+  leftchild->lev = lev+1;
+  leftchild->type = cgi.expansion->children[type][0];
+  return leftchild;
+  }
+
+void mycell::gchildren() {
+  mycell *child = gleftchild();
+  if(child->rightsibling) return;  
+  bool first = true;
+  for(int c: cgi.expansion->children[type]) {
+    if(first) {
+      first = false;
+      continue;
+      }
+    child->rightsibling = new mycell();
+    child->rightsibling->leftsibling = child;
+    child = child->rightsibling;
+    child->leftparent = child->rightparent = this;
+    child->type = c;
+    child->lev = lev + 1;
+    }
+  child->rightsibling = grightsibling()->gleftchild();
+  child->rightsibling->leftsibling = child;
+  }
+#endif
+
+vector<mycell*> allchildren(mycell *m, int dir=0) {
+  m->build();
+  vector<mycell*> res;
+  if(m->lev == 0) {
+    mycell *f = mroot->leftchild;
+    int origtype = cgi.expansion->rootid;
+    for(int i: cgi.expansion->children[origtype]) {
+      ignore(i);
+      res.push_back(f);
+      f = f->rightsibling;
+      }      
+    return res;
+    }
+  auto m1 = m->gleftchild();
+  while(true) {
+    m1->build();
+    bool isright = m1->rightparent == m;
+    bool isleft = m1->leftparent == m;
+    if(!isright && !isleft) return res;
+    if(dir > 0 && !isright) ;
+    else if(dir < 0 && !isleft) ;
+    else res.push_back(m1);
+    m1 = m1->grightsibling();
+    }
+  }
+
+vector<mycell*> allneighbors(mycell *m) {
+  auto ret = allchildren(m);
+  if(m->gleftsibling() != m) {
+    ret.push_back(m->leftsibling);
+    ret.push_back(m->grightsibling());
+    }
+  if(m->leftparent) {
+    ret.push_back(m->leftparent);
+    if(m->rightparent != m->leftparent)
+      ret.push_back(m->rightparent);
+    }
+  return ret;
+  }
+
+mycell *find_mycell_by_path(const string& s) {
+  mycell *at = mroot;
+  for(char c: s) {
+    at = at->gleftchild();
+    while(c > '0') c--, at = at->grightsibling();
+    }
+  return at;
+  }
+
+int childindex(mycell *c) {
+  mycell *p = c->rightparent->leftchild;
+  int i = 0;
+  while(p != c) p = p->grightsibling(), i++;
+  return i;
+  }
+
+string get_path(mycell *c) {
+  string s;
+  while(c != mroot) s += '0' + childindex(c), c = c->rightparent;
+  reverse(s.begin(), s.end());
+  return s;
+  }
+
+#ifndef BUILD_ON_HR
+cell *mycell::ascell() {
+  if(lev == 0) return croot();
+  auto m = this; m->build();
+  auto c = rightparent->ascell();
+  int childid = 0;
+  while(m != m->rightparent->leftchild) 
+    childid++, m = m->gleftsibling();
+  if(lev == 1) return createMov(croot(), childid);
+  cell *c2 = ts::child_number(c, childid, celldist);
+
+  return c2;
+  }
+#endif
+
+}
diff --git a/rogueviz/dhrg/paths.cpp b/rogueviz/dhrg/paths.cpp
new file mode 100644
index 00000000..a74e06d7
--- /dev/null
+++ b/rogueviz/dhrg/paths.cpp
@@ -0,0 +1,55 @@
+// computePath to convert HyperRogue's cell to DHRG's mycell (mycell.cpp)
+
+namespace dhrg {
+
+// compute the path to the given cell
+
+hrmap *mapptr;
+
+string computePath(cell *c) {
+  string s = "";
+  int d = celldist(c);
+  while(d > 0) {
+    char z = '0';
+    cell *c1 = ts::left_parent(c, celldist);
+    int id = neighborId(c1, c);
+    d--;
+    if(d == 0) z = '0'+id;
+    else while(true) {
+      id--;
+      if(id<0) id += c1->type;
+      if(celldist(createMov(c1, id)) <= d) break;
+      if(z >= 'a') { printf("<%d> ", celldist(createMov(c1, id))); }
+      z++;
+      }
+    s += z; c = c1;
+    }
+  reverse(s.begin(), s.end());
+  return s;
+  }
+
+int recycle_counter = 0;
+
+void recycle_compute_map() {
+  if(mapptr) {
+    delete mapptr;
+    mapptr = NULL;
+    }
+  }
+
+string computePath(const transmatrix& T) {
+  if(recycle_counter >= 1000)
+    recycle_compute_map();
+  if(!mapptr) mapptr = bt::in() ? bt::new_map() : new hrmap_hyperbolic;
+  recycle_counter++;
+  dynamicval<hrmap*> dv (currentmap, mapptr);
+  cell *c = mapptr->gamestart();
+  hyperpoint T0 = T * C0;
+  // call HyperRogue's function virtualRebase: 
+  // a point in the hyperbolic plane is given by cell c and point T0 relative to c
+  // change c and T0 so that the same point is specified, but with minimal T0
+  virtualRebase(c, T0);
+  return computePath(c);
+  }
+
+}
diff --git a/rogueviz/dhrg/readgraph.cpp b/rogueviz/dhrg/readgraph.cpp
new file mode 100644
index 00000000..1092e54d
--- /dev/null
+++ b/rogueviz/dhrg/readgraph.cpp
@@ -0,0 +1,81 @@
+namespace dhrg {
+  double graph_R, graph_alpha, graph_T;
+  vector<pair<double, double> > coords;
+  
+  rogueviz::edgetype *any;
+  
+  int N;
+               
+  void fixedges() {
+    using namespace rogueviz;
+    for(int i=N; i<isize(vdata); i++) if(vdata[i].m) vdata[i].m->dead = true;
+    for(int i=0; i<isize(vdata); i++) vdata[i].edges.clear();
+    vdata.resize(N);
+    for(auto e: edgeinfos) {
+      e->orig = NULL;
+      addedge(e->i, e->j, e);
+      }
+    }
+  
+  void tst() {}
+
+  void read_graph(string fn, bool subdiv, bool doRebase, bool doStore) {
+
+    any = rogueviz::add_edgetype("embedded edges");
+    rogueviz::fname = fn;
+    fhstream f(fn + "-coordinates.txt", "rt");
+    if(!f.f) {
+      printf("Missing file: %s-coordinates.txt\n", rogueviz::fname.c_str());
+      exit(1);
+      }
+    printf("Reading coordinates...\n");
+    string ignore;
+    if(!scan(f, ignore, ignore, ignore, ignore, N, graph_R, graph_alpha, graph_T)) {
+      printf("Error: incorrect format of the first line\n"); exit(1);
+      }
+    rogueviz::vdata.reserve(N);
+    while(true) {
+      string s = scan<string>(f);
+      if(s == "D11.11") tst();
+      if(s == "" || s == "#ROGUEVIZ_ENDOFDATA") break;
+      int id = rogueviz::getid(s);
+      rogueviz::vertexdata& vd(rogueviz::vdata[id]);
+      vd.name = s;
+      vd.cp = rogueviz::colorpair(rogueviz::dftcolor);
+      
+      double r, alpha;
+      if(!scan(f, r, alpha)) { printf("Error: incorrect format of r/alpha\n"); exit(1); }
+      coords.push_back(make_pair(r, alpha));
+  
+      transmatrix h = spin(alpha * degree) * xpush(r);
+      
+      rogueviz::createViz(id, currentmap->gamestart(), h);
+      }
+    
+    fhstream g(fn + "-links.txt", "rt");
+    if(!g.f) {
+      println(hlog, "Missing file: ", rogueviz::fname, "-links.txt");
+      exit(1);
+      }
+    println(hlog, "Reading links...");
+    int qlink = 0;
+    while(true) {
+      int i = rogueviz::readLabel(g), j = rogueviz::readLabel(g);
+      if(i == -1 || j == -1) break;
+      addedge(i, j, 1, subdiv, any);
+      qlink++;
+      }
+  
+    if(doRebase) {
+      printf("Rebasing...\n");
+      for(int i=0; i<isize(rogueviz::vdata); i++) {
+        if(i % 10000 == 0) printf("%d/%d\n", i, isize(rogueviz::vdata));
+        if(rogueviz::vdata[i].m) virtualRebase(rogueviz::vdata[i].m);
+        }
+      printf("Done.\n");
+      }
+    
+    if(doStore) rogueviz::storeall();
+    }
+  
+  }
diff --git a/rogueviz/dhrg/regular.cpp b/rogueviz/dhrg/regular.cpp
new file mode 100644
index 00000000..9b0538e5
--- /dev/null
+++ b/rogueviz/dhrg/regular.cpp
@@ -0,0 +1,113 @@
+// find the value of D(G) algorithmically (see the paper)
+
+namespace dhrg {
+
+// c2-c1
+
+int cycle_minus(cell *c2, cell *c1) {
+  int acc = 0;
+  while(c1 != c2) c1 = ts::right_of(c1, celldist), acc++;
+  return acc;
+  }
+
+// c2 to the right from c1
+int unlimited_distance(cell *c2, cell *c1) {
+  int at_least = cycle_minus(c2, c1);
+  int steps = 0;
+  while(true) {
+    steps += 2;
+    if(steps >= at_least) return at_least;
+    c1 = ts::right_parent(c1, celldist);
+    c2 = ts::left_parent(c2, celldist);
+    int ndist = steps + cycle_minus(c2, c1);
+    if(ndist < at_least) at_least = ndist;
+    }
+  }
+
+int gettypeof(cell *c) { return type_in_reduced(*(cgi.expansion), c, celldist); }
+
+vector<bool> grow_forever;
+
+set<vector<int>> checked;
+
+bool err = false;
+
+int my_sibling_limit;
+
+void find_sibling_limit(cell *c2, cell *c1) {
+  if(err) return;
+  if(celldist(c2) != celldist(c1)) {
+    printf("not the same ring %d/%d\n", celldist(c1), celldist(c2));
+    c1->item = itSilver;
+    c2->item = itGold;
+    err = true;
+    return;
+    }
+  vector<int> signature;
+  cell *cx = c1;
+  cell *cy = c1;
+  bool gf = false;
+  while(cx != c2) { 
+    int t = gettypeof(cx);
+    if(cx != c1 && grow_forever[t]) gf = true;
+    signature.push_back(t); cy = cx; cx = ts::right_of(cx, celldist);
+    }
+  signature.push_back(gettypeof(cx)); signature.push_back(unlimited_distance(cy, c1) - unlimited_distance(c2, c1));
+  if(checked.count(signature)) return;
+  checked.insert(signature);
+  // for(int v: signature) printf("%d ", v); 
+  int cm = cycle_minus(c2, c1);
+  int ud = c1 == c2 ? -1 : 2 + unlimited_distance(ts::left_parent(c2, celldist), ts::right_parent(c1, celldist));
+  // printf(": %d/%d {%p/%p} [%d]\n", cm, ud, c1, c2, my_sibling_limit);
+  if(cm < ud && cm > my_sibling_limit) { my_sibling_limit = cm; }
+  if(gf) return;
+  int t1 = gettypeof(c1);
+  int t2 = gettypeof(c2);
+  for(int i1=0; i1<isize(cgi.expansion->children[t1]); i1++) 
+  for(int i2=0; i2<isize(cgi.expansion->children[t2]); i2++) 
+    if(c1 != c2 || i1 <= i2+1)
+      find_sibling_limit(ts::child_number(c2, i2+1, celldist), ts::child_number(c1, i1, celldist));
+  }
+
+void correct_sibling_limit() {
+  my_sibling_limit = 0;
+  if(S3 < 4) {
+    grow_forever.clear();
+    grow_forever.resize(cgi.expansion->N, true);
+    
+    while(true) {
+      bool changed = false;
+      for(int i=0; i<cgi.expansion->N; i++) if(grow_forever[i]) {
+        grow_forever[i] = false;
+        if(isize(cgi.expansion->children[i]) == 0)
+          throw hr_exception("Error: our algorithm does not work if some vertices have no tree children");
+        if(isize(cgi.expansion->children[i]) > 1) 
+          for(int c: cgi.expansion->children[i])
+            if(grow_forever[c] || c == i) grow_forever[i] = true;
+        if(!grow_forever[i]) changed = true;
+        }
+      if(!changed) break;
+      }
+     
+    print(hlog, "The following grow forever:"); for(int i=0; i<cgi.expansion->N; i++) if(grow_forever[i]) print(hlog, " ", i); println(hlog);
+  
+    cell *root = currentmap->gamestart();
+    my_sibling_limit = 0;
+    forCellCM(c1, root) forCellCM(c2, root) find_sibling_limit(c2, c1);        
+    }
+  println(hlog, "The correct value of sibling_limit is ", my_sibling_limit);
+  cgi.expansion->sibling_limit = my_sibling_limit;
+  }
+
+void regular_info() {
+  indenter_finish im("regular_info");
+
+  cgi.expansion->get_descendants(0);
+  println(hlog, "growth = ", cgi.expansion->get_growth());
+  
+  // int typecount = expansion.N;
+  
+  correct_sibling_limit();
+  }
+
+}
diff --git a/rogueviz/dhrg/routing.cpp b/rogueviz/dhrg/routing.cpp
new file mode 100644
index 00000000..9bb72b3b
--- /dev/null
+++ b/rogueviz/dhrg/routing.cpp
@@ -0,0 +1,150 @@
+
+//	$(VARIANT) -nogui -dhrg embedded-graphs/facebook_combined_result -contll -iterate 99 -contll -esave > $@
+
+namespace dhrg {
+
+struct pairdata { ld success, route_length; };
+
+#define NOYET 127
+
+vector<vector<char> > actual;
+
+vector<pairdata> pairs;
+
+vector<int> last_goal;
+
+void prepare_pairs() {
+  pairs.resize(N);
+  actual.resize(N);
+  for(int i=0; i<N; i++) actual[i].resize(N, NOYET);
+  for(int i=0; i<N; i++) {
+    vector<int> bfsqueue;
+    auto& p = actual[i];
+    auto visit = [&] (int j, int d) {
+      if(p[j] == NOYET) {
+        p[j] = d;
+        bfsqueue.push_back(j);
+        }
+      };
+    visit(i, 0);
+    for(int k=0; k<isize(bfsqueue); k++) {
+      int a = bfsqueue[k];
+      for(auto ed: rogueviz::vdata[a].edges) {
+        int b = ed.second->i ^ ed.second->j ^ a;
+        visit(b, p[a] + 1);
+        }
+      }
+    }
+  last_goal.clear();
+  last_goal.resize(N, -1);
+  }
+
+void route_from(int src, int goal, const vector<ld>& distances_from_goal) {
+  if(last_goal[src] == goal) return;
+  if(src == goal) {
+    pairs[src].success = 1;
+    pairs[src].route_length = 0;
+    }
+  else {
+    ld bestd = distances_from_goal[src] - 1e-5;
+
+    /* iprintf("route_from goal=%d a=%d, bestd = %f\n", goal, src, bestd + 1e-5);
+    indent += 2; */
+
+    vector<int> candidates;
+    for(auto ed: rogueviz::vdata[src].edges) {
+      int e = ed.second->i ^ ed.second->j ^ src;
+      ld d = e == goal ? -1 : distances_from_goal[e];
+      if(d < bestd) bestd = d, candidates.clear();
+      if(d == bestd) candidates.push_back(e);
+      }
+    pairs[src].success = pairs[src].route_length = 0;
+    for(int c: candidates) {
+      route_from(c, goal, distances_from_goal);
+      // iprintf("candidate = %d\n", c);
+      pairs[src].success += pairs[c].success / isize(candidates);
+      pairs[src].route_length += (1 + pairs[c].route_length) / isize(candidates);
+      }
+    // iprintf("success = %f, route = %f\n", pairs[src].success, pairs[src].route_length);
+    // indent -= 2;
+    }
+  last_goal[src] = goal;
+  }
+
+struct iddata {
+  ld tot, suc, routedist, bestdist;
+  iddata() { tot = suc = routedist = bestdist = 0; }
+  } datas[6];
+
+template<class T> void greedy_routing(int id, const T& distance_function) {
+  for(int goal=0; goal<N; goal++) {
+    vector<ld> distances_from_goal(N);
+    for(int src=0; src<N; src++)
+      distances_from_goal[src] = distance_function(goal, src);
+    for(int src=0; src<N; src++)
+      route_from(src, goal, distances_from_goal);
+      
+    auto& d = datas[id];
+
+    for(int j=0; j<N; j++) if(j != goal){
+      d.tot++;
+      ld p = pairs[j].success;
+      d.suc += p;
+      d.routedist += pairs[j].route_length;
+      d.bestdist += p * actual[goal][j];
+      }
+    }
+  }
+
+void routing_test(string s) {
+
+  vector<string> reps;
+  vector<ld> values;
+
+  auto report = [&] (string s, ld val) {
+    values.push_back(val);
+    reps.push_back(s);
+    println(hlog, "REPORT ", s, " = ", val);
+    };
+
+  if(1) {
+    FILE *f = fopen(("embout/" + s).c_str(), "rb");
+    char buf[999999];
+    int siz = fread(buf, 1, 999999, f);
+    
+    for(int i=0; i<siz; i++) if(is(buf+i, "Embedded Log-likelihood: "))
+      report("bfkl", atof(buf+i+25));
+
+    fclose(f);
+    }
+
+  if(1) {
+    dhrg_init(); read_graph_full("data/sime-" + s);
+    origcoords();
+    prepare_pairs();
+    greedy_routing(0, [] (int i, int j) { return hdist(vertexcoords[i], vertexcoords[j]); });
+    greedy_routing(1, [] (int i, int j) { return quickdist(vertices[i], vertices[j], 0); });
+
+    embedder_loop(20);
+
+    greedy_routing(2, [] (int i, int j) { return quickdist(vertices[i], vertices[j], 0); });
+
+    cellcoords();
+    greedy_routing(3, [] (int i, int j) { return hdist(vertexcoords[i], vertexcoords[j]); });
+
+    greedy_routing(4, [] (int i, int j) { return hdist(vertexcoords[i], vertexcoords[j]) + 100 * quickdist(vertices[i], vertices[j], 0); });
+    }
+
+  for(int id: {0,1,2,3,4}) {
+    string caps[5] = {"ec", "ed", "od", "oc", "tie"};
+    string cap = caps[id];
+    auto& d = datas[id];
+    report("suc_" + cap, d.suc / d.tot);
+    report("str_" + cap, d.routedist / d.bestdist);
+    }
+
+  println(hlog, "HDR;", separated(";", reps), ";N;GROWTH\n");
+  println(hlog, "RES;", separated(";", values), ";", N, ";", cgi.expansion->get_growth());
+  }
+
+}
diff --git a/rogueviz/dhrg/segment.cpp b/rogueviz/dhrg/segment.cpp
new file mode 100644
index 00000000..f4e908b2
--- /dev/null
+++ b/rogueviz/dhrg/segment.cpp
@@ -0,0 +1,261 @@
+// algorithms to compute distances between mycells, and an implementation of tally counter
+
+namespace dhrg {
+
+int segmentcount;
+
+struct qtybox {
+  int minv, maxv;
+  int qty[BOXSIZE];
+  int& operator [] (int i) { 
+    if(i>=BOXSIZE)
+      throw hr_exception("not enough memory reserved -- increase BOXSIZE constant");
+    return qty[i]; 
+    }
+  qtybox() { for(int i=0; i<BOXSIZE; i++) qty[i] = 0; minv = BOXSIZE; maxv = 0; }
+  };
+  
+struct segment {
+  mycell *left;
+  mycell *right;
+  segment *nextleft;
+  segment *nextright;
+  segment *parent;
+  segment *nextchild;
+  segment *firstchild;
+  qtybox qty;
+  int seen;
+  segment() { seen = -1; segmentcount++; }
+  ~segment() { segmentcount--; }
+  };
+
+segment *getsegment(mycell *pleft, mycell *pright, int setid, bool docreate = true) {
+  pleft->build(); pright->build();
+  segment *c1 = pleft->byleft[setid];
+  while(c1) {
+    if(c1->right == pright) {
+      return c1;
+      }
+    c1 = c1 -> nextleft;
+    }
+  if(!docreate) return c1;
+  segment *p = new segment;
+  p->left = pleft;
+  p->right = pright;
+  p->nextleft = pleft->byleft[setid]; pleft->byleft[setid] = p;
+  p->nextright = pright->byright[setid]; pright->byright[setid] = p;
+  p->nextchild = NULL; p->firstchild = NULL;
+  if(pleft->leftparent) {
+    p->parent = getsegment(pleft->leftparent, pright->rightparent, setid);
+    p->nextchild = p->parent->firstchild;
+    p->parent->firstchild = p;
+    }
+  else 
+    p->parent = NULL;
+  #ifdef LONG_BRACKETS
+  if(pleft != pright) {
+    pleft = pleft->grightsibling();
+    int slen = 0;
+    while(pleft != pright) {
+      segmentlist *n = new segmentlist;
+      slen++;
+      n->s = p;
+      n->next = pleft->bracketing;
+      pleft->bracketing = n;
+      pleft = pleft->grightsibling();
+      }
+    if(slen > 10) throw hr_exception("bad bracket");
+    }
+  #endif
+  return p;
+  }
+
+ll tally[MAXDIST];
+
+ll *whichtally = tally;
+
+vector<segment*> acknowledged;
+
+void tallybox(qtybox& box, int d, int mul) {
+  for(int i=box.minv; i<box.maxv; i++)
+    whichtally[i+d] += box[i] * mul;
+  }
+
+void acknowledge(segment *p, int d) {
+  if(!p) return;
+  if(p->seen == -1) {
+    p->seen = d;
+    acknowledged.emplace_back(p);
+    }
+  else if(p->seen > d) 
+    p->seen = d;
+  }
+
+void acknowledgments(int mul) {
+  for(segment* p: acknowledged) {
+    tallybox(p->qty, p->seen, mul);
+    segment *p2 = p->parent;
+    int dist = 1;
+    while(p2) {
+      if(p2->seen != -1) {
+        tallybox(p->qty, p2->seen+dist, -mul);
+        break;
+        }
+      p2=p2->parent; dist++;
+      }
+    p->seen = -1;
+    }
+  acknowledged.clear();
+  }
+
+void trim(qtybox& b) {
+  while(b.minv < b.maxv && !b[b.minv]) b.minv++;
+  if(b.minv == b.maxv) b.minv = BOXSIZE, b.maxv = 0;
+  else while(!b[b.maxv-1]) b.maxv--;
+  }
+
+void add_to_set(mycell *c, int mul, int setid) {
+  // assume mul does not equal 0 !
+  segment *p = getsegment(c, c, setid);
+  int d = 0;
+  while(p) {    
+    p->qty[d] += mul;
+    if(p->qty[d]) {
+      p->qty.minv = min(p->qty.minv, d);
+      p->qty.maxv = max(p->qty.maxv, d+1);
+      }
+    else trim(p->qty);
+    d++; p = p->parent;
+    }
+  }
+
+int segmentdist(segment *p1, segment *p2, int d);
+
+void build_ack(mycell *c, int setid) {
+  segment *p = getsegment(c, c, setid);
+  int d = 0;
+  
+  #ifdef LONG_BRACKETS
+  segmentlist *brackets = c->bracketing;
+  while(brackets) {
+    acknowledge(brackets->s, 0);
+    brackets = brackets->next;
+    }
+  #else
+  acknowledge(getsegment(c->gleftsibling(), c->grightsibling(), setid, false), 0); // larger!
+  #endif
+
+  while(p) {
+    if(p->left != p->right) {
+      mycell *cc = p->left->grightsibling();
+      while(cc != p->right) {
+        acknowledge(getsegment(cc, cc, setid, false), d);
+        cc = cc->grightsibling();
+        }
+      }
+    
+    int sl = cgi.expansion->sibling_limit;
+
+    mycell *mright = p->right;    
+    mycell *mright0 = mright;
+    for(int u=0; u<=sl; u++) {
+      mright->build();
+      segment *mrightbyleft = mright->byleft[setid];
+      while(mrightbyleft) {
+        if(u == 3 && p->right->rightparent == mright->leftparent)
+          acknowledge(mrightbyleft, d+2);
+        else if(u > 3)
+          acknowledge(mrightbyleft, segmentdist(getsegment(mright0, mright0, setid), mrightbyleft, d));
+        else 
+          acknowledge(mrightbyleft, d+u);
+        mrightbyleft = mrightbyleft->nextleft;
+        }
+      mright=mright->grightsibling();
+      }
+    
+    mycell *mleft = p->left;
+    mycell *mleft0 = mleft;
+    for(int u=0; u<=sl; u++) {
+      mleft->build();
+      segment *mleftbyright = mleft->byright[setid];
+      while(mleftbyright) {
+        if(u == 3 && p->left -> leftparent == mleft->rightparent)
+          acknowledge(mleftbyright, d+2);
+        else if(u > 3)
+          acknowledge(mleftbyright, segmentdist(mleftbyright, getsegment(mleft0, mleft0, setid), d));
+        else 
+          acknowledge(mleftbyright, d+u);
+        mleftbyright = mleftbyright->nextright;
+        }
+      mleft=mleft->gleftsibling();
+      }
+ 
+    // go forth
+    d++; p = p->parent;
+    }
+  }
+
+void add_to_tally(mycell *c, int mul, int setid) {
+  build_ack(c, setid);
+  acknowledgments(mul);
+  }
+
+bool is_insegment(mycell *c, segment *p) {
+  mycell *p1 = p->left;
+  mycell *p2 = p->right;
+  while(true) {
+    if(c == p1) return true;
+    if(p1 == p2) return false;
+    p1 = p1->grightsibling();
+    }
+  }
+
+int segmentdist(segment *p1, segment *p2, int d) {
+
+  if(p1->left == p1->right && is_insegment(p1->left, p2)) return d;
+  if(p2->left == p2->right && is_insegment(p2->left, p1)) return d;
+  
+  int best = 1000000;
+  
+  while(true) {
+    if(d >= best) return best;
+    mycell *mright;
+    
+    int sl = cgi.expansion->sibling_limit;
+    
+    mright = p1->right;
+    for(int u=0; u<=sl; u++) {
+      mright->build();
+      if(mright == p2->left) best = min(best, d+u);
+      mright = mright->grightsibling();
+      }
+
+    mright = p2->right;
+    for(int u=0; u<=sl; u++) {
+      mright->build();
+      if(mright == p1->left) best = min(best, d+u);
+      mright = mright->grightsibling();
+      }
+    
+    p1 = p1->parent;
+    p2 = p2->parent;
+    d += 2;
+    }
+  }
+
+int quickdist(mycell *c1, mycell *c2, int setid=0) {
+  int d = 0;
+  c1->build();
+  c2->build();
+  int d1 = c1->lev;
+  int d2 = c2->lev;
+  segment *p1 = getsegment(c1, c1, setid);
+  segment *p2 = getsegment(c2, c2, setid);
+  
+  while(d1>d2) { p1 = p1->parent; d1--; d++; }
+  while(d2>d1) { p2 = p2->parent; d2--; d++; }
+  
+  return segmentdist(p1, p2, d);
+  }
+  
+}
diff --git a/rogueviz/dhrg/tests.cpp b/rogueviz/dhrg/tests.cpp
new file mode 100644
index 00000000..2edbc56f
--- /dev/null
+++ b/rogueviz/dhrg/tests.cpp
@@ -0,0 +1,91 @@
+// check whether our distance and tallycounter algorithms work correctly
+
+namespace dhrg {
+
+void test_paths(int radius) {
+  celllister cl(croot(), radius, 1000000, NULL);
+  int N = isize(cl.lst);
+  printf("N = %d\n", N);
+  vector<mycell*> mycells;
+  for(cell *c: cl.lst)
+    mycells.push_back(find_mycell_by_path(computePath(c)));
+  
+  int ctable[MAXDIST];
+
+  for(int u=0; u<MAXDIST; u++)
+    ctable[u] = 0;
+  
+  int errorcount = 0;
+
+  for(int i=0; i<N; i++) {
+    cell *c = cl.lst[i];
+    mycell *mc = mycells[i];
+    cell *c1 = mc->ascell();
+    if(c != c1) {
+      printf("ascell error %s / %s\n", computePath(c).c_str(), computePath(c1).c_str());
+      c->item = itEmerald;
+      c1->item = itRuby;
+      errorcount++;
+      return;
+      }
+    }
+  
+  for(int i=0; i<N; i++)
+    for(int j=0; j<N; j++) {
+      int a = celldistance(cl.lst[i], cl.lst[j]);
+      int b = quickdist(mycells[i], mycells[j]);
+      if(a != b) {
+        printf("celldistance = %d\n", a);
+        printf("quickdist = %d\n", b);
+        celllister cl2(cl.lst[i], 10, 1500, cl.lst[j]);
+        if(cl2.listed(cl.lst[j]))
+          printf("actual distance = %d\n", cl2.getdist(cl.lst[j]));
+        cl.lst[i]->item = itDiamond;
+        cl.lst[j]->item = itSilver;
+        errorcount++;
+        return;
+        }
+      else ctable[a]++;
+      }
+  
+  for(mycell* mc1: mycells) for(mycell* mc2: mycells) {
+    add_to_set(mc1, 1, 0);
+    add_to_tally(mc2, 1, 0);
+    // int v = 0;
+    if(tally[quickdist(mc1, mc2)] != 1) {
+      printf("[%p] [%p]\n", mc1, mc2);
+      printf("quickdist = %d\n", quickdist(mc1, mc2));
+      for(int i=0; i<MAXDIST; i++) if(tally[i]) printf("in tally = %d\n", i);
+      mc1->ascell()->item = itDiamond;
+      mc2->ascell()->item = itSilver;
+      errorcount++;
+      add_to_tally(mc2, -1, 0);
+      add_to_set(mc1, -1, 0);
+      return;
+      }
+    add_to_tally(mc2, -1, 0);
+    add_to_set(mc1, -1, 0);
+    }
+
+  for(mycell* mc: mycells)
+    add_to_set(mc, 1, 0);
+
+  for(mycell* mc: mycells)
+    add_to_tally(mc, 1, 0);
+
+  for(int u=0; u<MAXDIST; u++) if(ctable[u] || tally[u]) {
+    printf("%d: %d/%d\n", u, ctable[u], (int) tally[u]);
+    if(ctable[u] != tally[u]) errorcount++;
+    }
+  
+  for(mycell* mc: mycells)
+    add_to_tally(mc, -1, 0);
+
+  for(mycell* mc: mycells)
+    add_to_set(mc, -1, 0);
+
+  if(errorcount)
+    printf("errors found: %d\n", errorcount);
+  }
+
+}
diff --git a/rogueviz/dhrg/visualize.cpp b/rogueviz/dhrg/visualize.cpp
new file mode 100644
index 00000000..3cd4a20f
--- /dev/null
+++ b/rogueviz/dhrg/visualize.cpp
@@ -0,0 +1,116 @@
+// a RogueViz dialog to manipulate DHRG embeddings
+
+int held_id = -1;
+
+void show_likelihood() {
+  cmode = sm::SIDE | sm::MAYDARK | sm::DIALOG_STRICT_X;
+  gamescreen();
+
+  dialog::init("DHRG information");
+  
+  ll bonus_tally[MAXDIST], bonus_edgetally[MAXDIST];
+
+  if(held_id >= 0) {
+    for(int i=0; i<MAXDIST; i++)
+      bonus_tally[i] = tally[i],
+      bonus_edgetally[i] = edgetally[i];
+    auto& mc = vertices[held_id];
+    add_to_set(mc, -1, 0);
+    add_to_tally(mc, -1, 0);
+    tallyedgesof(held_id, -1, mc);
+    for(int i=0; i<MAXDIST; i++)
+      bonus_tally[i] -= tally[i],
+      bonus_edgetally[i] -= edgetally[i];
+    tallyedgesof(held_id, +1, mc);
+    add_to_tally(mc, +1, 0);    
+    add_to_set(mc, +1, 0);
+    }
+  else 
+    for(int i=0; i<MAXDIST; i++) bonus_tally[i] = 0;
+
+  for(int u=0; u<MAXDIST; u++) if(tally[u] || bonus_tally[u]) {
+    char buf[20];
+    sprintf(buf, "%.6lf", lc_type == 'R' ? current_logistic.yes(u) : double(edgetally[u] * 1. / tally[u]));    
+    string s = its(u);
+    if(isize(s) == 1) s = "0" + s;
+    s += ": " + its(edgetally[u]) + " / " + its(tally[u]);
+    if(bonus_tally[u]) s += " [" + its(bonus_edgetally[u]) + "/" + its(bonus_tally[u]) + "]";
+    dialog::addSelItem(s, buf, 0);
+    }
+  
+  char letters[3] = {'O', 'R', 'M'};
+  string lltypes[3] = {"opt", "logistic", "mono."};
+  string clltype = "?";
+  char nextletter;
+  for(int i=0; i<3; i++) if(lc_type == letters[i])
+    clltype = lltypes[i], nextletter = letters[(i+1)%3];
+  
+  getcstat = '-';
+
+  dialog::addBreak(50);
+  dialog::addSelItem("loglikelihood (" + clltype + ")", fts(loglik_chosen()), 'l');
+  dialog::add_action([nextletter] () { lc_type = nextletter; });
+  
+  dialog::addSelItem("iterations of local search", its(iterations), 'i');
+  dialog::add_action([] () { 
+    embedder_loop(1); 
+    ts_vertices = ts_rbase;
+    place_rogueviz_vertices();
+    if(!stored) rogueviz::storeall(), stored = true;
+    else shmup::fixStorage();
+    });
+  
+  dialog::addItem("move", 'm');
+  dialog::add_action([] () { popScreen(); });
+
+  dialog::display();
+
+  if(held_id >= 0) {
+    if(!mousepressed) held_id = -1;
+    else if(mouseover) {    
+      rogueviz::vdata[held_id].m->base = mouseover;
+      shmup::fixStorage();
+      dhrg::fixedges(); 
+
+      auto& mc = vertices[held_id];
+      tallyedgesof(held_id, -1, mc);
+      add_to_set(mc, -1, 0);
+      add_to_tally(mc, -1, 0);
+
+      mc = find_mycell_by_path(computePath(rogueviz::vdata[held_id].m->base));
+
+      tallyedgesof(held_id, 1, mc);
+      add_to_tally(mc, 1, 0);
+      add_to_set(mc, 1, 0);
+      }
+    }
+  
+  keyhandler = [] (int sym, int uni) {
+
+    if(uni == '-' && held_id == -1) {
+      for(int i=0; i<N; i++) if(rogueviz::vdata[i].m->base == mouseover)
+        held_id = i;
+      return;
+      }
+
+    dialog::handleNavigation(sym, uni);    
+    };
+  }
+
+bool dhrg_animate(int sym, int uni) {
+  if((cmode & sm::NORMAL) && uni == '/') {
+    clearMessages();
+
+    if(ts_rbase > ts_vertices) {
+      dhrg_init(); graph_from_rv(); 
+      ts_vertices = ts_rbase;
+      place_rogueviz_vertices();
+      if(!stored) rogueviz::storeall(), stored = true;
+      else shmup::fixStorage();
+      }    
+    pushScreen(show_likelihood);
+
+    return true;
+    }
+  return false;
+  }
diff --git a/rogueviz/rogueviz-all.cpp b/rogueviz/rogueviz-all.cpp
index f61d0524..795cae0e 100644
--- a/rogueviz/rogueviz-all.cpp
+++ b/rogueviz/rogueviz-all.cpp
@@ -76,6 +76,7 @@
 
 // the following comments are read by mymake so that it knows that the files include other files there:
 // hidden dependencies: rogueviz/nilrider/
+// hidden dependencies: rogueviz/dhrg/
 
 //#endif