commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1370951 [2/2] - in /commons/proper/math/trunk/src: changes/ main/java/org/apache/commons/math3/analysis/differentiation/ main/java/org/apache/commons/math3/util/ test/java/org/apache/commons/math3/analysis/differentiation/
Date Wed, 08 Aug 2012 20:33:43 GMT
Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/package-info.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/package-info.java?rev=1370951&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/package-info.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/package-info.java Wed Aug  8 20:33:43 2012
@@ -0,0 +1,35 @@
+/*
+ * 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.
+ */
+/**
+ *
+ * <p>
+ *   This package holds the main interfaces and basic building block classes
+ *   dealing with differentiation.
+ *   The core class is {@link DerivativeStructure} which holds the value and
+ *   the differentials of a function. This class handles some arbitrary number
+ *   of free parameters and arbitrary derivation order. It is used both as
+ *   the input and the output type for the {@link UnivariateDifferential}
+ *   interface. Any differentiable function should implement this interface.
+ * </p>
+ * <p>
+ *   The {@link UnivariateDifferentiator} interface defines a way to differentation
+ *   a simple {@link org.apache.commons.math3.analysis.UnivariateFunction
+ *   univariate function} and get a {@link differential function}.
+ * </p>
+ *
+ */
+package org.apache.commons.math3.analysis.differentiation;

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/package-info.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/analysis/differentiation/package-info.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/FastMath.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/FastMath.java?rev=1370951&r1=1370950&r2=1370951&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/FastMath.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/util/FastMath.java Wed Aug  8 20:33:43 2012
@@ -18,6 +18,9 @@ package org.apache.commons.math3.util;
 
 import java.io.PrintStream;
 
+import org.apache.commons.math3.exception.NotPositiveException;
+import org.apache.commons.math3.exception.util.LocalizedFormats;
+
 /**
  * Faster, more accurate, portable alternative to {@link Math} and
  * {@link StrictMath} for large scale computation.
@@ -1145,7 +1148,7 @@ public class FastMath {
             /* Normalize the subnormal number. */
             bits <<= 1;
             while ( (bits & 0x0010000000000000L) == 0) {
-                --exp;
+                exp--;
                 bits <<= 1;
             }
         }
