001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.analysis.interpolation;
018
019 import java.util.ArrayList;
020 import java.util.Arrays;
021 import java.util.List;
022
023 import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
024 import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableVectorFunction;
025 import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
026 import org.apache.commons.math3.exception.MathArithmeticException;
027 import org.apache.commons.math3.exception.NoDataException;
028 import org.apache.commons.math3.exception.ZeroException;
029 import org.apache.commons.math3.exception.util.LocalizedFormats;
030 import org.apache.commons.math3.util.ArithmeticUtils;
031
032 /** Polynomial interpolator using both sample values and sample derivatives.
033 * <p>
034 * The interpolation polynomials match all sample points, including both values
035 * and provided derivatives. There is one polynomial for each component of
036 * the values vector. All polynomials have the same degree. The degree of the
037 * polynomials depends on the number of points and number of derivatives at each
038 * point. For example the interpolation polynomials for n sample points without
039 * any derivatives all have degree n-1. The interpolation polynomials for n
040 * sample points with the two extreme points having value and first derivative
041 * and the remaining points having value only all have degree n+1. The
042 * interpolation polynomial for n sample points with value, first and second
043 * derivative for all points all have degree 3n-1.
044 * </p>
045 *
046 * @version $Id: HermiteInterpolator.java 1410460 2012-11-16 16:49:38Z erans $
047 * @since 3.1
048 */
049 public class HermiteInterpolator implements UnivariateDifferentiableVectorFunction {
050
051 /** Sample abscissae. */
052 private final List<Double> abscissae;
053
054 /** Top diagonal of the divided differences array. */
055 private final List<double[]> topDiagonal;
056
057 /** Bottom diagonal of the divided differences array. */
058 private final List<double[]> bottomDiagonal;
059
060 /** Create an empty interpolator.
061 */
062 public HermiteInterpolator() {
063 this.abscissae = new ArrayList<Double>();
064 this.topDiagonal = new ArrayList<double[]>();
065 this.bottomDiagonal = new ArrayList<double[]>();
066 }
067
068 /** Add a sample point.
069 * <p>
070 * This method must be called once for each sample point. It is allowed to
071 * mix some calls with values only with calls with values and first
072 * derivatives.
073 * </p>
074 * <p>
075 * The point abscissae for all calls <em>must</em> be different.
076 * </p>
077 * @param x abscissa of the sample point
078 * @param value value and derivatives of the sample point
079 * (if only one row is passed, it is the value, if two rows are
080 * passed the first one is the value and the second the derivative
081 * and so on)
082 * @exception ZeroException if the abscissa difference between added point
083 * and a previous point is zero (i.e. the two points are at same abscissa)
084 * @exception MathArithmeticException if the number of derivatives is larger
085 * than 20, which prevents computation of a factorial
086 */
087 public void addSamplePoint(final double x, final double[] ... value)
088 throws ZeroException, MathArithmeticException {
089
090 for (int i = 0; i < value.length; ++i) {
091
092 final double[] y = value[i].clone();
093 if (i > 1) {
094 double inv = 1.0 / ArithmeticUtils.factorial(i);
095 for (int j = 0; j < y.length; ++j) {
096 y[j] *= inv;
097 }
098 }
099
100 // update the bottom diagonal of the divided differences array
101 final int n = abscissae.size();
102 bottomDiagonal.add(n - i, y);
103 double[] bottom0 = y;
104 for (int j = i; j < n; ++j) {
105 final double[] bottom1 = bottomDiagonal.get(n - (j + 1));
106 final double inv = 1.0 / (x - abscissae.get(n - (j + 1)));
107 if (Double.isInfinite(inv)) {
108 throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x);
109 }
110 for (int k = 0; k < y.length; ++k) {
111 bottom1[k] = inv * (bottom0[k] - bottom1[k]);
112 }
113 bottom0 = bottom1;
114 }
115
116 // update the top diagonal of the divided differences array
117 topDiagonal.add(bottom0.clone());
118
119 // update the abscissae array
120 abscissae.add(x);
121
122 }
123
124 }
125
126 /** Compute the interpolation polynomials.
127 * @return interpolation polynomials array
128 * @exception NoDataException if sample is empty
129 */
130 public PolynomialFunction[] getPolynomials()
131 throws NoDataException {
132
133 // safety check
134 checkInterpolation();
135
136 // iteration initialization
137 final PolynomialFunction zero = polynomial(0);
138 PolynomialFunction[] polynomials = new PolynomialFunction[topDiagonal.get(0).length];
139 for (int i = 0; i < polynomials.length; ++i) {
140 polynomials[i] = zero;
141 }
142 PolynomialFunction coeff = polynomial(1);
143
144 // build the polynomials by iterating on the top diagonal of the divided differences array
145 for (int i = 0; i < topDiagonal.size(); ++i) {
146 double[] tdi = topDiagonal.get(i);
147 for (int k = 0; k < polynomials.length; ++k) {
148 polynomials[k] = polynomials[k].add(coeff.multiply(polynomial(tdi[k])));
149 }
150 coeff = coeff.multiply(polynomial(-abscissae.get(i), 1.0));
151 }
152
153 return polynomials;
154
155 }
156
157 /** Interpolate value at a specified abscissa.
158 * <p>
159 * Calling this method is equivalent to call the {@link PolynomialFunction#value(double)
160 * value} methods of all polynomials returned by {@link #getPolynomials() getPolynomials},
161 * except it does not build the intermediate polynomials, so this method is faster and
162 * numerically more stable.
163 * </p>
164 * @param x interpolation abscissa
165 * @return interpolated value
166 * @exception NoDataException if sample is empty
167 */
168 public double[] value(double x)
169 throws NoDataException {
170
171 // safety check
172 checkInterpolation();
173
174 final double[] value = new double[topDiagonal.get(0).length];
175 double valueCoeff = 1;
176 for (int i = 0; i < topDiagonal.size(); ++i) {
177 double[] dividedDifference = topDiagonal.get(i);
178 for (int k = 0; k < value.length; ++k) {
179 value[k] += dividedDifference[k] * valueCoeff;
180 }
181 final double deltaX = x - abscissae.get(i);
182 valueCoeff *= deltaX;
183 }
184
185 return value;
186
187 }
188
189 /** Interpolate value at a specified abscissa.
190 * <p>
191 * Calling this method is equivalent to call the {@link
192 * PolynomialFunction#value(DerivativeStructure) value} methods of all polynomials
193 * returned by {@link #getPolynomials() getPolynomials}, except it does not build the
194 * intermediate polynomials, so this method is faster and numerically more stable.
195 * </p>
196 * @param x interpolation abscissa
197 * @return interpolated value
198 * @exception NoDataException if sample is empty
199 */
200 public DerivativeStructure[] value(final DerivativeStructure x)
201 throws NoDataException {
202
203 // safety check
204 checkInterpolation();
205
206 final DerivativeStructure[] value = new DerivativeStructure[topDiagonal.get(0).length];
207 Arrays.fill(value, x.getField().getZero());
208 DerivativeStructure valueCoeff = x.getField().getOne();
209 for (int i = 0; i < topDiagonal.size(); ++i) {
210 double[] dividedDifference = topDiagonal.get(i);
211 for (int k = 0; k < value.length; ++k) {
212 value[k] = value[k].add(valueCoeff.multiply(dividedDifference[k]));
213 }
214 final DerivativeStructure deltaX = x.subtract(abscissae.get(i));
215 valueCoeff = valueCoeff.multiply(deltaX);
216 }
217
218 return value;
219
220 }
221
222 /** Check interpolation can be performed.
223 * @exception NoDataException if interpolation cannot be performed
224 * because sample is empty
225 */
226 private void checkInterpolation() throws NoDataException {
227 if (abscissae.isEmpty()) {
228 throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
229 }
230 }
231
232 /** Create a polynomial from its coefficients.
233 * @param c polynomials coefficients
234 * @return polynomial
235 */
236 private PolynomialFunction polynomial(double ... c) {
237 return new PolynomialFunction(c);
238 }
239
240 }