commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r919829 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java
Date Sat, 06 Mar 2010 19:39:18 GMT
Author: luc
Date: Sat Mar  6 19:39:17 2010
New Revision: 919829

URL: http://svn.apache.org/viewvc?rev=919829&view=rev
Log:
added miscellaneous methods to get some feedback on the integration process
(step start, step size, number of evaluations ...)

Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java?rev=919829&r1=919828&r2=919829&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java
(original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobians.java
Sat Mar  6 19:39:17 2010
@@ -25,6 +25,7 @@
 import java.util.Collection;
 
 import org.apache.commons.math.MathRuntimeException;
+import org.apache.commons.math.MaxEvaluationsExceededException;
 import org.apache.commons.math.ode.DerivativeException;
 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.FirstOrderIntegrator;
@@ -53,6 +54,12 @@
     /** Raw equations to integrate. */
     private final ParameterizedODEWithJacobians ode;
 
+    /** Maximal number of evaluations allowed. */
+    private int maxEvaluations;
+
+    /** Number of evaluations already performed. */
+    private int evaluations;
+
     /** Build an enhanced integrator using internal differentiation to compute jacobians.
      * @param integrator underlying integrator to solve the compound problem
      * @param ode original problem (f in the equation y' = f(t, y))
@@ -74,6 +81,7 @@
         checkDimension(ode.getParametersDimension(), hP);
         this.integrator = integrator;
         this.ode = new FiniteDifferencesWrapper(ode, p, hY, hP);
+        setMaxEvaluations(-1);
     }
 
     /** Build an enhanced integrator using ODE builtin jacobian computation features.
@@ -86,6 +94,7 @@
                                              final ParameterizedODEWithJacobians ode) {
         this.integrator = integrator;
         this.ode = ode;
+        setMaxEvaluations(-1);
     }
 
     /** Add a step handler to this integrator.
@@ -204,7 +213,8 @@
         }
 
         // integrate the compound state variational equations
-        final double stopTime = integrator.integrate(new MappingWrapper(ode), t0, z, t, z);
+        evaluations = 0;
+        final double stopTime = integrator.integrate(new MappingWrapper(), t0, z, t, z);
 
         // dispatch the final compound state into the state and partial derivatives arrays
         System.arraycopy(z, 0, y, 0, n);
@@ -219,6 +229,62 @@
 
     }
 
+    /** Get the current value of the step start time t<sub>i</sub>.
+     * <p>This method can be called during integration (typically by
+     * the object implementing the {@link FirstOrderDifferentialEquations
+     * differential equations} problem) if the value of the current step that
+     * is attempted is needed.</p>
+     * <p>The result is undefined if the method is called outside of
+     * calls to <code>integrate</code>.</p>
+     * @return current value of the step start time t<sub>i</sub>
+     */
+    public double getCurrentStepStart() {
+        return integrator.getCurrentStepStart();
+    }
+
+    /** Get the current signed value of the integration stepsize.
+     * <p>This method can be called during integration (typically by
+     * the object implementing the {@link FirstOrderDifferentialEquations
+     * differential equations} problem) if the signed value of the current stepsize
+     * that is tried is needed.</p>
+     * <p>The result is undefined if the method is called outside of
+     * calls to <code>integrate</code>.</p>
+     * @return current signed value of the stepsize
+     */
+    public double getCurrentSignedStepsize() {
+        return integrator.getCurrentSignedStepsize();
+    }
+
+    /** Set the maximal number of differential equations function evaluations.
+     * <p>The purpose of this method is to avoid infinite loops which can occur
+     * for example when stringent error constraints are set or when lots of
+     * discrete events are triggered, thus leading to many rejected steps.</p>
+     * @param maxEvaluations maximal number of function evaluations (negative
+     * values are silently converted to maximal integer value, thus representing
+     * almost unlimited evaluations)
+     */
+    public void setMaxEvaluations(int maxEvaluations) {
+        this.maxEvaluations = (maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations;
+    }
+
+    /** Get the maximal number of functions evaluations.
+     * @return maximal number of functions evaluations
+     */
+    public int getMaxEvaluations() {
+        return maxEvaluations;
+    }
+
+    /** Get the number of evaluations of the differential equations function.
+     * <p>
+     * The number of evaluations corresponds to the last call to the
+     * <code>integrate</code> method. It is 0 if the method has not been called
yet.
+     * </p>
+     * @return number of evaluations of the differential equations function
+     */
+    public int getEvaluations() {
+        return evaluations;
+    }
+
     /** Check array dimensions.
      * @param expected expected dimension
      * @param array (may be null if expected is 0)
@@ -234,10 +300,7 @@
     }
 
     /** Wrapper class used to map state and jacobians into compound state. */
-    private static class MappingWrapper implements  FirstOrderDifferentialEquations {
-
-        /** Underlying ODE with jacobians. */
-        private final ParameterizedODEWithJacobians ode;
+    private class MappingWrapper implements  FirstOrderDifferentialEquations {
 
         /** Current state. */
         private final double[]   y;
@@ -252,11 +315,8 @@
         private final double[][] dFdP;
 
         /** Simple constructor.
-         * @param ode underlying ODE with jacobians
          */
-        public MappingWrapper(final ParameterizedODEWithJacobians ode) {
-
-            this.ode = ode;
+        public MappingWrapper() {
 
             final int n = ode.getDimension();
             final int k = ode.getParametersDimension();
@@ -283,6 +343,9 @@
 
             // compute raw ODE and its jacobians: dy/dt, d[dy/dt]/dy0 and d[dy/dt]/dp
             System.arraycopy(z,    0, y,    0, n);
+            if (++evaluations > maxEvaluations) {
+                throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations));
+            }
             ode.computeDerivatives(t, y, yDot);
             ode.computeJacobians(t, y, yDot, dFdY, dFdP);
 
@@ -325,7 +388,7 @@
     }
 
     /** Wrapper class to compute jacobians by finite differences for ODE which do not compute
them themselves. */
-    private static class FiniteDifferencesWrapper
+    private class FiniteDifferencesWrapper
         implements ParameterizedODEWithJacobians {
 
         /** Raw ODE without jacobians computation. */
@@ -365,6 +428,8 @@
 
         /** {@inheritDoc} */
         public void computeDerivatives(double t, double[] y, double[] yDot) throws DerivativeException
{
+            // this call to computeDerivatives has already been counted,
+            // we must not increment the counter again
             ode.computeDerivatives(t, y, yDot);
         }
 
@@ -386,6 +451,11 @@
             final int n = ode.getDimension();
             final int k = ode.getParametersDimension();
 
+            evaluations += n + k;
+            if (evaluations > maxEvaluations) {
+                throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations));
+            }
+
             // compute df/dy where f is the ODE and y is the state array
             for (int j = 0; j < n; ++j) {
                 final double savedYj = y[j];

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java?rev=919829&r1=919828&r2=919829&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java
(original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/jacobians/FirstOrderIntegratorWithJacobiansTest.java
Sat Mar  6 19:39:17 2010
@@ -106,7 +106,12 @@
             FirstOrderIntegratorWithJacobians extInt =
                 new FirstOrderIntegratorWithJacobians(integ, brusselator, new double[] {
b },
                                                       new double[] { hY, hY }, new double[]
{ hP });
+            extInt.setMaxEvaluations(5000);
             extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 } }, 20.0, z, dZdZ0, dZdP);
+            Assert.assertEquals(5000, extInt.getMaxEvaluations());
+            Assert.assertTrue(extInt.getEvaluations() > 2000);
+            Assert.assertTrue(extInt.getEvaluations() < 2500);
+            Assert.assertEquals(4 * integ.getEvaluations(), extInt.getEvaluations());
             residualsP0.addValue(dZdP[0][0] - brusselator.dYdP0());
             residualsP1.addValue(dZdP[1][0] - brusselator.dYdP1());
         }
