rogueviz:: SAG improvements

rogueviz/sag.cpp

@@ -7,14 +7,23 @@
 // see: https://www.youtube.com/watch?v=WSyygk_3j9o (SAG roguelikes)
 // see: https://www.youtube.com/watch?v=HWQkDkeEUeM (SAG programming languages)
 
+#include "dhrg/dhrg.h"
+#include <thread>
+
 namespace rogueviz {
 
 namespace sag {
 
+  int threads = 1;
+
+  int informat; /* format ID */
+
   bool turn(int delta);
   
   int sagpar = 0;
 
+  int best_cost = 1000000000;
+
   enum eSagmode { sagOff, sagHC, sagSA };
   
   eSagmode sagmode; // 0 - off, 1 - hillclimbing, 2 - SA
@@ -24,6 +33,10 @@ namespace sag {
   ld temperature = -4;
   const char *loadfname;
   
+  string auto_save;
+
+  bool auto_visualize = true;
+
   int vizsa_start;
   int vizsa_len;
   
@@ -45,6 +58,38 @@ namespace sag {
   /** if i in neighbors[j], sagcells[i] is a neighbor of sagcells[j] */
   vector<vector<int>> neighbors;
   
+  ld pdist(hyperpoint hi, hyperpoint hj);
+  
+  /** matrix for every sagcell */
+  vector<transmatrix> cell_matrix;
+  
+  /** precision of geometric distances */
+  int gdist_prec;
+  
+  /** the maximum value in sagdist +1 */
+  int max_sag_dist;
+
+  vector<edgeinfo> sagedges;  
+  vector<vector<int>> edges_yes, edges_no;
+
+  int logistic_cost; /* 0 = disable, 1 = enable */
+  dhrg::logistic lgsag(1, 1);
+  
+  vector<ld> loglik_tab_y, loglik_tab_n;
+
+  int ipturn = 100;
+  int numiter = 0;
+  
+  int hightemp = 10;
+  int lowtemp = -15;
+  
+  /* for the embedding method: */  
+  bool embedding;
+  dhrg::logistic lgemb(1, 1);
+  vector<hyperpoint> placement;  
+
+  void optimize_sag_loglik();
+
   void compute_dists() {
     int N = isize(sagcells);
 
@@ -55,16 +100,49 @@ namespace sag {
       for(cell *c1: adj_minefield_cells(sagcells[i]))
         if(ids.count(c1)) neighbors[i].push_back(ids[c1]);
 
-    sagdist.clear();
+    const ld ERROR = -17.3;
-    sagdist.resize(N);
+    transmatrix unknown = Id; unknown[0][0] = ERROR;
-    for(int i=0; i<N; i++) {
+    cell_matrix.clear();
-      auto &sdi = sagdist[i];
+    cell_matrix.resize(N, unknown);
-      sdi.resize(N, N);
+    vector<int> visited;
-      vector<int> q;
+
-      auto visit = [&] (int j, int dist) { if(sdi[j] < N) return; sdi[j] = dist; q.push_back(j); };
+    auto visit = [&] (int id, const transmatrix& T) {
-      visit(i, 0);
+      if(cell_matrix[id][0][0] != ERROR) return;
-      for(int j=0; j<isize(q); j++) for(int k: neighbors[q[j]]) visit(k, sdi[q[j]]+1);
+      cell_matrix[id] = T;
+      visited.push_back(id);
+      };
+      
+    visit(0, Id);
+    for(int i=0; i<isize(visited); i++) {
+      cell *c0 = sagcells[i];
+      const transmatrix& T0 = cell_matrix[i];
+      for(int d=0; d<c0->type; d++)
+        if(ids.count(c0->move(d)))
+          visit(ids[c0->move(d)], T0 * currentmap->adj(c0, d));
       }
+    
+    if(gdist_prec) {
+      for(int i=0; i<N; i++)
+      for(int j=0; j<N; j++)
+        sagdist[i][j] = (pdist(tC0(cell_matrix[i]), tC0(cell_matrix[j])) + .5) * gdist_prec;
+      }
+    
+    else {
+      sagdist.clear();
+      sagdist.resize(N);
+      for(int i=0; i<N; i++) {
+        auto &sdi = sagdist[i];
+        sdi.resize(N, N);
+        vector<int> q;
+        auto visit = [&] (int j, int dist) { if(sdi[j] < N) return; sdi[j] = dist; q.push_back(j); };
+        visit(i, 0);
+        for(int j=0; j<isize(q); j++) for(int k: neighbors[q[j]]) visit(k, sdi[q[j]]+1);
+        }
+      }
+    
+    max_sag_dist = 0;
+    for(auto& d: sagdist) for(auto& x: d) max_sag_dist = max(max_sag_dist, x);
+    max_sag_dist++;
     }
 
   bool legacy;
@@ -107,6 +185,16 @@ namespace sag {
   double costat(int vid, int sid) {
     if(vid < 0) return 0;
     double cost = 0;
+    
+    if(logistic_cost) {
+      auto &s = sagdist[sid];
+      for(auto j: edges_yes[vid])
+        cost += loglik_tab_y[s[sagid[j]]];
+      for(auto j: edges_no[vid])
+        cost += loglik_tab_n[s[sagid[j]]];
+      return -cost;
+      }
+    
     vertexdata& vd = vdata[vid];
     for(int j=0; j<isize(vd.edges); j++) {
       edgeinfo *ei = vd.edges[j].second;
@@ -178,18 +266,37 @@ namespace sag {
 
     sagnode[sid1] = t2; sagnode[sid2] = t1;
     sagid[t1] = sid2; if(t2 >= 0) sagid[t2] = sid1;
-    cost += 2*change;
+    cost += change;
     }
   
-  void create_sagnode() {
+  void prepare_graph() {
+    int DN = isize(sagid);
+
+    set<pair<int, int>> alledges;
+    for(auto e: sagedges) {
+      if(e.i == e.j) continue;
+      alledges.emplace(e.i, e.j);
+      alledges.emplace(e.j, e.i);
+      }
+    
+    edges_yes.clear(); edges_yes.resize(DN);
+    edges_no.clear(); edges_no.resize(DN);
+    
+    for(int i=0; i<DN; i++) for(int j=0; j<DN; j++) if(i != j) {
+      if(alledges.count({i, j}))
+        edges_yes[i].push_back(j);
+      else
+        edges_no[i].push_back(j);
+      }          
+
     sagnode.clear();
     sagnode.resize(isize(sagcells), -1);
-    int DN = isize(sagid);
     for(int i=0; i<DN; i++)
       sagnode[sagid[i]] = i;
     cost = 0;
     for(int i=0; i<DN; i++)
       cost += costat(i, sagid[i]);
+    cost /= 2;
     }
   
   void reassign() {
@@ -229,41 +336,33 @@ namespace sag {
     afterload: 
     if(sf) fclose(sf);
 
-    create_sagnode();
+    prepare_graph();
     reassign();
     }
   
-  vector<edgeinfo> sagedges;
-  
-  int ipturn = 100;
-  int numiter = 0;
-  
-  int hightemp = 10;
-  int lowtemp = -15;
-  
   void dofullsa(int satime) {
     sagmode = sagSA;
     int t1 = SDL_GetTicks();
+    int tl = -999999;
     
     while(true) {
       int t2 = SDL_GetTicks();
       double d = (t2-t1) / (1000. * satime);
       if(d > 1) break;
+
       temperature = hightemp - (d*(hightemp-lowtemp));
-      chgs.clear();
+      for(int i=0; i<10000; i++) {
-      for(int i=0; i<50000; i++) {
         numiter++;
         sag::saiter();
         }
       
-      print(hlog, format("it %8d temp %6.4f [1/e at %13.6f] cost = %f ", 
+      if(t2 - tl > 980) {
-        numiter, double(sag::temperature), (double) exp(sag::temperature),
+        tl = t2;
-        double(sag::cost)));
+        println(hlog, format("it %8d temp %6.4f [1/e at %13.6f] cost = %f ", 
+          numiter, double(sag::temperature), (double) exp(sag::temperature),
+          double(sag::cost)));
+        }
       
-      sort(chgs.begin(), chgs.end());
-      int cc = chgs.size() - 1;
-      println(hlog, format("%9.4f .. %9.4f .. %9.4f .. %9.4f .. %9.4f", 
-        double(chgs[0]), double(chgs[cc/4]), double(chgs[cc/2]), double(chgs[cc*3/4]), double(chgs[cc])));
       }
     
     temperature = -5;
@@ -290,7 +389,7 @@ namespace sag {
       (double) exp(-50 * exp(-sag::temperature)),
       (double) sag::cost));
 
-    reassign();
+    if(auto_visualize) reassign();
     }
   
   void save_sag_solution(const string& fname) {
@@ -300,36 +399,324 @@ namespace sag {
     fclose(f);
     }
 
-  void loglik() {
+  void compute_loglik_tab() {
-    int indist[30], pedge[30];
+    loglik_tab_y.resize(max_sag_dist);
-    for(int d=0; d<30; d++) indist[d] = 0, pedge[d] = 0;
+    loglik_tab_n.resize(max_sag_dist);
+    for(int i=0; i<max_sag_dist; i++) {
+      loglik_tab_y[i] = lgsag.lyes(i);
+      loglik_tab_n[i] = lgsag.lno(i);
+      }
+    }
+
+  void optimize_sag_loglik() {
+    vector<int> indist(max_sag_dist, 0);
+    
+    const int mul = 1;
 
     int N = isize(sagid);
     for(int i=0; i<N; i++)
-    for(int j=0; j<i; j++)
+    for(int j=0; j<i; j++) {
-      indist[sagdist[sagid[i]][sagid[j]]]++;
+      int d = sagdist[sagid[i]][sagid[j]];
+      indist[d]++;
+      }
+    
+    vector<int> pedge(max_sag_dist, 0);
       
     for(int i=0; i<isize(sagedges); i++) {
       edgeinfo& ei = sagedges[i];
+      if(int(sagdist[sagid[ei.i]][sagid[ei.j]] * mul) == 136) printf("E %d,%d\n", ei.i, ei.j);
       if(ei.i != ei.j)
       if(ei.weight >= sag_edge->visible_from)
-        pedge[sagdist[sagid[ei.i]][sagid[ei.j]]]++;
+        pedge[sagdist[sagid[ei.i]][sagid[ei.j]] * mul]++;
       }
     
-    for(int d=0; d<30; d++) 
+    for(int d=0; d<max_sag_dist; d++) 
       if(indist[d])
         printf("%2d: %7d/%7d %7.3lf\n", 
           d, pedge[d], indist[d], double(pedge[d] * 100. / indist[d]));
         
     ld loglik = 0;
-    for(int d=0; d<30; d++) {
+    for(int d=0; d<max_sag_dist; d++) {
       int p = pedge[d], pq = indist[d];
       int q = pq - p;
-      if(p && q)
+      if(p && q) {
         loglik += p * log(p) + q * log(q) - pq * log(pq);
+        println(hlog, tie(d, p, q), loglik);
+        }
       }
     
-    println(hlog, "loglikelihood = ", fts(loglik));
+    println(hlog, "loglikelihood best = ", fts(loglik));
+    
+    auto logisticf = [&] (dhrg::logistic&  l) {
+      ld loglik = 0;
+      for(int d=0; d<max_sag_dist; d++) {
+        int p = pedge[d], pq = indist[d];
+        if(p) loglik += p * l.lyes(d);
+        if(pq > p) loglik += (pq-p) * l.lno(d);
+        }
+      return loglik;
+      };
+    
+    dhrg::fast_loglik_cont(lgsag, logisticf, nullptr, 1, 1e-5);
+    println(hlog, "loglikelihood logistic = ", logisticf(lgsag), " R= ", lgsag.R, " T= ", lgsag.T);    
+    
+    if(logistic_cost) {
+      compute_loglik_tab();
+      prepare_graph();       
+      println(hlog, "cost = ", cost);
+      }
+    }
+
+  void disttable_add(ld dist, int qty0, int qty1) {
+    using namespace dhrg;
+    size_t i = dist * llcont_approx_prec;
+    constexpr array<ll, 2> zero = {0, 0};
+    while(disttable_approx.size() <= i) disttable_approx.push_back(zero);
+    disttable_approx[i][0] += qty0;
+    disttable_approx[i][1] += qty1;
+    }
+  
+  ld approx_01(hyperpoint h) {
+    ld d = 0;
+    if(h[0] > 1) {
+      ld z = log(h[0]);
+      d += z; h[1] *= h[0]; h[0] = 1; h[2] += z;
+      }
+    d += h[0];
+    if(h[1] > 1) {
+      ld z = log(h[1]);
+      d += z; h[1] = 1; h[2] -= z;
+      }
+    d += h[1];
+    d += abs(h[2]);
+    return d;
+    }
+
+  ld pdist(hyperpoint hi, hyperpoint hj) {
+    if(sol) {
+      hyperpoint h = rgpushxto0(hi) * hj;
+      
+      h[0] = abs(h[0]);
+      h[1] = abs(h[1]);
+      
+      ld d1 = approx_01(h);
+      ld d2 = approx_01(hyperpoint(h[1], h[0], -h[2], 1));
+
+      h = rgpushxto0(hj) * hi;
+      
+      h[0] = abs(h[0]);
+      h[1] = abs(h[1]);
+      
+      ld d3 = approx_01(h);
+      ld d4 = approx_01(hyperpoint(h[1], h[0], -h[2], 1));
+
+      return min(min(d1, d2), min(d3, d4));
+      }
+    return geo_dist(hi, hj);
+    };
+
+  ld pdist(int i, int j) {
+    return pdist(placement[i], placement[j]);
+    };
+
+  void prepare_embedding() {
+    map<int, transmatrix> maps;    
+    vector<int> visited;
+
+    auto visit = [&] (int id, const transmatrix& T) {
+      if(maps.count(id)) return;
+      maps[id] = T;
+      visited.push_back(id);
+      };
+      
+    visit(0, Id);
+    for(int i=0; i<isize(visited); i++) {
+      cell *c0 = sagcells[i];
+      transmatrix T0 = maps[i];
+      for(int d=0; d<c0->type; d++)
+        if(ids.count(c0->move(d)))
+          visit(ids[c0->move(d)], T0 * currentmap->adj(c0, d));
+      }
+    
+    int DN = isize(sagid);
+    placement.resize(DN);
+    for(int i=0; i<DN; i++) placement[i] = tC0(maps[sagid[i]]);    
+    }
+  
+  int embiter;
+  
+  void compute_loglik() {    
+    dhrg::llcont_approx_prec = 10;
+    
+    dhrg::disttable_approx.clear();
+    int DN = isize(sagid);
+    for(int i=0; i<DN; i++)
+    for(int j=0; j<i; j++)
+      disttable_add(pdist(i, j), 1, 0);
+
+    for(int i=0; i<isize(sagedges); i++) {
+      edgeinfo& ei = sagedges[i];
+      if(ei.i != ei.j)
+        disttable_add(pdist(ei.i, ei.j), -1, 1);
+      }
+    
+    dhrg::logisticfun lc = dhrg::loglik_cont_approx;
+
+    dhrg::fast_loglik_cont(lgemb, lc, nullptr, 1, 1e-5);
+    
+    println(hlog, "loglik = ", format("%.6f", lc(lgemb)), " R = ", lgemb.R, " T = ", lgemb.T, " iterations = ", embiter);
+    }
+
+  void reassign_embedding() {
+    int DN = isize(sagid);
+    for(int i=0; i<DN; i++) {
+      vdata[i].m->base = sagcells[0];
+      vdata[i].m->at = rgpushxto0(placement[i]);
+      virtualRebase(vdata[i].m);
+      forgetedges(i);
+      }
+    shmup::fixStorage();
+    }
+
+  void improve_embedding() {
+    embiter++;
+    if(placement.empty()) {
+      prepare_embedding();
+      compute_loglik();
+      }
+    ld eps = .1;
+    int DN = isize(sagid);
+
+    hyperpoint h = C0;
+    for(int i=0; i<WDIM; i++) h[i] += (hrandf() - 0.5) * eps;
+    h = normalize(h);
+    
+    auto nplacement = placement;
+    parallelize(DN, [&] (int a, int b) {
+      for(int i=a; i<b; i++) {
+
+        hyperpoint np = rgpushxto0(placement[i]) * h;
+        
+        ld change;
+        for(auto e: edges_yes[i]) change -= lgemb.lyes(pdist(placement[i], placement[e]));
+        for(auto e: edges_no[i]) change -= lgemb.lno(pdist(placement[i], placement[e]));
+        for(auto e: edges_yes[i]) change += lgemb.lyes(pdist(np, placement[e]));
+        for(auto e: edges_no[i]) change += lgemb.lno(pdist(np, placement[e]));
+  
+        if(change > 0) nplacement[i] = np;
+        }
+      return 0;
+      });
+    
+    placement = nplacement;
+    }
+  
+  int embturn = 1;
+  void embedding_iterate() {
+    int t1 = SDL_GetTicks();
+    for(int i=0; i<embturn; i++) {
+      improve_embedding();
+      }
+    int t2 = SDL_GetTicks();
+    int t = t2 - t1;
+    if(t < 50) embturn *= 2;
+    else if(t > 200) embturn = (embturn + 1) / 2;
+    else embturn = (embturn * 100 + (t-1)) / t;
+    
+    compute_loglik();
+    if(auto_visualize) reassign_embedding();
+    }
+
+  void save_embedding(const string& fname) {
+    fhstream f(fname, "wt");
+    for(int i=0; i<isize(sagid); i++) {
+      println(f, vdata[i].name);
+      for(int d=0; d<MDIM; d++)
+        println(f, format("%.20f", placement[i][d]));
+      }
+    }
+
+  void load_embedding(const string& fname) {
+    prepare_embedding();
+    fhstream f(fname, "rt");
+    if(informat == 2) {
+      /* H2 embedding */
+      while(!feof(f.f)) {
+        string lab = scan<string>(f);
+        int id;
+        if(!labeler.count(lab)) {
+          printf("unknown vertex: %s\n", lab.c_str());
+          continue;
+          }
+        else id = getid(lab);
+        ld alpha, r;
+        if(1) {
+          dynamicval<eGeometry> g(geometry, gNormal);
+          hyperpoint h;
+          for(int d=0; d<MDIM; d++) h[d] = scan<ld>(f);
+          alpha = atan2(h);
+          r = hdist0(h);
+          println(hlog, "read ", lab, " as ", h, " which is ", tie(alpha, r));
+          }
+        placement[id] = direct_exp(cspin(0, 2, alpha) * ctangent(0, r));
+        println(hlog, "dist = ", pdist(placement[id], C0), " expected: ", r);
+        }
+      }
+    else if(informat == 3) {
+      /* BFKL */
+      string ignore;
+      if(!scan(f, ignore, ignore, ignore, ignore, ignore, ignore, ignore, ignore)) {
+        printf("Error: incorrect format of the first line\n"); exit(1);
+        }      
+      while(true) {
+        string lab = scan<string>(f);
+        if(lab == "" || lab == "#ROGUEVIZ_ENDOFDATA") break;
+        ld r, alpha;
+        if(!scan(f, r, alpha)) { printf("Error: incorrect format of r/alpha\n"); exit(1); }
+        hyperpoint h = spin(alpha * degree) * xpush0(r);
+  
+        if(!labeler.count(lab)) {
+          printf("unknown vertex: %s\n", lab.c_str());
+          }
+        else {
+          int id = getid(lab);
+          placement[id] = h;
+          }      
+        }
+      }
+    else if(informat == 4) {
+      while(true) {
+        string lab = scan<string>(f);
+        if(lab == "") break;
+        ld r, alpha;
+        if(!scan(f, r, alpha)) { printf("Error: incorrect format of r/alpha\n"); exit(1); }
+        hyperpoint h = spin(alpha) * xpush0(r);
+  
+        if(!labeler.count(lab)) {
+          printf("unknown vertex: %s\n", lab.c_str());
+          }
+        else {
+          int id = getid(lab);
+          placement[id] = h;
+          }      
+        }
+      }
+    else {
+      while(!feof(f.f)) {
+        string lab = scan<string>(f);
+        int id;
+        if(!labeler.count(lab)) {
+          printf("unknown vertex: %s\n", lab.c_str());
+          continue;
+          }
+        else id = getid(lab);
+        hyperpoint h;
+        for(int d=0; d<MDIM; d++) h[d] = scan<ld>(f);
+        placement[id] = h;
+        }
+      }
+    reassign_embedding();
+    compute_loglik();
     }
   
   void read_hubs(const string& fname) {
@@ -369,6 +756,30 @@ namespace sag {
       printf("Failed to open SAG file: %s\n", fname);
       throw "failed to open SAG file";
       }
+    if(informat == 1) {
+      scanline(f);
+      set<pair<int, int> > edges;
+      
+      int all = 0, good = 0;
+      while(!feof(f.f)) {        
+        string l1 = scan<string>(f);
+        string l2 = scan<string>(f);
+        if(l1 == "") continue;
+        if(l2 == "") continue;
+        edgeinfo ei(sag_edge);
+        ei.i = getid(l1);
+        ei.j = getid(l2);
+        if(ei.i > ei.j) swap(ei.i, ei.j);
+        all++;
+        if(edges.count({ei.i, ei.j})) continue;
+        good++;
+        edges.emplace(ei.i, ei.j);
+        ei.weight = 1;
+        sagedges.push_back(ei);
+        }
+      println(hlog, "N = ", isize(vdata), " edges = ", good, "/", all);
+      return;
+      }
     while(!feof(f.f)) {
       string l1, l2;
       while(true) {
@@ -410,6 +821,31 @@ namespace sag {
         });
       dialog::addSelItem(XLAT("SAG mode"), sag::sagmodes[sag::sagmode], 'm');      
       dialog::add_action([] { sag::sagmode = sag::eSagmode( (1+sag::sagmode) % 3 ); });
+
+      dialog::addSelItem(XLAT("min temperature"), fts(sag::lowtemp), 'i');
+      dialog::add_action([] {
+        dialog::editNumber(sag::lowtemp, -20, 20, 1, 0, XLAT("min temperature"), "");
+        });
+
+      dialog::addSelItem(XLAT("max temperature"), fts(sag::hightemp), 'i');
+      dialog::add_action([] {
+        dialog::editNumber(sag::hightemp, -20, 20, 1, 0, XLAT("high temperature"), "");
+        });
+
+      dialog::addSelItem(XLAT("automatic cycle"), fts(sag::vizsa_len), 'c');
+      dialog::add_action([] {
+        dialog::editNumber(sag::vizsa_len, 5, 1800, 1, 0, XLAT("automatic cycle"), "");
+        });
+      
+      dialog::addBoolItem(XLAT("automatic"), sag::vizsa_start, 'a');
+      dialog::add_action([] { 
+        sag::vizsa_start = sag::vizsa_start ? 0 : SDL_GetTicks();
+        sag::sagmode = sagOff;
+        });
+
+      dialog::addBoolItem_action(XLAT("auto-visualize"), sag::auto_visualize, 'b');
+
+      dialog::addBoolItem_action(XLAT("continuous embedding"), sag::embedding, 'e'); 
       });
 
     weight_label = "min weight";
@@ -444,7 +880,7 @@ namespace sag {
       }
     sagid.resize(DN);
     for(int i=0; i<DN; i++) sagid[i] = i;
-    create_sagnode();
+    prepare_graph();
 
     for(int i=0; i<DN; i++) {
       int ii = i;
@@ -456,6 +892,47 @@ namespace sag {
     storeall();
     }
 
+ld compute_mAP() {
+  ld mAP = 0;
+  int DN = isize(sagid);
+  for(int i=0; i<DN; i++) {
+    vector<int> alldist;
+    for(int j=0; j<DN; j++) if(i != j) alldist.push_back(sagdist[sagid[i]][sagid[j]]);
+    sort(alldist.begin(), alldist.end());
+    int q = isize(edges_yes[i]);
+    int qmin = q-1, qmax = q+1;
+    int threshold = alldist[q-1];
+    while(qmin && alldist[qmin-1] == threshold) qmin--;
+    while(qmax < isize(alldist)-2 && alldist[qmax+1] == threshold) qmax++;
+    ld on_threshold = (q - qmin) / (qmax + 1. - qmin);
+    int good = 0, onthr = 0;
+    for(auto j: edges_yes[i]) {
+      int d = sagdist[sagid[i]][sagid[j]];
+      if(d < threshold) good++;
+      if(d == threshold) onthr++;
+      }
+    mAP += (good + onthr * on_threshold) / q / DN;
+    }
+  return mAP;
+  }
+
+int logid;
+
+void output_stats() {
+  if(auto_save != "" && cost < best_cost) {
+    println(hlog, "cost ", cost, " beats ", best_cost);
+    best_cost = cost;
+    sag::save_sag_solution(auto_save);
+    }
+  println(hlog, "solution: ", sagid);
+  int DN = isize(sagid);
+  ld mAP = compute_mAP();
+  dhrg::iddata routing_result;
+  dhrg::prepare_pairs(DN, [] (int i) { return edges_yes[i]; });
+  dhrg::greedy_routing(routing_result, [] (int i, int j) { return sagdist[sagid[i]][sagid[j]]; });
+  println(hlog, "CSV;", logid++, ";", isize(sagnode), ";", DN, ";", isize(sagedges), ";", lgsag.R, ";", lgsag.T, ";", cost, ";", mAP, ";", routing_result.suc / routing_result.tot, ";", routing_result.routedist / routing_result.bestdist);
+  }
+
 int readArgs() {
 #if CAP_COMMANDLINE
   using namespace arg;
@@ -470,9 +947,23 @@ int readArgs() {
   else if(argis("-sagminhi")) {
     shift_arg_formula(default_edgetype.visible_from_hi);
     }
+  else if(argis("-sag_gdist")) {
+    shift(); sag::gdist_prec = argi();
+    }
+  else if(argis("-sagrt")) {
+    shift(); sag::lgsag.R = argf();
+    shift(); sag::lgsag.T = argf();
+    }
+  else if(argis("-sag_use_loglik")) {  
+    shift(); sag::logistic_cost = argi();
+    if(sag::logistic_cost) compute_loglik_tab();
+    }
   else if(argis("-sagminhelp")) {
     shift_arg_formula(default_edgetype.visible_from_help);
     }
+  else if(argis("-sagformat")) {
+    shift(); informat = argi();
+    }
 
 // (1) configure edge weights
   else if(argis("-sag-edgepower")) {
@@ -494,6 +985,10 @@ int readArgs() {
     PHASE(3); 
     shift(); sag::read(args());
     }
+  else if(argis("-sagaviz")) {
+    PHASE(3); 
+    shift(); sag::auto_visualize = argi();
+    }
   else if(argis("-saghubs")) {
     println(hlog, "HUBS");
     PHASE(3); 
@@ -512,21 +1007,40 @@ int readArgs() {
     sag::vizsa_start = SDL_GetTicks();
     shift(); sag::vizsa_len = argi();
     }
+  else if(argis("-sagstats")) {
+    output_stats();
+    }
 // (5) save the positioning
   else if(argis("-sagsave")) {
     PHASE(3); shift(); sag::save_sag_solution(args());
     }
+  else if(argis("-sagsave-auto")) {
+    PHASE(3); shift(); auto_save = args();
+    }
 // (6) output loglikelihood
   else if(argis("-sagloglik")) {
-    sag::loglik();
+    sag::optimize_sag_loglik();
     }
   else if(argis("-sagmode")) {
     shift();
+    vizsa_start = 0;
     sagmode = (eSagmode) argi();
     if(sagmode == sagSA) {
       shift(); temperature = argf();
       }
     }
+  else if(argis("-sagembed")) {
+    sag::embedding = true;
+    }
+  else if(argis("-sagembedoff")) {
+    sag::embedding = false;
+    }
+  else if(argis("-sagsavee")) {
+    PHASE(3); shift(); sag::save_embedding(args());
+    }
+  else if(argis("-sagloade")) {
+    PHASE(3); shift(); sag::load_embedding(args());
+    }
   else return 1;
 #endif
   return 0;
@@ -536,12 +1050,20 @@ bool turn(int delta) {
   if(vizsa_start) {
     auto t = ticks;
     double d = (t-vizsa_start) / (1000. * vizsa_len);
+    if(d > 1 && logistic_cost == 2) {
+      vizsa_start = ticks;
+      optimize_sag_loglik();
+      output_stats();
+      }
     if(d > 1) sagmode = sagOff;
     else {
       temperature = hightemp - (d*(hightemp-lowtemp));
       sagmode = sagSA;
       }
     }
+  if(sagmode == sagOff && embedding) {
+    embedding_iterate();
+    }
   iterate();
   return false;
   // shmup::pc[0]->rebase();