commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r506592 - in /jakarta/commons/proper/math/trunk/src: java/org/apache/commons/math/ java/org/apache/commons/math/analysis/ test/org/apache/commons/math/
Date Mon, 12 Feb 2007 19:27:17 GMT
Author: luc
Date: Mon Feb 12 11:27:16 2007
New Revision: 506592

URL: http://svn.apache.org/viewvc?view=rev&rev=506592
Log:
Added and used a specialized exception for duplicate abscissas in sampled functions

Added:
    jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/DuplicateSampleAbscissaException.java
  (with props)
    jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/DuplicateSampleAbscissaExceptionTest.java
  (with props)
Modified:
    jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java
  (contents, props changed)
    jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java
  (contents, props changed)

Added: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/DuplicateSampleAbscissaException.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/DuplicateSampleAbscissaException.java?view=auto&rev=506592
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/DuplicateSampleAbscissaException.java
(added)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/DuplicateSampleAbscissaException.java
Mon Feb 12 11:27:16 2007
@@ -0,0 +1,47 @@
+/*
+ * 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.math;
+
+/**
+ * Exeption thrown when a sample contains several entries at the same abscissa.
+ * @version $Revision:$
+ */
+public class DuplicateSampleAbscissaException extends MathException  {
+    
+    /** Serializable version identifier */
+    private static final long serialVersionUID = -2271007547170169872L;
+
+    /**
+     * Construct an exception indicating the duplicate abscissa.
+     * @param abscissa duplicate abscissa
+     * @param i1 index of one entry having the duplicate abscissa
+     * @param i2 index of another entry having the duplicate abscissa
+     */
+    public DuplicateSampleAbscissaException(double abscissa, int i1, int i2) {
+        super("Abscissa {0} is duplicated at both indices {1} and {2}",
+              new Object[] { new Double(abscissa), new Integer(i1), new Integer(i2) });
+    }
+
+    /**
+     * Get the duplicate abscissa.
+     * @return duplicate abscissa
+     */
+    public double getDuplicateAbscissa() {
+        return ((Double) getArguments()[0]).doubleValue();
+    }
+    
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/DuplicateSampleAbscissaException.java
------------------------------------------------------------------------------
    svn:eol-style = native

Modified: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java?view=diff&rev=506592&r1=506591&r2=506592
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java
(original)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java
Mon Feb 12 11:27:16 2007
@@ -1,122 +1,122 @@
-/*
- * 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.math.analysis;
-
-import java.io.Serializable;
-import org.apache.commons.math.MathException;
-
-/**
- * Implements the <a href="
- * "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html">
- * Divided Difference Algorithm</a> for interpolation of real univariate
- * functions. For reference, see <b>Introduction to Numerical Analysis</b>,
- * ISBN 038795452X, chapter 2.
- * <p>
- * The actual code of Neville's evalution is in PolynomialFunctionLagrangeForm,
- * this class provides an easy-to-use interface to it.
- *
- * @version $Revision$ $Date$
- */
-public class DividedDifferenceInterpolator implements UnivariateRealInterpolator,
-    Serializable {
-
-    /** serializable version identifier */
-    static final long serialVersionUID = 107049519551235069L;
-
-    /**
-     * Computes an interpolating function for the data set.
-     *
-     * @param x the interpolating points array
-     * @param y the interpolating values array
-     * @return a function which interpolates the data set
-     * @throws MathException if arguments are invalid
-     */
-    public UnivariateRealFunction interpolate(double x[], double y[]) throws
-        MathException {
-
-        /**
-         * a[] and c[] are defined in the general formula of Newton form:
-         * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
-         *        a[n](x-c[0])(x-c[1])...(x-c[n-1])
-         */
-        double a[], c[];
-
-        PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
-
-        /**
-         * When used for interpolation, the Newton form formula becomes
-         * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
-         *        f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
-         * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
-         * <p>
-         * Note x[], y[], a[] have the same length but c[]'s size is one less.
-         */
-        c = new double[x.length-1];
-        for (int i = 0; i < c.length; i++) {
-            c[i] = x[i];
-        }
-        a = computeDividedDifference(x, y);
-
-        PolynomialFunctionNewtonForm p;
-        p = new PolynomialFunctionNewtonForm(a, c);
-        return p;
-    }
-
-    /**
-     * Returns a copy of the divided difference array.
-     * <p> 
-     * The divided difference array is defined recursively by <pre>
-     * f[x0] = f(x0)
-     * f[x0,x1,...,xk] = (f(x1,...,xk) - f(x0,...,x[k-1])) / (xk - x0)
-     * </pre><p>
-     * The computational complexity is O(N^2).
-     *
-     * @return a fresh copy of the divided difference array
-     * @throws MathException if any abscissas coincide
-     */
-    protected static double[] computeDividedDifference(double x[], double y[])
-        throws MathException {
-
-        int i, j, n;
-        double divdiff[], a[], denominator;
-
-        PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
-
-        n = x.length;
-        divdiff = new double[n];
-        for (i = 0; i < n; i++) {
-            divdiff[i] = y[i];      // initialization
-        }
-
-        a = new double [n];
-        a[0] = divdiff[0];
-        for (i = 1; i < n; i++) {
-            for (j = 0; j < n-i; j++) {
-                denominator = x[j+i] - x[j];
-                if (denominator == 0.0) {
-                    // This happens only when two abscissas are identical.
-                    throw new MathException
-                        ("Identical abscissas cause division by zero.");
-                }
-                divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
-            }
-            a[i] = divdiff[0];
-        }
-
-        return a;
-    }
-}
+/*
+ * 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.math.analysis;
+
+import java.io.Serializable;
+
+import org.apache.commons.math.DuplicateSampleAbscissaException;
+
+/**
+ * Implements the <a href="
+ * "http://mathworld.wolfram.com/NewtonsDividedDifferenceInterpolationFormula.html">
+ * Divided Difference Algorithm</a> for interpolation of real univariate
+ * functions. For reference, see <b>Introduction to Numerical Analysis</b>,
+ * ISBN 038795452X, chapter 2.
+ * <p>
+ * The actual code of Neville's evalution is in PolynomialFunctionLagrangeForm,
+ * this class provides an easy-to-use interface to it.
+ *
+ * @version $Revision$ $Date$
+ */
+public class DividedDifferenceInterpolator implements UnivariateRealInterpolator,
+    Serializable {
+
+    /** serializable version identifier */
+    private static final long serialVersionUID = 107049519551235069L;
+
+    /**
+     * Computes an interpolating function for the data set.
+     *
+     * @param x the interpolating points array
+     * @param y the interpolating values array
+     * @return a function which interpolates the data set
+     * @throws DuplicateSampleAbscissaException if arguments are invalid
+     */
+    public UnivariateRealFunction interpolate(double x[], double y[]) throws
+        DuplicateSampleAbscissaException {
+
+        /**
+         * a[] and c[] are defined in the general formula of Newton form:
+         * p(x) = a[0] + a[1](x-c[0]) + a[2](x-c[0])(x-c[1]) + ... +
+         *        a[n](x-c[0])(x-c[1])...(x-c[n-1])
+         */
+        double a[], c[];
+
+        PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
+
+        /**
+         * When used for interpolation, the Newton form formula becomes
+         * p(x) = f[x0] + f[x0,x1](x-x0) + f[x0,x1,x2](x-x0)(x-x1) + ... +
+         *        f[x0,x1,...,x[n-1]](x-x0)(x-x1)...(x-x[n-2])
+         * Therefore, a[k] = f[x0,x1,...,xk], c[k] = x[k].
+         * <p>
+         * Note x[], y[], a[] have the same length but c[]'s size is one less.
+         */
+        c = new double[x.length-1];
+        for (int i = 0; i < c.length; i++) {
+            c[i] = x[i];
+        }
+        a = computeDividedDifference(x, y);
+
+        PolynomialFunctionNewtonForm p;
+        p = new PolynomialFunctionNewtonForm(a, c);
+        return p;
+    }
+
+    /**
+     * Returns a copy of the divided difference array.
+     * <p> 
+     * The divided difference array is defined recursively by <pre>
+     * f[x0] = f(x0)
+     * f[x0,x1,...,xk] = (f(x1,...,xk) - f(x0,...,x[k-1])) / (xk - x0)
+     * </pre><p>
+     * The computational complexity is O(N^2).
+     *
+     * @return a fresh copy of the divided difference array
+     * @throws DuplicateSampleAbscissaException if any abscissas coincide
+     */
+    protected static double[] computeDividedDifference(double x[], double y[])
+        throws DuplicateSampleAbscissaException {
+
+        int i, j, n;
+        double divdiff[], a[], denominator;
+
+        PolynomialFunctionLagrangeForm.verifyInterpolationArray(x, y);
+
+        n = x.length;
+        divdiff = new double[n];
+        for (i = 0; i < n; i++) {
+            divdiff[i] = y[i];      // initialization
+        }
+
+        a = new double [n];
+        a[0] = divdiff[0];
+        for (i = 1; i < n; i++) {
+            for (j = 0; j < n-i; j++) {
+                denominator = x[j+i] - x[j];
+                if (denominator == 0.0) {
+                    // This happens only when two abscissas are identical.
+                    throw new DuplicateSampleAbscissaException(x[j], j, j+i);
+                }
+                divdiff[j] = (divdiff[j+1] - divdiff[j]) / denominator;
+            }
+            a[i] = divdiff[0];
+        }
+
+        return a;
+    }
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/DividedDifferenceInterpolator.java
            ('svn:executable' removed)

Modified: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java?view=diff&rev=506592&r1=506591&r2=506592
==============================================================================
--- jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java
(original)
+++ jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java
Mon Feb 12 11:27:16 2007
@@ -1,291 +1,295 @@
-/*
- * 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.math.analysis;
-
-import java.io.Serializable;
-import org.apache.commons.math.FunctionEvaluationException;
-
-/**
- * Implements the representation of a real polynomial function in
- * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
- * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
- * Analysis</b>, ISBN 038795452X, chapter 2.
- * <p>
- * The approximated function should be smooth enough for Lagrange polynomial
- * to work well. Otherwise, consider using splines instead.
- *
- * @version $Revision$ $Date$
- */
-public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction,
-    Serializable {
-
-    /** serializable version identifier */
-    static final long serialVersionUID = -3965199246151093920L;
-
-    /**
-     * The coefficients of the polynomial, ordered by degree -- i.e.
-     * coefficients[0] is the constant term and coefficients[n] is the 
-     * coefficient of x^n where n is the degree of the polynomial.
-     */
-    private double coefficients[];
-
-    /**
-     * Interpolating points (abscissas) and the function values at these points.
-     */
-    private double x[], y[];
-
-    /**
-     * Whether the polynomial coefficients are available.
-     */
-    private boolean coefficientsComputed;
-
-    /**
-     * Construct a Lagrange polynomial with the given abscissas and function
-     * values. The order of interpolating points are not important.
-     * <p>
-     * The constructor makes copy of the input arrays and assigns them.
-     * 
-     * @param x interpolating points
-     * @param y function values at interpolating points
-     * @throws IllegalArgumentException if input arrays are not valid
-     */
-    PolynomialFunctionLagrangeForm(double x[], double y[]) throws
-        IllegalArgumentException {
-
-        verifyInterpolationArray(x, y);
-        this.x = new double[x.length];
-        this.y = new double[y.length];
-        System.arraycopy(x, 0, this.x, 0, x.length);
-        System.arraycopy(y, 0, this.y, 0, y.length);
-        coefficientsComputed = false;
-    }
-
-    /**
-     * Calculate the function value at the given point.
-     *
-     * @param z the point at which the function value is to be computed
-     * @return the function value
-     * @throws FunctionEvaluationException if a runtime error occurs
-     * @see UnivariateRealFunction#value(double)
-     */
-    public double value(double z) throws FunctionEvaluationException {
-       return evaluate(x, y, z);
-    }
-
-    /**
-     * Returns the degree of the polynomial.
-     * 
-     * @return the degree of the polynomial
-     */
-    public int degree() {
-        return x.length - 1;
-    }
-
-    /**
-     * Returns a copy of the interpolating points array.
-     * <p>
-     * Changes made to the returned copy will not affect the polynomial.
-     * 
-     * @return a fresh copy of the interpolating points array
-     */
-    public double[] getInterpolatingPoints() {
-        double[] out = new double[x.length];
-        System.arraycopy(x, 0, out, 0, x.length);
-        return out;
-    }
-
-    /**
-     * Returns a copy of the interpolating values array.
-     * <p>
-     * Changes made to the returned copy will not affect the polynomial.
-     * 
-     * @return a fresh copy of the interpolating values array
-     */
-    public double[] getInterpolatingValues() {
-        double[] out = new double[y.length];
-        System.arraycopy(y, 0, out, 0, y.length);
-        return out;
-    }
-
-    /**
-     * Returns a copy of the coefficients array.
-     * <p>
-     * Changes made to the returned copy will not affect the polynomial.
-     * 
-     * @return a fresh copy of the coefficients array
-     */
-    public double[] getCoefficients() {
-        if (!coefficientsComputed) {
-            computeCoefficients();
-        }
-        double[] out = new double[coefficients.length];
-        System.arraycopy(coefficients, 0, out, 0, coefficients.length);
-        return out;
-    }
-
-    /**
-     * Evaluate the Lagrange polynomial using 
-     * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
-     * Neville's Algorithm</a>. It takes O(N^2) time.
-     * <p>
-     * This function is made public static so that users can call it directly
-     * without instantiating PolynomialFunctionLagrangeForm object.
-     *
-     * @param x the interpolating points array
-     * @param y the interpolating values array
-     * @param z the point at which the function value is to be computed
-     * @return the function value
-     * @throws FunctionEvaluationException if a runtime error occurs
-     * @throws IllegalArgumentException if inputs are not valid
-     */
-    public static double evaluate(double x[], double y[], double z) throws
-        FunctionEvaluationException, IllegalArgumentException {
-
-        int i, j, n, nearest = 0;
-        double value, c[], d[], tc, td, divider, w, dist, min_dist;
-
-        verifyInterpolationArray(x, y);
-
-        n = x.length;
-        c = new double[n];
-        d = new double[n];
-        min_dist = Double.POSITIVE_INFINITY;
-        for (i = 0; i < n; i++) {
-            // initialize the difference arrays
-            c[i] = y[i];
-            d[i] = y[i];
-            // find out the abscissa closest to z
-            dist = Math.abs(z - x[i]);
-            if (dist < min_dist) {
-                nearest = i;
-                min_dist = dist;
-            }
-        }
-
-        // initial approximation to the function value at z
-        value = y[nearest];
-
-        for (i = 1; i < n; i++) {
-            for (j = 0; j < n-i; j++) {
-                tc = x[j] - z;
-                td = x[i+j] - z;
-                divider = x[j] - x[i+j];
-                if (divider == 0.0) {
-                    // This happens only when two abscissas are identical.
-                    throw new FunctionEvaluationException(z, 
-                        "Identical abscissas cause division by zero: x[" +
-                        i + "] = x[" + (i+j) + "] = " + x[i]);
-                }
-                // update the difference arrays
-                w = (c[j+1] - d[j]) / divider;
-                c[j] = tc * w;
-                d[j] = td * w;
-            }
-            // sum up the difference terms to get the final value
-            if (nearest < 0.5*(n-i+1)) {
-                value += c[nearest];    // fork down
-            } else {
-                nearest--;
-                value += d[nearest];    // fork up
-            }
-        }
-
-        return value;
-    }
-
-    /**
-     * Calculate the coefficients of Lagrange polynomial from the
-     * interpolation data. It takes O(N^2) time.
-     * <p>
-     * Note this computation can be ill-conditioned. Use with caution
-     * and only when it is necessary.
-     *
-     * @throws ArithmeticException if any abscissas coincide
-     */
-    protected void computeCoefficients() throws ArithmeticException {
-        int i, j, n;
-        double c[], tc[], d, t;
-
-        n = degree() + 1;
-        coefficients = new double[n];
-        for (i = 0; i < n; i++) {
-            coefficients[i] = 0.0;
-        }
-
-        // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
-        c = new double[n+1];
-        c[0] = 1.0;
-        for (i = 0; i < n; i++) {
-            for (j = i; j > 0; j--) {
-                c[j] = c[j-1] - c[j] * x[i];
-            }
-            c[0] *= (-x[i]);
-            c[i+1] = 1;
-        }
-
-        tc = new double[n];
-        for (i = 0; i < n; i++) {
-            // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
-            d = 1;
-            for (j = 0; j < n; j++) {
-                if (i != j) {
-                    d *= (x[i] - x[j]);
-                }
-            }
-            if (d == 0.0) {
-                // This happens only when two abscissas are identical.
-                throw new ArithmeticException
-                    ("Identical abscissas cause division by zero.");
-            }
-            t = y[i] / d;
-            // Lagrange polynomial is the sum of n terms, each of which is a
-            // polynomial of degree n-1. tc[] are the coefficients of the i-th
-            // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
-            tc[n-1] = c[n];     // actually c[n] = 1
-            coefficients[n-1] += t * tc[n-1];
-            for (j = n-2; j >= 0; j--) {
-                tc[j] = c[j+1] + tc[j+1] * x[i];
-                coefficients[j] += t * tc[j];
-            }
-        }
-
-        coefficientsComputed = true;
-    }
-
-    /**
-     * Verifies that the interpolation arrays are valid.
-     * <p>
-     * The interpolating points must be distinct. However it is not
-     * verified here, it is checked in evaluate() and computeCoefficients().
-     * 
-     * @throws IllegalArgumentException if not valid
-     * @see #evaluate(double[], double[], double)
-     * @see #computeCoefficients()
-     */
-    protected static void verifyInterpolationArray(double x[], double y[]) throws
-        IllegalArgumentException {
-
-        if (x.length < 2 || y.length < 2) {
-            throw new IllegalArgumentException
-                ("Interpolation requires at least two points.");
-        }
-        if (x.length != y.length) {
-            throw new IllegalArgumentException
-                ("Abscissa and value arrays must have the same length.");
-        }
-    }
-}
+/*
+ * 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.math.analysis;
+
+import java.io.Serializable;
+
+import org.apache.commons.math.DuplicateSampleAbscissaException;
+import org.apache.commons.math.FunctionEvaluationException;
+
+/**
+ * Implements the representation of a real polynomial function in
+ * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
+ * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
+ * Analysis</b>, ISBN 038795452X, chapter 2.
+ * <p>
+ * The approximated function should be smooth enough for Lagrange polynomial
+ * to work well. Otherwise, consider using splines instead.
+ *
+ * @version $Revision$ $Date$
+ */
+public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction,
+    Serializable {
+
+    /** serializable version identifier */
+    static final long serialVersionUID = -3965199246151093920L;
+
+    /**
+     * The coefficients of the polynomial, ordered by degree -- i.e.
+     * coefficients[0] is the constant term and coefficients[n] is the 
+     * coefficient of x^n where n is the degree of the polynomial.
+     */
+    private double coefficients[];
+
+    /**
+     * Interpolating points (abscissas) and the function values at these points.
+     */
+    private double x[], y[];
+
+    /**
+     * Whether the polynomial coefficients are available.
+     */
+    private boolean coefficientsComputed;
+
+    /**
+     * Construct a Lagrange polynomial with the given abscissas and function
+     * values. The order of interpolating points are not important.
+     * <p>
+     * The constructor makes copy of the input arrays and assigns them.
+     * 
+     * @param x interpolating points
+     * @param y function values at interpolating points
+     * @throws IllegalArgumentException if input arrays are not valid
+     */
+    PolynomialFunctionLagrangeForm(double x[], double y[]) throws
+        IllegalArgumentException {
+
+        verifyInterpolationArray(x, y);
+        this.x = new double[x.length];
+        this.y = new double[y.length];
+        System.arraycopy(x, 0, this.x, 0, x.length);
+        System.arraycopy(y, 0, this.y, 0, y.length);
+        coefficientsComputed = false;
+    }
+
+    /**
+     * Calculate the function value at the given point.
+     *
+     * @param z the point at which the function value is to be computed
+     * @return the function value
+     * @throws FunctionEvaluationException if a runtime error occurs
+     * @see UnivariateRealFunction#value(double)
+     */
+    public double value(double z) throws FunctionEvaluationException {
+        try {
+            return evaluate(x, y, z);
+        } catch (DuplicateSampleAbscissaException e) {
+            throw new FunctionEvaluationException(z, e.getPattern(), e.getArguments(), e);
+        }
+    }
+
+    /**
+     * Returns the degree of the polynomial.
+     * 
+     * @return the degree of the polynomial
+     */
+    public int degree() {
+        return x.length - 1;
+    }
+
+    /**
+     * Returns a copy of the interpolating points array.
+     * <p>
+     * Changes made to the returned copy will not affect the polynomial.
+     * 
+     * @return a fresh copy of the interpolating points array
+     */
+    public double[] getInterpolatingPoints() {
+        double[] out = new double[x.length];
+        System.arraycopy(x, 0, out, 0, x.length);
+        return out;
+    }
+
+    /**
+     * Returns a copy of the interpolating values array.
+     * <p>
+     * Changes made to the returned copy will not affect the polynomial.
+     * 
+     * @return a fresh copy of the interpolating values array
+     */
+    public double[] getInterpolatingValues() {
+        double[] out = new double[y.length];
+        System.arraycopy(y, 0, out, 0, y.length);
+        return out;
+    }
+
+    /**
+     * Returns a copy of the coefficients array.
+     * <p>
+     * Changes made to the returned copy will not affect the polynomial.
+     * 
+     * @return a fresh copy of the coefficients array
+     */
+    public double[] getCoefficients() {
+        if (!coefficientsComputed) {
+            computeCoefficients();
+        }
+        double[] out = new double[coefficients.length];
+        System.arraycopy(coefficients, 0, out, 0, coefficients.length);
+        return out;
+    }
+
+    /**
+     * Evaluate the Lagrange polynomial using 
+     * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
+     * Neville's Algorithm</a>. It takes O(N^2) time.
+     * <p>
+     * This function is made public static so that users can call it directly
+     * without instantiating PolynomialFunctionLagrangeForm object.
+     *
+     * @param x the interpolating points array
+     * @param y the interpolating values array
+     * @param z the point at which the function value is to be computed
+     * @return the function value
+     * @throws DuplicateSampleAbscissaException if the sample has duplicate abscissas
+     * @throws IllegalArgumentException if inputs are not valid
+     */
+    public static double evaluate(double x[], double y[], double z) throws
+        DuplicateSampleAbscissaException, IllegalArgumentException {
+
+        int i, j, n, nearest = 0;
+        double value, c[], d[], tc, td, divider, w, dist, min_dist;
+
+        verifyInterpolationArray(x, y);
+
+        n = x.length;
+        c = new double[n];
+        d = new double[n];
+        min_dist = Double.POSITIVE_INFINITY;
+        for (i = 0; i < n; i++) {
+            // initialize the difference arrays
+            c[i] = y[i];
+            d[i] = y[i];
+            // find out the abscissa closest to z
+            dist = Math.abs(z - x[i]);
+            if (dist < min_dist) {
+                nearest = i;
+                min_dist = dist;
+            }
+        }
+
+        // initial approximation to the function value at z
+        value = y[nearest];
+
+        for (i = 1; i < n; i++) {
+            for (j = 0; j < n-i; j++) {
+                tc = x[j] - z;
+                td = x[i+j] - z;
+                divider = x[j] - x[i+j];
+                if (divider == 0.0) {
+                    // This happens only when two abscissas are identical.
+                    throw new DuplicateSampleAbscissaException(x[i], i, i+j);
+                }
+                // update the difference arrays
+                w = (c[j+1] - d[j]) / divider;
+                c[j] = tc * w;
+                d[j] = td * w;
+            }
+            // sum up the difference terms to get the final value
+            if (nearest < 0.5*(n-i+1)) {
+                value += c[nearest];    // fork down
+            } else {
+                nearest--;
+                value += d[nearest];    // fork up
+            }
+        }
+
+        return value;
+    }
+
+    /**
+     * Calculate the coefficients of Lagrange polynomial from the
+     * interpolation data. It takes O(N^2) time.
+     * <p>
+     * Note this computation can be ill-conditioned. Use with caution
+     * and only when it is necessary.
+     *
+     * @throws ArithmeticException if any abscissas coincide
+     */
+    protected void computeCoefficients() throws ArithmeticException {
+        int i, j, n;
+        double c[], tc[], d, t;
+
+        n = degree() + 1;
+        coefficients = new double[n];
+        for (i = 0; i < n; i++) {
+            coefficients[i] = 0.0;
+        }
+
+        // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
+        c = new double[n+1];
+        c[0] = 1.0;
+        for (i = 0; i < n; i++) {
+            for (j = i; j > 0; j--) {
+                c[j] = c[j-1] - c[j] * x[i];
+            }
+            c[0] *= (-x[i]);
+            c[i+1] = 1;
+        }
+
+        tc = new double[n];
+        for (i = 0; i < n; i++) {
+            // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
+            d = 1;
+            for (j = 0; j < n; j++) {
+                if (i != j) {
+                    d *= (x[i] - x[j]);
+                }
+            }
+            if (d == 0.0) {
+                // This happens only when two abscissas are identical.
+                throw new ArithmeticException
+                    ("Identical abscissas cause division by zero.");
+            }
+            t = y[i] / d;
+            // Lagrange polynomial is the sum of n terms, each of which is a
+            // polynomial of degree n-1. tc[] are the coefficients of the i-th
+            // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
+            tc[n-1] = c[n];     // actually c[n] = 1
+            coefficients[n-1] += t * tc[n-1];
+            for (j = n-2; j >= 0; j--) {
+                tc[j] = c[j+1] + tc[j+1] * x[i];
+                coefficients[j] += t * tc[j];
+            }
+        }
+
+        coefficientsComputed = true;
+    }
+
+    /**
+     * Verifies that the interpolation arrays are valid.
+     * <p>
+     * The interpolating points must be distinct. However it is not
+     * verified here, it is checked in evaluate() and computeCoefficients().
+     * 
+     * @throws IllegalArgumentException if not valid
+     * @see #evaluate(double[], double[], double)
+     * @see #computeCoefficients()
+     */
+    protected static void verifyInterpolationArray(double x[], double y[]) throws
+        IllegalArgumentException {
+
+        if (x.length < 2 || y.length < 2) {
+            throw new IllegalArgumentException
+                ("Interpolation requires at least two points.");
+        }
+        if (x.length != y.length) {
+            throw new IllegalArgumentException
+                ("Abscissa and value arrays must have the same length.");
+        }
+    }
+}

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: jakarta/commons/proper/math/trunk/src/java/org/apache/commons/math/analysis/PolynomialFunctionLagrangeForm.java
            ('svn:executable' removed)

Added: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/DuplicateSampleAbscissaExceptionTest.java
URL: http://svn.apache.org/viewvc/jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/DuplicateSampleAbscissaExceptionTest.java?view=auto&rev=506592
==============================================================================
--- jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/DuplicateSampleAbscissaExceptionTest.java
(added)
+++ jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/DuplicateSampleAbscissaExceptionTest.java
Mon Feb 12 11:27:16 2007
@@ -0,0 +1,38 @@
+/*
+ * 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.math;
+
+import java.util.Locale;
+
+import junit.framework.TestCase;
+
+/**
+ * @version $Revision:$
+ */
+public class DuplicateSampleAbscissaExceptionTest extends TestCase {
+    
+    public void testConstructor(){
+        DuplicateSampleAbscissaException ex = new DuplicateSampleAbscissaException(1.2, 10,
11);
+        assertNull(ex.getCause());
+        assertNotNull(ex.getMessage());
+        assertTrue(ex.getMessage().indexOf("1.2") > 0);
+        assertEquals(1.2, ex.getDuplicateAbscissa(), 0);
+        assertFalse(ex.getMessage().equals(ex.getMessage(Locale.FRENCH)));
+    }
+    
+}

Propchange: jakarta/commons/proper/math/trunk/src/test/org/apache/commons/math/DuplicateSampleAbscissaExceptionTest.java
------------------------------------------------------------------------------
    svn:eol-style = native



---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Mime
View raw message