View Javadoc
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 }