DHRG not exports some things

rogueviz/dhrg/dhrg-utils.cpp

@@ -3,8 +3,6 @@
 
 namespace dhrg {
 
-typedef long long ll;
-
 struct progressbar : indenter_finish {
   string name;
   static const int PBSIZE = 64;

rogueviz/dhrg/dhrg.cpp

@@ -1,41 +1,9 @@
 // see this paper: https://arxiv.org/abs/2109.11772
 
-#define DHRGVER "7.1"
+#include "dhrg.h"
-#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"
@@ -53,6 +21,10 @@ bool dhrg_animate(int sym, int uni);
 
 namespace dhrg {
 
+int M;
+vector<struct mycell*> vertices;
+vector<ld> disttable0, disttable1;
+
 void memoryInfo() {
   string s = "";
   

rogueviz/dhrg/dhrg.h

@@ -0,0 +1,76 @@
+#define DHRGVER "7.1"
+#include "../rogueviz.h"
+
+#define LONG_BRACKETS
+
+namespace dhrg {
+
+using namespace hr;
+
+typedef long long ll;
+
+#ifndef BOXSIZE
+static const int BOXSIZE = 32;
+#endif
+static const int MAXDIST = (2*BOXSIZE);
+static const int SETS = 4;
+
+struct segment;
+
+inline cell *croot() { return currentmap->gamestart(); }
+
+extern int M;
+extern vector<struct mycell*> vertices;
+
+extern vector<ld> disttable0, disttable1;
+
+void memoryInfo();
+
+void cellcoords();
+void origcoords();
+void build_disttable();
+
+void dhrg_init();
+bool dhrg_animate(int sym, int uni);
+
+/* implemented in loglik.cpp: */
+
+/* for logistic regression */
+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)); }
+  logistic() {}
+  logistic(ld _R, ld _T) : R(_R), T(_T) {}  
+  void setRT(ld _R, ld _T) { R = _R; T = _T; }
+  };
+
+extern ld llcont_approx_prec;
+extern vector<array<ll, 2>> disttable_approx;
+
+using logisticfun = std::function<ld(logistic&)>;
+
+extern logistic current_logistic;
+
+ld loglik_cont_approx(logistic& l = current_logistic);
+
+void fast_loglik_cont(logistic& l, const logisticfun& f, const char *name, ld start, ld eps);
+
+/* greedy routing */
+
+struct iddata {
+  ld tot, suc, routedist, bestdist;
+  iddata() { tot = suc = routedist = bestdist = 0; }
+  };
+
+using neighborhoodfun = std::function<vector<int> (int)>;
+using distfun = std::function<ld(int a, int b)>;
+
+void prepare_pairs(int N, const neighborhoodfun& nei);
+void greedy_routing(iddata& d, const distfun& distance_function);
+
+}
+
+

rogueviz/dhrg/loglik.cpp

@@ -2,10 +2,11 @@
 
 #include <thread>
 #define USE_THREADS
-int threads = 32;
 
 namespace dhrg {
 
+int threads = 32;
+
 ld llcont_approx_prec = 10000;
 
 // tally edges of the given vertex at the given index
@@ -74,15 +75,6 @@ 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);
@@ -379,7 +371,7 @@ void build_disttable_approx() {
         disttable_approx[i][j] += r[i][j];
   }
 
-ld loglik_cont_approx(logistic& l = current_logistic) {
+ld loglik_cont_approx(logistic& l) {
     
   ld llh = 0;
   int N = isize(disttable_approx);
@@ -389,46 +381,50 @@ ld loglik_cont_approx(logistic& l = current_logistic) {
     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) {
+using logisticfun = std::function<ld(logistic&)>;
 
-  indenter_finish im("fix_logistic_parameters");
+void fast_loglik_cont(logistic& l, const logisticfun& f, const char *name, ld start, ld eps) {
+
+  if(name) println(hlog, "fix_logistic_parameters");
+  indenter_finish im(name);
   ld cur = f(l);
-  println(hlog, format("%s = %20.10" PLDF " (R=%10.5" PLDF " T=%" PLDF ")\n", name, cur, l.R, l.T));
+  if(name) println(hlog, format("%s = %20.10" PLDF " (R=%10.5" PLDF " T=%" PLDF ")", name, cur, l.R, l.T));
 
   map<pair<double, double>, double> memo;
   auto ff = [&] () {
+    if(l.T < -5) exit(1);
     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);
     };
   
+  int steps = 0;
+
   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; }
+    while(true) { steps++; l.R += step; ld t = ff(); if(t <= cur || steps > 1000) break; cur = t; changed = true; }
     l.R -= step;
     
-    while(true) { l.R -= step; ld t = ff(); if(t <= cur) break; cur = t; changed = true; }
+    while(true) { steps++; l.R -= step; ld t = ff(); if(t <= cur || steps > 1000) break; cur = t; changed = true; }
     l.R += step;
 
-    while(true) { l.T += step; ld t = ff(); if(t <= cur) break; cur = t; changed = true; }
+    while(true) { steps++; l.T += step; ld t = ff(); if(t <= cur || steps > 1000) break; cur = t; changed = true; }
     l.T -= step;
     
-    while(true) { l.T -= step; ld t = ff(); if(t <= cur) break; cur = t; changed = true; }
+    while(true) { steps++; l.T -= step; ld t = ff(); if(t <= cur || l.T < 1e-3 || steps > 1000) 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));
+    if(name) println(hlog, format("%s = %20.10" PLDF " (R=%10.5" PLDF " T=%10.5" PLDF ")", name, cur, l.R, l.T));
     fflush(stdout);
     }
   }
 
