diff --git a/kohonen.cpp b/kohonen.cpp
index 896d09bb..2a4e9d81 100644
--- a/kohonen.cpp
+++ b/kohonen.cpp
@@ -1,1112 +1,1112 @@
-// Hyperbolic Rogue
-// Copyright (C) 2011-2017 Zeno and Tehora Rogue, see 'hyper.cpp' for details
-
-// Kohonen's self-organizing networks. 
-// This is a part of RogueViz, not a part of HyperRogue.
-
-namespace kohonen {
-
-int cols;
-
-typedef vector<double> kohvec;
-
-struct sample {
-  kohvec val;
-  string name;
-  };
-
-vector<sample> data;
-
-vector<int> samples_shown;
-
-int whattodraw[3] = {-2,-2,-2};
-
-struct neuron {
-  kohvec net;
-  cell *where;
-  double udist;
-  int lpbak;
-  int col;
-  int samples, csample, bestsample;
-  };
-
-vector<string> colnames;
-
-kohvec weights;
-
-vector<neuron> net;
-
-int neuronId(neuron& n) { return &n - &(net[0]); }
-
-void alloc(kohvec& k) { k.resize(cols); }
-
-bool neurons_indexed = false;
-
-int samples;
-
-template<class T> T sqr(T x) { return x*x; }
-
-vector<neuron*> whowon;
-
-void normalize() {
-  alloc(weights);
-  for(int k=0; k<cols; k++) {
-    double sum = 0, sqsum = 0;
-    for(sample& s: data)
-      sum += s.val[k],
-      sqsum += s.val[k] * s.val[k];
-    double variance = sqsum/samples - sqr(sum/samples);
-    weights[k] = 1 / sqrt(variance);
-    }
-  }
-
-double vnorm(kohvec& a, kohvec& b) {
-  double diff = 0;
-  for(int k=0; k<cols; k++) diff += sqr((a[k]-b[k]) * weights[k]);
-  return diff;
-  }
-
-void sominit(int);
-void uninit(int);
-
-void loadsamples(const char *fname) {
-  FILE *f = fopen(fname, "rt");
-  if(!f) {
-    fprintf(stderr, "Could not load samples\n");
-    return;
-    }
-  if(fscanf(f, "%d", &cols) != 1) { fclose(f); return; }
-  while(true) {
-    sample s;
-    bool shown = false;
-    alloc(s.val);
-    if(feof(f)) break;
-    for(int i=0; i<cols; i++)
-      if(fscanf(f, "%lf", &s.val[i]) != 1) { goto bigbreak; }
-    fgetc(f);
-    while(true) {
-      int c = fgetc(f);
-      if(c == -1 || c == 10 || c == 13) break;
-      if(c == '!' && s.name == "") shown = true;
-      else if(c != 32 && c != 9) s.name += c;
-      }
-    if(shown) samples_shown.push_back(size(data));
-    data.push_back(move(s));
-    }
-  bigbreak:
-  fclose(f);
-  samples = size(data);
-  normalize();
-  colnames.resize(cols);
-  for(int i=0; i<cols; i++) colnames[i] = "Column " + its(i);
-  uninit(0); sominit(1);
-  }
-                  
-int tmax = 30000;
-double distmul = 1;
-double learning_factor = .1;
-
-int qpct = 100;
-
-int t, lpct, cells;
-double maxdist;
-
-neuron& winner(int id) {
-  double bdiff = 1e20;
-  neuron *bcell = NULL;
-  for(neuron& n: net) {
-    double diff = vnorm(n.net, data[id].val);
-    if(diff < bdiff) bdiff = diff, bcell = &n;
-    }
-  return *bcell;
-  }
-
-void setindex(bool b) {
-  if(b == neurons_indexed) return;
-  neurons_indexed = b;
-  if(b) {
-    for(neuron& n: net) n.lpbak = n.where->landparam, n.where->landparam = neuronId(n);
-    }
-  else {
-    for(neuron& n: net) n.where->landparam = n.lpbak;
-    }
-  }
-
-neuron *getNeuron(cell *c) {
-  if(!c) return NULL;
-  setindex(true);
-  if(c->landparam < 0 || c->landparam >= cells) return NULL;
-  neuron& ret = net[c->landparam];
-  if(ret.where != c) return NULL;
-  return &ret;
-  }
-
-neuron *getNeuronSlow(cell *c) {
-  if(neurons_indexed) return getNeuron(c);
-  for(neuron& n: net) if(n.where == c) return &n;
-  return NULL;
-  }
-
-double maxudist;
-
-neuron *distfrom;
-
-bool noshow = false;
-
-void coloring() {
-  if(noshow) return;
-  setindex(false);
-  
-  bool besttofind = true;
-
-  for(int pid=0; pid<3; pid++) {
-    int c = whattodraw[pid];
-    
-    if(c == -5) {
-      if(besttofind) {
-        besttofind = false;
-        for(neuron& n: net) {
-          double bdiff = 1e20;
-          for(int i=0; i<size(samples_shown); i++) {
-            double diff = vnorm(n.net, data[samples_shown[i]].val);
-            if(diff < bdiff) bdiff = diff, n.bestsample = i;
-            }
-          }
-        }
-
-      for(int i=0; i<cells; i++) 
-        part(net[i].where->landparam, pid) = part(vdata[net[i].bestsample].cp.color1, pid+1);
-      }
-
-    else {
-      vector<double> listing;
-      for(neuron& n: net) switch(c) {
-        case -4:
-          listing.push_back(log(5+n.samples));
-          break;
-          
-        case -3:
-          if(distfrom) 
-            listing.push_back(vnorm(n.net, distfrom->net));
-          else
-            listing.push_back(0);
-          break;
-        
-        case -2:
-          listing.push_back(n.udist);
-          break;
-        
-        case -1: 
-          listing.push_back(-n.udist);
-          break;
-        
-        default:
-          listing.push_back(n.net[c]);
-          break;
-        }
-      
-      double minl = listing[0], maxl = listing[0];
-      for(double& d: listing) minl = min(minl, d), maxl = max(maxl, d);
-      if(maxl-minl < 1e-3) maxl = minl+1e-3;
-      
-      for(int i=0; i<cells; i++) 
-        part(net[i].where->landparam, pid) = (255 * (listing[i] - minl)) / (maxl - minl);
-      }
-    }
-  }
-
-void analyze() {
-
-  setindex(true);
-  
-  maxudist = 0;
-  for(neuron& n: net) {
-    int qty = 0;
-    double total = 0;
-    forCellEx(c2, n.where) {
-      neuron *n2 = getNeuron(c2);
-      if(!n2) continue;
-      qty++;
-      total += sqrt(vnorm(n.net, n2->net));
-      }
-    n.udist = total / qty;
-    maxudist = max(maxudist, n.udist);
-    }
-    
-  if(!noshow) {
-
-    whowon.resize(samples);
-    
-    for(neuron& n: net) n.samples = 0;
-    
-    for(int id=0; id<size(samples_shown); id++) {
-      int s = samples_shown[id];
-      auto& w = winner(s);
-      whowon[s] = &w;
-      w.samples++;
-      }
-    
-    for(int id=0; id<size(samples_shown); id++) {
-      int s = samples_shown[id];
-      auto& w = *whowon[s];
-      vdata[id].m->base = w.where;
-      vdata[id].m->at = 
-        spin(2*M_PI*w.csample / w.samples) * xpush(.25 * (w.samples-1) / w.samples);
-      w.csample++;
-      }
-    
-    shmup::fixStorage();  
-    setindex(false);
-    }
-  
-  coloring();
-  }
-
-// traditionally Gaussian blur is used in the Kohonen algoritm
-// but it does not seem to make much sense in hyperbolic geometry
-// especially wrapped one.
-// GAUSSIAN==1: use the Gaussian blur, on celldistance
-// GAUSSIAN==2: use the Gaussian blur, on true distance
-// GAUSSIAN==0: simulate the dispersion on our network
-
-int gaussian = 0;
-
-double mydistance(cell *c1, cell *c2) {
-  if(gaussian == 2) return hdist(tC0(shmup::ggmatrix(c1)), tC0(shmup::ggmatrix(c2)));
-  else return celldistance(c1, c2);
-  }
-
-struct cellcrawler {
-
-  struct cellcrawlerdata {
-    cellwalker orig;
-    int from, spin, dist;
-    cellwalker target;
-    cellcrawlerdata(const cellwalker& o, int fr, int sp) : orig(o), from(fr), spin(sp) {}
-    };
-  
-  vector<cellcrawlerdata> data;
-  
-  void store(const cellwalker& o, int from, int spin) {
-    if(eq(o.c->aitmp, sval)) return;
-    o.c->aitmp = sval;
-    data.emplace_back(o, from, spin);
-    }
-  
-  void build(const cellwalker& start) {
-    sval++;
-    data.clear();
-    store(start, 0, 0);
-    for(int i=0; i<size(data); i++) {
-      cellwalker cw0 = data[i].orig;
-      for(int j=0; j<cw0.c->type; j++) {
-        cellwalker cw = cw0;
-        cwspin(cw, j); cwstep(cw);
-        if(!getNeuron(cw.c)) continue;
-        store(cw, i, j);
-        }
-      }
-    if(gaussian) for(cellcrawlerdata& s: data)
-      s.dist = mydistance(s.orig.c, start.c);
-    }
-  
-  void sprawl(const cellwalker& start) {
-    data[0].target = start;
-    
-    for(int i=1; i<size(data); i++) {
-      cellcrawlerdata& s = data[i];
-      s.target = data[s.from].target;
-      if(!s.target.c) continue;
-      cwspin(s.target, s.spin);
-      if(cwstepcreates(s.target)) s.target.c = NULL;
-      else cwstep(s.target);
-      }
-    }
-  };
-
-cellcrawler scc[2]; // hex and non-hex
-
-double dispersion_end_at = 1.5;
-
-double dispersion_precision = .0001;
-int dispersion_each = 1;
-
-vector<vector<ld>> dispersion[2]; 
-
-int dispersion_count;
-
-void buildcellcrawler(cell *c) {
-  int sccid = c->type != 6;
-  
-  cellcrawler& cr =  scc[sccid];
-  cr.build(cellwalker(c,0));
-
-  if(!gaussian) {
-    vector<ld> curtemp;
-    vector<ld> newtemp;
-    vector<int> qty;
-    vector<pair<ld*, ld*> > pairs;
-    int N = size(net);
-    
-    curtemp.resize(N, 0);
-    newtemp.resize(N, 0);
-    qty.resize(N, 0);
-  
-    for(int i=0; i<N; i++) 
-    forCellEx(c2, net[i].where) {
-      neuron *nj = getNeuron(c2);
-      if(nj) {
-        pairs.emplace_back(&curtemp[i], &newtemp[neuronId(*nj)]);
-        qty[i]++;
-        }
-      }
-    
-    curtemp[neuronId(*getNeuron(c))] = 1;
-  
-    ld vmin = 0, vmax = 1;
-    int iter;
-    
-    auto &d = dispersion[sccid];
-    
-    d.clear();
-  
-    printf("Building dispersion...\n");
-    
-    for(iter=0; dispersion_count ? true : vmax > vmin * dispersion_end_at; iter++) {
-      if(iter % dispersion_each == 0) {
-        d.emplace_back(N);
-        auto& dispvec = d.back();
-        for(int i=0; i<N; i++) dispvec[i] = curtemp[neuronId(*getNeuron(cr.data[i].orig.c))] / vmax;
-        if(size(d) == dispersion_count) break;
-        }
-      double df = dispersion_precision * (iter+1);
-      double df0 = df / ceil(df);
-      for(int i=0; i<df; i++) {
-        for(auto& p: pairs) 
-          *p.second += *p.first;
-        for(int i=0; i<N; i++) {
-          curtemp[i] += (newtemp[i] / qty[i] - curtemp[i]) * df0;
-          newtemp[i] = 0;
-          }
-        }
-      vmin = vmax = curtemp[0];
-      for(int i=0; i<N; i++) 
-        if(curtemp[i] < vmin) vmin = curtemp[i];
-        else if(curtemp[i] > vmax) vmax = curtemp[i];
-      }
-  
-    dispersion_count = size(d);
-    printf("Dispersion count = %d\n", dispersion_count);
-    }
-  }
-
-bool finished() { return t == 0; }
-
-int krad;
-
-double ttpower = 1;
-
-void sominit(int);
-
-void step() {
-
-  if(t == 0) return;
-  sominit(2);
-  
-  double tt = (t-1.) / tmax;
-  tt = pow(tt, ttpower);
-
-  double sigma = maxdist * tt;
-  int dispid = int(dispersion_count * tt);
-
-  if(qpct) {
-    int pct = (int) ((qpct * (t+.0)) / tmax);
-    if(pct != lpct) {
-      printf("pct %d lpct %d\n", pct, lpct);
-      lpct = pct;
-      analyze();
-      
-      if(gaussian)
-        printf("t = %6d/%6d %3d%% sigma=%10.7lf maxudist=%10.7lf\n", t, tmax, pct, sigma, maxudist);
-      else
-        printf("t = %6d/%6d %3d%% dispid=%5d maxudist=%10.7lf\n", t, tmax, pct, dispid, maxudist);
-      }
-    }
-  int id = hrand(samples);
-  neuron& n = winner(id);
-  whowon.resize(samples);
-  whowon[id] = &n;
-    
-  /* 
-  for(neuron& n2: net) {
-    int d = celldistance(n.where, n2.where);
-    double nu = learning_factor; 
-//  nu *= exp(-t*(double)maxdist/perdist);
-//  nu *= exp(-t/t2);
-    nu *= exp(-sqr(d/sigma));
-    for(int k=0; k<cols; k++)
-      n2.net[k] += nu * (irisdata[id][k] - n2.net[k]);
-    } */
-    
-  int sccid = n.where->type != 6;
-  cellcrawler& s = scc[sccid];
-  s.sprawl(cellwalker(n.where, 0));
-
-  vector<double> fake(1,1);
-  auto it = gaussian ? fake.begin() : dispersion[sccid][dispid].begin();
-
-  for(auto& sd: s.data) {
-    neuron *n2 = getNeuron(sd.target.c);
-    if(!n2) continue;
-    double nu = learning_factor;
-    
-    if(gaussian)
-      nu *= exp(-sqr(sd.dist/sigma));
-    else
-      nu *= *(it++);
-
-    for(int k=0; k<cols; k++)
-      n2->net[k] += nu * (data[id].val[k] - n2->net[k]);
-    }
-  
-  t--;
-  if(t == 0) analyze();
-  }
-
-int initdiv = 1;
-
-int inited = 0;
-
-void uninit(int initto) {
-  if(inited > initto) inited = initto;
-  }
-
-void showsample(int id) {
-  for(int ii: samples_shown) 
-    if(ii == id) 
-      return;
-  int i = vdata.size();
-  samples_shown.push_back(id);
-  vdata.emplace_back();
-  auto& v = vdata.back();
-  v.name = data[id].name;
-  v.cp = dftcolor;
-  createViz(i, cwt.c, Id);
-  v.m->store();
-  }
-
-void showsample(string s) {
-  if(s == "") return;
-  for(int i=0; i<samples; i++) {
-    if(s[0] != '*' && data[i].name == s)
-      showsample(i);
-    if(s[0] == '*' && data[i].name.find(s.substr(1)) != string::npos)
-      showsample(i);
-    }
-  }
-
-void showbestsamples() {
-  vector<int> samplesbak;
-  for(auto& n: net) 
-    if(n.samples)
-      showsample(n.bestsample);
-  analyze();
-  for(auto& n: net) n.samples = 0;
-  for(int i=0; i<samples; i++) 
-    if(whowon[i])
-      whowon[i]->samples++;
-  }
-  
-void sominit(int initto) {
-
-  if(inited < 1 && initto >= 1) {
-    inited = 1;
-    if(!samples) {
-      fprintf(stderr, "Error: SOM without samples\n");
-      exit(1);
-      }
-    
-    init(); kind = kKohonen;
-    
-    /* if(geometry != gQuotient1) {
-      targetGeometry = gQuotient1;
-      restartGame('g');
-      }
-    if(!purehepta) restartGame('7'); */
-    
-    printf("Initializing SOM (1)\n");
-  
-    vector<cell*> allcells;
-    
-    if(krad) {
-      celllister cl(cwt.c, krad, 1000000, NULL);
-      allcells = cl.lst;
-      }
-    else allcells = currentmap->allcells();
-  
-    cells = size(allcells);
-    net.resize(cells);
-    for(int i=0; i<cells; i++) net[i].where = allcells[i], allcells[i]->landparam = i;
-    for(int i=0; i<cells; i++) {
-      net[i].where->land = laCanvas;
-      alloc(net[i].net);
-  
-      for(int k=0; k<cols; k++)
-      for(int z=0; z<initdiv; z++)
-        net[i].net[k] += data[hrand(samples)].val[k] / initdiv;
-      }
-      
-    for(neuron& n: net) for(int d=BARLEV; d>=7; d--) setdist(n.where, d, NULL);
-  
-    printf("samples = %d (%d) cells = %d\n", samples, size(samples_shown), cells);
-
-    if(!noshow) {
-      vdata.resize(size(samples_shown));
-      for(int i=0; i<size(samples_shown); i++) {
-        vdata[i].name = data[samples_shown[i]].name;
-        vdata[i].cp = dftcolor;
-        createViz(i, cwt.c, Id);
-        }
-      
-      storeall();
-      }
-
-    analyze();
-    }
-  
-  if(inited < 2 && initto >= 2) {
-    inited = 2;
-
-    printf("Initializing SOM (2)\n");
-
-    if(gaussian) {
-      printf("dist = %lf\n", mydistance(net[0].where, net[1].where));
-      cell *c1 = net[cells/2].where;
-      vector<double> mapdist;
-      for(neuron &n2: net) mapdist.push_back(mydistance(c1,n2.where));
-      sort(mapdist.begin(), mapdist.end());
-      maxdist = mapdist[size(mapdist)*5/6] * distmul;
-      printf("maxdist = %lf\n", maxdist);
-      }
-      
-    dispersion_count = 0;  
-    cell *c1 = currentmap->gamestart();
-    cell *c2 = createMov(c1, 0);
-    buildcellcrawler(c1);
-    if(c1->type != c2->type) buildcellcrawler(c2);
-  
-    lpct = -46130;
-    }
-  }
-
-void describe(cell *c) {
-  if(cmode & sm::HELP) return;
-  neuron *n = getNeuronSlow(c);
-  if(!n) return;
-  help += "cell number: " + its(neuronId(*n)) + " (" + its(n->samples) + ")\n";
-  help += "parameters:"; for(int k=0; k<cols; k++) help += " " + fts(n->net[k]); 
-  help += ", u-matrix = " + fts(n->udist);
-  help += "\n";
-  vector<pair<double, int>> v;
-  for(int s=0; s<samples; s++) if(whowon[s] == n) v.emplace_back(vnorm(n->net, data[s].val), s);
-  random_shuffle(v.begin(), v.end());
-  sort(v.begin(), v.end(), [] (auto a, auto b) { return a.first < b.first; });
-  
-  for(int i=0; i<size(v) && i<20; i++) {
-    int s = v[i].second;
-    help += "sample "+its(s)+":"; 
-    for(int k=0; k<cols; k++) help += " " + fts(data[s].val[k]); 
-    help += " "; help += data[s].name; help += "\n";
-    }
-  }
-
-namespace levelline {
-
-  struct levelline {
-    int column, qty;
-    unsigned int color;
-    vector<double> values;
-    bool modified;
-    };
-  
-  vector<levelline> levellines;
-  
-  bool on;
-  
-  void create() {
-    int xlalpha = int(pow(ld(.5), ggamma) * 255);
-    for(int i=0; i<cols; i++) {
-      levellines.emplace_back();
-      levelline& lv = levellines.back();
-      lv.column = i;    
-      lv.color = ((hrandpos() & 0xFFFFFF) << 8) | xlalpha;
-      lv.qty = 0;
-      }
-    }
-  
-  void build() {
-    if(levellines.size() == 0) create();
-    on = false;
-    for(auto& lv: levellines) {
-      if(!lv.qty) { lv.values.clear(); continue; }
-      on = true;
-      if(!lv.modified) continue;
-      lv.modified = false;
-      vector<double> sample;
-      for(int j=0; j<=1024; j++) sample.push_back(data[hrand(samples)].val[lv.column]);
-      sort(sample.begin(), sample.end());
-      lv.values.clear();
-      lv.values.push_back(-1e10);
-      for(int j=0; j<=1024; j+=1024 >> (lv.qty)) lv.values.push_back(sample[j]);
-      lv.values.push_back(1e10);
-      }
-    }
-  
-  void draw() {
-    if(!on) return;
-    for(auto& g: gmatrix) {
-      cell *c1 = g.first;
-      transmatrix T = g.second;
-      neuron *n1 = getNeuron(c1);
-      if(!n1) continue;
-      for(int i=0; i<c1->type; i++) {
-        cell *c2 = c1->mov[i];
-        if(!c2) continue;
-        cell *c3 = c1->mov[i ? i-1 : c1->type-1];
-        if(!c3) continue;
-  
-        if(!gmatrix.count(c2)) continue;
-        if(!gmatrix.count(c3)) continue;
-        double d2 = hdist(tC0(T), tC0(gmatrix[c2]));
-        double d3 = hdist(tC0(T), tC0(gmatrix[c3]));
-        
-        neuron *n2 = getNeuron(c2);
-        if(!n2) continue;
-        neuron *n3 = getNeuron(c3);
-        if(!n3) continue;
-              
-        for(auto& l: levellines) {
-          auto val1 = n1->net[l.column];
-          auto val2 = n2->net[l.column];
-          auto val3 = n3->net[l.column];
-          auto v1 = lower_bound(l.values.begin(), l.values.end(), val1);
-          auto v2 = lower_bound(l.values.begin(), l.values.end(), val2);
-          auto v3 = lower_bound(l.values.begin(), l.values.end(), val3);
-          auto draw = [&] () {
-            auto vmid = *v1;
-            queueline(
-              tC0(T * ddspin(c1,i) * xpush(d2 * (vmid-val1) / (val2-val1))), 
-              tC0(T * ddspin(c1,i-1) * xpush(d3 * (vmid-val1) / (val3-val1))), 
-              l.color);
-            };
-          while(v1 < v2 && v1 < v3) {
-            draw();
-            v1++;
-            }
-          while(v1 > v2 && v1 > v3) {
-            v1--;
-            draw();
-            }
-          }
-        }      
-      }
-    setindex(false);
-    }
-
-  void show() {
-    if(levellines.size() == 0) create();
-    gamescreen(0);
-    cmode = vid.xres > vid.yres * 1.4 ? sm::SIDE : 0;
-    dialog::init("level lines");
-    char nx = 'a';
-    for(auto &l : levellines) {
-      dialog::addSelItem(colnames[l.column], its(l.qty), nx++);
-      dialog::lastItem().colorv = l.color >> 8;
-      }
-    dialog::addItem("exit menu", '0');
-    dialog::addItem("shift+letter to change color", 0);
-    dialog::display();
-    keyhandler = [] (int sym, int uni) {
-      dialog::handleNavigation(sym, uni);
-      if(uni >= 'a' && uni - 'a' + size(levellines)) {
-        auto& l = levellines[uni - 'a'];
-        dialog::editNumber(l.qty, 0, 10, 1, 0, colnames[l.column], 
-          XLAT("Controls the number of level lines."));
-        dialog::reaction = [&l] () {
-          l.modified = true;
-          build();
-          };
-        }
-      else if(uni >= 'A' && uni - 'A' + size(levellines)) {
-        auto& l = levellines[uni - 'A'];
-        dialog::openColorDialog(l.color, NULL);
-        }
-      else if(doexiton(sym, uni)) popScreen();
-      };
-    }
-
-  
-  }
-
-void ksave(const char *fname) {
-  sominit(1);
-  FILE *f = fopen(fname, "wt");
-  if(!f) {
-    fprintf(stderr, "Could not save the network\n");
-    return;
-    }
-  fprintf(f, "%d %d\n", cells, t);
-  for(neuron& n: net) {
-    for(int k=0; k<cols; k++)
-      fprintf(f, "%.9lf ", n.net[k]);
-    fprintf(f, "\n");
-    }
-  fclose(f);
-  }
-
-void kload(const char *fname) {
-  sominit(1);
-  int xcells;
-  FILE *f = fopen(fname, "rt");
-  if(!f) {
-    fprintf(stderr, "Could not load the network\n");
-    return;
-    }
-  if(fscanf(f, "%d%d\n", &xcells, &t) != 2) return;
-  if(xcells != cells) {
-    fprintf(stderr, "Error: bad number of cells\n");
-    exit(1);
-    }
-  for(neuron& n: net) {
-    for(int k=0; k<cols; k++) if(fscanf(f, "%lf", &n.net[k]) != 1) return;
-    }
-  fclose(f);
-  analyze();
-  }
-
-void ksavew(const char *fname) {
-  sominit(1);
-  FILE *f = fopen(fname, "wt");
-  if(!f) {
-    fprintf(stderr, "Could not save the weights\n");
-    return;
-    }
-  for(int i=0; i<cols; i++)
-    fprintf(f, "%s=%.9lf\n", colnames[i].c_str(), weights[i]);
-  fclose(f);
-  }
-
-void kloadw(const char *fname) {
-  sominit(1);
-  FILE *f = fopen(fname, "rt");
-  if(!f) {
-    fprintf(stderr, "Could not load the weights\n");
-    return;
-    }
-  for(int i=0; i<cols; i++) {
-    string s1, s2;
-    char kind = 0;
-    while(true) {
-      int c = fgetc(f);
-      if(c == 10 || c == 13 || c == -1) {
-        if(s1 == "" && !kind && c != -1) continue;
-        if(s1 != "") colnames[i] = s1;
-        if(kind == '=') weights[i] = atof(s2.c_str());
-        if(kind == '*') weights[i] *= atof(s2.c_str());
-        if(kind == '/') weights[i] /= atof(s2.c_str());
-        if(c == -1) break;
-        goto nexti;
-        }
-      else if(c == '=' || c == '/' || c == '*') kind = c;
-      else (kind?s2:s1) += c;
-      }
-    nexti: ;
-    }
-  fclose(f);
-  analyze();
-  }
-
-
-unsigned lastprogress;
-void progress(string s) {
-  if(SDL_GetTicks() >= lastprogress + (noGUI ? 500 : 100)) {
-    if(noGUI)
-      printf("%s\n", s.c_str());
-    else {
-      clearMessages();
-      addMessage(s);
-      mainloopiter();
-      }
-    lastprogress = SDL_GetTicks();
-    } 
-  }
-
-void kclassify(const char *fname_classify) {
-
-  sominit(1);
-  vector<double> bdiffs(samples, 1e20);
-  vector<int> bids(samples, 0);
-
-  printf("Classifying...\n");
-
-  for(neuron& n: net) n.samples = 0;
-
-  for(int s=0; s<samples; s++) {
-    for(int n=0; n<cells; n++) {
-      double diff = vnorm(net[n].net, data[s].val);
-      if(diff < bdiffs[s]) bdiffs[s] = diff, bids[s] = n, whowon[s] = &net[n];
-      }
-    if(!(s % 128))
-      progress("Classifying: " + its(s) + "/" + its(samples));
-    }
-  
-  vector<double> bdiffn(cells, 1e20);
-
-  printf("Finding samples...\n");
-
-  for(int s=0; s<samples; s++) {
-    int n = bids[s];
-    double diff = bdiffs[s];
-    if(diff < bdiffn[n]) bdiffn[n] = diff, net[n].bestsample = s;
-    }
-
-  for(int s=0; s<samples; s++) net[bids[s]].samples++;
-  
-  if(fname_classify != NULL) {  
-    printf("Listing classification...\n");  
-    FILE *f = fopen(fname_classify, "wt");  
-    if(!f) {
-      printf("Failed to open file\n");
-      }
-    else {
-      for(int s=0; s<samples; s++)
-        fprintf(f, "%s;%d\n", data[s].name.c_str(), bids[s]);
-      fclose(f);
-      }
-    }
-  coloring();
-  }
-
-void klistsamples(const char *fname_samples, bool best, bool colorformat) {
-  if(fname_samples != NULL) { 
-    printf("Listing samples...\n");
-    FILE *f = fopen(fname_samples, "wt");  
-    if(!f) {
-      printf("Failed to open file\n");
-      }
-    else {
-      auto klistsample = [f, colorformat] (int id, int neu) {
-        if(colorformat) {
-          fprintf(f, "%s;+#%d\n", data[id].name.c_str(), neu);
-          }
-        else {
-          for(int k=0; k<cols; k++)
-            fprintf(f, "%.4lf ", data[id].val[k]);
-          fprintf(f, "!%s\n", data[id].name.c_str());
-          }
-        };
-      if(!colorformat) fprintf(f, "%d\n", cols);
-      if(best)
-        for(int n=0; n<cells; n++) {
-          if(!net[n].samples) { if(!colorformat) fprintf(f, "\n"); continue; }
-          klistsample(net[n].bestsample, n);
-          }
-      else
-        for(int i=0; i<size(samples_shown); i++) 
-          klistsample(samples_shown[i], neuronId(*(whowon[i])));  
-      fclose(f);
-      }
-    }
-  }
-
-void neurondisttable(const char *fname) {
-  FILE *f = fopen(fname, "wt");
-  if(!f) {
-    printf("Could not open file: %s\n", fname);
-    return;
-    }
-  int neurons = size(net);
-  fprintf(f, "%d\n", neurons);
-  for(int i=0; i<neurons; i++) {
-    for(int j=0; j<neurons; j++) fprintf(f, "%3d", celldistance(net[i].where, net[j].where));
-    // todo: build the table correctly for gaussian=2
-    fprintf(f, "\n");
-    }
-  fclose(f);
-  }
-
-void steps() {
-  if(!kohonen::finished()) {
-    unsigned int t = SDL_GetTicks();
-    while(SDL_GetTicks() < t+20) kohonen::step();
-    setindex(false);
-    }
-  }
-
-void showMenu() {
-  string parts[3] = {"red", "green", "blue"};
-  for(int i=0; i<3; i++) {
-    string c;
-    if(whattodraw[i] == -1) c = "u-matrix";
-    else if(whattodraw[i] == -2) c = "u-matrix reversed";
-    else if(whattodraw[i] == -3) c = "distance from marked ('m')";
-    else if(whattodraw[i] == -4) c = "number of samples";
-    else if(whattodraw[i] == -5) c = "best sample's color";
-    else if(whattodraw[i] == -6) c = "sample names to colors";
-    else c = colnames[whattodraw[i]];
-    dialog::addSelItem(XLAT("coloring (%1)", parts[i]), c, '1'+i);
-    }
-  dialog::addItem("coloring (all)", '0');
-  dialog::addItem("level lines", '4');
-  }
-
-bool handleMenu(int sym, int uni) {
-  if(uni >= '1' && uni <= '3') {
-    int i = uni - '1';
-    whattodraw[i]++;
-    if(whattodraw[i] == cols) whattodraw[i] = -5;
-    coloring();
-    return true;
-    }
-  if(uni == '0') {
-    for(char x: {'1','2','3'}) handleMenu(x, x);
-    return true;
-    }
-  if(uni == '4') {
-    pushScreen(levelline::show);
-    return true;
-    }
-  return false;
-  }
-
-int readArgs() {
-  using namespace arg;
-
-  // #1: load the samples
-  
-  if(argis("-som")) {
-    PHASE(3);
-    shift(); kohonen::loadsamples(args());
-    }
-
-  // #2: set parameters
-
-  else if(argis("-somkrad")) {
-    gaussian = 0; uninit(0);
-    }
-  else if(argis("-somsim")) {
-    gaussian = 0; uninit(1);
-    }
-  else if(argis("-somcgauss")) {
-    gaussian = 1; uninit(1);
-    }
-  else if(argis("-somggauss")) {
-    gaussian = 2; uninit(1);
-    }
-  else if(argis("-sompct")) {
-    shift(); qpct = argi();
-    }
-  else if(argis("-sompower")) {
-    shift(); ttpower = argf();
-    }
-  else if(argis("-somparam")) {
-    shift(); (gaussian ? distmul : dispersion_end_at) = argf();
-    if(dispersion_end_at <= 1) {
-      fprintf(stderr, "Dispersion parameter illegal\n");
-      dispersion_end_at = 1.5;
-      }
-    uninit(1);
-    }
-  else if(argis("-sominitdiv")) {
-    shift(); initdiv = argi(); uninit(0);
-    }
-  else if(argis("-somtmax")) {
-    shift(); t = (t*1./tmax) * argi();
-    tmax = argi();
-    }
-  else if(argis("-somlearn")) {
-    // this one can be changed at any moment
-    shift(); learning_factor = argf();
-    }
-
-  else if(argis("-somrun")) {
-    t = tmax; sominit(1);
-    }
-
-  // #3: load the neuron data (usually without #2)
-  else if(argis("-somload")) {
-    PHASE(3);
-    shift(); kohonen::kload(args());
-    }
-
-  // #4: run, stop etc.
-  else if(argis("-somrunto")) {
-    int i = argi();
-    shift(); while(t > i) {
-      if(t % 128 == 0) progress("Steps left: " + its(t));
-      kohonen::step();
-      }
-    }
-  else if(argis("-somstop")) {
-    t = 0;
-    }
-  else if(argis("-somnoshow")) {
-    noshow = true;
-    }
-  else if(argis("-somfinish")) {
-    while(!finished()) {
-      kohonen::step();
-      if(t % 128 == 0) progress("Steps left: " + its(t));
-      }
-    }
-
-  // #5 save data, classify etc.
-  else if(argis("-somsave")) {
-    PHASE(3);
-    shift(); kohonen::ksave(args());
-    }
-  else if(argis("-somsavew")) {
-    PHASE(3);
-    shift(); kohonen::ksavew(args());
-    }
-  else if(argis("-somloadw")) {
-    PHASE(3);
-    shift(); kohonen::kloadw(args());
-    }
-  else if(argis("-somclassify")) {
-    PHASE(3);
-    shift(); kohonen::kclassify(args());
-    }
-  else if(argis("-somlistshown")) {
-    PHASE(3);
-    shift(); kohonen::klistsamples(args(), false, false);
-    }
-  else if(argis("-somlistbest")) {
-    PHASE(3);
-    shift(); kohonen::klistsamples(args(), true, false);
-    }
-  else if(argis("-somlistbestc")) {
-    PHASE(3);
-    shift(); kohonen::klistsamples(args(), true, true);
-    }
-  else if(argis("-somndist")) {
-    PHASE(3);
-    shift(); kohonen::neurondisttable(args());
-    }
-  else if(argis("-somshowbest")) {
-    showbestsamples();
-    }
-
-  else return 1;
-  return 0;
-  }
-
-auto hooks = addHook(hooks_args, 100, readArgs) + addHook(hooks_frame, 50, levelline::draw);
-}
-
-void mark(cell *c) {
-  using namespace kohonen;
-  distfrom = getNeuronSlow(c);
-  coloring();
-  }
-
+// Hyperbolic Rogue
+// Copyright (C) 2011-2017 Zeno and Tehora Rogue, see 'hyper.cpp' for details
+
+// Kohonen's self-organizing networks. 
+// This is a part of RogueViz, not a part of HyperRogue.
+
+namespace kohonen {
+
+int cols;
+
+typedef vector<double> kohvec;
+
+struct sample {
+  kohvec val;
+  string name;
+  };
+
+vector<sample> data;
+
+vector<int> samples_shown;
+
+int whattodraw[3] = {-2,-2,-2};
+
+struct neuron {
+  kohvec net;
+  cell *where;
+  double udist;
+  int lpbak;
+  int col;
+  int samples, csample, bestsample;
+  };
+
+vector<string> colnames;
+
+kohvec weights;
+
+vector<neuron> net;
+
+int neuronId(neuron& n) { return &n - &(net[0]); }
+
+void alloc(kohvec& k) { k.resize(cols); }
+
+bool neurons_indexed = false;
+
+int samples;
+
+template<class T> T sqr(T x) { return x*x; }
+
+vector<neuron*> whowon;
+
+void normalize() {
+  alloc(weights);
+  for(int k=0; k<cols; k++) {
+    double sum = 0, sqsum = 0;
+    for(sample& s: data)
+      sum += s.val[k],
+      sqsum += s.val[k] * s.val[k];
+    double variance = sqsum/samples - sqr(sum/samples);
+    weights[k] = 1 / sqrt(variance);
+    }
+  }
+
+double vnorm(kohvec& a, kohvec& b) {
+  double diff = 0;
+  for(int k=0; k<cols; k++) diff += sqr((a[k]-b[k]) * weights[k]);
+  return diff;
+  }
+
+void sominit(int);
+void uninit(int);
+
+void loadsamples(const char *fname) {
+  FILE *f = fopen(fname, "rt");
+  if(!f) {
+    fprintf(stderr, "Could not load samples\n");
+    return;
+    }
+  if(fscanf(f, "%d", &cols) != 1) { fclose(f); return; }
+  while(true) {
+    sample s;
+    bool shown = false;
+    alloc(s.val);
+    if(feof(f)) break;
+    for(int i=0; i<cols; i++)
+      if(fscanf(f, "%lf", &s.val[i]) != 1) { goto bigbreak; }
+    fgetc(f);
+    while(true) {
+      int c = fgetc(f);
+      if(c == -1 || c == 10 || c == 13) break;
+      if(c == '!' && s.name == "") shown = true;
+      else if(c != 32 && c != 9) s.name += c;
+      }
+    if(shown) samples_shown.push_back(size(data));
+    data.push_back(move(s));
+    }
+  bigbreak:
+  fclose(f);
+  samples = size(data);
+  normalize();
+  colnames.resize(cols);
+  for(int i=0; i<cols; i++) colnames[i] = "Column " + its(i);
+  uninit(0); sominit(1);
+  }
+                  
+int tmax = 30000;
+double distmul = 1;
+double learning_factor = .1;
+
+int qpct = 100;
+
+int t, lpct, cells;
+double maxdist;
+
+neuron& winner(int id) {
+  double bdiff = 1e20;
+  neuron *bcell = NULL;
+  for(neuron& n: net) {
+    double diff = vnorm(n.net, data[id].val);
+    if(diff < bdiff) bdiff = diff, bcell = &n;
+    }
+  return *bcell;
+  }
+
+void setindex(bool b) {
+  if(b == neurons_indexed) return;
+  neurons_indexed = b;
+  if(b) {
+    for(neuron& n: net) n.lpbak = n.where->landparam, n.where->landparam = neuronId(n);
+    }
+  else {
+    for(neuron& n: net) n.where->landparam = n.lpbak;
+    }
+  }
+
+neuron *getNeuron(cell *c) {
+  if(!c) return NULL;
+  setindex(true);
+  if(c->landparam < 0 || c->landparam >= cells) return NULL;
+  neuron& ret = net[c->landparam];
+  if(ret.where != c) return NULL;
+  return &ret;
+  }
+
+neuron *getNeuronSlow(cell *c) {
+  if(neurons_indexed) return getNeuron(c);
+  for(neuron& n: net) if(n.where == c) return &n;
+  return NULL;
+  }
+
+double maxudist;
+
+neuron *distfrom;
+
+bool noshow = false;
+
+void coloring() {
+  if(noshow) return;
+  setindex(false);
+  
+  bool besttofind = true;
+
+  for(int pid=0; pid<3; pid++) {
+    int c = whattodraw[pid];
+    
+    if(c == -5) {
+      if(besttofind) {
+        besttofind = false;
+        for(neuron& n: net) {
+          double bdiff = 1e20;
+          for(int i=0; i<size(samples_shown); i++) {
+            double diff = vnorm(n.net, data[samples_shown[i]].val);
+            if(diff < bdiff) bdiff = diff, n.bestsample = i;
+            }
+          }
+        }
+
+      for(int i=0; i<cells; i++) 
+        part(net[i].where->landparam, pid) = part(vdata[net[i].bestsample].cp.color1, pid+1);
+      }
+
+    else {
+      vector<double> listing;
+      for(neuron& n: net) switch(c) {
+        case -4:
+          listing.push_back(log(5+n.samples));
+          break;
+          
+        case -3:
+          if(distfrom) 
+            listing.push_back(vnorm(n.net, distfrom->net));
+          else
+            listing.push_back(0);
+          break;
+        
+        case -2:
+          listing.push_back(n.udist);
+          break;
+        
+        case -1: 
+          listing.push_back(-n.udist);
+          break;
+        
+        default:
+          listing.push_back(n.net[c]);
+          break;
+        }
+      
+      double minl = listing[0], maxl = listing[0];
+      for(double& d: listing) minl = min(minl, d), maxl = max(maxl, d);
+      if(maxl-minl < 1e-3) maxl = minl+1e-3;
+      
+      for(int i=0; i<cells; i++) 
+        part(net[i].where->landparam, pid) = (255 * (listing[i] - minl)) / (maxl - minl);
+      }
+    }
+  }
+
+void analyze() {
+
+  setindex(true);
+  
+  maxudist = 0;
+  for(neuron& n: net) {
+    int qty = 0;
+    double total = 0;
+    forCellEx(c2, n.where) {
+      neuron *n2 = getNeuron(c2);
+      if(!n2) continue;
+      qty++;
+      total += sqrt(vnorm(n.net, n2->net));
+      }
+    n.udist = total / qty;
+    maxudist = max(maxudist, n.udist);
+    }
+    
+  if(!noshow) {
+
+    whowon.resize(samples);
+    
+    for(neuron& n: net) n.samples = 0;
+    
+    for(int id=0; id<size(samples_shown); id++) {
+      int s = samples_shown[id];
+      auto& w = winner(s);
+      whowon[s] = &w;
+      w.samples++;
+      }
+    
+    for(int id=0; id<size(samples_shown); id++) {
+      int s = samples_shown[id];
+      auto& w = *whowon[s];
+      vdata[id].m->base = w.where;
+      vdata[id].m->at = 
+        spin(2*M_PI*w.csample / w.samples) * xpush(.25 * (w.samples-1) / w.samples);
+      w.csample++;
+      }
+    
+    shmup::fixStorage();  
+    setindex(false);
+    }
+  
+  coloring();
+  }
+
+// traditionally Gaussian blur is used in the Kohonen algoritm
+// but it does not seem to make much sense in hyperbolic geometry
+// especially wrapped one.
+// GAUSSIAN==1: use the Gaussian blur, on celldistance
+// GAUSSIAN==2: use the Gaussian blur, on true distance
+// GAUSSIAN==0: simulate the dispersion on our network
+
+int gaussian = 0;
+
+double mydistance(cell *c1, cell *c2) {
+  if(gaussian == 2) return hdist(tC0(shmup::ggmatrix(c1)), tC0(shmup::ggmatrix(c2)));
+  else return celldistance(c1, c2);
+  }
+
+struct cellcrawler {
+
+  struct cellcrawlerdata {
+    cellwalker orig;
+    int from, spin, dist;
+    cellwalker target;
+    cellcrawlerdata(const cellwalker& o, int fr, int sp) : orig(o), from(fr), spin(sp) {}
+    };
+  
+  vector<cellcrawlerdata> data;
+  
+  void store(const cellwalker& o, int from, int spin) {
+    if(eq(o.c->aitmp, sval)) return;
+    o.c->aitmp = sval;
+    data.emplace_back(o, from, spin);
+    }
+  
+  void build(const cellwalker& start) {
+    sval++;
+    data.clear();
+    store(start, 0, 0);
+    for(int i=0; i<size(data); i++) {
+      cellwalker cw0 = data[i].orig;
+      for(int j=0; j<cw0.c->type; j++) {
+        cellwalker cw = cw0;
+        cwspin(cw, j); cwstep(cw);
+        if(!getNeuron(cw.c)) continue;
+        store(cw, i, j);
+        }
+      }
+    if(gaussian) for(cellcrawlerdata& s: data)
+      s.dist = mydistance(s.orig.c, start.c);
+    }
+  
+  void sprawl(const cellwalker& start) {
+    data[0].target = start;
+    
+    for(int i=1; i<size(data); i++) {
+      cellcrawlerdata& s = data[i];
+      s.target = data[s.from].target;
+      if(!s.target.c) continue;
+      cwspin(s.target, s.spin);
+      if(cwstepcreates(s.target)) s.target.c = NULL;
+      else cwstep(s.target);
+      }
+    }
+  };
+
+cellcrawler scc[2]; // hex and non-hex
+
+double dispersion_end_at = 1.5;
+
+double dispersion_precision = .0001;
+int dispersion_each = 1;
+
+vector<vector<ld>> dispersion[2]; 
+
+int dispersion_count;
+
+void buildcellcrawler(cell *c) {
+  int sccid = c->type != 6;
+  
+  cellcrawler& cr =  scc[sccid];
+  cr.build(cellwalker(c,0));
+
+  if(!gaussian) {
+    vector<ld> curtemp;
+    vector<ld> newtemp;
+    vector<int> qty;
+    vector<pair<ld*, ld*> > pairs;
+    int N = size(net);
+    
+    curtemp.resize(N, 0);
+    newtemp.resize(N, 0);
+    qty.resize(N, 0);
+  
+    for(int i=0; i<N; i++) 
+    forCellEx(c2, net[i].where) {
+      neuron *nj = getNeuron(c2);
+      if(nj) {
+        pairs.emplace_back(&curtemp[i], &newtemp[neuronId(*nj)]);
+        qty[i]++;
+        }
+      }
+    
+    curtemp[neuronId(*getNeuron(c))] = 1;
+  
+    ld vmin = 0, vmax = 1;
+    int iter;
+    
+    auto &d = dispersion[sccid];
+    
+    d.clear();
+  
+    printf("Building dispersion...\n");
+    
+    for(iter=0; dispersion_count ? true : vmax > vmin * dispersion_end_at; iter++) {
+      if(iter % dispersion_each == 0) {
+        d.emplace_back(N);
+        auto& dispvec = d.back();
+        for(int i=0; i<N; i++) dispvec[i] = curtemp[neuronId(*getNeuron(cr.data[i].orig.c))] / vmax;
+        if(size(d) == dispersion_count) break;
+        }
+      double df = dispersion_precision * (iter+1);
+      double df0 = df / ceil(df);
+      for(int i=0; i<df; i++) {
+        for(auto& p: pairs) 
+          *p.second += *p.first;
+        for(int i=0; i<N; i++) {
+          curtemp[i] += (newtemp[i] / qty[i] - curtemp[i]) * df0;
+          newtemp[i] = 0;
+          }
+        }
+      vmin = vmax = curtemp[0];
+      for(int i=0; i<N; i++) 
+        if(curtemp[i] < vmin) vmin = curtemp[i];
+        else if(curtemp[i] > vmax) vmax = curtemp[i];
+      }
+  
+    dispersion_count = size(d);
+    printf("Dispersion count = %d\n", dispersion_count);
+    }
+  }
+
+bool finished() { return t == 0; }
+
+int krad;
+
+double ttpower = 1;
+
+void sominit(int);
+
+void step() {
+
+  if(t == 0) return;
+  sominit(2);
+  
+  double tt = (t-1.) / tmax;
+  tt = pow(tt, ttpower);
+
+  double sigma = maxdist * tt;
+  int dispid = int(dispersion_count * tt);
+
+  if(qpct) {
+    int pct = (int) ((qpct * (t+.0)) / tmax);
+    if(pct != lpct) {
+      printf("pct %d lpct %d\n", pct, lpct);
+      lpct = pct;
+      analyze();
+      
+      if(gaussian)
+        printf("t = %6d/%6d %3d%% sigma=%10.7lf maxudist=%10.7lf\n", t, tmax, pct, sigma, maxudist);
+      else
+        printf("t = %6d/%6d %3d%% dispid=%5d maxudist=%10.7lf\n", t, tmax, pct, dispid, maxudist);
+      }
+    }
+  int id = hrand(samples);
+  neuron& n = winner(id);
+  whowon.resize(samples);
+  whowon[id] = &n;
+    
+  /* 
+  for(neuron& n2: net) {
+    int d = celldistance(n.where, n2.where);
+    double nu = learning_factor; 
+//  nu *= exp(-t*(double)maxdist/perdist);
+//  nu *= exp(-t/t2);
+    nu *= exp(-sqr(d/sigma));
+    for(int k=0; k<cols; k++)
+      n2.net[k] += nu * (irisdata[id][k] - n2.net[k]);
+    } */
+    
+  int sccid = n.where->type != 6;
+  cellcrawler& s = scc[sccid];
+  s.sprawl(cellwalker(n.where, 0));
+
+  vector<double> fake(1,1);
+  auto it = gaussian ? fake.begin() : dispersion[sccid][dispid].begin();
+
+  for(auto& sd: s.data) {
+    neuron *n2 = getNeuron(sd.target.c);
+    if(!n2) continue;
+    double nu = learning_factor;
+    
+    if(gaussian)
+      nu *= exp(-sqr(sd.dist/sigma));
+    else
+      nu *= *(it++);
+
+    for(int k=0; k<cols; k++)
+      n2->net[k] += nu * (data[id].val[k] - n2->net[k]);
+    }
+  
+  t--;
+  if(t == 0) analyze();
+  }
+
+int initdiv = 1;
+
+int inited = 0;
+
+void uninit(int initto) {
+  if(inited > initto) inited = initto;
+  }
+
+void showsample(int id) {
+  for(int ii: samples_shown) 
+    if(ii == id) 
+      return;
+  int i = vdata.size();
+  samples_shown.push_back(id);
+  vdata.emplace_back();
+  auto& v = vdata.back();
+  v.name = data[id].name;
+  v.cp = dftcolor;
+  createViz(i, cwt.c, Id);
+  v.m->store();
+  }
+
+void showsample(string s) {
+  if(s == "") return;
+  for(int i=0; i<samples; i++) {
+    if(s[0] != '*' && data[i].name == s)
+      showsample(i);
+    if(s[0] == '*' && data[i].name.find(s.substr(1)) != string::npos)
+      showsample(i);
+    }
+  }
+
+void showbestsamples() {
+  vector<int> samplesbak;
+  for(auto& n: net) 
+    if(n.samples)
+      showsample(n.bestsample);
+  analyze();
+  for(auto& n: net) n.samples = 0;
+  for(int i=0; i<samples; i++) 
+    if(whowon[i])
+      whowon[i]->samples++;
+  }
+  
+void sominit(int initto) {
+
+  if(inited < 1 && initto >= 1) {
+    inited = 1;
+    if(!samples) {
+      fprintf(stderr, "Error: SOM without samples\n");
+      exit(1);
+      }
+    
+    init(); kind = kKohonen;
+    
+    /* if(geometry != gQuotient1) {
+      targetGeometry = gQuotient1;
+      restartGame('g');
+      }
+    if(!purehepta) restartGame('7'); */
+    
+    printf("Initializing SOM (1)\n");
+  
+    vector<cell*> allcells;
+    
+    if(krad) {
+      celllister cl(cwt.c, krad, 1000000, NULL);
+      allcells = cl.lst;
+      }
+    else allcells = currentmap->allcells();
+  
+    cells = size(allcells);
+    net.resize(cells);
+    for(int i=0; i<cells; i++) net[i].where = allcells[i], allcells[i]->landparam = i;
+    for(int i=0; i<cells; i++) {
+      net[i].where->land = laCanvas;
+      alloc(net[i].net);
+  
+      for(int k=0; k<cols; k++)
+      for(int z=0; z<initdiv; z++)
+        net[i].net[k] += data[hrand(samples)].val[k] / initdiv;
+      }
+      
+    for(neuron& n: net) for(int d=BARLEV; d>=7; d--) setdist(n.where, d, NULL);
+  
+    printf("samples = %d (%d) cells = %d\n", samples, size(samples_shown), cells);
+
+    if(!noshow) {
+      vdata.resize(size(samples_shown));
+      for(int i=0; i<size(samples_shown); i++) {
+        vdata[i].name = data[samples_shown[i]].name;
+        vdata[i].cp = dftcolor;
+        createViz(i, cwt.c, Id);
+        }
+      
+      storeall();
+      }
+
+    analyze();
+    }
+  
+  if(inited < 2 && initto >= 2) {
+    inited = 2;
+
+    printf("Initializing SOM (2)\n");
+
+    if(gaussian) {
+      printf("dist = %lf\n", mydistance(net[0].where, net[1].where));
+      cell *c1 = net[cells/2].where;
+      vector<double> mapdist;
+      for(neuron &n2: net) mapdist.push_back(mydistance(c1,n2.where));
+      sort(mapdist.begin(), mapdist.end());
+      maxdist = mapdist[size(mapdist)*5/6] * distmul;
+      printf("maxdist = %lf\n", maxdist);
+      }
+      
+    dispersion_count = 0;  
+    cell *c1 = currentmap->gamestart();
+    cell *c2 = createMov(c1, 0);
+    buildcellcrawler(c1);
+    if(c1->type != c2->type) buildcellcrawler(c2);
+  
+    lpct = -46130;
+    }
+  }
+
+void describe(cell *c) {
+  if(cmode & sm::HELP) return;
+  neuron *n = getNeuronSlow(c);
+  if(!n) return;
+  help += "cell number: " + its(neuronId(*n)) + " (" + its(n->samples) + ")\n";
+  help += "parameters:"; for(int k=0; k<cols; k++) help += " " + fts(n->net[k]); 
+  help += ", u-matrix = " + fts(n->udist);
+  help += "\n";
+  vector<pair<double, int>> v;
+  for(int s=0; s<samples; s++) if(whowon[s] == n) v.emplace_back(vnorm(n->net, data[s].val), s);
+  random_shuffle(v.begin(), v.end());
+  sort(v.begin(), v.end(), [] (auto a, auto b) { return a.first < b.first; });
+  
+  for(int i=0; i<size(v) && i<20; i++) {
+    int s = v[i].second;
+    help += "sample "+its(s)+":"; 
+    for(int k=0; k<cols; k++) help += " " + fts(data[s].val[k]); 
+    help += " "; help += data[s].name; help += "\n";
+    }
+  }
+
+namespace levelline {
+
+  struct levelline {
+    int column, qty;
+    unsigned int color;
+    vector<double> values;
+    bool modified;
+    };
+  
+  vector<levelline> levellines;
+  
+  bool on;
+  
+  void create() {
+    int xlalpha = int(pow(ld(.5), ggamma) * 255);
+    for(int i=0; i<cols; i++) {
+      levellines.emplace_back();
+      levelline& lv = levellines.back();
+      lv.column = i;    
+      lv.color = ((hrandpos() & 0xFFFFFF) << 8) | xlalpha;
+      lv.qty = 0;
+      }
+    }
+  
+  void build() {
+    if(levellines.size() == 0) create();
+    on = false;
+    for(auto& lv: levellines) {
+      if(!lv.qty) { lv.values.clear(); continue; }
+      on = true;
+      if(!lv.modified) continue;
+      lv.modified = false;
+      vector<double> sample;
+      for(int j=0; j<=1024; j++) sample.push_back(data[hrand(samples)].val[lv.column]);
+      sort(sample.begin(), sample.end());
+      lv.values.clear();
+      lv.values.push_back(-1e10);
+      for(int j=0; j<=1024; j+=1024 >> (lv.qty)) lv.values.push_back(sample[j]);
+      lv.values.push_back(1e10);
+      }
+    }
+  
+  void draw() {
+    if(!on) return;
+    for(auto& g: gmatrix) {
+      cell *c1 = g.first;
+      transmatrix T = g.second;
+      neuron *n1 = getNeuron(c1);
+      if(!n1) continue;
+      for(int i=0; i<c1->type; i++) {
+        cell *c2 = c1->mov[i];
+        if(!c2) continue;
+        cell *c3 = c1->mov[i ? i-1 : c1->type-1];
+        if(!c3) continue;
+  
+        if(!gmatrix.count(c2)) continue;
+        if(!gmatrix.count(c3)) continue;
+        double d2 = hdist(tC0(T), tC0(gmatrix[c2]));
+        double d3 = hdist(tC0(T), tC0(gmatrix[c3]));
+        
+        neuron *n2 = getNeuron(c2);
+        if(!n2) continue;
+        neuron *n3 = getNeuron(c3);
+        if(!n3) continue;
+              
+        for(auto& l: levellines) {
+          auto val1 = n1->net[l.column];
+          auto val2 = n2->net[l.column];
+          auto val3 = n3->net[l.column];
+          auto v1 = lower_bound(l.values.begin(), l.values.end(), val1);
+          auto v2 = lower_bound(l.values.begin(), l.values.end(), val2);
+          auto v3 = lower_bound(l.values.begin(), l.values.end(), val3);
+          auto draw = [&] () {
+            auto vmid = *v1;
+            queueline(
+              tC0(T * ddspin(c1,i) * xpush(d2 * (vmid-val1) / (val2-val1))), 
+              tC0(T * ddspin(c1,i-1) * xpush(d3 * (vmid-val1) / (val3-val1))), 
+              l.color);
+            };
+          while(v1 < v2 && v1 < v3) {
+            draw();
+            v1++;
+            }
+          while(v1 > v2 && v1 > v3) {
+            v1--;
+            draw();
+            }
+          }
+        }      
+      }
+    setindex(false);
+    }
+
+  void show() {
+    if(levellines.size() == 0) create();
+    gamescreen(0);
+    cmode = vid.xres > vid.yres * 1.4 ? sm::SIDE : 0;
+    dialog::init("level lines");
+    char nx = 'a';
+    for(auto &l : levellines) {
+      dialog::addSelItem(colnames[l.column], its(l.qty), nx++);
+      dialog::lastItem().colorv = l.color >> 8;
+      }
+    dialog::addItem("exit menu", '0');
+    dialog::addItem("shift+letter to change color", 0);
+    dialog::display();
+    keyhandler = [] (int sym, int uni) {
+      dialog::handleNavigation(sym, uni);
+      if(uni >= 'a' && uni - 'a' + size(levellines)) {
+        auto& l = levellines[uni - 'a'];
+        dialog::editNumber(l.qty, 0, 10, 1, 0, colnames[l.column], 
+          XLAT("Controls the number of level lines."));
+        dialog::reaction = [&l] () {
+          l.modified = true;
+          build();
+          };
+        }
+      else if(uni >= 'A' && uni - 'A' + size(levellines)) {
+        auto& l = levellines[uni - 'A'];
+        dialog::openColorDialog(l.color, NULL);
+        }
+      else if(doexiton(sym, uni)) popScreen();
+      };
+    }
+
+  
+  }
+
+void ksave(const char *fname) {
+  sominit(1);
+  FILE *f = fopen(fname, "wt");
+  if(!f) {
+    fprintf(stderr, "Could not save the network\n");
+    return;
+    }
+  fprintf(f, "%d %d\n", cells, t);
+  for(neuron& n: net) {
+    for(int k=0; k<cols; k++)
+      fprintf(f, "%.9lf ", n.net[k]);
+    fprintf(f, "\n");
+    }
+  fclose(f);
+  }
+
+void kload(const char *fname) {
+  sominit(1);
+  int xcells;
+  FILE *f = fopen(fname, "rt");
+  if(!f) {
+    fprintf(stderr, "Could not load the network\n");
+    return;
+    }
+  if(fscanf(f, "%d%d\n", &xcells, &t) != 2) return;
+  if(xcells != cells) {
+    fprintf(stderr, "Error: bad number of cells\n");
+    exit(1);
+    }
+  for(neuron& n: net) {
+    for(int k=0; k<cols; k++) if(fscanf(f, "%lf", &n.net[k]) != 1) return;
+    }
+  fclose(f);
+  analyze();
+  }
+
+void ksavew(const char *fname) {
+  sominit(1);
+  FILE *f = fopen(fname, "wt");
+  if(!f) {
+    fprintf(stderr, "Could not save the weights\n");
+    return;
+    }
+  for(int i=0; i<cols; i++)
+    fprintf(f, "%s=%.9lf\n", colnames[i].c_str(), weights[i]);
+  fclose(f);
+  }
+
+void kloadw(const char *fname) {
+  sominit(1);
+  FILE *f = fopen(fname, "rt");
+  if(!f) {
+    fprintf(stderr, "Could not load the weights\n");
+    return;
+    }
+  for(int i=0; i<cols; i++) {
+    string s1, s2;
+    char kind = 0;
+    while(true) {
+      int c = fgetc(f);
+      if(c == 10 || c == 13 || c == -1) {
+        if(s1 == "" && !kind && c != -1) continue;
+        if(s1 != "") colnames[i] = s1;
+        if(kind == '=') weights[i] = atof(s2.c_str());
+        if(kind == '*') weights[i] *= atof(s2.c_str());
+        if(kind == '/') weights[i] /= atof(s2.c_str());
+        if(c == -1) break;
+        goto nexti;
+        }
+      else if(c == '=' || c == '/' || c == '*') kind = c;
+      else (kind?s2:s1) += c;
+      }
+    nexti: ;
+    }
+  fclose(f);
+  analyze();
+  }
+
+
+unsigned lastprogress;
+void progress(string s) {
+  if(SDL_GetTicks() >= lastprogress + (noGUI ? 500 : 100)) {
+    if(noGUI)
+      printf("%s\n", s.c_str());
+    else {
+      clearMessages();
+      addMessage(s);
+      mainloopiter();
+      }
+    lastprogress = SDL_GetTicks();
+    } 
+  }
+
+void kclassify(const char *fname_classify) {
+
+  sominit(1);
+  vector<double> bdiffs(samples, 1e20);
+  vector<int> bids(samples, 0);
+
+  printf("Classifying...\n");
+
+  for(neuron& n: net) n.samples = 0;
+
+  for(int s=0; s<samples; s++) {
+    for(int n=0; n<cells; n++) {
+      double diff = vnorm(net[n].net, data[s].val);
+      if(diff < bdiffs[s]) bdiffs[s] = diff, bids[s] = n, whowon[s] = &net[n];
+      }
+    if(!(s % 128))
+      progress("Classifying: " + its(s) + "/" + its(samples));
+    }
+  
+  vector<double> bdiffn(cells, 1e20);
+
+  printf("Finding samples...\n");
+
+  for(int s=0; s<samples; s++) {
+    int n = bids[s];
+    double diff = bdiffs[s];
+    if(diff < bdiffn[n]) bdiffn[n] = diff, net[n].bestsample = s;
+    }
+
+  for(int s=0; s<samples; s++) net[bids[s]].samples++;
+  
+  if(fname_classify != NULL) {  
+    printf("Listing classification...\n");  
+    FILE *f = fopen(fname_classify, "wt");  
+    if(!f) {
+      printf("Failed to open file\n");
+      }
+    else {
+      for(int s=0; s<samples; s++)
+        fprintf(f, "%s;%d\n", data[s].name.c_str(), bids[s]);
+      fclose(f);
+      }
+    }
+  coloring();
+  }
+
+void klistsamples(const char *fname_samples, bool best, bool colorformat) {
+  if(fname_samples != NULL) { 
+    printf("Listing samples...\n");
+    FILE *f = fopen(fname_samples, "wt");  
+    if(!f) {
+      printf("Failed to open file\n");
+      }
+    else {
+      auto klistsample = [f, colorformat] (int id, int neu) {
+        if(colorformat) {
+          fprintf(f, "%s;+#%d\n", data[id].name.c_str(), neu);
+          }
+        else {
+          for(int k=0; k<cols; k++)
+            fprintf(f, "%.4lf ", data[id].val[k]);
+          fprintf(f, "!%s\n", data[id].name.c_str());
+          }
+        };
+      if(!colorformat) fprintf(f, "%d\n", cols);
+      if(best)
+        for(int n=0; n<cells; n++) {
+          if(!net[n].samples) { if(!colorformat) fprintf(f, "\n"); continue; }
+          klistsample(net[n].bestsample, n);
+          }
+      else
+        for(int i=0; i<size(samples_shown); i++) 
+          klistsample(samples_shown[i], neuronId(*(whowon[i])));  
+      fclose(f);
+      }
+    }
+  }
+
+void neurondisttable(const char *fname) {
+  FILE *f = fopen(fname, "wt");
+  if(!f) {
+    printf("Could not open file: %s\n", fname);
+    return;
+    }
+  int neurons = size(net);
+  fprintf(f, "%d\n", neurons);
+  for(int i=0; i<neurons; i++) {
+    for(int j=0; j<neurons; j++) fprintf(f, "%3d", celldistance(net[i].where, net[j].where));
+    // todo: build the table correctly for gaussian=2
+    fprintf(f, "\n");
+    }
+  fclose(f);
+  }
+
+void steps() {
+  if(!kohonen::finished()) {
+    unsigned int t = SDL_GetTicks();
+    while(SDL_GetTicks() < t+20) kohonen::step();
+    setindex(false);
+    }
+  }
+
+void showMenu() {
+  string parts[3] = {"red", "green", "blue"};
+  for(int i=0; i<3; i++) {
+    string c;
+    if(whattodraw[i] == -1) c = "u-matrix";
+    else if(whattodraw[i] == -2) c = "u-matrix reversed";
+    else if(whattodraw[i] == -3) c = "distance from marked ('m')";
+    else if(whattodraw[i] == -4) c = "number of samples";
+    else if(whattodraw[i] == -5) c = "best sample's color";
+    else if(whattodraw[i] == -6) c = "sample names to colors";
+    else c = colnames[whattodraw[i]];
+    dialog::addSelItem(XLAT("coloring (%1)", parts[i]), c, '1'+i);
+    }
+  dialog::addItem("coloring (all)", '0');
+  dialog::addItem("level lines", '4');
+  }
+
+bool handleMenu(int sym, int uni) {
+  if(uni >= '1' && uni <= '3') {
+    int i = uni - '1';
+    whattodraw[i]++;
+    if(whattodraw[i] == cols) whattodraw[i] = -5;
+    coloring();
+    return true;
+    }
+  if(uni == '0') {
+    for(char x: {'1','2','3'}) handleMenu(x, x);
+    return true;
+    }
+  if(uni == '4') {
+    pushScreen(levelline::show);
+    return true;
+    }
+  return false;
+  }
+
+int readArgs() {
+  using namespace arg;
+
+  // #1: load the samples
+  
+  if(argis("-som")) {
+    PHASE(3);
+    shift(); kohonen::loadsamples(args());
+    }
+
+  // #2: set parameters
+
+  else if(argis("-somkrad")) {
+    gaussian = 0; uninit(0);
+    }
+  else if(argis("-somsim")) {
+    gaussian = 0; uninit(1);
+    }
+  else if(argis("-somcgauss")) {
+    gaussian = 1; uninit(1);
+    }
+  else if(argis("-somggauss")) {
+    gaussian = 2; uninit(1);
+    }
+  else if(argis("-sompct")) {
+    shift(); qpct = argi();
+    }
+  else if(argis("-sompower")) {
+    shift(); ttpower = argf();
+    }
+  else if(argis("-somparam")) {
+    shift(); (gaussian ? distmul : dispersion_end_at) = argf();
+    if(dispersion_end_at <= 1) {
+      fprintf(stderr, "Dispersion parameter illegal\n");
+      dispersion_end_at = 1.5;
+      }
+    uninit(1);
+    }
+  else if(argis("-sominitdiv")) {
+    shift(); initdiv = argi(); uninit(0);
+    }
+  else if(argis("-somtmax")) {
+    shift(); t = (t*1./tmax) * argi();
+    tmax = argi();
+    }
+  else if(argis("-somlearn")) {
+    // this one can be changed at any moment
+    shift(); learning_factor = argf();
+    }
+
+  else if(argis("-somrun")) {
+    t = tmax; sominit(1);
+    }
+
+  // #3: load the neuron data (usually without #2)
+  else if(argis("-somload")) {
+    PHASE(3);
+    shift(); kohonen::kload(args());
+    }
+
+  // #4: run, stop etc.
+  else if(argis("-somrunto")) {
+    int i = argi();
+    shift(); while(t > i) {
+      if(t % 128 == 0) progress("Steps left: " + its(t));
+      kohonen::step();
+      }
+    }
+  else if(argis("-somstop")) {
+    t = 0;
+    }
+  else if(argis("-somnoshow")) {
+    noshow = true;
+    }
+  else if(argis("-somfinish")) {
+    while(!finished()) {
+      kohonen::step();
+      if(t % 128 == 0) progress("Steps left: " + its(t));
+      }
+    }
+
+  // #5 save data, classify etc.
+  else if(argis("-somsave")) {
+    PHASE(3);
+    shift(); kohonen::ksave(args());
+    }
+  else if(argis("-somsavew")) {
+    PHASE(3);
+    shift(); kohonen::ksavew(args());
+    }
+  else if(argis("-somloadw")) {
+    PHASE(3);
+    shift(); kohonen::kloadw(args());
+    }
+  else if(argis("-somclassify")) {
+    PHASE(3);
+    shift(); kohonen::kclassify(args());
+    }
+  else if(argis("-somlistshown")) {
+    PHASE(3);
+    shift(); kohonen::klistsamples(args(), false, false);
+    }
+  else if(argis("-somlistbest")) {
+    PHASE(3);
+    shift(); kohonen::klistsamples(args(), true, false);
+    }
+  else if(argis("-somlistbestc")) {
+    PHASE(3);
+    shift(); kohonen::klistsamples(args(), true, true);
+    }
+  else if(argis("-somndist")) {
+    PHASE(3);
+    shift(); kohonen::neurondisttable(args());
+    }
+  else if(argis("-somshowbest")) {
+    showbestsamples();
+    }
+
+  else return 1;
+  return 0;
+  }
+
+auto hooks = addHook(hooks_args, 100, readArgs) + addHook(hooks_frame, 50, levelline::draw);
+}
+
+void mark(cell *c) {
+  using namespace kohonen;
+  distfrom = getNeuronSlow(c);
+  coloring();
+  }
+