commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r1061508 [1/2] - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/ode/ main/java/org/apache/commons/math/ode/events/ main/java/org/apache/commons/math/ode/nonstiff/ main/java/org/apache/commons/math/ode/sampling/ site/xd...
Date Thu, 20 Jan 2011 20:57:13 GMT
Author: luc
Date: Thu Jan 20 20:57:11 2011
New Revision: 1061508

URL: http://svn.apache.org/viewvc?rev=1061508&view=rev
Log:
separate discrete event detection from adaptive step size handling in ODE solvers,
thus improving robustness, maintainability and speed
JIRA: MATH-484

Removed:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/CombinedEventsManager.java
Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ODEIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/EventState.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblem4.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblemAbstract.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblemHandler.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaIntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ClassicalRungeKuttaStepInterpolatorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince54IntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince54StepInterpolatorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853IntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolatorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/EulerIntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GillIntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GillStepInterpolatorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolatorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54IntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/HighamHall54StepInterpolatorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/MidpointIntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/MidpointStepInterpolatorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesIntegratorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/ThreeEighthesStepInterpolatorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/sampling/DummyStepInterpolatorTest.java
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolatorTest.java

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/AbstractIntegrator.java Thu Jan 20 20:57:11 2011
@@ -20,15 +20,21 @@ package org.apache.commons.math.ode;
 import java.util.ArrayList;
 import java.util.Collection;
 import java.util.Collections;
+import java.util.List;
+import java.util.SortedSet;
+import java.util.TreeSet;
 
+import org.apache.commons.math.ConvergenceException;
 import org.apache.commons.math.MaxEvaluationsExceededException;
 import org.apache.commons.math.exception.MathUserException;
 import org.apache.commons.math.exception.util.LocalizedFormats;
-import org.apache.commons.math.ode.events.CombinedEventsManager;
+import org.apache.commons.math.ode.events.EventException;
 import org.apache.commons.math.ode.events.EventHandler;
 import org.apache.commons.math.ode.events.EventState;
+import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
 import org.apache.commons.math.ode.sampling.StepHandler;
 import org.apache.commons.math.util.FastMath;
+import org.apache.commons.math.util.MathUtils;
 
 /**
  * Base class managing common boilerplate for all integrators.
@@ -46,8 +52,17 @@ public abstract class AbstractIntegrator
     /** Current stepsize. */
     protected double stepSize;
 
-    /** Events handlers manager. */
-    protected CombinedEventsManager eventsHandlersManager;
+    /** Indicator for last step. */
+    protected boolean isLastStep;
+
+    /** Indicator that a state or derivative reset was triggered by some event. */
+    protected boolean resetOccurred;
+
+    /** Events states. */
+    protected Collection<EventState> eventsStates;
+
+    /** Initialization indicator of events states. */
+    protected boolean statesInitialized;
 
     /** Name of the method. */
     private final String name;
@@ -69,7 +84,8 @@ public abstract class AbstractIntegrator
         stepHandlers = new ArrayList<StepHandler>();
         stepStart = Double.NaN;
         stepSize  = Double.NaN;
-        eventsHandlersManager = new CombinedEventsManager();
+        eventsStates = new ArrayList<EventState>();
+        statesInitialized = false;
         setMaxEvaluations(-1);
         resetEvaluations();
     }
@@ -101,22 +117,25 @@ public abstract class AbstractIntegrator
     }
 
     /** {@inheritDoc} */
-    public void addEventHandler(final EventHandler function,
+    public void addEventHandler(final EventHandler handler,
                                 final double maxCheckInterval,
                                 final double convergence,
                                 final int maxIterationCount) {
-        eventsHandlersManager.addEventHandler(function, maxCheckInterval,
-                                              convergence, maxIterationCount);
+        eventsStates.add(new EventState(handler, maxCheckInterval, convergence, maxIterationCount));
     }
 
     /** {@inheritDoc} */
     public Collection<EventHandler> getEventHandlers() {
-        return eventsHandlersManager.getEventsHandlers();
+        final List<EventHandler> list = new ArrayList<EventHandler>();
+        for (EventState state : eventsStates) {
+            list.add(state.getEventHandler());
+        }
+        return Collections.unmodifiableCollection(list);
     }
 
     /** {@inheritDoc} */
     public void clearEventHandlers() {
-        eventsHandlersManager.clearEventsHandlers();
+        eventsStates.clear();
     }
 
     /** Check if one of the step handlers requires dense output.
@@ -185,6 +204,109 @@ public abstract class AbstractIntegrator
         equations.computeDerivatives(t, y, yDot);
     }
 
+    /** Accept a step, triggering events and step handlers.
+     * @param interpolator step interpolator
+     * @param handlers step handlers
+     * @param y state vector at step end time, must be reset if an event
+     * asks for resetting or if an events stops integration during the step
+     * @param yDot placeholder array where to put the time derivative of the state vector
+     * @param tEnd final integration time
+     * @return time at end of step
+     * @exception IntegratorException if the value of one event state cannot be evaluated
+     */
+    protected double acceptStep(final AbstractStepInterpolator interpolator,
+                                final Collection<StepHandler> handlers,
+                                final double[] y,
+                                final double[] yDot, final double tEnd)
+        throws IntegratorException {
+
+        try {
+            double previousT = interpolator.getGlobalPreviousTime();
+            final double currentT = interpolator.getGlobalCurrentTime();
+            resetOccurred = false;
+
+            // initialize the events states if needed
+            if (! statesInitialized) {
+                for (EventState state : eventsStates) {
+                    state.reinitializeBegin(interpolator);
+                }
+                statesInitialized = true;
+            }
+
+            // find all events that occur during the step
+            SortedSet<EventState> occuringEvents = new TreeSet<EventState>();
+            for (final EventState state : eventsStates) {
+                if (state.evaluateStep(interpolator)) {
+                    // the event occurs during the current step
+                    occuringEvents.add(state);
+                }
+            }
+
+            // handle the events chronologically
+            for (final EventState state : occuringEvents) {
+
+                // restrict the interpolator to the first part of the step, up to the event
+                final double eventT = state.getEventTime();
+                interpolator.setSoftBounds(previousT, eventT);
+
+                // trigger the event
+                interpolator.setInterpolatedTime(eventT);
+                final double[] eventY = interpolator.getInterpolatedState();
+                state.stepAccepted(eventT, eventY);
+                isLastStep = state.stop();
+
+                // handle the first part of the step, up to the event
+                for (final StepHandler handler : stepHandlers) {
+                    handler.handleStep(interpolator, isLastStep);
+                }
+
+                if (isLastStep) {
+                    // the event asked to stop integration
+                    System.arraycopy(eventY, 0, y, 0, y.length);
+                    return eventT;
+                }
+
+                if (state.reset(eventT, eventY)) {
+                    // some event handler has triggered changes that
+                    // invalidate the derivatives, we need to recompute them
+                    System.arraycopy(eventY, 0, y, 0, y.length);
+                    computeDerivatives(eventT, y, yDot);
+                    resetOccurred = true;
+                    return eventT;
+                }
+
+                // prepare handling of the remaining part of the step
+                previousT = eventT;
+                interpolator.setSoftBounds(eventT, currentT);
+
+            }
+
+            interpolator.setInterpolatedTime(currentT);
+            final double[] currentY = interpolator.getInterpolatedState();
+            for (final EventState state : eventsStates) {
+                state.stepAccepted(currentT, currentY);
+                isLastStep = isLastStep || state.stop();
+            }
+            isLastStep = isLastStep || MathUtils.equals(currentT, tEnd, 1);
+
+            // handle the remaining part of the step, after all events if any
+            for (StepHandler handler : stepHandlers) {
+                handler.handleStep(interpolator, isLastStep);
+            }
+
+            return currentT;
+        } catch (EventException se) {
+            final Throwable cause = se.getCause();
+            if ((cause != null) && (cause instanceof MathUserException)) {
+                throw (MathUserException) cause;
+            }
+            throw new IntegratorException(se);
+        } catch (ConvergenceException ce) {
+            throw new IntegratorException(ce);
+        }
+
+    }
+
     /** Perform some sanity checks on the integration parameters.
      * @param ode differential equations set
      * @param t0 start time
@@ -216,60 +338,4 @@ public abstract class AbstractIntegrator
 
     }
 
-    /** Add an event handler for end time checking.
-     * <p>This method can be used to simplify handling of integration end time.
-     * It leverages the nominal stop condition with the exceptional stop
-     * conditions.</p>
-     * @param startTime integration start time
-     * @param endTime desired end time
-     * @param manager manager containing the user-defined handlers
-     * @return a new manager containing all the user-defined handlers plus a
-     * dedicated manager triggering a stop event at entTime
-     */
-    protected CombinedEventsManager addEndTimeChecker(final double startTime,
-                                                      final double endTime,
-                                                      final CombinedEventsManager manager) {
-        CombinedEventsManager newManager = new CombinedEventsManager();
-        for (final EventState state : manager.getEventsStates()) {
-            newManager.addEventHandler(state.getEventHandler(),
-                                       state.getMaxCheckInterval(),
-                                       state.getConvergence(),
-                                       state.getMaxIterationCount());
-        }
-        newManager.addEventHandler(new EndTimeChecker(endTime),
-                                   Double.POSITIVE_INFINITY,
-                                   FastMath.ulp(FastMath.max(FastMath.abs(startTime), FastMath.abs(endTime))),
-                                   100);
-        return newManager;
-    }
-
-    /** Specialized event handler to stop integration. */
-    private static class EndTimeChecker implements EventHandler {
-
-        /** Desired end time. */
-        private final double endTime;
-
-        /** Build an instance.
-         * @param endTime desired time
-         */
-        public EndTimeChecker(final double endTime) {
-            this.endTime = endTime;
-        }
-
-        /** {@inheritDoc} */
-        public int eventOccurred(double t, double[] y, boolean increasing) {
-            return STOP;
-        }
-
-        /** {@inheritDoc} */
-        public double g(double t, double[] y) {
-            return t - endTime;
-        }
-
-        /** {@inheritDoc} */
-        public void resetState(double t, double[] y) {
-        }
-
-    }
-
 }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ODEIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ODEIntegrator.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ODEIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/ODEIntegrator.java Thu Jan 20 20:57:11 2011