@@ -1165,9 +1168,8 @@ public class FastMath {
                 xa = aa;
                 xb = ab;
 
-                final double[] lnCoef_last = LN_QUICK_COEF[LN_QUICK_COEF.length - 1];
-                double ya = lnCoef_last[0];
-                double yb = lnCoef_last[1];
+                double ya = LN_QUICK_COEF[LN_QUICK_COEF.length-1][0];
+                double yb = LN_QUICK_COEF[LN_QUICK_COEF.length-1][1];
 
                 for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
                     /* Multiply a = y * x */
@@ -1179,9 +1181,8 @@ public class FastMath {
                     yb = aa - ya + ab;
 
                     /* Add  a = y + lnQuickCoef */
-                    final double[] lnCoef_i = LN_QUICK_COEF[i];
-                    aa = ya + lnCoef_i[0];
-                    ab = yb + lnCoef_i[1];
+                    aa = ya + LN_QUICK_COEF[i][0];
+                    ab = yb + LN_QUICK_COEF[i][1];
                     /* Split y = a */
                     tmp = aa * HEX_40000000;
                     ya = aa + tmp - tmp;
@@ -1201,7 +1202,7 @@ public class FastMath {
         }
 
         // lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm < ln(2)
-        final double[] lnm = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
+        double lnm[] = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
 
         /*
     double epsilon = x / Double.longBitsToDouble(bits & 0xfffffc0000000000L);
@@ -1212,7 +1213,7 @@ public class FastMath {
         // y is the most significant 10 bits of the mantissa
         //double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L);
         //double epsilon = (x - y) / y;
-        final double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));
+        double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));
 
         double lnza = 0.0;
         double lnzb = 0.0;
@@ -1226,15 +1227,14 @@ public class FastMath {
             double xb = ab;
 
             /* Need a more accurate epsilon, so adjust the division. */
-            final double numer = bits & 0x3ffffffffffL;
-            final double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
+            double numer = bits & 0x3ffffffffffL;
+            double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
             aa = numer - xa*denom - xb * denom;
             xb += aa / denom;
 
             /* Remez polynomial evaluation */
-            final double[] lnCoef_last = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1];
-            double ya = lnCoef_last[0];
-            double yb = lnCoef_last[1];
+            double ya = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1][0];
+            double yb = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1][1];
 
             for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) {
                 /* Multiply a = y * x */
@@ -1246,9 +1246,8 @@ public class FastMath {
                 yb = aa - ya + ab;
 
                 /* Add  a = y + lnHiPrecCoef */
-                final double[] lnCoef_i = LN_HI_PREC_COEF[i];
-                aa = ya + lnCoef_i[0];
-                ab = yb + lnCoef_i[1];
+                aa = ya + LN_HI_PREC_COEF[i][0];
+                ab = yb + LN_HI_PREC_COEF[i][1];
                 /* Split y = a */
                 tmp = aa * HEX_40000000;
                 ya = aa + tmp - tmp;
@@ -1582,6 +1581,34 @@ public class FastMath {
 
 
     /**
+     * Raise a double to an int power.
+     *
+     * @param d Number to raise.
+     * @param e Exponent.
+     * @return d<sup>e</sup>
+     */
+    public static double pow(double d, int e) {
+        if (e == 0) {
+            return 1.0;
+        } else if (e < 0) {
+            e = -e;
+            d = 1.0 / d;
+        }
+
+        double result = 1;
+        double d2p    = d;
+        while (e != 0) {
+            if ((e & 0x1) != 0) {
+                result *= d2p;
+            }
+            d2p *= d2p;
+            e = e >> 1;
+        }
+
+        return result;
+    }
+
+    /**
      *  Computes sin(x) - x, where |x| < 1/16.
      *  Use a Remez polynomial approximation.
      *  @param x a number smaller than 1/16

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DSCompilerTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DSCompilerTest.java?rev=1370951&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DSCompilerTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DSCompilerTest.java Wed Aug  8 20:33:43 2012
@@ -0,0 +1,466 @@
+/*
+ * 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.differentiation;
+
+import java.lang.reflect.Field;
+import java.util.HashMap;
+import java.util.Map;
+
+import org.apache.commons.math3.util.ArithmeticUtils;
+import org.junit.Assert;
+import org.junit.Test;
+
+
+/**
+ * Test for class {@link DSCompiler}.
+ */
+public class DSCompilerTest {
+
+    @Test
+    public void testSize() {
+        for (int i = 0; i < 6; ++i) {
+            for (int j = 0; j < 6; ++j) {
+                long expected = ArithmeticUtils.binomialCoefficient(i + j, i);
+                Assert.assertEquals(expected, DSCompiler.getCompiler(i, j).getSize());
+                Assert.assertEquals(expected, DSCompiler.getCompiler(j, i).getSize());
+            }
+        }
+    }
+
+    @Test
+    public void testIndices() {
+
+        DSCompiler c = DSCompiler.getCompiler(0, 0);
+        checkIndices(c.getPartialDerivativeOrders(0), new int[0]);
+
+        c = DSCompiler.getCompiler(0, 1);
+        checkIndices(c.getPartialDerivativeOrders(0), new int[0]);
+
+        c = DSCompiler.getCompiler(1, 0);
+        checkIndices(c.getPartialDerivativeOrders(0), 0);
+
+        c = DSCompiler.getCompiler(1, 1);
+        checkIndices(c.getPartialDerivativeOrders(0), 0);
+        checkIndices(c.getPartialDerivativeOrders(1), 1);
+
+        c = DSCompiler.getCompiler(1, 2);
+        checkIndices(c.getPartialDerivativeOrders(0), 0);
+        checkIndices(c.getPartialDerivativeOrders(1), 1);
+        checkIndices(c.getPartialDerivativeOrders(2), 2);
+
+        c = DSCompiler.getCompiler(2, 1);
+        checkIndices(c.getPartialDerivativeOrders(0), 0, 0);
+        checkIndices(c.getPartialDerivativeOrders(1), 1, 0);
+        checkIndices(c.getPartialDerivativeOrders(2), 0, 1);
+
+        c = DSCompiler.getCompiler(1, 3);
+        checkIndices(c.getPartialDerivativeOrders(0), 0);
+        checkIndices(c.getPartialDerivativeOrders(1), 1);
+        checkIndices(c.getPartialDerivativeOrders(2), 2);
+        checkIndices(c.getPartialDerivativeOrders(3), 3);
+
+        c = DSCompiler.getCompiler(2, 2);
+        checkIndices(c.getPartialDerivativeOrders(0), 0, 0);
+        checkIndices(c.getPartialDerivativeOrders(1), 1, 0);
+        checkIndices(c.getPartialDerivativeOrders(2), 2, 0);
+        checkIndices(c.getPartialDerivativeOrders(3), 0, 1);
+        checkIndices(c.getPartialDerivativeOrders(4), 1, 1);
+        checkIndices(c.getPartialDerivativeOrders(5), 0, 2);
+
+        c = DSCompiler.getCompiler(3, 1);
+        checkIndices(c.getPartialDerivativeOrders(0), 0, 0, 0);
+        checkIndices(c.getPartialDerivativeOrders(1), 1, 0, 0);
+        checkIndices(c.getPartialDerivativeOrders(2), 0, 1, 0);
+        checkIndices(c.getPartialDerivativeOrders(3), 0, 0, 1);
+
+        c = DSCompiler.getCompiler(1, 4);
+        checkIndices(c.getPartialDerivativeOrders(0), 0);
+        checkIndices(c.getPartialDerivativeOrders(1), 1);
+        checkIndices(c.getPartialDerivativeOrders(2), 2);
+        checkIndices(c.getPartialDerivativeOrders(3), 3);
+        checkIndices(c.getPartialDerivativeOrders(4), 4);
+
+        c = DSCompiler.getCompiler(2, 3);
+        checkIndices(c.getPartialDerivativeOrders(0), 0, 0);
+        checkIndices(c.getPartialDerivativeOrders(1), 1, 0);
+        checkIndices(c.getPartialDerivativeOrders(2), 2, 0);
+        checkIndices(c.getPartialDerivativeOrders(3), 3, 0);
+        checkIndices(c.getPartialDerivativeOrders(4), 0, 1);
+        checkIndices(c.getPartialDerivativeOrders(5), 1, 1);
+        checkIndices(c.getPartialDerivativeOrders(6), 2, 1);
+        checkIndices(c.getPartialDerivativeOrders(7), 0, 2);
+        checkIndices(c.getPartialDerivativeOrders(8), 1, 2);
+        checkIndices(c.getPartialDerivativeOrders(9), 0, 3);
+
+        c = DSCompiler.getCompiler(3, 2);
+        checkIndices(c.getPartialDerivativeOrders(0), 0, 0, 0);
+        checkIndices(c.getPartialDerivativeOrders(1), 1, 0, 0);
+        checkIndices(c.getPartialDerivativeOrders(2), 2, 0, 0);
+        checkIndices(c.getPartialDerivativeOrders(3), 0, 1, 0);
+        checkIndices(c.getPartialDerivativeOrders(4), 1, 1, 0);
+        checkIndices(c.getPartialDerivativeOrders(5), 0, 2, 0);
+        checkIndices(c.getPartialDerivativeOrders(6), 0, 0, 1);
+        checkIndices(c.getPartialDerivativeOrders(7), 1, 0, 1);
+        checkIndices(c.getPartialDerivativeOrders(8), 0, 1, 1);
+        checkIndices(c.getPartialDerivativeOrders(9), 0, 0, 2);
+
+        c = DSCompiler.getCompiler(4, 1);
+        checkIndices(c.getPartialDerivativeOrders(0), 0, 0, 0, 0);
+        checkIndices(c.getPartialDerivativeOrders(1), 1, 0, 0, 0);
+        checkIndices(c.getPartialDerivativeOrders(2), 0, 1, 0, 0);
+        checkIndices(c.getPartialDerivativeOrders(3), 0, 0, 1, 0);
+        checkIndices(c.getPartialDerivativeOrders(4), 0, 0, 0, 1);
+
+    }
+
+    @Test
+    public void testSymmetry() {
+        for (int i = 0; i < 6; ++i) {
+            for (int j = 0; j < 6; ++j) {
+                DSCompiler c = DSCompiler.getCompiler(i, j);
+                for (int k = 0; k < c.getSize(); ++k) {
+                    Assert.assertEquals(k, c.getPartialDerivativeIndex(c.getPartialDerivativeOrders(k)));
+                }
+            }
+        }
+    }
+
+    @Test public void testMultiplicationRules()
+        throws SecurityException, NoSuchFieldException, IllegalArgumentException, IllegalAccessException {
+
+        Map<String,String> referenceRules = new HashMap<String, String>();
+        referenceRules.put("(f*g)",          "f * g");
+        referenceRules.put("d(f*g)/dx",      "f * dg/dx + df/dx * g");
+        referenceRules.put("d(f*g)/dy",      referenceRules.get("d(f*g)/dx").replaceAll("x", "y"));
+        referenceRules.put("d(f*g)/dz",      referenceRules.get("d(f*g)/dx").replaceAll("x", "z"));
+        referenceRules.put("d(f*g)/dt",      referenceRules.get("d(f*g)/dx").replaceAll("x", "t"));
+        referenceRules.put("d2(f*g)/dx2",    "f * d2g/dx2 + 2 * df/dx * dg/dx + d2f/dx2 * g");
+        referenceRules.put("d2(f*g)/dy2",    referenceRules.get("d2(f*g)/dx2").replaceAll("x", "y"));
+        referenceRules.put("d2(f*g)/dz2",    referenceRules.get("d2(f*g)/dx2").replaceAll("x", "z"));
+        referenceRules.put("d2(f*g)/dt2",    referenceRules.get("d2(f*g)/dx2").replaceAll("x", "t"));
+        referenceRules.put("d2(f*g)/dxdy",   "f * d2g/dxdy + df/dy * dg/dx + df/dx * dg/dy + d2f/dxdy * g");
+        referenceRules.put("d2(f*g)/dxdz",   referenceRules.get("d2(f*g)/dxdy").replaceAll("y", "z"));
+        referenceRules.put("d2(f*g)/dxdt",   referenceRules.get("d2(f*g)/dxdy").replaceAll("y", "t"));
+        referenceRules.put("d2(f*g)/dydz",   referenceRules.get("d2(f*g)/dxdz").replaceAll("x", "y"));
+        referenceRules.put("d2(f*g)/dydt",   referenceRules.get("d2(f*g)/dxdt").replaceAll("x", "y"));
+        referenceRules.put("d2(f*g)/dzdt",   referenceRules.get("d2(f*g)/dxdt").replaceAll("x", "z"));
+        referenceRules.put("d3(f*g)/dx3",    "f * d3g/dx3 +" +
+                                             " 3 * df/dx * d2g/dx2 +" +
+                                             " 3 * d2f/dx2 * dg/dx +" +
+                                             " d3f/dx3 * g");
+        referenceRules.put("d3(f*g)/dy3",   referenceRules.get("d3(f*g)/dx3").replaceAll("x", "y"));
+        referenceRules.put("d3(f*g)/dz3",   referenceRules.get("d3(f*g)/dx3").replaceAll("x", "z"));
+        referenceRules.put("d3(f*g)/dt3",   referenceRules.get("d3(f*g)/dx3").replaceAll("x", "t"));
+        referenceRules.put("d3(f*g)/dx2dy",  "f * d3g/dx2dy +" +
+                                             " df/dy * d2g/dx2 +" +
+                                             " 2 * df/dx * d2g/dxdy +" +
+                                             " 2 * d2f/dxdy * dg/dx +" +
+                                             " d2f/dx2 * dg/dy +" +
+                                             " d3f/dx2dy * g");
+        referenceRules.put("d3(f*g)/dxdy2",  "f * d3g/dxdy2 +" +
+                                             " 2 * df/dy * d2g/dxdy +" +
+                                             " d2f/dy2 * dg/dx +" +
+                                             " df/dx * d2g/dy2 +" +
+                                             " 2 * d2f/dxdy * dg/dy +" +
+                                             " d3f/dxdy2 * g");
+        referenceRules.put("d3(f*g)/dx2dz",   referenceRules.get("d3(f*g)/dx2dy").replaceAll("y", "z"));
+        referenceRules.put("d3(f*g)/dy2dz",   referenceRules.get("d3(f*g)/dx2dz").replaceAll("x", "y"));
+        referenceRules.put("d3(f*g)/dxdz2",   referenceRules.get("d3(f*g)/dxdy2").replaceAll("y", "z"));
+        referenceRules.put("d3(f*g)/dydz2",   referenceRules.get("d3(f*g)/dxdz2").replaceAll("x", "y"));
+        referenceRules.put("d3(f*g)/dx2dt",   referenceRules.get("d3(f*g)/dx2dz").replaceAll("z", "t"));
+        referenceRules.put("d3(f*g)/dy2dt",   referenceRules.get("d3(f*g)/dx2dt").replaceAll("x", "y"));
+        referenceRules.put("d3(f*g)/dz2dt",   referenceRules.get("d3(f*g)/dx2dt").replaceAll("x", "z"));
+        referenceRules.put("d3(f*g)/dxdt2",   referenceRules.get("d3(f*g)/dxdy2").replaceAll("y", "t"));
+        referenceRules.put("d3(f*g)/dydt2",   referenceRules.get("d3(f*g)/dxdt2").replaceAll("x", "y"));
+        referenceRules.put("d3(f*g)/dzdt2",   referenceRules.get("d3(f*g)/dxdt2").replaceAll("x", "z"));
+        referenceRules.put("d3(f*g)/dxdydz", "f * d3g/dxdydz +" +
+                                             " df/dz * d2g/dxdy +" +
+                                             " df/dy * d2g/dxdz +" +
+                                             " d2f/dydz * dg/dx +" +
+                                             " df/dx * d2g/dydz +" +
+                                             " d2f/dxdz * dg/dy +" +
+                                             " d2f/dxdy * dg/dz +" +
+                                             " d3f/dxdydz * g");
+        referenceRules.put("d3(f*g)/dxdydt", referenceRules.get("d3(f*g)/dxdydz").replaceAll("z", "t"));
+        referenceRules.put("d3(f*g)/dxdzdt", referenceRules.get("d3(f*g)/dxdydt").replaceAll("y", "z"));
+        referenceRules.put("d3(f*g)/dydzdt", referenceRules.get("d3(f*g)/dxdzdt").replaceAll("x", "y"));
+
+        Field multFieldArrayField = DSCompiler.class.getDeclaredField("multIndirection");
+        multFieldArrayField.setAccessible(true);
+        for (int i = 0; i < 5; ++i) {
+            for (int j = 0; j < 4; ++j) {
+                DSCompiler compiler = DSCompiler.getCompiler(i, j);
+                int[][][] multIndirection = (int[][][]) multFieldArrayField.get(compiler);
+                for (int k = 0; k < multIndirection.length; ++k) {
+                    String product = ordersToString(compiler.getPartialDerivativeOrders(k),
+                                                    "(f*g)", "x", "y", "z", "t");
+                    StringBuilder rule = new StringBuilder();
+                    for (int[] term : multIndirection[k]) {
+                        if (rule.length() > 0) {
+                            rule.append(" + ");
+                        }
+                        if (term[0] > 1) {
+                            rule.append(term[0]).append(" * ");
+                        }
+                        rule.append(ordersToString(compiler.getPartialDerivativeOrders(term[1]),
+                                                   "f", "x", "y", "z", "t"));
+                        rule.append(" * ");
+                        rule.append(ordersToString(compiler.getPartialDerivativeOrders(term[2]),
+                                                   "g", "x", "y", "z", "t"));
+                    }
+                    Assert.assertEquals(product, referenceRules.get(product), rule.toString());
+                }
+            }
+        }
+    }
+
+    @Test public void testCompositionRules()
+        throws SecurityException, NoSuchFieldException, IllegalArgumentException, IllegalAccessException {
+
+        // the following reference rules have all been computed independently from the library,
+        // using only pencil and paper and some search and replace to handle symmetries
+        Map<String,String> referenceRules = new HashMap<String, String>();
+        referenceRules.put("(f(g))",              "(f(g))");
+        referenceRules.put("d(f(g))/dx",          "d(f(g))/dg * dg/dx");
+        referenceRules.put("d(f(g))/dy",          referenceRules.get("d(f(g))/dx").replaceAll("x", "y"));
+        referenceRules.put("d(f(g))/dz",          referenceRules.get("d(f(g))/dx").replaceAll("x", "z"));
+        referenceRules.put("d(f(g))/dt",          referenceRules.get("d(f(g))/dx").replaceAll("x", "t"));
+        referenceRules.put("d2(f(g))/dx2",        "d2(f(g))/dg2 * dg/dx * dg/dx + d(f(g))/dg * d2g/dx2");
+        referenceRules.put("d2(f(g))/dy2",        referenceRules.get("d2(f(g))/dx2").replaceAll("x", "y"));
+        referenceRules.put("d2(f(g))/dz2",        referenceRules.get("d2(f(g))/dx2").replaceAll("x", "z"));
+        referenceRules.put("d2(f(g))/dt2",        referenceRules.get("d2(f(g))/dx2").replaceAll("x", "t"));
+        referenceRules.put("d2(f(g))/dxdy",       "d2(f(g))/dg2 * dg/dx * dg/dy + d(f(g))/dg * d2g/dxdy");
+        referenceRules.put("d2(f(g))/dxdz",       referenceRules.get("d2(f(g))/dxdy").replaceAll("y", "z"));
+        referenceRules.put("d2(f(g))/dxdt",       referenceRules.get("d2(f(g))/dxdy").replaceAll("y", "t"));
+        referenceRules.put("d2(f(g))/dydz",       referenceRules.get("d2(f(g))/dxdz").replaceAll("x", "y"));
+        referenceRules.put("d2(f(g))/dydt",       referenceRules.get("d2(f(g))/dxdt").replaceAll("x", "y"));
+        referenceRules.put("d2(f(g))/dzdt",       referenceRules.get("d2(f(g))/dxdt").replaceAll("x", "z"));
+        referenceRules.put("d3(f(g))/dx3",        "d3(f(g))/dg3 * dg/dx * dg/dx * dg/dx +" +
+                                                  " 3 * d2(f(g))/dg2 * dg/dx * d2g/dx2 +" +
+                                                  " d(f(g))/dg * d3g/dx3");
+        referenceRules.put("d3(f(g))/dy3",        referenceRules.get("d3(f(g))/dx3").replaceAll("x", "y"));
+        referenceRules.put("d3(f(g))/dz3",        referenceRules.get("d3(f(g))/dx3").replaceAll("x", "z"));
+        referenceRules.put("d3(f(g))/dt3",        referenceRules.get("d3(f(g))/dx3").replaceAll("x", "t"));
+        referenceRules.put("d3(f(g))/dxdy2",      "d3(f(g))/dg3 * dg/dx * dg/dy * dg/dy +" +
+                                                  " 2 * d2(f(g))/dg2 * dg/dy * d2g/dxdy +" +
+                                                  " d2(f(g))/dg2 * dg/dx * d2g/dy2 +" +
+                                                  " d(f(g))/dg * d3g/dxdy2");
+        referenceRules.put("d3(f(g))/dxdz2",      referenceRules.get("d3(f(g))/dxdy2").replaceAll("y", "z"));
+        referenceRules.put("d3(f(g))/dxdt2",      referenceRules.get("d3(f(g))/dxdy2").replaceAll("y", "t"));
+        referenceRules.put("d3(f(g))/dydz2",      referenceRules.get("d3(f(g))/dxdz2").replaceAll("x", "y"));
+        referenceRules.put("d3(f(g))/dydt2",      referenceRules.get("d3(f(g))/dxdt2").replaceAll("x", "y"));
+        referenceRules.put("d3(f(g))/dzdt2",      referenceRules.get("d3(f(g))/dxdt2").replaceAll("x", "z"));
+        referenceRules.put("d3(f(g))/dx2dy",      "d3(f(g))/dg3 * dg/dx * dg/dx * dg/dy +" +
+                                                  " 2 * d2(f(g))/dg2 * dg/dx * d2g/dxdy +" +
+                                                  " d2(f(g))/dg2 * d2g/dx2 * dg/dy +" +
+                                                  " d(f(g))/dg * d3g/dx2dy");
+        referenceRules.put("d3(f(g))/dx2dz",      referenceRules.get("d3(f(g))/dx2dy").replaceAll("y", "z"));
+        referenceRules.put("d3(f(g))/dx2dt",      referenceRules.get("d3(f(g))/dx2dy").replaceAll("y", "t"));
+        referenceRules.put("d3(f(g))/dy2dz",      referenceRules.get("d3(f(g))/dx2dz").replaceAll("x", "y"));
+        referenceRules.put("d3(f(g))/dy2dt",      referenceRules.get("d3(f(g))/dx2dt").replaceAll("x", "y"));
+        referenceRules.put("d3(f(g))/dz2dt",      referenceRules.get("d3(f(g))/dx2dt").replaceAll("x", "z"));
+        referenceRules.put("d3(f(g))/dxdydz",     "d3(f(g))/dg3 * dg/dx * dg/dy * dg/dz +" +
+                                                  " d2(f(g))/dg2 * dg/dy * d2g/dxdz +" +
+                                                  " d2(f(g))/dg2 * dg/dx * d2g/dydz +" +
+                                                  " d2(f(g))/dg2 * d2g/dxdy * dg/dz +" +
+                                                  " d(f(g))/dg * d3g/dxdydz");
+        referenceRules.put("d3(f(g))/dxdydt",     referenceRules.get("d3(f(g))/dxdydz").replaceAll("z", "t"));
+        referenceRules.put("d3(f(g))/dxdzdt",     referenceRules.get("d3(f(g))/dxdydt").replaceAll("y", "z"));
+        referenceRules.put("d3(f(g))/dydzdt",     referenceRules.get("d3(f(g))/dxdzdt").replaceAll("x", "y"));
+        referenceRules.put("d4(f(g))/dx4",        "d4(f(g))/dg4 * dg/dx * dg/dx * dg/dx * dg/dx +" +
+                                                  " 6 * d3(f(g))/dg3 * dg/dx * dg/dx * d2g/dx2 +" +
+                                                  " 3 * d2(f(g))/dg2 * d2g/dx2 * d2g/dx2 +" +
+                                                  " 4 * d2(f(g))/dg2 * dg/dx * d3g/dx3 +" +
+                                                  " d(f(g))/dg * d4g/dx4");
+        referenceRules.put("d4(f(g))/dy4",        referenceRules.get("d4(f(g))/dx4").replaceAll("x", "y"));
+        referenceRules.put("d4(f(g))/dz4",        referenceRules.get("d4(f(g))/dx4").replaceAll("x", "z"));
+        referenceRules.put("d4(f(g))/dt4",        referenceRules.get("d4(f(g))/dx4").replaceAll("x", "t"));
+        referenceRules.put("d4(f(g))/dx3dy",      "d4(f(g))/dg4 * dg/dx * dg/dx * dg/dx * dg/dy +" +
+                                                  " 3 * d3(f(g))/dg3 * dg/dx * dg/dx * d2g/dxdy +" +
+                                                  " 3 * d3(f(g))/dg3 * dg/dx * d2g/dx2 * dg/dy +" +
+                                                  " 3 * d2(f(g))/dg2 * d2g/dx2 * d2g/dxdy +" +
+                                                  " 3 * d2(f(g))/dg2 * dg/dx * d3g/dx2dy +" +
+                                                  " d2(f(g))/dg2 * d3g/dx3 * dg/dy +" +
+                                                  " d(f(g))/dg * d4g/dx3dy");
+        referenceRules.put("d4(f(g))/dx3dz",      referenceRules.get("d4(f(g))/dx3dy").replaceAll("y", "z"));
+        referenceRules.put("d4(f(g))/dx3dt",      referenceRules.get("d4(f(g))/dx3dy").replaceAll("y", "t"));
+        referenceRules.put("d4(f(g))/dxdy3",      "d4(f(g))/dg4 * dg/dx * dg/dy * dg/dy * dg/dy +" +
+                                                  " 3 * d3(f(g))/dg3 * dg/dy * dg/dy * d2g/dxdy +" +
+                                                  " 3 * d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dy2 +" +
+                                                  " 3 * d2(f(g))/dg2 * d2g/dxdy * d2g/dy2 +" +
+                                                  " 3 * d2(f(g))/dg2 * dg/dy * d3g/dxdy2 +" +
+                                                  " d2(f(g))/dg2 * dg/dx * d3g/dy3 +" +
+                                                  " d(f(g))/dg * d4g/dxdy3");
+        referenceRules.put("d4(f(g))/dxdz3",      referenceRules.get("d4(f(g))/dxdy3").replaceAll("y", "z"));
+        referenceRules.put("d4(f(g))/dxdt3",      referenceRules.get("d4(f(g))/dxdy3").replaceAll("y", "t"));
+        referenceRules.put("d4(f(g))/dy3dz",      referenceRules.get("d4(f(g))/dx3dz").replaceAll("x", "y"));
+        referenceRules.put("d4(f(g))/dy3dt",      referenceRules.get("d4(f(g))/dx3dt").replaceAll("x", "y"));
+        referenceRules.put("d4(f(g))/dydz3",      referenceRules.get("d4(f(g))/dxdz3").replaceAll("x", "y"));
+        referenceRules.put("d4(f(g))/dydt3",      referenceRules.get("d4(f(g))/dxdt3").replaceAll("x", "y"));
+        referenceRules.put("d4(f(g))/dz3dt",      referenceRules.get("d4(f(g))/dx3dt").replaceAll("x", "z"));
+        referenceRules.put("d4(f(g))/dzdt3",      referenceRules.get("d4(f(g))/dxdt3").replaceAll("x", "z"));
+        referenceRules.put("d4(f(g))/dx2dy2",     "d4(f(g))/dg4 * dg/dx * dg/dx * dg/dy * dg/dy +" +
+                                                  " 4 * d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dxdy +" +
+                                                  " d3(f(g))/dg3 * dg/dx * dg/dx * d2g/dy2 +" +
+                                                  " 2 * d2(f(g))/dg2 * d2g/dxdy * d2g/dxdy +" +
+                                                  " 2 * d2(f(g))/dg2 * dg/dx * d3g/dxdy2 +" +
+                                                  " d3(f(g))/dg3 * d2g/dx2 * dg/dy * dg/dy +" +
+                                                  " 2 * d2(f(g))/dg2 * dg/dy * d3g/dx2dy +" +
+                                                  " d2(f(g))/dg2 * d2g/dx2 * d2g/dy2 +" +
+                                                  " d(f(g))/dg * d4g/dx2dy2");
+        referenceRules.put("d4(f(g))/dx2dz2",     referenceRules.get("d4(f(g))/dx2dy2").replaceAll("y", "z"));
+        referenceRules.put("d4(f(g))/dx2dt2",     referenceRules.get("d4(f(g))/dx2dy2").replaceAll("y", "t"));
+        referenceRules.put("d4(f(g))/dy2dz2",     referenceRules.get("d4(f(g))/dx2dz2").replaceAll("x", "y"));
+        referenceRules.put("d4(f(g))/dy2dt2",     referenceRules.get("d4(f(g))/dx2dt2").replaceAll("x", "y"));
+        referenceRules.put("d4(f(g))/dz2dt2",     referenceRules.get("d4(f(g))/dx2dt2").replaceAll("x", "z"));
+
+        referenceRules.put("d4(f(g))/dx2dydz",    "d4(f(g))/dg4 * dg/dx * dg/dx * dg/dy * dg/dz +" +
+                                                  " 2 * d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dxdz +" +
+                                                  " d3(f(g))/dg3 * dg/dx * dg/dx * d2g/dydz +" +
+                                                  " 2 * d3(f(g))/dg3 * dg/dx * d2g/dxdy * dg/dz +" +
+                                                  " 2 * d2(f(g))/dg2 * d2g/dxdy * d2g/dxdz +" +
+                                                  " 2 * d2(f(g))/dg2 * dg/dx * d3g/dxdydz +" +
+                                                  " d3(f(g))/dg3 * d2g/dx2 * dg/dy * dg/dz +" +
+                                                  " d2(f(g))/dg2 * dg/dy * d3g/dx2dz +" +
+                                                  " d2(f(g))/dg2 * d2g/dx2 * d2g/dydz +" +
+                                                  " d2(f(g))/dg2 * d3g/dx2dy * dg/dz +" +
+                                                  " d(f(g))/dg * d4g/dx2dydz");
+        referenceRules.put("d4(f(g))/dx2dydt",    referenceRules.get("d4(f(g))/dx2dydz").replaceAll("z", "t"));
+        referenceRules.put("d4(f(g))/dx2dzdt",    referenceRules.get("d4(f(g))/dx2dydt").replaceAll("y", "z"));
+        referenceRules.put("d4(f(g))/dxdy2dz",    "d4(f(g))/dg4 * dg/dx * dg/dy * dg/dy * dg/dz +" +
+                                                  " d3(f(g))/dg3 * dg/dy * dg/dy * d2g/dxdz +" +
+                                                  " 2 * d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dydz +" +
+                                                  " 2 * d3(f(g))/dg3 * dg/dy * d2g/dxdy * dg/dz +" +
+                                                  " 2 * d2(f(g))/dg2 * d2g/dxdy * d2g/dydz +" +
+                                                  " 2 * d2(f(g))/dg2 * dg/dy * d3g/dxdydz +" +
+                                                  " d3(f(g))/dg3 * dg/dx * d2g/dy2 * dg/dz +" +
+                                                  " d2(f(g))/dg2 * d2g/dy2 * d2g/dxdz +" +
+                                                  " d2(f(g))/dg2 * dg/dx * d3g/dy2dz +" +
+                                                  " d2(f(g))/dg2 * d3g/dxdy2 * dg/dz +" +
+                                                  " d(f(g))/dg * d4g/dxdy2dz");
+        referenceRules.put("d4(f(g))/dxdy2dt",    referenceRules.get("d4(f(g))/dxdy2dz").replaceAll("z", "t"));
+        referenceRules.put("d4(f(g))/dy2dzdt",    referenceRules.get("d4(f(g))/dx2dzdt").replaceAll("x", "y"));
+        referenceRules.put("d4(f(g))/dxdydz2",    "d4(f(g))/dg4 * dg/dx * dg/dy * dg/dz * dg/dz +" +
+                                                  " 2 * d3(f(g))/dg3 * dg/dy * dg/dz * d2g/dxdz +" +
+                                                  " 2 * d3(f(g))/dg3 * dg/dx * dg/dz * d2g/dydz +" +
+                                                  " d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dz2 +" +
+                                                  " 2 * d2(f(g))/dg2 * d2g/dxdz * d2g/dydz +" +
+                                                  " d2(f(g))/dg2 * dg/dy * d3g/dxdz2 +" +
+                                                  " d2(f(g))/dg2 * dg/dx * d3g/dydz2 +" +
+                                                  " d3(f(g))/dg3 * d2g/dxdy * dg/dz * dg/dz +" +
+                                                  " 2 * d2(f(g))/dg2 * dg/dz * d3g/dxdydz +" +
+                                                  " d2(f(g))/dg2 * d2g/dxdy * d2g/dz2 +" +
+                                                  " d(f(g))/dg * d4g/dxdydz2");
+        referenceRules.put("d4(f(g))/dxdz2dt",    referenceRules.get("d4(f(g))/dxdy2dt").replaceAll("y", "z"));
+        referenceRules.put("d4(f(g))/dydz2dt",    referenceRules.get("d4(f(g))/dxdz2dt").replaceAll("x", "y"));
+        referenceRules.put("d4(f(g))/dxdydt2",    referenceRules.get("d4(f(g))/dxdydz2").replaceAll("z", "t"));
+        referenceRules.put("d4(f(g))/dxdzdt2",    referenceRules.get("d4(f(g))/dxdydt2").replaceAll("y", "z"));
+        referenceRules.put("d4(f(g))/dydzdt2",    referenceRules.get("d4(f(g))/dxdzdt2").replaceAll("x", "y"));
+        referenceRules.put("d4(f(g))/dxdydzdt",   "d4(f(g))/dg4 * dg/dx * dg/dy * dg/dz * dg/dt +" +
+                                                  " d3(f(g))/dg3 * dg/dy * dg/dz * d2g/dxdt +" +
+                                                  " d3(f(g))/dg3 * dg/dx * dg/dz * d2g/dydt +" +
+                                                  " d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dzdt +" +
+                                                  " d3(f(g))/dg3 * dg/dy * d2g/dxdz * dg/dt +" +
+                                                  " d2(f(g))/dg2 * d2g/dxdz * d2g/dydt +" +
+                                                  " d2(f(g))/dg2 * dg/dy * d3g/dxdzdt +" +
+                                                  " d3(f(g))/dg3 * dg/dx * d2g/dydz * dg/dt +" +
+                                                  " d2(f(g))/dg2 * d2g/dydz * d2g/dxdt +" +
+                                                  " d2(f(g))/dg2 * dg/dx * d3g/dydzdt +" +
+                                                  " d3(f(g))/dg3 * d2g/dxdy * dg/dz * dg/dt +" +
+                                                  " d2(f(g))/dg2 * dg/dz * d3g/dxdydt +" +
+                                                  " d2(f(g))/dg2 * d2g/dxdy * d2g/dzdt +" +
+                                                  " d2(f(g))/dg2 * d3g/dxdydz * dg/dt +" +
+                                                  " d(f(g))/dg * d4g/dxdydzdt");
+
+        Field compFieldArrayField = DSCompiler.class.getDeclaredField("compIndirection");
+        compFieldArrayField.setAccessible(true);
+        for (int i = 0; i < 5; ++i) {
+            for (int j = 0; j < 5; ++j) {
+                DSCompiler compiler = DSCompiler.getCompiler(i, j);
+                int[][][] compIndirection = (int[][][]) compFieldArrayField.get(compiler);
+                for (int k = 0; k < compIndirection.length; ++k) {
+                    String product = ordersToString(compiler.getPartialDerivativeOrders(k),
+                                                    "(f(g))", "x", "y", "z", "t");
+                    StringBuilder rule = new StringBuilder();
+                    for (int[] term : compIndirection[k]) {
+                        if (rule.length() > 0) {
+                            rule.append(" + ");
+                        }
+                        if (term[0] > 1) {
+                            rule.append(term[0]).append(" * ");
+                        }
+                        rule.append(orderToString(term[1], "(f(g))", "g"));
+                        for (int l = 2; l < term.length; ++l) {
+                            rule.append(" * ");
+                            rule.append(ordersToString(compiler.getPartialDerivativeOrders(term[l]),
+                                                       "g", "x", "y", "z", "t"));
+                        }
+                    }
+                    Assert.assertEquals(product, referenceRules.get(product), rule.toString());
+                }
+            }
+        }
+    }
+
+    private void checkIndices(int[] indices, int ... expected) {
+        Assert.assertEquals(expected.length, indices.length);
+        for (int i = 0; i < expected.length; ++i) {
+            Assert.assertEquals(expected[i], indices[i]);
+        }
+    }
+
+    private String orderToString(int order, String functionName, String parameterName) {
+        if (order == 0) {
+            return functionName;
+        } else if (order == 1) {
+            return "d" + functionName + "/d" + parameterName;
+        } else {
+            return "d" + order + functionName + "/d" + parameterName + order;
+        }
+    }
+
+    private String ordersToString(int[] orders, String functionName, String ... parametersNames) {
+
+        int sumOrders = 0;
+        for (int order : orders) {
+            sumOrders += order;
+        }
+
+        if (sumOrders == 0) {
+            return functionName;
+        }
+
+        StringBuilder builder = new StringBuilder();
+        builder.append('d');
+        if (sumOrders > 1) {
+            builder.append(sumOrders);
+        }
+        builder.append(functionName).append('/');
+        for (int i = 0; i < orders.length; ++i) {
+            if (orders[i] > 0) {
+                builder.append('d').append(parametersNames[i]);
+                if (orders[i] > 1) {
+                    builder.append(orders[i]);
+                }
+            }
+        }
+        return builder.toString();
+
+    }
+
+}

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DSCompilerTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DSCompilerTest.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"

Added: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java?rev=1370951&view=auto
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java (added)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java Wed Aug  8 20:33:43 2012
@@ -0,0 +1,539 @@
+/*
+ * 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.differentiation;
+
+import java.util.Arrays;
+import java.util.List;
+
+import org.apache.commons.math3.exception.DimensionMismatchException;
+import org.apache.commons.math3.exception.NumberIsTooLargeException;
+import org.apache.commons.math3.util.ArithmeticUtils;
+import org.apache.commons.math3.util.FastMath;
+import org.junit.Assert;
+import org.junit.Test;
+
+/**
+ * Test for class {@link DerivativeStructure}.
+ */
+public class DerivativeStructureTest {
+
+    @Test(expected=NumberIsTooLargeException.class)
+    public void testWrongVariableIndex() {
+        new DerivativeStructure(3, 1, 3, 1.0);
+    }
+
+    @Test(expected=DimensionMismatchException.class)
+    public void testMissingOrders() {
+        new DerivativeStructure(3, 1, 0, 1.0).getPartialDerivative(0, 1);
+    }
+
+    @Test(expected=NumberIsTooLargeException.class)
+    public void testTooLargeOrder() {
+        new DerivativeStructure(3, 1, 0, 1.0).getPartialDerivative(1, 1, 2);
+    }
+
+    @Test
+    public void testVariableWithoutDerivative0() {
+        DerivativeStructure v = new DerivativeStructure(1, 0, 0, 1.0);
+        Assert.assertEquals(1.0, v.getValue(), 1.0e-15);
+    }
+
+    @Test(expected=NumberIsTooLargeException.class)
+    public void testVariableWithoutDerivative1() {
+        DerivativeStructure v = new DerivativeStructure(1, 0, 0, 1.0);
+        Assert.assertEquals(1.0, v.getPartialDerivative(1), 1.0e-15);
+    }
+
+    @Test
+    public void testVariable() {
+        for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
+            checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0),
+                      1.0, 1.0, 0.0, 0.0);
+            checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0),
+                      2.0, 0.0, 1.0, 0.0);
+            checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0),
+                      3.0, 0.0, 0.0, 1.0);
+        }
+    }
+
+    @Test
+    public void testConstant() {
+        for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
+            checkF0F1(new DerivativeStructure(3, maxOrder, FastMath.PI),
+                      FastMath.PI, 0.0, 0.0, 0.0);
+        }
+    }
+
+    @Test
+    public void testPrimitiveAdd() {
+        for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
+            checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0).add(5), 6.0, 1.0, 0.0, 0.0);
+            checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0).add(5), 7.0, 0.0, 1.0, 0.0);
+            checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0).add(5), 8.0, 0.0, 0.0, 1.0);
+        }
+    }
+
+    @Test
+    public void testAdd() {
+        for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
+            DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
+            DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 2.0);
+            DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 3.0);
+            DerivativeStructure xyz = x.add(y.add(z));
+            checkF0F1(xyz, x.getValue() + y.getValue() + z.getValue(), 1.0, 1.0, 1.0);
+        }
+    }
+
+    @Test
+    public void testPrimitiveSubtract() {
+        for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
+            checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0).subtract(5), -4.0, 1.0, 0.0, 0.0);
+            checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0).subtract(5), -3.0, 0.0, 1.0, 0.0);
+            checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0).subtract(5), -2.0, 0.0, 0.0, 1.0);
+        }
+    }
+
+    @Test
+    public void testSubtract() {
+        for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
+            DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
+            DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 2.0);
+            DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 3.0);
+            DerivativeStructure xyz = x.subtract(y.subtract(z));
+            checkF0F1(xyz, x.getValue() - (y.getValue() - z.getValue()), 1.0, -1.0, 1.0);
+        }
+    }
+
+    @Test
+    public void testPrimitiveMultiply() {
+        for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
+            checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0).multiply(5),  5.0, 5.0, 0.0, 0.0);
+            checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0).multiply(5), 10.0, 0.0, 5.0, 0.0);
+            checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0).multiply(5), 15.0, 0.0, 0.0, 5.0);
+        }
+    }
+
+    @Test
+    public void testMultiply() {
+        for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
+            DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
+            DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 2.0);
+            DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 3.0);
+            DerivativeStructure xyz = x.multiply(y.multiply(z));
+            for (int i = 0; i <= maxOrder; ++i) {
+                for (int j = 0; j <= maxOrder; ++j) {
+                    for (int k = 0; k <= maxOrder; ++k) {
+                        if (i + j + k <= maxOrder) {
+                            Assert.assertEquals((i == 0 ? x.getValue() : (i == 1 ? 1.0 : 0.0)) *
+                                                (j == 0 ? y.getValue() : (j == 1 ? 1.0 : 0.0)) *
+                                                (k == 0 ? z.getValue() : (k == 1 ? 1.0 : 0.0)),
+                                                xyz.getPartialDerivative(i, j, k),
+                                                1.0e-15);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    @Test
+    public void testNegate() {
+        for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
+            checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0).negate(), -1.0, -1.0, 0.0, 0.0);
+            checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0).negate(), -2.0, 0.0, -1.0, 0.0);
+            checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0).negate(), -3.0, 0.0, 0.0, -1.0);
+        }
+    }
+
+    @Test
+    public void testReciprocal() {
+        for (double x = 0.1; x < 1.2; x += 0.1) {
+            DerivativeStructure r = new DerivativeStructure(1, 6, 0, x).reciprocal();
+            Assert.assertEquals(1 / x, r.getValue(), 1.0e-15);
+            for (int i = 1; i < r.getOrder(); ++i) {
+                double expected = ArithmeticUtils.pow(-1, i) * ArithmeticUtils.factorial(i) /
+                                  FastMath.pow(x, i + 1);
+                Assert.assertEquals(expected, r.getPartialDerivative(i), 1.0e-15 * FastMath.abs(expected));
+            }
+        }
+    }
+
+    @Test
+    public void testPower() {
+        for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
+            for (int n = 0; n < 10; ++n) {
+
+                DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
+                DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 2.0);
+                DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 3.0);
+                List<DerivativeStructure> list = Arrays.asList(x, y, z,
+                                                               x.add(y).add(z),
+                                                               x.multiply(y).multiply(z));
+
+                if (n == 0) {
+                    for (DerivativeStructure ds : list) {
+                        checkEquals(ds.getField().getOne(), ds.pow(n), 1.0e-15);
+                    }
+                } else if (n == 1) {
+                    for (DerivativeStructure ds : list) {
+                        checkEquals(ds, ds.pow(n), 1.0e-15);
+                    }
+                } else {
+                    for (DerivativeStructure ds : list) {
+                        DerivativeStructure p = ds.getField().getOne();
+                        for (int i = 0; i < n; ++i) {
+                            p = p.multiply(ds);
+                        }
+                        checkEquals(p, ds.pow(n), 1.0e-15);
+                    }
+                }
+            }
+        }
+    }
+
+    @Test
+    public void testExpression() {
+        double epsilon = 2.5e-13;
+        for (double x = 0; x < 2; x += 0.2) {
+            DerivativeStructure dsX = new DerivativeStructure(3, 5, 0, x);
+            for (double y = 0; y < 2; y += 0.2) {
+                DerivativeStructure dsY = new DerivativeStructure(3, 5, 1, y);
+                for (double z = 0; z >- 2; z -= 0.2) {
+                    DerivativeStructure dsZ = new DerivativeStructure(3, 5, 2, z);
+
+                    // f(x, y, z) = x + 5 x y - 2 z + (8 z x - y)^3
+                    DerivativeStructure ds =
+                            new DerivativeStructure(1, dsX,
+                                                    5, dsX.multiply(dsY),
+                                                    -2, dsZ,
+                                                    1, new DerivativeStructure(8, dsZ.multiply(dsX),
+                                                                               -1, dsY).pow(3));
+                    double f = x + 5 * x * y - 2 * z + FastMath.pow(8 * z * x - y, 3);
+                    Assert.assertEquals(f, ds.getValue(),
+                                        FastMath.abs(epsilon * f));
+
+                    // df/dx = 1 + 5 y + 24 (8 z x - y)^2 z
+                    double dfdx = 1 + 5 * y + 24 * z * FastMath.pow(8 * z * x - y, 2);
+                    Assert.assertEquals(dfdx, ds.getPartialDerivative(1, 0, 0),
+                                        FastMath.abs(epsilon * dfdx));
+
+                    // df/dxdy = 5 + 48 z*(y - 8 z x)
+                    double dfdxdy = 5 + 48 * z * (y - 8 * z * x);
+                    Assert.assertEquals(dfdxdy, ds.getPartialDerivative(1, 1, 0),
+                                        FastMath.abs(epsilon * dfdxdy));
+
+                    // df/dxdydz = 48 (y - 16 z x)
+                    double dfdxdydz = 48 * (y - 16 * z * x);
+                    Assert.assertEquals(dfdxdydz, ds.getPartialDerivative(1, 1, 1),
+                                        FastMath.abs(epsilon * dfdxdydz));
+
+                }
+                
+            }
+        }
+    }
+
+    @Test
+    public void testCompositionOneVariableX() {
+        double epsilon = 1.0e-13;
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.1) {
+                DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
+                for (double y = 0.1; y < 1.2; y += 0.1) {
+                    DerivativeStructure dsY = new DerivativeStructure(1, maxOrder, y);
+                    DerivativeStructure f = dsX.divide(dsY).sqrt();
+                    double f0 = FastMath.sqrt(x / y);
+                    Assert.assertEquals(f0, f.getValue(), FastMath.abs(epsilon * f0));
+                    if (f.getOrder() > 0) {
+                        double f1 = 1 / (2 * FastMath.sqrt(x * y));
+                        Assert.assertEquals(f1, f.getPartialDerivative(1), FastMath.abs(epsilon * f1));
+                        if (f.getOrder() > 1) {
+                            double f2 = -f1 / (2 * x); 
+                            Assert.assertEquals(f2, f.getPartialDerivative(2), FastMath.abs(epsilon * f2));
+                            if (f.getOrder() > 2) {
+                                double f3 = (f0 + x / (2 * y * f0)) / (4 * x * x * x); 
+                                Assert.assertEquals(f3, f.getPartialDerivative(3), FastMath.abs(epsilon * f3));
+                            }
+                        }
+                    }
+                }
+            }
+        }        
+    }
+
+    @Test
+    public void testTrigo() {
+        double epsilon = 2.0e-12;
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.1) {
+                DerivativeStructure dsX = new DerivativeStructure(3, maxOrder, 0, x);
+                for (double y = 0.1; y < 1.2; y += 0.1) {
+                    DerivativeStructure dsY = new DerivativeStructure(3, maxOrder, 1, y);
+                    for (double z = 0.1; z < 1.2; z += 0.1) {
+                        DerivativeStructure dsZ = new DerivativeStructure(3, maxOrder, 2, z);
+                        DerivativeStructure f = dsX.divide(dsY.cos().add(dsZ.tan())).sin();
+                        double a = FastMath.cos(y) + FastMath.tan(z);
+                        double f0 = FastMath.sin(x / a);
+                        Assert.assertEquals(f0, f.getValue(), FastMath.abs(epsilon * f0));
+                        if (f.getOrder() > 0) {
+                            double dfdx = FastMath.cos(x / a) / a;
+                            Assert.assertEquals(dfdx, f.getPartialDerivative(1, 0, 0), FastMath.abs(epsilon * dfdx));
+                            double dfdy =  x * FastMath.sin(y) * dfdx / a;
+                            Assert.assertEquals(dfdy, f.getPartialDerivative(0, 1, 0), FastMath.abs(epsilon * dfdy));
+                            double cz = FastMath.cos(z);
+                            double cz2 = cz * cz;
+                            double dfdz = -x * dfdx / (a * cz2);
+                            Assert.assertEquals(dfdz, f.getPartialDerivative(0, 0, 1), FastMath.abs(epsilon * dfdz));
+                            if (f.getOrder() > 1) {
+                                double df2dx2 = -(f0 / (a * a));
+                                Assert.assertEquals(df2dx2, f.getPartialDerivative(2, 0, 0), FastMath.abs(epsilon * df2dx2));
+                                double df2dy2 = x * FastMath.cos(y) * dfdx / a -
+                                                x * x * FastMath.sin(y) * FastMath.sin(y) * f0 / (a * a * a * a) +
+                                                2 * FastMath.sin(y) * dfdy / a;
+                                Assert.assertEquals(df2dy2, f.getPartialDerivative(0, 2, 0), FastMath.abs(epsilon * df2dy2));
+                                double c4 = cz2 * cz2;
+                                double df2dz2 = x * (2 * a * (1 - a * cz * FastMath.sin(z)) * dfdx - x * f0 / a ) / (a * a * a * c4);
+                                Assert.assertEquals(df2dz2, f.getPartialDerivative(0, 0, 2), FastMath.abs(epsilon * df2dz2));
+                                double df2dxdy = dfdy / x  - x * FastMath.sin(y) * f0 / (a * a * a);
+                                Assert.assertEquals(df2dxdy, f.getPartialDerivative(1, 1, 0), FastMath.abs(epsilon * df2dxdy));
+                            }
+                        }
+                    }
+                }
+            }        
+        }
+    }
+
+    @Test
+    public void testSqrtDefinition() {
+        double[] epsilon = new double[] { 5.0e-16, 5.0e-16, 2.0e-15, 5.0e-14, 2.0e-12 };
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.001) {
+                DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
+                DerivativeStructure sqrt1 = dsX.pow(0.5);
+                DerivativeStructure sqrt2 = dsX.sqrt();
+                DerivativeStructure zero = sqrt1.subtract(sqrt2);
+                for (int n = 0; n <= maxOrder; ++n) {
+                    Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
+                }
+            }
+        }
+    }
+
+    @Test
+    public void testSqrtPow2() {
+        double[] epsilon = new double[] { 1.0e-16, 3.0e-16, 2.0e-15, 6.0e-14, 6.0e-12 };
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.001) {
+                DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
+                DerivativeStructure rebuiltX = dsX.multiply(dsX).sqrt();
+                DerivativeStructure zero = rebuiltX.subtract(dsX);
+                for (int n = 0; n <= maxOrder; ++n) {
+                    Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
+                }
+            }
+        }
+    }
+
+    @Test
+    public void testCbrtDefinition() {
+        double[] epsilon = new double[] { 4.0e-16, 9.0e-16, 6.0e-15, 2.0e-13, 4.0e-12 };
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.001) {
+                DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
+                DerivativeStructure cbrt1 = dsX.pow(1.0 / 3.0);
+                DerivativeStructure cbrt2 = dsX.cbrt();
+                DerivativeStructure zero = cbrt1.subtract(cbrt2);
+                for (int n = 0; n <= maxOrder; ++n) {
+                    Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
+                }
+            }
+        }
+    }
+
+    @Test
+    public void testCbrtPow3() {
+        double[] epsilon = new double[] { 1.0e-16, 5.0e-16, 8.0e-15, 3.0e-13, 4.0e-11 };
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.001) {
+                DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
+                DerivativeStructure rebuiltX = dsX.multiply(dsX.multiply(dsX)).cbrt();
+                DerivativeStructure zero = rebuiltX.subtract(dsX);
+                for (int n = 0; n <= maxOrder; ++n) {
+                    Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
+                }
+            }
+        }
+    }
+
+    @Test
+    public void testExp() {
+        double[] epsilon = new double[] { 1.0e-16, 1.0e-16, 1.0e-16, 1.0e-16, 1.0e-16 };
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.001) {
+                double refExp = FastMath.exp(x);
+                DerivativeStructure exp = new DerivativeStructure(1, maxOrder, 0, x).exp();
+                for (int n = 0; n <= maxOrder; ++n) {
+                    Assert.assertEquals(refExp, exp.getPartialDerivative(n), epsilon[n]);
+                }
+            }
+        }
+    }
+
+    @Test
+    public void testLog() {
+        double[] epsilon = new double[] { 1.0e-16, 1.0e-16, 3.0e-14, 7.0e-13, 3.0e-11 };
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.001) {
+                DerivativeStructure log = new DerivativeStructure(1, maxOrder, 0, x).log();
+                Assert.assertEquals(FastMath.log(x), log.getValue(), epsilon[0]);
+                for (int n = 1; n <= maxOrder; ++n) {
+                    double refDer = -ArithmeticUtils.factorial(n - 1) / FastMath.pow(-x, n);
+                    Assert.assertEquals(refDer, log.getPartialDerivative(n), epsilon[n]);
+                }
+            }
+        }
+    }
+
+    @Test
+    public void testLogExp() {
+        double[] epsilon = new double[] { 2.0e-16, 2.0e-16, 3.0e-16, 2.0e-15, 6.0e-15 };
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.001) {
+                DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
+                DerivativeStructure rebuiltX = dsX.exp().log();
+                DerivativeStructure zero = rebuiltX.subtract(dsX);
+                for (int n = 0; n <= maxOrder; ++n) {
+                    Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
+                }
+            }
+        }
+    }
+
+    @Test
+    public void testTangentDefinition() {
+        double[] epsilon = new double[] { 5.0e-16, 2.0e-15, 3.0e-14, 5.0e-13, 2.0e-11 };
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.001) {
+                DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
+                DerivativeStructure tan1 = dsX.sin().divide(dsX.cos());
+                DerivativeStructure tan2 = dsX.tan();
+                DerivativeStructure zero = tan1.subtract(tan2);
+                for (int n = 0; n <= maxOrder; ++n) {
+                    Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
+                }
+            }
+        }
+    }
+
+    @Test
+    public void testCompositionOneVariableY() {
+        double epsilon = 1.0e-13;
+        for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
+            for (double x = 0.1; x < 1.2; x += 0.1) {
+                DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, x);
+                for (double y = 0.1; y < 1.2; y += 0.1) {
+                    DerivativeStructure dsY = new DerivativeStructure(1, maxOrder, 0, y);
+                    DerivativeStructure f = dsX.divide(dsY).sqrt();
+                    double f0 = FastMath.sqrt(x / y);
+                    Assert.assertEquals(f0, f.getValue(), FastMath.abs(epsilon * f0));
+                    if (f.getOrder() > 0) {
+                        double f1 = -x / (2 * y * y * f0);
+                        Assert.assertEquals(f1, f.getPartialDerivative(1), FastMath.abs(epsilon * f1));
+                        if (f.getOrder() > 1) {
+                            double f2 = (f0 - x / (4 * y * f0)) / (y * y); 
+                            Assert.assertEquals(f2, f.getPartialDerivative(2), FastMath.abs(epsilon * f2));
+                            if (f.getOrder() > 2) {
+                                double f3 = (x / (8 * y * f0) - 2 * f0) / (y * y * y); 
+                                Assert.assertEquals(f3, f.getPartialDerivative(3), FastMath.abs(epsilon * f3));
+                            }
+                        }
+                    }
+                }
+            }
+        }        
+    }
+
+    @Test
+    public void testField() {
+        for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
+            DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
+            checkF0F1(x.getField().getZero(), 0.0, 0.0, 0.0, 0.0);
+            checkF0F1(x.getField().getOne(), 1.0, 0.0, 0.0, 0.0);
+            Assert.assertEquals(maxOrder, x.getField().getZero().getOrder());
+            Assert.assertEquals(3, x.getField().getZero().getFreeParameters());
+            Assert.assertEquals(DerivativeStructure.class, x.getField().getRuntimeClass());
+        }
+    }
+
+    private void checkF0F1(DerivativeStructure ds, double value, double...derivatives) {
+
+        // check dimension
+        Assert.assertEquals(derivatives.length, ds.getFreeParameters());
+
+        // check value, directly and also as 0th order derivative
+        Assert.assertEquals(value, ds.getValue(), 1.0e-15);
+        Assert.assertEquals(value, ds.getPartialDerivative(new int[ds.getFreeParameters()]), 1.0e-15);
+
+        // check first order derivatives
+        for (int i = 0; i < derivatives.length; ++i) {
+            int[] orders = new int[derivatives.length];
+            orders[i] = 1;
+            Assert.assertEquals(derivatives[i], ds.getPartialDerivative(orders), 1.0e-15);
+        }
+
+    }
+
+    private void checkEquals(DerivativeStructure ds1, DerivativeStructure ds2, double epsilon) {
+
+        // check dimension
+        Assert.assertEquals(ds1.getFreeParameters(), ds2.getFreeParameters());
+        Assert.assertEquals(ds1.getOrder(), ds2.getOrder());
+
+        int[] derivatives = new int[ds1.getFreeParameters()];
+        int sum = 0;
+        while (true) {
+
+            if (sum <= ds1.getOrder()) {
+                Assert.assertEquals(ds1.getPartialDerivative(derivatives),
+                                    ds2.getPartialDerivative(derivatives),
+                                    epsilon);
+            }
+
+            boolean increment = true;
+            sum = 0;
+            for (int i = derivatives.length - 1; i >= 0; --i) {
+                if (increment) {
+                    if (derivatives[i] == ds1.getOrder()) {
+                        derivatives[i] = 0;
+                    } else {
+                        derivatives[i]++;
+                        increment = false;
+                    }
+                }
+                sum += derivatives[i];
+            }
+            if (increment) {
+                return;
+            }
+
+        }
+
+    }
+
+}

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java
------------------------------------------------------------------------------
    svn:keywords = "Author Date Id Revision"



Mime
View raw message