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