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