@@ -73,10 +73,8 @@ public interface ODEIntegrator  {
      * @see #getEventHandlers()
      * @see #clearEventHandlers()
      */
-    void addEventHandler(EventHandler handler,
-                                         double maxCheckInterval,
-                                         double convergence,
-                                         int maxIterationCount);
+    void addEventHandler(EventHandler handler, double maxCheckInterval,
+                         double convergence, int maxIterationCount);
 
     /** Get all the event handlers that have been added to the integrator.
      * @return an unmodifiable collection of the added events handlers

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/EventState.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/EventState.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/EventState.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/events/EventState.java Thu Jan 20 20:57:11 2011
@@ -18,10 +18,10 @@
 package org.apache.commons.math.ode.events;
 
 import org.apache.commons.math.ConvergenceException;
-import org.apache.commons.math.exception.MathUserException;
-import org.apache.commons.math.MathRuntimeException;
 import org.apache.commons.math.analysis.UnivariateRealFunction;
 import org.apache.commons.math.analysis.solvers.BrentSolver;
+import org.apache.commons.math.exception.MathInternalError;
+import org.apache.commons.math.exception.MathUserException;
 import org.apache.commons.math.ode.sampling.StepInterpolator;
 import org.apache.commons.math.util.FastMath;
 
@@ -33,13 +33,12 @@ import org.apache.commons.math.util.Fast
  * of one handler during one integration step, with references to the
  * state at the end of the preceding step. This information is used to
  * decide if the handler should trigger an event or not during the
- * proposed step (and hence the step should be reduced to ensure the
- * event occurs at a bound rather than inside the step).</p>
+ * proposed step.</p>
  *
  * @version $Revision$ $Date$
  * @since 1.2
  */
-public class EventState {
+public class EventState implements Comparable<EventState> {
 
     /** Event handler. */
     private final EventHandler handler;
@@ -187,8 +186,7 @@ public class EventState {
     /** Evaluate the impact of the proposed step on the event handler.
      * @param interpolator step interpolator for the proposed step
      * @return true if the event handler triggers an event before
-     * the end of the proposed step (this implies the step should be
-     * rejected)
+     * the end of the proposed step
      * @exception MathUserException if the interpolator fails to
      * compute the switching function somewhere within the step
      * @exception EventException if the switching function
@@ -202,16 +200,21 @@ public class EventState {
 
             forward = interpolator.isForward();
             final double t1 = interpolator.getCurrentTime();
-            final int    n  = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(t1 - t0) / maxCheckInterval));
-            final double h  = (t1 - t0) / n;
+            if (FastMath.abs(t1 - t0) < convergence) {
+                // we cannot do anything on such a small step, don't trigger any events
+                return false;
+            }
+            final double start = forward ? (t0 + convergence) : t0 - convergence;
+            final double dt    = t1 - start;
+            final int    n     = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / maxCheckInterval));
+            final double h     = dt / n;
 
             double ta = t0;
             double ga = g0;
-            double tb = t0 + (interpolator.isForward() ? convergence : -convergence);
             for (int i = 0; i < n; ++i) {
 
                 // evaluate handler value at the end of the substep
-                tb += h;
+                final double tb = start + (i + 1) * h;
                 interpolator.setInterpolatedTime(tb);
                 final double gb = handler.g(tb, interpolator.getInterpolatedState());
 
@@ -219,7 +222,22 @@ public class EventState {
                 if (g0Positive ^ (gb >= 0)) {
                     // there is a sign change: an event is expected during this step
 
-                    if (ga * gb > 0) {
+                    // variation direction, with respect to the integration direction
+                    increasing = gb >= ga;
+
+                    final UnivariateRealFunction f = new UnivariateRealFunction() {
+                        public double value(final double t) throws MathUserException {
+                            try {
+                                interpolator.setInterpolatedTime(t);
+                                return handler.g(t, interpolator.getInterpolatedState());
+                            } catch (EventException e) {
+                                throw new MathUserException(e);
+                            }
+                        }
+                    };
+                    final BrentSolver solver = new BrentSolver(convergence);
+
+                    if (ga * gb >= 0) {
                         // this is a corner case:
                         // - there was an event near ta,
                         // - there is another event between ta and tb
@@ -230,51 +248,33 @@ public class EventState {
                         final double epsilon = (forward ? 0.25 : -0.25) * convergence;
                         for (int k = 0; (k < 4) && (ga * gb > 0); ++k) {
                             ta += epsilon;
-                            interpolator.setInterpolatedTime(ta);
-                            ga = handler.g(ta, interpolator.getInterpolatedState());
+                            ga = f.value(ta);
                         }
                         if (ga * gb > 0) {
                             // this should never happen
-                            throw MathRuntimeException.createInternalError(null);
+                            throw new MathInternalError();
                         }
                     }
 
-                    // variation direction, with respect to the integration direction
-                    increasing = gb >= ga;
-
-                    final UnivariateRealFunction f = new UnivariateRealFunction() {
-                        public double value(final double t) throws MathUserException {
-                            try {
-                                interpolator.setInterpolatedTime(t);
-                                return handler.g(t, interpolator.getInterpolatedState());
-                            } catch (EventException e) {
-                                throw new MathUserException(e);
-                            }
-                        }
-                    };
-                    final BrentSolver solver = new BrentSolver(convergence);
                     final double root = (ta <= tb) ?
-                        solver.solve(maxIterationCount, f, ta, tb) :
-                        solver.solve(maxIterationCount, f, tb, ta);
-                    if ((FastMath.abs(root - ta) <= convergence) &&
-                         (FastMath.abs(root - previousEventTime) <= convergence)) {
+                                        solver.solve(maxIterationCount, f, ta, tb) :
+                                        solver.solve(maxIterationCount, f, tb, ta);
+
+                    if ((!Double.isNaN(previousEventTime)) &&
+                        (FastMath.abs(root - ta) <= convergence) &&
+                        (FastMath.abs(root - previousEventTime) <= convergence)) {
                         // we have either found nothing or found (again ?) a past event, we simply ignore it
                         ta = tb;
                         ga = gb;
                     } else if (Double.isNaN(previousEventTime) ||
                                (FastMath.abs(previousEventTime - root) > convergence)) {
                         pendingEventTime = root;
-                        if (pendingEvent && (FastMath.abs(t1 - pendingEventTime) <= convergence)) {
-                            // we were already waiting for this event which was
-                            // found during a previous call for a step that was
-                            // rejected, this step must now be accepted since it
-                            // properly ends exactly at the event occurrence
-                            return false;
-                        }
-                        // either we were not waiting for the event or it has
-                        // moved in such a way the step cannot be accepted
                         pendingEvent = true;
                         return true;
+                    } else {
+                        // no sign change: there is no event for now
+                        ta = tb;
+                        ga = gb;
                     }
 
                 } else {
@@ -323,7 +323,7 @@ public class EventState {
         t0 = t;
         g0 = handler.g(t, y);
 
-        if (pendingEvent) {
+        if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) {
             // force the sign to its value "just after the event"
             previousEventTime = t;
             g0Positive        = increasing;
@@ -354,7 +354,7 @@ public class EventState {
     public boolean reset(final double t, final double[] y)
         throws EventException {
 
-        if (! pendingEvent) {
+        if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) {
             return false;
         }
 
@@ -369,4 +369,21 @@ public class EventState {
 
     }
 
+    /** Compare the instance with another event state.
+     * <p>
+     * Event state ordering is based on occurrence time within the last
+     * evaluated step. If no event occurs during the step, a time arbitrarily
+     * set to positive infinity is used.
+     * </p>
+     * @param state other event state to compare the instance to
+     * @return a negative integer, zero, or a positive integer as the event
+     * occurs before, simultaneous, or after the specified event of the
+     * specified state.
+     */
+    public int compareTo(final EventState state) {
+        final double instanceTime = pendingEvent ? pendingEventTime : Double.POSITIVE_INFINITY;
+        final double otherTime = state.pendingEvent ? state.pendingEventTime : Double.POSITIVE_INFINITY;
+        return Double.compare(instanceTime, otherTime);
+    }
+
 }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegrator.java Thu Jan 20 20:57:11 2011
@@ -21,7 +21,6 @@ import org.apache.commons.math.exception
 import org.apache.commons.math.linear.Array2DRowRealMatrix;
 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.IntegratorException;
-import org.apache.commons.math.ode.events.CombinedEventsManager;
 import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;
 import org.apache.commons.math.ode.sampling.StepHandler;
 import org.apache.commons.math.util.FastMath;
@@ -202,19 +201,16 @@ public class AdamsBashforthIntegrator ex
             System.arraycopy(y0, 0, y, 0, n);
         }
         final double[] yDot = new double[n];
-        final double[] yTmp = new double[y0.length];
 
         // set up an interpolator sharing the integrator arrays
         final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator();
         interpolator.reinitialize(y, forward);
-        final NordsieckStepInterpolator interpolatorTmp = new NordsieckStepInterpolator();
-        interpolatorTmp.reinitialize(yTmp, forward);
 
         // set up integration control objects
         for (StepHandler handler : stepHandlers) {
             handler.reset();
         }
-        CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
+        statesInitialized = false;
 
         // compute the initial Nordsieck vector using the configured starter integrator
         start(t0, y, t);
@@ -226,14 +222,12 @@ public class AdamsBashforthIntegrator ex
         double hNew = stepSize;
         interpolator.rescale(hNew);
 
-        boolean lastStep = false;
-        while (!lastStep) {
+        // main integration loop
+        isLastStep = false;
+        do {
 
-            // shift all data
-            interpolator.shift();
-
-            double error = 0;
-            for (boolean loop = true; loop;) {
+            double error = 10;
+            while (error >= 1.0) {
 
                 stepSize = hNew;
 
@@ -249,92 +243,51 @@ public class AdamsBashforthIntegrator ex
                 }
                 error = FastMath.sqrt(error / mainSetDimension);
 
-                if (error <= 1.0) {
-
-                    // predict a first estimate of the state at step end
-                    final double stepEnd = stepStart + stepSize;
-                    interpolator.setInterpolatedTime(stepEnd);
-                    System.arraycopy(interpolator.getInterpolatedState(), 0, yTmp, 0, y0.length);
-
-                    // evaluate the derivative
-                    computeDerivatives(stepEnd, yTmp, yDot);
-
-                    // update Nordsieck vector
-                    final double[] predictedScaled = new double[y0.length];
-                    for (int j = 0; j < y0.length; ++j) {
-                        predictedScaled[j] = stepSize * yDot[j];
-                    }
-                    final Array2DRowRealMatrix nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);
-                    updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
-
-                    // discrete events handling
-                    interpolatorTmp.reinitialize(stepEnd, stepSize, predictedScaled, nordsieckTmp);
-                    interpolatorTmp.storeTime(stepStart);
-                    interpolatorTmp.shift();
-                    interpolatorTmp.storeTime(stepEnd);
-                    if (manager.evaluateStep(interpolatorTmp)) {
-                        final double dt = manager.getEventTime() - stepStart;
-                        if (FastMath.abs(dt) <= FastMath.ulp(stepStart)) {
-                            // we cannot simply truncate the step, reject the current computation
-                            // and let the loop compute another state with the truncated step.
-                            // it is so small (much probably exactly 0 due to limited accuracy)
-                            // that the code above would fail handling it.
-                            // So we set up an artificial 0 size step by copying states
-                            interpolator.storeTime(stepStart);
-                            System.arraycopy(y, 0, yTmp, 0, y0.length);
-                            hNew     = 0;
-                            stepSize = 0;
-                            loop     = false;
-                        } else {
-                            // reject the step to match exactly the next switch time
-                            hNew = dt;
-                            interpolator.rescale(hNew);
-                        }
-                    } else {
-                        // accept the step
-                        scaled    = predictedScaled;
-                        nordsieck = nordsieckTmp;
-                        interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
-                        loop = false;
-                    }
-
-                } else {
+                if (error >= 1.0) {
                     // reject the step and attempt to reduce error by stepsize control
                     final double factor = computeStepGrowShrinkFactor(error);
                     hNew = filterStep(stepSize * factor, forward, false);
                     interpolator.rescale(hNew);
-                }
-
-            }
 
-            // the step has been accepted (may have been truncated)
-            final double nextStep = stepStart + stepSize;
-            System.arraycopy(yTmp, 0, y, 0, n);
-            interpolator.storeTime(nextStep);
-            manager.stepAccepted(nextStep, y);
-            lastStep = manager.stop();
-
-            // provide the step data to the step handler
-            for (StepHandler handler : stepHandlers) {
-                interpolator.setInterpolatedTime(nextStep);
-                handler.handleStep(interpolator, lastStep);
+                }
             }
-            stepStart = nextStep;
 
-            if (!lastStep && manager.reset(stepStart, y)) {
+            // predict a first estimate of the state at step end
+            final double stepEnd = stepStart + stepSize;
+            interpolator.shift();
+            interpolator.setInterpolatedTime(stepEnd);
+            System.arraycopy(interpolator.getInterpolatedState(), 0, y, 0, y0.length);
 
-                // some events handler has triggered changes that
-                // invalidate the derivatives, we need to restart from scratch
-                start(stepStart, y, t);
-                interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
+            // evaluate the derivative
+            computeDerivatives(stepEnd, y, yDot);
 
+            // update Nordsieck vector
+            final double[] predictedScaled = new double[y0.length];
+            for (int j = 0; j < y0.length; ++j) {
+                predictedScaled[j] = stepSize * yDot[j];
             }
-
-            if (! lastStep) {
-                // in some rare cases we may get here with stepSize = 0, for example
-                // when an event occurs at integration start, reducing the first step
-                // to zero; we have to reset the step to some safe non zero value
-                stepSize = filterStep(stepSize, forward, true);
+            final Array2DRowRealMatrix nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);
+            updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
+            interpolator.reinitialize(stepEnd, stepSize, predictedScaled, nordsieckTmp);
+
+            // discrete events handling
+            interpolator.storeTime(stepEnd);
+            stepStart = acceptStep(interpolator, stepHandlers, y, yDot, t);
+            scaled    = predictedScaled;
+            nordsieck = nordsieckTmp;
+            interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
+
+            if (!isLastStep) {
+
+                // prepare next step
+                interpolator.storeTime(stepStart);
+
+                if (resetOccurred) {
+                    // some events handler has triggered changes that
+                    // invalidate the derivatives, we need to restart from scratch
+                    start(stepStart, y, t);
+                    interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
+                }
 
                 // stepsize control for next step
                 final double  factor     = computeStepGrowShrinkFactor(error);
@@ -342,14 +295,21 @@ public class AdamsBashforthIntegrator ex
                 final double  nextT      = stepStart + scaledH;
                 final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
                 hNew = filterStep(scaledH, forward, nextIsLast);
+
+                final double  filteredNextT      = stepStart + hNew;
+                final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t);
+                if (filteredNextIsLast) {
+                    hNew = t - stepStart;
+                }
+
                 interpolator.rescale(hNew);
+
             }
 
-        }
+        } while (!isLastStep);
 
-        final double stopTime  = stepStart;
-        stepStart = Double.NaN;
-        stepSize  = Double.NaN;
+        final double stopTime = stepStart;
+        resetInternalState();
         return stopTime;
 
     }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/AdamsMoultonIntegrator.java Thu Jan 20 20:57:11 2011
