added a lazy evaluation for step interpolators

moved the possible exceptions from setInterpolatedTime (which now does almost nothing) to getInterpolatedState and getInterpolatedDerivatives (which now do the evaluations)
renamed an internal function to make sure people overriding it don't forget to set the derivatives


git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@782431 13f79535-47bb-0310-9956-ffa450edef68

src/java/org/apache/commons/math/ode/ContinuousOutputModel.java

@@ -233,7 +233,6 @@ public class ContinuousOutputModel
    */
   public void setInterpolatedTime(final double time) {
 
-    try {
       // initialize the search with the complete steps table
       int iMin = 0;
       final StepInterpolator sMin = steps.get(iMin);
@@ -318,17 +317,15 @@ public class ContinuousOutputModel
 
       steps.get(index).setInterpolatedTime(time);
 
-    } catch (DerivativeException de) {
-      throw new MathRuntimeException(de, "unexpected exception caught");
-    }
-
   }
 
   /**
    * Get the state vector of the interpolated point.
    * @return state vector at time {@link #getInterpolatedTime}
+   * @throws DerivativeException if this call induces an automatic
+   * step finalization that throws one
    */
-  public double[] getInterpolatedState() {
+  public double[] getInterpolatedState() throws DerivativeException {
     return steps.get(index).getInterpolatedState();
   }
 

src/java/org/apache/commons/math/ode/sampling/AbstractStepInterpolator.java

@@ -74,6 +74,10 @@ public abstract class AbstractStepInterpolator
   /** integration direction. */
   private boolean forward;
 
+  /** indicator for dirty state. */
+  private boolean dirtyState;
+
+
   /** Simple constructor.
    * This constructor builds an instance that is not usable yet, the
    * {@link #reinitialize} method should be called before using the
@@ -94,6 +98,7 @@ public abstract class AbstractStepInterpolator
     interpolatedDerivatives = null;
     finalized               = false;
     this.forward            = true;
+    this.dirtyState         = true;
   }
 
   /** Simple constructor.
@@ -114,6 +119,7 @@ public abstract class AbstractStepInterpolator
 
     finalized         = false;
     this.forward      = forward;
+    this.dirtyState   = true;
 
   }
 
@@ -151,8 +157,9 @@ public abstract class AbstractStepInterpolator
       interpolatedDerivatives = null;
     }
 
-    finalized = interpolator.finalized;
+    finalized  = interpolator.finalized;
-    forward   = interpolator.forward;
+    forward    = interpolator.forward;
+    dirtyState = interpolator.dirtyState;
 
   }
 
@@ -174,6 +181,7 @@ public abstract class AbstractStepInterpolator
 
     finalized         = false;
     this.forward      = forward;
+    this.dirtyState   = true;
 
   }
 
@@ -210,14 +218,12 @@ public abstract class AbstractStepInterpolator
    */
   public void storeTime(final double t) {
 
-    currentTime      = t;
+    currentTime = t;
-    h                = currentTime - previousTime;
+    h           = currentTime - previousTime;
-    interpolatedTime = t;
+    setInterpolatedTime(t);
-    System.arraycopy(currentState, 0, interpolatedState, 0,
-                     currentState.length);
 
     // the step is not finalized anymore
-    finalized = false;
+    finalized  = false;
 
   }
 
@@ -237,12 +243,9 @@ public abstract class AbstractStepInterpolator
   }
     
   /** {@inheritDoc} */
-  public void setInterpolatedTime(final double time)
+  public void setInterpolatedTime(final double time) {
-      throws DerivativeException {
       interpolatedTime = time;
-      final double oneMinusThetaH = currentTime - interpolatedTime;
+      dirtyState       = true;
-      final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
-      computeInterpolatedState(theta, oneMinusThetaH);
   }
 
   /** {@inheritDoc} */
@@ -250,7 +253,7 @@ public abstract class AbstractStepInterpolator
     return forward;
   }
 
-  /** Compute the state at the interpolated time.
+  /** Compute the state and derivatives at the interpolated time.
    * This is the main processing method that should be implemented by
    * the derived classes to perform the interpolation.
    * @param theta normalized interpolation abscissa within the step
@@ -260,18 +263,38 @@ public abstract class AbstractStepInterpolator
    * @throws DerivativeException this exception is propagated to the caller if the
    * underlying user function triggers one
    */
-  protected abstract void computeInterpolatedState(double theta,
+  protected abstract void computeInterpolatedStateAndDerivatives(double theta,
-                                                   double oneMinusThetaH)
+                                                                 double oneMinusThetaH)
     throws DerivativeException;
     
   /** {@inheritDoc} */
-  public double[] getInterpolatedState() {
+  public double[] getInterpolatedState() throws DerivativeException {
-    return interpolatedState;
+
+      // lazy evaluation of the state
+      if (dirtyState) {
+          final double oneMinusThetaH = currentTime - interpolatedTime;
+          final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
+          computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
+          dirtyState = false;
+      }
+
+      return interpolatedState;
+
   }
 
   /** {@inheritDoc} */
-  public double[] getInterpolatedDerivatives() {
+  public double[] getInterpolatedDerivatives() throws DerivativeException {
-    return interpolatedDerivatives;
+
+      // lazy evaluation of the state
+      if (dirtyState) {
+          final double oneMinusThetaH = currentTime - interpolatedTime;
+          final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
+          computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
+          dirtyState = false;
+      }
+
+      return interpolatedDerivatives;
+
   }
 
   /**
@@ -398,6 +421,7 @@ public abstract class AbstractStepInterpolator
     currentTime   = in.readDouble();
     h             = in.readDouble();
     forward       = in.readBoolean();
+    dirtyState    = true;
 
     if (dimension < 0) {
         currentState = null;

src/java/org/apache/commons/math/ode/sampling/DummyStepInterpolator.java

@@ -17,11 +17,10 @@
 
 package org.apache.commons.math.ode.sampling;
 
+import java.io.IOException;
 import java.io.ObjectInput;
 import java.io.ObjectOutput;
-import java.io.IOException;
 
-import org.apache.commons.math.MathRuntimeException;
 import org.apache.commons.math.ode.DerivativeException;
 import org.apache.commons.math.ode.nonstiff.EmbeddedRungeKuttaIntegrator;
 
@@ -93,7 +92,7 @@ public class DummyStepInterpolator
    * underlying user function triggers one
    */
   @Override
-  protected void computeInterpolatedState(final double theta, final double oneMinusThetaH)
+  protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH)
     throws DerivativeException {
       System.arraycopy(currentState, 0, interpolatedState, 0, currentState.length);
   }
@@ -120,12 +119,8 @@ public class DummyStepInterpolator
     // read the base class 
     final double t = readBaseExternal(in);
 
-    try {
+    // we can now set the interpolated time and state
-      // we can now set the interpolated time and state
+    setInterpolatedTime(t);
-      setInterpolatedTime(t);
-    } catch (DerivativeException e) {
-      throw MathRuntimeException.createIOException(e);
-    }
 
   }
 

src/java/org/apache/commons/math/ode/sampling/NordsieckStepInterpolator.java

@@ -22,12 +22,9 @@ import java.io.ObjectInput;
 import java.io.ObjectOutput;
 import java.util.Arrays;
 
-import org.apache.commons.math.MathRuntimeException;
 import org.apache.commons.math.linear.RealMatrix;
 import org.apache.commons.math.linear.RealMatrixImpl;
 import org.apache.commons.math.linear.RealMatrixPreservingVisitor;
-import org.apache.commons.math.ode.DerivativeException;
-import org.apache.commons.math.ode.nonstiff.AdamsIntegrator;
 
 /**
  * This class implements an interpolator for integrators using Nordsieck representation.
@@ -35,7 +32,8 @@ import org.apache.commons.math.ode.nonstiff.AdamsIntegrator;
  * <p>This interpolator computes dense output around the current point.
  * The interpolation equation is based on Taylor series formulas.
  *
- * @see AdamsIntegrator
+ * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
+ * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
  * @version $Revision$ $Date$
  * @since 2.0
  */
@@ -48,6 +46,14 @@ public class NordsieckStepInterpolator extends AbstractStepInterpolator {
     /** Step size used in the first scaled derivative and Nordsieck vector. */
     private double scalingH;
 
+    /** Reference time for all arrays.
+     * <p>Sometimes, the reference time is the same as previousTime,
+     * sometimes it is the same as currentTime, so we use a separate
+     * field to avoid any confusion.
+     * </p>
+     */
+    private double referenceTime;
+
     /** First scaled derivative. */
     private double[] scaled;
 
@@ -101,34 +107,29 @@ public class NordsieckStepInterpolator extends AbstractStepInterpolator {
     /** Reinitialize the instance
      * <p>Beware that all arrays <em>must</em> be references to integrator
      * arrays, in order to ensure proper update without copy.</p>
+     * @param referenceTime time at which all arrays are defined
      * @param scalingH step size used in the scaled and nordsieck arrays
      * @param scaled reference to the integrator array holding the first
      * scaled derivative
      * @param nordsieck reference to the integrator matrix holding the
      * nordsieck vector
      */
-    public void reinitialize(final double scalingH, final double[] scaled,
+    public void reinitialize(final double referenceTime, final double scalingH,
-                             final RealMatrix nordsieck) {
+                             final double[] scaled, final RealMatrix nordsieck) {
-        this.scalingH  = scalingH;
+        this.referenceTime = referenceTime;
-        this.scaled    = scaled;
+        this.scalingH      = scalingH;
-        this.nordsieck = nordsieck;
+        this.scaled        = scaled;
-    }
+        this.nordsieck     = nordsieck;
+
+        // make sure the state and derivatives will depend on the new arrays
+        setInterpolatedTime(getInterpolatedTime());
 
-    /** Store the current step time.
-     * @param t current time
-     */
-    @Override
-    public void storeTime(final double t) {
-      currentTime      = t;
-      h                = currentTime - previousTime;
-      interpolatedTime = t;
-      computeInterpolatedState(1.0, 0.0);
     }
 
     /** {@inheritDoc} */
     @Override
-    protected void computeInterpolatedState(final double theta, final double oneMinusThetaH) {
+    protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH) {
-        final double x = theta * h;
+        final double x = interpolatedTime - referenceTime;
         nordsieck.walkInOptimizedOrder(new StateEstimator(x, x / scalingH));
     }
 
@@ -254,13 +255,9 @@ public class NordsieckStepInterpolator extends AbstractStepInterpolator {
             nordsieck = null;
         }
 
-        try {
+        if (hasScaled && hasNordsieck) {
-            if (hasScaled && hasNordsieck) {
+            // we can now set the interpolated time and state
-                // we can now set the interpolated time and state
+            setInterpolatedTime(t);
-                setInterpolatedTime(t);
-            }
-        } catch (DerivativeException e) {
-            throw MathRuntimeException.createIOException(e);
         }
 
     }

src/java/org/apache/commons/math/ode/sampling/StepInterpolator.java

@@ -82,11 +82,8 @@ public interface StepInterpolator
    * specific state must be preserved, a copy of the instance must be
    * created using {@link #copy()}.</p>
    * @param time time of the interpolated point
-   * @throws DerivativeException if this call induces an automatic
-   * step finalization that throws one
    */
-  public void setInterpolatedTime(double time)
+  public void setInterpolatedTime(double time);
-    throws DerivativeException;
 
   /**
    * Get the state vector of the interpolated point.
@@ -95,8 +92,11 @@ public interface StepInterpolator
    * to be preserved across several calls.</p>
    * @return state vector at time {@link #getInterpolatedTime}
    * @see #getInterpolatedDerivatives()
+   * @throws DerivativeException if this call induces an automatic
+   * step finalization that throws one
    */
-  public double[] getInterpolatedState();
+  public double[] getInterpolatedState()
+      throws DerivativeException;
 
   /**
    * Get the derivatives of the state vector of the interpolated point.
@@ -105,9 +105,12 @@ public interface StepInterpolator
    * to be preserved across several calls.</p>
    * @return derivatives of the state vector at time {@link #getInterpolatedTime}
    * @see #getInterpolatedState()
+   * @throws DerivativeException if this call induces an automatic
+   * step finalization that throws one
    * @since 2.0
    */
-  public double[] getInterpolatedDerivatives();
+  public double[] getInterpolatedDerivatives()
+      throws DerivativeException;
 
   /** Check if the natural integration direction is forward.
    * <p>This method provides the integration direction as specified by