| /* |
| * 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.math4.legacy.analysis.interpolation; |
| |
| import java.util.ArrayList; |
| import java.util.Arrays; |
| import java.util.List; |
| |
| import org.apache.commons.math4.legacy.analysis.differentiation.DerivativeStructure; |
| import org.apache.commons.math4.legacy.analysis.differentiation.UnivariateDifferentiableVectorFunction; |
| import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunction; |
| import org.apache.commons.math4.legacy.exception.MathArithmeticException; |
| import org.apache.commons.math4.legacy.exception.NoDataException; |
| import org.apache.commons.math4.legacy.exception.ZeroException; |
| import org.apache.commons.math4.legacy.exception.util.LocalizedFormats; |
| import org.apache.commons.numbers.combinatorics.Factorial; |
| |
| /** 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<>(); |
| this.topDiagonal = new ArrayList<>(); |
| this.bottomDiagonal = new ArrayList<>(); |
| } |
| |
| /** 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 { |
| |
| if (value.length > 20) { |
| throw new MathArithmeticException(LocalizedFormats.NUMBER_TOO_LARGE, value.length, 20); |
| } |
| for (int i = 0; i < value.length; ++i) { |
| final double[] y = value[i].clone(); |
| if (i > 1) { |
| double inv = 1.0 / Factorial.value(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 |
| */ |
| @Override |
| 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 |
| */ |
| @Override |
| 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); |
| } |
| |
| } |