HermiteInterpolator.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math3.analysis.interpolation;

  18. import java.util.ArrayList;
  19. import java.util.Arrays;
  20. import java.util.List;

  21. import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
  22. import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableVectorFunction;
  23. import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
  24. import org.apache.commons.math3.exception.MathArithmeticException;
  25. import org.apache.commons.math3.exception.NoDataException;
  26. import org.apache.commons.math3.exception.ZeroException;
  27. import org.apache.commons.math3.exception.util.LocalizedFormats;
  28. import org.apache.commons.math3.util.CombinatoricsUtils;

  29. /** Polynomial interpolator using both sample values and sample derivatives.
  30.  * <p>
  31.  * The interpolation polynomials match all sample points, including both values
  32.  * and provided derivatives. There is one polynomial for each component of
  33.  * the values vector. All polynomials have the same degree. The degree of the
  34.  * polynomials depends on the number of points and number of derivatives at each
  35.  * point. For example the interpolation polynomials for n sample points without
  36.  * any derivatives all have degree n-1. The interpolation polynomials for n
  37.  * sample points with the two extreme points having value and first derivative
  38.  * and the remaining points having value only all have degree n+1. The
  39.  * interpolation polynomial for n sample points with value, first and second
  40.  * derivative for all points all have degree 3n-1.
  41.  * </p>
  42.  *
  43.  * @since 3.1
  44.  */
  45. public class HermiteInterpolator implements UnivariateDifferentiableVectorFunction {

  46.     /** Sample abscissae. */
  47.     private final List<Double> abscissae;

  48.     /** Top diagonal of the divided differences array. */
  49.     private final List<double[]> topDiagonal;

  50.     /** Bottom diagonal of the divided differences array. */
  51.     private final List<double[]> bottomDiagonal;

  52.     /** Create an empty interpolator.
  53.      */
  54.     public HermiteInterpolator() {
  55.         this.abscissae      = new ArrayList<Double>();
  56.         this.topDiagonal    = new ArrayList<double[]>();
  57.         this.bottomDiagonal = new ArrayList<double[]>();
  58.     }

  59.     /** Add a sample point.
  60.      * <p>
  61.      * This method must be called once for each sample point. It is allowed to
  62.      * mix some calls with values only with calls with values and first
  63.      * derivatives.
  64.      * </p>
  65.      * <p>
  66.      * The point abscissae for all calls <em>must</em> be different.
  67.      * </p>
  68.      * @param x abscissa of the sample point
  69.      * @param value value and derivatives of the sample point
  70.      * (if only one row is passed, it is the value, if two rows are
  71.      * passed the first one is the value and the second the derivative
  72.      * and so on)
  73.      * @exception ZeroException if the abscissa difference between added point
  74.      * and a previous point is zero (i.e. the two points are at same abscissa)
  75.      * @exception MathArithmeticException if the number of derivatives is larger
  76.      * than 20, which prevents computation of a factorial
  77.      */
  78.     public void addSamplePoint(final double x, final double[] ... value)
  79.         throws ZeroException, MathArithmeticException {

  80.         for (int i = 0; i < value.length; ++i) {

  81.             final double[] y = value[i].clone();
  82.             if (i > 1) {
  83.                 double inv = 1.0 / CombinatoricsUtils.factorial(i);
  84.                 for (int j = 0; j < y.length; ++j) {
  85.                     y[j] *= inv;
  86.                 }
  87.             }

  88.             // update the bottom diagonal of the divided differences array
  89.             final int n = abscissae.size();
  90.             bottomDiagonal.add(n - i, y);
  91.             double[] bottom0 = y;
  92.             for (int j = i; j < n; ++j) {
  93.                 final double[] bottom1 = bottomDiagonal.get(n - (j + 1));
  94.                 final double inv = 1.0 / (x - abscissae.get(n - (j + 1)));
  95.                 if (Double.isInfinite(inv)) {
  96.                     throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x);
  97.                 }
  98.                 for (int k = 0; k < y.length; ++k) {
  99.                     bottom1[k] = inv * (bottom0[k] - bottom1[k]);
  100.                 }
  101.                 bottom0 = bottom1;
  102.             }

  103.             // update the top diagonal of the divided differences array
  104.             topDiagonal.add(bottom0.clone());

  105.             // update the abscissae array
  106.             abscissae.add(x);

  107.         }

  108.     }

  109.     /** Compute the interpolation polynomials.
  110.      * @return interpolation polynomials array
  111.      * @exception NoDataException if sample is empty
  112.      */
  113.     public PolynomialFunction[] getPolynomials()
  114.         throws NoDataException {

  115.         // safety check
  116.         checkInterpolation();

  117.         // iteration initialization
  118.         final PolynomialFunction zero = polynomial(0);
  119.         PolynomialFunction[] polynomials = new PolynomialFunction[topDiagonal.get(0).length];
  120.         for (int i = 0; i < polynomials.length; ++i) {
  121.             polynomials[i] = zero;
  122.         }
  123.         PolynomialFunction coeff = polynomial(1);

  124.         // build the polynomials by iterating on the top diagonal of the divided differences array
  125.         for (int i = 0; i < topDiagonal.size(); ++i) {
  126.             double[] tdi = topDiagonal.get(i);
  127.             for (int k = 0; k < polynomials.length; ++k) {
  128.                 polynomials[k] = polynomials[k].add(coeff.multiply(polynomial(tdi[k])));
  129.             }
  130.             coeff = coeff.multiply(polynomial(-abscissae.get(i), 1.0));
  131.         }

  132.         return polynomials;

  133.     }

  134.     /** Interpolate value at a specified abscissa.
  135.      * <p>
  136.      * Calling this method is equivalent to call the {@link PolynomialFunction#value(double)
  137.      * value} methods of all polynomials returned by {@link #getPolynomials() getPolynomials},
  138.      * except it does not build the intermediate polynomials, so this method is faster and
  139.      * numerically more stable.
  140.      * </p>
  141.      * @param x interpolation abscissa
  142.      * @return interpolated value
  143.      * @exception NoDataException if sample is empty
  144.      */
  145.     public double[] value(double x)
  146.         throws NoDataException {

  147.         // safety check
  148.         checkInterpolation();

  149.         final double[] value = new double[topDiagonal.get(0).length];
  150.         double valueCoeff = 1;
  151.         for (int i = 0; i < topDiagonal.size(); ++i) {
  152.             double[] dividedDifference = topDiagonal.get(i);
  153.             for (int k = 0; k < value.length; ++k) {
  154.                 value[k] += dividedDifference[k] * valueCoeff;
  155.             }
  156.             final double deltaX = x - abscissae.get(i);
  157.             valueCoeff *= deltaX;
  158.         }

  159.         return value;

  160.     }

  161.     /** Interpolate value at a specified abscissa.
  162.      * <p>
  163.      * Calling this method is equivalent to call the {@link
  164.      * PolynomialFunction#value(DerivativeStructure) value} methods of all polynomials
  165.      * returned by {@link #getPolynomials() getPolynomials}, except it does not build the
  166.      * intermediate polynomials, so this method is faster and numerically more stable.
  167.      * </p>
  168.      * @param x interpolation abscissa
  169.      * @return interpolated value
  170.      * @exception NoDataException if sample is empty
  171.      */
  172.     public DerivativeStructure[] value(final DerivativeStructure x)
  173.         throws NoDataException {

  174.         // safety check
  175.         checkInterpolation();

  176.         final DerivativeStructure[] value = new DerivativeStructure[topDiagonal.get(0).length];
  177.         Arrays.fill(value, x.getField().getZero());
  178.         DerivativeStructure valueCoeff = x.getField().getOne();
  179.         for (int i = 0; i < topDiagonal.size(); ++i) {
  180.             double[] dividedDifference = topDiagonal.get(i);
  181.             for (int k = 0; k < value.length; ++k) {
  182.                 value[k] = value[k].add(valueCoeff.multiply(dividedDifference[k]));
  183.             }
  184.             final DerivativeStructure deltaX = x.subtract(abscissae.get(i));
  185.             valueCoeff = valueCoeff.multiply(deltaX);
  186.         }

  187.         return value;

  188.     }

  189.     /** Check interpolation can be performed.
  190.      * @exception NoDataException if interpolation cannot be performed
  191.      * because sample is empty
  192.      */
  193.     private void checkInterpolation() throws NoDataException {
  194.         if (abscissae.isEmpty()) {
  195.             throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
  196.         }
  197.     }

  198.     /** Create a polynomial from its coefficients.
  199.      * @param c polynomials coefficients
  200.      * @return polynomial
  201.      */
  202.     private PolynomialFunction polynomial(double ... c) {
  203.         return new PolynomialFunction(c);
  204.     }

  205. }