Mailing-List: contact commons-dev-help@jakarta.apache.org; run by ezmlm
Precedence: bulk
Reply-To: "Jakarta Commons Developers List" <commons-dev@jakarta.apache.org>
Received-SPF: neutral (herse.apache.org: local policy)
Message-ID: <4560D631.9000702@free.fr>
Date: Sun, 19 Nov 2006 23:09:53 +0100
From: Luc Maisonobe <Luc.Maisonobe@free.fr>
User-Agent: Thunderbird 1.5.0.7 (X11/20060918)
MIME-Version: 1.0
To: Jakarta Commons Developers List <commons-dev@jakarta.apache.org>
Subject: Re: [math] patching a mantissa issue
References: <4560C82B.4040900@free.fr>
 <8a81b4af0611191319t7fe3efd1tbd3a3301b3d8b670@mail.gmail.com>
In-Reply-To: <8a81b4af0611191319t7fe3efd1tbd3a3301b3d8b670@mail.gmail.com>
Content-Type: multipart/mixed;
 boundary="------------060003030606040000030105"

--------------060003030606040000030105
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: quoted-printable

Phil Steitz a =E9crit :
> Luc,
>=20
> I am about to commit the code.  Actually was working on that when I
> saw this.  Once I have the code in, submit the patch against the code
> in svn.

Here is the patch. I hope it is not too large to be submitted through=20
this list. If anybody is upset with this, please accept my apologies.

The commit message I suggest is the following one:

fixed a problem when switching functions triggered derivatives=20
discontinuities

Luc

--------------060003030606040000030105
Content-Type: text/x-patch;
 name="mantissa.patch"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="mantissa.patch"

Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java	(revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/geometry/Rotation.java	(working copy)
@@ -99,22 +99,6 @@
   }
 
   /** Build a rotation from the quaternion coordinates.
-   * @param q0 scalar part of the quaternion
-   * @param q1 first coordinate of the vectorial part of the quaternion
-   * @param q2 second coordinate of the vectorial part of the quaternion
-   * @param q3 third coordinate of the vectorial part of the quaternion
-   * @deprecated since Mantissa 6.3, this method as been deprecated as it
-   * does not properly handles non-normalized quaternions, it should be
-   * replaced by {@link #Rotation(double, double, double, double, boolean)}
-   */
-  public Rotation(double q0, double q1, double q2, double q3) {
-    this.q0 = q0;
-    this.q1 = q1;
-    this.q2 = q2;
-    this.q3 = q3;
-  }
-
-  /** Build a rotation from the quaternion coordinates.
    * <p>A rotation can be built from a <em>normalized</em> quaternion,
    * i.e. a quaternion for which q<sub>0</sub><sup>2</sup> +
    * q<sub>1</sub><sup>2</sup> + q<sub>2</sub><sup>2</sup> +
@@ -120,9 +104,6 @@
    * q<sub>1</sub><sup>2</sup> + q<sub>2</sub><sup>2</sup> +
    * q<sub>3</sub><sup>2</sup> = 1. If the quaternion is not normalized,
    * the constructor can normalize it in a preprocessing step.</p>
-   * <p>This method replaces the {@link #Rotation(double, double,
-   * double, double) constructor using only 4 doubles} which was deprecated
-   * as of version 6.3 of Mantissa.</p>
    * @param q0 scalar part of the quaternion
    * @param q1 first coordinate of the vectorial part of the quaternion
    * @param q2 second coordinate of the vectorial part of the quaternion
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java	(revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/MessagesResources_fr.java	(working copy)
@@ -91,7 +91,7 @@
     { "dimensions mismatch: ODE problem has dimension {0},"
     + " state vector has dimension {1}",
       "incompatibilit\u00e9 de dimensions entre le probl\u00e8me ODE ({0}),"
-    + " et le vecteur d'\u00e9tat ({1})" },
+    + " et le vecteur d''\u00e9tat ({1})" },
     { "too small integration interval: length = {0}",
       "intervalle d''int\u00e9gration trop petit : {0}" },
 
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java	(revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegrator.java	(working copy)
@@ -880,7 +880,11 @@
         interpolator.storeTime(currentT);
         handler.handleStep(interpolator, lastStep);
 
-        switchesHandler.reset(currentT, y);
+        if (switchesHandler.reset(currentT, y) && ! lastStep) {
+          // some switching function has triggered changes that
+          // invalidate the derivatives, we need to recompute them
+          firstStepAlreadyComputed = false;
+        }
 
         int optimalIter;
         if (k == 1) {
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java	(revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaFehlbergIntegrator.java	(working copy)
@@ -311,7 +311,11 @@
         System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
       }
 
-      switchesHandler.reset(currentT, y);
+      if (switchesHandler.reset(currentT, y) && ! lastStep) {
+        // some switching function has triggered changes that
+        // invalidate the derivatives, we need to recompute them
+        equations.computeDerivatives(currentT, y, yDotK[0]);
+      }
 
       if (! lastStep) {
         // stepsize control for next step
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java	(revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/RungeKuttaIntegrator.java	(working copy)
@@ -234,7 +234,11 @@
         System.arraycopy(yDotK[stages - 1], 0, yDotK[0], 0, y0.length);
       }
 
-      switchesHandler.reset(currentT, y);
+      if (switchesHandler.reset(currentT, y) && ! lastStep) {
+        // some switching function has triggered changes that
+        // invalidate the derivatives, we need to recompute them
+        equations.computeDerivatives(currentT, y, yDotK[0]);
+      }
 
       if (needUpdate) {
         // a switching function has changed the step
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java	(revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunction.java	(working copy)
@@ -22,18 +22,20 @@
  * <p>A switching function allows to handle discrete events in
  * integration problems. These events occur for example when the
  * integration process should be stopped as some value is reached
- * (G-stop facility), or when the derivatives has
- * discontinuities. These events are traditionally defined as
- * occurring when a <code>g</code> function sign changes.</p>
+ * (G-stop facility), or when the derivatives have
+ * discontinuities, or simply when the user wants to monitor some
+ * states boundaries crossings. These events are traditionally defined
+ * as occurring when a <code>g</code> function sign changes, hence
+ * the name <em>switching functions</em>.</p>
  *
  * <p>Since events are only problem-dependent and are triggered by the
  * independant <i>time</i> variable and the state vector, they can
- * occur at virtually any time. 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
+ * 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 feature is available on all integrators, including
  * fixed step ones).</p>
@@ -52,14 +54,23 @@
    */
   public static final int STOP = 0;
 
