    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    }