@@ -120,7 +125,7 @@
     public void testAnalyticalDifferentiation()
         throws IntegratorException, DerivativeException {
         FirstOrderIntegrator integ =
-            new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
+            new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-4, 1.0e-4);
         SummaryStatistics residualsP0 = new SummaryStatistics();
         SummaryStatistics residualsP1 = new SummaryStatistics();
         for (double b = 2.88; b < 3.08; b += 0.001) {
@@ -131,7 +136,12 @@
             double[][] dZdP  = new double[2][1];
             FirstOrderIntegratorWithJacobians extInt =
                 new FirstOrderIntegratorWithJacobians(integ, brusselator);
+            extInt.setMaxEvaluations(5000);
             extInt.integrate(0, z, new double[][] { { 0.0 }, { 1.0 } }, 20.0, z, dZdZ0, dZdP);
+            Assert.assertEquals(5000, extInt.getMaxEvaluations());
+            Assert.assertTrue(extInt.getEvaluations() > 510);
+            Assert.assertTrue(extInt.getEvaluations() < 610);
+            Assert.assertEquals(integ.getEvaluations(), extInt.getEvaluations());
             residualsP0.addValue(dZdP[0][0] - brusselator.dYdP0());
             residualsP1.addValue(dZdP[1][0] - brusselator.dYdP1());
         }
@@ -177,7 +187,7 @@
         double[][] dydy0 = new double[2][2];
         double[][] dydp  = new double[2][3];
         double t = 18 * Math.PI;
-        FirstOrderIntegratorWithJacobians extInt =
+        final FirstOrderIntegratorWithJacobians extInt =
             new FirstOrderIntegratorWithJacobians(integ, circle);
         extInt.addStepHandler(new StepHandlerWithJacobians() {
             
@@ -194,6 +204,8 @@
                 double[]   y     = interpolator.getInterpolatedY();
                 double[][] dydy0 = interpolator.getInterpolatedDyDy0();
                 double[][] dydp  = interpolator.getInterpolatedDyDp();
+                Assert.assertEquals(interpolator.getPreviousTime(), extInt.getCurrentStepStart(),
1.0e-10);
+                Assert.assertTrue(extInt.getCurrentSignedStepsize() < 0.5);
                 for (int i = 0; i < y.length; ++i) {
                     Assert.assertEquals(circle.exactY(t)[i], y[i], 1.0e-10);
                 }



Mime
View raw message