-  /** Reset indicator.
+  /** 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 through the {@link #resetState
+   * vector (which will be retrieved thanks to the {@link #resetState
    * resetState} method).</p>
    */
-  public static final int RESET = 1;
+  public static final 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
+   * FirstOrderDifferentialEquations#computeDerivatives} method).</p>
+   */
+  public static final int RESET_DERIVATIVES = 2;
 
   /** Continue indicator.
    * <p>This value should be used as the return value of the {@link
@@ -66,7 +77,7 @@
    * #eventOccurred eventOccurred} method when the integration should go
    * on after the event ending the current step.</p>
    */
-  public static final int CONTINUE = 2;
+  public static final int CONTINUE = 3;
 
   /** Compute the value of the switching function.
 
@@ -73,8 +84,8 @@
    * <p>Discrete events are generated when the sign of this function
    * changes, the integrator will take care to change the stepsize in
    * such a way these events occur exactly at step boundaries. This
-   * function must be continuous, as the integrator will need to find
-   * its roots to locate the events.</p>
+   * function must be continuous (at least in its roots neighborhood),
+   * as the integrator will need to find its roots to locate the events.</p>
 
    * @param t current value of the independant <i>time</i> variable
    * @param y array containing the current value of the state vector
@@ -88,16 +99,27 @@
    * ending exactly on a sign change of the function, just before the
    * step handler itself is called. It allows the user to update his
    * internal data to acknowledge the fact the event has been handled
-   * (for example setting a flag to switch the derivatives computation
-   * in case of discontinuity), and it allows to direct the integrator
-   * to either stop or continue integration, possibly with a reset
-   * state.</p>
+   * (for example setting a flag in the {@link
+   * FirstOrderDifferentialEquations 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>
 
-   * <p>If {@link #STOP} is returned, the step handler will be called
-   * with the <code>isLast</code> flag of the {@link
-   * StepHandler#handleStep handleStep} method set to true. If {@link
-   * #RESET} is returned, the {@link #resetState resetState} method
-   * will be called once the step handler has finished its task.</p>
+   * <ul>
+   *   <li>if {@link #STOP} is returned, the step handler will be called
+   *   with the <code>isLast</code> flag of the {@link
+   *   StepHandler#handleStep 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>
 
    * @param t current value of the independant <i>time</i> variable
    * @param y array containing the current value of the state vector
@@ -102,8 +124,8 @@
    * @param t current value of the independant <i>time</i> variable
    * @param y array containing the current value of the state vector
    * @return indication of what the integrator should do next, this
-   * value must be one of {@link #STOP}, {@link #RESET} or {@link
-   * #CONTINUE}
+   * value must be one of {@link #STOP}, {@link #RESET_STATE},
+   * {@link #RESET_DERIVATIVES} or {@link #CONTINUE}
    */
   public int eventOccurred(double t, double[] y);
   
@@ -111,11 +133,11 @@
 
    * <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}
