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;
18
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;
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 * 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;
60
61 /**
62 * Uses the original Akima algorithm.
63 */
64 public AkimaSplineInterpolator() {
65 this(false);
66 }
67
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 }
76
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 }
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 }