rogueviz:: added DHRG

This commit is contained in:
Zeno Rogue 2022-07-12 11:12:35 +02:00
parent e04a6bc49a
commit 9a916f1ef3
19 changed files with 2994 additions and 0 deletions

View File

@ -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);
}
}

View File

@ -0,0 +1,2 @@
#define BOXSIZE 128
#include "dhrg.cpp"

View File

@ -0,0 +1,2 @@
#define BOXSIZE 64
#include "dhrg.cpp"

View File

@ -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)++; }
};
}

276
rogueviz/dhrg/dhrg.cpp Normal file
View File

@ -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
}

146
rogueviz/dhrg/dynamic.cpp Normal file
View File

@ -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()));
}
}

340
rogueviz/dhrg/embedder.cpp Normal file
View File

@ -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();
}
}

View File

@ -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;
}
}
}

View File

@ -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());
}
}

434
rogueviz/dhrg/loglik.cpp Normal file
View File

@ -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);
}
}
}

284
rogueviz/dhrg/mycell.cpp Normal file
View File

@ -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
}

55
rogueviz/dhrg/paths.cpp Normal file
View File

@ -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);
}
}

View File

@ -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();
}
}

113
rogueviz/dhrg/regular.cpp Normal file
View File

@ -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();
}
}

150
rogueviz/dhrg/routing.cpp Normal file
View File

@ -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());
}
}

261
rogueviz/dhrg/segment.cpp Normal file
View File

@ -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);
}
}

91
rogueviz/dhrg/tests.cpp Normal file
View File

@ -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);
}
}

116
rogueviz/dhrg/visualize.cpp Normal file
View File

@ -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;
}

View File

@ -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