+   * #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} indicator, this function will never be called, and it is
+   * #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 independant <i>time</i> variable
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java	(revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchingFunctionsHandler.java	(working copy)
@@ -166,11 +166,16 @@
    * beginning of the next step
    * @param y array were to put the desired state vector at the beginning
    * of the next step
+   * @return true if the integrator should reset the derivatives too
    */
-  public void reset(double t, double[] y) {
+  public boolean reset(double t, double[] y) {
+    boolean resetDerivatives = false;
     for (Iterator iter = functions.iterator(); iter.hasNext();) {
-      ((SwitchState) iter.next()).reset(t, y);
+      if (((SwitchState) iter.next()).reset(t, y)) {
+        resetDerivatives = true;
+      }
     }
+    return resetDerivatives;
   }
 
   /** Switching functions. */
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java	(revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/src/org/spaceroots/mantissa/ode/SwitchState.java	(working copy)
@@ -239,15 +239,23 @@
    * beginning of the next step
    * @param y array were to put the desired state vector at the beginning
    * of the next step
+   * @return true if the integrator should reset the derivatives too
    */
-  public void reset(double t, double[] y) {
-    if (pendingEvent) {
-      if (nextAction == SwitchingFunction.RESET) {
-        function.resetState(t, y);
-      }
-      pendingEvent      = false;
-      pendingEventTime  = Double.NaN;
+  public boolean reset(double t, double[] y) {
+
+    if (! pendingEvent) {
+      return false;
+    }
+
+    if (nextAction == SwitchingFunction.RESET_STATE) {
+      function.resetState(t, y);
     }
+    pendingEvent      = false;
+    pendingEventTime  = Double.NaN;
+
+    return (nextAction == SwitchingFunction.RESET_STATE)
+        || (nextAction == SwitchingFunction.RESET_DERIVATIVES);
+
   }
 
   /** Get the value of the g function at the specified time.
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java	(revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/DormandPrince853IntegratorTest.java	(working copy)
@@ -305,6 +305,17 @@
 
   }
 
+  public void testUnstableDerivative()
+  throws DerivativeException, IntegratorException {
+    final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
+    FirstOrderIntegrator integ =
+      new DormandPrince853Integrator(0.1, 10, 1.0e-12, 0.0);
+    integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12);
+    double[] y = { Double.NaN };
+    integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
+    assertEquals(8.0, y[0], 1.0e-12);
+  }
+
   public static Test suite() {
     return new TestSuite(DormandPrince853IntegratorTest.class);
   }
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java	(revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GillIntegratorTest.java	(working copy)
@@ -196,6 +196,16 @@
                     pb.getFinalTime(), new double[pb.getDimension()]);
   }
 
+  public void testUnstableDerivative()
+  throws DerivativeException, IntegratorException {
+    final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
+    FirstOrderIntegrator integ = new GillIntegrator(0.3);
+    integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12);
+    double[] y = { Double.NaN };
+    integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
+    assertEquals(8.0, y[0], 1.0e-12);
+  }
+
   public static Test suite() {
     return new TestSuite(GillIntegratorTest.class);
   }
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java	(revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/GraggBulirschStoerIntegratorTest.java	(working copy)
@@ -256,6 +256,17 @@
                     pb.getFinalTime(), new double[pb.getDimension()]);
   }
 
+  public void testUnstableDerivative()
+    throws DerivativeException, IntegratorException {
+    final StepProblem stepProblem = new StepProblem(0.0, 1.0, 2.0);
+    FirstOrderIntegrator integ =
+      new GraggBulirschStoerIntegrator(0.1, 10, 1.0e-12, 0.0);
+    integ.addSwitchingFunction(stepProblem, 1.0, 1.0e-12);
+    double[] y = { Double.NaN };
+    integ.integrate(stepProblem, 0.0, new double[] { 0.0 }, 10.0, y);
+    assertEquals(8.0, y[0], 1.0e-12);
+  }
+
   public static Test suite() {
     return new TestSuite(GraggBulirschStoerIntegratorTest.class);
   }
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/StepProblem.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/StepProblem.java	(revision 0)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/StepProblem.java	(revision 0)
@@ -0,0 +1,42 @@
+package org.spaceroots.mantissa.ode;
+
+
+public class StepProblem
+  implements FirstOrderDifferentialEquations, SwitchingFunction {
+
+  public StepProblem(double rateBefore, double rateAfter,
+                     double switchTime) {
+    this.rateAfter  = rateAfter;
+    this.switchTime = switchTime;
+    setRate(rateBefore);
+  }
+
+  public void computeDerivatives(double t, double[] y, double[] yDot) {
+    yDot[0] = rate;
+  }
+
+  public int getDimension() {
+    return 1;
+  }
+
+  public void setRate(double rate) {
+    this.rate = rate;
+  }
+
+  public int eventOccurred(double t, double[] y) {
+    setRate(rateAfter);
+    return RESET_DERIVATIVES;
+  }
+
+  public double g(double t, double[] y) {
+    return t - switchTime;
+  }
+
+  public void resetState(double t, double[] y) {
+  }
+
+  private double rate;
+  private double rateAfter;
+  private double switchTime;
+
+}
Index: /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java
===================================================================
--- /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java	(revision 476939)
+++ /home/luc/sources/jakarta/commons-math/src/mantissa/tests-src/org/spaceroots/mantissa/ode/TestProblem4.java	(working copy)
@@ -104,7 +104,7 @@
     public int eventOccurred(double t, double[] y) {
       // this sign change is needed because the state will be reset soon
       sign = -sign;
-      return SwitchingFunction.RESET;
+      return SwitchingFunction.RESET_STATE;
     }
   
     public void resetState(double t, double[] y) {


--------------060003030606040000030105
Content-Type: text/plain; charset=us-ascii

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