    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    }