hyperrogue/expansion.cpp

// HyperRogue -- expansion_analyzer
// Copyright (C) 2011-2018 Zeno Rogue, see 'hyper.cpp' for details

namespace hr {

int subtype(cell *c) {
  return patterns::getpatterninfo(c, patterns::PAT_NONE, 0).id;
  }

bignum& bignum::operator +=(const bignum& b) {
  int K = isize(b.digits);
  if(K > isize(digits)) digits.resize(K);
  int carry = 0;
  for(int i=0; i<K || carry; i++) {
    if(i >= isize(digits)) digits.push_back(0);
    digits[i] += carry;
    if(i < K) digits[i] += b.digits[i];
    if(digits[i] >= BASE) {
      digits[i] -= BASE;
      carry = 1;
      }
    else carry = 0;
    }
  return *this;
  }

void bignum::addmul(const bignum& b, int factor) {
  int K = isize(b.digits);
  if(K > isize(digits)) digits.resize(K);
  int carry = 0;
  for(int i=0; i<K || (carry > 0 && carry < -1); i++) {
    if(i >= isize(digits)) digits.push_back(0);
    long long l = digits[i];
    l += carry;
    if(i < K) l += b.digits[i] * factor;
    carry = 0;
    if(l >= BASE) carry = l / BASE;
    if(l < 0) carry = -(BASE-1-l) / BASE;
    l -= carry * BASE;
    digits[i] = l;
    }
  if(carry < 0) digits.back() -= BASE;
  }

void operator ++(bignum &b, int) {
  int i = 0;
  while(true) {
    if(isize(b.digits) == i) { b.digits.push_back(1); break; }
    else if(b.digits[i] == bignum::BASE-1) {
      b.digits[i] = 0;
      i++;
      }
    else {
      b.digits[i]++;
      break;
      }      
    }
  }

string bignum::get_str(int max_length) {
  if(digits.empty()) return "0";
  string ret = its(digits.back());
  for(int i=isize(digits)-2; i>=0; i--) {
    if(isize(ret) > max_length && i) {
      ret += XLAT(" (%1 more digits)", its(9 * (i+1)));
      return ret;
      }

    ret += " ";
    string val = its(digits[i]);
    while(isize(val) < 9) val = "0" + val;
    ret += val;
    }
  return ret;
  }

void canonicize(vector<int>& t) {
  for(int i=2; i<isize(t); i++)
    if((t[i] & 3) == 1 && (t[i-1] & 3) != 1)
      std::rotate(t.begin()+1, t.begin()+i, t.end());
  }

vector<int> expansion_analyzer::gettype(cell *c) {
  vector<int> res;
  res.push_back(subtype(c) * 4 + 2);
  int d = celldist(c);
  for(int i=0; i<c->type; i++) {
    cell *c1 = c->cmove(i);
    res.push_back(subtype(c1) * 4 + celldist(c1) - d);
    }
  canonicize(res);
  return res;
  }

int expansion_analyzer::sample_id(cell *c) {
  auto t = gettype(c);
  if(codeid.count(t)) return codeid[t];
  auto &cit = codeid[t];
  cit = isize(samples);
  samples.push_back(c);
  return cit;
  }
  
void expansion_analyzer::preliminary_grouping() {
  samples.clear();
  codeid.clear();
  children.clear();
  sample_id(currentmap->gamestart());
  for(int i=0; i<isize(samples); i++) {
    children.emplace_back();
    auto c = samples[i];
    int d = celldist(c);
    for(int k=0; k<c->type; k++) {
      cell *c1 = c->cmove(k);
      if(celldist(c1) != d+1) continue;
      cell *c2 = c->cmove((k+1) % c->type);
      if(celldist(c2) != d+1) continue;
      children.back().push_back(sample_id(c1));
      }
    }
  N = isize(samples);
  rootid = 0;
  diskid = N;
  children.push_back(children[rootid]);
  children[diskid].push_back(diskid);
  N++;
  }

void expansion_analyzer::reduce_grouping() {
  int old_N = N;
  vector<int> grouping;
  grouping.resize(N);
  int nogroups = 1;
  for(int i=0; i<N; i++) grouping[i] = 0;
  while(true) {
    vector< pair<vector<int>, int > > childgroups(N);
    for(int i=0; i<N; i++) {
      childgroups[i].second = i;
      for(int j: children[i])
        childgroups[i].first.push_back(grouping[j]);
      sort(childgroups[i].first.begin(), childgroups[i].first.end());
      }
    sort(childgroups.begin(), childgroups.end());
    int newgroups = 0;
    for(int i=0; i<N; i++) {
      if(i == 0 || childgroups[i].first != childgroups[i-1].first) newgroups++;
      grouping[childgroups[i].second] = newgroups - 1;
      }
    if(nogroups == newgroups) break;
    nogroups = newgroups;
    }
  
  vector<int> groupsample(nogroups);
  for(int i=0; i<N; i++) groupsample[grouping[i]] = i;
  vector<vector<int>> newchildren(nogroups);
  for(int i=0; i<nogroups; i++) 
    for(int j: children[groupsample[i]])
      newchildren[i].push_back(grouping[j]);
  children = move(newchildren);
  for(auto& p: codeid) p.second = grouping[p.second];
  N = nogroups;
  rootid = grouping[rootid];
  diskid = grouping[diskid];
  printf("%d -> %d\n", old_N, N);
  for(int g=0; g<old_N; g++) if(grouping[g] != g) descendants.clear();
  }

template<class T> int size_upto(vector<T>& v, int s) {
  int res = isize(v);
  if(res < s) v.resize(s);
  return res;
  }

bignum& expansion_analyzer::get_descendants(int level) {
  if(!N) preliminary_grouping(), reduce_grouping();
  return get_descendants(level, rootid);
  }

bignum& expansion_analyzer::get_descendants(int level, int type) {
  auto& pd = descendants;
  size_upto(pd, level+1);
  for(int d=0; d<=level; d++)
  for(int i=size_upto(pd[d], N); i<N; i++)
    if(d == 0) pd[d][i].be(1);
    else for(int j: children[i]) 
      pd[d][i] += pd[d-1][j];
  return pd[level][type];
  }

bool expansion_analyzer::verify(int id) {
  if(id < isize(coef)) return false;
  long long res = 0;
  for(int t=0; t<isize(coef); t++)
    res += coef[t] * get_descendants(id-t-1).approx_ll();
  return res == get_descendants(id).approx_ll();
  }

int expansion_analyzer::valid(int v, int step) {
  if(step < 0) return 0;
  if(get_descendants(step+v+v+5).approx_int() >= bignum::BASE) return 0;
  ld matrix[100][128];
  for(int i=0; i<v; i++)
  for(int j=0; j<v+1; j++)
    matrix[i][j] = get_descendants(step+i+j).approx();
  
  for(int k=0; k<v; k++) {
    int nextrow = k;
    while(nextrow < v && std::abs(matrix[nextrow][k]) < 1e-6)
      nextrow++;
    if(nextrow == v) return 1;
    if(nextrow != k) {
      // printf("swap %d %d\n", k, nextrow);
      for(int l=0; l<=v; l++) swap(matrix[k][l], matrix[nextrow][l]);
      // display();
      }
    ld divv = 1. / matrix[k][k];
    for(int k1=k; k1<=v; k1++) matrix[k][k1] *= divv;
    // printf("divide %d\n", k);
    // display();
    for(int k1=k+1; k1<v; k1++) if(matrix[k1][k] != 0) {
      ld coef = -matrix[k1][k];
      for(int k2=k; k2<=v; k2++) matrix[k1][k2] += matrix[k][k2] * coef;
      }
    // printf("zeros below %d\n", k);
    // display();
    }
  
  for(int k=v-1; k>=0; k--)
  for(int l=k-1; l>=0; l--)
    if(matrix[l][k]) matrix[l][v] -= matrix[l][k] * matrix[k][v];
  
  coef.resize(v);
  for(int i=0; i<v; i++) coef[i] = int(floor(matrix[v-1-i][v] + .5));
    
  for(int t=step+v; t<step+v+v+5; t++) if(!verify(t)) return 2;
  tested_to = step+v+v+5;
  while(tested_to < step+v+v+100 && get_descendants(tested_to).approx_ll() < bignum::BASE2) {
    if(!verify(tested_to)) return 2;
    tested_to++;
    }
  
  valid_from = step+v;
  return 3;
  }

void expansion_analyzer::find_coefficients() {
  if(coefficients_known) return;
  if(!N) preliminary_grouping(), reduce_grouping();
  for(int v=1; v<25; v++) 
  for(int step=0; step<1000; step++) { 
    int val = valid(v, step);
    if(val == 0) break;
    if(val == 3) { coefficients_known = 2; return; }
    }
  coefficients_known = 1;
  }

ld growth;

ld expansion_analyzer::get_growth() {
  if(growth >= 1) return growth;
  if(!N) preliminary_grouping(), reduce_grouping();
  vector<ld> eigen(N, 1);
  ld total;
  for(int iter=0; iter<100000; iter++) {
    total = 0;
    vector<ld> neweigen(N, 0);
    for(int i=0; i<N; i++) {
      for(int j: children[i]) neweigen[i] += eigen[j];
      total += neweigen[i];
      }
    for(int i=0; i<N; i++) eigen[i] = .1 * eigen[i] + .9 * neweigen[i] / total;
    // for(int i=0; i<N; i++) printf("%lf ", eigen[i]); printf("total = %lf\n", total);
    }
  return growth = total;
  }

void expansion_analyzer::reset() {
  N = 0;
  growth = 0;
  coefficients_known = 0;
  samples.clear();
  codeid.clear();
  children.clear();
  coef.clear();
  descendants.clear();
  }

template<class T> int type_in(expansion_analyzer& ea, cell *c, const T& f) {
  if(!ea.N) ea.preliminary_grouping(), ea.reduce_grouping();
  vector<int> res;
  res.push_back(subtype(c) * 4 + 2);
  int d = f(c);
  for(int i=0; i<c->type; i++) {
    cell *c1 = c->cmove(i);
    res.push_back(subtype(c1) * 4 + f(c1) - d);
    }
  
  canonicize(res);
  if(ea.codeid.count(res)) return ea.codeid[res];
  int ret = ea.N++;
  ea.codeid[res] = ret;
  
  vector<int> rec(MAX_EDGE, -1);
  
  ea.children.emplace_back();
  for(int k=0; k<c->type; k++) {
    cell *c1 = c->cmove(k);
    if(f(c1) != d+1) continue;
    cell *c2 = c->cmove((k+1) % c->type);
    if(f(c2) != d+1) continue;
    auto ti = rec[k] = type_in(ea, c1, f);
    // note: ea.children[ret].push_back(type_in(...)) would not work because of invalidation
    ea.children[ret].push_back(ti);
    }

  /*
  printf("extra type created: [%d, d=%d]", ret, d); for(auto i: ea.children[ret]) printf(" %d", i); printf("\n");
  printf("  list:");
  for(int k=0; k<c->type; k++) {
    cell *c1 = c->cmove(k);
    printf(" %d/%d/%d", k, rec[k], f(c1) - d);
    }
  printf("\n");
  */

  return ret;
  }

template<class T> int type_in_quick(expansion_analyzer& ea, cell *c, const T& f) {
  vector<int> res;
  res.push_back(subtype(c) * 4 + 2);
  int d = f(c);
  for(int i=0; i<c->type; i++) {
    cell *c1 = c->cmove(i);
    int dd = f(c1) - d;
    if(dd < -1 || dd > 1) return -1;
    res.push_back(subtype(c1) * 4 + dd);
    }
  
  canonicize(res);
  if(ea.codeid.count(res)) return ea.codeid[res];
  return -1;
  }

bool sizes_known() {
  if(bounded) return false;
  // Castle Anthrax is infinite
  if(binarytiling) return false;
  // not implemented
  if(archimedean) return false;
  return true;  
  }

string expansion_analyzer::approximate_descendants(int d, int max_length) {
  auto t = SDL_GetTicks();
  while(isize(descendants) <= d && SDL_GetTicks() < t + 100)
    get_descendants(isize(descendants));
  if(isize(descendants) > d) 
    return get_descendants(d).get_str(max_length);
  int v = isize(descendants) - 1;
  bignum& b = get_descendants(v);
  if(b.digits.empty()) return "0";
  ld log_10 = log(b.digits.back()) / log(10) + 9 * (isize(b.digits) - 1) + (d - v) * log(get_growth()) / log(10);
  int more_digits = int(log_10);
  return XLAT("about ") + fts(pow(10, log_10 - more_digits)) + "E" + its(more_digits);
  }

enum eDistanceFrom { dfPlayer, dfStart, dfWorld };
string dfnames[3] = { "player", "start", "land" };

eDistanceFrom distance_from = dfPlayer;

enum eNumberCoding { ncNone, ncDistance, ncType, ncDebug };
string ncnames[4] = { "NO", "distance", "type", "debug" };
eNumberCoding number_coding = ncDistance;

bool mod_allowed() {
  return cheater || autocheat || archimedean || tour::on;
  }

int curr_dist(cell *c) {
  switch(distance_from) {
    case dfPlayer:
      return c->cpdist < INFD ? c->cpdist : celldistance(cwt.at, c);
    case dfStart:
      return celldist(c);
    case dfWorld:
      if(!mod_allowed() && !among(c->land, laOcean, laIvoryTower, laEndorian, laDungeon, laTemple, laWhirlpool))
        return 0;
      if(c->master->alt) return celldistAlt(c);
      return inmirror(c) ? (c->landparam & 255) : c->landparam;
    }
  return 0;    
  }

int position;

template<class T> int type_in_reduced(expansion_analyzer& ea, cell *c, const T& f) {
  int a = ea.N;
  int t = type_in(ea, c, f);
  if(expansion.N != a) {
    expansion.reduce_grouping();
    t = type_in(ea, c, f);
    }
  return t;
  }

bool viewdists = false, use_color_codes = true, use_analyzer = true;

int first_distance = 0, scrolltime = 0;
bool scrolling_distances = false;

color_t distribute_color(int id) {
  color_t v = 0xFFFFFF;
  for(int z=0; z<24; z++) if(id & (1<<z)) part(v, (id%3)) &=~ 1<<(7-(z/3));
  return v; 
  }

void do_viewdist(cell *c, const transmatrix& V, color_t& wcol, color_t& fcol) {
  if(behindsphere(V)) return;

  int cd = (use_color_codes || number_coding == ncDistance || number_coding == ncDebug) ? curr_dist(c) : 0;
  
  if(use_color_codes) {
    int dc = distcolors[cd&7];
    wcol = gradient(wcol, dc, 0, .4, 1);
    fcol = gradient(fcol, dc, 0, .4, 1);
    }
  
  string label = "";
  int dc;
 
  switch(number_coding) {
    case ncDistance: { 
      label = its(cd);
      dc = distcolors[cd&7];
      break;
      }
    case ncType: {
      int t = type_in_quick(expansion, c, curr_dist);
      if(t >= 0) label = its(t), dc = distribute_color(t);
      break;
      }
    case ncDebug: {
      int d = celldistance(c, distance_from == dfPlayer ? cwt.at : currentmap->gamestart());
      dc = (d != cd) ? 0xFF0000 : 0x00FF00;
      label = its(d);
      }
    case ncNone: ;
    }

  // string label = its(fieldpattern::getriverdistleft(c)) + its(fieldpattern::getriverdistright(c));
  /* queuepolyat(V, shFloor[ct6], darkena(gradient(0, distcolors[cd&7], 0, .25, 1), fd, 0xC0),
    PPR::TEXT); */
  if(label != "")
    queuestr(V, (isize(label) > 1 ? .6 : 1), label, 0xFF000000 + dc, 1);
  }

void viewdist_configure_dialog() {
  dialog::init("");
  cmode |= sm::SIDE | sm::MAYDARK | sm::EXPANSION;
  gamescreen(0);
  
  dialog::addSelItem("which distance", dfnames[distance_from], 'c');
  dialog::add_action([] () { distance_from = mod_allowed() ? eDistanceFrom(2 - distance_from) : eDistanceFrom((distance_from + 1) % 3); });
             
  dialog::addSelItem("number codes", ncnames[number_coding], 'n');
  dialog::add_action([] () { number_coding = eNumberCoding((number_coding + 1) % (mod_allowed() ? 4 : 2)); });

  dialog::addBoolItem("color codes", use_color_codes, 'u');
  dialog::add_action([] () { use_color_codes = !use_color_codes; });

  dialog::addBoolItem("use analyzer", use_analyzer, 'a');
  dialog::add_action([] () { use_analyzer = !use_analyzer; });

  dialog::addSelItem("display distances from", its(first_distance), 'd');
  dialog::add_action([] () { 
    scrolling_distances = false;
    dialog::editNumber(first_distance, 0, 3000, 0, 1, "display distances from", "");
    });
  
  for(auto& lp: linepatterns::patterns) {
    using namespace linepatterns;
    int id = 0;
    if(among(lp.id, patTriTree, patTriRings, patTriOther)) {
      dialog::addColorItem(XLAT(lp.lpname), lp.color, '1'+(id++));
      dialog::add_action([&lp] () {
        dialog::openColorDialog(lp.color, NULL);
        dialog::dialogflags |= sm::MAYDARK | sm::SIDE | sm::EXPANSION;
        });
      }
    }
  
  if(!mod_allowed()) 
    dialog::addInfo(XLAT("note: enable the cheat mode for additional options"));
  else 
    dialog::addBreak(100);
  
  if(distance_from)
    dialog::addInfo("numbers show the descendants of current player position");
  else
    dialog::addBreak(100);

  dialog::addBack();  
  dialog::display();
  }

bool is_descendant(cell *c) {
  if(c == cwt.at) return true;
  if(curr_dist(c) < curr_dist(cwt.at)) return false;
  return is_descendant(chosenDown(c, -1, 0, curr_dist));
  }

const int scrollspeed = 100;
    
void expansion_analyzer::view_distances_dialog() {
  static int lastticks;
  if(scrolling_distances && !bounded) {
    scrolltime += SDL_GetTicks() - lastticks;
    first_distance += scrolltime / scrollspeed;
    scrolltime %= scrollspeed;
    }
  lastticks = SDL_GetTicks();
  
  distcolors[0] = forecolor;
  dialog::init("");
  cmode |= sm::DIALOG_STRICT_X | sm::EXPANSION;
  
  int maxlen = bounded ? 128 : 16 + first_distance;
  vector<bignum> qty(maxlen);
  
  bool really_use_analyzer = use_analyzer && sizes_known();
  
  if(really_use_analyzer) {
    int t = type_in_reduced(expansion, cwt.at, curr_dist);
    for(int r=0; r<maxlen; r++)
      qty[r] = expansion.get_descendants(r, t);
    }
  else {
    if(distance_from == dfPlayer) {
      celllister cl(cwt.at, bounded ? maxlen-1 : gamerange(), 100000, NULL);
      for(int d: cl.dists)
        if(d >= 0 && d < maxlen) qty[d]++;
      }
    else {
      celllister cl(cwt.at, bounded ? maxlen-1 : gamerange(), 100000, NULL);
      for(cell *c: cl.lst) if(is_descendant(c)) qty[curr_dist(c)]++;
      }
    if(sizes_known()) {
      find_coefficients();
      if(gamerange()+1 >= valid_from && coefficients_known == 2) {
        for(int i=gamerange()+1; i<maxlen; i++)
          for(int j=0; j<isize(coef); j++) {
            qty[i].addmul(qty[i-1-j], coef[j]);
            }
        }
      }
    }
  
  dialog::addBreak(100 - 100 * scrolltime / scrollspeed);

  for(int i=first_distance; i<maxlen; i++) if(!qty[i].digits.empty())
    dialog::addInfo(its(i) + ": " + qty[i].get_str(100), distcolors[i&7]);
  
  dialog::addBreak(100 * scrolltime / scrollspeed);

  if(binarytiling){
    dialog::addBreak(200);
    dialog::addInfo("a(d) ~ 2ᵈ");
    }
  else if(sizes_known()) {
    if(euclid) {
      dialog::addBreak(200);
      dialog::addInfo("a(d) = " + its(get_descendants(10).approx_int() - get_descendants(9).approx_int()) + "d", forecolor);
      }
    else {
      dialog::addBreak(100);
      
      if(coefficients_known == 2) {
        string fmt = "a(d+" + its(isize(coef)) + ") = ";
        bool first = true;
        for(int i=0; i<isize(coef); i++) if(coef[i]) {
          if(first && coef[i] == 1) ;
          else if(first) fmt += its(coef[i]);
          else if(coef[i] == 1) fmt += " + ";
          else if(coef[i] == -1) fmt += " - ";
          else if(coef[i] > 1) fmt += " + " + its(coef[i]);
          else if(coef[i] < -1) fmt += " - " + its(-coef[i]);
          fmt += "a(d";
          if(i != isize(coef) - 1)
            fmt += "+" + its(isize(coef) - 1 - i);
          fmt += ")";
          first = false;
          }
        fmt += " (d>" + its(valid_from-1) + ")";
        dialog::addHelp(fmt);
        }
      else dialog::addBreak(100);
      
      char buf[20];
      snprintf(buf, 20, "%.8lf", (double) get_growth());
      dialog::addInfo("a(d) = Θ(" + string(buf) + "...ᵈ)", forecolor);
      }
    }
  
  dialog::addItem("scroll", 'S');
  dialog::addItem("configure", 'C');
  dialog::display();
  }

void enable_viewdists() {
  first_distance = 0;
  scrolltime = 0;
  viewdists = true;
  if(!mod_allowed()) {
    number_coding = ncDistance;
    distance_from = dfPlayer;
    }
  }

bool expansion_handleKey(int sym, int uni) {
  if((cmode & sm::EXPANSION) && (cmode & sm::NORMAL)) {
    if(uni == 'S') scrolling_distances = !scrolling_distances;
    else if(uni == 'C') pushScreen(viewdist_configure_dialog);
    else if(sym == SDLK_ESCAPE) first_distance = 0, viewdists = false;
    else return false;
    return true;
    }
  return false;
  }

int expansion_hook = addHook(hooks_handleKey, 0, expansion_handleKey);

#if !CAP_MINI
void compute_coefficients() {
  printf("%s %s\n", gp::operation_name().c_str(), ginf[geometry].name);
  start_game();
  
    printf("  sizes:");
    for(int i=0; i<10; i++) printf(" %d", expansion.get_descendants(i).approx_int());
    
    printf("  N = %d\n", expansion.N);

  expansion.find_coefficients();      
  if(expansion.coefficients_known == 2) {
    printf("  coefficients:"); for(int x: expansion.coef) printf(" %d", x);
    printf(" (tested on %d to %d)\n", expansion.valid_from, expansion.tested_to);
    }
  }

int readArgs() {
  using namespace arg;
           
  if(0) ;
  else if(argis("-vap")) { 
    PHASEFROM(2); 
    start_game();
    while(true) {
      string s = expansion.approximate_descendants(10000, 100);
      printf("s = %s\n", s.c_str());
      if(isize(expansion.descendants) >= 10000) break;
      }
    }
  else if(argis("-csizes")) { 
    PHASEFROM(2); 
    start_game();
    expansion.get_growth();
    for(int i=0; i<30; i++)
      printf("%s / %s\n", expansion.get_descendants(i).get_str(1000).c_str(), expansion.get_descendants(i, expansion.diskid).get_str(1000).c_str());  
    }
  else if(argis("-csolve")) { 
    PHASEFROM(2); 
    start_game();
    printf("preliminary_grouping...\n");
    expansion.preliminary_grouping();
    printf("N = %d\n", expansion.N);
    printf("reduce_grouping...\n");
    expansion.reduce_grouping();
    printf("N = %d\n", expansion.N);
    printf("growth = %lf\n", (double) expansion.get_growth());
    expansion.find_coefficients();      
    if(expansion.coefficients_known == 2) {
      printf("coefficients:"); for(int x: expansion.coef) printf(" %d", x);
      printf(", valid from %d to %d\n", expansion.valid_from, expansion.tested_to);
      }
    }
  else if(argis("-csolve_tab")) {
    for(eGeometry geo: {gNormal, gOctagon, g45, g46, g47}) {
      set_geometry(geo);
      set_variation(eVariation::pure);
      compute_coefficients();
      set_variation(eVariation::bitruncated);
      compute_coefficients();
      for(int x=1; x<9; x++)
      for(int y=0; y<=x; y++) {
        if(x == 1 && y == 0) continue;
        if(x == 1 && y == 1 && S3 == 3) continue;
        if(x+y > 10) continue;
        stop_game();
        gp::param = gp::loc(x, y);
        need_reset_geometry = true;
        set_variation(eVariation::goldberg);
        compute_coefficients();
        }
      }
    }
  else return 1;
  return 0;
  }

auto ea_hook = addHook(hooks_args, 100, readArgs);
#endif

expansion_analyzer expansion;

}