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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.analysis.interpolation;
23
24 import java.lang.reflect.Array;
25
26 import org.hipparchus.Field;
27 import org.hipparchus.CalculusFieldElement;
28 import org.hipparchus.analysis.polynomials.FieldPolynomialFunction;
29 import org.hipparchus.analysis.polynomials.FieldPolynomialSplineFunction;
30 import org.hipparchus.analysis.polynomials.PolynomialFunction;
31 import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
32 import org.hipparchus.exception.LocalizedCoreFormats;
33 import org.hipparchus.exception.MathIllegalArgumentException;
34 import org.hipparchus.exception.NullArgumentException;
35 import org.hipparchus.util.FastMath;
36 import org.hipparchus.util.MathArrays;
37 import org.hipparchus.util.Precision;
38
39 /**
40 * Computes a cubic spline interpolation for the data set using the Akima
41 * algorithm, as originally formulated by Hiroshi Akima in his 1970 paper
42 * "A New Method of Interpolation and Smooth Curve Fitting Based on Local Procedures."
43 * J. ACM 17, 4 (October 1970), 589-602. DOI=10.1145/321607.321609
44 * http://doi.acm.org/10.1145/321607.321609
45 * <p>
46 * This implementation is based on the Akima implementation in the CubicSpline
47 * class in the Math.NET Numerics library. The method referenced is
48 * CubicSpline.InterpolateAkimaSorted
49 * </p>
50 * <p>
51 * The {@link #interpolate(double[], double[]) interpolate} method returns a
52 * {@link PolynomialSplineFunction} consisting of n cubic polynomials, defined
53 * over the subintervals determined by the x values, {@code x[0] < x[i] ... < x[n]}.
54 * The Akima algorithm requires that {@code n >= 5}.
55 * </p>
56 */
57 public class AkimaSplineInterpolator
58 implements UnivariateInterpolator, FieldUnivariateInterpolator {
59
60 /** The minimum number of points that are needed to compute the function. */
61 private static final int MINIMUM_NUMBER_POINTS = 5;
62
63 /** Weight modifier to avoid overshoots. */
64 private final boolean useModifiedWeights;
65
66 /** Simple constructor.
67 * <p>
68 * This constructor is equivalent to call {@link #AkimaSplineInterpolator(boolean)
69 * AkimaSplineInterpolator(false)}, i.e. to use original Akima weights
70 * </p>
71 * @since 2.1
72 */
73 public AkimaSplineInterpolator() {
74 this(false);
75 }
76
77 /** Simple constructor.
78 * <p>
79 * The weight modification is described in <a
80 * href="https://blogs.mathworks.com/cleve/2019/04/29/makima-piecewise-cubic-interpolation/">
81 * Makima Piecewise Cubic Interpolation</a>. It attempts to avoid overshoots
82 * near near constant slopes sub-samples.
83 * </p>
84 * @param useModifiedWeights if true, use modified weights to avoid overshoots
85 * @since 2.1
86 */
87 public AkimaSplineInterpolator(final boolean useModifiedWeights) {
88 this.useModifiedWeights = useModifiedWeights;
89 }
90
91 /**
92 * Computes an interpolating function for the data set.
93 *
94 * @param xvals the arguments for the interpolation points
95 * @param yvals the values for the interpolation points
96 * @return a function which interpolates the data set
97 * @throws MathIllegalArgumentException if {@code xvals} and {@code yvals} have
98 * different sizes.
99 * @throws MathIllegalArgumentException if {@code xvals} is not sorted in
100 * strict increasing order.
101 * @throws MathIllegalArgumentException if the size of {@code xvals} is smaller
102 * than 5.
103 */
104 @Override
105 public PolynomialSplineFunction interpolate(double[] xvals,
106 double[] yvals)
107 throws MathIllegalArgumentException {
108 if (xvals == null ||
109 yvals == null) {
110 throw new NullArgumentException();
111 }
112
113 MathArrays.checkEqualLength(xvals, yvals);
114
115 if (xvals.length < MINIMUM_NUMBER_POINTS) {
116 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_POINTS,
117 xvals.length,
118 MINIMUM_NUMBER_POINTS, true);
119 }
120
121 MathArrays.checkOrder(xvals);
122
123 final int numberOfDiffAndWeightElements = xvals.length - 1;
124
125 final double[] differences = new double[numberOfDiffAndWeightElements];
126 final double[] weights = new double[numberOfDiffAndWeightElements];
127
128 for (int i = 0; i < differences.length; i++) {
129 differences[i] = (yvals[i + 1] - yvals[i]) / (xvals[i + 1] - xvals[i]);
130 }
131
132 for (int i = 1; i < weights.length; i++) {
133 weights[i] = FastMath.abs(differences[i] - differences[i - 1]);
134 if (useModifiedWeights) {
135 // modify weights to avoid overshoots near constant slopes sub-samples
136 weights[i] += FastMath.abs(differences[i] + differences[i - 1]);
137 }
138 }
139
140 // Prepare Hermite interpolation scheme.
141 final double[] firstDerivatives = new double[xvals.length];
142
143 for (int i = 2; i < firstDerivatives.length - 2; i++) {
144 final double wP = weights[i + 1];
145 final double wM = weights[i - 1];
146 if (Precision.equals(wP, 0.0) &&
147 Precision.equals(wM, 0.0)) {
148 final double xv = xvals[i];
149 final double xvP = xvals[i + 1];
150 final double xvM = xvals[i - 1];
151 firstDerivatives[i] = (((xvP - xv) * differences[i - 1]) + ((xv - xvM) * differences[i])) / (xvP - xvM);
152 } else {
153 firstDerivatives[i] = ((wP * differences[i - 1]) + (wM * differences[i])) / (wP + wM);
154 }
155 }
156
157 firstDerivatives[0] = differentiateThreePoint(xvals, yvals, 0, 0, 1, 2);
158 firstDerivatives[1] = differentiateThreePoint(xvals, yvals, 1, 0, 1, 2);
159 firstDerivatives[xvals.length - 2] = differentiateThreePoint(xvals, yvals, xvals.length - 2,
160 xvals.length - 3, xvals.length - 2,
161 xvals.length - 1);
162 firstDerivatives[xvals.length - 1] = differentiateThreePoint(xvals, yvals, xvals.length - 1,
163 xvals.length - 3, xvals.length - 2,
164 xvals.length - 1);
165
166 return interpolateHermiteSorted(xvals, yvals, firstDerivatives);
167
168 }
169
170 /**
171 * Computes an interpolating function for the data set.
172 *
173 * @param xvals the arguments for the interpolation points
174 * @param yvals the values for the interpolation points
175 * @param <T> the type of the field elements
176 * @return a function which interpolates the data set
177 * @throws MathIllegalArgumentException if {@code xvals} and {@code yvals} have
178 * different sizes.
179 * @throws MathIllegalArgumentException if {@code xvals} is not sorted in
180 * strict increasing order.
181 * @throws MathIllegalArgumentException if the size of {@code xvals} is smaller
182 * than 5.
183 * @since 1.5
184 */
185 @Override
186 public <T extends CalculusFieldElement<T>> FieldPolynomialSplineFunction<T> interpolate(final T[] xvals,
187 final T[] yvals)
188 throws MathIllegalArgumentException {
189 if (xvals == null ||
190 yvals == null) {
191 throw new NullArgumentException();
192 }
193
194 MathArrays.checkEqualLength(xvals, yvals);
195
196 if (xvals.length < MINIMUM_NUMBER_POINTS) {
197 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_POINTS,
198 xvals.length,
199 MINIMUM_NUMBER_POINTS, true);
200 }
201
202 MathArrays.checkOrder(xvals);
203
204 final Field<T> field = xvals[0].getField();
205 final int numberOfDiffAndWeightElements = xvals.length - 1;
206
207 final T[] differences = MathArrays.buildArray(field, numberOfDiffAndWeightElements);
208 final T[] weights = MathArrays.buildArray(field, numberOfDiffAndWeightElements);
209
210 for (int i = 0; i < differences.length; i++) {
211 differences[i] = yvals[i + 1].subtract(yvals[i]).divide(xvals[i + 1].subtract(xvals[i]));
212 }
213
214 for (int i = 1; i < weights.length; i++) {
215 weights[i] = FastMath.abs(differences[i].subtract(differences[i - 1]));
216 }
217
218 // Prepare Hermite interpolation scheme.
219 final T[] firstDerivatives = MathArrays.buildArray(field, xvals.length);
220
221 for (int i = 2; i < firstDerivatives.length - 2; i++) {
222 final T wP = weights[i + 1];
223 final T wM = weights[i - 1];
224 if (Precision.equals(wP.getReal(), 0.0) &&
225 Precision.equals(wM.getReal(), 0.0)) {
226 final T xv = xvals[i];
227 final T xvP = xvals[i + 1];
228 final T xvM = xvals[i - 1];
229 firstDerivatives[i] = xvP.subtract(xv).multiply(differences[i - 1]).
230 add(xv.subtract(xvM).multiply(differences[i])).
231 divide(xvP.subtract(xvM));
232 } else {
233 firstDerivatives[i] = wP.multiply(differences[i - 1]).
234 add(wM.multiply(differences[i])).
235 divide(wP.add(wM));
236 }
237 }
238
239 firstDerivatives[0] = differentiateThreePoint(xvals, yvals, 0, 0, 1, 2);
240 firstDerivatives[1] = differentiateThreePoint(xvals, yvals, 1, 0, 1, 2);
241 firstDerivatives[xvals.length - 2] = differentiateThreePoint(xvals, yvals, xvals.length - 2,
242 xvals.length - 3, xvals.length - 2,
243 xvals.length - 1);
244 firstDerivatives[xvals.length - 1] = differentiateThreePoint(xvals, yvals, xvals.length - 1,
245 xvals.length - 3, xvals.length - 2,
246 xvals.length - 1);
247
248 return interpolateHermiteSorted(xvals, yvals, firstDerivatives);
249
250 }
251
252 /**
253 * Three point differentiation helper, modeled off of the same method in the
254 * Math.NET CubicSpline class. This is used by both the Apache Math and the
255 * Math.NET Akima Cubic Spline algorithms
256 *
257 * @param xvals x values to calculate the numerical derivative with
258 * @param yvals y values to calculate the numerical derivative with
259 * @param indexOfDifferentiation index of the elemnt we are calculating the derivative around
260 * @param indexOfFirstSample index of the first element to sample for the three point method
261 * @param indexOfSecondsample index of the second element to sample for the three point method
262 * @param indexOfThirdSample index of the third element to sample for the three point method
263 * @return the derivative
264 */
265 private double differentiateThreePoint(double[] xvals, double[] yvals,
266 int indexOfDifferentiation,
267 int indexOfFirstSample,
268 int indexOfSecondsample,
269 int indexOfThirdSample) {
270 final double x0 = yvals[indexOfFirstSample];
271 final double x1 = yvals[indexOfSecondsample];
272 final double x2 = yvals[indexOfThirdSample];
273
274 final double t = xvals[indexOfDifferentiation] - xvals[indexOfFirstSample];
275 final double t1 = xvals[indexOfSecondsample] - xvals[indexOfFirstSample];
276 final double t2 = xvals[indexOfThirdSample] - xvals[indexOfFirstSample];
277
278 final double a = (x2 - x0 - (t2 / t1 * (x1 - x0))) / (t2 * t2 - t1 * t2);
279 final double b = (x1 - x0 - a * t1 * t1) / t1;
280
281 return (2 * a * t) + b;
282 }
283
284 /**
285 * Three point differentiation helper, modeled off of the same method in the
286 * Math.NET CubicSpline class. This is used by both the Apache Math and the
287 * Math.NET Akima Cubic Spline algorithms
288 *
289 * @param xvals x values to calculate the numerical derivative with
290 * @param yvals y values to calculate the numerical derivative with
291 * @param <T> the type of the field elements
292 * @param indexOfDifferentiation index of the elemnt we are calculating the derivative around
293 * @param indexOfFirstSample index of the first element to sample for the three point method
294 * @param indexOfSecondsample index of the second element to sample for the three point method
295 * @param indexOfThirdSample index of the third element to sample for the three point method
296 * @return the derivative
297 * @since 1.5
298 */
299 private <T extends CalculusFieldElement<T>> T differentiateThreePoint(T[] xvals, T[] yvals,
300 int indexOfDifferentiation,
301 int indexOfFirstSample,
302 int indexOfSecondsample,
303 int indexOfThirdSample) {
304 final T x0 = yvals[indexOfFirstSample];
305 final T x1 = yvals[indexOfSecondsample];
306 final T x2 = yvals[indexOfThirdSample];
307
308 final T t = xvals[indexOfDifferentiation].subtract(xvals[indexOfFirstSample]);
309 final T t1 = xvals[indexOfSecondsample].subtract(xvals[indexOfFirstSample]);
310 final T t2 = xvals[indexOfThirdSample].subtract(xvals[indexOfFirstSample]);
311
312 final T a = x2.subtract(x0).subtract(t2.divide(t1).multiply(x1.subtract(x0))).
313 divide(t2.multiply(t2).subtract(t1.multiply(t2)));
314 final T b = x1.subtract(x0).subtract(a.multiply(t1).multiply(t1)).divide(t1);
315
316 return a.multiply(t).multiply(2).add(b);
317 }
318
319 /**
320 * Creates a Hermite cubic spline interpolation from the set of (x,y) value
321 * pairs and their derivatives. This is modeled off of the
322 * InterpolateHermiteSorted method in the Math.NET CubicSpline class.
323 *
324 * @param xvals x values for interpolation
325 * @param yvals y values for interpolation
326 * @param firstDerivatives first derivative values of the function
327 * @return polynomial that fits the function
328 */
329 private PolynomialSplineFunction interpolateHermiteSorted(double[] xvals,
330 double[] yvals,
331 double[] firstDerivatives) {
332 MathArrays.checkEqualLength(xvals, yvals);
333 MathArrays.checkEqualLength(xvals, firstDerivatives);
334
335 final int minimumLength = 2;
336 if (xvals.length < minimumLength) {
337 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_POINTS,
338 xvals.length, minimumLength,
339 true);
340 }
341
342 final int size = xvals.length - 1;
343 final PolynomialFunction[] polynomials = new PolynomialFunction[size];
344 final double[] coefficients = new double[4];
345
346 for (int i = 0; i < polynomials.length; i++) {
347 final double w = xvals[i + 1] - xvals[i];
348 final double w2 = w * w;
349
350 final double yv = yvals[i];
351 final double yvP = yvals[i + 1];
352
353 final double fd = firstDerivatives[i];
354 final double fdP = firstDerivatives[i + 1];
355
356 coefficients[0] = yv;
357 coefficients[1] = firstDerivatives[i];
358 coefficients[2] = (3 * (yvP - yv) / w - 2 * fd - fdP) / w;
359 coefficients[3] = (2 * (yv - yvP) / w + fd + fdP) / w2;
360 polynomials[i] = new PolynomialFunction(coefficients);
361 }
362
363 return new PolynomialSplineFunction(xvals, polynomials);
364
365 }
366 /**
367 * Creates a Hermite cubic spline interpolation from the set of (x,y) value
368 * pairs and their derivatives. This is modeled off of the
369 * InterpolateHermiteSorted method in the Math.NET CubicSpline class.
370 *
371 * @param xvals x values for interpolation
372 * @param yvals y values for interpolation
373 * @param firstDerivatives first derivative values of the function
374 * @param <T> the type of the field elements
375 * @return polynomial that fits the function
376 * @since 1.5
377 */
378 private <T extends CalculusFieldElement<T>> FieldPolynomialSplineFunction<T> interpolateHermiteSorted(T[] xvals,
379 T[] yvals,
380 T[] firstDerivatives) {
381 MathArrays.checkEqualLength(xvals, yvals);
382 MathArrays.checkEqualLength(xvals, firstDerivatives);
383
384 final int minimumLength = 2;
385 if (xvals.length < minimumLength) {
386 throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_POINTS,
387 xvals.length, minimumLength,
388 true);
389 }
390
391 final Field<T> field = xvals[0].getField();
392 final int size = xvals.length - 1;
393 @SuppressWarnings("unchecked")
394 final FieldPolynomialFunction<T>[] polynomials =
395 (FieldPolynomialFunction<T>[]) Array.newInstance(FieldPolynomialFunction.class, size);
396 final T[] coefficients = MathArrays.buildArray(field, 4);
397
398 for (int i = 0; i < polynomials.length; i++) {
399 final T w = xvals[i + 1].subtract(xvals[i]);
400 final T w2 = w.multiply(w);
401
402 final T yv = yvals[i];
403 final T yvP = yvals[i + 1];
404
405 final T fd = firstDerivatives[i];
406 final T fdP = firstDerivatives[i + 1];
407
408 coefficients[0] = yv;
409 coefficients[1] = firstDerivatives[i];
410 final T ratio = yvP.subtract(yv).divide(w);
411 coefficients[2] = ratio.multiply(+3).subtract(fd.add(fd)).subtract(fdP).divide(w);
412 coefficients[3] = ratio.multiply(-2).add(fd).add(fdP).divide(w2);
413 polynomials[i] = new FieldPolynomialFunction<>(coefficients);
414 }
415
416 return new FieldPolynomialSplineFunction<>(xvals, polynomials);
417
418 }
419 }