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 */
017package org.apache.commons.math3.analysis.interpolation;
018
019import java.util.ArrayList;
020import java.util.Arrays;
021import java.util.List;
022
023import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
024import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableVectorFunction;
025import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
026import org.apache.commons.math3.exception.MathArithmeticException;
027import org.apache.commons.math3.exception.NoDataException;
028import org.apache.commons.math3.exception.ZeroException;
029import org.apache.commons.math3.exception.util.LocalizedFormats;
030import org.apache.commons.math3.util.CombinatoricsUtils;
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 * @since 3.1
047 */
048public class HermiteInterpolator implements UnivariateDifferentiableVectorFunction {
049
050    /** Sample abscissae. */
051    private final List<Double> abscissae;
052
053    /** Top diagonal of the divided differences array. */
054    private final List<double[]> topDiagonal;
055
056    /** Bottom diagonal of the divided differences array. */
057    private final List<double[]> bottomDiagonal;
058
059    /** Create an empty interpolator.
060     */
061    public HermiteInterpolator() {
062        this.abscissae      = new ArrayList<Double>();
063        this.topDiagonal    = new ArrayList<double[]>();
064        this.bottomDiagonal = new ArrayList<double[]>();
065    }
066
067    /** Add a sample point.
068     * <p>
069     * This method must be called once for each sample point. It is allowed to
070     * mix some calls with values only with calls with values and first
071     * derivatives.
072     * </p>
073     * <p>
074     * The point abscissae for all calls <em>must</em> be different.
075     * </p>
076     * @param x abscissa of the sample point
077     * @param value value and derivatives of the sample point
078     * (if only one row is passed, it is the value, if two rows are
079     * passed the first one is the value and the second the derivative
080     * and so on)
081     * @exception ZeroException if the abscissa difference between added point
082     * and a previous point is zero (i.e. the two points are at same abscissa)
083     * @exception MathArithmeticException if the number of derivatives is larger
084     * than 20, which prevents computation of a factorial
085     */
086    public void addSamplePoint(final double x, final double[] ... value)
087        throws ZeroException, MathArithmeticException {
088
089        for (int i = 0; i < value.length; ++i) {
090
091            final double[] y = value[i].clone();
092            if (i > 1) {
093                double inv = 1.0 / CombinatoricsUtils.factorial(i);
094                for (int j = 0; j < y.length; ++j) {
095                    y[j] *= inv;
096                }
097            }
098
099            // update the bottom diagonal of the divided differences array
100            final int n = abscissae.size();
101            bottomDiagonal.add(n - i, y);
102            double[] bottom0 = y;
103            for (int j = i; j < n; ++j) {
104                final double[] bottom1 = bottomDiagonal.get(n - (j + 1));
105                final double inv = 1.0 / (x - abscissae.get(n - (j + 1)));
106                if (Double.isInfinite(inv)) {
107                    throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x);
108                }
109                for (int k = 0; k < y.length; ++k) {
110                    bottom1[k] = inv * (bottom0[k] - bottom1[k]);
111                }
112                bottom0 = bottom1;
113            }
114
115            // update the top diagonal of the divided differences array
116            topDiagonal.add(bottom0.clone());
117
118            // update the abscissae array
119            abscissae.add(x);
120
121        }
122
123    }
124
125    /** Compute the interpolation polynomials.
126     * @return interpolation polynomials array
127     * @exception NoDataException if sample is empty
128     */
129    public PolynomialFunction[] getPolynomials()
130        throws NoDataException {
131
132        // safety check
133        checkInterpolation();
134
135        // iteration initialization
136        final PolynomialFunction zero = polynomial(0);
137        PolynomialFunction[] polynomials = new PolynomialFunction[topDiagonal.get(0).length];
138        for (int i = 0; i < polynomials.length; ++i) {
139            polynomials[i] = zero;
140        }
141        PolynomialFunction coeff = polynomial(1);
142
143        // build the polynomials by iterating on the top diagonal of the divided differences array
144        for (int i = 0; i < topDiagonal.size(); ++i) {
145            double[] tdi = topDiagonal.get(i);
146            for (int k = 0; k < polynomials.length; ++k) {
147                polynomials[k] = polynomials[k].add(coeff.multiply(polynomial(tdi[k])));
148            }
149            coeff = coeff.multiply(polynomial(-abscissae.get(i), 1.0));
150        }
151
152        return polynomials;
153
154    }
155
156    /** Interpolate value at a specified abscissa.
157     * <p>
158     * Calling this method is equivalent to call the {@link PolynomialFunction#value(double)
159     * value} methods of all polynomials returned by {@link #getPolynomials() getPolynomials},
160     * except it does not build the intermediate polynomials, so this method is faster and
161     * numerically more stable.
162     * </p>
163     * @param x interpolation abscissa
164     * @return interpolated value
165     * @exception NoDataException if sample is empty
166     */
167    public double[] value(double x)
168        throws NoDataException {
169
170        // safety check
171        checkInterpolation();
172
173        final double[] value = new double[topDiagonal.get(0).length];
174        double valueCoeff = 1;
175        for (int i = 0; i < topDiagonal.size(); ++i) {
176            double[] dividedDifference = topDiagonal.get(i);
177            for (int k = 0; k < value.length; ++k) {
178                value[k] += dividedDifference[k] * valueCoeff;
179            }
180            final double deltaX = x - abscissae.get(i);
181            valueCoeff *= deltaX;
182        }
183
184        return value;
185
186    }
187
188    /** Interpolate value at a specified abscissa.
189     * <p>
190     * Calling this method is equivalent to call the {@link
191     * PolynomialFunction#value(DerivativeStructure) value} methods of all polynomials
192     * returned by {@link #getPolynomials() getPolynomials}, except it does not build the
193     * intermediate polynomials, so this method is faster and numerically more stable.
194     * </p>
195     * @param x interpolation abscissa
196     * @return interpolated value
197     * @exception NoDataException if sample is empty
198     */
199    public DerivativeStructure[] value(final DerivativeStructure x)
200        throws NoDataException {
201
202        // safety check
203        checkInterpolation();
204
205        final DerivativeStructure[] value = new DerivativeStructure[topDiagonal.get(0).length];
206        Arrays.fill(value, x.getField().getZero());
207        DerivativeStructure valueCoeff = x.getField().getOne();
208        for (int i = 0; i < topDiagonal.size(); ++i) {
209            double[] dividedDifference = topDiagonal.get(i);
210            for (int k = 0; k < value.length; ++k) {
211                value[k] = value[k].add(valueCoeff.multiply(dividedDifference[k]));
212            }
213            final DerivativeStructure deltaX = x.subtract(abscissae.get(i));
214            valueCoeff = valueCoeff.multiply(deltaX);
215        }
216
217        return value;
218
219    }
220
221    /** Check interpolation can be performed.
222     * @exception NoDataException if interpolation cannot be performed
223     * because sample is empty
224     */
225    private void checkInterpolation() throws NoDataException {
226        if (abscissae.isEmpty()) {
227            throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
228        }
229    }
230
231    /** Create a polynomial from its coefficients.
232     * @param c polynomials coefficients
233     * @return polynomial
234     */
235    private PolynomialFunction polynomial(double ... c) {
236        return new PolynomialFunction(c);
237    }
238
239}