AkimaSplineInterpolator.java
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You 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 org.apache.commons.math4.legacy.analysis.interpolation;
import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunction;
import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
import org.apache.commons.math4.legacy.exception.NullArgumentException;
import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
import org.apache.commons.math4.core.jdkmath.JdkMath;
import org.apache.commons.math4.legacy.core.MathArrays;
import org.apache.commons.numbers.core.Precision;
/**
* Computes a cubic spline interpolation for the data set using the Akima
* algorithm, as originally formulated by Hiroshi Akima in his 1970 paper
* "A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures."
* J. ACM 17, 4 (October 1970), 589-602. DOI=10.1145/321607.321609
* http://doi.acm.org/10.1145/321607.321609
* <p>
* This implementation is based on the Akima implementation in the CubicSpline
* class in the Math.NET Numerics library. The method referenced is
* "CubicSpline.InterpolateAkimaSorted".
* </p>
* <p>
* When the {@link #AkimaSplineInterpolator(boolean) argument to the constructor}
* is {@code true}, the weights are modified in order to
* <a href="https://nl.mathworks.com/help/matlab/ref/makima.html">smooth out
* spurious oscillations</a>.
* </p>
* <p>
* The {@link #interpolate(double[], double[]) interpolate} method returns a
* {@link PolynomialSplineFunction} consisting of n cubic polynomials, defined
* over the subintervals determined by the x values, {@code x[0] < x[i] ... < x[n]}.
* The Akima algorithm requires that {@code n >= 5}.
* </p>
*/
public class AkimaSplineInterpolator
implements UnivariateInterpolator {
/** The minimum number of points that are needed to compute the function. */
private static final int MINIMUM_NUMBER_POINTS = 5;
/** Whether to use the "modified Akima" interpolation. */
private final boolean useModified;
/**
* Uses the original Akima algorithm.
*/
public AkimaSplineInterpolator() {
this(false);
}
/**
* @param useModified Whether to use the
* <a href="https://nl.mathworks.com/help/matlab/ref/makima.html">modified Akima</a>
* interpolation.
*/
public AkimaSplineInterpolator(boolean useModified) {
this.useModified = useModified;
}
/**
* Computes an interpolating function for the data set.
*
* @param xvals the arguments for the interpolation points
* @param yvals the values for the interpolation points
* @return a function which interpolates the data set
* @throws DimensionMismatchException if {@code xvals} and {@code yvals} have
* different sizes.
* @throws NonMonotonicSequenceException if {@code xvals} is not sorted in
* strict increasing order.
* @throws NumberIsTooSmallException if the size of {@code xvals} is smaller
* than 5.
*/
@Override
public PolynomialSplineFunction interpolate(double[] xvals,
double[] yvals)
throws DimensionMismatchException,
NumberIsTooSmallException,
NonMonotonicSequenceException {
if (xvals == null ||
yvals == null) {
throw new NullArgumentException();
}
if (xvals.length != yvals.length) {
throw new DimensionMismatchException(xvals.length, yvals.length);
}
if (xvals.length < MINIMUM_NUMBER_POINTS) {
throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS,
xvals.length,
MINIMUM_NUMBER_POINTS, true);
}
MathArrays.checkOrder(xvals);
final int numberOfDiffAndWeightElements = xvals.length - 1;
final double[] differences = new double[numberOfDiffAndWeightElements];
final double[] weights = new double[numberOfDiffAndWeightElements];
for (int i = 0; i < differences.length; i++) {
differences[i] = (yvals[i + 1] - yvals[i]) / (xvals[i + 1] - xvals[i]);
}
if (useModified) {
for (int i = 1; i < weights.length; i++) {
final double a = differences[i];
final double b = differences[i - 1];
weights[i] = JdkMath.abs(a - b) + 0.5 * JdkMath.abs(a + b);
}
} else {
for (int i = 1; i < weights.length; i++) {
weights[i] = JdkMath.abs(differences[i] - differences[i - 1]);
}
}
// Prepare Hermite interpolation scheme.
final double[] firstDerivatives = new double[xvals.length];
for (int i = 2; i < firstDerivatives.length - 2; i++) {
final double wP = weights[i + 1];
final double wM = weights[i - 1];
if (Precision.equals(wP, 0.0) &&
Precision.equals(wM, 0.0)) {
final double xv = xvals[i];
final double xvP = xvals[i + 1];
final double xvM = xvals[i - 1];
firstDerivatives[i] = (((xvP - xv) * differences[i - 1]) + ((xv - xvM) * differences[i])) / (xvP - xvM);
} else {
firstDerivatives[i] = ((wP * differences[i - 1]) + (wM * differences[i])) / (wP + wM);
}
}
firstDerivatives[0] = differentiateThreePoint(xvals, yvals, 0, 0, 1, 2);
firstDerivatives[1] = differentiateThreePoint(xvals, yvals, 1, 0, 1, 2);
firstDerivatives[xvals.length - 2] = differentiateThreePoint(xvals, yvals, xvals.length - 2,
xvals.length - 3, xvals.length - 2,
xvals.length - 1);
firstDerivatives[xvals.length - 1] = differentiateThreePoint(xvals, yvals, xvals.length - 1,
xvals.length - 3, xvals.length - 2,
xvals.length - 1);
return interpolateHermiteSorted(xvals, yvals, firstDerivatives);
}
/**
* Three point differentiation helper, modeled off of the same method in the
* Math.NET CubicSpline class. This is used by both the Apache Math and the
* Math.NET Akima Cubic Spline algorithms
*
* @param xvals x values to calculate the numerical derivative with
* @param yvals y values to calculate the numerical derivative with
* @param indexOfDifferentiation index of the elemnt we are calculating the derivative around
* @param indexOfFirstSample index of the first element to sample for the three point method
* @param indexOfSecondsample index of the second element to sample for the three point method
* @param indexOfThirdSample index of the third element to sample for the three point method
* @return the derivative
*/
private double differentiateThreePoint(double[] xvals, double[] yvals,
int indexOfDifferentiation,
int indexOfFirstSample,
int indexOfSecondsample,
int indexOfThirdSample) {
final double x0 = yvals[indexOfFirstSample];
final double x1 = yvals[indexOfSecondsample];
final double x2 = yvals[indexOfThirdSample];
final double t = xvals[indexOfDifferentiation] - xvals[indexOfFirstSample];
final double t1 = xvals[indexOfSecondsample] - xvals[indexOfFirstSample];
final double t2 = xvals[indexOfThirdSample] - xvals[indexOfFirstSample];
final double a = (x2 - x0 - (t2 / t1 * (x1 - x0))) / (t2 * t2 - t1 * t2);
final double b = (x1 - x0 - a * t1 * t1) / t1;
return (2 * a * t) + b;
}
/**
* Creates a Hermite cubic spline interpolation from the set of (x,y) value
* pairs and their derivatives. This is modeled off of the
* InterpolateHermiteSorted method in the Math.NET CubicSpline class.
*
* @param xvals x values for interpolation
* @param yvals y values for interpolation
* @param firstDerivatives first derivative values of the function
* @return polynomial that fits the function
*/
private PolynomialSplineFunction interpolateHermiteSorted(double[] xvals,
double[] yvals,
double[] firstDerivatives) {
if (xvals.length != yvals.length) {
throw new DimensionMismatchException(xvals.length, yvals.length);
}
if (xvals.length != firstDerivatives.length) {
throw new DimensionMismatchException(xvals.length,
firstDerivatives.length);
}
final int minimumLength = 2;
if (xvals.length < minimumLength) {
throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS,
xvals.length, minimumLength,
true);
}
final int size = xvals.length - 1;
final PolynomialFunction[] polynomials = new PolynomialFunction[size];
final double[] coefficients = new double[4];
for (int i = 0; i < polynomials.length; i++) {
final double w = xvals[i + 1] - xvals[i];
final double w2 = w * w;
final double yv = yvals[i];
final double yvP = yvals[i + 1];
final double fd = firstDerivatives[i];
final double fdP = firstDerivatives[i + 1];
coefficients[0] = yv;
coefficients[1] = firstDerivatives[i];
coefficients[2] = (3 * (yvP - yv) / w - 2 * fd - fdP) / w;
coefficients[3] = (2 * (yv - yvP) / w + fd + fdP) / w2;
polynomials[i] = new PolynomialFunction(coefficients);
}
return new PolynomialSplineFunction(xvals, polynomials);
}
}