@@ -24,7 +24,6 @@ import org.apache.commons.math.linear.Ar
 import org.apache.commons.math.linear.RealMatrixPreservingVisitor;
 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.IntegratorException;
-import org.apache.commons.math.ode.events.CombinedEventsManager;
 import org.apache.commons.math.ode.sampling.NordsieckStepInterpolator;
 import org.apache.commons.math.ode.sampling.StepHandler;
 import org.apache.commons.math.util.FastMath;
@@ -220,19 +219,18 @@ public class AdamsMoultonIntegrator exte
         }
         final double[] yDot = new double[y0.length];
         final double[] yTmp = new double[y0.length];
+        final double[] predictedScaled = new double[y0.length];
+        Array2DRowRealMatrix nordsieckTmp = null;
 
         // set up two interpolators sharing the integrator arrays
         final NordsieckStepInterpolator interpolator = new NordsieckStepInterpolator();
         interpolator.reinitialize(y, forward);
-        final NordsieckStepInterpolator interpolatorTmp = new NordsieckStepInterpolator();
-        interpolatorTmp.reinitialize(yTmp, forward);
 
         // set up integration control objects
         for (StepHandler handler : stepHandlers) {
             handler.reset();
         }
-        CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
-
+        statesInitialized = false;
 
         // compute the initial Nordsieck vector using the configured starter integrator
         start(t0, y, t);
