import java.util.*;
import FloatPoint;
import java.awt.*;
import java.lang.*;
import FourByFourSystem;

public class FunctionCurve extends Object implements Cloneable {

  Vector points, dpoints;
  double moments[][];
  int depth, p, u;
  boolean reg;
  FloatPoint origin;

  static final int VECTOR_SIZE = 33;
  static final int DRAWING_SCALER = 5;
  static final int DRAWING_TOLERANCE = 5;

  static final double[][] pca = 
  { {0}, 
    {.5,.5},
    {-.0625,.5625,.5625,-.0625},
    {.0117188,-.0976562,.585938,.585938,-.0976562,.0117188} };

  static final double[][] uca = 
  { {0}, 
    {.25,.25}, 
    {-.03125,.28125,.28125,-.03125} };

/*************************/
  public FunctionCurve() {
    points = new Vector(VECTOR_SIZE);
    dpoints = new Vector(VECTOR_SIZE*DRAWING_SCALER);
    depth = 1;
    origin = new FloatPoint(0,0);
  }

/*************************/
  public void setOrigin(int x, int y) {
    origin.x = x;
    origin.y = y;
  }

/*************************/
  public FloatPoint screenToGraph(FloatPoint p) {
    return new FloatPoint(p.x - origin.x, p.y - origin.y);
  }

/*************************/
  public FloatPoint graphToScreen(FloatPoint p) {
    return new FloatPoint(p.x + origin.x, p.y + origin.y);
  }

/*************************/
  public void addPoint(FloatPoint p) {
    p = screenToGraph(p);
    for(int i=0; i<points.size(); i++) {
      if(((FloatPoint)pAt(i)).x > p.x) {
	points.insertElementAt(p, i);
        return;
      }
      if(((FloatPoint)pAt(i)).x == p.x) {
	points.setElementAt(p, i);
	return;
      }
    }
    points.addElement(p);
  }

/************************/
  public FloatPoint pAt(int i) { return (FloatPoint)points.elementAt(i); }
  public FloatPoint dAt(int i) { return (FloatPoint)dpoints.elementAt(i); }

/************************/
  void calcMoments() {
    double a,b;
    moments = new double[u][points.size()];
    for(int x=0; x<u; x++)
      for(int y=0; y<(points.size()); y++) {
	if(y==0) {
	  b = pAt(1).x - pAt(0).x;
	  a = 0;
	}
	else
	  if(y==points.size()-1) {
	    a = pAt(points.size()-1).x - pAt(points.size()-2).x;
	    b = 0;
	  }
	  else {
	    a = pAt(y).x - pAt(y-1).x;
	    b = pAt(y+1).x - pAt(y).x;
	  }
	switch(x) {
	case 0 : moments[x][y] = a/2 + b/2;  break;
	case 1 : moments[x][y] = b*b/6 - a*a/6;  break;
	case 2 : moments[x][y] = a*a*a/12 + b*b*b/12;  break;
	case 3 : moments[x][y] = b*b*b*b/20 - a*a*a*a/20;  break;
	}
      }
  }

/***************************/
  public boolean canDrop() {
    if( ((points.size()-1)%(depth*2) != 0) || ((points.size()/depth+1) < (2*p-1)) )
      return false;
    return true;
  }
  
/***************************/
  public void drop() {    
    if(depth==1 && !reg)
      calcMoments();

    FloatPoint ua[] = new FloatPoint[u], pa[] = new FloatPoint[p];
    
    for(int i=0; i<p; i++)
      pa[i] = pAt(2*depth*i);
    
    if(!reg)
      for(int i=0; i<u; i++)
	ua[i] = pAt(2*depth*i);

    FloatPoint d;
    for(int i=depth; i<(points.size()-1); i += 2*depth) {
      d = pAt(i);
      if(reg)
	d.y -= regPredict(pa, pca[p/2]);
      else {
	d.y -= neville(pa, d.x);
	FloatPoint uc[] = new FloatPoint[u];
	for(int order=0; order<(u/2+1); order++)
	  for(int a=0; a<u; a++) {
	    for(int b=0; b<u; b++) 
	      if(b==a)
		uc[b] = new FloatPoint(ua[b].x, 1.0);
	      else
		uc[b] = new FloatPoint(ua[b].x, 0.0);
	    moments[order][points.indexOf(ua[a])] += moments[order][i]*neville(uc, d.x);
	  }
      }
      
      if((i>=(p*depth-depth)) && i<((points.size()-1)-(p*depth-depth))) {
	for(int a=0; a<(p-1); a++)	
	  pa[a] = pa[a+1];	
	pa[p-1] = pAt(points.indexOf(pa[p-2]) + 2*depth);
      }
      
      if(!reg && (i>=(u*depth-depth)) && i<((points.size()-1)-(u*depth-depth)) && (u!=0)) {
	for(int a=0; a<(u-1); a++)	
	  ua[a] = ua[a+1];	
	ua[u-1] = pAt(points.indexOf(ua[u-2]) + 2*depth);
      }
    }  

    for(int i=0; i<u; i++)
      ua[i] = pAt(2*depth*i);  
    double uc[] = new double[4];
    
    for(int i=depth; i<(points.size()-1); i += 2*depth) {
      d = pAt(i);
      
      if(reg)
	for(int a=0; a<u; a++)
	  ua[a].y += d.y * uca[u/2][a];
      else {
	uc = calcUpdateCoeffs(d, ua);
	for(int z=0; z<u; z++)
	  ua[z].y += d.y*uc[z];
      }
      
      if((i>=(u*depth-depth)) && i<((points.size()-1)-(u*depth-depth)) && (u!=0)) {
	for(int a=0; a<(u-1); a++)	
	  ua[a] = ua[a+1];	
	ua[u-1] = pAt(points.indexOf(ua[u-2]) + 2*depth);
      }
    }
    
    depth *= 2;
  }
  
/*****************************/
  public void lift() {
    if(depth == 1)
      return;
    depth /= 2;
   
    FloatPoint[] pa = new FloatPoint[p], ua = new FloatPoint[u];
    
    for(int i=0; i<u; i++)
      ua[i] = pAt(2*depth*i);  

    FloatPoint d;
    for(int i=depth; i<(points.size()-1); i += 2*depth) {
      d = pAt(i);
      double uc[] = new double[u];

      if(reg)
	for(int a=0; a<u; a++)
	  ua[a].y -= d.y * uca[u/2][a];
      else {
	uc = calcUpdateCoeffs(d, ua);
	for(int z=0; z<u; z++)
	  ua[z].y -= d.y*uc[z];
      }
      
      if((i>=(u*depth-depth)) && i<((points.size()-1)-(u*depth-depth)) && (u!=0)) {
	for(int a=0; a<(u-1); a++)	
	  ua[a] = ua[a+1];	
	ua[u-1] = pAt(points.indexOf(ua[u-2]) + 2*depth);
      }
    }

    for(int i=0; i<p; i++)
      pa[i] = pAt(2*depth*i);
    
    if(!reg)
      for(int i=0; i<u; i++)
	ua[i] = pAt(2*depth*i);
    
    for(int i=depth; i<points.size(); i += 2*depth) {
      d = pAt(i);
      if(reg)
	d.y += regPredict(pa, pca[p/2]);
      else {
	FloatPoint uc[] = new FloatPoint[u];
	for(int order=0; order<(u/2+1); order++)
	  for(int a=0; a<u; a++) {
	    for(int b=0; b<u; b++) 
	      if(b==a)
		uc[b] = new FloatPoint(ua[b].x, 1.0);
	      else
		uc[b] = new FloatPoint(ua[b].x, 0.0);
	    moments[order][points.indexOf(ua[a])] -= moments[order][i]*neville(uc, d.x);
	  }
	d.y += neville(pa, d.x);
      }

      if((i>=(p*depth-depth)) && i<((points.size()-1)-(p*depth-depth))) {
	for(int a=0; a<(p-1); a++)
	  pa[a] = pa[a+1];
	pa[p-1] = pAt(points.indexOf(pa[p-2]) + 2*depth);
      }

      if(!reg && (i>=(u*depth-depth)) && i<((points.size()-1)-(u*depth-depth)) && (u!=0)) {
	for(int a=0; a<(u-1); a++)	
	  ua[a] = ua[a+1];	
	ua[u-1] = pAt(points.indexOf(ua[u-2]) + 2*depth);
      } 
    }
  }

/*********************************/
  public void drawCurve(Graphics g) {
    FloatPoint a=null,b=null,p1, p2;
    try
      a = (FloatPoint)dpoints.firstElement();
    catch (Exception e) return;
    p1 = graphToScreen(a);
    for(int i=1; i<dpoints.size(); i++) {
      b = dAt(i);
      p2 = graphToScreen(b);
      g.drawLine((int)p1.x,(int)p1.y,(int)p2.x,(int)p2.y);
      a = b;
      p1 = p2;
    }
  }

/********************************/
  public void smoothScalingCurve() {
    dpoints.removeAllElements();
    for(int i=0; i<points.size(); i += depth)
      dpoints.addElement(points.elementAt(i));
    smoothCurve();
  }

/********************************/
  public void smoothWaveletCurve() {
    dpoints.removeAllElements();
    if(depth==1)
      return;
    for(int i=depth/2; i<points.size(); i+=depth)
      dpoints.addElement(points.elementAt(i));
    smoothCurve();
  }


/********************************/
  public void smoothCurve() {
           
    if((dpoints.size()) <= 1 || (p==2))
      return;
    
    boolean done = false;
    FloatPoint a,b;
    int origsize, o;
    while(!done) {
      o = p;
      while(dpoints.size() < o)
	o -= 2;
      FloatPoint[] pa = new FloatPoint[o];
      origsize = dpoints.size();
      done = true;
      for(int i=0; i<o; i++)
	pa[i] = dAt(i);
      a = dAt(0);
      for(int i=1; i<origsize; i++) {
	b = (FloatPoint)dpoints.elementAt(dpoints.indexOf(a) + 1);
	if( (b.x - a.x) > DRAWING_TOLERANCE) {
	  double x = a.x + (b.x-a.x)/2;
	  dpoints.insertElementAt(new FloatPoint(x, neville(pa,x)), dpoints.indexOf(b));
	  done = false;
	}
	a = b;
	if( (i >= (o/2)) && (i < (origsize-o/2)) ) {
	  for(int z=0; z<(o-1); z++)
	    pa[z] = pa[z+1];
	  pa[o-1] = (FloatPoint)dpoints.elementAt(dpoints.indexOf(pa[o-2]) + 1);
	}
      }
    }   
  }   
  
/**********************/
  public void drawPoint(Graphics g, FloatPoint p) { 
    FloatPoint ps = graphToScreen(p);
    g.fillOval((int)ps.x-3,(int)ps.y-3,6,6); 
  }
  //need a const here for the point width;

/**********************/
  public void drawPoints(Graphics g) {
    for(int i=0; i<points.size(); i+= depth)
      drawPoint(g, pAt(i));
  }

/**********************/		
  public void clear() {
    points.removeAllElements();
    dpoints.removeAllElements();
    depth = 1;
  }

/**********************/		
  public FloatPoint findPoint(int x, int y) {
    FloatPoint p = screenToGraph(new FloatPoint((double)x,(double)y));
    for( int i=0; i < points.size(); i++ ) {
      FloatPoint n = pAt(i);
      if( (int)n.x - 3 < p.x && p.x < (int)n.x + 3 && (int)n.y - 3 < p.y && p.y < (int)n.y + 3 )
      //same const as in drawPoint here;
	return n;
    }
    return null;
  }

/**********************/		
  public void movePoint(FloatPoint p, int x, int y) {
    p.x = x - origin.x;
    p.y = y - origin.y;
  }

/**********************/		
  public void killPoint(FloatPoint p) { points.removeElement(p); }

/**********************/
  public double[] calcUpdateCoeffs(FloatPoint d, FloatPoint[] ua) {
    double uc[] = new double[u];
    switch(u) {
    case 0: break;
    case 2:
      int ai = points.indexOf(ua[0]);
      int bi = points.indexOf(ua[1]);
      double as = d.x - ua[0].x;
      double bs = ua[1].x - d.x;
      double ma0 = moments[0][ai];
      double ma1 = moments[1][ai] - as*moments[0][ai];
      double mb0 = moments[0][bi];
      double mb1 = bs*moments[0][bi] + moments[1][bi];
      double mw0 = moments[0][points.indexOf(d)];
      double mw1 = moments[1][points.indexOf(d)];
      uc[1] = (mw1*ma0-mw0*ma1)/(mb1*ma0-mb0*ma1); 
      uc[0] = (mw0-uc[1]*mb0)/ma0;
      break;
    case 4:
      int ind[] = new int[4];
      double sep[] = new double[4];
      double w[] = new double[4];
      for(int z=0; z<4; z++) {
	ind[z] = points.indexOf(ua[z]);
	sep[z] = ua[z].x - d.x;
	w[z] = moments[z][points.indexOf(d)];
      }
      double A[][] = new double[4][4];
      for(int r=0; r<4; r++)
	for(int c=0; c<4; c++)
	  switch(r) {
	  case 0:
	    A[r][c] = moments[0][ind[c]];
	    break;
	  case 1:
	    A[r][c] = moments[1][ind[c]] + sep[c]*moments[0][ind[c]];
	    break;
	  case 2:
	    A[r][c] = moments[2][ind[c]] + 2*sep[c]*moments[1][ind[c]];
	    A[r][c] += sep[c]*sep[c]*moments[0][ind[c]];
	    break;				   
	  case 3:
	    A[r][c] = moments[3][ind[c]] + 3*sep[c]*moments[2][ind[c]];
	    A[r][c] += 3*sep[c]*sep[c]*moments[1][ind[c]];
	    A[r][c] += sep[c]*sep[c]*sep[c]*moments[0][ind[c]];
	    break;
	  }
      FourByFourSystem.linSolve(A, uc, w);
      break;
    }
    return uc;
  }

/**********************/		
  double regPredict(FloatPoint[] pa, double[] ca) {
    double y=0;
    for(int i=0; i<pa.length; i++)
      y += pa[i].y * ca[i];
    return y;
  }

/**********************/		
  void regUpdate(FloatPoint[] ua, double[] ca, int d) {
    for(int i=0; i<ua.length; i++)
      ua[i].y += d * ca[i];
  }

/**********************/		
  boolean overlapping(FloatPoint p, int x) {
    x -= origin.x;

    if(p == points.firstElement()) {
      FloatPoint b = (FloatPoint)points.elementAt(points.indexOf(p) + 1);
      return (x >= (int)b.x);
    }
    FloatPoint a = (FloatPoint)points.elementAt(points.indexOf(p) - 1);
    if(p == points.lastElement())
      return (x <= (int)a.x);
    FloatPoint b = (FloatPoint)points.elementAt(points.indexOf(p) + 1);
    return ((x <= (int)a.x) || (x >= (int)b.x)); 
  }

/**********************/		
/*Given (p.length) Points this method uses the Neville algorithm to find 
  the value at xx of the polynomial of degree (pa.length-1) interpolating 
  the points pa.*/

  double neville(FloatPoint[] pa, double xx) {
    double[] vy = new double[pa.length];
    
    for (int i=0; i<pa.length; i++) {
      vy[i] = pa[i].y;
      for (int j=i-1; j>=0; j--)
	try
	  vy[j] = vy[j+1] + (vy[j+1] - vy[j]) * (xx - pa[i].x) / (pa[i].x - pa[j].x);
	catch(ArithmeticException e); //possible division by zero;
    }
    return vy[0];
  }
 
/******************************/
  public Object clone() {
    FunctionCurve n = new FunctionCurve();
    n.setOrigin((int)origin.x, (int)origin.y);
    for(int i=0; i<points.size(); i++)
      (n.points).addElement(pAt(i).clone());
    try {
      n.moments = new double[moments.length][];
      for(int a=0; a<moments.length; a++) {
	n.moments[a] = new double[moments[a].length];
	System.arraycopy(moments[a],0,n.moments[a],0,moments[a].length);
      }
    }
    catch(NullPointerException e);
    n.depth = depth;
    n.p = p;
    n.u = u;
    n.reg = reg;
    return n;
  }  
}