com.actelion.research.util.SmoothingSplineInterpolator Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of openchemlib Show documentation
Show all versions of openchemlib Show documentation
Open Source Chemistry Library
package com.actelion.research.util;
/**
* Smoothing Spline Interpolator based on the algorithm described at
* http://www.qmw.ac.uk/~ugte133/book/11_tsd/splines.pdf
*
*
* The Smoothing Spline is used to minimize Sum(sqr((Si-yi)/sigmai)) + lambda*Sum(sqr(S''i))
* If lambda=0 (default), this is equivalent to the cubic spline interpolation
* If lambda=Infinity, this is equivalent to the least square fitting
*
* FastSpline spline = new SmoothingSplineInterpolator().interpolate(X, Y);
*
* @author freyssj
*/
public class SmoothingSplineInterpolator {
private double lambda = 0;
private double[] sigma = null;
private double residuals;
private double smoothing;
/**
*
*/
private static void quincunx(int n, double[] u, double[] v, double[] w, double[] q) {
//Factorisation
u[0] = 0;
v[0] = 0;
w[0] = 0;
for (int j = 1; j <= n-1; j++) {
u[j] = u[j] - (j-2>=0? u[j-2]*sqr(w[j-2]): 0) - u[j-1]*sqr(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[0] = 0;
for (int j = 1; j <= n-1; j++) {
q[j] = q[j] - v[j-1] * q[j-1] - (j-2>=0? w[j-2]*q[j-2] :0);
}
for (int j = 1; j <= n-1; j++) {
q[j] = q[j] / u[j];
}
//Back Substitution
q[n-1] = q[n-1];
q[n-2] = q[n-2] - v[n-2] * q[n-2+1];
for (int j = n-3; j >=1; j--) {
q[j] = q[j] - v[j] * q[j+1] - w[j] * q[j+2];
}
}
/**
* @see org.apache.commons.math.analysis.UnivariateRealInterpolator#interpolate(double[], double[])
*/
public FastSpline interpolate(double[] x, double[] y) {
if (x.length != y.length) throw new IllegalArgumentException("Dataset arrays must have same length.");
if (x.length < 3) throw new IllegalArgumentException("At least 3 datapoints are required to compute a spline interpolant");
if (sigma!=null && sigma.length != x.length) throw new IllegalArgumentException("Sigma and dataset arrays must have same length");
for (int i = 0; i < x.length-1; i++) {
if(x[i+1]<=x[i]) throw new IllegalArgumentException("the X must be strictly increasing");
}
int n = x.length;
double[] h = new double[n-1];
double[] r = new double[n];
double[] f = new double[n];
double[] p = new double[n];
double[] q = new double[n];
double[] u = new double[n];
double[] v = new double[n];
double[] w = new double[n];
for (int i = 0; i < n-1; i++) {
h[i] = x[i + 1] - x[i];
r[i] = 3d / h[i];
}
for (int i = 1; i < n-1; i++) {
f[i] = -(r[i-1] + r[i]);
p[i] = 2d * (x[i + 1] - x[i - 1]);
q[i] = 3d * (y[i + 1] - y[i]) / h[i]
- 3d * (y[i] - y[i - 1]) / h[i - 1];
}
r[n-1] = 0;
f[n-1] = 0;
for (int i = 1; i <= n-1; i++) {
u[i] = sqr(r[i-1])*sigma(i-1) + sqr(f[i])*sigma(i) + sqr(r[i])*sigma(i + 1);
u[i] = lambda * u[i] + p[i];
if(i u=p, v=h, w=0
quincunx(n, u, v, w, q);
//Spline Parameters
residuals = smoothing = 0;
double a[] = new double[n];
double b[] = new double[n];
double c[] = new double[n];
double d[] = new double[n];
d[0] = y[0] - lambda * r[0] * q[1] * sigma(0);
d[1] = y[1] - lambda * (f[1] * q[1] + r[1] * q[2]) * sigma(1);
a[0] = q[1] / (3d * h[0]);
b[0] = 0;
c[0] = (d[1] - d[0])/h[0] - q[1] * h[0]/3d;
r[0] = 0;
for (int j = 1; j < n-1; j++) {
a[j] = (q[j + 1]-q[j]) / (3d * h[j]);
b[j] = q[j];
c[j] = (q[j] + q[j-1]) * h[j-1] + c[j-1];
d[j] = r[j-1] * q[j-1] + f[j] * q[j] + r[j] * q[j+1];
d[j] = y[j] - lambda * d[j] * sigma(j);
}
for (int j = 0; j < n-1; j++) {
if(sigma(j)>0) residuals += (d[j]-y[j]) * (d[j]-y[j]) / (sigma(j) * sigma(j));
smoothing += a[j]*a[j];
}
FastSpline.Polynome polynomials[] = new FastSpline.Polynome[n-1];
for (int i = 0; i < polynomials.length; i++) {
polynomials[i] = new FastSpline.Polynome(new double[] {d[i], c[i], b[i], a[i]});
}
return new FastSpline(x, polynomials);
}
public double getResiduals() {
return residuals;
}
public double getSmoothing() {
return smoothing;
}
public double getLambda() {
return lambda;
}
public double[] getSigma() {
return sigma;
}
public void setLambda(double d) {
lambda = d;
}
public void setSigma(double[] ds) {
sigma = ds;
}
private static final double sqr(double v) {
return v*v;
}
private final double sigma(int index) {
if(sigma==null || index>=sigma.length) return 1;
return sigma[index];
}
}