SmoothingCubicSpline.java

/*
 * Class:        SmoothingCubicSpline
 * Description:  smoothing cubic spline algorithm of Schoenberg
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since
 *
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */
package umontreal.ssj.functionfit;
 
import umontreal.ssj.functions.*;
import umontreal.ssj.functions.Polynomial;
import umontreal.ssj.functions.MathFunctionWithFirstDerivative;
import umontreal.ssj.functions.MathFunctionWithDerivative;
import umontreal.ssj.functions.MathFunctionWithIntegral;
import umontreal.ssj.functions.SquareMathFunction;
import umontreal.ssj.functions.MathFunctionUtil;

public class SmoothingCubicSpline
      implements MathFunction, MathFunctionWithFirstDerivative, MathFunctionWithDerivative, MathFunctionWithIntegral {
 
   private Polynomial[] splineVector;
   private double[] x, y, weight;
   private double rho;

   public SmoothingCubicSpline(double[] x, double[] y, double[] w, double rho) {
      if (x.length != y.length)
         throw new IllegalArgumentException("x.length != y.length");
      if (w != null && x.length != w.length)
         throw new IllegalArgumentException("x.length != w.length");
      if (rho < 0 || rho > 1)
         throw new IllegalArgumentException("rho not in [0, 1]");
 
      splineVector = new Polynomial[x.length + 1];
      this.rho = rho;
      this.x = x.clone();
      this.y = y.clone();
 
      weight = new double[x.length];
      if (w == null) {
         for (int i = 0; i < weight.length; i++)
            weight[i] = 1.0;
      } else {
         for (int i = 0; i < weight.length; i++)
            weight[i] = w[i];
      }
 
      resolve();
   }

   public SmoothingCubicSpline(double[] x, double[] y, double rho) {
      this(x, y, null, rho);
   }

   public double evaluate(double z) {
      int i = getFitPolynomialIndex(z);
      if (i == 0)
         return splineVector[i].evaluate(z - x[0]);
      else
         return splineVector[i].evaluate(z - x[i - 1]);
   }

   public double integral(double a, double b) {
      int iA = getFitPolynomialIndex(a);
      int iB = getFitPolynomialIndex(b);
      double retour;
      int i = 1;
 
      // Calcule l'integrale
      if (iA == iB) { // les deux valeurs sont sur le meme polynome
         retour = splineVector[iB].integral(a - x[iB], b - x[iB]);
      } else {
         if (iA == 0)
            retour = splineVector[iA].integral(a - x[iA], 0);
         else
            retour = splineVector[iA].integral(a - x[iA], x[iA + 1] - x[iA]);
         for (i = iA + 1; i < iB; i++) {
            retour += splineVector[i].integral(0, x[i + 1] - x[i]);
         }
         retour += splineVector[iB].integral(0, b - x[iB]);
      }
      return retour;
   }

   public double derivative(double z) {
      int i = getFitPolynomialIndex(z);
      if (i == 0)
         return splineVector[i].derivative(z - x[0]);
      else
         return splineVector[i].derivative(z - x[i - 1]);
   }

   public double derivative(double z, int n) {
      int i = getFitPolynomialIndex(z);
      if (i == 0)
         return splineVector[i].derivative(z - x[0], n);
      else
         return splineVector[i].derivative(z - x[i - 1], n);
   }

   public double[] getX() {
      return x.clone();
   }

   public double[] getY() {
      return y.clone();
   }

   public double[] getWeights() {
      return weight;
   }

   public double getRho() {
      return rho;
   }

   public Polynomial[] getSplinePolynomials() {
      return splineVector.clone();
   }

   public int getFitPolynomialIndex(double x) {
      // Algorithme de recherche binaire legerement modifie
      int j = this.x.length - 1;
      if (x > this.x[j])
         return j + 1;
      int tmp = 0;
      int i = 0;
 
      while (i + 1 != j) {
         if (x > this.x[tmp]) {
            i = tmp;
            tmp = i + (j - i) / 2;
         } else {
            j = tmp;
            tmp = i + (j - i) / 2;
         }
         if (j == 0) // le point est < a x_0, on sort
            i--;
      }
      return i + 1;
   }
 
   private void resolve() {
      /*
       * taken from D.S.G Pollock's paper, "Smoothing with Cubic Splines", Queen Mary,
       * University of London (1993) http://www.qmw.ac.uk/~ugte133/PAPERS/SPLINES.PDF
       */
 
      double[] h = new double[x.length];
      double[] r = new double[x.length];
      double[] u = new double[x.length];
      double[] v = new double[x.length];
      double[] w = new double[x.length];
      double[] q = new double[x.length + 1];
      double[] sigma = new double[weight.length];
 
      for (int i = 0; i < weight.length; i++) {
         if (weight[i] <= 0.0)
            sigma[i] = 1.0e100;
         else
            sigma[i] = 1.0 / Math.sqrt(weight[i]);
      }
 
      double mu;
      int n = x.length - 1;
 
      if (rho <= 0)
         mu = 1.0e100; // arbitrary large number to avoid 1/0
      else
         mu = 2 * (1 - rho) / (3 * rho);
 
      h[0] = x[1] - x[0];
      r[0] = 3 / h[0];
      for (int i = 1; i < n; i++) {
         h[i] = x[i + 1] - x[i];
         r[i] = 3 / h[i];
         q[i] = 3 * (y[i + 1] - y[i]) / h[i] - 3 * (y[i] - y[i - 1]) / h[i - 1];
      }
 
      for (int i = 1; i < n; i++) {
         u[i] = r[i - 1] * r[i - 1] * sigma[i - 1] + (r[i - 1] + r[i]) * (r[i - 1] + r[i]) * sigma[i]
               + r[i] * r[i] * sigma[i + 1];
         u[i] = mu * u[i] + 2 * (x[i + 1] - x[i - 1]);
         v[i] = -(r[i - 1] + r[i]) * r[i] * sigma[i] - r[i] * (r[i] + r[i + 1]) * sigma[i + 1];
         v[i] = mu * v[i] + h[i];
         w[i] = mu * r[i] * r[i + 1] * sigma[i + 1];
      }
      q = Quincunx(u, v, w, q);
 
      // extrapolation a gauche
      double[] params = new double[4];
      double dd;
      params[0] = y[0] - mu * r[0] * q[1] * sigma[0];
      dd = y[1] - mu * ((-r[0] - r[1]) * q[1] + r[1] * q[2]) * sigma[1];
      params[1] = (dd - params[0]) / h[0] - q[1] * h[0] / 3;
      splineVector[0] = new Polynomial(params);
 
      // premier polynome
      params[0] = y[0] - mu * r[0] * q[1] * sigma[0];
      dd = y[1] - mu * ((-r[0] - r[1]) * q[1] + r[1] * q[2]) * sigma[1];
      params[3] = q[1] / (3 * h[0]);
      params[2] = 0;
      params[1] = (dd - params[0]) / h[0] - q[1] * h[0] / 3;
      splineVector[1] = new Polynomial(params);
 
      // les polynomes suivants
      int j;
      for (j = 1; j < n; j++) {
         params[3] = (q[j + 1] - q[j]) / (3 * h[j]);
         params[2] = q[j];
         params[1] = (q[j] + q[j - 1]) * h[j - 1] + splineVector[j].getCoefficient(1);
         params[0] = r[j - 1] * q[j - 1] + (-r[j - 1] - r[j]) * q[j] + r[j] * q[j + 1];
         params[0] = y[j] - mu * params[0] * sigma[j];
         splineVector[j + 1] = new Polynomial(params);
      }
 
      // extrapolation a droite
      j = n;
      params[3] = 0;
      params[2] = 0;
      params[1] = splineVector[j].derivative(x[x.length - 1] - x[x.length - 2]);
      params[0] = splineVector[j].evaluate(x[x.length - 1] - x[x.length - 2]);
      splineVector[n + 1] = new Polynomial(params);
   }
 
   private static double[] Quincunx(double[] u, double[] v, double[] w, double[] q) {
      u[0] = 0;
      v[1] = v[1] / u[1];
      w[1] = w[1] / u[1];
      int j;
 
      for (j = 2; j < u.length - 1; j++) {
         u[j] = u[j] - u[j - 2] * w[j - 2] * w[j - 2] - u[j - 1] * v[j - 1] * v[j - 1];
         v[j] = (v[j] - u[j - 1] * v[j - 1] * w[j - 1]) / u[j];
         w[j] = w[j] / u[j];
      }
 
      // forward substitution
      q[1] = q[1] - v[0] * q[0];
      for (j = 2; j < u.length - 1; j++) {
         q[j] = q[j] - v[j - 1] * q[j - 1] - w[j - 2] * q[j - 2];
      }
      for (j = 1; j < u.length - 1; j++) {
         q[j] = q[j] / u[j];
      }
 
      // back substitution
      q[u.length - 1] = 0;
      for (j = u.length - 3; j > 0; j--) {
         q[j] = q[j] - v[j] * q[j + 1] - w[j] * q[j + 2];
      }
      return q;
   }
}