@@ -242,14 +240,11 @@ public class AdamsMoultonIntegrator exte
         double hNew = stepSize;
         interpolator.rescale(hNew);
 
-        boolean lastStep = false;
-        while (!lastStep) {
+        isLastStep = false;
+        do {
 
-            // shift all data
-            interpolator.shift();
-
-            double error = 0;
-            for (boolean loop = true; loop;) {
+            double error = 10;
+            while (error >= 1.0) {
 
                 stepSize = hNew;
 
@@ -262,96 +257,56 @@ public class AdamsMoultonIntegrator exte
                 computeDerivatives(stepEnd, yTmp, yDot);
 
                 // update Nordsieck vector
-                final double[] predictedScaled = new double[y0.length];
                 for (int j = 0; j < y0.length; ++j) {
                     predictedScaled[j] = stepSize * yDot[j];
                 }
-                final Array2DRowRealMatrix nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);
+                nordsieckTmp = updateHighOrderDerivativesPhase1(nordsieck);
                 updateHighOrderDerivativesPhase2(scaled, predictedScaled, nordsieckTmp);
 
                 // apply correction (C in the PECE sequence)
                 error = nordsieckTmp.walkInOptimizedOrder(new Corrector(y, predictedScaled, yTmp));
 
-                if (error <= 1.0) {
-
-                    // evaluate a final estimate of the derivative (second E in the PECE sequence)
-                    computeDerivatives(stepEnd, yTmp, yDot);
-
-                    // update Nordsieck vector
-                    final double[] correctedScaled = new double[y0.length];
-                    for (int j = 0; j < y0.length; ++j) {
-                        correctedScaled[j] = stepSize * yDot[j];
-                    }
-                    updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp);
-
-                    // discrete events handling
-                    interpolatorTmp.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp);
-                    interpolatorTmp.storeTime(stepStart);
-                    interpolatorTmp.shift();
-                    interpolatorTmp.storeTime(stepEnd);
-                    if (manager.evaluateStep(interpolatorTmp)) {
-                        final double dt = manager.getEventTime() - stepStart;
-                        if (FastMath.abs(dt) <= FastMath.ulp(stepStart)) {
-                            // we cannot simply truncate the step, reject the current computation
-                            // and let the loop compute another state with the truncated step.
-                            // it is so small (much probably exactly 0 due to limited accuracy)
-                            // that the code above would fail handling it.
-                            // So we set up an artificial 0 size step by copying states
-                            interpolator.storeTime(stepStart);
-                            System.arraycopy(y, 0, yTmp, 0, y0.length);
-                            hNew     = 0;
-                            stepSize = 0;
-                            loop     = false;
-                        } else {
-                            // reject the step to match exactly the next switch time
-                            hNew = dt;
-                            interpolator.rescale(hNew);
-                        }
-                    } else {
-                        // accept the step
-                        scaled    = correctedScaled;
-                        nordsieck = nordsieckTmp;
-                        interpolator.reinitialize(stepEnd, stepSize, scaled, nordsieck);
-                        loop = false;
-                    }
-
-                } else {
+                if (error >= 1.0) {
                     // reject the step and attempt to reduce error by stepsize control
                     final double factor = computeStepGrowShrinkFactor(error);
                     hNew = filterStep(stepSize * factor, forward, false);
                     interpolator.rescale(hNew);
                 }
-
             }
 
-            // the step has been accepted (may have been truncated)
-            final double nextStep = stepStart + stepSize;
-            System.arraycopy(yTmp, 0, y, 0, n);
-            interpolator.storeTime(nextStep);
-            manager.stepAccepted(nextStep, y);
-            lastStep = manager.stop();
-
-            // provide the step data to the step handler
-            for (StepHandler handler : stepHandlers) {
-                interpolator.setInterpolatedTime(nextStep);
-                handler.handleStep(interpolator, lastStep);
+            // evaluate a final estimate of the derivative (second E in the PECE sequence)
+            final double stepEnd = stepStart + stepSize;
+            computeDerivatives(stepEnd, yTmp, yDot);
+
+            // update Nordsieck vector
+            final double[] correctedScaled = new double[y0.length];
+            for (int j = 0; j < y0.length; ++j) {
+                correctedScaled[j] = stepSize * yDot[j];
             }
-            stepStart = nextStep;
-
-            if (!lastStep && manager.reset(stepStart, y)) {
+            updateHighOrderDerivativesPhase2(predictedScaled, correctedScaled, nordsieckTmp);
 
-                // some events handler has triggered changes that
-                // invalidate the derivatives, we need to restart from scratch
-                start(stepStart, y, t);
-                interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
-
-            }
+            // discrete events handling
+            System.arraycopy(yTmp, 0, y, 0, n);
+            interpolator.reinitialize(stepEnd, stepSize, correctedScaled, nordsieckTmp);
+            interpolator.storeTime(stepStart);
+            interpolator.shift();
+            interpolator.storeTime(stepEnd);
+            stepStart = acceptStep(interpolator, stepHandlers, y, yDot, t);
+            scaled    = correctedScaled;
+            nordsieck = nordsieckTmp;
+
+            if (!isLastStep) {
+
+                // prepare next step
+                interpolator.storeTime(stepStart);
+
+                if (resetOccurred) {
+                    // some events handler has triggered changes that
+                    // invalidate the derivatives, we need to restart from scratch
+                    start(stepStart, y, t);
+                    interpolator.reinitialize(stepStart, stepSize, scaled, nordsieck);
 
-            if (! lastStep) {
-                // in some rare cases we may get here with stepSize = 0, for example
-                // when an event occurs at integration start, reducing the first step
-                // to zero; we have to reset the step to some safe non zero value
-                stepSize = filterStep(stepSize, forward, true);
+                }
 
                 // stepsize control for next step
                 final double  factor     = computeStepGrowShrinkFactor(error);
@@ -359,10 +314,17 @@ public class AdamsMoultonIntegrator exte
                 final double  nextT      = stepStart + scaledH;
                 final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
                 hNew = filterStep(scaledH, forward, nextIsLast);
+
+                final double  filteredNextT      = stepStart + hNew;
+                final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t);
+                if (filteredNextIsLast) {
+                    hNew = t - stepStart;
+                }
+
                 interpolator.rescale(hNew);
             }
 
-        }
+        } while (!isLastStep);
 
         final double stopTime  = stepStart;
         stepStart = Double.NaN;

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/DormandPrince853StepInterpolator.java Thu Jan 20 20:57:11 2011
@@ -21,7 +21,6 @@ import java.io.IOException;
 import java.io.ObjectInput;
 import java.io.ObjectOutput;
 
-import org.apache.commons.math.MathRuntimeException;
 import org.apache.commons.math.exception.MathUserException;
 import org.apache.commons.math.ode.AbstractIntegrator;
 import org.apache.commons.math.ode.sampling.StepInterpolator;
@@ -396,6 +395,7 @@ class DormandPrince853StepInterpolator
 
     double s;
     final double[] yTmp = new double[currentState.length];
+    final double pT = getGlobalPreviousTime();
 
     // k14
     for (int j = 0; j < currentState.length; ++j) {
@@ -404,7 +404,7 @@ class DormandPrince853StepInterpolator
           K14_11 * yDotK[10][j] + K14_12 * yDotK[11][j] + K14_13 * yDotK[12][j];
       yTmp[j] = currentState[j] + h * s;
     }
