#include "../hyper.h" // This program generates the error table for Solv approxiations. #define D3 1 #define D2 0 #if CAP_FIELD namespace hr { ld solerror(hyperpoint ok, hyperpoint chk) { return geo_dist(chk, ok); } ld minz = -1e-9, maxz = 1e-9; int max_iter = 999999; bool isok; hyperpoint iterative_solve(hyperpoint xp, hyperpoint candidate, ld minerr, bool debug = false) { transmatrix T = Id; T[0][1] = 8; T[2][2] = 5; auto f = [&] (hyperpoint x) { return nisot::numerical_exp(x); }; // T * x; }; auto ver = f(candidate); ld err = solerror(xp, ver); auto at = candidate; ld eps = 1e-6; hyperpoint c[6]; for(int a=0; a<3; a++) c[a] = point3(a==0, a==1, a==2); for(int a=0; a<3; a++) c[3+a] = point3(-(a==0), -(a==1), -(a==2)); int iter = 0; while(err > minerr) { again: iter++; if(iter > max_iter) { isok = false; return at; } // cands.push_back(at); if(debug) println(hlog, "\n\nf(", at, "?) = ", ver, " (error ", err, ")"); array pnear; for(int a=0; a<3; a++) { auto x = at + c[a] * eps; if(debug) println(hlog, "f(", x, ") = ", f(x), " = y + ", f(x)-ver, "imp ", err - solerror(xp, f(x)) ); auto y = at - c[a] * eps; if(debug) println(hlog, "f(", y, ") = ", f(y), " = y + ", f(y)-ver, "imp ", err - solerror(xp, f(y)) ); pnear[a] = (f(x) - ver) / eps; // (direct_exp(at + c[a] * eps, prec) - ver) / eps; } transmatrix U = Id; for(int a=0; a<3; a++) for(int b=0; b<3; b++) U[a][b] = pnear[b][a]; hyperpoint diff = (xp - ver); hyperpoint bonus = inverse(U) * diff; ld lbonus = hypot_d(3, bonus); if(lbonus > 0.1) bonus = bonus * 0.1 / hypot_d(3, bonus); if(false && lbonus > 1000) { int best = -1; ld besti = err; for(int a=0; a<6; a++) { auto x = at + c[a] * eps; auto nerr = solerror(xp, f(x)); if(nerr < besti) best = a, besti = nerr; } if(best == -1) { println(hlog, "local best"); for(int a=0; a<1000000; a++) { auto x = at; for(int i=0; i<3; i++) x[i] += (hrand(1000000) - hrand(1000000)) * 1e-5; auto nerr = solerror(xp, f(x)); if(nerr < besti) { println(hlog, "moved to ", x); at = x; goto again; } } break; } bonus = c[best] * 1e-3; } int fixes = 0; if(debug) println(hlog, "\nU = ", U, "\ndiff = ", diff, "\nbonus = ", bonus, " of ", lbonus, "\n"); nextfix: hyperpoint next = at + bonus; hyperpoint nextver = f(next); ld nexterr = solerror(xp, nextver); if(debug) println(hlog, "f(", next, ") = ", nextver, ", imp = ", err - nexterr); if(nexterr < err) { // println(hlog, "reduced error ", err, " to ", nexterr); at = next; ver = nextver; err = nexterr; continue; } else { bonus /= 2; fixes++; if(fixes > 10) { if(err > 999) { for(ld s = 1; abs(s) > 1e-9; s *= 0.5) for(int k=0; k<27; k++) { int kk = k; next = at; for(int i=0; i<3; i++) { if(kk%3 == 1) next[i] += s; if(kk%3 == 2) next[i] -= s; kk /= 3; } // next = at + c[k] * s; nextver = f(next); nexterr = solerror(xp, nextver); // println(hlog, "f(", next, ") = ", nextver, ", error = ", nexterr); if(nexterr < err) { at = next; ver = nextver; err = nexterr; goto nextiter; } } println(hlog, "cannot improve error ", err); exit(1); } if(debug) println(hlog, "fixes = ", fixes, " : break"); isok = false; return at; } goto nextfix; } nextiter: ; } if(debug) println(hlog, "\n\nsolution found: f(", at, ") = ", ver, " (error ", err, ")"); isok = true; return at; } EX void geodesic_step_euler(hyperpoint& at, hyperpoint& velocity) { auto acc = nisot::christoffel(at, velocity, velocity); at = at + velocity + acc / 2; velocity += acc; } EX void geodesic_step_poor(hyperpoint& at, hyperpoint& velocity) { auto acc = nisot::christoffel(at, velocity, velocity); at = at + velocity; velocity += acc; } EX void geodesic_step_midpoint(hyperpoint& at, hyperpoint& velocity) { // y(n+1) = y(n) + f(y(n) + 1/2 f(y(n))) auto acc = nisot::christoffel(at, velocity, velocity); auto at2 = at + velocity / 2; auto velocity2 = velocity + acc / 2; auto acc2 = nisot::christoffel(at2, velocity2, velocity2); at = at + velocity + acc2 / 2; velocity = velocity + acc2; } auto& chr = nisot::get_acceleration; EX bool invalid_any(const hyperpoint h) { return isnan(h[0]) || isnan(h[1]) || isnan(h[2]) || isinf(h[0]) || isinf(h[1]) || isinf(h[2]) || abs(h[0]) > 1e20 || abs(h[1]) > 1e20 || abs(h[2]) > 1e20; } EX void geodesic_step_rk4(hyperpoint& at, hyperpoint& vel) { auto acc1 = chr(at, vel); auto acc2 = chr(at + vel/2, vel + acc1/2); auto acc3 = chr(at + vel/2 + acc1/4, vel + acc2/2); auto acc4 = chr(at + vel + acc2/2, vel + acc3); at += vel + (acc1+acc2+acc3)/6; vel += (acc1+2*acc2+2*acc3+acc4)/6; } template hyperpoint numerical_exp(hyperpoint v, int steps, const T& gstep) { hyperpoint at = point31(0, 0, 0); v /= steps; v[3] = 0; for(int i=0; i, map> maxerr; bool scatterplot; void queueline1(hyperpoint a, hyperpoint b, color_t c) { queueline(shiftless(a), shiftless(b), c); } void draw_graph() { vid.linewidth *= 2; queueline1(pt(0, 950), pt(1500, 950), 0xFF); queueline1(pt(150, 0), pt(150, 1000), 0xFF); vid.linewidth /= 2; for(int i=1; i<=9; i++) { queueline1(pt(x_to_scr(i), 950), pt(x_to_scr(i), 960), 0xFF); queuestr(x_to_scr(i), 980, 0, 60, its(i), 0, 0, 8); } for(int i=-8; i<=2; i++) { ld v = pow(10, i); queueline1(pt(140, y_to_scr(v)), pt(150, y_to_scr(v)), 0xFF); queuestr(70, y_to_scr(v), 0, 60, "1e"+its(i), 0, 0, 8); vid.linewidth /= 2; queueline1(pt(1100, y_to_scr(v)), pt(150, y_to_scr(v)), 0xFF); vid.linewidth *= 2; } vid.linewidth *= 2; for(auto& [id, graph]: maxerr) { auto& [name, col] = id; ld last = 1e-9; ld lastx = 0; for(auto [x, y]: graph) { if(scatterplot) { curvepoint(pt(x_to_scr(x)+2, y_to_scr(y))); curvepoint(pt(x_to_scr(x)-2, y_to_scr(y))); queuecurve(shiftless(Id), col, 0, PPR::LINE); curvepoint(pt(x_to_scr(x), y_to_scr(y)+2)); curvepoint(pt(x_to_scr(x), y_to_scr(y)-2)); queuecurve(shiftless(Id), col, 0, PPR::LINE); } if(y_to_scr(y) > y_to_scr(last) - x_to_scr(lastx) + x_to_scr(x)) continue; if(y > 100) y = 100; last = y; lastx = x; ld xx = x; if(xx > 9) xx = 9; if(!scatterplot) curvepoint(pt(x_to_scr(x), y_to_scr(y))); if(xx == 9) break; } if(!scatterplot) { queuestr(1100, y_to_scr(last), 0, 60, name, col >> 8, 0, 0); queuecurve(shiftless(Id), col, 0, PPR::LINE); } } vid.linewidth /= 2; drawqueue(); } void draw_sol_diffeq_graph() { } void make_graph(string fname) { start_game(); flat_model_enabler fme; shot::shotx = 1500; shot::shoty = 1000; shot::format = shot::screenshot_format::svg; svg::divby = 1; shot::take(fname, draw_graph); } void sol_diffeq_graph() { auto& s = sn::get_tabled(); s.load(); for(int x=0; x quantiles(vector data) { sort(data.begin(), data.end()); if(isize(data) <= 20) return data; vector q; for(int i=0; i<=20; i++) q.push_back(data[(isize(data)-1)*i/20]); return q; } auto smax(auto& tab, ld& i, ld x) { if(x) tab[i] = max(tab[i], x); } ld median(vector v) { sort(v.begin(), v.end()); return v[isize(v)/2]; } void sol_table_test() { // auto& length_good = maxerr[{"length/good", 0x408040FF}]; // auto& angle_good = maxerr[{"angle/good", 0x404080FF}]; // auto& length_good2 = maxerr[{"length/mid", 0x808040FF}]; // auto& angle_good2 = maxerr[{"angle/mid", 0x804080FF}]; // auto& length_bad = maxerr[{"length/bad", 0xC08040FF}]; // auto& angle_bad = maxerr[{"angle/bad", 0xC04080FF}]; // map wins; auto& s = sn::get_tabled(); s.load(); map maxerr; int good = 0, bad = 0; vector length_errors; vector angle_errors; vector split; vector lerrs[4][4][4], aerrs[4][4][4]; for(int a: {16, 32, 48, 60}) println(hlog, "xy_", a, " : ", sn::ix_to_x(a / (s.PRECX-1.))); for(int a: {16, 32, 48, 60}) println(hlog, "z_", a, " : ", sn::iz_to_z(a / (s.PRECZ-1.))); FILE *g = fopen("solv-error-data.csv", "wt"); for(ld x=0; x 0) a0++; else b0++; } bool bad_region = x > s.PRECX/2 && y > s.PRECY/2 && z < s.PRECZ/2; bool bad_break = bad_region && a0 && b0; auto ax = sn::ix_to_x(x / (s.PRECX-1.)); auto ay = sn::ix_to_x(y / (s.PRECY-1.)); auto az = sn::iz_to_z(z / (s.PRECZ-1.)); hyperpoint h = point31(ax, ay, az); hyperpoint v = inverse_exp(shiftless(h), bad_break ? pfNO_INTERPOLATION : pNORMAL); // println(hlog, "looking for ", h); // println(hlog, "exp(", v, ") = ", nisot::numerical_exp(v)); hyperpoint v1 = iterative_solve(h, v, 1e-9, false); // println(hlog, "exp(", v1, ") = ", nisot::numerical_exp(v1)); hyperpoint h2 = nisot::numerical_exp(v1); if(sqhypot_d(3, h-h2) > 1e-6) { bad++; continue; } else good++; ld dv = hypot_d(3, v); ld dv1 = hypot_d(3, v1); ld lerr = abs(dv - dv1); ld aerr = asin(hypot_d(3, v^v1) / dv / dv1); ld d = hypot_d(3, v1); if(dv == 0 || dv1 == 0) continue; if(invalid_any(v1) || invalid_any(v)) { println(hlog, "invalid"); continue; } if(isnan(aerr)) println(hlog, "v = ", v, " v1 = ", v1, "aerr"); else fprintf(g, "%lf;%lf;%lf;%lf;%lf;%lf;%lf;%lf;%d\n", x, y, z, ax, ay, az, lerr, aerr, bad_break ); lerrs[zp][yp][xp].push_back(lerr); aerrs[zp][yp][xp].push_back(aerr); } fclose(g); /* if(d >= 3 && d <= 3.1 && !bad_region) { println(hlog, tie(x,y,z), " : ", lerr); split.push_back(lerr); } if(bad_break) smax(length_bad, d, lerr), smax(angle_bad, d, aerr), 0; else if(bad_region) smax(length_good2, d, lerr), smax(angle_good2, d, aerr), 0; else smax(length_good, d, lerr), smax(angle_good, d, aerr), 0; length_errors.push_back(lerr); ld cross = hypot_d(3, v^v1) / dv / dv1; angle_errors.push_back(cross); } // println(hlog, quantiles(length_errors)); println(hlog, quantiles(split)); */ // for(auto p: angle_good) println(hlog, p); // make_graph("sol-la-errors.svg"); FILE *f = fopen("devmods/graph.tex", "wt"); fprintf(f, "\\documentclass{article}\n\\begin{document}\n"); fprintf(f, "\\small\\setlength{\\tabcolsep}{3pt}\n"); fprintf(f, "\\begin{tabular}{|c|cccc|cccc|cccc|cccc|}\n\\hline\n"); for(int z=0; z<4; z++) { fprintf(f, " & "); fprintf(f, "\\multicolumn{4}{|c%s}{$z_%d$}", z==3?"|":"", z); } fprintf(f, "|\\\\\n"); for(int z=0; z<4; z++) { for(int x=0; x<4; x++) { fprintf(f, " & "); fprintf(f, "$x_%d$", x); } } fprintf(f, "\\\\\n\\hline"); for(int y=0; y<4; y++) { fprintf(f, "$y_%d$ ", y); for(int z=0; z<4; z++) { for(int x=0; x<4; x++) { fprintf(f, " & "); fprintf(f, "%4.2g", log10(median(lerrs[z][y][x]))); } } fprintf(f, "\\\\\n"); } fprintf(f, "\\hline \n"); for(int y=0; y<4; y++) { fprintf(f, "$y_%d$ ", y); for(int z=0; z<4; z++) { for(int x=0; x<4; x++) { fprintf(f, " & "); fprintf(f, "%4.2g", log10(median(aerrs[z][y][x]))); } } fprintf(f, "\\\\\n"); } fprintf(f, "\\hline\n"); fprintf(f, "\\end{tabular}\n"); fprintf(f, "\\end{document}\n"); fclose(f); } int readArgs() { using namespace arg; if(0) ; else if(argis("-sol-diff-graph")) { sol_diffeq_graph(); } else if(argis("-sol-tabletest")) { sol_table_test(); } else if(argis("-sol-numerics")) { sol_numerics_out(); } else return 1; return 0; } auto nhook = addHook(hooks_args, 100, readArgs); } #endif