AkimaSplineInterpolator.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math4.legacy.analysis.interpolation;

  19. import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunction;
  20. import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialSplineFunction;
  21. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  22. import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
  23. import org.apache.commons.math4.legacy.exception.NullArgumentException;
  24. import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
  25. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  26. import org.apache.commons.math4.core.jdkmath.JdkMath;
  27. import org.apache.commons.math4.legacy.core.MathArrays;
  28. import org.apache.commons.numbers.core.Precision;

  30. /**
  31.  * Computes a cubic spline interpolation for the data set using the Akima
  32.  * algorithm, as originally formulated by Hiroshi Akima in his 1970 paper
  33.  * "A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures."
  34.  * J. ACM 17, 4 (October 1970), 589-602. DOI=10.1145/321607.321609
  35.  * http://doi.acm.org/10.1145/321607.321609
  36.  * <p>
  37.  * This implementation is based on the Akima implementation in the CubicSpline
  38.  * class in the Math.NET Numerics library. The method referenced is
  39.  * "CubicSpline.InterpolateAkimaSorted".
  40.  * </p>
  41.  * <p>
  42.  * When the {@link #AkimaSplineInterpolator(boolean) argument to the constructor}
  43.  * is {@code true}, the weights are modified in order to
  44.  * <a href="https://nl.mathworks.com/help/matlab/ref/makima.html">smooth out
  45.  * spurious oscillations</a>.
  46.  * </p>
  47.  * <p>
  48.  * The {@link #interpolate(double[], double[]) interpolate} method returns a
  49.  * {@link PolynomialSplineFunction} consisting of n cubic polynomials, defined
  50.  * over the subintervals determined by the x values, {@code x[0] < x[i] ... < x[n]}.
  51.  * The Akima algorithm requires that {@code n >= 5}.
  52.  * </p>
  53.  */
  54. public class AkimaSplineInterpolator
  55.     implements UnivariateInterpolator {
  56.     /** The minimum number of points that are needed to compute the function. */
  57.     private static final int MINIMUM_NUMBER_POINTS = 5;
  58.     /** Whether to use the "modified Akima" interpolation. */
  59.     private final boolean useModified;

  61.     /**
  62.      * Uses the original Akima algorithm.
  63.      */
  64.     public AkimaSplineInterpolator() {
  65.         this(false);
  66.     }

  68.     /**
  69.      * @param useModified Whether to use the
  70.      * <a href="https://nl.mathworks.com/help/matlab/ref/makima.html">modified Akima</a>
  71.      * interpolation.
  72.      */
  73.     public AkimaSplineInterpolator(boolean useModified) {
  74.         this.useModified = useModified;
  75.     }

  77.     /**
  78.      * Computes an interpolating function for the data set.
  79.      *
  80.      * @param xvals the arguments for the interpolation points
  81.      * @param yvals the values for the interpolation points
  82.      * @return a function which interpolates the data set
  83.      * @throws DimensionMismatchException if {@code xvals} and {@code yvals} have
  84.      *         different sizes.
  85.      * @throws NonMonotonicSequenceException if {@code xvals} is not sorted in
  86.      *         strict increasing order.
  87.      * @throws NumberIsTooSmallException if the size of {@code xvals} is smaller
  88.      *         than 5.
  89.      */
  90.     @Override
  91.     public PolynomialSplineFunction interpolate(double[] xvals,
  92.                                                 double[] yvals)
  93.         throws DimensionMismatchException,
  94.                NumberIsTooSmallException,
  95.                NonMonotonicSequenceException {
  96.         if (xvals == null ||
  97.             yvals == null) {
  98.             throw new NullArgumentException();
  99.         }

  101.         if (xvals.length != yvals.length) {
  102.             throw new DimensionMismatchException(xvals.length, yvals.length);
  103.         }

  105.         if (xvals.length < MINIMUM_NUMBER_POINTS) {
  106.             throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS,
  107.                                                 xvals.length,
  108.                                                 MINIMUM_NUMBER_POINTS, true);
  109.         }

  111.         MathArrays.checkOrder(xvals);

  113.         final int numberOfDiffAndWeightElements = xvals.length - 1;

  115.         final double[] differences = new double[numberOfDiffAndWeightElements];
  116.         final double[] weights = new double[numberOfDiffAndWeightElements];

  118.         for (int i = 0; i < differences.length; i++) {
  119.             differences[i] = (yvals[i + 1] - yvals[i]) / (xvals[i + 1] - xvals[i]);
  120.         }

  122.         if (useModified) {
  123.             for (int i = 1; i < weights.length; i++) {
  124.                 final double a = differences[i];
  125.                 final double b = differences[i - 1];
  126.                 weights[i] = JdkMath.abs(a - b) + 0.5 * JdkMath.abs(a + b);
  127.             }
  128.         } else {
  129.             for (int i = 1; i < weights.length; i++) {
  130.                 weights[i] = JdkMath.abs(differences[i] - differences[i - 1]);
  131.             }
  132.         }

  134.         // Prepare Hermite interpolation scheme.
  135.         final double[] firstDerivatives = new double[xvals.length];

  137.         for (int i = 2; i < firstDerivatives.length - 2; i++) {
  138.             final double wP = weights[i + 1];
  139.             final double wM = weights[i - 1];
  140.             if (Precision.equals(wP, 0.0) &&
  141.                 Precision.equals(wM, 0.0)) {
  142.                 final double xv = xvals[i];
  143.                 final double xvP = xvals[i + 1];
  144.                 final double xvM = xvals[i - 1];
  145.                 firstDerivatives[i] = (((xvP - xv) * differences[i - 1]) + ((xv - xvM) * differences[i])) / (xvP - xvM);
  146.             } else {
  147.                 firstDerivatives[i] = ((wP * differences[i - 1]) + (wM * differences[i])) / (wP + wM);
  148.             }
  149.         }

  151.         firstDerivatives[0] = differentiateThreePoint(xvals, yvals, 0, 0, 1, 2);
  152.         firstDerivatives[1] = differentiateThreePoint(xvals, yvals, 1, 0, 1, 2);
  153.         firstDerivatives[xvals.length - 2] = differentiateThreePoint(xvals, yvals, xvals.length - 2,
  154.                                                                      xvals.length - 3, xvals.length - 2,
  155.                                                                      xvals.length - 1);
  156.         firstDerivatives[xvals.length - 1] = differentiateThreePoint(xvals, yvals, xvals.length - 1,
  157.                                                                      xvals.length - 3, xvals.length - 2,
  158.                                                                      xvals.length - 1);

  160.         return interpolateHermiteSorted(xvals, yvals, firstDerivatives);
  161.     }

  163.     /**
  164.      * Three point differentiation helper, modeled off of the same method in the
  165.      * Math.NET CubicSpline class. This is used by both the Apache Math and the
  166.      * Math.NET Akima Cubic Spline algorithms
  167.      *
  168.      * @param xvals x values to calculate the numerical derivative with
  169.      * @param yvals y values to calculate the numerical derivative with
  170.      * @param indexOfDifferentiation index of the elemnt we are calculating the derivative around
  171.      * @param indexOfFirstSample index of the first element to sample for the three point method
  172.      * @param indexOfSecondsample index of the second element to sample for the three point method
  173.      * @param indexOfThirdSample index of the third element to sample for the three point method
  174.      * @return the derivative
  175.      */
  176.     private double differentiateThreePoint(double[] xvals, double[] yvals,
  177.                                            int indexOfDifferentiation,
  178.                                            int indexOfFirstSample,
  179.                                            int indexOfSecondsample,
  180.                                            int indexOfThirdSample) {
  181.         final double x0 = yvals[indexOfFirstSample];
  182.         final double x1 = yvals[indexOfSecondsample];
  183.         final double x2 = yvals[indexOfThirdSample];

  185.         final double t = xvals[indexOfDifferentiation] - xvals[indexOfFirstSample];
  186.         final double t1 = xvals[indexOfSecondsample] - xvals[indexOfFirstSample];
  187.         final double t2 = xvals[indexOfThirdSample] - xvals[indexOfFirstSample];

  189.         final double a = (x2 - x0 - (t2 / t1 * (x1 - x0))) / (t2 * t2 - t1 * t2);
  190.         final double b = (x1 - x0 - a * t1 * t1) / t1;

  192.         return (2 * a * t) + b;
  193.     }

  195.     /**
  196.      * Creates a Hermite cubic spline interpolation from the set of (x,y) value
  197.      * pairs and their derivatives. This is modeled off of the
  198.      * InterpolateHermiteSorted method in the Math.NET CubicSpline class.
  199.      *
  200.      * @param xvals x values for interpolation
  201.      * @param yvals y values for interpolation
  202.      * @param firstDerivatives first derivative values of the function
  203.      * @return polynomial that fits the function
  204.      */
  205.     private PolynomialSplineFunction interpolateHermiteSorted(double[] xvals,
  206.                                                               double[] yvals,
  207.                                                               double[] firstDerivatives) {
  208.         if (xvals.length != yvals.length) {
  209.             throw new DimensionMismatchException(xvals.length, yvals.length);
  210.         }

  212.         if (xvals.length != firstDerivatives.length) {
  213.             throw new DimensionMismatchException(xvals.length,
  214.                                                  firstDerivatives.length);
  215.         }

  217.         final int minimumLength = 2;
  218.         if (xvals.length < minimumLength) {
  219.             throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS,
  220.                                                 xvals.length, minimumLength,
  221.                                                 true);
  222.         }

  224.         final int size = xvals.length - 1;
  225.         final PolynomialFunction[] polynomials = new PolynomialFunction[size];
  226.         final double[] coefficients = new double[4];

  228.         for (int i = 0; i < polynomials.length; i++) {
  229.             final double w = xvals[i + 1] - xvals[i];
  230.             final double w2 = w * w;

  232.             final double yv = yvals[i];
  233.             final double yvP = yvals[i + 1];

  235.             final double fd = firstDerivatives[i];
  236.             final double fdP = firstDerivatives[i + 1];

  238.             coefficients[0] = yv;
  239.             coefficients[1] = firstDerivatives[i];
  240.             coefficients[2] = (3 * (yvP - yv) / w - 2 * fd - fdP) / w;
  241.             coefficients[3] = (2 * (yv - yvP) / w + fd + fdP) / w2;
  242.             polynomials[i] = new PolynomialFunction(coefficients);
  243.         }

  245.         return new PolynomialSplineFunction(xvals, polynomials);
  246.     }
  247. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT