001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.math4.legacy.analysis.interpolation; 018 019import java.util.ArrayList; 020import java.util.Arrays; 021import java.util.List; 022 023import org.apache.commons.math4.legacy.analysis.differentiation.DerivativeStructure; 024import org.apache.commons.math4.legacy.analysis.differentiation.UnivariateDifferentiableVectorFunction; 025import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunction; 026import org.apache.commons.math4.legacy.exception.MathArithmeticException; 027import org.apache.commons.math4.legacy.exception.NoDataException; 028import org.apache.commons.math4.legacy.exception.ZeroException; 029import org.apache.commons.math4.legacy.exception.util.LocalizedFormats; 030import org.apache.commons.numbers.combinatorics.Factorial; 031 032/** Polynomial interpolator using both sample values and sample derivatives. 033 * <p> 034 * The interpolation polynomials match all sample points, including both values 035 * and provided derivatives. There is one polynomial for each component of 036 * the values vector. All polynomials have the same degree. The degree of the 037 * polynomials depends on the number of points and number of derivatives at each 038 * point. For example the interpolation polynomials for n sample points without 039 * any derivatives all have degree n-1. The interpolation polynomials for n 040 * sample points with the two extreme points having value and first derivative 041 * and the remaining points having value only all have degree n+1. The 042 * interpolation polynomial for n sample points with value, first and second 043 * derivative for all points all have degree 3n-1. 044 * </p> 045 * 046 * @since 3.1 047 */ 048public class HermiteInterpolator implements UnivariateDifferentiableVectorFunction { 049 050 /** Sample abscissae. */ 051 private final List<Double> abscissae; 052 053 /** Top diagonal of the divided differences array. */ 054 private final List<double[]> topDiagonal; 055 056 /** Bottom diagonal of the divided differences array. */ 057 private final List<double[]> bottomDiagonal; 058 059 /** Create an empty interpolator. 060 */ 061 public HermiteInterpolator() { 062 this.abscissae = new ArrayList<>(); 063 this.topDiagonal = new ArrayList<>(); 064 this.bottomDiagonal = new ArrayList<>(); 065 } 066 067 /** Add a sample point. 068 * <p> 069 * This method must be called once for each sample point. It is allowed to 070 * mix some calls with values only with calls with values and first 071 * derivatives. 072 * </p> 073 * <p> 074 * The point abscissae for all calls <em>must</em> be different. 075 * </p> 076 * @param x abscissa of the sample point 077 * @param value value and derivatives of the sample point 078 * (if only one row is passed, it is the value, if two rows are 079 * passed the first one is the value and the second the derivative 080 * and so on) 081 * @exception ZeroException if the abscissa difference between added point 082 * and a previous point is zero (i.e. the two points are at same abscissa) 083 * @exception MathArithmeticException if the number of derivatives is larger 084 * than 20, which prevents computation of a factorial 085 */ 086 public void addSamplePoint(final double x, final double[] ... value) 087 throws ZeroException, MathArithmeticException { 088 089 if (value.length > 20) { 090 throw new MathArithmeticException(LocalizedFormats.NUMBER_TOO_LARGE, value.length, 20); 091 } 092 for (int i = 0; i < value.length; ++i) { 093 final double[] y = value[i].clone(); 094 if (i > 1) { 095 double inv = 1.0 / Factorial.value(i); 096 for (int j = 0; j < y.length; ++j) { 097 y[j] *= inv; 098 } 099 } 100 101 // update the bottom diagonal of the divided differences array 102 final int n = abscissae.size(); 103 bottomDiagonal.add(n - i, y); 104 double[] bottom0 = y; 105 for (int j = i; j < n; ++j) { 106 final double[] bottom1 = bottomDiagonal.get(n - (j + 1)); 107 final double inv = 1.0 / (x - abscissae.get(n - (j + 1))); 108 if (Double.isInfinite(inv)) { 109 throw new ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x); 110 } 111 for (int k = 0; k < y.length; ++k) { 112 bottom1[k] = inv * (bottom0[k] - bottom1[k]); 113 } 114 bottom0 = bottom1; 115 } 116 117 // update the top diagonal of the divided differences array 118 topDiagonal.add(bottom0.clone()); 119 120 // update the abscissae array 121 abscissae.add(x); 122 } 123 } 124 125 /** Compute the interpolation polynomials. 126 * @return interpolation polynomials array 127 * @exception NoDataException if sample is empty 128 */ 129 public PolynomialFunction[] getPolynomials() 130 throws NoDataException { 131 132 // safety check 133 checkInterpolation(); 134 135 // iteration initialization 136 final PolynomialFunction zero = polynomial(0); 137 PolynomialFunction[] polynomials = new PolynomialFunction[topDiagonal.get(0).length]; 138 for (int i = 0; i < polynomials.length; ++i) { 139 polynomials[i] = zero; 140 } 141 PolynomialFunction coeff = polynomial(1); 142 143 // build the polynomials by iterating on the top diagonal of the divided differences array 144 for (int i = 0; i < topDiagonal.size(); ++i) { 145 double[] tdi = topDiagonal.get(i); 146 for (int k = 0; k < polynomials.length; ++k) { 147 polynomials[k] = polynomials[k].add(coeff.multiply(polynomial(tdi[k]))); 148 } 149 coeff = coeff.multiply(polynomial(-abscissae.get(i), 1.0)); 150 } 151 152 return polynomials; 153 } 154 155 /** Interpolate value at a specified abscissa. 156 * <p> 157 * Calling this method is equivalent to call the {@link PolynomialFunction#value(double) 158 * value} methods of all polynomials returned by {@link #getPolynomials() getPolynomials}, 159 * except it does not build the intermediate polynomials, so this method is faster and 160 * numerically more stable. 161 * </p> 162 * @param x interpolation abscissa 163 * @return interpolated value 164 * @exception NoDataException if sample is empty 165 */ 166 @Override 167 public double[] value(double x) throws NoDataException { 168 169 // safety check 170 checkInterpolation(); 171 172 final double[] value = new double[topDiagonal.get(0).length]; 173 double valueCoeff = 1; 174 for (int i = 0; i < topDiagonal.size(); ++i) { 175 double[] dividedDifference = topDiagonal.get(i); 176 for (int k = 0; k < value.length; ++k) { 177 value[k] += dividedDifference[k] * valueCoeff; 178 } 179 final double deltaX = x - abscissae.get(i); 180 valueCoeff *= deltaX; 181 } 182 183 return value; 184 } 185 186 /** Interpolate value at a specified abscissa. 187 * <p> 188 * Calling this method is equivalent to call the {@link 189 * PolynomialFunction#value(DerivativeStructure) value} methods of all polynomials 190 * returned by {@link #getPolynomials() getPolynomials}, except it does not build the 191 * intermediate polynomials, so this method is faster and numerically more stable. 192 * </p> 193 * @param x interpolation abscissa 194 * @return interpolated value 195 * @exception NoDataException if sample is empty 196 */ 197 @Override 198 public DerivativeStructure[] value(final DerivativeStructure x) 199 throws NoDataException { 200 201 // safety check 202 checkInterpolation(); 203 204 final DerivativeStructure[] value = new DerivativeStructure[topDiagonal.get(0).length]; 205 Arrays.fill(value, x.getField().getZero()); 206 DerivativeStructure valueCoeff = x.getField().getOne(); 207 for (int i = 0; i < topDiagonal.size(); ++i) { 208 double[] dividedDifference = topDiagonal.get(i); 209 for (int k = 0; k < value.length; ++k) { 210 value[k] = value[k].add(valueCoeff.multiply(dividedDifference[k])); 211 } 212 final DerivativeStructure deltaX = x.subtract(abscissae.get(i)); 213 valueCoeff = valueCoeff.multiply(deltaX); 214 } 215 216 return value; 217 } 218 219 /** Check interpolation can be performed. 220 * @exception NoDataException if interpolation cannot be performed 221 * because sample is empty 222 */ 223 private void checkInterpolation() throws NoDataException { 224 if (abscissae.isEmpty()) { 225 throw new NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE); 226 } 227 } 228 229 /** Create a polynomial from its coefficients. 230 * @param c polynomials coefficients 231 * @return polynomial 232 */ 233 private PolynomialFunction polynomial(double ... c) { 234 return new PolynomialFunction(c); 235 } 236}