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.  *      https://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.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.analysis.interpolation;

  24. import java.util.ArrayList;
  25. import java.util.Arrays;
  26. import java.util.List;

  28. import org.hipparchus.analysis.differentiation.Derivative;
  29. import org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction;
  30. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  31. import org.hipparchus.exception.LocalizedCoreFormats;
  32. import org.hipparchus.exception.MathIllegalArgumentException;
  33. import org.hipparchus.exception.MathRuntimeException;
  34. import org.hipparchus.exception.NullArgumentException;
  35. import org.hipparchus.util.CombinatoricsUtils;
  36. import org.hipparchus.util.MathArrays;
  37. import org.hipparchus.util.MathUtils;

  39. /** Polynomial interpolator using both sample values and sample derivatives.
  40.  * <p>
  41.  * The interpolation polynomials match all sample points, including both values
  42.  * and provided derivatives. There is one polynomial for each component of
  43.  * the values vector. All polynomials have the same degree. The degree of the
  44.  * polynomials depends on the number of points and number of derivatives at each
  45.  * point. For example the interpolation polynomials for n sample points without
  46.  * any derivatives all have degree n-1. The interpolation polynomials for n
  47.  * sample points with the two extreme points having value and first derivative
  48.  * and the remaining points having value only all have degree n+1. The
  49.  * interpolation polynomial for n sample points with value, first and second
  50.  * derivative for all points all have degree 3n-1.
  51.  * </p>
  52.  *
  53.  */
  54. public class HermiteInterpolator implements UnivariateDifferentiableVectorFunction {

  56.     /** Sample abscissae. */
  57.     private final List<Double> abscissae;

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

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

  65.     /** Create an empty interpolator.
  66.      */
  67.     public HermiteInterpolator() {
  68.         this.abscissae      = new ArrayList<>();
  69.         this.topDiagonal    = new ArrayList<>();
  70.         this.bottomDiagonal = new ArrayList<>();
  71.     }

  73.     /** Add a sample point.
  74.      * <p>
  75.      * This method must be called once for each sample point. It is allowed to
  76.      * mix some calls with values only with calls with values and first
  77.      * derivatives.
  78.      * </p>
  79.      * <p>
  80.      * The point abscissae for all calls <em>must</em> be different.
  81.      * </p>
  82.      * @param x abscissa of the sample point
  83.      * @param value value and derivatives of the sample point
  84.      * (if only one row is passed, it is the value, if two rows are
  85.      * passed the first one is the value and the second the derivative
  86.      * and so on)
  87.      * @exception MathIllegalArgumentException if the abscissa difference between added point
  88.      * and a previous point is zero (i.e. the two points are at same abscissa)
  89.      * @exception MathRuntimeException if the number of derivatives is larger
  90.      * than 20, which prevents computation of a factorial
  91.      */
  92.     public void addSamplePoint(final double x, final double[] ... value)
  93.         throws MathRuntimeException {

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

  97.             final double[] y = value[i].clone();
  98.             if (i > 1) {
  99.                 double inv = 1.0 / CombinatoricsUtils.factorial(i);
  100.                 for (int j = 0; j < y.length; ++j) {
  101.                     y[j] *= inv;
  102.                 }
  103.             }

  105.             // update the bottom diagonal of the divided differences array
  106.             final int n = abscissae.size();
  107.             bottomDiagonal.add(n - i, y);
  108.             double[] bottom0 = y;
  109.             for (int j = i; j < n; ++j) {
  110.                 final double[] bottom1 = bottomDiagonal.get(n - (j + 1));
  111.                 final double inv = 1.0 / (x - abscissae.get(n - (j + 1)));
  112.                 if (Double.isInfinite(inv)) {
  113.                     throw new MathIllegalArgumentException(LocalizedCoreFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x);
  114.                 }
  115.                 for (int k = 0; k < y.length; ++k) {
  116.                     bottom1[k] = inv * (bottom0[k] - bottom1[k]);
  117.                 }
  118.                 bottom0 = bottom1;
  119.             }

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

  124.             // update the abscissae array
  125.             abscissae.add(x);

  127.         }

  129.     }

  131.     /** Compute the interpolation polynomials.
  132.      * @return interpolation polynomials array
  133.      * @exception MathIllegalArgumentException if sample is empty
  134.      */
  135.     public PolynomialFunction[] getPolynomials()
  136.         throws MathIllegalArgumentException {

  138.         // safety check
  139.         checkInterpolation();

  141.         // iteration initialization
  142.         final PolynomialFunction zero = polynomial(0);
  143.         PolynomialFunction[] polynomials = new PolynomialFunction[topDiagonal.get(0).length];
  144.         for (int i = 0; i < polynomials.length; ++i) {
  145.             polynomials[i] = zero;
  146.         }
  147.         PolynomialFunction coeff = polynomial(1);

  149.         // build the polynomials by iterating on the top diagonal of the divided differences array
  150.         for (int i = 0; i < topDiagonal.size(); ++i) {
  151.             double[] tdi = topDiagonal.get(i);
  152.             for (int k = 0; k < polynomials.length; ++k) {
  153.                 polynomials[k] = polynomials[k].add(coeff.multiply(polynomial(tdi[k])));
  154.             }
  155.             coeff = coeff.multiply(polynomial(-abscissae.get(i), 1.0));
  156.         }

  158.         return polynomials;

  160.     }

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

  176.         // safety check
  177.         checkInterpolation();

  179.         final double[] value = new double[topDiagonal.get(0).length];
  180.         double valueCoeff = 1;
  181.         for (int i = 0; i < topDiagonal.size(); ++i) {
  182.             double[] dividedDifference = topDiagonal.get(i);
  183.             for (int k = 0; k < value.length; ++k) {
  184.                 value[k] += dividedDifference[k] * valueCoeff;
  185.             }
  186.             final double deltaX = x - abscissae.get(i);
  187.             valueCoeff *= deltaX;
  188.         }

  190.         return value;

  192.     }

  194.     /** {@inheritDoc}. */
  195.     @Override
  196.     public <T extends Derivative<T>> T[] value(T x)
  197.         throws MathIllegalArgumentException {

  199.         // safety check
  200.         checkInterpolation();

  202.         final T[] value = MathArrays.buildArray(x.getField(), topDiagonal.get(0).length);
  203.         Arrays.fill(value, x.getField().getZero());
  204.         T valueCoeff = x.getField().getOne();
  205.         for (int i = 0; i < topDiagonal.size(); ++i) {
  206.             double[] dividedDifference = topDiagonal.get(i);
  207.             for (int k = 0; k < value.length; ++k) {
  208.                 value[k] = value[k].add(valueCoeff.multiply(dividedDifference[k]));
  209.             }
  210.             final T deltaX = x.subtract(abscissae.get(i));
  211.             valueCoeff = valueCoeff.multiply(deltaX);
  212.         }

  214.         return value;

  216.     }


  219.     /** Interpolate value and first derivatives at a specified abscissa.
  220.      * @param x interpolation abscissa
  221.      * @param order maximum derivation order
  222.      * @return interpolated value and derivatives (value in row 0,
  223.      * 1<sup>st</sup> derivative in row 1, ... n<sup>th</sup> derivative in row n)
  224.      * @exception MathIllegalArgumentException if sample is empty
  225.      * @throws NullArgumentException if x is null
  226.      */
  227.     public double[][] derivatives(double x, int order)
  228.         throws MathIllegalArgumentException, NullArgumentException {

  230.         // safety check
  231.         MathUtils.checkNotNull(x);
  232.         if (abscissae.isEmpty()) {
  233.             throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_INTERPOLATION_SAMPLE);
  234.         }

  236.         final double[] tj = new double[order + 1];
  237.         tj[0] = 0.0;
  238.         for (int i = 0; i < order; ++i) {
  239.             tj[i + 1] = tj[i] + 1;
  240.         }

  242.         final double[][] derivatives = new double[order + 1][topDiagonal.get(0).length];
  243.         final double[] valueCoeff = new double[order + 1];
  244.         valueCoeff[0] = 1.0;
  245.         for (int i = 0; i < topDiagonal.size(); ++i) {
  246.             double[] dividedDifference = topDiagonal.get(i);
  247.             final double deltaX = x - abscissae.get(i);
  248.             for (int j = order; j >= 0; --j) {
  249.                 for (int k = 0; k < derivatives[j].length; ++k) {
  250.                     derivatives[j][k] += dividedDifference[k] * valueCoeff[j];
  251.                 }
  252.                 valueCoeff[j] *= deltaX;
  253.                 if (j > 0) {
  254.                     valueCoeff[j] += tj[j] * valueCoeff[j - 1];
  255.                 }
  256.             }
  257.         }

  259.         return derivatives;

  261.     }

  263.     /** Check interpolation can be performed.
  264.      * @exception MathIllegalArgumentException if interpolation cannot be performed
  265.      * because sample is empty
  266.      */
  267.     private void checkInterpolation() throws MathIllegalArgumentException {
  268.         if (abscissae.isEmpty()) {
  269.             throw new MathIllegalArgumentException(LocalizedCoreFormats.EMPTY_INTERPOLATION_SAMPLE);
  270.         }
  271.     }

  273.     /** Create a polynomial from its coefficients.
  274.      * @param c polynomials coefficients
  275.      * @return polynomial
  276.      */
  277.     private PolynomialFunction polynomial(double ... c) {
  278.         return new PolynomialFunction(c);
  279.     }

  281. }
Created with JaCoCo 0.8.12.202403310830