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