package umontreal.iro.lecuyer.functionfit;

import umontreal.iro.lecuyer.util.Misc;
import umontreal.iro.lecuyer.util.RootFinder;
import umontreal.iro.lecuyer.functions.MathFunction;
import umontreal.iro.lecuyer.functions.MathFunctionWithIntegral;
import umontreal.iro.lecuyer.functions.MathFunctionWithDerivative;
import umontreal.iro.lecuyer.functions.MathFunctionWithFirstDerivative;
import umontreal.iro.lecuyer.functions.MathFunctionUtil;
import cern.colt.matrix.linalg.Algebra;
import cern.colt.matrix.DoubleMatrix2D;
import cern.colt.matrix.DoubleFactory2D;
import java.util.Arrays;
import java.util.Random;
import java.io.*;



/**
 * Represents a B-spline with control points at
 * 
 * (Xi, Yi).
 * Let 
 * Pi = (Xi, Yi), for 
 * i = 0,…, n - 1, be a control point and
 * let tj, for 
 * j = 0,…, m - 1 be a knot.
 * A B-spline of degree p = m - n - 1 is a parametric curve defined as
 * 

 * 
 * 

 * P(t) = ∑i=0n-1Ni, p(t)Pi, for tp <= t <= tm-p-1.
 * 

 * Here,
 * 

 * 

 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
 * 
Ni, p(t)
=
((t - ti)/(ti+p - ti))Ni, p-1(t) + ((ti+p+1 - t)/(ti+p+1 - ti+1))Ni+1, p-1(t)

 *  
Ni, 0(t)
=
{

 *  
 * 

 * 
 * 

 * This class provides methods to evaluate
 * 
 * P(t) = (X(t), Y(t)) at any value of t,
 * for a B-spline of any degree p >= 1.
 * Note that the
 * {@link #evaluate(double) evaluate} method
 * of this class can be slow, since it
 * uses a root finder to determine
 * the value of t* for which X(t*) = x
 * before it computes Y(t*).
 * 
 */
public class BSpline implements MathFunction

