Author: luc
Date: Mon Feb 25 10:59:10 2013
New Revision: 1449657
URL: http://svn.apache.org/r1449657
Log:
Added Hermite interpolator for ExtendFieldElement instances.
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolator.java
(with props)
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolatorTest.java
(with props)
Modified:
commons/proper/math/trunk/src/changes/changes.xml
Modified: commons/proper/math/trunk/src/changes/changes.xml
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/changes/changes.xml?rev=1449657&r1=1449656&r2=1449657&view=diff
==============================================================================
--- commons/proper/math/trunk/src/changes/changes.xml (original)
+++ commons/proper/math/trunk/src/changes/changes.xml Mon Feb 25 10:59:10 2013
@@ -56,6 +56,9 @@ This is a minor release: It combines bug
way such as to allow drop-in replacement of the v3.1[.1] JAR file.
">
<action dev="luc" type="add" >
+ Added Hermite interpolator for ExtendFieldElement instances.
+ </action>
+ <action dev="luc" type="add" >
Added ExtendFieldElement interface to represent anything that is
real number like, implemented by both Decimal64, Dfp and
DerivativeStructure.
</action>
Added:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolator.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolator.java?rev=1449657&view=auto
==============================================================================
---
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolator.java
(added)
+++
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolator.java
Mon Feb 25 10:59:10 2013
@@ -0,0 +1,152 @@
+/*
+ * 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.FieldElement;
+import org.apache.commons.math3.exception.MathArithmeticException;
+import org.apache.commons.math3.exception.NoDataException;
+import org.apache.commons.math3.exception.ZeroException;
+import org.apache.commons.math3.exception.util.LocalizedFormats;
+import org.apache.commons.math3.util.MathArrays;
+
+/** 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>
+ *
+ * @version $Id$
+ * @since 3.2
+ */
+public class FieldHermiteInterpolator<T extends FieldElement<T>> {
+
+ /** Sample abscissae. */
+ private final List<T> abscissae;
+
+ /** Top diagonal of the divided differences array. */
+ private final List<T[]> topDiagonal;
+
+ /** Bottom diagonal of the divided differences array. */
+ private final List<T[]> bottomDiagonal;
+
+ /** Create an empty interpolator.
+ */
+ public FieldHermiteInterpolator() {
+ this.abscissae = new ArrayList<T>();
+ this.topDiagonal = new ArrayList<T[]>();
+ this.bottomDiagonal = new ArrayList<T[]>();
+ }
+
+ /** 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 T x, final T[] ... value)
+ throws ZeroException, MathArithmeticException {
+
+ T factorial = x.getField().getOne();
+ for (int i = 0; i < value.length; ++i) {
+
+ final T[] y = value[i].clone();
+ if (i > 1) {
+ factorial = factorial.multiply(i);
+ final T inv = factorial.reciprocal();
+ for (int j = 0; j < y.length; ++j) {
+ y[j] = y[j].multiply(inv);
+ }
+ }
+
+ // update the bottom diagonal of the divided differences array
+ final int n = abscissae.size();
+ bottomDiagonal.add(n - i, y);
+ T[] bottom0 = y;
+ for (int j = i; j < n; ++j) {
+ final T[] bottom1 = bottomDiagonal.get(n - (j + 1));
+ if (x.equals(abscissae.get(n - (j + 1)))) {
+ throw new
ZeroException(LocalizedFormats.DUPLICATED_ABSCISSA_DIVISION_BY_ZERO, x);
+ }
+ final T inv = x.subtract(abscissae.get(n - (j +
1))).reciprocal();
+ for (int k = 0; k < y.length; ++k) {
+ bottom1[k] = inv.multiply(bottom0[k].subtract(bottom1[k]));
+ }
+ bottom0 = bottom1;
+ }
+
+ // update the top diagonal of the divided differences array
+ topDiagonal.add(bottom0.clone());
+
+ // update the abscissae array
+ abscissae.add(x);
+
+ }
+
+ }
+
+ /** Interpolate value at a specified abscissa.
+ * @param x interpolation abscissa
+ * @return interpolated value
+ * @exception NoDataException if sample is empty
+ */
+ public T[] value(T x) throws NoDataException {
+
+ // safety check
+ if (abscissae.isEmpty()) {
+ throw new
NoDataException(LocalizedFormats.EMPTY_INTERPOLATION_SAMPLE);
+ }
+
+ final T[] value = MathArrays.buildArray(x.getField(),
topDiagonal.get(0).length);
+ T valueCoeff = x.getField().getOne();
+ for (int i = 0; i < topDiagonal.size(); ++i) {
+ T[] dividedDifference = topDiagonal.get(i);
+ for (int k = 0; k < value.length; ++k) {
+ value[k] =
value[k].add(dividedDifference[k].multiply(valueCoeff));
+ }
+ final T deltaX = x.subtract(abscissae.get(i));
+ valueCoeff = valueCoeff.multiply(deltaX);
+ }
+
+ return value;
+
+ }
+
+}
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolator.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolator.java
------------------------------------------------------------------------------
svn:keywords = "Author Date Id Revision"
Added:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolatorTest.java
URL:
http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolatorTest.java?rev=1449657&view=auto
==============================================================================
---
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolatorTest.java
(added)
+++
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolatorTest.java
Mon Feb 25 10:59:10 2013
@@ -0,0 +1,245 @@
+/*
+ * 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.dfp.Dfp;
+import org.apache.commons.math3.dfp.DfpField;
+import org.apache.commons.math3.exception.NoDataException;
+import org.apache.commons.math3.fraction.BigFraction;
+import org.apache.commons.math3.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+public class FieldHermiteInterpolatorTest {
+
+ @Test
+ public void testZero() {
+ FieldHermiteInterpolator<BigFraction> interpolator = new
FieldHermiteInterpolator<BigFraction>();
+ interpolator.addSamplePoint(new BigFraction(0), new BigFraction[] {
new BigFraction(0) });
+ for (int x = -10; x < 10; x++) {
+ BigFraction y = interpolator.value(new BigFraction(x))[0];
+ Assert.assertEquals(BigFraction.ZERO, y);
+ }
+ }
+
+ @Test
+ public void testQuadratic() {
+ FieldHermiteInterpolator<BigFraction> interpolator = new
FieldHermiteInterpolator<BigFraction>();
+ interpolator.addSamplePoint(new BigFraction(0), new BigFraction[] {
new BigFraction(2) });
+ interpolator.addSamplePoint(new BigFraction(1), new BigFraction[] {
new BigFraction(0) });
+ interpolator.addSamplePoint(new BigFraction(2), new BigFraction[] {
new BigFraction(0) });
+ for (double x = -10; x < 10; x += 1.0) {
+ BigFraction y = interpolator.value(new BigFraction(x))[0];
+ Assert.assertEquals((x - 1) * (x - 2), y.doubleValue(), 1.0e-15);
+ }
+ }
+
+ @Test
+ public void testMixedDerivatives() {
+ FieldHermiteInterpolator<BigFraction> interpolator = new
FieldHermiteInterpolator<BigFraction>();
+ interpolator.addSamplePoint(new BigFraction(0), new BigFraction[] {
new BigFraction(1) }, new BigFraction[] { new BigFraction(2) });
+ interpolator.addSamplePoint(new BigFraction(1), new BigFraction[] {
new BigFraction(4) });
+ interpolator.addSamplePoint(new BigFraction(2), new BigFraction[] {
new BigFraction(5) }, new BigFraction[] { new BigFraction(2) });
+ Assert.assertEquals(new BigFraction(1), interpolator.value(new
BigFraction(0))[0]);
+ Assert.assertEquals(new BigFraction(4), interpolator.value(new
BigFraction(1))[0]);
+ Assert.assertEquals(new BigFraction(5), interpolator.value(new
BigFraction(2))[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);
+ }
+
+ DfpField field = new DfpField(30);
+ Dfp step = field.getOne().divide(field.newDfp(10));
+ FieldHermiteInterpolator<Dfp> interpolator = new
FieldHermiteInterpolator<Dfp>();
+ for (int j = 0; j < 1 + maxDegree; ++j) {
+ Dfp x = field.newDfp(j).multiply(step);
+ Dfp[] values = new Dfp[p.length];
+ for (int k = 0; k < p.length; ++k) {
+ values[k] = field.newDfp(p[k].value(x.getReal()));
+ }
+ interpolator.addSamplePoint(x, values);
+ }
+
+ for (int j = 0; j < 20; ++j) {
+ Dfp x = field.newDfp(j).multiply(step);
+ Dfp[] values = interpolator.value(x);
+ Assert.assertEquals(p.length, values.length);
+ for (int k = 0; k < p.length; ++k) {
+ Assert.assertEquals(p[k].value(x.getReal()),
+ values[k].getReal(),
+ 1.0e-8 *
FastMath.abs(p[k].value(x.getReal())));
+ }
+ }
+
+ }
+
+ }
+
+ @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);
+ }
+
+ DfpField field = new DfpField(30);
+ Dfp step = field.getOne().divide(field.newDfp(10));
+ FieldHermiteInterpolator<Dfp> interpolator = new
FieldHermiteInterpolator<Dfp>();
+ for (int j = 0; j < 1 + maxDegree / 2; ++j) {
+ Dfp x = field.newDfp(j).multiply(step);
+ Dfp[] values = new Dfp[p.length];
+ Dfp[] derivatives = new Dfp[p.length];
+ for (int k = 0; k < p.length; ++k) {
+ values[k] = field.newDfp(p[k].value(x.getReal()));
+ derivatives[k] =
field.newDfp(pPrime[k].value(x.getReal()));
+ }
+ interpolator.addSamplePoint(x, values, derivatives);
+ }
+
+ Dfp h = step.divide(field.newDfp(100000));
+ for (int j = 0; j < 20; ++j) {
+ Dfp x = field.newDfp(j).multiply(step);
+ Dfp[] y = interpolator.value(x);
+ Dfp[] yP = interpolator.value(x.add(h));
+ Dfp[] yM = interpolator.value(x.subtract(h));
+ Assert.assertEquals(p.length, y.length);
+ for (int k = 0; k < p.length; ++k) {
+ Assert.assertEquals(p[k].value(x.getReal()),
+ y[k].getReal(),
+ 1.0e-8 *
FastMath.abs(p[k].value(x.getReal())));
+ Assert.assertEquals(pPrime[k].value(x.getReal()),
+
yP[k].subtract(yM[k]).divide(h.multiply(2)).getReal(),
+ 4.0e-8 *
FastMath.abs(p[k].value(x.getReal())));
+ }
+ System.out.println();
+ }
+
+ }
+ }
+
+ @Test
+ public void testSine() {
+ DfpField field = new DfpField(30);
+ FieldHermiteInterpolator<Dfp> interpolator = new
FieldHermiteInterpolator<Dfp>();
+ for (Dfp x = field.getZero(); x.getReal() < FastMath.PI; x =
x.add(0.5)) {
+ interpolator.addSamplePoint(x, new Dfp[] { x.sin() });
+ }
+ for (Dfp x = field.newDfp(0.1); x.getReal() < 2.9; x = x.add(0.01)) {
+ Dfp y = interpolator.value(x)[0];
+ Assert.assertEquals( x.sin().getReal(), y.getReal(), 3.5e-5);
+ }
+ }
+
+ @Test
+ public void testSquareRoot() {
+ DfpField field = new DfpField(30);
+ FieldHermiteInterpolator<Dfp> interpolator = new
FieldHermiteInterpolator<Dfp>();
+ for (Dfp x = field.getOne(); x.getReal() < 3.6; x = x.add(0.5)) {
+ interpolator.addSamplePoint(x, new Dfp[] { x.sqrt() });
+ }
+ for (Dfp x = field.newDfp(1.1); x.getReal() < 3.5; x = x.add(0.01)) {
+ Dfp y = interpolator.value(x)[0];
+ Assert.assertEquals(x.sqrt().getReal(), y.getReal(), 1.5e-4);
+ }
+ }
+
+ @Test
+ public void testWikipedia() {
+ // this test corresponds to the example from Wikipedia page:
+ // http://en.wikipedia.org/wiki/Hermite_interpolation
+ FieldHermiteInterpolator<BigFraction> interpolator = new
FieldHermiteInterpolator<BigFraction>();
+ interpolator.addSamplePoint(new BigFraction(-1),
+ new BigFraction[] { new BigFraction( 2) },
+ new BigFraction[] { new BigFraction(-8) },
+ new BigFraction[] { new BigFraction(56) });
+ interpolator.addSamplePoint(new BigFraction( 0),
+ new BigFraction[] { new BigFraction( 1) },
+ new BigFraction[] { new BigFraction( 0) },
+ new BigFraction[] { new BigFraction( 0) });
+ interpolator.addSamplePoint(new BigFraction( 1),
+ new BigFraction[] { new BigFraction( 2) },
+ new BigFraction[] { new BigFraction( 8) },
+ new BigFraction[] { new BigFraction(56) });
+ for (BigFraction x = new BigFraction(-1); x.doubleValue() <= 1.0; x =
x.add(new BigFraction(1, 8))) {
+ BigFraction y = interpolator.value(x)[0];
+ BigFraction x2 = x.multiply(x);
+ BigFraction x4 = x2.multiply(x2);
+ BigFraction x8 = x4.multiply(x4);
+ Assert.assertEquals(x8.add(new BigFraction(1)), y);
+ }
+ }
+
+ @Test
+ public void testOnePointParabola() {
+ FieldHermiteInterpolator<BigFraction> interpolator = new
FieldHermiteInterpolator<BigFraction>();
+ interpolator.addSamplePoint(new BigFraction(0),
+ new BigFraction[] { new BigFraction(1) },
+ new BigFraction[] { new BigFraction(1) },
+ new BigFraction[] { new BigFraction(2) });
+ for (BigFraction x = new BigFraction(-1); x.doubleValue() <= 1.0; x =
x.add(new BigFraction(1, 8))) {
+ BigFraction y = interpolator.value(x)[0];
+
Assert.assertEquals(BigFraction.ONE.add(x.multiply(BigFraction.ONE.add(x))), y);
+ }
+ }
+
+ 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=NoDataException.class)
+ public void testEmptySample() {
+ new FieldHermiteInterpolator<BigFraction>().value(BigFraction.ZERO);
+ }
+
+ @Test(expected=IllegalArgumentException.class)
+ public void testDuplicatedAbscissa() {
+ FieldHermiteInterpolator<BigFraction> interpolator = new
FieldHermiteInterpolator<BigFraction>();
+ interpolator.addSamplePoint(new BigFraction(1), new BigFraction[] {
new BigFraction(0) });
+ interpolator.addSamplePoint(new BigFraction(1), new BigFraction[] {
new BigFraction(1) });
+ }
+
+}
+
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolatorTest.java
------------------------------------------------------------------------------
svn:eol-style = native
Propchange:
commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/interpolation/FieldHermiteInterpolatorTest.java
------------------------------------------------------------------------------
svn:keywords = "Author Date Id Revision"
