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