,
MathFunctionWithIntegral, MathFunctionWithDerivative, MathFunctionWithFirstDerivative{
   // Formula taken from http://www.ibiblio.org/e-notes/Splines/Basis.htm
   // and http://en.wikipedia.org/wiki/De_Boor_algorithm
   private double[] x;     //x original
   private double[] y;     //y original

   private int degree;

   //variables sur lesquelles on travaille
   private double[] myX;
   private double[] myY;
   private double[] knots;



   /**
    * Constructs a new uniform B-spline of degree degree
    *    with control points at (x[i], y[i]).
    *    The knots of the resulting B-spline are set uniformly from
    *     x[0] to x[n-1].
    * 
    * @param x the values of X.
    * 
    *    @param y the values of Y.
    * 
    *    @param degree the degree of the B-spline.
    * 
    * 
    */
   public BSpline (final double[] x, final double[] y, final int degree) {
      if (x.length != y.length)
         throw new IllegalArgumentException("The arrays x and y must share the same length");
      this.degree = degree;
      this.x = x.clone();
      this.y = y.clone();
      init(x, y, null);
   }


   /**
    * Constructs a new uniform B-spline
    *    with control points at (x[i], y[i]), and
    *    knot vector given by the array knots.
    * 
    * @param x the values of X.
    * 
    *    @param y the values of Y.
    * 
    *    @param knots the knots of the B-spline.
    * 
    */
   public BSpline (final double[] x, final double[] y, final double[] knots) {
      if (x.length != y.length)
         throw new IllegalArgumentException("The arrays x and y must share the same length");
      if (knots.length <= x.length + 1)
         throw new IllegalArgumentException("The number of knots must be at least n+2");

      this.x = x.clone();
      this.y = y.clone();
      this.knots = knots.clone();
      init(x, y, knots);
   }


   /**
    * Returns the Xi coordinates for this spline.
    * 
    * @return the Xi coordinates.
    * 
    */
   public double[] getX() {
      return myX.clone ();
   }


   /**
    * Returns the Yi coordinates for this spline.
    * 
    * @return the Yi coordinates.
    * 
    */
   public double[] getY() {
      return myY.clone ();
   }


   /**
    * Returns the knot maximal value.
    * 
    * @return the Yi coordinates.
    * 
    */
   public double getMaxKnot() {
      return knots[knots.length-1];
   }


   /**
    * Returns the knot minimal value.
    * 
    * @return the Yi coordinates.
    * 
    */
   public double getMinKnot() {
      return knots[0];
   }


   /**
    * Returns an array containing the knot vector 
    * (t0, tm-1).
    * 
    * @return the knot vector.
    * 
    */
   public double[] getKnots() {
      return knots.clone ();
   }


   /**
    * Returns a B-spline curve of degree degree interpolating the
    *    
    * (xi, yi) points.
    *    This method uses the uniformly spaced method for interpolating
    *    points with a B-spline curve, and a uniformed clamped knot vector,
    *    as described in http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/.
    *    
    * @param x the values of X.
    * 
    *    @param y the values of Y.
    * 
    *    @param degree the degree of the B-spline.
    * 
    *    @return the B-spline curve.
    * 
    */
   public static BSpline createInterpBSpline (double[] x, double[] y,
                                              int degree)  {
      if (x.length != y.length)
         throw new IllegalArgumentException("The arrays x and y must share the same length");
      if (x.length <= degree)
         throw new IllegalArgumentException("The arrays length must be greater than degree");

      int n = x.length-1;
      //compute t : parameters vector uniformly from 0 to 1
      double[] t = new double[x.length];
      for(int i =0; idegree smoothing
    *    
    * (xi, yi), for 
    * i = 0,…, n points. The
    *    precision depends on the parameter h: 
    * 1 <= degree <= h < n, which
    *    represents the number of control points used by the new B-spline curve,
    *    minimizing the quadratic error
    *    

    * 
    * 

    * L = ∑i=0n()2.
    * 

    * This method uses the uniformly spaced method for interpolating
    *    points with a B-spline curve and a uniformed clamped knot vector,
    *    as described in http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/.
    * 
    * @param x the values of X.
    * 
    *    @param y the values of Y.
    * 
    *    @param degree the degree of the B-spline.
    * 
    *    @param h the desired number of control points.
    * 
    *    @return the B-spline curve.
    * 
    */
   public static BSpline createApproxBSpline (double[] x, double[] y,
                                              int degree, int h)  {
      if (x.length != y.length)
         throw new IllegalArgumentException("The arrays x and y must share the same length");
      if (x.length <= degree)
         throw new IllegalArgumentException("The arrays length must be greater than degree");

      //compute t : parameters vector uniformly from 0 to 1
      double[] t = new double[x.length];
      for(int i =0; iderivative method,
    *    is strongly recommended if one wants to compute many derivative values.
    * 
    * @return the derivative B-spline of the current variable.
    * 
    */
   public BSpline derivativeBSpline()  {
      double xTemp[] = new double[this.myX.length-1];
      double yTemp[] = new double[this.myY.length-1];
      for(int i = 0; i xTemp[j])
               k++;
         }
         while(xTemp2[k] != 0)
            k++;
         xTemp2[k] = xTemp[i];
         yTemp2[k] = yTemp[i];
      }

      return new BSpline(xTemp2, yTemp2, newKnots);
   }


   /**
    * Returns the ith derivative B-spline object of the current variable;
    *    i must be less than the degree of the original B-spline.
    *    Using this function and the returned object, instead of the
    *   derivative method,
    *    is strongly recommended if one wants to compute many derivative values.
    * 
    * @param i the degree of the derivative.
    * 
    *    @return the ith derivative.
    * 
    */
   public BSpline derivativeBSpline (int i)  {
      BSpline bs = this;
      while(i > 0) {
         i--;
         bs = bs.derivativeBSpline();
      }
      return bs;
   }

   public double evaluate(final double u) {
      final MathFunction xFunction = new MathFunction () {
         public double evaluate (double t) {
            return evalX(t) - u;
         }
      };
      final double t = RootFinder.brentDekker (0, 1-1.0E-6, xFunction, 1e-6);
      return evalY(t);
   }

   public double integral (double a, double b) {
      return MathFunctionUtil.simpsonIntegral (this, a, b, 500);
   }

   public double derivative(double u) {
      return derivativeBSpline().evaluate(u);
   }

   public double derivative(double u, int n) {
      return derivativeBSpline(n).evaluate(u);
   }

   private void init(double[] x, double[] y, double [] initialKnots) {
      if(initialKnots == null) {
      //Cree un vecteur de noeud uniforme entre 0 et 1
         knots = new double[x.length+degree+1];
         for(int i = degree; i < this.knots.length-degree; i++)
            this.knots[i]= (double)(i-degree)/(knots.length - (2.0*degree) -1);
         for(int i = this.knots.length-degree; i < this.knots.length; i++)
            this.knots[i]=this.knots[i-1];
         for(int i = degree; i > 0; i--)
            this.knots[i-1]=this.knots[i];

      // cree notre vecteur interne de Points de controle
      // ici, aucune modification a faire sur les tableaux originaux
         myX = x;
         myY = y;
      }
      else {
         degree = initialKnots.length - x.length -1;

      // on adapte le tableau des noeuds a notre algorithme
      // le tableau knot necessite d'avoir degree fois la meme valeur en debut et en fin de tableau
      // on adapte la taille des tableau X et Y en consequence afin de continuer a respecter la condition :
      // x.length + degree + 1 = this.knots.length
      // Cette modification n'influence pas le resultat et permet de fermer notre courbe

         //Compute the number of values wich need to be added
         int iBorneInf = 1, iBorneSup = initialKnots.length-2;
         while(initialKnots[iBorneInf] == initialKnots[0])
            iBorneInf++;
         if(iBorneInf <= degree)
            iBorneInf = degree-iBorneInf+1;
         else
            iBorneInf=0;//on a alors iBorneInf valeurs a rajouter en debut de tableau

         while(initialKnots[iBorneSup] == initialKnots[initialKnots.length-1])
            iBorneSup--;
         if(iBorneSup >= initialKnots.length-1-degree)
            iBorneSup = degree+1-(initialKnots.length-1-iBorneSup);
         else
            iBorneSup = 0; //on a alors iBorneSup valeurs a rajouter en fin de tableau

         //add computed values
         this.knots = new double[initialKnots.length + iBorneInf + iBorneSup];
         myX = new double[x.length + iBorneInf + iBorneSup];
         myY = new double[y.length + iBorneInf + iBorneSup];
         for(int i = 0; i