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 }