//Kohonen Self-Organising-Map
//adapted from
//<http://www.patol.com/java/fill/index.html>
class SOM{
  int W,H,NGEONEURON,j;
  double  px,py,pz,theta,phi,momentum,mindist;
  geoNeuron [] gn;
  double [][] r;
  double [] xv,yv,zv;
  boolean stable;
  SOM(int size){
    theta = 0.5;
    phi   = 0.5;
    momentum = 0.999;
    W = size;
    H = W;
    NGEONEURON = W * H;
    gn = new geoNeuron[NGEONEURON];
    for(int x = 0; x < W; x++)
      for(int y = 0; y < H; y++){
        gn[x + y * W] = new geoNeuron((double)x/(double)(W-1), (double)y/(double)(H-1), 0.0);
      }
    r = new double[NGEONEURON][NGEONEURON];
    makeR(theta);
  }
  void run(){
    px = Math.random();
    py = Math.random();
    pz = Math.random();
    // SEARCH FOR MAXIMAL
    mindist = 0.0;
    j = -1;
    for(int i = 0; i < NGEONEURON;i++){
      double d = (px - gn[i].wx) * (px - gn[i].wx) + (py - gn[i].wy) * (py - gn[i].wy) + (pz - gn[i].wz) * (pz - gn[i].wz);
      if(d > mindist){
        mindist = d;
        j = i;
      }
    }
    // UPDATE WEIGHTS
    for(int i=0; i < NGEONEURON;i++){
      gn[i].wx += (phi * r[i][j] * (px - gn[i].wx));
      gn[i].wy += (phi * r[i][j] * (py - gn[i].wy));
      gn[i].wz += (phi * r[i][j] * (pz - gn[i].wz));
    }
    // DECREASE LEARNING PARAMETERS
    phi *= momentum;
    theta *= momentum;
    if (theta < 0.1) stable = true;
    // RE-COMPUTE r MATRIX
    makeR(theta);
  }
  void makeR(double th){
    for(int i = 0; i < NGEONEURON; i++){
      r[i][i] = 1.0;
      for(int j = i + 1; j < NGEONEURON; j++){
        r[i][j] = Math.exp( -1.0 * ( gn[i].dist(gn[j]) * gn[i].dist(gn[j])) / (2.0 * th * th));
        r[j][i] = r[i][j];
      }
    }
  }
  void addVector(double x, double y, double z){
    double [] xtemp = new double[xv.length+1];
    double [] ytemp = new double[yv.length+1];
    double [] ztemp = new double[zv.length+1];
    System.arraycopy(xv, 0, xtemp, 0, xv.length);
    System.arraycopy(yv, 0, ytemp, 0, yv.length);
    System.arraycopy(zv, 0, ztemp, 0, zv.length);
    xtemp[xv.length] = x;
    ytemp[yv.length] = y;
    ztemp[zv.length] = z;
    xv = xtemp;
    yv = ytemp;
    zv = ztemp;
  }
  void draw(){
    /*
    stroke(0);
     //DRAW PATH
     fill(255);
     stroke(255,0,0);
     for(int i = 0; i < gn.length; i++){
     pushMatrix();
     translate(-w+(float)(gn[i].wx)*width, -h+(float)(gn[i].wy)*height, (float)(gn[i].wz)*width);
     ellipse(0,0,4,4);
     popMatrix();
     }
     for(int i = 0; i < H; i++)
     for(int j = 1; j < W; j++){
     line(-w+(float)(gn[(j-1)+(i*H)].wx)*width, -h+(float)(gn[(j-1)+(i*H)].wy)*height, (float)(gn[(j-1)+(i*H)].wz)*height, -w+(float)(gn[(j)+(i*H)].wx)*width, -h+(float)(gn[(j)+(i*H)].wy)*height, (float)(gn[(j)+(i*H)].wz)*height);
     line(-w+(float)(gn[(i)+((j-1)*H)].wx)*width, -h+(float)(gn[(i)+((j-1)*H)].wy)*height, (float)(gn[(i)+((j-1)*H)].wz)*height, -w+(float)(gn[(i)+(j*H)].wx)*width, -h+(float)(gn[(i)+(j*H)].wy)*height, (float)(gn[(i)+(j*H)].wz)*height);
     }
     */
    noStroke();
    fill(255,0,0);
    //float size = buffer.width / W;
    for(int i = 1; i < H; i++)
      for(int j = 1; j < W; j++){
        beginShape(QUADS);
        //texture(buffer);
        vertex(-w+(float)(gn[(j-1)+((i-1)*H)].wx)*width, -h+(float)(gn[(j-1)+((i-1)*H)].wy)*height, (float)(gn[(j-1)+((i-1)*H)].wz)*height);//, (i-1)*size, (j-1)*size);
        vertex(-w+(float)(gn[(j)+((i-1)*H)].wx)*width, -h+(float)(gn[(j)+((i-1)*H)].wy)*height, (float)(gn[(j)+((i-1)*H)].wz)*height);//, (i)*size, (j-1)*size);
        vertex(-w+(float)(gn[(j)+(i*H)].wx)*width, -h+(float)(gn[(j)+(i*H)].wy)*height, (float)(gn[(j)+((i)*H)].wz)*height);//, (i)*size, (j)*size);
        vertex(-w+(float)(gn[(j-1)+(i*H)].wx)*width, -h+(float)(gn[(j-1)+(i*H)].wy)*height, (float)(gn[(j-1)+((i)*H)].wz)*height);//, (i-1)*size, (j)*size);
        endShape();
      }
    /*
    stroke(255,0,255);
     noFill();
     for (int i = 0; i < xv.length; i++){
     pushMatrix();
     translate(-w+(float)xv[i]*width, -h+(float)yv[i]*height, (float)zv[i]*height);
     ellipse(0,0,15,15);
     popMatrix();
     }
     fill(255);
     pushMatrix();
     translate(-w+(float)(px)*width, -h+(float)(py)*height, (float)(pz)*height);
     ellipse(0,0,10,10);
     popMatrix();
     */
  }
  Vec3 [] netImage(){
    Vec3 [] vectorImage = new Vec3[gn.length];
    for(int i = 0; i < gn.length; i++)
      vectorImage[i] = new Vec3((float)gn[i].wx,(float)gn[i].wy,(float)gn[i].wz);
    return vectorImage;
  }
}
void drawNetImage(int k){
  int H = netSize;
  int W = H;
  //float size = buffer.width / W;
  float step = (1.0 / frames) * k;
  for(int i = 1; i < H; i++)
    for(int j = 1; j < W; j++){
      beginShape(QUADS);
      //texture(buffer);
      vertex(-w+(animTemp[0][(j-1)+((i-1)*H)].lerp(animTemp[1][(j-1)+((i-1)*H)], step).x)*width, -h+(animTemp[0][(j-1)+((i-1)*H)].lerp(animTemp[1][(j-1)+((i-1)*H)], step).y)*height, (animTemp[0][(j-1)+((i-1)*H)].lerp(animTemp[1][(j-1)+((i-1)*H)], step).z)*height);//, (i-1)*size, (j-1)*size);
      vertex(-w+(animTemp[0][(j)+((i-1)*H)].lerp(animTemp[1][(j)+((i-1)*H)], step).x)*width, -h+(animTemp[0][(j)+((i-1)*H)].lerp(animTemp[1][(j)+((i-1)*H)], step).y)*height, (animTemp[0][(j)+((i-1)*H)].lerp(animTemp[1][(j)+((i-1)*H)], step).z)*height);//, (i)*size, (j-1)*size);
      vertex(-w+(animTemp[0][(j)+(i*H)].lerp(animTemp[1][(j)+(i*H)], step).x)*width, -h+(animTemp[0][(j)+(i*H)].lerp(animTemp[1][(j)+(i*H)], step).y)*height, (animTemp[0][(j)+((i)*H)].lerp(animTemp[1][(j)+((i)*H)], step).z)*height);//, (i)*size, (j)*size);
      vertex(-w+(animTemp[0][(j-1)+(i*H)].lerp(animTemp[1][(j-1)+(i*H)], step).x)*width, -h+(animTemp[0][(j-1)+(i*H)].lerp(animTemp[1][(j-1)+(i*H)], step).y)*height, (animTemp[0][(j-1)+((i)*H)].lerp(animTemp[1][(j-1)+((i)*H)], step).z)*height);//, (i-1)*size, (j)*size);
      endShape();
    }
}