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.util.ArrayList;
25  import java.util.Arrays;
26  import java.util.List;
27  
28  import org.hipparchus.analysis.differentiation.Derivative;
29  import org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction;
30  import org.hipparchus.analysis.polynomials.PolynomialFunction;
31  import org.hipparchus.exception.LocalizedCoreFormats;
32  import org.hipparchus.exception.MathIllegalArgumentException;
33  import org.hipparchus.exception.MathRuntimeException;
34  import org.hipparchus.exception.NullArgumentException;
35  import org.hipparchus.util.CombinatoricsUtils;
36  import org.hipparchus.util.MathArrays;
37  import org.hipparchus.util.MathUtils;
38  
39  /** Polynomial interpolator using both sample values and sample derivatives.
40   * <p>
41   * The interpolation polynomials match all sample points, including both values
42   * and provided derivatives. There is one polynomial for each component of
43   * the values vector. All polynomials have the same degree. The degree of the
44   * polynomials depends on the number of points and number of derivatives at each
45   * point. For example the interpolation polynomials for n sample points without
46   * any derivatives all have degree n-1. The interpolation polynomials for n
47   * sample points with the two extreme points having value and first derivative
48   * and the remaining points having value only all have degree n+1. The
49   * interpolation polynomial for n sample points with value, first and second
50   * derivative for all points all have degree 3n-1.
51   * </p>
52   *
53   */
54  public class HermiteInterpolator implements UnivariateDifferentiableVectorFunction {
55  
56      /** Sample abscissae. */
57      private final List<Double> abscissae;
58  
59      /** Top diagonal of the divided differences array. */
60      private final List<double[]> topDiagonal;
61  
62      /** Bottom diagonal of the divided differences array. */
63      private final List<double[]> bottomDiagonal;
64  
65      /** Create an empty interpolator.
66       */
67      public HermiteInterpolator() {
68          this.abscissae      = new ArrayList<>();
69          this.topDiagonal    = new ArrayList<>();
70          this.bottomDiagonal = new ArrayList<>();
71      }
72  
73      /** Add a sample point.
74       * <p>
75       * This method must be called once for each sample point. It is allowed to
76       * mix some calls with values only with calls with values and first
77       * derivatives.
78       * </p>
79       * <p>
80       * The point abscissae for all calls <em>must</em> be different.
81       * </p>
82       * @param x abscissa of the sample point
83       * @param value value and derivatives of the sample point
84       * (if only one row is passed, it is the value, if two rows are
85       * passed the first one is the value and the second the derivative
86       * and so on)
87       * @exception MathIllegalArgumentException if the abscissa difference between added point
88       * and a previous point is zero (i.e. the two points are at same abscissa)
89       * @exception MathRuntimeException if the number of derivatives is larger
90       * than 20, which prevents computation of a factorial
91       */
92      public void addSamplePoint(final double x, final double[] ... value)
93          throws MathIllegalArgumentException, MathRuntimeException {
94  
95          for (int i = 0; i < value.length; ++i) {
96  
97              final double[] y = value[i].clone();
98              if (i > 1) {
99                  double inv = 1.0 / CombinatoricsUtils.factorial(i);
100                 for (int j = 0; j < y.length; ++j) {
101                     y[j] *= inv;
102                 }
103             }
104 
105             // update the bottom diagonal of the divided differences array
106             final int n = abscissae.size();
107             bottomDiagonal.add(n - i, y);
108             double[] bottom0 = y;
109             for (int j = i; j < n; ++j) {
110                 final double[] bottom1 = bottomDiagonal.get(n - (j + 1));
111                 final double inv = 1.0 / (x - abscissae.get(n - (j + 1)));
112                 if (Double.isInfinite(inv)) {
113                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x);
114                 }
115                 for (int k = 0; k < y.length; ++k) {
116                     bottom1[k] = inv * (bottom0[k] - bottom1[k]);
117                 }
118                 bottom0 = bottom1;
119             }
120 
121             // update the top diagonal of the divided differences array
122             topDiagonal.add(bottom0.clone());
123 
124             // update the abscissae array
125             abscissae.add(x);
126 
127         }
128 
129     }
130 
131     /** Compute the interpolation polynomials.
132      * @return interpolation polynomials array
133      * @exception MathIllegalArgumentException if sample is empty
134      */
135     public PolynomialFunction[] getPolynomials()
136         throws MathIllegalArgumentException {
137 
138         // safety check
139         checkInterpolation();
140 
141         // iteration initialization
142         final PolynomialFunction zero = polynomial(0);
143         PolynomialFunction[] polynomials = new PolynomialFunction[topDiagonal.get(0).length];
144         for (int i = 0; i < polynomials.length; ++i) {
145             polynomials[i] = zero;
146         }
147         PolynomialFunction coeff = polynomial(1);
148 
149         // build the polynomials by iterating on the top diagonal of the divided differences array
150         for (int i = 0; i < topDiagonal.size(); ++i) {
151             double[] tdi = topDiagonal.get(i);
152             for (int k = 0; k < polynomials.length; ++k) {
153                 polynomials[k] = polynomials[k].add(coeff.multiply(polynomial(tdi[k])));
154             }
155             coeff = coeff.multiply(polynomial(-abscissae.get(i), 1.0));
156         }
157 
158         return polynomials;
159 
160     }
161 
162     /** Interpolate value at a specified abscissa.
163      * <p>
164      * Calling this method is equivalent to call the {@link PolynomialFunction#value(double)
165      * value} methods of all polynomials returned by {@link #getPolynomials() getPolynomials},
166      * except it does not build the intermediate polynomials, so this method is faster and
167      * numerically more stable.
168      * </p>
169      * @param x interpolation abscissa
170      * @return interpolated value
171      * @exception MathIllegalArgumentException if sample is empty
172      */
173     @Override
174     public double[] value(double x) throws MathIllegalArgumentException {
175 
176         // safety check
177         checkInterpolation();
178 
179         final double[] value = new double[topDiagonal.get(0).length];
180         double valueCoeff = 1;
181         for (int i = 0; i < topDiagonal.size(); ++i) {
182             double[] dividedDifference = topDiagonal.get(i);
183             for (int k = 0; k < value.length; ++k) {
184                 value[k] += dividedDifference[k] * valueCoeff;
185             }
186             final double deltaX = x - abscissae.get(i);
187             valueCoeff *= deltaX;
188         }
189 
190         return value;
191 
192     }
193 
194     /** {@inheritDoc}. */
195     @Override
196     public <T extends Derivative<T>> T[] value(T x)
197         throws MathIllegalArgumentException {
198 
199         // safety check
200         checkInterpolation();
201 
202         final T[] value = MathArrays.buildArray(x.getField(), topDiagonal.get(0).length);
203         Arrays.fill(value, x.getField().getZero());
204         T valueCoeff = x.getField().getOne();
205         for (int i = 0; i < topDiagonal.size(); ++i) {
206             double[] dividedDifference = topDiagonal.get(i);
207             for (int k = 0; k < value.length; ++k) {
208                 value[k] = value[k].add(valueCoeff.multiply(dividedDifference[k]));
209             }
210             final T deltaX = x.subtract(abscissae.get(i));
211             valueCoeff = valueCoeff.multiply(deltaX);
212         }
213 
214         return value;
215 
216     }
217 
218 
219     /** Interpolate value and first derivatives at a specified abscissa.
220      * @param x interpolation abscissa
221      * @param order maximum derivation order
222      * @return interpolated value and derivatives (value in row 0,
223      * 1<sup>st</sup> derivative in row 1, ... n<sup>th</sup> derivative in row n)
224      * @exception MathIllegalArgumentException if sample is empty
225      * @throws NullArgumentException if x is null
226      */
227     public double[][] derivatives(double x, int order)
228         throws MathIllegalArgumentException, NullArgumentException {
229 
230         // safety check
231         MathUtils.checkNotNull(x);
232         if (abscissae.isEmpty()) {
233             throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_INTERPOLATION_SAMPLE);
234         }
235 
236         final double[] tj = new double[order + 1];
237         tj[0] = 0.0;
238         for (int i = 0; i < order; ++i) {
239             tj[i + 1] = tj[i] + 1;
240         }
241 
242         final double[][] derivatives = new double[order + 1][topDiagonal.get(0).length];
243         final double[] valueCoeff = new double[order + 1];
244         valueCoeff[0] = 1.0;
245         for (int i = 0; i < topDiagonal.size(); ++i) {
246             double[] dividedDifference = topDiagonal.get(i);
247             final double deltaX = x - abscissae.get(i);
248             for (int j = order; j >= 0; --j) {
249                 for (int k = 0; k < derivatives[j].length; ++k) {
250                     derivatives[j][k] += dividedDifference[k] * valueCoeff[j];
251                 }
252                 valueCoeff[j] *= deltaX;
253                 if (j > 0) {
254                     valueCoeff[j] += tj[j] * valueCoeff[j - 1];
255                 }
256             }
257         }
258 
259         return derivatives;
260 
261     }
262 
263     /** Check interpolation can be performed.
264      * @exception MathIllegalArgumentException if interpolation cannot be performed
265      * because sample is empty
266      */
267     private void checkInterpolation() throws MathIllegalArgumentException {
268         if (abscissae.isEmpty()) {
269             throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_INTERPOLATION_SAMPLE);
270         }
271     }
272 
273     /** Create a polynomial from its coefficients.
274      * @param c polynomials coefficients
275      * @return polynomial
276      */
277     private PolynomialFunction polynomial(double ... c) {
278         return new PolynomialFunction(c);
279     }
280 
281 }