package umontreal.iro.lecuyer.functionfit;

import java.io.Serializable;
import umontreal.iro.lecuyer.functions.Polynomial;

import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.impl.DenseDoubleMatrix2D;
import cern.colt.matrix.linalg.Algebra;


/**
 * Represents a polynomial obtained by the least squares method on a set of points.
 * More specifically, let 
 * (x0, y0),…,(xn, yn) be a set of points and
 * p(x) the constructed polynomial of degree m. The constructed polynomial
 * minimizes the square error 

 * 
 * 

 * E2 = ∑i=0n[yi - p(xi)]2.
 * 

 * 
 */
public class LeastSquares extends Polynomial implements Serializable {
   private static final long serialVersionUID = -4997132164503234983L;
   private static final Algebra alg = new Algebra ();
   private double[] x;
   private double[] y;



   /**
    * Constructs a new least squares polynomial with points (x[0],
    *  y[0]),..., (x[n], y[n]). The constructed polynomial has degree
    *  degree.
    * 
    * @param x the x coordinates of the points.
    * 
    *    @param y the y coordinates of the points.
    * 
    *    @param degree the degree of the polynomial.
    * 
    *    @exception NullPointerException if x or y are null.
    * 
    *    @exception IllegalArgumentException if the lengths of x and y are different,
    *                or if less than degree + 1 points are specified.
    * 
    */
   public LeastSquares (double[] x, double[] y, int degree) {
      super (getCoefficients (x, y, degree));
      this.x = x.clone ();
      this.y = y.clone ();
   }


   /**
    * Computes and returns the coefficients of the fitting polynomial of
    *   degree degree.
    *  The coordinates of the given points are  (x[i], y[i]).
    * 
    * @param x the x coordinates of the points.
    * 
    *    @param y the y coordinates of the points.
    * 
    *    @param degree the degree of the polynomial.
    * 
    *    @return the coefficients of the fitting polynomial.
    * 
    */
   public static double[] getCoefficients (double[] x, double[] y,
                                           int degree) {
      if (x.length != y.length)
         throw new IllegalArgumentException ("Length of x and y not equal");
      if (x.length < degree + 1)
         throw new IllegalArgumentException ("Not enough points");

      final double[] xSums = new double[2 * degree + 1];
      final double[] xySums = new double[degree + 1];
      xSums[0] = x.length;
      for (int i = 0; i < x.length; i++) {
         double xv = x[i];
         xySums[0] += y[i];
         for (int j = 1; j <= 2 * degree; j++) {
            xSums[j] += xv;
            if (j <= degree)
               xySums[j] += xv * y[i];
            xv *= x[i];
         }
      }
      final DoubleMatrix2D A = new DenseDoubleMatrix2D (degree + 1, degree + 1);
      final DoubleMatrix2D B = new DenseDoubleMatrix2D (degree + 1, 1);
      for (int i = 0; i <= degree; i++) {
         for (int j = 0; j <= degree; j++) {
            final int d = i + j;
            A.setQuick (i, j, xSums[d]);
         }
         B.setQuick (i, 0, xySums[i]);
      }
      final DoubleMatrix2D aVec = alg.solve (A, B);
      return aVec.viewColumn (0).toArray ();
   }


   /**
    * Returns the x coordinates of the fitted points.
    * 
    * @return the x coordinates of the fitted points.
    * 
    */
   public double[] getX() {
      return x;
   }


   /**
    * Returns the y coordinates of the fitted points.
    * 
    * @return the y coordinates of the fitted points.
    * 
    */
   public double[] getY() {
      return y;
   }


   /**
    * Calls {@link umontreal.iro.lecuyer.functionfit.PolInterp#toString(double[], double[]) toString}
    *  with the associated points.
    * 
    * @return a string containing the points.
    * 
    */
   public String toString() {
      return PolInterp.toString (x, y);
   }

   @Override



   public LeastSquares clone() {
      final LeastSquares ls = (LeastSquares) super.clone ();
      ls.x = x.clone ();
      ls.y = y.clone ();
      return ls;
   }
}