    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    }