From 1f27ae428a6653e8eba6d62583d7c7710a0e0659 Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Tue, 27 Jan 2026 10:31:01 +0100
Subject: [PATCH] extra tools to compare embeddings

---
 rogueviz/embeddings/compare.cpp    | 348 +++++++++++++++++++++++++++++
 rogueviz/embeddings/embeddings.cpp |   3 +
 rogueviz/embeddings/embeddings.h   |   2 +
 rogueviz/embeddings/maprank.cpp    |  72 ++++++
 rogueviz/embeddings/mdl.cpp        | 169 ++++++++++++++
 5 files changed, 594 insertions(+)
 create mode 100644 rogueviz/embeddings/compare.cpp
 create mode 100644 rogueviz/embeddings/maprank.cpp
 create mode 100644 rogueviz/embeddings/mdl.cpp

diff --git a/rogueviz/embeddings/compare.cpp b/rogueviz/embeddings/compare.cpp
new file mode 100644
index 00000000..81c3e664
--- /dev/null
+++ b/rogueviz/embeddings/compare.cpp
@@ -0,0 +1,348 @@
+#include "../som/kohonen.h"
+
+namespace rogueviz {
+
+namespace embeddings {
+
+int count_directed_edges() {
+  int M = 0;
+  for(auto& v: directed_edges) M += isize(v);
+  return M;
+  }
+
+void full_evaluation() {
+
+  println(hlog, "evaluating embedding: " + current->name());
+
+  try {
+    ld maxradius = 0;  
+    for(int i=0; i<isize(rogueviz::vdata); i++) maxradius = max(maxradius, current->zero_distance(i));
+    println(hlog, "maximum radius = ", maxradius);
+    }
+  catch(hr_exception&) {}
+
+  bool symmetric = count_directed_edges() == 2 * isize(edgeinfos);
+  if(symmetric) full_routing();
+  continuous_ranks();
+  analyze_mdl_symmetric(symmetric);
+  }
+
+struct abstract_embedding : embedding {
+  hyperpoint as_hyperpoint(int id) override { throw hr_exception("abstract embedding"); }
+  pair<cell*, hyperpoint>  as_location(int id) override { throw hr_exception("abstract embedding"); }
+  };
+
+/* convert a hyperbolic embedding into a landscape embedding */
+
+using rogueviz::kohonen::kohvec;
+using rogueviz::kohonen::columns;
+
+int landscape_dim;
+
+vector<kohvec> landscape_weights;
+
+ld dist(const kohvec& a, const kohvec& b) {
+  ld res = 0;
+  for(int i=0; i<columns; i++) res += (a[i]-b[i]) * (a[i] - b[i]);
+  return res;
+  }
+
+struct landscape_embedding : abstract_embedding {
+  virtual string name() override { return "landscape"; }
+  eGeometryClass get_geometry() override { return gcEuclid; }
+  int get_dimension() override { return isize(landscape_weights[0]); }
+  ld distance(int i, int j) override {
+    return dist(landscape_weights[i], landscape_weights[j]);
+    }
+  };
+
+void landscape_build(int dim) {
+  rogueviz::embeddings::init_landscape(columns = landscape_dim = dim);
+  int N = isize(rogueviz::vdata);
+  landscape_weights.resize(N);
+  if(1) {
+    progressbar pb(N, "compute landscape");
+    for(int i=0; i<N; i++) {
+      landscape_weights[i] = rogueviz::embeddings::get_landscape_at(current->as_location(i).first);
+      pb++;
+      }
+    }
+  println(hlog, "delta size = ", isize(rogueviz::embeddings::delta_at));
+  enable_embedding(make_shared<landscape_embedding>());
+  }
+
+void landscape_output(string shape, ld target_dist, string fname) {
+  ld max_dist = 0;
+  int N = isize(rogueviz::vdata);
+  for(int i=0; i<N; i++) {
+    ld tot = 0;
+    for(int d=0; d<columns; d++) tot += pow(landscape_weights[i][d], 2);
+    tot = sqrt(tot);
+    max_dist = max(max_dist, tot);
+    }
+  println(hlog, "max_dist = ", max_dist, " to ", target_dist);
+  fhstream f(fname, "w");
+  for(int i=0; i<N; i++) {
+    ld tot = 0;
+    println(f, rogueviz::vdata[i].name);
+    for(int d=0; d<columns; d++) {
+      auto w = landscape_weights[i][d] * target_dist / max_dist;
+      if(d) print(f, " ");
+      print(f, w);
+      tot += w * w;
+      }
+    if(shape == "euclid") println(f);
+    else if(shape == "poincare") println(f);
+    else if(shape == "lorentz") println(f, " ", sqrt(tot+1));
+    else if(shape == "horosphere") {
+      tot /= 2;
+      ld shift = -asinh(target_dist)/2;
+      println(f, " ", cosh(shift) * tot + sinh(shift) * (1+tot), " ", cosh(shift) * (1+tot) + sinh(shift) * tot);
+      }
+    }
+  }
+
+/* read a distance table */
+
+struct disttable_embedding : abstract_embedding {
+  vector<vector<ld>> disttable;
+  virtual string name() override { return "disttable"; }
+  eGeometryClass get_geometry() override { return gcHyperbolic; }
+  int get_dimension() override { return 1; }
+  ld zero_distance(int id) override { throw hr_exception("unknown"); }
+  ld distance(int i, int j) override { return disttable[i][j]; }
+  };
+
+struct zdisttable_embedding : disttable_embedding {
+  vector<ld> dist0;
+  virtual string name() override { return "zdisttable"; }
+  ld zero_distance(int id) override { return dist0[id]; }
+  };
+
+void read_disttable(string nodelist_fname, string disttable_fname) {
+  auto de = make_shared<disttable_embedding> ();
+  int N = isize(rogueviz::vdata);
+
+  fhstream g(nodelist_fname, "rt");
+  vector<int> our_index;
+  while(true) {
+    int i = rogueviz::readLabel(g);
+    if(i == -1) break;
+    our_index.push_back(i);
+    }
+  println(hlog, "read ", isize(our_index), " node labels from ", nodelist_fname);
+  if(isize(our_index) != N) { println(hlog, "wrong number of labels, ", N, " expected"); return; }
+
+  ld maxdist = 0;
+
+  de->disttable.resize(N);
+  for(int i=0; i<N; i++) de->disttable[i].resize(N, 0);
+  fhstream g1(disttable_fname, "rt");
+  println(hlog, "reading distance table from ", disttable_fname);
+  for(int i=0; i<N; i++)
+  for(int j=0; j<i; j++) {
+    ld dist = scan<ld>(g1);
+    maxdist = max(dist, maxdist);
+    de->disttable[our_index[i]][our_index[j]] = dist;
+    de->disttable[our_index[j]][our_index[i]] = dist;
+    }
+
+  println(hlog, "maximum distance = ", maxdist);
+
+  enable_embedding(de);
+  }
+
+/* read a 3D Mercator embeding */
+
+struct mercator_coord {
+  ld distance = -1;
+  array<ld, 3> h;
+  };
+
+struct dmercator_embedding : abstract_embedding {
+  vector<mercator_coord> mcs;
+  virtual string name() override { return "dmercator"; }
+  eGeometryClass get_geometry() override { return gcHyperbolic; }
+  int get_dimension() override { return 3; }
+  hyperpoint as_hyperpoint(int id) {
+    auto& mc = mcs[id];
+    hyperpoint h;
+    for(int i=0; i<3; i++) 
+      h[i] = sinh(mc.distance) * mc.h[i];
+    h[3] = cosh(mc.distance);
+    return h;
+    }
+
+  ld zero_distance(int id) override { return mcs[id].distance; }
+
+  ld distance(int a, int b) {
+    if(a == b) return ld(0);
+    ld da = mcs[a].distance;
+    ld db = mcs[b].distance;
+
+    ld cosphi = 0;
+    for(int i=0; i<3; i++) cosphi += mcs[a].h[i] * mcs[b].h[i];
+  
+    ld co = sinh(da) * sinh(db) * (1 - cosphi);  
+    ld v = cosh(da - db) + co;
+    if(v < 1) return ld(0);
+  
+    return acosh(v);
+    };
+
+  };
+
+
+void read_dmercator(string fname, bool only_read = false) {
+  auto me = make_shared<dmercator_embedding> ();
+  int N = isize(rogueviz::vdata);
+
+  me->mcs.resize(N);
+
+  fhstream g(fname, "rt");
+
+  while(true) {
+    string s = scan<string>(g);
+    if(s == "") break;
+    if(s == "#") { scanline_noblank(g); continue; }
+    int id = rogueviz::getid(s);
+    if(id == -1) { printf("got -1\n"); exit(1); }
+    ld kappa = scan<ld> (g); hr::ignore(kappa);
+    auto& mc = me->mcs[id];
+    mc.distance = scan<ld> (g);
+    for(int i=0; i<3; i++) mc.h[i] = scan<ld> (g);
+    ld radius = mc.h[0] * mc.h[0] + mc.h[1] * mc.h[1] + mc.h[2] * mc.h[2];
+    radius = sqrt(radius);
+    for(int i=0; i<3; i++) mc.h[i] /= radius;
+    }
+
+  for(auto& p: me->mcs) if(p.distance < 0) throw hr_exception("mercator read incorrectly");
+
+  enable_embedding(me);
+  }
+
+/* read a high-dimensional Euclidean embeding */
+
+struct euclid_embedding : zdisttable_embedding {
+  vector<vector<ld>> coords;
+  virtual string name() override { return "euclid"; }
+  eGeometryClass get_geometry() override { return gcEuclid; }
+  int get_dimension() override { return isize(coords[0]); }
+  };
+
+void read_euclid(string fname, int dim) {
+  auto ee = make_shared<euclid_embedding> ();
+  int N = isize(rogueviz::vdata);
+
+  println(hlog, "Open file ", fname);
+  fhstream g(fname, "rt");
+  if(!g.f) return file_error(fname);
+  println(hlog, "Reading Euclid (dim = ", dim, ")...");
+  ee->coords.resize(N);
+  ee->dist0.resize(N);
+  while(true) {
+    int i = rogueviz::readLabel(g);
+    if(i == -1) break;
+    vector<ld> co(dim);
+    for(int i=0; i<dim; i++) co[i] = scan<ld>(g);
+
+    ld d = 0; for(int i=0; i<dim; i++) d += co[i] * co[i];
+    ee->dist0[i] = sqrt(d);
+
+    ee->coords[i] = std::move(co);
+    }
+
+  ee->disttable.resize(N);
+  for(int i=0; i<N; i++) ee->disttable[i].resize(N, 0);
+  for(int i=0; i<N; i++)
+  for(int j=0; j<=i; j++) {
+    ld d = 0; for(int k=0; k<dim; k++) d += pow(ee->coords[i][k]-ee->coords[j][k], 2);
+    ee->disttable[i][j] = ee->disttable[j][i] = sqrt(d);
+    }
+
+  enable_embedding(ee);
+  }
+
+/* read a high-dimensional Poincare embeding */
+
+struct poincare_embedding : zdisttable_embedding {
+  vector<vector<ld>> coords;
+  virtual string name() override { return "poincare"; }
+  eGeometryClass get_geometry() override { return gcHyperbolic; }
+  int get_dimension() override { return isize(coords[0])-1; }
+  };
+
+void read_poincare(string fname, int dim) {
+  auto pe = make_shared<poincare_embedding>();
+  int N = isize(rogueviz::vdata);
+
+  println(hlog, "Open file ", fname);
+  fhstream g(fname, "rt");
+  if(!g.f) return file_error(fname);
+  println(hlog, "Reading Poincare (dim = ", dim, ")...");
+  pe->coords.resize(N);
+  pe->dist0.resize(N);
+
+  while(true) {
+    int i = rogueviz::readLabel(g);
+    if(i == -1) break;
+    vector<ld> co(dim+1);
+    for(int i=0; i<dim; i++) co[i] = scan<ld>(g);
+
+    ld d = 0; for(int i=0; i<dim; i++) d += co[i] * co[i];
+    hyperpoint h; h[0] = sqrt(d); h[1] = 0; auto h1 = perspective_to_space(h, 1);
+
+    for(int i=0; i<dim; i++) co[i] = co[i] * h1[0] / d; co[dim] = h1[2];
+
+    pe->dist0[i] = acos_auto_clamp(h1[2]);
+    pe->coords[i] = std::move(co);
+    }
+
+  pe->disttable.resize(N);
+  for(int i=0; i<N; i++) pe->disttable[i].resize(N, 0);
+  for(int i=0; i<N; i++)
+  for(int j=0; j<=i; j++) {
+    ld d = 0; for(int k=0; k<dim; k++) d += pow(pe->coords[i][k]-pe->coords[j][k], 2);
+    d -= pow(pe->coords[i][dim]-pe->coords[j][dim], 2);
+    pe->disttable[i][j] = pe->disttable[j][i] = acos_auto_clamp(d);
+    }
+
+  enable_embedding(pe);
+  }
+
+int compare_args() {
+  using namespace arg;
+
+  if(argis("-fulleval")) {
+    full_evaluation();
+    }
+  else if(argis("-el-disttable")) {
+    shift(); string nodelist_fname = args();
+    shift(); string disttable_fname = args();
+    read_disttable(nodelist_fname, disttable_fname);
+    }
+  else if(argis("-el-dmercator")) {
+    shift(); read_dmercator(args());
+    }
+  else if(argis("-el-poincare")) {
+    shift(); string fname = args();
+    shift(); int dim = argi();
+    read_poincare(fname, dim);
+    }
+  else if(argis("-el-euclid")) {
+    shift(); string fname = args();
+    shift(); int dim = argi();
+    read_euclid(fname, dim);
+    }
+  else if(argis("-build-landscape")) {
+    shift(); landscape_build(argi());
+    }
+  else return 1;
+
+  return 0;
+  }
+
+auto acmp = addHook(hooks_args, 50, compare_args);
+
+}
+}
diff --git a/rogueviz/embeddings/embeddings.cpp b/rogueviz/embeddings/embeddings.cpp
index 88dea487..c182eefa 100644
--- a/rogueviz/embeddings/embeddings.cpp
+++ b/rogueviz/embeddings/embeddings.cpp
@@ -118,3 +118,6 @@ int a = arg::add3("-edgelist", [] { arg::shift(); read_edgelist(arg::args()); })
 #include "loglik.cpp"
 #include "routing.cpp"
 
+#include "maprank.cpp"
+#include "mdl.cpp"
+#include "compare.cpp"
diff --git a/rogueviz/embeddings/embeddings.h b/rogueviz/embeddings/embeddings.h
index 4d7d7d87..afbc8ede 100644
--- a/rogueviz/embeddings/embeddings.h
+++ b/rogueviz/embeddings/embeddings.h
@@ -81,6 +81,8 @@ std::vector<int> path(int src);
 int get_actual(int src);
 void prepare_goal(int goal);
 
+int count_directed_edges();
+
 void read_edgelist(const string& fn);
 void read_polar(const string& fn);
 }
diff --git a/rogueviz/embeddings/maprank.cpp b/rogueviz/embeddings/maprank.cpp
new file mode 100644
index 00000000..cf30cf56
--- /dev/null
+++ b/rogueviz/embeddings/maprank.cpp
@@ -0,0 +1,72 @@
+// compute mAP/MeanRank rankings, for general continuous distances
+
+namespace rogueviz {
+
+namespace embeddings {
+
+void continuous_ranks() {
+  ld tot_ranks = 0, tot_rby = 0, tot_map = 0, tot_n = 0;
+
+  int N = isize(rogueviz::vdata);
+
+  if(1) {
+    progressbar pb(N/threads, "continuous ranks");
+
+    std::mutex lock; 
+    parallelize(N, [&] (int a, int b) {
+      ld ltot_ranks = 0, ltot_rby = 0, ltot_map = 0, ltot_n = 0;
+      for(int i=a; i<b; i++) {
+        int qty = isize(directed_edges[i]);
+        if(!qty) continue;
+        vector<pair<ld, int>> dists(N);
+        // println(hlog, i, ": ", dists);
+        for(int j=0; j<N; j++) dists[j] = { current->distance(i, j), 0};
+
+        for(auto e: directed_edges[i]) dists[e].second = 1;
+
+/*      for(auto& d: dists) if(d.second == 0) d.first += 0.5;
+        vector<int> ba(MAXDIST, 0), bb(MAXDIST, 0);
+        for(auto& d: dists) { if(d.second == 1) ba[d.first]++; bb[d.first]++; } */
+
+/*    bb[0]--; print(hlog, i, ": "); for(int d=0; d<MAXDIST; d++) if(bb[d]) print(hlog, d, " : ", ba[d],"/",bb[d], " ");
+    println(hlog); */
+
+        // dists[i].first -= 2;
+        sort(dists.begin(), dists.end());
+
+        vector<int> ranks;
+        int seen = 0;
+        int bad = 0;
+        ld for_map = 0;
+        // note: i itself is in first position
+        for(int j=0; seen<qty; j++)
+          if(dists[j].second) {
+            ltot_ranks += bad;
+            ranks.push_back(bad+seen);
+            // if(bad+seen+1 <= qty) for_map++; // 
+            for_map += (seen+1.) / (bad+seen);
+            // printf("%d/%d ", seen+1, bad+seen);
+            seen++;
+            }
+          else
+            bad++;
+        if(a == 0) pb++;
+        ltot_rby += qty;
+        ltot_map += for_map / qty;
+// println(hlog, "contribution: ", for_map / qty);
+        ltot_n++;
+        }
+      std::unique_lock<std::mutex> lk(lock);
+      tot_rby += ltot_rby;
+      tot_map += ltot_map;
+      tot_ranks += ltot_ranks;
+      tot_n += ltot_n;
+      return 0;
+      });
+    }
+  println(hlog, "MeanRank = ", tot_ranks / tot_rby, " MAP = ", tot_map / tot_n);
+  println(hlog, "data: ", tie(tot_rby, tot_map, tot_ranks, tot_n, N));
+  }
+
+
+}}
diff --git a/rogueviz/embeddings/mdl.cpp b/rogueviz/embeddings/mdl.cpp
new file mode 100644
index 00000000..c73f5473
--- /dev/null
+++ b/rogueviz/embeddings/mdl.cpp
@@ -0,0 +1,169 @@
+namespace rogueviz {
+
+namespace embeddings {
+
+int tal_limit = 0;
+
+void analyze_mdl_symmetric(bool symmetric) {
+  int dim = current->get_dimension();
+  int N = isize(rogueviz::vdata);
+
+  if(dim == 1) return;
+
+  vector<ld> center_distances(N);
+  for(int id=0; id<N; id++) center_distances[id] = current->zero_distance(id);
+
+  build_disttable_approx();
+  logistic cont;
+  cont.setRT(0, 1);
+  if(symmetric) fast_loglik_cont(cont, loglik_cont_approx, "lcont", 1, 1e-6);
+  if(symmetric) println(hlog, "loglikelihood = ", format("%.1f", -loglik_cont_approx(cont)));
+  ld N1 = N * (N-1.);
+  ld M = count_directed_edges();
+  if(symmetric) { N1 /= 2; M /= 2; }
+  ld max_radius = 0;
+  for(auto& d: center_distances) max_radius = max(max_radius, d);
+  ld M1 = N1 - M;
+  ld entropy = M * log(M/N1) + M1 * log(M1/N1);
+  println(hlog, "N = ", N, " M = ", M, " radius = ", max_radius, " entropy = ", format("%.1f", -entropy));
+  
+  auto entropy_for_exponent = [&] (ld alpha) {
+    // integral of exp(alpha*t) from 0 to max_radius
+    ld integral = alpha ? (exp(alpha*max_radius)-1) / alpha : max_radius;
+    ld entropy = 0;
+    for(auto& d: center_distances) entropy += -log(exp(alpha*d) / integral);
+    return entropy;
+    };
+
+  // we minimize entropy
+  ld alpha_min = -3, alpha_max = +3;
+  for(int it=0; it<100; it++) {
+    ld a = (alpha_min * 2 + alpha_max) / 3;
+    ld b = (alpha_min + 2 * alpha_max) / 3;
+    if(entropy_for_exponent(a) < entropy_for_exponent(b)) alpha_max = b; else alpha_min = a;
+    }
+
+  ld alpha = alpha_min;
+  println(hlog, "distance entropy for zero alpha = ", entropy_for_exponent(0));
+  println(hlog, "obtained alpha = ", alpha);
+  println(hlog, "distance entropy = ", entropy_for_exponent(alpha));
+
+  vector<vector<ld>> bad;
+  vector<vector<ld>> good;
+  vector<logistic> logistics(N, cont);
+
+  if(!symmetric) {
+    if(tal_limit == 0) { println(hlog, "not computing TAL -- the graph is asymmetric"); return; }
+    if(N >= tal_limit) { println(hlog, "too large graph to compute TAL"); return; }
+
+    println(hlog, "computing TAL");
+    bad.resize(N);
+    good.resize(N);
+    vector<int> last_id(N, -1);
+    ld total_loglik = 0;
+    for(int i=0; i<N; i++) {
+      for(int j: directed_edges[i]) last_id[j] = i;
+      for(int j=0; j<N; j++) {
+        if(last_id[j] == i)
+          good[i].push_back(current->distance(i, j));
+        else
+          bad[i].push_back(current->distance(i, j));
+        }
+      sort(good[i].begin(), good[i].end());
+      sort(bad[i].begin(), bad[i].end());
+
+      auto fun = [&] (logistic& l) {
+        ld res = 0;
+        for(auto& g: good[i]) res += l.lyes(g);
+        for(auto& b: bad[i]) res += l.lno(b);
+        return res;
+        };
+
+      if(good[i].size() && bad[i].size() && good[i].back() >= bad[i][0]) {
+        fast_loglik_cont(logistics[i], fun, nullptr, 1, 1e-6);
+        total_loglik += -fun(logistics[i]);
+        // println(hlog, i, ": R=", logistics[i].R, " T=", logistics[i].T, " f = ", -fun(logistics[i]));
+        }
+      }
+    println(hlog, "total asymmetric loglikelihood = ", total_loglik);
+    }
+
+  auto llcont_eps = [] (ld eps) {
+    return [eps] (logistic& cont) {
+      ld res = 0;
+      int DS = isize(disttable_approx);
+      for(int i=0; i<DS; i++) {
+        ld d = (i + .5) / llcont_approx_prec;
+        if(disttable_approx[i][0])
+          res += (cont.lno(d - eps) + cont.lno(d + eps)) / 2 * disttable_approx[i][0];
+        if(disttable_approx[i][1])
+          res += (cont.lyes(d - eps) + cont.lyes(d + eps)) / 2 * disttable_approx[i][1];
+        }
+      return res;
+      };
+    };
+  auto compute_for_epsilon = [&] (ld eps, bool report) {
+    ld dist_entropy = 0;
+    ld sphere_entropy = 0;
+    ld integral = alpha ? (exp(alpha*max_radius)-1) / alpha : max_radius;
+    for(auto& d: center_distances) {
+      dist_entropy += alpha ? -log( (exp(alpha*(d+eps))-exp(alpha*(d-eps))) / alpha / integral) : - log((2*eps)/integral);
+      if(d > eps) sphere_entropy += log(sinh(d) / eps) * (dim-1); // 4pi is O(1)
+      }
+    ld extra_info = 0;
+    if(symmetric) extra_info = -llcont_eps(eps)(cont);
+    if(!symmetric) {
+      for(int i=0; i<N; i++) {
+        auto last = logistics[i];
+        auto fun = [&] (logistic& l) {
+          ld res = 0;
+          for(auto& g: good[i]) res += l.lyes(g+eps) + l.lyes(g-eps);
+          for(auto& b: bad[i]) res += l.lno(b+eps) + l.lno(b-eps);
+          return res/2;
+          };
+        if(good[i].size() && bad[i].size() && good[i].back() >= bad[i][0] - 2 * eps) {
+          fast_loglik_cont(last, fun, nullptr, 0.1, 1e-4);
+          extra_info -= fun(last);
+          }
+        }
+      }
+    if(report) {
+      println(hlog, "FOR epsilon = ", eps);
+      println(hlog, "  dist_entropy = ", dist_entropy);
+      println(hlog, "  sphere_entropy = ", sphere_entropy);
+      println(hlog, "  negloglik = ", extra_info);
+      println(hlog, "  control = ", format("%.1f", dist_entropy + sphere_entropy + extra_info));
+      }
+    return dist_entropy + sphere_entropy + extra_info;
+    };
+
+  ld last_rep = HUGE_VAL;
+  bool rep = !symmetric;
+  ld step = symmetric ? 0.1 : 1;
+  for(double eps=step; eps<=10; eps+=(symmetric ? 0.1 : 1)) {
+    ld cur_rep = compute_for_epsilon(eps, rep);
+    if(eps >= 9) { if(!rep) compute_for_epsilon(eps, true); break; }
+    if(cur_rep > last_rep) {
+      ld eps_min = eps-2*step, eps_max = eps;
+      int iterations = 1;
+      ld last_control = 0;
+      while(iterations--) {
+        for(int it=0; it< (symmetric ? 100 : 10); it++) {
+          ld a = (eps_min * 2 + eps_max) / 3;
+          ld b = (eps_min + 2 * eps_max) / 3;
+          if(compute_for_epsilon(a, rep) < compute_for_epsilon(b, false)) eps_max = b; else eps_min = a;
+          }
+        ld control = compute_for_epsilon(eps_min, true);
+        if(abs(control - last_control) < .1) break;
+        last_control = control;
+        if(iterations) fast_loglik_cont(cont, llcont_eps(eps_min), "lcont_eps", 1, 1e-6);  
+        eps_min = eps_min - 1; eps_max = eps_max + 1;
+        }
+      break;
+      }
+    last_rep = cur_rep;
+    }
+  }
+
+}
+}