001    /*
002     * Licensed to the Apache Software Foundation (ASF) under one or more
003     * contributor license agreements.  See the NOTICE file distributed with
004     * this work for additional information regarding copyright ownership.
005     * The ASF licenses this file to You under the Apache License, Version 2.0
006     * (the "License"); you may not use this file except in compliance with
007     * the License.  You may obtain a copy of the License at
008     *
009     *      http://www.apache.org/licenses/LICENSE-2.0
010     *
011     * Unless required by applicable law or agreed to in writing, software
012     * distributed under the License is distributed on an "AS IS" BASIS,
013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014     * See the License for the specific language governing permissions and
015     * limitations under the License.
016     */
017    
018    package org.apache.commons.math.ode.sampling;
019    
020    import java.io.IOException;
021    import java.io.ObjectInput;
022    import java.io.ObjectOutput;
023    
024    import org.apache.commons.math.MathRuntimeException;
025    import org.apache.commons.math.ode.DerivativeException;
026    
027    /** This abstract class represents an interpolator over the last step
028     * during an ODE integration.
029     *
030     * <p>The various ODE integrators provide objects extending this class
031     * to the step handlers. The handlers can use these objects to
032     * retrieve the state vector at intermediate times between the
033     * previous and the current grid points (dense output).</p>
034     *
035     * @see org.apache.commons.math.ode.FirstOrderIntegrator
036     * @see org.apache.commons.math.ode.SecondOrderIntegrator
037     * @see StepHandler
038     *
039     * @version $Revision: 811685 $ $Date: 2009-09-05 13:36:48 -0400 (Sat, 05 Sep 2009) $
040     * @since 1.2
041     *
042     */
043    
044    public abstract class AbstractStepInterpolator
045      implements StepInterpolator {
046    
047      /** previous time */
048      protected double previousTime;
049    
050      /** current time */
051      protected double currentTime;
052    
053      /** current time step */
054      protected double h;
055    
056      /** current state */
057      protected double[] currentState;
058    
059      /** interpolated time */
060      protected double interpolatedTime;
061    
062      /** interpolated state */
063      protected double[] interpolatedState;
064    
065      /** interpolated derivatives */
066      protected double[] interpolatedDerivatives;
067    
068      /** indicate if the step has been finalized or not. */
069      private boolean finalized;
070    
071      /** integration direction. */
072      private boolean forward;
073    
074      /** indicator for dirty state. */
075      private boolean dirtyState;
076    
077    
078      /** Simple constructor.
079       * This constructor builds an instance that is not usable yet, the
080       * {@link #reinitialize} method should be called before using the
081       * instance in order to initialize the internal arrays. This
082       * constructor is used only in order to delay the initialization in
083       * some cases. As an example, the {@link
084       * org.apache.commons.math.ode.nonstiff.EmbeddedRungeKuttaIntegrator}
085       * class uses the prototyping design pattern to create the step
086       * interpolators by cloning an uninitialized model and latter
087       * initializing the copy.
088       */
089      protected AbstractStepInterpolator() {
090        previousTime            = Double.NaN;
091        currentTime             = Double.NaN;
092        h                       = Double.NaN;
093        interpolatedTime        = Double.NaN;
094        currentState            = null;
095        interpolatedState       = null;
096        interpolatedDerivatives = null;
097        finalized               = false;
098        this.forward            = true;
099        this.dirtyState         = true;
100      }
101    
102      /** Simple constructor.
103       * @param y reference to the integrator array holding the state at
104       * the end of the step
105       * @param forward integration direction indicator
106       */
107      protected AbstractStepInterpolator(final double[] y, final boolean forward) {
108    
109        previousTime      = Double.NaN;
110        currentTime       = Double.NaN;
111        h                 = Double.NaN;
112        interpolatedTime  = Double.NaN;
113    
114        currentState            = y;
115        interpolatedState       = new double[y.length];
116        interpolatedDerivatives = new double[y.length];
117    
118        finalized         = false;
119        this.forward      = forward;
120        this.dirtyState   = true;
121    
122      }
123    
124      /** Copy constructor.
125    
126       * <p>The copied interpolator should have been finalized before the
127       * copy, otherwise the copy will not be able to perform correctly
128       * any derivative computation and will throw a {@link
129       * NullPointerException} later. Since we don't want this constructor
130       * to throw the exceptions finalization may involve and since we
131       * don't want this method to modify the state of the copied
132       * interpolator, finalization is <strong>not</strong> done
133       * automatically, it remains under user control.</p>
134    
135       * <p>The copy is a deep copy: its arrays are separated from the
136       * original arrays of the instance.</p>
137    
138       * @param interpolator interpolator to copy from.
139    
140       */
141      protected AbstractStepInterpolator(final AbstractStepInterpolator interpolator) {
142    
143        previousTime      = interpolator.previousTime;
144        currentTime       = interpolator.currentTime;
145        h                 = interpolator.h;
146        interpolatedTime  = interpolator.interpolatedTime;
147    
148        if (interpolator.currentState != null) {
149          currentState            = interpolator.currentState.clone();
150          interpolatedState       = interpolator.interpolatedState.clone();
151          interpolatedDerivatives = interpolator.interpolatedDerivatives.clone();
152        } else {
153          currentState            = null;
154          interpolatedState       = null;
155          interpolatedDerivatives = null;
156        }
157    
158        finalized  = interpolator.finalized;
159        forward    = interpolator.forward;
160        dirtyState = interpolator.dirtyState;
161    
162      }
163    
164      /** Reinitialize the instance
165       * @param y reference to the integrator array holding the state at
166       * the end of the step
167       * @param isForward integration direction indicator
168       */
169      protected void reinitialize(final double[] y, final boolean isForward) {
170    
171        previousTime      = Double.NaN;
172        currentTime       = Double.NaN;
173        h                 = Double.NaN;
174        interpolatedTime  = Double.NaN;
175    
176        currentState            = y;
177        interpolatedState       = new double[y.length];
178        interpolatedDerivatives = new double[y.length];
179    
180        finalized         = false;
181        this.forward      = isForward;
182        this.dirtyState   = true;
183    
184      }
185    
186      /** {@inheritDoc} */
187       public StepInterpolator copy() throws DerivativeException {
188    
189         // finalize the step before performing copy
190         finalizeStep();
191    
192         // create the new independent instance
193         return doCopy();
194    
195       }
196    
197       /** Really copy the finalized instance.
198        * <p>This method is called by {@link #copy()} after the
199        * step has been finalized. It must perform a deep copy
200        * to have an new instance completely independent for the
201        * original instance.
202        * @return a copy of the finalized instance
203        */
204       protected abstract StepInterpolator doCopy();
205    
206      /** Shift one step forward.
207       * Copy the current time into the previous time, hence preparing the
208       * interpolator for future calls to {@link #storeTime storeTime}
209       */
210      public void shift() {
211        previousTime = currentTime;
212      }
213    
214      /** Store the current step time.
215       * @param t current time
216       */
217      public void storeTime(final double t) {
218    
219        currentTime = t;
220        h           = currentTime - previousTime;
221        setInterpolatedTime(t);
222    
223        // the step is not finalized anymore
224        finalized  = false;
225    
226      }
227    
228      /** {@inheritDoc} */
229      public double getPreviousTime() {
230        return previousTime;
231      }
232    
233      /** {@inheritDoc} */
234      public double getCurrentTime() {
235        return currentTime;
236      }
237    
238      /** {@inheritDoc} */
239      public double getInterpolatedTime() {
240        return interpolatedTime;
241      }
242    
243      /** {@inheritDoc} */
244      public void setInterpolatedTime(final double time) {
245          interpolatedTime = time;
246          dirtyState       = true;
247      }
248    
249      /** {@inheritDoc} */
250      public boolean isForward() {
251        return forward;
252      }
253    
254      /** Compute the state and derivatives at the interpolated time.
255       * This is the main processing method that should be implemented by
256       * the derived classes to perform the interpolation.
257       * @param theta normalized interpolation abscissa within the step
258       * (theta is zero at the previous time step and one at the current time step)
259       * @param oneMinusThetaH time gap between the interpolated time and
260       * the current time
261       * @throws DerivativeException this exception is propagated to the caller if the
262       * underlying user function triggers one
263       */
264      protected abstract void computeInterpolatedStateAndDerivatives(double theta,
265                                                                     double oneMinusThetaH)
266        throws DerivativeException;
267    
268      /** {@inheritDoc} */
269      public double[] getInterpolatedState() throws DerivativeException {
270    
271          // lazy evaluation of the state
272          if (dirtyState) {
273              final double oneMinusThetaH = currentTime - interpolatedTime;
274              final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
275              computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
276              dirtyState = false;
277          }
278    
279          return interpolatedState;
280    
281      }
282    
283      /** {@inheritDoc} */
284      public double[] getInterpolatedDerivatives() throws DerivativeException {
285    
286          // lazy evaluation of the state
287          if (dirtyState) {
288              final double oneMinusThetaH = currentTime - interpolatedTime;
289              final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
290              computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
291              dirtyState = false;
292          }
293    
294          return interpolatedDerivatives;
295    
296      }
297    
298      /**
299       * Finalize the step.
300    
301       * <p>Some embedded Runge-Kutta integrators need fewer functions
302       * evaluations than their counterpart step interpolators. These
303       * interpolators should perform the last evaluations they need by
304       * themselves only if they need them. This method triggers these
305       * extra evaluations. It can be called directly by the user step
306       * handler and it is called automatically if {@link
307       * #setInterpolatedTime} is called.</p>
308    
309       * <p>Once this method has been called, <strong>no</strong> other
310       * evaluation will be performed on this step. If there is a need to
311       * have some side effects between the step handler and the
312       * differential equations (for example update some data in the
313       * equations once the step has been done), it is advised to call
314       * this method explicitly from the step handler before these side
315       * effects are set up. If the step handler induces no side effect,
316       * then this method can safely be ignored, it will be called
317       * transparently as needed.</p>
318    
319       * <p><strong>Warning</strong>: since the step interpolator provided
320       * to the step handler as a parameter of the {@link
321       * StepHandler#handleStep handleStep} is valid only for the duration
322       * of the {@link StepHandler#handleStep handleStep} call, one cannot
323       * simply store a reference and reuse it later. One should first
324       * finalize the instance, then copy this finalized instance into a
325       * new object that can be kept.</p>
326    
327       * <p>This method calls the protected <code>doFinalize</code> method
328       * if it has never been called during this step and set a flag
329       * indicating that it has been called once. It is the <code>
330       * doFinalize</code> method which should perform the evaluations.
331       * This wrapping prevents from calling <code>doFinalize</code> several
332       * times and hence evaluating the differential equations too often.
333       * Therefore, subclasses are not allowed not reimplement it, they
334       * should rather reimplement <code>doFinalize</code>.</p>
335    
336       * @throws DerivativeException this exception is propagated to the
337       * caller if the underlying user function triggers one
338       */
339      public final void finalizeStep()
340        throws DerivativeException {
341        if (! finalized) {
342          doFinalize();
343          finalized = true;
344        }
345      }
346    
347      /**
348       * Really finalize the step.
349       * The default implementation of this method does nothing.
350       * @throws DerivativeException this exception is propagated to the
351       * caller if the underlying user function triggers one
352       */
353      protected void doFinalize()
354        throws DerivativeException {
355      }
356    
357      /** {@inheritDoc} */
358      public abstract void writeExternal(ObjectOutput out)
359        throws IOException;
360    
361      /** {@inheritDoc} */
362      public abstract void readExternal(ObjectInput in)
363        throws IOException, ClassNotFoundException;
364    
365      /** Save the base state of the instance.
366       * This method performs step finalization if it has not been done
367       * before.
368       * @param out stream where to save the state
369       * @exception IOException in case of write error
370       */
371      protected void writeBaseExternal(final ObjectOutput out)
372        throws IOException {
373    
374        if (currentState == null) {
375            out.writeInt(-1);
376        } else {
377            out.writeInt(currentState.length);
378        }
379        out.writeDouble(previousTime);
380        out.writeDouble(currentTime);
381        out.writeDouble(h);
382        out.writeBoolean(forward);
383    
384        if (currentState != null) {
385            for (int i = 0; i < currentState.length; ++i) {
386                out.writeDouble(currentState[i]);
387            }
388        }
389    
390        out.writeDouble(interpolatedTime);
391    
392        // we do not store the interpolated state,
393        // it will be recomputed as needed after reading
394    
395        // finalize the step (and don't bother saving the now true flag)
396        try {
397          finalizeStep();
398        } catch (DerivativeException e) {
399          throw MathRuntimeException.createIOException(e);
400        }
401    
402      }
403    
404      /** Read the base state of the instance.
405       * This method does <strong>neither</strong> set the interpolated
406       * time nor state. It is up to the derived class to reset it
407       * properly calling the {@link #setInterpolatedTime} method later,
408       * once all rest of the object state has been set up properly.
409       * @param in stream where to read the state from
410       * @return interpolated time be set later by the caller
411       * @exception IOException in case of read error
412       */
413      protected double readBaseExternal(final ObjectInput in)
414        throws IOException {
415    
416        final int dimension = in.readInt();
417        previousTime  = in.readDouble();
418        currentTime   = in.readDouble();
419        h             = in.readDouble();
420        forward       = in.readBoolean();
421        dirtyState    = true;
422    
423        if (dimension < 0) {
424            currentState = null;
425        } else {
426            currentState  = new double[dimension];
427            for (int i = 0; i < currentState.length; ++i) {
428                currentState[i] = in.readDouble();
429            }
430        }
431    
432        // we do NOT handle the interpolated time and state here
433        interpolatedTime        = Double.NaN;
434        interpolatedState       = (dimension < 0) ? null : new double[dimension];
435        interpolatedDerivatives = (dimension < 0) ? null : new double[dimension];
436    
437        finalized = true;
438    
439        return in.readDouble();
440    
441      }
442    
443    }