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}