diff --git a/NOTICE.txt b/NOTICE.txt index be2d74481..8c313273f 100644 --- a/NOTICE.txt +++ b/NOTICE.txt @@ -54,7 +54,13 @@ All rights reserved The LocalizedFormatsTest class in the unit tests is an adapted version of the OrekitMessagesTest class from the orekit library distributed under the terms of the Apache 2 licence. Original source copyright: -Copyright 2010 CS Communication & Systèmes +Copyright 2010 CS Systèmes d'Information +=============================================================================== + +The HermiteInterpolator class and its corresponding test have been imported from +the orekit library distributed under the terms of the Apache 2 licence. Original +source copyright: +Copyright 2010-2012 CS Systèmes d'Information =============================================================================== The complete text of licenses and disclaimers associated with the the original diff --git a/src/changes/changes.xml b/src/changes/changes.xml index 6919d6dc8..f23c413f9 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -52,6 +52,11 @@ If the output is not quite correct, check for invisible trailing spaces! + + A new HermiteInterpolator class allow interpolation of vector-valued + functions using both values and derivatives of the function at sample + points. + Parameterized "CurveFitter" class (package "o.a.c.m.optimization.fitting") with the type of the fitting function. Updated subclasses "PolynomialFitter", diff --git a/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java b/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java new file mode 100644 index 000000000..397014bc7 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolator.java @@ -0,0 +1,254 @@ +/* + * 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.math3.analysis.interpolation; + +import java.util.ArrayList; +import java.util.List; + +import org.apache.commons.math3.analysis.DifferentiableUnivariateVectorFunction; +import org.apache.commons.math3.analysis.UnivariateVectorFunction; +import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; +import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.MathIllegalStateException; +import org.apache.commons.math3.exception.util.LocalizedFormats; +import org.apache.commons.math3.util.ArithmeticUtils; + +/** Polynomial interpolator using both sample values and sample derivatives. + *

+ * 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 polynomial 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. + *

+ *

+ * This class has been imported from the Orekit space flight dynamics library + * also distributed under the terms of the Apache License V2. Original copyright + * is: Copyright 2002-2012 CS Systèmes d'Information. + *

+ * @version $Id$ + * @since 3.1 + */ +public class HermiteInterpolator implements DifferentiableUnivariateVectorFunction { + + /** Sample abscissae. */ + private final List abscissae; + + /** Top diagonal of the divided differences array. */ + private final List topDiagonal; + + /** Bottom diagonal of the divided differences array. */ + private final List bottomDiagonal; + + /** Create an empty interpolator. + */ + public HermiteInterpolator() { + this.abscissae = new ArrayList(); + this.topDiagonal = new ArrayList(); + this.bottomDiagonal = new ArrayList(); + } + + /** Add a sample point. + *

+ * 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. + *

+ *

+ * The point abscissae for all calls must be different. + *

+ * @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 MathIllegalArgumentException if the abscissa is equals to a previously + * added sample point + */ + public void addSamplePoint(final double x, final double[] ... value) + throws MathIllegalArgumentException { + + for (int i = 0; i < value.length; ++i) { + + final double[] y = value[i].clone(); + if (i > 1) { + double inv = 1.0 / ArithmeticUtils.factorial(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 MathIllegalArgumentException(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 MathIllegalStateException if sample is empty + */ + public PolynomialFunction[] getPolynomials() + throws MathIllegalStateException { + + // 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. + *

+ * 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. + *

+ * @param x interpolation abscissa + * @return interpolated value + * @exception MathIllegalStateException if sample is empty + */ + public double[] value(double x) + throws MathIllegalStateException { + + // 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 first derivative at a specified abscissa. + *

+ * Calling this method is equivalent to call the {@link PolynomialFunction#value(double) + * value} methods of the derivatives of all polynomials returned by {@link + * #getPolynomials() getPolynomials}, except it builds neither the intermediate + * polynomials nor their derivatives, so this method is faster and numerically more stable. + *

+ * @param x interpolation abscissa + * @return interpolated derivative + * @exception MathIllegalStateException if sample is empty + */ + public double[] derivative(double x) + throws MathIllegalStateException { + + // safety check + checkInterpolation(); + + final double[] derivative = new double[topDiagonal.get(0).length]; + double valueCoeff = 1; + double derivativeCoeff = 0; + for (int i = 0; i < topDiagonal.size(); ++i) { + double[] dividedDifference = topDiagonal.get(i); + for (int k = 0; k < derivative.length; ++k) { + derivative[k] += dividedDifference[k] * derivativeCoeff; + } + final double deltaX = x - abscissae.get(i); + derivativeCoeff = valueCoeff + derivativeCoeff * deltaX; + valueCoeff *= deltaX; + } + + return derivative; + + } + + /** {@inheritDoc}} */ + public UnivariateVectorFunction derivative() { + return new UnivariateVectorFunction() { + + /** {@inheritDoc}} */ + public double[] value(double x) { + return derivative(x); + } + + }; + } + + /** Check interpolation can be performed. + * @exception MathIllegalStateException if interpolation cannot be performed + * because sample is empty + */ + private void checkInterpolation() throws MathIllegalStateException { + if (abscissae.isEmpty()) { + throw new MathIllegalStateException(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); + } + +} diff --git a/src/main/java/org/apache/commons/math3/exception/util/LocalizedFormats.java b/src/main/java/org/apache/commons/math3/exception/util/LocalizedFormats.java index db6c5463e..1ca563575 100644 --- a/src/main/java/org/apache/commons/math3/exception/util/LocalizedFormats.java +++ b/src/main/java/org/apache/commons/math3/exception/util/LocalizedFormats.java @@ -88,9 +88,10 @@ public enum LocalizedFormats implements Localizable { DIMENSIONS_MISMATCH("dimensions mismatch"), /* keep */ DISCRETE_CUMULATIVE_PROBABILITY_RETURNED_NAN("Discrete cumulative probability function returned NaN for argument {0}"), DISTRIBUTION_NOT_LOADED("distribution not loaded"), - DUPLICATED_ABSCISSA("Abscissa {0} is duplicated at both indices {1} and {2}"), + DUPLICATED_ABSCISSA_DIVISION_BY_ZERO("duplicated abscissa {0} causes division by zero"), ELITISM_RATE("elitism rate ({0})"), EMPTY_CLUSTER_IN_K_MEANS("empty cluster in k-means"), + EMPTY_INTERPOLATION_SAMPLE("sample for interpolation is empty"), EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY("empty polynomials coefficients array"), /* keep */ EMPTY_SELECTED_COLUMN_INDEX_ARRAY("empty selected column index array"), EMPTY_SELECTED_ROW_INDEX_ARRAY("empty selected row index array"), @@ -112,7 +113,6 @@ public enum LocalizedFormats implements Localizable { GCD_OVERFLOW_32_BITS("overflow: gcd({0}, {1}) is 2^31"), GCD_OVERFLOW_64_BITS("overflow: gcd({0}, {1}) is 2^63"), HOLE_BETWEEN_MODELS_TIME_RANGES("{0} wide hole between models time ranges"), - IDENTICAL_ABSCISSAS_DIVISION_BY_ZERO("identical abscissas x[{0}] == x[{1}] == {2} cause division by zero"), ILL_CONDITIONED_OPERATOR("condition number {1} is too high "), INDEX_LARGER_THAN_MAX("the index specified: {0} is larger than the current maximal index {1}"), INDEX_NOT_POSITIVE("index ({0}) is not positive"), diff --git a/src/main/resources/assets/org/apache/commons/math3/exception/util/LocalizedFormats_fr.properties b/src/main/resources/assets/org/apache/commons/math3/exception/util/LocalizedFormats_fr.properties index 5123c2baa..de0a771f2 100644 --- a/src/main/resources/assets/org/apache/commons/math3/exception/util/LocalizedFormats_fr.properties +++ b/src/main/resources/assets/org/apache/commons/math3/exception/util/LocalizedFormats_fr.properties @@ -60,9 +60,10 @@ DIMENSIONS_MISMATCH_SIMPLE = {0} != {1} DIMENSIONS_MISMATCH = dimensions incoh\u00e9rentes DISCRETE_CUMULATIVE_PROBABILITY_RETURNED_NAN = Discr\u00e8tes fonction de probabilit\u00e9 cumulative retourn\u00e9 NaN \u00e0 l''argument de {0} DISTRIBUTION_NOT_LOADED = aucune distribution n''a \u00e9t\u00e9 charg\u00e9e -DUPLICATED_ABSCISSA = Abscisse {0} dupliqu\u00e9e aux indices {1} et {2} +DUPLICATED_ABSCISSA_DIVISION_BY_ZERO = la duplication de l''abscisse {0} engendre une division par z\u00e9ro ELITISM_RATE = proportion d''\u00e9litisme ({0}) EMPTY_CLUSTER_IN_K_MEANS = groupe vide dans l''algorithme des k-moyennes +EMPTY_INTERPOLATION_SAMPLE = \u00e9chantillon d''interpolation vide EMPTY_POLYNOMIALS_COEFFICIENTS_ARRAY = tableau de coefficients polyn\u00f4miaux vide EMPTY_SELECTED_COLUMN_INDEX_ARRAY = tableau des indices de colonnes s\u00e9lectionn\u00e9es vide EMPTY_SELECTED_ROW_INDEX_ARRAY = tableau des indices de lignes s\u00e9lectionn\u00e9es vide @@ -84,7 +85,6 @@ FUNCTION_NOT_POLYNOMIAL = la fonction n''est pas p\u00f4lynomiale GCD_OVERFLOW_32_BITS = d\u00e9passement de capacit\u00e9 : le PGCD de {0} et {1} vaut 2^31 GCD_OVERFLOW_64_BITS = d\u00e9passement de capacit\u00e9 : le PGCD de {0} et {1} vaut 2^63 HOLE_BETWEEN_MODELS_TIME_RANGES = trou de longueur {0} entre les domaines temporels des mod\u00e8les -IDENTICAL_ABSCISSAS_DIVISION_BY_ZERO = les abscisses identiques x[{0}] == x[{1}] == {2} engendrent une division par z\u00e9ro ILL_CONDITIONED_OPERATOR = le conditionnement {1} est trop \u00e9lev\u00e9 INDEX_LARGER_THAN_MAX = l''index sp\u00e9cifi\u00e9 ({0}) d\u00e9passe l''index maximal courant ({1}) INDEX_NOT_POSITIVE = l''indice ({0}) n''est pas positif diff --git a/src/site/xdoc/userguide/analysis.xml b/src/site/xdoc/userguide/analysis.xml index 4904d4a67..f0e761dc6 100644 --- a/src/site/xdoc/userguide/analysis.xml +++ b/src/site/xdoc/userguide/analysis.xml @@ -436,6 +436,27 @@ System.out println("f(" + interpolationX + ") = " + interpolatedY); It has been described in William Dudziak's MS thesis.

+

+ Hermite interpolation + is an interpolation method that can use derivatives in addition to function values at sample points. The HermiteInterpolator + class implements this method for vector-valued functions. The sampling points can have any spacing (there are + no requirements for a regular grid) and some points may provide derivatives while others don't provide them + (or provide derivatives to a smaller order). Points are added one at a time, as shown in the following example: +

+ HermiteInterpolator interpolator = new HermiteInterpolator; +// at x = 0, we provide both value and first derivative +interpolator.addSamplePoint(0.0, new double[] { 1.0 }, new double[] { 2.0 }); +// at x = 1, we provide only function value +interpolator.addSamplePoint(1.0, new double[] { 4.0 }); +// at x = 2, we provide both value and first derivative +interpolator.addSamplePoint(2.0, new double[] { 5.0 }, new double[] { 2.0 }); +// should print "value at x = 0.5: 2.5625" +System.out.println("value at x = 0.5: " + interpolator.value(0.5)[0]); +// should print "derivative at x = 0.5: 3.5" +System.out.println("derivative at x = 0.5: " + interpolator.derivative(0.5)[0]); +// should print "interpolation polynomial: 1 + 2 x + 4 x^2 - 4 x^3 + x^4" +System.out.println("interpolation polynomial: " + interpolator.getPolynomials()[0]);

A BivariateGridInterpolator is used to find a bivariate real-valued diff --git a/src/test/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolatorTest.java b/src/test/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolatorTest.java new file mode 100644 index 000000000..46228e1f8 --- /dev/null +++ b/src/test/java/org/apache/commons/math3/analysis/interpolation/HermiteInterpolatorTest.java @@ -0,0 +1,246 @@ +/* + * 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.math3.analysis.interpolation; + +import java.util.Random; + +import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; +import org.apache.commons.math3.util.FastMath; +import org.junit.Assert; +import org.junit.Test; + +public class HermiteInterpolatorTest { + + @Test + public void testZero() { + HermiteInterpolator interpolator = new HermiteInterpolator(); + interpolator.addSamplePoint(0.0, new double[] { 0.0 }); + for (double x = -10; x < 10; x += 1.0) { + Assert.assertEquals(0.0, interpolator.value(x)[0], 1.0e-15); + Assert.assertEquals(0.0, interpolator.derivative(x)[0], 1.0e-15); + } + checkPolynomial(new PolynomialFunction(new double[] { 0.0 }), + interpolator.getPolynomials()[0]); + } + + @Test + public void testQuadratic() { + HermiteInterpolator interpolator = new HermiteInterpolator(); + interpolator.addSamplePoint(0.0, new double[] { 2.0 }); + interpolator.addSamplePoint(1.0, new double[] { 0.0 }); + interpolator.addSamplePoint(2.0, new double[] { 0.0 }); + for (double x = -10; x < 10; x += 1.0) { + Assert.assertEquals((x - 1.0) * (x - 2.0), interpolator.value(x)[0], 1.0e-15); + Assert.assertEquals(2 * x - 3.0, interpolator.derivative(x)[0], 1.0e-15); + } + checkPolynomial(new PolynomialFunction(new double[] { 2.0, -3.0, 1.0 }), + interpolator.getPolynomials()[0]); + } + + @Test + public void testMixedDerivatives() { + HermiteInterpolator interpolator = new HermiteInterpolator(); + interpolator.addSamplePoint(0.0, new double[] { 1.0 }, new double[] { 2.0 }); + interpolator.addSamplePoint(1.0, new double[] { 4.0 }); + interpolator.addSamplePoint(2.0, new double[] { 5.0 }, new double[] { 2.0 }); + Assert.assertEquals(4, interpolator.getPolynomials()[0].degree()); + Assert.assertEquals(1.0, interpolator.value(0.0)[0], 1.0e-15); + Assert.assertEquals(2.0, interpolator.derivative(0.0)[0], 1.0e-15); + Assert.assertEquals(4.0, interpolator.value(1.0)[0], 1.0e-15); + Assert.assertEquals(5.0, interpolator.value(2.0)[0], 1.0e-15); + Assert.assertEquals(2.0, interpolator.derivative(2.0)[0], 1.0e-15); + checkPolynomial(new PolynomialFunction(new double[] { 1.0, 2.0, 4.0, -4.0, 1.0 }), + interpolator.getPolynomials()[0]); + } + + @Test + public void testRandomPolynomialsValuesOnly() { + + Random random = new Random(0x42b1e7dbd361a932l); + + for (int i = 0; i < 100; ++i) { + + int maxDegree = 0; + PolynomialFunction[] p = new PolynomialFunction[5]; + for (int k = 0; k < p.length; ++k) { + int degree = random.nextInt(7); + p[k] = randomPolynomial(degree, random); + maxDegree = FastMath.max(maxDegree, degree); + } + + HermiteInterpolator interpolator = new HermiteInterpolator(); + for (int j = 0; j < 1 + maxDegree; ++j) { + double x = 0.1 * j; + double[] values = new double[p.length]; + for (int k = 0; k < p.length; ++k) { + values[k] = p[k].value(x); + } + interpolator.addSamplePoint(x, values); + } + + for (double x = 0; x < 2; x += 0.1) { + double[] values = interpolator.value(x); + Assert.assertEquals(p.length, values.length); + for (int k = 0; k < p.length; ++k) { + Assert.assertEquals(p[k].value(x), values[k], 1.0e-8 * FastMath.abs(p[k].value(x))); + } + } + + PolynomialFunction[] result = interpolator.getPolynomials(); + for (int k = 0; k < p.length; ++k) { + checkPolynomial(p[k], result[k]); + } + + } + } + + @Test + public void testRandomPolynomialsFirstDerivative() { + + Random random = new Random(0x570803c982ca5d3bl); + + for (int i = 0; i < 100; ++i) { + + int maxDegree = 0; + PolynomialFunction[] p = new PolynomialFunction[5]; + PolynomialFunction[] pPrime = new PolynomialFunction[5]; + for (int k = 0; k < p.length; ++k) { + int degree = random.nextInt(7); + p[k] = randomPolynomial(degree, random); + pPrime[k] = p[k].polynomialDerivative(); + maxDegree = FastMath.max(maxDegree, degree); + } + + HermiteInterpolator interpolator = new HermiteInterpolator(); + for (int j = 0; j < 1 + maxDegree / 2; ++j) { + double x = 0.1 * j; + double[] values = new double[p.length]; + double[] derivatives = new double[p.length]; + for (int k = 0; k < p.length; ++k) { + values[k] = p[k].value(x); + derivatives[k] = pPrime[k].value(x); + } + interpolator.addSamplePoint(x, values, derivatives); + } + + for (double x = 0; x < 2; x += 0.1) { + double[] values = interpolator.value(x); + double[] derivatives = interpolator.derivative(x); + Assert.assertEquals(p.length, values.length); + for (int k = 0; k < p.length; ++k) { + Assert.assertEquals(p[k].value(x), values[k], 1.0e-8 * FastMath.abs(p[k].value(x))); + Assert.assertEquals(pPrime[k].value(x), derivatives[k], 4.0e-8 * FastMath.abs(p[k].value(x))); + } + } + + PolynomialFunction[] result = interpolator.getPolynomials(); + for (int k = 0; k < p.length; ++k) { + checkPolynomial(p[k], result[k]); + } + + } + } + + @Test + public void testSine() { + HermiteInterpolator interpolator = new HermiteInterpolator(); + for (double x = 0; x < FastMath.PI; x += 0.5) { + interpolator.addSamplePoint(x, new double[] { FastMath.sin(x) }); + } + for (double x = 0.1; x <= 2.9; x += 0.01) { + Assert.assertEquals(FastMath.sin(x), interpolator.value(x)[0], 3.5e-5); + Assert.assertEquals(FastMath.cos(x), interpolator.derivative(x)[0], 1.3e-4); + } + } + + @Test + public void testSquareRoot() { + HermiteInterpolator interpolator = new HermiteInterpolator(); + for (double x = 1.0; x < 3.6; x += 0.5) { + interpolator.addSamplePoint(x, new double[] { FastMath.sqrt(x) }); + } + for (double x = 1.1; x < 3.5; x += 0.01) { + Assert.assertEquals(FastMath.sqrt(x), interpolator.value(x)[0], 1.5e-4); + Assert.assertEquals(0.5 / FastMath.sqrt(x), interpolator.derivative(x)[0], 8.5e-4); + } + } + + @Test + public void testWikipedia() { + // this test corresponds to the example from Wikipedia page: + // http://en.wikipedia.org/wiki/Hermite_interpolation + HermiteInterpolator interpolator = new HermiteInterpolator(); + interpolator.addSamplePoint(-1, new double[] { 2 }, new double[] { -8 }, new double[] { 56 }); + interpolator.addSamplePoint( 0, new double[] { 1 }, new double[] { 0 }, new double[] { 0 }); + interpolator.addSamplePoint( 1, new double[] { 2 }, new double[] { 8 }, new double[] { 56 }); + for (double x = -1.0; x <= 1.0; x += 0.125) { + double x2 = x * x; + double x4 = x2 * x2; + double x8 = x4 * x4; + Assert.assertEquals(x8 + 1, interpolator.value(x)[0], 1.0e-15); + Assert.assertEquals(8 * x4 * x2 * x, interpolator.derivative(x)[0], 1.0e-15); + } + checkPolynomial(new PolynomialFunction(new double[] { 1, 0, 0, 0, 0, 0, 0, 0, 1 }), + interpolator.getPolynomials()[0]); + } + + @Test + public void testOnePointParabola() { + HermiteInterpolator interpolator = new HermiteInterpolator(); + interpolator.addSamplePoint(0, new double[] { 1 }, new double[] { 1 }, new double[] { 2 }); + for (double x = -1.0; x <= 1.0; x += 0.125) { + Assert.assertEquals(1 + x * (1 + x), interpolator.value(x)[0], 1.0e-15); + Assert.assertEquals(1 + 2 * x, interpolator.derivative(x)[0], 1.0e-15); + } + checkPolynomial(new PolynomialFunction(new double[] { 1, 1, 1 }), + interpolator.getPolynomials()[0]); + } + + private PolynomialFunction randomPolynomial(int degree, Random random) { + double[] coeff = new double[ 1 + degree]; + for (int j = 0; j < degree; ++j) { + coeff[j] = random.nextDouble(); + } + return new PolynomialFunction(coeff); + } + + @Test(expected=IllegalStateException.class) + public void testEmptySample() { + new HermiteInterpolator().value(0.0); + } + + @Test(expected=IllegalArgumentException.class) + public void testDuplicatedAbscissa() { + HermiteInterpolator interpolator = new HermiteInterpolator(); + interpolator.addSamplePoint(1.0, new double[] { 0.0 }); + interpolator.addSamplePoint(1.0, new double[] { 1.0 }); + } + + private void checkPolynomial(PolynomialFunction expected, PolynomialFunction result) { + Assert.assertTrue(result.degree() >= expected.degree()); + double[] cE = expected.getCoefficients(); + double[] cR = result.getCoefficients(); + for (int i = 0; i < cE.length; ++i) { + Assert.assertEquals(cE[i], cR[i], 1.0e-8 * FastMath.abs(cE[i])); + } + for (int i = cE.length; i < cR.length; ++i) { + Assert.assertEquals(0.0, cR[i], 1.0e-9); + } + } + +} +