AbstractStepInterpolator.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  17. package org.apache.commons.math3.ode.sampling;

  18. import java.io.IOException;
  19. import java.io.ObjectInput;
  20. import java.io.ObjectOutput;

  21. import org.apache.commons.math3.exception.MaxCountExceededException;
  22. import org.apache.commons.math3.ode.EquationsMapper;

  23. /** This abstract class represents an interpolator over the last step
  24.  * during an ODE integration.
  25.  *
  26.  * <p>The various ODE integrators provide objects extending this class
  27.  * to the step handlers. The handlers can use these objects to
  28.  * retrieve the state vector at intermediate times between the
  29.  * previous and the current grid points (dense output).</p>
  30.  *
  31.  * @see org.apache.commons.math3.ode.FirstOrderIntegrator
  32.  * @see org.apache.commons.math3.ode.SecondOrderIntegrator
  33.  * @see StepHandler
  34.  *
  35.  * @since 1.2
  36.  *
  37.  */

  38. public abstract class AbstractStepInterpolator
  39.   implements StepInterpolator {

  40.   /** current time step */
  41.   protected double h;

  42.   /** current state */
  43.   protected double[] currentState;

  44.   /** interpolated time */
  45.   protected double interpolatedTime;

  46.   /** interpolated state */
  47.   protected double[] interpolatedState;

  48.   /** interpolated derivatives */
  49.   protected double[] interpolatedDerivatives;

  50.   /** interpolated primary state */
  51.   protected double[] interpolatedPrimaryState;

  52.   /** interpolated primary derivatives */
  53.   protected double[] interpolatedPrimaryDerivatives;

  54.   /** interpolated secondary state */
  55.   protected double[][] interpolatedSecondaryState;

  56.   /** interpolated secondary derivatives */
  57.   protected double[][] interpolatedSecondaryDerivatives;

  58.   /** global previous time */
  59.   private double globalPreviousTime;

  60.   /** global current time */
  61.   private double globalCurrentTime;

  62.   /** soft previous time */
  63.   private double softPreviousTime;

  64.   /** soft current time */
  65.   private double softCurrentTime;

  66.   /** indicate if the step has been finalized or not. */
  67.   private boolean finalized;

  68.   /** integration direction. */
  69.   private boolean forward;

  70.   /** indicator for dirty state. */
  71.   private boolean dirtyState;

  72.   /** Equations mapper for the primary equations set. */
  73.   private EquationsMapper primaryMapper;

  74.   /** Equations mappers for the secondary equations sets. */
  75.   private EquationsMapper[] secondaryMappers;

  76.   /** Simple constructor.
  77.    * This constructor builds an instance that is not usable yet, the
  78.    * {@link #reinitialize} method should be called before using the
  79.    * instance in order to initialize the internal arrays. This
  80.    * constructor is used only in order to delay the initialization in
  81.    * some cases. As an example, the {@link
  82.    * org.apache.commons.math3.ode.nonstiff.EmbeddedRungeKuttaIntegrator}
  83.    * class uses the prototyping design pattern to create the step
  84.    * interpolators by cloning an uninitialized model and latter
  85.    * initializing the copy.
  86.    */
  87.   protected AbstractStepInterpolator() {
  88.     globalPreviousTime = Double.NaN;
  89.     globalCurrentTime  = Double.NaN;
  90.     softPreviousTime   = Double.NaN;
  91.     softCurrentTime    = Double.NaN;
  92.     h                  = Double.NaN;
  93.     interpolatedTime   = Double.NaN;
  94.     currentState       = null;
  95.     finalized          = false;
  96.     this.forward       = true;
  97.     this.dirtyState    = true;
  98.     primaryMapper      = null;
  99.     secondaryMappers   = null;
  100.     allocateInterpolatedArrays(-1);
  101.   }

  102.   /** Simple constructor.
  103.    * @param y reference to the integrator array holding the state at
  104.    * the end of the step
  105.    * @param forward integration direction indicator
  106.    * @param primaryMapper equations mapper for the primary equations set
  107.    * @param secondaryMappers equations mappers for the secondary equations sets
  108.    */
  109.   protected AbstractStepInterpolator(final double[] y, final boolean forward,
  110.                                      final EquationsMapper primaryMapper,
  111.                                      final EquationsMapper[] secondaryMappers) {

  112.     globalPreviousTime    = Double.NaN;
  113.     globalCurrentTime     = Double.NaN;
  114.     softPreviousTime      = Double.NaN;
  115.     softCurrentTime       = Double.NaN;
  116.     h                     = Double.NaN;
  117.     interpolatedTime      = Double.NaN;
  118.     currentState          = y;
  119.     finalized             = false;
  120.     this.forward          = forward;
  121.     this.dirtyState       = true;
  122.     this.primaryMapper    = primaryMapper;
  123.     this.secondaryMappers = (secondaryMappers == null) ? null : secondaryMappers.clone();
  124.     allocateInterpolatedArrays(y.length);

  125.   }

  126.   /** Copy constructor.

  127.    * <p>The copied interpolator should have been finalized before the
  128.    * copy, otherwise the copy will not be able to perform correctly
  129.    * any derivative computation and will throw a {@link
  130.    * NullPointerException} later. Since we don't want this constructor
  131.    * to throw the exceptions finalization may involve and since we
  132.    * don't want this method to modify the state of the copied
  133.    * interpolator, finalization is <strong>not</strong> done
  134.    * automatically, it remains under user control.</p>

  135.    * <p>The copy is a deep copy: its arrays are separated from the
  136.    * original arrays of the instance.</p>

  137.    * @param interpolator interpolator to copy from.

  138.    */
  139.   protected AbstractStepInterpolator(final AbstractStepInterpolator interpolator) {

  140.     globalPreviousTime = interpolator.globalPreviousTime;
  141.     globalCurrentTime  = interpolator.globalCurrentTime;
  142.     softPreviousTime   = interpolator.softPreviousTime;
  143.     softCurrentTime    = interpolator.softCurrentTime;
  144.     h                  = interpolator.h;
  145.     interpolatedTime   = interpolator.interpolatedTime;

  146.     if (interpolator.currentState == null) {
  147.         currentState     = null;
  148.         primaryMapper    = null;
  149.         secondaryMappers = null;
  150.         allocateInterpolatedArrays(-1);
  151.     } else {
  152.       currentState                     = interpolator.currentState.clone();
  153.       interpolatedState                = interpolator.interpolatedState.clone();
  154.       interpolatedDerivatives          = interpolator.interpolatedDerivatives.clone();
  155.       interpolatedPrimaryState         = interpolator.interpolatedPrimaryState.clone();
  156.       interpolatedPrimaryDerivatives   = interpolator.interpolatedPrimaryDerivatives.clone();
  157.       interpolatedSecondaryState       = new double[interpolator.interpolatedSecondaryState.length][];
  158.       interpolatedSecondaryDerivatives = new double[interpolator.interpolatedSecondaryDerivatives.length][];
  159.       for (int i = 0; i < interpolatedSecondaryState.length; ++i) {
  160.           interpolatedSecondaryState[i]       = interpolator.interpolatedSecondaryState[i].clone();
  161.           interpolatedSecondaryDerivatives[i] = interpolator.interpolatedSecondaryDerivatives[i].clone();
  162.       }
  163.     }

  164.     finalized        = interpolator.finalized;
  165.     forward          = interpolator.forward;
  166.     dirtyState       = interpolator.dirtyState;
  167.     primaryMapper    = interpolator.primaryMapper;
  168.     secondaryMappers = (interpolator.secondaryMappers == null) ?
  169.                        null : interpolator.secondaryMappers.clone();

  170.   }

  171.   /** Allocate the various interpolated states arrays.
  172.    * @param dimension total dimension (negative if arrays should be set to null)
  173.    */
  174.   private void allocateInterpolatedArrays(final int dimension) {
  175.       if (dimension < 0) {
  176.           interpolatedState                = null;
  177.           interpolatedDerivatives          = null;
  178.           interpolatedPrimaryState         = null;
  179.           interpolatedPrimaryDerivatives   = null;
  180.           interpolatedSecondaryState       = null;
  181.           interpolatedSecondaryDerivatives = null;
  182.       } else {
  183.           interpolatedState                = new double[dimension];
  184.           interpolatedDerivatives          = new double[dimension];
  185.           interpolatedPrimaryState         = new double[primaryMapper.getDimension()];
  186.           interpolatedPrimaryDerivatives   = new double[primaryMapper.getDimension()];
  187.           if (secondaryMappers == null) {
  188.               interpolatedSecondaryState       = null;
  189.               interpolatedSecondaryDerivatives = null;
  190.           } else {
  191.               interpolatedSecondaryState       = new double[secondaryMappers.length][];
  192.               interpolatedSecondaryDerivatives = new double[secondaryMappers.length][];
  193.               for (int i = 0; i < secondaryMappers.length; ++i) {
  194.                   interpolatedSecondaryState[i]       = new double[secondaryMappers[i].getDimension()];
  195.                   interpolatedSecondaryDerivatives[i] = new double[secondaryMappers[i].getDimension()];
  196.               }
  197.           }
  198.       }
  199.   }

  200.   /** Reinitialize the instance
  201.    * @param y reference to the integrator array holding the state at the end of the step
  202.    * @param isForward integration direction indicator
  203.    * @param primary equations mapper for the primary equations set
  204.    * @param secondary equations mappers for the secondary equations sets
  205.    */
  206.   protected void reinitialize(final double[] y, final boolean isForward,
  207.                               final EquationsMapper primary,
  208.                               final EquationsMapper[] secondary) {

  209.     globalPreviousTime    = Double.NaN;
  210.     globalCurrentTime     = Double.NaN;
  211.     softPreviousTime      = Double.NaN;
  212.     softCurrentTime       = Double.NaN;
  213.     h                     = Double.NaN;
  214.     interpolatedTime      = Double.NaN;
  215.     currentState          = y;
  216.     finalized             = false;
  217.     this.forward          = isForward;
  218.     this.dirtyState       = true;
  219.     this.primaryMapper    = primary;
  220.     this.secondaryMappers = secondary.clone();
  221.     allocateInterpolatedArrays(y.length);

  222.   }

  223.   /** {@inheritDoc} */
  224.    public StepInterpolator copy() throws MaxCountExceededException {

  225.      // finalize the step before performing copy
  226.      finalizeStep();

  227.      // create the new independent instance
  228.      return doCopy();

  229.    }

  230.    /** Really copy the finalized instance.
  231.     * <p>This method is called by {@link #copy()} after the
  232.     * step has been finalized. It must perform a deep copy
  233.     * to have an new instance completely independent for the
  234.     * original instance.
  235.     * @return a copy of the finalized instance
  236.     */
  237.    protected abstract StepInterpolator doCopy();

  238.   /** Shift one step forward.
  239.    * Copy the current time into the previous time, hence preparing the
  240.    * interpolator for future calls to {@link #storeTime storeTime}
  241.    */
  242.   public void shift() {
  243.     globalPreviousTime = globalCurrentTime;
  244.     softPreviousTime   = globalPreviousTime;
  245.     softCurrentTime    = globalCurrentTime;
  246.   }

  247.   /** Store the current step time.
  248.    * @param t current time
  249.    */
  250.   public void storeTime(final double t) {

  251.     globalCurrentTime = t;
  252.     softCurrentTime   = globalCurrentTime;
  253.     h                 = globalCurrentTime - globalPreviousTime;
  254.     setInterpolatedTime(t);

  255.     // the step is not finalized anymore
  256.     finalized  = false;

  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 #getPreviousTime()},
  263.    * it does not change any other property
  264.    * </p>
  265.    * @param softPreviousTime start of the restricted step
  266.    * @since 2.2
  267.    */
  268.   public void setSoftPreviousTime(final double softPreviousTime) {
  269.       this.softPreviousTime = softPreviousTime;
  270.   }

  271.   /** Restrict step range to a limited part of the global step.
  272.    * <p>
  273.    * This method can be used to restrict a step and make it appear
  274.    * as if the original step was smaller. Calling this method
  275.    * <em>only</em> changes the value returned by {@link #getCurrentTime()},
  276.    * it does not change any other property
  277.    * </p>
  278.    * @param softCurrentTime end of the restricted step
  279.    * @since 2.2
  280.    */
  281.   public void setSoftCurrentTime(final double softCurrentTime) {
  282.       this.softCurrentTime  = softCurrentTime;
  283.   }

  284.   /**
  285.    * Get the previous global grid point time.
  286.    * @return previous global grid point time
  287.    */
  288.   public double getGlobalPreviousTime() {
  289.     return globalPreviousTime;
  290.   }

  291.   /**
  292.    * Get the current global grid point time.
  293.    * @return current global grid point time
  294.    */
  295.   public double getGlobalCurrentTime() {
  296.     return globalCurrentTime;
  297.   }

  298.   /**
  299.    * Get the previous soft grid point time.
  300.    * @return previous soft grid point time
  301.    * @see #setSoftPreviousTime(double)
  302.    */
  303.   public double getPreviousTime() {
  304.     return softPreviousTime;
  305.   }

  306.   /**
  307.    * Get the current soft grid point time.
  308.    * @return current soft grid point time
  309.    * @see #setSoftCurrentTime(double)
  310.    */
  311.   public double getCurrentTime() {
  312.     return softCurrentTime;
  313.   }

  314.   /** {@inheritDoc} */
  315.   public double getInterpolatedTime() {
  316.     return interpolatedTime;
  317.   }

  318.   /** {@inheritDoc} */
  319.   public void setInterpolatedTime(final double time) {
  320.       interpolatedTime = time;
  321.       dirtyState       = true;
  322.   }

  323.   /** {@inheritDoc} */
  324.   public boolean isForward() {
  325.     return forward;
  326.   }

  327.   /** Compute the state and derivatives at the interpolated time.
  328.    * This is the main processing method that should be implemented by
  329.    * the derived classes to perform the interpolation.
  330.    * @param theta normalized interpolation abscissa within the step
  331.    * (theta is zero at the previous time step and one at the current time step)
  332.    * @param oneMinusThetaH time gap between the interpolated time and
  333.    * the current time
  334.    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  335.    */
  336.   protected abstract void computeInterpolatedStateAndDerivatives(double theta,
  337.                                                                  double oneMinusThetaH)
  338.       throws MaxCountExceededException;

  339.   /** Lazy evaluation of complete interpolated state.
  340.    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  341.    */
  342.   private void evaluateCompleteInterpolatedState()
  343.       throws MaxCountExceededException {
  344.       // lazy evaluation of the state
  345.       if (dirtyState) {
  346.           final double oneMinusThetaH = globalCurrentTime - interpolatedTime;
  347.           final double theta = (h == 0) ? 0 : (h - oneMinusThetaH) / h;
  348.           computeInterpolatedStateAndDerivatives(theta, oneMinusThetaH);
  349.           dirtyState = false;
  350.       }
  351.   }

  352.   /** {@inheritDoc} */
  353.   public double[] getInterpolatedState() throws MaxCountExceededException {
  354.       evaluateCompleteInterpolatedState();
  355.       primaryMapper.extractEquationData(interpolatedState,
  356.                                         interpolatedPrimaryState);
  357.       return interpolatedPrimaryState;
  358.   }

  359.   /** {@inheritDoc} */
  360.   public double[] getInterpolatedDerivatives() throws MaxCountExceededException {
  361.       evaluateCompleteInterpolatedState();
  362.       primaryMapper.extractEquationData(interpolatedDerivatives,
  363.                                         interpolatedPrimaryDerivatives);
  364.       return interpolatedPrimaryDerivatives;
  365.   }

  366.   /** {@inheritDoc} */
  367.   public double[] getInterpolatedSecondaryState(final int index) throws MaxCountExceededException {
  368.       evaluateCompleteInterpolatedState();
  369.       secondaryMappers[index].extractEquationData(interpolatedState,
  370.                                                   interpolatedSecondaryState[index]);
  371.       return interpolatedSecondaryState[index];
  372.   }

  373.   /** {@inheritDoc} */
  374.   public double[] getInterpolatedSecondaryDerivatives(final int index) throws MaxCountExceededException {
  375.       evaluateCompleteInterpolatedState();
  376.       secondaryMappers[index].extractEquationData(interpolatedDerivatives,
  377.                                                   interpolatedSecondaryDerivatives[index]);
  378.       return interpolatedSecondaryDerivatives[index];
  379.   }

  380.   /**
  381.    * Finalize the step.

  382.    * <p>Some embedded Runge-Kutta integrators need fewer functions
  383.    * evaluations than their counterpart step interpolators. These
  384.    * interpolators should perform the last evaluations they need by
  385.    * themselves only if they need them. This method triggers these
  386.    * extra evaluations. It can be called directly by the user step
  387.    * handler and it is called automatically if {@link
  388.    * #setInterpolatedTime} is called.</p>

  389.    * <p>Once this method has been called, <strong>no</strong> other
  390.    * evaluation will be performed on this step. If there is a need to
  391.    * have some side effects between the step handler and the
  392.    * differential equations (for example update some data in the
  393.    * equations once the step has been done), it is advised to call
  394.    * this method explicitly from the step handler before these side
  395.    * effects are set up. If the step handler induces no side effect,
  396.    * then this method can safely be ignored, it will be called
  397.    * transparently as needed.</p>

  398.    * <p><strong>Warning</strong>: since the step interpolator provided
  399.    * to the step handler as a parameter of the {@link
  400.    * StepHandler#handleStep handleStep} is valid only for the duration
  401.    * of the {@link StepHandler#handleStep handleStep} call, one cannot
  402.    * simply store a reference and reuse it later. One should first
  403.    * finalize the instance, then copy this finalized instance into a
  404.    * new object that can be kept.</p>

  405.    * <p>This method calls the protected <code>doFinalize</code> method
  406.    * if it has never been called during this step and set a flag
  407.    * indicating that it has been called once. It is the <code>
  408.    * doFinalize</code> method which should perform the evaluations.
  409.    * This wrapping prevents from calling <code>doFinalize</code> several
  410.    * times and hence evaluating the differential equations too often.
  411.    * Therefore, subclasses are not allowed not reimplement it, they
  412.    * should rather reimplement <code>doFinalize</code>.</p>

  413.    * @exception MaxCountExceededException if the number of functions evaluations is exceeded

  414.    */
  415.   public final void finalizeStep() throws MaxCountExceededException {
  416.     if (! finalized) {
  417.       doFinalize();
  418.       finalized = true;
  419.     }
  420.   }

  421.   /**
  422.    * Really finalize the step.
  423.    * The default implementation of this method does nothing.
  424.    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  425.    */
  426.   protected void doFinalize() throws MaxCountExceededException {
  427.   }

  428.   /** {@inheritDoc} */
  429.   public abstract void writeExternal(ObjectOutput out)
  430.     throws IOException;

  431.   /** {@inheritDoc} */
  432.   public abstract void readExternal(ObjectInput in)
  433.     throws IOException, ClassNotFoundException;

  434.   /** Save the base state of the instance.
  435.    * This method performs step finalization if it has not been done
  436.    * before.
  437.    * @param out stream where to save the state
  438.    * @exception IOException in case of write error
  439.    */
  440.   protected void writeBaseExternal(final ObjectOutput out)
  441.     throws IOException {

  442.     if (currentState == null) {
  443.         out.writeInt(-1);
  444.     } else {
  445.         out.writeInt(currentState.length);
  446.     }
  447.     out.writeDouble(globalPreviousTime);
  448.     out.writeDouble(globalCurrentTime);
  449.     out.writeDouble(softPreviousTime);
  450.     out.writeDouble(softCurrentTime);
  451.     out.writeDouble(h);
  452.     out.writeBoolean(forward);
  453.     out.writeObject(primaryMapper);
  454.     out.write(secondaryMappers.length);
  455.     for (final EquationsMapper  mapper : secondaryMappers) {
  456.         out.writeObject(mapper);
  457.     }

  458.     if (currentState != null) {
  459.         for (int i = 0; i < currentState.length; ++i) {
  460.             out.writeDouble(currentState[i]);
  461.         }
  462.     }

  463.     out.writeDouble(interpolatedTime);

  464.     // we do not store the interpolated state,
  465.     // it will be recomputed as needed after reading

  466.     try {
  467.         // finalize the step (and don't bother saving the now true flag)
  468.         finalizeStep();
  469.     } catch (MaxCountExceededException mcee) {
  470.         final IOException ioe = new IOException(mcee.getLocalizedMessage());
  471.         ioe.initCause(mcee);
  472.         throw ioe;
  473.     }

  474.   }

  475.   /** Read the base state of the instance.
  476.    * This method does <strong>neither</strong> set the interpolated
  477.    * time nor state. It is up to the derived class to reset it
  478.    * properly calling the {@link #setInterpolatedTime} method later,
  479.    * once all rest of the object state has been set up properly.
  480.    * @param in stream where to read the state from
  481.    * @return interpolated time to be set later by the caller
  482.    * @exception IOException in case of read error
  483.    * @exception ClassNotFoundException if an equation mapper class
  484.    * cannot be found
  485.    */
  486.   protected double readBaseExternal(final ObjectInput in)
  487.     throws IOException, ClassNotFoundException {

  488.     final int dimension = in.readInt();
  489.     globalPreviousTime  = in.readDouble();
  490.     globalCurrentTime   = in.readDouble();
  491.     softPreviousTime    = in.readDouble();
  492.     softCurrentTime     = in.readDouble();
  493.     h                   = in.readDouble();
  494.     forward             = in.readBoolean();
  495.     primaryMapper       = (EquationsMapper) in.readObject();
  496.     secondaryMappers    = new EquationsMapper[in.read()];
  497.     for (int i = 0; i < secondaryMappers.length; ++i) {
  498.         secondaryMappers[i] = (EquationsMapper) in.readObject();
  499.     }
  500.     dirtyState          = true;

  501.     if (dimension < 0) {
  502.         currentState = null;
  503.     } else {
  504.         currentState  = new double[dimension];
  505.         for (int i = 0; i < currentState.length; ++i) {
  506.             currentState[i] = in.readDouble();
  507.         }
  508.     }

  509.     // we do NOT handle the interpolated time and state here
  510.     interpolatedTime = Double.NaN;
  511.     allocateInterpolatedArrays(dimension);

  512.     finalized = true;

  513.     return in.readDouble();

  514.   }

  515. }