You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2013/02/25 11:59:11 UTC

svn commit: r1449657 - in /commons/proper/math/trunk/src: changes/ main/java/org/apache/commons/math3/analysis/interpolation/ test/java/org/apache/commons/math3/analysis/interpolation/

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"