-    integrator.computeDerivatives(previousTime + C14 * h, yTmp, yDotKLast[0]);
+    integrator.computeDerivatives(pT + C14 * h, yTmp, yDotKLast[0]);
 
     // k15
     for (int j = 0; j < currentState.length; ++j) {
@@ -414,7 +414,7 @@ class DormandPrince853StepInterpolator
          K15_14 * yDotKLast[0][j];
      yTmp[j] = currentState[j] + h * s;
     }
-    integrator.computeDerivatives(previousTime + C15 * h, yTmp, yDotKLast[1]);
+    integrator.computeDerivatives(pT + C15 * h, yTmp, yDotKLast[1]);
 
     // k16
     for (int j = 0; j < currentState.length; ++j) {
@@ -424,7 +424,7 @@ class DormandPrince853StepInterpolator
           K16_14 * yDotKLast[0][j] +  K16_15 * yDotKLast[1][j];
       yTmp[j] = currentState[j] + h * s;
     }
-    integrator.computeDerivatives(previousTime + C16 * h, yTmp, yDotKLast[2]);
+    integrator.computeDerivatives(pT + C16 * h, yTmp, yDotKLast[2]);
 
   }
 
@@ -437,7 +437,9 @@ class DormandPrince853StepInterpolator
       // save the local attributes
       finalizeStep();
     } catch (MathUserException e) {
-      throw MathRuntimeException.createIOException(e);
+        IOException ioe = new IOException(e.getLocalizedMessage());
+        ioe.initCause(e);
+        throw ioe;
     }
     final int dimension = (currentState == null) ? -1 : currentState.length;
     out.writeInt(dimension);

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/EmbeddedRungeKuttaIntegrator.java Thu Jan 20 20:57:11 2011
@@ -20,7 +20,6 @@ package org.apache.commons.math.ode.nons
 import org.apache.commons.math.exception.MathUserException;
 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.IntegratorException;
-import org.apache.commons.math.ode.events.CombinedEventsManager;
 import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
 import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
 import org.apache.commons.math.ode.sampling.StepHandler;
@@ -206,11 +205,12 @@ public abstract class EmbeddedRungeKutta
       System.arraycopy(y0, 0, y, 0, y0.length);
     }
     final double[][] yDotK = new double[stages][y0.length];
-    final double[] yTmp = new double[y0.length];
+    final double[] yTmp    = new double[y0.length];
+    final double[] yDotTmp = new double[y0.length];
 
     // set up an interpolator sharing the integrator arrays
     AbstractStepInterpolator interpolator;
