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    }