commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From l..@apache.org
Subject svn commit: r919844 - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/ode/jacobians/ test/java/org/apache/commons/math/ode/jacobians/
Date Sat, 06 Mar 2010 20:38:55 GMT
Author: luc
Date: Sat Mar  6 20:38:55 2010
New Revision: 919844

URL: http://svn.apache.org/viewvc?rev=919844&view=rev
Log:
added event handling to ODE integrators with jacobians

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java
  (with props)
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

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java?rev=919844&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java
(added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java
Sat Mar  6 20:38:55 2010
@@ -0,0 +1,203 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.ode.jacobians;
+
+import org.apache.commons.math.ode.events.EventException;
+
+/** This interface represents a handler for discrete events triggered
+ * during ODE integration.
+ *
+ * <p>Some events can be triggered at discrete times as an ODE problem
+ * is solved. This occurs for example when the integration process
+ * should be stopped as some state is reached (G-stop facility) when the
+ * precise date is unknown a priori, or when the derivatives have
+ * discontinuities, or simply when the user wants to monitor some
+ * states boundaries crossings.
+ * </p>
+ *
+ * <p>These events are defined as occurring when a <code>g</code>
+ * switching function sign changes.</p>
+ *
+ * <p>Since events are only problem-dependent and are triggered by the
+ * independent <i>time</i> variable and the state vector, they can
+ * occur at virtually any time, unknown in advance. The integrators will
+ * take care to avoid sign changes inside the steps, they will reduce
+ * the step size when such an event is detected in order to put this
+ * event exactly at the end of the current step. This guarantees that
+ * step interpolation (which always has a one step scope) is relevant
+ * even in presence of discontinuities. This is independent from the
+ * stepsize control provided by integrators that monitor the local
+ * error (this event handling feature is available for all integrators,
+ * including fixed step ones).</p>
+ *
+ * @version $Revision$ $Date$
+ * @since 2.1
+ */
+
+public interface EventHandlerWithJacobians  {
+
+    /** Stop indicator.
+     * <p>This value should be used as the return value of the {@link
+     * #eventOccurred eventOccurred} method when the integration should be
+     * stopped after the event ending the current step.</p>
+     */
+    int STOP = 0;
+
+    /** Reset state indicator.
+     * <p>This value should be used as the return value of the {@link
+     * #eventOccurred eventOccurred} method when the integration should
+     * go on after the event ending the current step, with a new state
+     * vector (which will be retrieved thanks to the {@link #resetState
+     * resetState} method).</p>
+     */
+    int RESET_STATE = 1;
+
+    /** Reset derivatives indicator.
+     * <p>This value should be used as the return value of the {@link
+     * #eventOccurred eventOccurred} method when the integration should
+     * go on after the event ending the current step, with a new derivatives
+     * vector (which will be retrieved thanks to the {@link
+     * org.apache.commons.math.ode.FirstOrderDifferentialEquations#computeDerivatives}
+     * method).</p>
+     */
+    int RESET_DERIVATIVES = 2;
+
+    /** Continue indicator.
+     * <p>This value should be used as the return value of the {@link
+     * #eventOccurred eventOccurred} method when the integration should go
+     * on after the event ending the current step.</p>
+     */
+    int CONTINUE = 3;
+
+    /** Compute the value of the switching function.
+
+     * <p>The discrete events are generated when the sign of this
+     * switching function changes. The integrator will take care to change
+     * the stepsize in such a way these events occur exactly at step boundaries.
+     * The switching function must be continuous in its roots neighborhood
+     * (but not necessarily smooth), as the integrator will need to find its
+     * roots to locate precisely the events.</p>
+
+     * @param t current value of the independent <i>time</i> variable
+     * @param y array containing the current value of the state vector
+     * @param dydy0 array containing the current value of the jacobian of
+     * the state vector with respect to initial state
+     * @param dydp array containing the current value of the jacobian of
+     * the state vector with respect to parameters
+     * @return value of the g switching function
+     * @exception EventException if the switching function cannot be evaluated
+     */
+    double g(double t, double[] y, double[][] dydy0, double[][] dydp)
+        throws EventException;
+
+    /** Handle an event and choose what to do next.
+
+     * <p>This method is called when the integrator has accepted a step
+     * ending exactly on a sign change of the function, just <em>before</em>
+     * the step handler itself is called (see below for scheduling). It
+     * allows the user to update his internal data to acknowledge the fact
+     * the event has been handled (for example setting a flag in the {@link
+     * org.apache.commons.math.ode.jacobians.ParameterizedODEWithJacobians
+     * differential equations} to switch the derivatives computation in
+     * case of discontinuity), or to direct the integrator to either stop
+     * or continue integration, possibly with a reset state or derivatives.</p>
+
+     * <ul>
+     *   <li>if {@link #STOP} is returned, the step handler will be called
+     *   with the <code>isLast</code> flag of the {@link
+     *   org.apache.commons.math.ode.jacobians.StepHandlerWithJacobians#handleStep(
+     *   StepInterpolatorWithJacobians, boolean) handleStep} method set to true and
+     *   the integration will be stopped,</li>
+     *   <li>if {@link #RESET_STATE} is returned, the {@link #resetState
+     *   resetState} method will be called once the step handler has
+     *   finished its task, and the integrator will also recompute the
+     *   derivatives,</li>
+     *   <li>if {@link #RESET_DERIVATIVES} is returned, the integrator
+     *   will recompute the derivatives,
+     *   <li>if {@link #CONTINUE} is returned, no specific action will
+     *   be taken (apart from having called this method) and integration
+     *   will continue.</li>
+     * </ul>
+
+     * <p>The scheduling between this method and the {@link
+     * org.apache.commons.math.ode.jacobians.StepHandlerWithJacobians
+     * StepHandlerWithJacobians} method {@link
+     * org.apache.commons.math.ode.jacobians.StepHandlerWithJacobians#handleStep(
+     * StepInterpolatorWithJacobians, boolean) handleStep(interpolator, isLast)}
+     * is to call this method first and <code>handleStep</code> afterwards. This
+     * scheduling allows the integrator to pass <code>true</code> as the
+     * <code>isLast</code> parameter to the step handler to make it aware the
step
+     * will be the last one if this method returns {@link #STOP}. As the
+     * interpolator may be used to navigate back throughout the last step (as {@link
+     * org.apache.commons.math.ode.sampling.StepNormalizer StepNormalizer}
+     * does for example), user code called by this method and user
+     * code called by step handlers may experience apparently out of order values
+     * of the independent time variable. As an example, if the same user object
+     * implements both this {@link EventHandlerWithJacobians EventHandler} interface and
the
+     * {@link org.apache.commons.math.ode.sampling.FixedStepHandler FixedStepHandler}
+     * interface, a <em>forward</em> integration may call its
+     * <code>eventOccurred</code> method with t = 10 first and call its
+     * <code>handleStep</code> method with t = 9 afterwards. Such out of order
+     * calls are limited to the size of the integration step for {@link
+     * org.apache.commons.math.ode.sampling.StepHandler variable step handlers} and
+     * to the size of the fixed step for {@link
+     * org.apache.commons.math.ode.sampling.FixedStepHandler fixed step handlers}.</p>
+
+     * @param t current value of the independent <i>time</i> variable
+     * @param y array containing the current value of the state vector
+     * @param dydy0 array containing the current value of the jacobian of
+     * the state vector with respect to initial state
+     * @param dydp array containing the current value of the jacobian of
+     * the state vector with respect to parameters
+     * @param increasing if true, the value of the switching function increases
+     * when times increases around event (note that increase is measured with respect
+     * to physical time, not with respect to integration which may go backward in time)
+     * @return indication of what the integrator should do next, this
+     * value must be one of {@link #STOP}, {@link #RESET_STATE},
+     * {@link #RESET_DERIVATIVES} or {@link #CONTINUE}
+     * @exception EventException if the event occurrence triggers an error
+     */
+    int eventOccurred(double t, double[] y, double[][] dydy0, double[][] dydp,
+                      boolean increasing) throws EventException;
+
+    /** Reset the state prior to continue the integration.
+
+     * <p>This method is called after the step handler has returned and
+     * before the next step is started, but only when {@link
+     * #eventOccurred} has itself returned the {@link #RESET_STATE}
+     * indicator. It allows the user to reset the state vector for the
+     * next step, without perturbing the step handler of the finishing
+     * step. If the {@link #eventOccurred} never returns the {@link
+     * #RESET_STATE} indicator, this function will never be called, and it is
+     * safe to leave its body empty.</p>
+
+     * @param t current value of the independent <i>time</i> variable
+     * @param y array containing the current value of the state vector
+     * the new state should be put in the same array
+     * @param dydy0 array containing the current value of the jacobian of
+     * the state vector with respect to initial state, the new jacobian
+     * should be put in the same array
+     * @param dydp array containing the current value of the jacobian of
+     * the state vector with respect to parameters, the new jacobian
+     * should be put in the same array
+     * @exception EventException if the state cannot be reseted
+     */
+    void resetState(double t, double[] y, double[][] dydy0, double[][] dydp)
+    throws EventException;
+
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java
------------------------------------------------------------------------------
    svn:eol-style = native

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/ode/jacobians/EventHandlerWithJacobians.java
------------------------------------------------------------------------------
    svn:keywords = Author Date Id Revision

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=919844&r1=919843&r2=919844&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 20:38:55 2010
@@ -30,6 +30,8 @@
 import org.apache.commons.math.ode.FirstOrderDifferentialEquations;
 import org.apache.commons.math.ode.FirstOrderIntegrator;
 import org.apache.commons.math.ode.IntegratorException;
+import org.apache.commons.math.ode.events.EventException;
+import org.apache.commons.math.ode.events.EventHandler;
 import org.apache.commons.math.ode.sampling.StepHandler;
 import org.apache.commons.math.ode.sampling.StepInterpolator;
 
@@ -129,23 +131,50 @@
      * @see #getStepHandlers()
      */
     public void clearStepHandlers() {
-
-        // preserve the handlers we did not add ourselves
-        final Collection<StepHandler> otherHandlers = new ArrayList<StepHandler>();
-        for (final StepHandler handler : integrator.getStepHandlers()) {
-            if (!(handler instanceof StepHandlerWrapper)) {
-                otherHandlers.add(handler);
-            }
-        }
-
-        // clear all handlers
         integrator.clearStepHandlers();
+    }
+
+    /** Add an event handler to the integrator.
+     * @param handler event handler
+     * @param maxCheckInterval maximal time interval between switching
+     * function checks (this interval prevents missing sign changes in
+     * case the integration steps becomes very large)
+     * @param convergence convergence threshold in the event time search
+     * @param maxIterationCount upper limit of the iteration count in
+     * the event time search
+     * @see #getEventHandlers()
+     * @see #clearEventHandlers()
+     */
+    public void addEventHandler(EventHandlerWithJacobians handler,
+                                double maxCheckInterval,
+                                double convergence,
+                                int maxIterationCount) {
+        integrator.addEventHandler(new EventHandlerWrapper(handler),
+                                   maxCheckInterval, convergence, maxIterationCount);
+    }
 
-        // put back the preserved handlers
-        for (final StepHandler handler : otherHandlers) {
-            integrator.addStepHandler(handler);
+    /** Get all the event handlers that have been added to the integrator.
+     * @return an unmodifiable collection of the added events handlers
+     * @see #addEventHandler(EventHandlerWithJacobians, double, double, int)
+     * @see #clearEventHandlers()
+     */
+    public Collection<EventHandlerWithJacobians> getEventHandlers() {
+        final Collection<EventHandlerWithJacobians> handlers =
+            new ArrayList<EventHandlerWithJacobians>();
+        for (final EventHandler handler : integrator.getEventHandlers()) {
+            if (handler instanceof EventHandlerWrapper) {
+                handlers.add(((EventHandlerWrapper) handler).getHandler());
+            }
         }
+        return handlers;
+    }
 
+    /** Remove all the event handlers that have been added to the integrator.
+     * @see #addEventHandler(EventHandlerWithJacobians, double, double, int)
+     * @see #getEventHandlers()
+     */
+    public void clearEventHandlers() {
+        integrator.clearEventHandlers();
     }
 
     /** Integrate the differential equations and the variational equations up to the given
time.
@@ -217,16 +246,37 @@
         final double stopTime = integrator.integrate(new MappingWrapper(), t0, z, t, z);
 
         // dispatch the final compound state into the state and partial derivatives arrays
+        dispatchCompoundState(z, y, dYdY0, dYdP);
+
+        return stopTime;
+
+    }
+
+    /** Dispatch a compound state array into state and jacobians arrays.
+     * @param z compound state
+     * @param y raw state array to fill
+     * @param dydy0 jacobian array to fill
+     * @param dydp jacobian array to fill
+     */
+    private static void dispatchCompoundState(final double[] z, final double[] y,
+                                              final double[][] dydy0, final double[][] dydp)
{
+
+        final int n = y.length;
+        final int k = dydp[0].length;
+
+        // state
         System.arraycopy(z, 0, y, 0, n);
+
+        // jacobian with respect to initial state
         for (int i = 0; i < n; ++i) {
-            System.arraycopy(z, n * (i + 1), dYdY0[i], 0, n);
+            System.arraycopy(z, n * (i + 1), dydy0[i], 0, n);
         }
+
+        // jacobian with respect to parameters
         for (int i = 0; i < n; ++i) {
-            System.arraycopy(z, n * (n + 1) + i * k, dYdP[i], 0, k);
+            System.arraycopy(z, n * (n + 1) + i * k, dydp[i], 0, k);
         }
 
-        return stopTime;
-
     }
 
     /** Get the current value of the step start time t<sub>i</sub>.
@@ -766,4 +816,62 @@
         }
 
     }
+
+    /** Wrapper for event handlers. */
+    private class EventHandlerWrapper implements EventHandler {
+
+        /** Underlying event handler with jacobians. */
+        private final EventHandlerWithJacobians handler;
+
+        /** State array. */
+        private double[] y;
+
+        /** Jacobian with respect to initial state dy/dy0. */
+        private double[][] dydy0;
+
+        /** Jacobian with respect to parameters dy/dp. */
+        private double[][] dydp;
+
+        /** Simple constructor.
+         * @param handler underlying event handler with jacobians
+         */
+        public EventHandlerWrapper(final EventHandlerWithJacobians handler) {
+            this.handler = handler;
+            final int n = ode.getDimension();
+            final int k = ode.getParametersDimension();
+            y        = new double[n];
+            dydy0    = new double[n][n];
+            dydp     = new double[n][k];
+        }
+
+        /** Get the underlying event handler with jacobians.
+         * @return underlying event handler with jacobians
+         */
+        public EventHandlerWithJacobians getHandler() {
+            return handler;
+        }
+
+        /** {@inheritDoc} */
+        public int eventOccurred(double t, double[] z, boolean increasing)
+            throws EventException {
+            dispatchCompoundState(z, y, dydy0, dydp);
+            return handler.eventOccurred(t, y, dydy0, dydp, increasing);
+        }
+
+        /** {@inheritDoc} */
+        public double g(double t, double[] z)
+            throws EventException {
+            dispatchCompoundState(z, y, dydy0, dydp);
+            return handler.g(t, y, dydy0, dydp);
+        }
+
+        /** {@inheritDoc} */
+        public void resetState(double t, double[] z)
+            throws EventException {
+            dispatchCompoundState(z, y, dydy0, dydp);
+            handler.resetState(t, y, dydy0, dydp);
+        }
+
+    }
+
 }

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=919844&r1=919843&r2=919844&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 20:38:55 2010
@@ -242,6 +242,39 @@
         extInt.integrate(0, y, circle.exactDyDp(0), t, y, dydy0, dydp);
     }
 
+    @Test
+    public void testEventHandler() throws IntegratorException, DerivativeException {
+        FirstOrderIntegrator integ =
+            new DormandPrince54Integrator(1.0e-8, 100.0, 1.0e-10, 1.0e-10);
+        double[] y = new double[] { 0.0, 1.0 };
+        final Circle circle = new Circle(y, 1.0, 1.0, 0.1);
+        double[][] dydy0 = new double[2][2];
+        double[][] dydp  = new double[2][3];
+        double t = 18 * Math.PI;
+        final FirstOrderIntegratorWithJacobians extInt =
+            new FirstOrderIntegratorWithJacobians(integ, circle);
+        extInt.addEventHandler(new EventHandlerWithJacobians() {
+
+            public int eventOccurred(double t, double[] y, double[][] dydy0,
+                                     double[][] dydp, boolean increasing) {
+                Assert.assertEquals(0.1, y[1], 1.0e-11);
+                Assert.assertTrue(!increasing);
+                return STOP;
+            }
+
+            public double g(double t, double[] y, double[][] dydy0,
+                            double[][] dydp) {
+                return y[1] - 0.1;
+            }
+
+            public void resetState(double t, double[] y, double[][] dydy0,
+                                   double[][] dydp) {
+            }
+        }, 10.0, 1.0e-10, 1000);
+        double stopTime = extInt.integrate(0, y, circle.exactDyDp(0), t, y, dydy0, dydp);
+        Assert.assertTrue(stopTime < 5.0 * Math.PI);
+    }
+
     private static class Brusselator implements ParameterizedODEWithJacobians {
 
         private double b;



Mime
View raw message