FieldHermiteInterpolator.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.math4.legacy.analysis.interpolation;

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

  22. import org.apache.commons.math4.legacy.core.FieldElement;
  23. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  24. import org.apache.commons.math4.legacy.exception.MathArithmeticException;
  25. import org.apache.commons.math4.legacy.exception.NoDataException;
  26. import org.apache.commons.math4.legacy.exception.NullArgumentException;
  27. import org.apache.commons.math4.legacy.exception.ZeroException;
  28. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  29. import org.apache.commons.math4.legacy.core.MathArrays;

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

  51.     /** Sample abscissae. */
  52.     private final List<T> abscissae;

  54.     /** Top diagonal of the divided differences array. */
  55.     private final List<T[]> topDiagonal;

  57.     /** Bottom diagonal of the divided differences array. */
  58.     private final List<T[]> bottomDiagonal;

  60.     /** Create an empty interpolator.
  61.      */
  62.     public FieldHermiteInterpolator() {
  63.         this.abscissae      = new ArrayList<>();
  64.         this.topDiagonal    = new ArrayList<>();
  65.         this.bottomDiagonal = new ArrayList<>();
  66.     }

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

  94.         NullArgumentException.check(x);
  95.         T factorial = x.getField().getOne();
  96.         for (int i = 0; i < value.length; ++i) {

  98.             final T[] y = value[i].clone();
  99.             if (i > 1) {
  100.                 factorial = factorial.multiply(i);
  101.                 final T inv = factorial.reciprocal();
  102.                 for (int j = 0; j < y.length; ++j) {
  103.                     y[j] = y[j].multiply(inv);
  104.                 }
  105.             }

  107.             // update the bottom diagonal of the divided differences array
  108.             final int n = abscissae.size();
  109.             bottomDiagonal.add(n - i, y);
  110.             T[] bottom0 = y;
  111.             for (int j = i; j < n; ++j) {
  112.                 final T[] bottom1 = bottomDiagonal.get(n - (j + 1));
  113.                 if (x.equals(abscissae.get(n - (j + 1)))) {
  114.                     throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x);
  115.                 }
  116.                 final T inv = x.subtract(abscissae.get(n - (j + 1))).reciprocal();
  117.                 for (int k = 0; k < y.length; ++k) {
  118.                     bottom1[k] = inv.multiply(bottom0[k].subtract(bottom1[k]));
  119.                 }
  120.                 bottom0 = bottom1;
  121.             }

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

  126.             // update the abscissae array
  127.             abscissae.add(x);
  128.         }
  129.     }

  131.     /** Interpolate value at a specified abscissa.
  132.      * @param x interpolation abscissa
  133.      * @return interpolated value
  134.      * @exception NoDataException if sample is empty
  135.      * @throws NullArgumentException if x is null
  136.      */
  137.     public T[] value(T x) throws NoDataException, NullArgumentException {

  139.         // safety check
  140.         NullArgumentException.check(x);
  141.         if (abscissae.isEmpty()) {
  142.             throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
  143.         }

  145.         final T[] value = MathArrays.buildArray(x.getField(), topDiagonal.get(0).length);
  146.         T valueCoeff = x.getField().getOne();
  147.         for (int i = 0; i < topDiagonal.size(); ++i) {
  148.             T[] dividedDifference = topDiagonal.get(i);
  149.             for (int k = 0; k < value.length; ++k) {
  150.                 value[k] = value[k].add(dividedDifference[k].multiply(valueCoeff));
  151.             }
  152.             final T deltaX = x.subtract(abscissae.get(i));
  153.             valueCoeff = valueCoeff.multiply(deltaX);
  154.         }

  156.         return value;
  157.     }

  159.     /** Interpolate value and first derivatives at a specified abscissa.
  160.      * @param x interpolation abscissa
  161.      * @param order maximum derivation order
  162.      * @return interpolated value and derivatives (value in row 0,
  163.      * 1<sup>st</sup> derivative in row 1, ... n<sup>th</sup> derivative in row n)
  164.      * @exception NoDataException if sample is empty
  165.      * @throws NullArgumentException if x is null
  166.      */
  167.     public T[][] derivatives(T x, int order) throws NoDataException, NullArgumentException {

  169.         // safety check
  170.         NullArgumentException.check(x);
  171.         if (abscissae.isEmpty()) {
  172.             throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
  173.         }

  175.         final T zero = x.getField().getZero();
  176.         final T one  = x.getField().getOne();
  177.         final T[] tj = MathArrays.buildArray(x.getField(), order + 1);
  178.         tj[0] = zero;
  179.         for (int i = 0; i < order; ++i) {
  180.             tj[i + 1] = tj[i].add(one);
  181.         }

  183.         final T[][] derivatives =
  184.                 MathArrays.buildArray(x.getField(), order + 1, topDiagonal.get(0).length);
  185.         final T[] valueCoeff = MathArrays.buildArray(x.getField(), order + 1);
  186.         valueCoeff[0] = x.getField().getOne();
  187.         for (int i = 0; i < topDiagonal.size(); ++i) {
  188.             T[] dividedDifference = topDiagonal.get(i);
  189.             final T deltaX = x.subtract(abscissae.get(i));
  190.             for (int j = order; j >= 0; --j) {
  191.                 for (int k = 0; k < derivatives[j].length; ++k) {
  192.                     derivatives[j][k] =
  193.                             derivatives[j][k].add(dividedDifference[k].multiply(valueCoeff[j]));
  194.                 }
  195.                 valueCoeff[j] = valueCoeff[j].multiply(deltaX);
  196.                 if (j > 0) {
  197.                     valueCoeff[j] = valueCoeff[j].add(tj[j].multiply(valueCoeff[j - 1]));
  198.                 }
  199.             }
  200.         }

  202.         return derivatives;
  203.     }
  204. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT