1 files changed, 239 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java
new file mode 100644
index 0000000..15ed322
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java
@@ -0,0 +1,239 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math3.analysis.interpolation;
+
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.List;
+
+import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
+import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableVectorFunction;
+import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
+import org.apache.commons.math3.exception.MathArithmeticException;
+import org.apache.commons.math3.exception.NoDataException;
+import org.apache.commons.math3.exception.ZeroException;
+import org.apache.commons.math3.exception.util.LocalizedFormats;
+import org.apache.commons.math3.util.CombinatoricsUtils;
+
+/** Polynomial interpolator using both sample values and sample derivatives.
+ * <p>
+ * The interpolation polynomials match all sample points, including both values
+ * and provided derivatives. There is one polynomial for each component of
+ * the values vector. All polynomials have the same degree. The degree of the
+ * polynomials depends on the number of points and number of derivatives at each
+ * point. For example the interpolation polynomials for n sample points without
+ * any derivatives all have degree n-1. The interpolation polynomials for n
+ * sample points with the two extreme points having value and first derivative
+ * and the remaining points having value only all have degree n+1. The
+ * interpolation polynomial for n sample points with value, first and second
+ * derivative for all points all have degree 3n-1.
+ * </p>
+ *
+ * @since 3.1
+ */
+public class HermiteInterpolator implements UnivariateDifferentiableVectorFunction {
+
+    /** Sample abscissae. */
+    private final List<Double> abscissae;
+
+    /** Top diagonal of the divided differences array. */
+    private final List<double[]> topDiagonal;
+
+    /** Bottom diagonal of the divided differences array. */
+    private final List<double[]> bottomDiagonal;
+
+    /** Create an empty interpolator.
+     */
+    public HermiteInterpolator() {
+        this.abscissae      = new ArrayList<Double>();
+        this.topDiagonal    = new ArrayList<double[]>();
+        this.bottomDiagonal = new ArrayList<double[]>();
+    }
+
+    /** Add a sample point.
+     * <p>
+     * This method must be called once for each sample point. It is allowed to
+     * mix some calls with values only with calls with values and first
+     * derivatives.
+     * </p>
+     * <p>
+     * The point abscissae for all calls <em>must</em> be different.
+     * </p>
+     * @param x abscissa of the sample point
+     * @param value value and derivatives of the sample point
+     * (if only one row is passed, it is the value, if two rows are
+     * passed the first one is the value and the second the derivative
+     * and so on)
+     * @exception ZeroException if the abscissa difference between added point
+     * and a previous point is zero (i.e. the two points are at same abscissa)
+     * @exception MathArithmeticException if the number of derivatives is larger
+     * than 20, which prevents computation of a factorial
+     */
+    public void addSamplePoint(final double x, final double[] ... value)
+        throws ZeroException, MathArithmeticException {
+
+        for (int i = 0; i < value.length; ++i) {
+
+            final double[] y = value[i].clone();
+            if (i > 1) {
+                double inv = 1.0 / CombinatoricsUtils.factorial(i);
+                for (int j = 0; j < y.length; ++j) {
+                    y[j] *= inv;
+                }
+            }
+
+            // update the bottom diagonal of the divided differences array
+            final int n = abscissae.size();
+            bottomDiagonal.add(n - i, y);
+            double[] bottom0 = y;
+            for (int j = i; j < n; ++j) {
+                final double[] bottom1 = bottomDiagonal.get(n - (j + 1));
+                final double inv = 1.0 / (x - abscissae.get(n - (j + 1)));
+                if (Double.isInfinite(inv)) {
+                    throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x);
+                }
+                for (int k = 0; k < y.length; ++k) {
+                    bottom1[k] = inv * (bottom0[k] - bottom1[k]);
+                }
+                bottom0 = bottom1;
+            }
+
+            // update the top diagonal of the divided differences array
+            topDiagonal.add(bottom0.clone());
+
+            // update the abscissae array
+            abscissae.add(x);
+
+        }
+
+    }
+
+    /** Compute the interpolation polynomials.
+     * @return interpolation polynomials array
+     * @exception NoDataException if sample is empty
+     */
+    public PolynomialFunction[] getPolynomials()
+        throws NoDataException {
+
+        // safety check
+        checkInterpolation();
+
+        // iteration initialization
+        final PolynomialFunction zero = polynomial(0);
+        PolynomialFunction[] polynomials = new PolynomialFunction[topDiagonal.get(0).length];
+        for (int i = 0; i < polynomials.length; ++i) {
+            polynomials[i] = zero;
+        }
+        PolynomialFunction coeff = polynomial(1);
+
+        // build the polynomials by iterating on the top diagonal of the divided differences array
+        for (int i = 0; i < topDiagonal.size(); ++i) {
+            double[] tdi = topDiagonal.get(i);
+            for (int k = 0; k < polynomials.length; ++k) {
+                polynomials[k] = polynomials[k].add(coeff.multiply(polynomial(tdi[k])));
+            }
+            coeff = coeff.multiply(polynomial(-abscissae.get(i), 1.0));
+        }
+
+        return polynomials;
+
+    }
+
+    /** Interpolate value at a specified abscissa.
+     * <p>
+     * Calling this method is equivalent to call the {@link PolynomialFunction#value(double)
+     * value} methods of all polynomials returned by {@link #getPolynomials() getPolynomials},
+     * except it does not build the intermediate polynomials, so this method is faster and
+     * numerically more stable.
+     * </p>
+     * @param x interpolation abscissa
+     * @return interpolated value
+     * @exception NoDataException if sample is empty
+     */
+    public double[] value(double x)
+        throws NoDataException {
+
+        // safety check
+        checkInterpolation();
+
+        final double[] value = new double[topDiagonal.get(0).length];
+        double valueCoeff = 1;
+        for (int i = 0; i < topDiagonal.size(); ++i) {
+            double[] dividedDifference = topDiagonal.get(i);
+            for (int k = 0; k < value.length; ++k) {
+                value[k] += dividedDifference[k] * valueCoeff;
+            }
+            final double deltaX = x - abscissae.get(i);
+            valueCoeff *= deltaX;
+        }
+
+        return value;
+
+    }
+
+    /** Interpolate value at a specified abscissa.
+     * <p>
+     * Calling this method is equivalent to call the {@link
+     * PolynomialFunction#value(DerivativeStructure) value} methods of all polynomials
+     * returned by {@link #getPolynomials() getPolynomials}, except it does not build the
+     * intermediate polynomials, so this method is faster and numerically more stable.
+     * </p>
+     * @param x interpolation abscissa
+     * @return interpolated value
+     * @exception NoDataException if sample is empty
+     */
+    public DerivativeStructure[] value(final DerivativeStructure x)
+        throws NoDataException {
+
+        // safety check
+        checkInterpolation();
+
+        final DerivativeStructure[] value = new DerivativeStructure[topDiagonal.get(0).length];
+        Arrays.fill(value, x.getField().getZero());
+        DerivativeStructure valueCoeff = x.getField().getOne();
+        for (int i = 0; i < topDiagonal.size(); ++i) {
+            double[] dividedDifference = topDiagonal.get(i);
+            for (int k = 0; k < value.length; ++k) {
+                value[k] = value[k].add(valueCoeff.multiply(dividedDifference[k]));
+            }
+            final DerivativeStructure deltaX = x.subtract(abscissae.get(i));
+            valueCoeff = valueCoeff.multiply(deltaX);
+        }
+
+        return value;
+
+    }
+
+    /** Check interpolation can be performed.
+     * @exception NoDataException if interpolation cannot be performed
+     * because sample is empty
+     */
+    private void checkInterpolation() throws NoDataException {
+        if (abscissae.isEmpty()) {
+            throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
+        }
+    }
+
+    /** Create a polynomial from its coefficients.
+     * @param c polynomials coefficients
+     * @return polynomial
+     */
+    private PolynomialFunction polynomial(double ... c) {
+        return new PolynomialFunction(c);
+    }
+
+}