-    if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
+    if (requiresDenseOutput() || (! eventsStates.isEmpty())) {
       final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
       rki.reinitialize(this, yTmp, yDotK, forward);
       interpolator = rki;
@@ -226,16 +226,17 @@ public abstract class EmbeddedRungeKutta
     for (StepHandler handler : stepHandlers) {
         handler.reset();
     }
-    CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
-    boolean lastStep = false;
+    statesInitialized = false;
 
     // main integration loop
-    while (!lastStep) {
+    isLastStep = false;
+    do {
 
       interpolator.shift();
 
-      double error = 0;
-      for (boolean loop = true; loop;) {
+      // iterate over step size, ensuring local normalized error is smaller than 1
+      double error = 10;
+      while (error >= 1.0) {
 
         if (firstTime || !fsal) {
           // first stage
@@ -248,11 +249,11 @@ public abstract class EmbeddedRungeKutta
               for (int i = 0; i < scale.length; ++i) {
                 scale[i] = scalAbsoluteTolerance + scalRelativeTolerance * FastMath.abs(y[i]);
               }
-            } else {
+          } else {
               for (int i = 0; i < scale.length; ++i) {
                 scale[i] = vecAbsoluteTolerance[i] + vecRelativeTolerance[i] * FastMath.abs(y[i]);
               }
-            }
+          }
           hNew = initializeStep(equations, forward, getOrder(), scale,
                                 stepStart, y, yDotK[0], yTmp, yDotK[1]);
           firstTime = false;
@@ -286,83 +287,49 @@ public abstract class EmbeddedRungeKutta
 
         // estimate the error at the end of the step
         error = estimateError(yDotK, y, yTmp, stepSize);
-        if (error <= 1.0) {
-
-          // discrete events handling
-          interpolator.storeTime(stepStart + stepSize);
-          if (manager.evaluateStep(interpolator)) {
-              final double dt = manager.getEventTime() - stepStart;
-              if (FastMath.abs(dt) <= FastMath.ulp(stepStart)) {
-                  // we cannot simply truncate the step, reject the current computation
-                  // and let the loop compute another state with the truncated step.
-                  // it is so small (much probably exactly 0 due to limited accuracy)
-                  // that the code above would fail handling it.
-                  // So we set up an artificial 0 size step by copying states
-                  interpolator.storeTime(stepStart);
-                  System.arraycopy(y, 0, yTmp, 0, y0.length);
-                  hNew     = 0;
-                  stepSize = 0;
-                  loop     = false;
-              } else {
-                  // reject the step to match exactly the next switch time
-                  hNew = dt;
-              }
-          } else {
-            // accept the step
-            loop = false;
-          }
-
-        } else {
+        if (error >= 1.0) {
           // reject the step and attempt to reduce error by stepsize control
           final double factor =
               FastMath.min(maxGrowth,
-                       FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
+                           FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
           hNew = filterStep(stepSize * factor, forward, false);
         }
 
       }
 
-      // the step has been accepted
-      final double nextStep = stepStart + stepSize;
+      // local error is small enough: accept the step, trigger events and step handlers
+      interpolator.storeTime(stepStart + stepSize);
       System.arraycopy(yTmp, 0, y, 0, y0.length);
-      manager.stepAccepted(nextStep, y);
-      lastStep = manager.stop();
+      System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
+      stepStart = acceptStep(interpolator, stepHandlers, y, yDotTmp, t);
 
-      // provide the step data to the step handler
-      interpolator.storeTime(nextStep);
-      for (StepHandler handler : stepHandlers) {
-          handler.handleStep(interpolator, lastStep);
-      }
-      stepStart = nextStep;
+      if (!isLastStep) {
 
-      if (fsal) {
-        // save the last evaluation for the next step
-        System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
-      }
+          // prepare next step
+          interpolator.storeTime(stepStart);
 
-      if (manager.reset(stepStart, y) && ! lastStep) {
-        // some event handler has triggered changes that
-        // invalidate the derivatives, we need to recompute them
-        computeDerivatives(stepStart, y, yDotK[0]);
-      }
+          if (fsal) {
+              // save the last evaluation for the next step
+              System.arraycopy(yDotTmp, 0, yDotK[0], 0, y0.length);
+          }
+
+          // stepsize control for next step
+          final double factor =
+              FastMath.min(maxGrowth, FastMath.max(minReduction, safety * FastMath.pow(error, exp)));
+          final double  scaledH    = stepSize * factor;
+          final double  nextT      = stepStart + scaledH;
+          final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
+          hNew = filterStep(scaledH, forward, nextIsLast);
+
+          final double  filteredNextT      = stepStart + hNew;
+          final boolean filteredNextIsLast = forward ? (filteredNextT >= t) : (filteredNextT <= t);
+          if (filteredNextIsLast) {
+              hNew = t - stepStart;
+          }
 
-      if (! lastStep) {
-        // in some rare cases we may get here with stepSize = 0, for example
-        // when an event occurs at integration start, reducing the first step
-        // to zero; we have to reset the step to some safe non zero value
-          stepSize = filterStep(stepSize, forward, true);
-
-        // stepsize control for next step
-        final double factor = FastMath.min(maxGrowth,
-                                       FastMath.max(minReduction,
-                                                safety * FastMath.pow(error, exp)));
-        final double  scaledH    = stepSize * factor;
-        final double  nextT      = stepStart + scaledH;
-        final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
-        hNew = filterStep(scaledH, forward, nextIsLast);
       }
 
-    }
+    } while (!isLastStep);
 
     final double stopTime = stepStart;
     resetInternalState();

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerIntegrator.java Thu Jan 20 20:57:11 2011
@@ -170,7 +170,7 @@ public class GraggBulirschStoerIntegrato
                                       final double scalRelativeTolerance) {
     super(METHOD_NAME, minStep, maxStep,
           scalAbsoluteTolerance, scalRelativeTolerance);
-    denseOutput = requiresDenseOutput() || (! eventsHandlersManager.isEmpty());
+    denseOutput = requiresDenseOutput() || (! eventsStates.isEmpty());
     setStabilityCheck(true, -1, -1, -1);
     setStepsizeControl(-1, -1, -1, -1);
     setOrderControl(-1, -1, -1);
@@ -193,7 +193,7 @@ public class GraggBulirschStoerIntegrato
                                       final double[] vecRelativeTolerance) {
     super(METHOD_NAME, minStep, maxStep,
           vecAbsoluteTolerance, vecRelativeTolerance);
-    denseOutput = requiresDenseOutput() || (! eventsHandlersManager.isEmpty());
+    denseOutput = requiresDenseOutput() || (! eventsStates.isEmpty());
     setStabilityCheck(true, -1, -1, -1);
     setStepsizeControl(-1, -1, -1, -1);
     setOrderControl(-1, -1, -1);
@@ -339,7 +339,7 @@ public class GraggBulirschStoerIntegrato
   public void addStepHandler (final StepHandler handler) {
 
     super.addStepHandler(handler);
-    denseOutput = requiresDenseOutput() || (! eventsHandlersManager.isEmpty());
+    denseOutput = requiresDenseOutput() || (! eventsStates.isEmpty());
 
     // reinitialize the arrays
     initializeArrays();
@@ -353,7 +353,7 @@ public class GraggBulirschStoerIntegrato
                               final double convergence,
                               final int maxIterationCount) {
     super.addEventHandler(function, maxCheckInterval, convergence, maxIterationCount);
-    denseOutput = requiresDenseOutput() || (! eventsHandlersManager.isEmpty());
+    denseOutput = requiresDenseOutput() || (! eventsStates.isEmpty());
 
     // reinitialize the arrays
     initializeArrays();
@@ -556,7 +556,7 @@ public class GraggBulirschStoerIntegrato
   @Override
   public double integrate(final FirstOrderDifferentialEquations equations,
                           final double t0, final double[] y0, final double t, final double[] y)
-  throws MathUserException, IntegratorException {
+      throws MathUserException, IntegratorException {
 
     sanityChecks(equations, t0, y0, t, y);
     setEquations(equations);
@@ -594,10 +594,9 @@ public class GraggBulirschStoerIntegrato
       System.arraycopy(y0, 0, y, 0, y0.length);
     }
 
-    double[] yDot1      = null;
+    double[] yDot1      = new double[y0.length];
     double[][] yMidDots = null;
     if (denseOutput) {
-      yDot1    = new double[y0.length];
       yMidDots = new double[1 + 2 * sequence.length][];
       for (int j = 0; j < yMidDots.length; ++j) {
         yMidDots[j] = new double[y0.length];
@@ -614,13 +613,13 @@ public class GraggBulirschStoerIntegrato
     // initial order selection
     final double tol =
         (vecRelativeTolerance == null) ? scalRelativeTolerance : vecRelativeTolerance[0];
-    final double log10R = FastMath.log(FastMath.max(1.0e-10, tol)) / FastMath.log(10.0);
+    final double log10R = FastMath.log10(FastMath.max(1.0e-10, tol));
     int targetIter = FastMath.max(1,
                               FastMath.min(sequence.length - 2,
                                        (int) FastMath.floor(0.5 - 0.6 * log10R)));
     // set up an interpolator sharing the integrator arrays
     AbstractStepInterpolator interpolator = null;
-    if (denseOutput || (! eventsHandlersManager.isEmpty())) {
+    if (denseOutput || (! eventsStates.isEmpty())) {
       interpolator = new GraggBulirschStoerStepInterpolator(y, yDot0,
                                                             y1, yDot1,
                                                             yMidDots, forward);
@@ -635,13 +634,14 @@ public class GraggBulirschStoerIntegrato
     boolean previousRejected = false;
     boolean firstTime        = true;
     boolean newStep          = true;
-    boolean lastStep         = false;
     boolean firstStepAlreadyComputed = false;
     for (StepHandler handler : stepHandlers) {
         handler.reset();
     }
+    statesInitialized = false;
     costPerTimeUnit[0] = 0;
-    while (! lastStep) {
+    isLastStep = false;
+    do {
 
       double error;
       boolean reject = false;
@@ -656,15 +656,9 @@ public class GraggBulirschStoerIntegrato
         }
 
         if (firstTime) {
-
           hNew = initializeStep(equations, forward,
                                 2 * targetIter + 1, scale,
                                 stepStart, y, yDot0, yTmp, yTmpDot);
-
-          if (! forward) {
-            hNew = -hNew;
-          }
-
         }
 
         newStep = false;
@@ -679,7 +673,7 @@ public class GraggBulirschStoerIntegrato
         stepSize = t - stepStart;
       }
       final double nextT = stepStart + stepSize;
-      lastStep = forward ? (nextT >= t) : (nextT <= t);
+      isLastStep = forward ? (nextT >= t) : (nextT <= t);
 
       // iterate over several substep sizes
       int k = -1;
@@ -804,7 +798,7 @@ public class GraggBulirschStoerIntegrato
                 break;
 
               default :
-                if ((firstTime || lastStep) && (error <= 1.0)) {
+                if ((firstTime || isLastStep) && (error <= 1.0)) {
                   loop = false;
                 }
                 break;
@@ -816,6 +810,11 @@ public class GraggBulirschStoerIntegrato
         }
       }
 
+      if (! reject) {
+          // derivatives at end of step
+          computeDerivatives(stepStart + stepSize, y1, yDot1);
+      }
+
       // dense output handling
       double hInt = getMaxStep();
       if (denseOutput && ! reject) {
@@ -825,9 +824,6 @@ public class GraggBulirschStoerIntegrato
           extrapolate(0, j, diagonal, yMidDots[0]);
         }
 
-        // derivative at end of step
-        computeDerivatives(stepStart + stepSize, y1, yDot1);
-
         final int mu = 2 * k - mudif + 3;
 
         for (int l = 0; l < mu; ++l) {
@@ -880,52 +876,21 @@ public class GraggBulirschStoerIntegrato
             }
           }
 
-          // Discrete events handling
-          if (!reject) {
-            interpolator.storeTime(stepStart + stepSize);
-            if (eventsHandlersManager.evaluateStep(interpolator)) {
-                final double dt = eventsHandlersManager.getEventTime() - stepStart;
-                if (FastMath.abs(dt) > FastMath.ulp(stepStart)) {
-                    // reject the step to match exactly the next switch time
-                    hNew = FastMath.abs(dt);
-                    reject = true;
-                }
-            }
-          }
-
-        }
-
-        if (!reject) {
-          // we will reuse the slope for the beginning of next step
-          firstStepAlreadyComputed = true;
-          System.arraycopy(yDot1, 0, yDot0, 0, y0.length);
         }
 
       }
 
       if (! reject) {
 
-        // store end of step state
-        final double nextStep = stepStart + stepSize;
-        System.arraycopy(y1, 0, y, 0, y0.length);
+        // Discrete events handling
+        interpolator.storeTime(stepStart + stepSize);
+        stepStart = acceptStep(interpolator, stepHandlers, y1, yDot1, t);
 
-        eventsHandlersManager.stepAccepted(nextStep, y);
-        if (eventsHandlersManager.stop()) {
-          lastStep = true;
-        }
-
-        // provide the step data to the step handler
-        interpolator.storeTime(nextStep);
-        for (StepHandler handler : stepHandlers) {
-            handler.handleStep(interpolator, lastStep);
-        }
-        stepStart = nextStep;
-
-        if (eventsHandlersManager.reset(stepStart, y) && ! lastStep) {
-          // some switching function has triggered changes that
-          // invalidate the derivatives, we need to recompute them
-          firstStepAlreadyComputed = false;
-        }
+        // prepare next step
+        interpolator.storeTime(stepStart);
+        System.arraycopy(y1, 0, y, 0, y0.length);
+        System.arraycopy(yDot1, 0, yDot0, 0, y0.length);
+        firstStepAlreadyComputed = true;
 
         int optimalIter;
         if (k == 1) {
@@ -987,15 +952,17 @@ public class GraggBulirschStoerIntegrato
       firstTime = false;
 
       if (reject) {
-        lastStep = false;
+        isLastStep = false;
         previousRejected = true;
       } else {
         previousRejected = false;
       }
 
-    }
+    } while (!isLastStep);
 
-    return stepStart;
+    final double stopTime = stepStart;
+    resetInternalState();
+    return stopTime;
 
   }
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/GraggBulirschStoerStepInterpolator.java Thu Jan 20 20:57:11 2011
@@ -21,7 +21,6 @@ import java.io.IOException;
 import java.io.ObjectInput;
 import java.io.ObjectOutput;
 
-import org.apache.commons.math.exception.MathUserException;
 import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
 import org.apache.commons.math.ode.sampling.StepInterpolator;
 import org.apache.commons.math.util.FastMath;
@@ -310,8 +309,7 @@ class GraggBulirschStoerStepInterpolator
   /** {@inheritDoc} */
   @Override
   protected void computeInterpolatedStateAndDerivatives(final double theta,
-                                          final double oneMinusThetaH)
-    throws MathUserException {
+                                          final double oneMinusThetaH) {
 
     final int dimension = currentState.length;
 

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/nonstiff/RungeKuttaIntegrator.java Thu Jan 20 20:57:11 2011
@@ -22,7 +22,6 @@ import org.apache.commons.math.exception
 import org.apache.commons.math.ode.AbstractIntegrator;
 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.IntegratorException;
-import org.apache.commons.math.ode.events.CombinedEventsManager;
 import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
 import org.apache.commons.math.ode.sampling.DummyStepInterpolator;
 import org.apache.commons.math.ode.sampling.StepHandler;
@@ -112,11 +111,12 @@ public abstract class RungeKuttaIntegrat
     for (int i = 0; i < stages; ++i) {
       yDotK [i] = new double[y0.length];
     }
-    final double[] yTmp = new double[y0.length];
+    final double[] yTmp    = new double[y0.length];
+    final double[] yDotTmp = new double[y0.length];
 
     // set up an interpolator sharing the integrator arrays
     AbstractStepInterpolator interpolator;
-    if (requiresDenseOutput() || (! eventsHandlersManager.isEmpty())) {
+    if (requiresDenseOutput() || (! eventsStates.isEmpty())) {
       final RungeKuttaStepInterpolator rki = (RungeKuttaStepInterpolator) prototype.copy();
       rki.reinitialize(this, yTmp, yDotK, forward);
       interpolator = rki;
@@ -131,90 +131,61 @@ public abstract class RungeKuttaIntegrat
     for (StepHandler handler : stepHandlers) {
         handler.reset();
     }
-    CombinedEventsManager manager = addEndTimeChecker(t0, t, eventsHandlersManager);
-    boolean lastStep = false;
+    statesInitialized = false;
 
     // main integration loop
-    while (!lastStep) {
+    isLastStep = false;
+    do {
 
       interpolator.shift();
 
-      for (boolean loop = true; loop;) {
+      // first stage
+      computeDerivatives(stepStart, y, yDotK[0]);
 
-        // first stage
-        computeDerivatives(stepStart, y, yDotK[0]);
-
-        // next stages
-        for (int k = 1; k < stages; ++k) {
+      // next stages
+      for (int k = 1; k < stages; ++k) {
 
           for (int j = 0; j < y0.length; ++j) {
-            double sum = a[k-1][0] * yDotK[0][j];
-            for (int l = 1; l < k; ++l) {
-              sum += a[k-1][l] * yDotK[l][j];
-            }
-            yTmp[j] = y[j] + stepSize * sum;
+              double sum = a[k-1][0] * yDotK[0][j];
+              for (int l = 1; l < k; ++l) {
+                  sum += a[k-1][l] * yDotK[l][j];
+              }
+              yTmp[j] = y[j] + stepSize * sum;
           }
 
           computeDerivatives(stepStart + c[k-1] * stepSize, yTmp, yDotK[k]);
 
-        }
+      }
 
-        // estimate the state at the end of the step
-        for (int j = 0; j < y0.length; ++j) {
+      // estimate the state at the end of the step
+      for (int j = 0; j < y0.length; ++j) {
           double sum    = b[0] * yDotK[0][j];
           for (int l = 1; l < stages; ++l) {
-            sum    += b[l] * yDotK[l][j];
+              sum    += b[l] * yDotK[l][j];
           }
           yTmp[j] = y[j] + stepSize * sum;
-        }
-
-        // discrete events handling
-        interpolator.storeTime(stepStart + stepSize);
-        if (manager.evaluateStep(interpolator)) {
-            final double dt = manager.getEventTime() - stepStart;
-            if (FastMath.abs(dt) <= FastMath.ulp(stepStart)) {
-                // we cannot simply truncate the step, reject the current computation
-                // and let the loop compute another state with the truncated step.
-                // it is so small (much probably exactly 0 due to limited accuracy)
-                // that the code above would fail handling it.
-                // So we set up an artificial 0 size step by copying states
-                interpolator.storeTime(stepStart);
-                System.arraycopy(y, 0, yTmp, 0, y0.length);
-                stepSize = 0;
-                loop     = false;
-            } else {
-                // reject the step to match exactly the next switch time
-                stepSize = dt;
-            }
-        } else {
-          loop = false;
-        }
-
       }
 
-      // the step has been accepted
-      final double nextStep = stepStart + stepSize;
+      // discrete events handling
+      interpolator.storeTime(stepStart + stepSize);
       System.arraycopy(yTmp, 0, y, 0, y0.length);
-      manager.stepAccepted(nextStep, y);
-      lastStep = manager.stop();
+      System.arraycopy(yDotK[stages - 1], 0, yDotTmp, 0, y0.length);
+      stepStart = acceptStep(interpolator, stepHandlers, y, yDotTmp, t);
 
-      // provide the step data to the step handler
-      interpolator.storeTime(nextStep);
-      for (StepHandler handler : stepHandlers) {
-          handler.handleStep(interpolator, lastStep);
-      }
-      stepStart = nextStep;
+      if (!isLastStep) {
 
-      if (manager.reset(stepStart, y) && ! lastStep) {
-        // some events handler has triggered changes that
-        // invalidate the derivatives, we need to recompute them
-        computeDerivatives(stepStart, y, yDotK[0]);
-      }
+          // prepare next step
+          interpolator.storeTime(stepStart);
 
-      // make sure step size is set to default before next step
-      stepSize = forward ? step : -step;
+          // stepsize control for next step
+          final double  nextT      = stepStart + stepSize;
+          final boolean nextIsLast = forward ? (nextT >= t) : (nextT <= t);
+          if (nextIsLast) {
+              stepSize = t - stepStart;
+          }
+      }
 
-    }
+    } while (!isLastStep);
 
     final double stopTime = stepStart;
     stepStart = Double.NaN;

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java Thu Jan 20 20:57:11 2011
@@ -21,7 +21,6 @@ import java.io.IOException;
 import java.io.ObjectInput;
 import java.io.ObjectOutput;
 
-import org.apache.commons.math.MathRuntimeException;
 import org.apache.commons.math.exception.MathUserException;
 
 /** This abstract class represents an interpolator over the last step
@@ -44,11 +43,17 @@ import org.apache.commons.math.exception
 public abstract class AbstractStepInterpolator
   implements StepInterpolator {
 
-  /** previous time */
-  protected double previousTime;
+  /** global previous time */
+  private double globalPreviousTime;
 
-  /** current time */
-  protected double currentTime;
+  /** global current time */
+  private double globalCurrentTime;
+
+  /** soft previous time */
+  private double softPreviousTime;
+
+  /** soft current time */
+  private double softCurrentTime;
 
   /** current time step */
   protected double h;
@@ -87,8 +92,10 @@ public abstract class AbstractStepInterp
    * initializing the copy.
    */
   protected AbstractStepInterpolator() {
-    previousTime            = Double.NaN;
-    currentTime             = Double.NaN;
+    globalPreviousTime      = Double.NaN;
+    globalCurrentTime       = Double.NaN;
+    softPreviousTime        = Double.NaN;
+    softCurrentTime         = Double.NaN;
     h                       = Double.NaN;
     interpolatedTime        = Double.NaN;
     currentState            = null;
@@ -106,10 +113,12 @@ public abstract class AbstractStepInterp
    */
   protected AbstractStepInterpolator(final double[] y, final boolean forward) {
 
-    previousTime      = Double.NaN;
-    currentTime       = Double.NaN;
-    h                 = Double.NaN;
-    interpolatedTime  = Double.NaN;
+    globalPreviousTime = Double.NaN;
+    globalCurrentTime  = Double.NaN;
+    softPreviousTime   = Double.NaN;
+    softCurrentTime    = Double.NaN;
+    h                  = Double.NaN;
+    interpolatedTime   = Double.NaN;
 
     currentState            = y;
     interpolatedState       = new double[y.length];
@@ -140,10 +149,12 @@ public abstract class AbstractStepInterp
    */
   protected AbstractStepInterpolator(final AbstractStepInterpolator interpolator) {
 
-    previousTime      = interpolator.previousTime;
-    currentTime       = interpolator.currentTime;
-    h                 = interpolator.h;
-    interpolatedTime  = interpolator.interpolatedTime;
+    globalPreviousTime = interpolator.globalPreviousTime;
+    globalCurrentTime  = interpolator.globalCurrentTime;
+    softPreviousTime   = interpolator.softPreviousTime;
+    softCurrentTime    = interpolator.softCurrentTime;
+    h                  = interpolator.h;
+    interpolatedTime   = interpolator.interpolatedTime;
 
     if (interpolator.currentState != null) {
       currentState            = interpolator.currentState.clone();
@@ -168,10 +179,12 @@ public abstract class AbstractStepInterp
    */
   protected void reinitialize(final double[] y, final boolean isForward) {
 
-    previousTime      = Double.NaN;
-    currentTime       = Double.NaN;
-    h                 = Double.NaN;
-    interpolatedTime  = Double.NaN;
+    globalPreviousTime = Double.NaN;
+    globalCurrentTime  = Double.NaN;
+    softPreviousTime   = Double.NaN;
+    softCurrentTime    = Double.NaN;
+    h                  = Double.NaN;
+    interpolatedTime   = Double.NaN;
 
     currentState            = y;
     interpolatedState       = new double[y.length];
@@ -208,7 +221,9 @@ public abstract class AbstractStepInterp
    * interpolator for future calls to {@link #storeTime storeTime}
    */
   public void shift() {
-    previousTime = currentTime;
+    globalPreviousTime = globalCurrentTime;
+    softPreviousTime   = globalPreviousTime;
+    softCurrentTime    = globalCurrentTime;
   }
 
   /** Store the current step time.
@@ -216,8 +231,9 @@ public abstract class AbstractStepInterp
    */
   public void storeTime(final double t) {
 
-    currentTime = t;
-    h           = currentTime - previousTime;
+    globalCurrentTime = t;
+    softCurrentTime   = globalCurrentTime;
+    h                 = globalCurrentTime - globalPreviousTime;
     setInterpolatedTime(t);
 
     // the step is not finalized anymore
@@ -225,14 +241,53 @@ public abstract class AbstractStepInterp
 
   }
 
-  /** {@inheritDoc} */
+  /** Restrict step range to a limited part of the global step.
+   * <p>
+   * This method can be used to restrict a step and make it appear
+   * as if the original step was smaller. Calling this method
+   * <em>only</em> changes the value returned by {@link #getPreviousTime()}
+   * and {@link #getCurrentTime()}, it does not change any
+   * </p>
+   * @param softPreviousTime start of the restricted step
+   * @param softCurrentTime end of the restricted step
+   */
+  public void setSoftBounds(final double softPreviousTime, final double softCurrentTime) {
+      this.softPreviousTime = softPreviousTime;
+      this.softCurrentTime  = softCurrentTime;
+  }
+
+  /**
+   * Get the previous global grid point time.
+   * @return previous global grid point time
+   */
+  public double getGlobalPreviousTime() {
+    return globalPreviousTime;
+  }
+
+  /**
+   * Get the current global grid point time.
+   * @return current global grid point time
+   */
+  public double getGlobalCurrentTime() {
+    return globalCurrentTime;
+  }
+
+  /**
+   * Get the previous soft grid point time.
+   * @return previous soft grid point time
+   * @see #setSoftBounds(double, double)
+   */
   public double getPreviousTime() {
-    return previousTime;
+    return softPreviousTime;
   }
 
-  /** {@inheritDoc} */
+  /**
+   * Get the current soft grid point time.
+   * @return current soft grid point time
+   * @see #setSoftBounds(double, double)
+   */
   public double getCurrentTime() {
-    return currentTime;
+    return softCurrentTime;
   }
 
   /** {@inheritDoc} */
@@ -270,7 +325,7 @@ public abstract class AbstractStepInterp
 
       // lazy evaluation of the state
       if (dirtyState) {
-          final double oneMinusThetaH = currentTime - interpolatedTime;
+          final double oneMinusThetaH = globalCurrentTime - interpolatedTime;
           final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
           computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
           dirtyState = false;
@@ -285,7 +340,7 @@ public abstract class AbstractStepInterp
 
       // lazy evaluation of the state
       if (dirtyState) {
-          final double oneMinusThetaH = currentTime - interpolatedTime;
+          final double oneMinusThetaH = globalCurrentTime - interpolatedTime;
           final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
           computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
           dirtyState = false;
@@ -376,8 +431,10 @@ public abstract class AbstractStepInterp
     } else {
         out.writeInt(currentState.length);
     }
-    out.writeDouble(previousTime);
-    out.writeDouble(currentTime);
+    out.writeDouble(globalPreviousTime);
+    out.writeDouble(globalCurrentTime);
+    out.writeDouble(softPreviousTime);
+    out.writeDouble(softCurrentTime);
     out.writeDouble(h);
     out.writeBoolean(forward);
 
@@ -396,7 +453,9 @@ public abstract class AbstractStepInterp
     try {
       finalizeStep();
     } catch (MathUserException e) {
-      throw MathRuntimeException.createIOException(e);
+        IOException ioe = new IOException(e.getLocalizedMessage());
+        ioe.initCause(e);
+        throw ioe;
     }
 
   }
@@ -414,11 +473,13 @@ public abstract class AbstractStepInterp
     throws IOException {
 
     final int dimension = in.readInt();
-    previousTime  = in.readDouble();
-    currentTime   = in.readDouble();
-    h             = in.readDouble();
-    forward       = in.readBoolean();
-    dirtyState    = true;
+    globalPreviousTime  = in.readDouble();
+    globalCurrentTime   = in.readDouble();
+    softPreviousTime    = in.readDouble();
+    softCurrentTime     = in.readDouble();
+    h                   = in.readDouble();
+    forward             = in.readBoolean();
+    dirtyState          = true;
 
     if (dimension < 0) {
         currentState = null;

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Thu Jan 20 20:57:11 2011
@@ -169,6 +169,10 @@ The <action> type attribute can be add,u
       </action>
     </release>
     <release version="2.2" date="TBD" description="TBD">
+      <action dev="luc" type="fix" issue="MATH-484">
+          separate discrete event detection from adaptive step size handling in ODE solvers,
+          thus improving robustness, maintainability and speed
+      </action>
       <action dev="luc" type="fix" issue="MATH-467">
         Fixed an awkward statement that triggered a false positive warning.
       </action>

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblem4.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblem4.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblem4.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblem4.java Thu Jan 20 20:57:11 2011
@@ -80,6 +80,20 @@ public TestProblem4 copy() {
     return new EventHandler[] { new Bounce(), new Stop() };
   }
 
+  /**
+   * Get the theoretical events times.
+   * @return theoretical events times
+   */
+  public double[] getTheoreticalEventsTimes() {
+      return new double[] {
+          1 * FastMath.PI - a,
+          2 * FastMath.PI - a,
+          3 * FastMath.PI - a,
+          4 * FastMath.PI - a,
+          12.0
+      };
+  }
+
   @Override
   public void doComputeDerivatives(double t, double[] y, double[] yDot) {
     yDot[0] =  y[1];

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblemAbstract.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblemAbstract.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblemAbstract.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblemAbstract.java Thu Jan 20 20:57:11 2011
@@ -159,6 +159,14 @@ public abstract class TestProblemAbstrac
   }
 
   /**
+   * Get the theoretical events times.
+   * @return theoretical events times
+   */
+  public double[] getTheoreticalEventsTimes() {
+      return new double[0];
+  }
+
+  /**
    * Get the number of calls.
    * @return nuber of calls
    */

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblemHandler.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblemHandler.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblemHandler.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/TestProblemHandler.java Thu Jan 20 20:57:11 2011
@@ -80,7 +80,13 @@ public class TestProblemHandler
         // multistep integrators do not handle the first steps themselves
         // so we have to make sure the integrator we look at has really started its work
         if (!Double.isNaN(expectedStepStart)) {
-            maxTimeError = FastMath.max(maxTimeError, FastMath.abs(start - expectedStepStart));
+            // the step should either start at the end of the integrator step
+            // or at an event if the step is split into several substeps
+            double stepError = FastMath.max(maxTimeError, FastMath.abs(start - expectedStepStart));
+            for (double eventTime : problem.getTheoreticalEventsTimes()) {
+                stepError = FastMath.min(stepError, FastMath.abs(start - eventTime));
+            }
+            maxTimeError = FastMath.max(maxTimeError, stepError);
         }
         expectedStepStart = start + integrator.getCurrentSignedStepsize();
     }

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java?rev=1061508&r1=1061507&r2=1061508&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/ode/nonstiff/AdamsBashforthIntegratorTest.java Thu Jan 20 20:57:11 2011
@@ -146,9 +146,9 @@ public class AdamsBashforthIntegratorTes
             integ.integrate(pb, pb.getInitialTime(), pb.getInitialState(),
                             pb.getFinalTime(), new double[pb.getDimension()]);
             if (nSteps < 4) {
-                assertTrue(integ.getEvaluations() > 160);
+                assertTrue(integ.getEvaluations() > 150);
             } else {
-                assertTrue(integ.getEvaluations() < 80);
+                assertTrue(integ.getEvaluations() < 70);
             }
         }
 



Mime
View raw message