SplineInterpolator.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.lang.reflect.Array;

  26. import org.hipparchus.Field;
  27. import org.hipparchus.CalculusFieldElement;
  28. import org.hipparchus.analysis.polynomials.FieldPolynomialFunction;
  29. import org.hipparchus.analysis.polynomials.FieldPolynomialSplineFunction;
  30. import org.hipparchus.analysis.polynomials.PolynomialFunction;
  31. import org.hipparchus.analysis.polynomials.PolynomialSplineFunction;
  32. import org.hipparchus.exception.LocalizedCoreFormats;
  33. import org.hipparchus.exception.MathIllegalArgumentException;
  34. import org.hipparchus.util.MathArrays;
  35. import org.hipparchus.util.MathUtils;

  37. /**
  38.  * Computes a natural (also known as "free", "unclamped") cubic spline interpolation for the data set.
  39.  * <p>
  40.  * The {@link #interpolate(double[], double[])} method returns a {@link PolynomialSplineFunction}
  41.  * consisting of n cubic polynomials, defined over the subintervals determined by the x values,
  42.  * {@code x[0] < x[i] ... < x[n].}  The x values are referred to as "knot points."</p>
  43.  * <p>
  44.  * The value of the PolynomialSplineFunction at a point x that is greater than or equal to the smallest
  45.  * knot point and strictly less than the largest knot point is computed by finding the subinterval to which
  46.  * x belongs and computing the value of the corresponding polynomial at <code>x - x[i] </code> where
  47.  * <code>i</code> is the index of the subinterval.  See {@link PolynomialSplineFunction} for more details.
  48.  * </p>
  49.  * <p>
  50.  * The interpolating polynomials satisfy:
  51.  * </p>
  52.  * <ol>
  53.  * <li>The value of the PolynomialSplineFunction at each of the input x values equals the
  54.  *  corresponding y value.</li>
  55.  * <li>Adjacent polynomials are equal through two derivatives at the knot points (i.e., adjacent polynomials
  56.  *  "match up" at the knot points, as do their first and second derivatives).</li>
  57.  * </ol>
  58.  * <p>
  59.  * The cubic spline interpolation algorithm implemented is as described in R.L. Burden, J.D. Faires,
  60.  * <u>Numerical Analysis</u>, 4th Ed., 1989, PWS-Kent, ISBN 0-53491-585-X, pp 126-131.
  61.  * </p>
  62.  *
  63.  */
  64. public class SplineInterpolator implements UnivariateInterpolator, FieldUnivariateInterpolator {

  66.     /** Empty constructor.
  67.      * <p>
  68.      * This constructor is not strictly necessary, but it prevents spurious
  69.      * javadoc warnings with JDK 18 and later.
  70.      * </p>
  71.      * @since 3.0
  72.      */
  73.     public SplineInterpolator() { // NOPMD - unnecessary constructor added intentionally to make javadoc happy
  74.         // nothing to do
  75.     }

  77.     /**
  78.      * Computes an interpolating function for the data set.
  79.      * @param x the arguments for the interpolation points
  80.      * @param y the values for the interpolation points
  81.      * @return a function which interpolates the data set
  82.      * @throws MathIllegalArgumentException if {@code x} and {@code y}
  83.      * have different sizes.
  84.      * @throws MathIllegalArgumentException if {@code x} is not sorted in
  85.      * strict increasing order.
  86.      * @throws MathIllegalArgumentException if the size of {@code x} is smaller
  87.      * than 3.
  88.      */
  89.     @Override
  90.     public PolynomialSplineFunction interpolate(double[] x, double[] y)
  91.         throws MathIllegalArgumentException {

  93.         MathUtils.checkNotNull(x);
  94.         MathUtils.checkNotNull(y);
  95.         MathArrays.checkEqualLength(x, y);
  96.         if (x.length < 3) {
  97.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_POINTS,
  98.                                                    x.length, 3, true);
  99.         }

  101.         // Number of intervals.  The number of data points is n + 1.
  102.         final int n = x.length - 1;

  104.         MathArrays.checkOrder(x);

  106.         // Differences between knot points
  107.         final double[] h = new double[n];
  108.         for (int i = 0; i < n; i++) {
  109.             h[i] = x[i + 1] - x[i];
  110.         }

  112.         final double[] mu = new double[n];
  113.         final double[] z = new double[n + 1];
  114.         mu[0] = 0d;
  115.         z[0] = 0d;
  116.         double g;
  117.         for (int i = 1; i < n; i++) {
  118.             g = 2d * (x[i+1]  - x[i - 1]) - h[i - 1] * mu[i -1];
  119.             mu[i] = h[i] / g;
  120.             z[i] = (3d * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1])+ y[i - 1] * h[i]) /
  121.                     (h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g;
  122.         }

  124.         // cubic spline coefficients --  b is linear, c quadratic, d is cubic (original y's are constants)
  125.         final double[] b = new double[n];
  126.         final double[] c = new double[n + 1];
  127.         final double[] d = new double[n];

  129.         z[n] = 0d;
  130.         c[n] = 0d;

  132.         for (int j = n -1; j >=0; j--) {
  133.             c[j] = z[j] - mu[j] * c[j + 1];
  134.             b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2d * c[j]) / 3d;
  135.             d[j] = (c[j + 1] - c[j]) / (3d * h[j]);
  136.         }

  138.         final PolynomialFunction[] polynomials = new PolynomialFunction[n];
  139.         final double[] coefficients = new double[4];
  140.         for (int i = 0; i < n; i++) {
  141.             coefficients[0] = y[i];
  142.             coefficients[1] = b[i];
  143.             coefficients[2] = c[i];
  144.             coefficients[3] = d[i];
  145.             polynomials[i] = new PolynomialFunction(coefficients);
  146.         }

  148.         return new PolynomialSplineFunction(x, polynomials);
  149.     }

  151.     /**
  152.      * Computes an interpolating function for the data set.
  153.      * @param x the arguments for the interpolation points
  154.      * @param y the values for the interpolation points
  155.      * @param <T> the type of the field elements
  156.      * @return a function which interpolates the data set
  157.      * @throws MathIllegalArgumentException if {@code x} and {@code y}
  158.      * have different sizes.
  159.      * @throws MathIllegalArgumentException if {@code x} is not sorted in
  160.      * strict increasing order.
  161.      * @throws MathIllegalArgumentException if the size of {@code x} is smaller
  162.      * than 3.
  163.      * @since 1.5
  164.      */
  165.     @Override
  166.     public <T extends CalculusFieldElement<T>> FieldPolynomialSplineFunction<T> interpolate(
  167.                     T[] x, T[] y)
  168.         throws MathIllegalArgumentException {

  170.         MathUtils.checkNotNull(x);
  171.         MathUtils.checkNotNull(y);
  172.         MathArrays.checkEqualLength(x, y);
  173.         if (x.length < 3) {
  174.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_OF_POINTS,
  175.                                                    x.length, 3, true);
  176.         }

  178.         // Number of intervals.  The number of data points is n + 1.
  179.         final int n = x.length - 1;

  181.         MathArrays.checkOrder(x);

  183.         // Differences between knot points
  184.         final Field<T> field = x[0].getField();
  185.         final T[] h = MathArrays.buildArray(field, n);
  186.         for (int i = 0; i < n; i++) {
  187.             h[i] = x[i + 1].subtract(x[i]);
  188.         }

  190.         final T[] mu = MathArrays.buildArray(field, n);
  191.         final T[] z = MathArrays.buildArray(field, n + 1);
  192.         mu[0] = field.getZero();
  193.         z[0]  = field.getZero();
  194.         for (int i = 1; i < n; i++) {
  195.             final T g = x[i+1].subtract(x[i - 1]).multiply(2).subtract(h[i - 1].multiply(mu[i -1]));
  196.             mu[i] = h[i].divide(g);
  197.             z[i] =          y[i + 1].multiply(h[i - 1]).
  198.                    subtract(y[i].multiply(x[i + 1].subtract(x[i - 1]))).
  199.                         add(y[i - 1].multiply(h[i])).
  200.                    multiply(3).
  201.                    divide(h[i - 1].multiply(h[i])).
  202.                    subtract(h[i - 1].multiply(z[i - 1])).
  203.                    divide(g);
  204.         }

  206.         // cubic spline coefficients --  b is linear, c quadratic, d is cubic (original y's are constants)
  207.         final T[] b = MathArrays.buildArray(field, n);
  208.         final T[] c = MathArrays.buildArray(field, n + 1);
  209.         final T[] d = MathArrays.buildArray(field, n);

  211.         z[n] = field.getZero();
  212.         c[n] = field.getZero();

  214.         for (int j = n -1; j >=0; j--) {
  215.             c[j] = z[j].subtract(mu[j].multiply(c[j + 1]));
  216.             b[j] = y[j + 1].subtract(y[j]).divide(h[j]).
  217.                    subtract(h[j].multiply(c[j + 1].add(c[j]).add(c[j])).divide(3));
  218.             d[j] = c[j + 1].subtract(c[j]).divide(h[j].multiply(3));
  219.         }

  221.         @SuppressWarnings("unchecked")
  222.         final FieldPolynomialFunction<T>[] polynomials =
  223.                         (FieldPolynomialFunction<T>[]) Array.newInstance(FieldPolynomialFunction.class, n);
  224.         final T[] coefficients = MathArrays.buildArray(field, 4);
  225.         for (int i = 0; i < n; i++) {
  226.             coefficients[0] = y[i];
  227.             coefficients[1] = b[i];
  228.             coefficients[2] = c[i];
  229.             coefficients[3] = d[i];
  230.             polynomials[i] = new FieldPolynomialFunction<>(coefficients);
  231.         }

  233.         return new FieldPolynomialSplineFunction<>(x, polynomials);
  234.     }

  236. }
Created with JaCoCo 0.8.12.202403310830