001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math4.legacy.analysis.interpolation;
018
019import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunction;
020import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialSplineFunction;
021import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
022import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
023import org.apache.commons.math4.legacy.exception.NullArgumentException;
024import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
025import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
026import org.apache.commons.math4.core.jdkmath.JdkMath;
027import org.apache.commons.math4.legacy.core.MathArrays;
028import org.apache.commons.numbers.core.Precision;
029
030/**
031 * Computes a cubic spline interpolation for the data set using the Akima
032 * algorithm, as originally formulated by Hiroshi Akima in his 1970 paper
033 * "A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures."
034 * J. ACM 17, 4 (October 1970), 589-602. DOI=10.1145/321607.321609
035 * http://doi.acm.org/10.1145/321607.321609
036 * <p>
037 * This implementation is based on the Akima implementation in the CubicSpline
038 * class in the Math.NET Numerics library. The method referenced is
039 * "CubicSpline.InterpolateAkimaSorted".
040 * </p>
041 * <p>
042 * When the {@link #AkimaSplineInterpolator(boolean) argument to the constructor}
043 * is {@code true}, the weights are modified in order to
044 * <a href="https://nl.mathworks.com/help/matlab/ref/makima.html">smooth out
045 * spurious oscillations</a>.
046 * </p>
047 * <p>
048 * The {@link #interpolate(double[], double[]) interpolate} method returns a
049 * {@link PolynomialSplineFunction} consisting of n cubic polynomials, defined
050 * over the subintervals determined by the x values, {@code x[0] < x[i] ... < x[n]}.
051 * The Akima algorithm requires that {@code n >= 5}.
052 * </p>
053 */
054public class AkimaSplineInterpolator
055    implements UnivariateInterpolator {
056    /** The minimum number of points that are needed to compute the function. */
057    private static final int MINIMUM_NUMBER_POINTS = 5;
058    /** Whether to use the "modified Akima" interpolation. */
059    private final boolean useModified;
060
061    /**
062     * Uses the original Akima algorithm.
063     */
064    public AkimaSplineInterpolator() {
065        this(false);
066    }
067
068    /**
069     * @param useModified Whether to use the
070     * <a href="https://nl.mathworks.com/help/matlab/ref/makima.html">modified Akima</a>
071     * interpolation.
072     */
073    public AkimaSplineInterpolator(boolean useModified) {
074        this.useModified = useModified;
075    }
076
077    /**
078     * Computes an interpolating function for the data set.
079     *
080     * @param xvals the arguments for the interpolation points
081     * @param yvals the values for the interpolation points
082     * @return a function which interpolates the data set
083     * @throws DimensionMismatchException if {@code xvals} and {@code yvals} have
084     *         different sizes.
085     * @throws NonMonotonicSequenceException if {@code xvals} is not sorted in
086     *         strict increasing order.
087     * @throws NumberIsTooSmallException if the size of {@code xvals} is smaller
088     *         than 5.
089     */
090    @Override
091    public PolynomialSplineFunction interpolate(double[] xvals,
092                                                double[] yvals)
093        throws DimensionMismatchException,
094               NumberIsTooSmallException,
095               NonMonotonicSequenceException {
096        if (xvals == null ||
097            yvals == null) {
098            throw new NullArgumentException();
099        }
100
101        if (xvals.length != yvals.length) {
102            throw new DimensionMismatchException(xvals.length, yvals.length);
103        }
104
105        if (xvals.length < MINIMUM_NUMBER_POINTS) {
106            throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS,
107                                                xvals.length,
108                                                MINIMUM_NUMBER_POINTS, true);
109        }
110
111        MathArrays.checkOrder(xvals);
112
113        final int numberOfDiffAndWeightElements = xvals.length - 1;
114
115        final double[] differences = new double[numberOfDiffAndWeightElements];
116        final double[] weights = new double[numberOfDiffAndWeightElements];
117
118        for (int i = 0; i < differences.length; i++) {
119            differences[i] = (yvals[i + 1] - yvals[i]) / (xvals[i + 1] - xvals[i]);
120        }
121
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        }
133
134        // Prepare Hermite interpolation scheme.
135        final double[] firstDerivatives = new double[xvals.length];
136
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        }
150
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);
159
160        return interpolateHermiteSorted(xvals, yvals, firstDerivatives);
161    }
162
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];
184
185        final double t = xvals[indexOfDifferentiation] - xvals[indexOfFirstSample];
186        final double t1 = xvals[indexOfSecondsample] - xvals[indexOfFirstSample];
187        final double t2 = xvals[indexOfThirdSample] - xvals[indexOfFirstSample];
188
189        final double a = (x2 - x0 - (t2 / t1 * (x1 - x0))) / (t2 * t2 - t1 * t2);
190        final double b = (x1 - x0 - a * t1 * t1) / t1;
191
192        return (2 * a * t) + b;
193    }
194
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        }
211
212        if (xvals.length != firstDerivatives.length) {
213            throw new DimensionMismatchException(xvals.length,
214                                                 firstDerivatives.length);
215        }
216
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        }
223
224        final int size = xvals.length - 1;
225        final PolynomialFunction[] polynomials = new PolynomialFunction[size];
226        final double[] coefficients = new double[4];
227
228        for (int i = 0; i < polynomials.length; i++) {
229            final double w = xvals[i + 1] - xvals[i];
230            final double w2 = w * w;
231
232            final double yv = yvals[i];
233            final double yvP = yvals[i + 1];
234
235            final double fd = firstDerivatives[i];
236            final double fdP = firstDerivatives[i + 1];
237
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        }
244
245        return new PolynomialSplineFunction(xvals, polynomials);
246    }
247}