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.math3.analysis.interpolation; 18 19 import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; 20 import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction; 21 import org.apache.commons.math3.exception.DimensionMismatchException; 22 import org.apache.commons.math3.exception.NonMonotonicSequenceException; 23 import org.apache.commons.math3.exception.NullArgumentException; 24 import org.apache.commons.math3.exception.NumberIsTooSmallException; 25 import org.apache.commons.math3.exception.util.LocalizedFormats; 26 import org.apache.commons.math3.util.FastMath; 27 import org.apache.commons.math3.util.MathArrays; 28 import org.apache.commons.math3.util.Precision; 29 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 * The {@link #interpolate(double[], double[]) interpolate} method returns a 43 * {@link PolynomialSplineFunction} consisting of n cubic polynomials, defined 44 * over the subintervals determined by the x values, {@code x[0] < x[i] ... < x[n]}. 45 * The Akima algorithm requires that {@code n >= 5}. 46 * </p> 47 */ 48 public class AkimaSplineInterpolator 49 implements UnivariateInterpolator { 50 /** The minimum number of points that are needed to compute the function. */ 51 private static final int MINIMUM_NUMBER_POINTS = 5; 52 53 /** 54 * Computes an interpolating function for the data set. 55 * 56 * @param xvals the arguments for the interpolation points 57 * @param yvals the values for the interpolation points 58 * @return a function which interpolates the data set 59 * @throws DimensionMismatchException if {@code xvals} and {@code yvals} have 60 * different sizes. 61 * @throws NonMonotonicSequenceException if {@code xvals} is not sorted in 62 * strict increasing order. 63 * @throws NumberIsTooSmallException if the size of {@code xvals} is smaller 64 * than 5. 65 */ interpolate(double[] xvals, double[] yvals)66 public PolynomialSplineFunction interpolate(double[] xvals, 67 double[] yvals) 68 throws DimensionMismatchException, 69 NumberIsTooSmallException, 70 NonMonotonicSequenceException { 71 if (xvals == null || 72 yvals == null) { 73 throw new NullArgumentException(); 74 } 75 76 if (xvals.length != yvals.length) { 77 throw new DimensionMismatchException(xvals.length, yvals.length); 78 } 79 80 if (xvals.length < MINIMUM_NUMBER_POINTS) { 81 throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS, 82 xvals.length, 83 MINIMUM_NUMBER_POINTS, true); 84 } 85 86 MathArrays.checkOrder(xvals); 87 88 final int numberOfDiffAndWeightElements = xvals.length - 1; 89 90 final double[] differences = new double[numberOfDiffAndWeightElements]; 91 final double[] weights = new double[numberOfDiffAndWeightElements]; 92 93 for (int i = 0; i < differences.length; i++) { 94 differences[i] = (yvals[i + 1] - yvals[i]) / (xvals[i + 1] - xvals[i]); 95 } 96 97 for (int i = 1; i < weights.length; i++) { 98 weights[i] = FastMath.abs(differences[i] - differences[i - 1]); 99 } 100 101 // Prepare Hermite interpolation scheme. 102 final double[] firstDerivatives = new double[xvals.length]; 103 104 for (int i = 2; i < firstDerivatives.length - 2; i++) { 105 final double wP = weights[i + 1]; 106 final double wM = weights[i - 1]; 107 if (Precision.equals(wP, 0.0) && 108 Precision.equals(wM, 0.0)) { 109 final double xv = xvals[i]; 110 final double xvP = xvals[i + 1]; 111 final double xvM = xvals[i - 1]; 112 firstDerivatives[i] = (((xvP - xv) * differences[i - 1]) + ((xv - xvM) * differences[i])) / (xvP - xvM); 113 } else { 114 firstDerivatives[i] = ((wP * differences[i - 1]) + (wM * differences[i])) / (wP + wM); 115 } 116 } 117 118 firstDerivatives[0] = differentiateThreePoint(xvals, yvals, 0, 0, 1, 2); 119 firstDerivatives[1] = differentiateThreePoint(xvals, yvals, 1, 0, 1, 2); 120 firstDerivatives[xvals.length - 2] = differentiateThreePoint(xvals, yvals, xvals.length - 2, 121 xvals.length - 3, xvals.length - 2, 122 xvals.length - 1); 123 firstDerivatives[xvals.length - 1] = differentiateThreePoint(xvals, yvals, xvals.length - 1, 124 xvals.length - 3, xvals.length - 2, 125 xvals.length - 1); 126 127 return interpolateHermiteSorted(xvals, yvals, firstDerivatives); 128 } 129 130 /** 131 * Three point differentiation helper, modeled off of the same method in the 132 * Math.NET CubicSpline class. This is used by both the Apache Math and the 133 * Math.NET Akima Cubic Spline algorithms 134 * 135 * @param xvals x values to calculate the numerical derivative with 136 * @param yvals y values to calculate the numerical derivative with 137 * @param indexOfDifferentiation index of the elemnt we are calculating the derivative around 138 * @param indexOfFirstSample index of the first element to sample for the three point method 139 * @param indexOfSecondsample index of the second element to sample for the three point method 140 * @param indexOfThirdSample index of the third element to sample for the three point method 141 * @return the derivative 142 */ differentiateThreePoint(double[] xvals, double[] yvals, int indexOfDifferentiation, int indexOfFirstSample, int indexOfSecondsample, int indexOfThirdSample)143 private double differentiateThreePoint(double[] xvals, double[] yvals, 144 int indexOfDifferentiation, 145 int indexOfFirstSample, 146 int indexOfSecondsample, 147 int indexOfThirdSample) { 148 final double x0 = yvals[indexOfFirstSample]; 149 final double x1 = yvals[indexOfSecondsample]; 150 final double x2 = yvals[indexOfThirdSample]; 151 152 final double t = xvals[indexOfDifferentiation] - xvals[indexOfFirstSample]; 153 final double t1 = xvals[indexOfSecondsample] - xvals[indexOfFirstSample]; 154 final double t2 = xvals[indexOfThirdSample] - xvals[indexOfFirstSample]; 155 156 final double a = (x2 - x0 - (t2 / t1 * (x1 - x0))) / (t2 * t2 - t1 * t2); 157 final double b = (x1 - x0 - a * t1 * t1) / t1; 158 159 return (2 * a * t) + b; 160 } 161 162 /** 163 * Creates a Hermite cubic spline interpolation from the set of (x,y) value 164 * pairs and their derivatives. This is modeled off of the 165 * InterpolateHermiteSorted method in the Math.NET CubicSpline class. 166 * 167 * @param xvals x values for interpolation 168 * @param yvals y values for interpolation 169 * @param firstDerivatives first derivative values of the function 170 * @return polynomial that fits the function 171 */ interpolateHermiteSorted(double[] xvals, double[] yvals, double[] firstDerivatives)172 private PolynomialSplineFunction interpolateHermiteSorted(double[] xvals, 173 double[] yvals, 174 double[] firstDerivatives) { 175 if (xvals.length != yvals.length) { 176 throw new DimensionMismatchException(xvals.length, yvals.length); 177 } 178 179 if (xvals.length != firstDerivatives.length) { 180 throw new DimensionMismatchException(xvals.length, 181 firstDerivatives.length); 182 } 183 184 final int minimumLength = 2; 185 if (xvals.length < minimumLength) { 186 throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS, 187 xvals.length, minimumLength, 188 true); 189 } 190 191 final int size = xvals.length - 1; 192 final PolynomialFunction[] polynomials = new PolynomialFunction[size]; 193 final double[] coefficients = new double[4]; 194 195 for (int i = 0; i < polynomials.length; i++) { 196 final double w = xvals[i + 1] - xvals[i]; 197 final double w2 = w * w; 198 199 final double yv = yvals[i]; 200 final double yvP = yvals[i + 1]; 201 202 final double fd = firstDerivatives[i]; 203 final double fdP = firstDerivatives[i + 1]; 204 205 coefficients[0] = yv; 206 coefficients[1] = firstDerivatives[i]; 207 coefficients[2] = (3 * (yvP - yv) / w - 2 * fd - fdP) / w; 208 coefficients[3] = (2 * (yv - yvP) / w + fd + fdP) / w2; 209 polynomials[i] = new PolynomialFunction(coefficients); 210 } 211 212 return new PolynomialSplineFunction(xvals, polynomials); 213 214 } 215 } 216