-
 }

rogueviz/dhrg/routing.cpp

@@ -14,7 +14,12 @@ vector<pairdata> pairs;
 vector<int> last_goal;
 vector<int> next_stop;
 
-void prepare_pairs() {
+int gr_N;
+
+using neighborhoodfun = std::function<vector<int> (int)>;
+
+void prepare_pairs(int N, const neighborhoodfun& nei) {
+  gr_N = N;
   pairs.resize(N);
   actual.resize(N);
   for(int i=0; i<N; i++) actual[i].resize(N, NOYET);
@@ -30,18 +35,27 @@ void prepare_pairs() {
     visit(i, 0);
     for(int k=0; k<isize(bfsqueue); k++) {
       int a = bfsqueue[k];
-      for(auto ed: rogueviz::vdata[a].edges) {
+      for(auto b: nei(a))
-        int b = ed.second->i ^ ed.second->j ^ a;
         visit(b, p[a] + 1);
       }
     }
-    }
   last_goal.clear();
   last_goal.resize(N, -1);
   next_stop.clear();
   next_stop.resize(N, -1);
   }
 
+void prepare_pairs() {
+  prepare_pairs(N, [] (int a) {
+    vector<int> res;
+    for(auto ed: rogueviz::vdata[a].edges) {
+      int b = ed.second->i ^ ed.second->j ^ a;
+      res.push_back(b);
+      }
+    return res;
+    });
+  }
+
 void route_from(int src, int goal, const vector<ld>& distances_from_goal) {
   if(last_goal[src] == goal) return;
   if(src == goal) {
@@ -78,21 +92,14 @@ void route_from(int src, int goal, const vector<ld>& distances_from_goal) {
   last_goal[src] = goal;
   }
 
-struct iddata {
+void greedy_routing_to(iddata& d, int goal, const distfun& distance_function) {
-  ld tot, suc, routedist, bestdist;
+  vector<ld> distances_from_goal(gr_N);
-  iddata() { tot = suc = routedist = bestdist = 0; }
+  for(int src=0; src<gr_N; src++)
-  } datas[6];
-
-template<class T> void greedy_routing_to(int id, int goal, const T& distance_function) {
-  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++)
+  for(int src=0; src<gr_N; src++)
     route_from(src, goal, distances_from_goal);
 
-  auto& d = datas[id];
+  for(int j=0; j<gr_N; j++) if(j != goal){
-
-  for(int j=0; j<N; j++) if(j != goal){
     d.tot++;
     ld p = pairs[j].success;
     d.suc += p;
@@ -101,12 +108,13 @@ template<class T> void greedy_routing_to(int id, int goal, const T& distance_fun
     }
   }
 
-template<class T> void greedy_routing(int id, const T& distance_function) {
+void greedy_routing(iddata& d, const distfun& distance_function) {
-  for(int goal=0; goal<N; goal++) greedy_routing_to(id, goal, distance_function);
+  for(int goal=0; goal<gr_N; goal++) greedy_routing_to(d, goal, distance_function);
   }
 
 void routing_test(string s) {
 
+  iddata datas[6];
   vector<string> reps;
   vector<ld> values;
 
@@ -131,17 +139,17 @@ void routing_test(string s) {
     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(datas[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); });
+    greedy_routing(datas[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); });
+    greedy_routing(datas[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(datas[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); });
+    greedy_routing(datas[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}) {
@@ -157,8 +165,9 @@ void routing_test(string s) {
   }
 
 int current_goal;
+iddata prepared;
 void prepare_goal(int goal) {
-  greedy_routing_to(0, current_goal = goal, [] (int i, int j) { return hdist(vertexcoords[i], vertexcoords[j]); });
+  greedy_routing_to(prepared, current_goal = goal, [] (int i, int j) { return hdist(vertexcoords[i], vertexcoords[j]); });
   }
 
 vector<int> path(int src) {