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.  */

  18. package org.apache.commons.math4.legacy.ode.sampling;

  20. import java.io.IOException;
  21. import java.io.ObjectInput;
  22. import java.io.ObjectOutput;

  24. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  25. import org.apache.commons.math4.legacy.ode.EquationsMapper;

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

  43. public abstract class AbstractStepInterpolator
  44.   implements StepInterpolator {

  46.   /** current time step. */
  47.   protected double h;

  49.   /** current state. */
  50.   protected double[] currentState;

  52.   /** interpolated time. */
  53.   protected double interpolatedTime;

  55.   /** interpolated state. */
  56.   protected double[] interpolatedState;

  58.   /** interpolated derivatives. */
  59.   protected double[] interpolatedDerivatives;

  61.   /** interpolated primary state. */
  62.   protected double[] interpolatedPrimaryState;

  64.   /** interpolated primary derivatives. */
  65.   protected double[] interpolatedPrimaryDerivatives;

  67.   /** interpolated secondary state. */
  68.   protected double[][] interpolatedSecondaryState;

  70.   /** interpolated secondary derivatives. */
  71.   protected double[][] interpolatedSecondaryDerivatives;

  73.   /** global previous time. */
  74.   private double globalPreviousTime;

  76.   /** global current time. */
  77.   private double globalCurrentTime;

  79.   /** soft previous time. */
  80.   private double softPreviousTime;

  82.   /** soft current time. */
  83.   private double softCurrentTime;

  85.   /** indicate if the step has been finalized or not. */
  86.   private boolean finalized;

  88.   /** integration direction. */
  89.   private boolean forward;

  91.   /** indicator for dirty state. */
  92.   private boolean dirtyState;

  94.   /** Equations mapper for the primary equations set. */
  95.   private EquationsMapper primaryMapper;

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

  100.   /** Simple constructor.
  101.    * This constructor builds an instance that is not usable yet, the
  102.    * {@link #reinitialize} method should be called before using the
  103.    * instance in order to initialize the internal arrays. This
  104.    * constructor is used only in order to delay the initialization in
  105.    * some cases. As an example, the {@link
  106.    * org.apache.commons.math4.legacy.ode.nonstiff.EmbeddedRungeKuttaIntegrator}
  107.    * class uses the prototyping design pattern to create the step
  108.    * interpolators by cloning an uninitialized model and latter
  109.    * initializing the copy.
  110.    */
  111.   protected AbstractStepInterpolator() {
  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       = null;
  119.     finalized          = false;
  120.     this.forward       = true;
  121.     this.dirtyState    = true;
  122.     primaryMapper      = null;
  123.     secondaryMappers   = null;
  124.     allocateInterpolatedArrays(-1);
  125.   }

  127.   /** Simple constructor.
  128.    * @param y reference to the integrator array holding the state at
  129.    * the end of the step
  130.    * @param forward integration direction indicator
  131.    * @param primaryMapper equations mapper for the primary equations set
  132.    * @param secondaryMappers equations mappers for the secondary equations sets
  133.    */
  134.   protected AbstractStepInterpolator(final double[] y, final boolean forward,
  135.                                      final EquationsMapper primaryMapper,
  136.                                      final EquationsMapper[] secondaryMappers) {

  138.     globalPreviousTime    = Double.NaN;
  139.     globalCurrentTime     = Double.NaN;
  140.     softPreviousTime      = Double.NaN;
  141.     softCurrentTime       = Double.NaN;
  142.     h                     = Double.NaN;
  143.     interpolatedTime      = Double.NaN;
  144.     currentState          = y;
  145.     finalized             = false;
  146.     this.forward          = forward;
  147.     this.dirtyState       = true;
  148.     this.primaryMapper    = primaryMapper;
  149.     this.secondaryMappers = (secondaryMappers == null) ? null : secondaryMappers.clone();
  150.     allocateInterpolatedArrays(y.length);
  151.   }

  153.   /** Copy constructor.

  155.    * <p>The copied interpolator should have been finalized before the
  156.    * copy, otherwise the copy will not be able to perform correctly
  157.    * any derivative computation and will throw a {@link
  158.    * NullPointerException} later. Since we don't want this constructor
  159.    * to throw the exceptions finalization may involve and since we
  160.    * don't want this method to modify the state of the copied
  161.    * interpolator, finalization is <strong>not</strong> done
  162.    * automatically, it remains under user control.</p>

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

  167.    * @param interpolator interpolator to copy from.

  169.    */
  170.   protected AbstractStepInterpolator(final AbstractStepInterpolator interpolator) {

  172.     globalPreviousTime = interpolator.globalPreviousTime;
  173.     globalCurrentTime  = interpolator.globalCurrentTime;
  174.     softPreviousTime   = interpolator.softPreviousTime;
  175.     softCurrentTime    = interpolator.softCurrentTime;
  176.     h                  = interpolator.h;
  177.     interpolatedTime   = interpolator.interpolatedTime;

  179.     if (interpolator.currentState == null) {
  180.         currentState     = null;
  181.         primaryMapper    = null;
  182.         secondaryMappers = null;
  183.         allocateInterpolatedArrays(-1);
  184.     } else {
  185.       currentState                     = interpolator.currentState.clone();
  186.       interpolatedState                = interpolator.interpolatedState.clone();
  187.       interpolatedDerivatives          = interpolator.interpolatedDerivatives.clone();
  188.       interpolatedPrimaryState         = interpolator.interpolatedPrimaryState.clone();
  189.       interpolatedPrimaryDerivatives   = interpolator.interpolatedPrimaryDerivatives.clone();
  190.       interpolatedSecondaryState       = new double[interpolator.interpolatedSecondaryState.length][];
  191.       interpolatedSecondaryDerivatives = new double[interpolator.interpolatedSecondaryDerivatives.length][];
  192.       for (int i = 0; i < interpolatedSecondaryState.length; ++i) {
  193.           interpolatedSecondaryState[i]       = interpolator.interpolatedSecondaryState[i].clone();
  194.           interpolatedSecondaryDerivatives[i] = interpolator.interpolatedSecondaryDerivatives[i].clone();
  195.       }
  196.     }

  198.     finalized        = interpolator.finalized;
  199.     forward          = interpolator.forward;
  200.     dirtyState       = interpolator.dirtyState;
  201.     primaryMapper    = interpolator.primaryMapper;
  202.     secondaryMappers = (interpolator.secondaryMappers == null) ?
  203.                        null : interpolator.secondaryMappers.clone();
  204.   }

  206.   /** Allocate the various interpolated states arrays.
  207.    * @param dimension total dimension (negative if arrays should be set to null)
  208.    */
  209.   private void allocateInterpolatedArrays(final int dimension) {
  210.       if (dimension < 0) {
  211.           interpolatedState                = null;
  212.           interpolatedDerivatives          = null;
  213.           interpolatedPrimaryState         = null;
  214.           interpolatedPrimaryDerivatives   = null;
  215.           interpolatedSecondaryState       = null;
  216.           interpolatedSecondaryDerivatives = null;
  217.       } else {
  218.           interpolatedState                = new double[dimension];
  219.           interpolatedDerivatives          = new double[dimension];
  220.           interpolatedPrimaryState         = new double[primaryMapper.getDimension()];
  221.           interpolatedPrimaryDerivatives   = new double[primaryMapper.getDimension()];
  222.           if (secondaryMappers == null) {
  223.               interpolatedSecondaryState       = null;
  224.               interpolatedSecondaryDerivatives = null;
  225.           } else {
  226.               interpolatedSecondaryState       = new double[secondaryMappers.length][];
  227.               interpolatedSecondaryDerivatives = new double[secondaryMappers.length][];
  228.               for (int i = 0; i < secondaryMappers.length; ++i) {
  229.                   interpolatedSecondaryState[i]       = new double[secondaryMappers[i].getDimension()];
  230.                   interpolatedSecondaryDerivatives[i] = new double[secondaryMappers[i].getDimension()];
  231.               }
  232.           }
  233.       }
  234.   }

  236.   /** Reinitialize the instance.
  237.    * @param y reference to the integrator array holding the state at the end of the step
  238.    * @param isForward integration direction indicator
  239.    * @param primary equations mapper for the primary equations set
  240.    * @param secondary equations mappers for the secondary equations sets
  241.    */
  242.   protected void reinitialize(final double[] y, final boolean isForward,
  243.                               final EquationsMapper primary,
  244.                               final EquationsMapper[] secondary) {

  246.     globalPreviousTime    = Double.NaN;
  247.     globalCurrentTime     = Double.NaN;
  248.     softPreviousTime      = Double.NaN;
  249.     softCurrentTime       = Double.NaN;
  250.     h                     = Double.NaN;
  251.     interpolatedTime      = Double.NaN;
  252.     currentState          = y;
  253.     finalized             = false;
  254.     this.forward          = isForward;
  255.     this.dirtyState       = true;
  256.     this.primaryMapper    = primary;
  257.     this.secondaryMappers = secondary.clone();
  258.     allocateInterpolatedArrays(y.length);
  259.   }

  261.   /** {@inheritDoc} */
  262.   @Override
  263.    public StepInterpolator copy() throws MaxCountExceededException {

  265.      // finalize the step before performing copy
  266.      finalizeStep();

  268.      // create the new independent instance
  269.      return doCopy();
  270.    }

  272.    /** Really copy the finalized instance.
  273.     * <p>This method is called by {@link #copy()} after the
  274.     * step has been finalized. It must perform a deep copy
  275.     * to have an new instance completely independent for the
  276.     * original instance.
  277.     * @return a copy of the finalized instance
  278.     */
  279.    protected abstract StepInterpolator doCopy();

  281.   /** Shift one step forward.
  282.    * Copy the current time into the previous time, hence preparing the
  283.    * interpolator for future calls to {@link #storeTime storeTime}
  284.    */
  285.   public void shift() {
  286.     globalPreviousTime = globalCurrentTime;
  287.     softPreviousTime   = globalPreviousTime;
  288.     softCurrentTime    = globalCurrentTime;
  289.   }

  291.   /** Store the current step time.
  292.    * @param t current time
  293.    */
  294.   public void storeTime(final double t) {

  296.     globalCurrentTime = t;
  297.     softCurrentTime   = globalCurrentTime;
  298.     h                 = globalCurrentTime - globalPreviousTime;
  299.     setInterpolatedTime(t);

  301.     // the step is not finalized anymore
  302.     finalized  = false;
  303.   }

  305.   /** Restrict step range to a limited part of the global step.
  306.    * <p>
  307.    * This method can be used to restrict a step and make it appear
  308.    * as if the original step was smaller. Calling this method
  309.    * <em>only</em> changes the value returned by {@link #getPreviousTime()},
  310.    * it does not change any other property
  311.    * </p>
  312.    * @param softPreviousTime start of the restricted step
  313.    * @since 2.2
  314.    */
  315.   public void setSoftPreviousTime(final double softPreviousTime) {
  316.       this.softPreviousTime = softPreviousTime;
  317.   }

  319.   /** Restrict step range to a limited part of the global step.
  320.    * <p>
  321.    * This method can be used to restrict a step and make it appear
  322.    * as if the original step was smaller. Calling this method
  323.    * <em>only</em> changes the value returned by {@link #getCurrentTime()},
  324.    * it does not change any other property
  325.    * </p>
  326.    * @param softCurrentTime end of the restricted step
  327.    * @since 2.2
  328.    */
  329.   public void setSoftCurrentTime(final double softCurrentTime) {
  330.       this.softCurrentTime  = softCurrentTime;
  331.   }

  333.   /**
  334.    * Get the previous global grid point time.
  335.    * @return previous global grid point time
  336.    */
  337.   public double getGlobalPreviousTime() {
  338.     return globalPreviousTime;
  339.   }

  341.   /**
  342.    * Get the current global grid point time.
  343.    * @return current global grid point time
  344.    */
  345.   public double getGlobalCurrentTime() {
  346.     return globalCurrentTime;
  347.   }

  349.   /**
  350.    * Get the previous soft grid point time.
  351.    * @return previous soft grid point time
  352.    * @see #setSoftPreviousTime(double)
  353.    */
  354.   @Override
  355. public double getPreviousTime() {
  356.     return softPreviousTime;
  357.   }

  359.   /**
  360.    * Get the current soft grid point time.
  361.    * @return current soft grid point time
  362.    * @see #setSoftCurrentTime(double)
  363.    */
  364.   @Override
  365. public double getCurrentTime() {
  366.     return softCurrentTime;
  367.   }

  369.   /** {@inheritDoc} */
  370.   @Override
  371.   public double getInterpolatedTime() {
  372.     return interpolatedTime;
  373.   }

  375.   /** {@inheritDoc} */
  376.   @Override
  377.   public void setInterpolatedTime(final double time) {
  378.       interpolatedTime = time;
  379.       dirtyState       = true;
  380.   }

  382.   /** {@inheritDoc} */
  383.   @Override
  384.   public boolean isForward() {
  385.     return forward;
  386.   }

  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;

  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.   }

  415.   /** {@inheritDoc} */
  416.   @Override
  417.   public double[] getInterpolatedState() throws MaxCountExceededException {
  418.       evaluateCompleteInterpolatedState();
  419.       primaryMapper.extractEquationData(interpolatedState,
  420.                                         interpolatedPrimaryState);
  421.       return interpolatedPrimaryState;
  422.   }

  424.   /** {@inheritDoc} */
  425.   @Override
  426.   public double[] getInterpolatedDerivatives() throws MaxCountExceededException {
  427.       evaluateCompleteInterpolatedState();
  428.       primaryMapper.extractEquationData(interpolatedDerivatives,
  429.                                         interpolatedPrimaryDerivatives);
  430.       return interpolatedPrimaryDerivatives;
  431.   }

  433.   /** {@inheritDoc} */
  434.   @Override
  435.   public double[] getInterpolatedSecondaryState(final int index) throws MaxCountExceededException {
  436.       evaluateCompleteInterpolatedState();
  437.       secondaryMappers[index].extractEquationData(interpolatedState,
  438.                                                   interpolatedSecondaryState[index]);
  439.       return interpolatedSecondaryState[index];
  440.   }

  442.   /** {@inheritDoc} */
  443.   @Override
  444.   public double[] getInterpolatedSecondaryDerivatives(final int index) throws MaxCountExceededException {
  445.       evaluateCompleteInterpolatedState();
  446.       secondaryMappers[index].extractEquationData(interpolatedDerivatives,
  447.                                                   interpolatedSecondaryDerivatives[index]);
  448.       return interpolatedSecondaryDerivatives[index];
  449.   }

  451.   /**
  452.    * Finalize the step.

  454.    * <p>Some embedded Runge-Kutta integrators need fewer functions
  455.    * evaluations than their counterpart step interpolators. These
  456.    * interpolators should perform the last evaluations they need by
  457.    * themselves only if they need them. This method triggers these
  458.    * extra evaluations. It can be called directly by the user step
  459.    * handler and it is called automatically if {@link
  460.    * #setInterpolatedTime} is called.</p>

  462.    * <p>Once this method has been called, <strong>no</strong> other
  463.    * evaluation will be performed on this step. If there is a need to
  464.    * have some side effects between the step handler and the
  465.    * differential equations (for example update some data in the
  466.    * equations once the step has been done), it is advised to call
  467.    * this method explicitly from the step handler before these side
  468.    * effects are set up. If the step handler induces no side effect,
  469.    * then this method can safely be ignored, it will be called
  470.    * transparently as needed.</p>

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

  480.    * <p>This method calls the protected <code>doFinalize</code> method
  481.    * if it has never been called during this step and set a flag
  482.    * indicating that it has been called once. It is the <code>
  483.    * doFinalize</code> method which should perform the evaluations.
  484.    * This wrapping prevents from calling <code>doFinalize</code> several
  485.    * times and hence evaluating the differential equations too often.
  486.    * Therefore, subclasses are not allowed not reimplement it, they
  487.    * should rather reimplement <code>doFinalize</code>.</p>

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

  491.    */
  492.   public final void finalizeStep() throws MaxCountExceededException {
  493.     if (! finalized) {
  494.       doFinalize();
  495.       finalized = true;
  496.     }
  497.   }

  499.   /**
  500.    * Really finalize the step.
  501.    * The default implementation of this method does nothing.
  502.    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  503.    */
  504.   protected void doFinalize() throws MaxCountExceededException {
  505.   }

  507.   /** {@inheritDoc} */
  508.   @Override
  509.   public abstract void writeExternal(ObjectOutput out)
  510.     throws IOException;

  512.   /** {@inheritDoc} */
  513.   @Override
  514.   public abstract void readExternal(ObjectInput in)
  515.     throws IOException, ClassNotFoundException;

  517.   /** Save the base state of the instance.
  518.    * This method performs step finalization if it has not been done
  519.    * before.
  520.    * @param out stream where to save the state
  521.    * @exception IOException in case of write error
  522.    */
  523.   protected void writeBaseExternal(final ObjectOutput out)
  524.     throws IOException {

  526.     if (currentState == null) {
  527.         out.writeInt(-1);
  528.     } else {
  529.         out.writeInt(currentState.length);
  530.     }
  531.     out.writeDouble(globalPreviousTime);
  532.     out.writeDouble(globalCurrentTime);
  533.     out.writeDouble(softPreviousTime);
  534.     out.writeDouble(softCurrentTime);
  535.     out.writeDouble(h);
  536.     out.writeBoolean(forward);
  537.     out.writeObject(primaryMapper);
  538.     out.write(secondaryMappers.length);
  539.     for (final EquationsMapper  mapper : secondaryMappers) {
  540.         out.writeObject(mapper);
  541.     }

  543.     if (currentState != null) {
  544.         for (int i = 0; i < currentState.length; ++i) {
  545.             out.writeDouble(currentState[i]);
  546.         }
  547.     }

  549.     out.writeDouble(interpolatedTime);

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

  554.     try {
  555.         // finalize the step (and don't bother saving the now true flag)
  556.         finalizeStep();
  557.     } catch (MaxCountExceededException mcee) {
  558.         final IOException ioe = new IOException(mcee.getLocalizedMessage());
  559.         ioe.initCause(mcee);
  560.         throw ioe;
  561.     }
  562.   }

  564.   /** Read the base state of the instance.
  565.    * This method does <strong>neither</strong> set the interpolated
  566.    * time nor state. It is up to the derived class to reset it
  567.    * properly calling the {@link #setInterpolatedTime} method later,
  568.    * once all rest of the object state has been set up properly.
  569.    * @param in stream where to read the state from
  570.    * @return interpolated time to be set later by the caller
  571.    * @exception IOException in case of read error
  572.    * @exception ClassNotFoundException if an equation mapper class
  573.    * cannot be found
  574.    */
  575.   protected double readBaseExternal(final ObjectInput in)
  576.     throws IOException, ClassNotFoundException {

  578.     final int dimension = in.readInt();
  579.     globalPreviousTime  = in.readDouble();
  580.     globalCurrentTime   = in.readDouble();
  581.     softPreviousTime    = in.readDouble();
  582.     softCurrentTime     = in.readDouble();
  583.     h                   = in.readDouble();
  584.     forward             = in.readBoolean();
  585.     primaryMapper       = (EquationsMapper) in.readObject();
  586.     secondaryMappers    = new EquationsMapper[in.read()];
  587.     for (int i = 0; i < secondaryMappers.length; ++i) {
  588.         secondaryMappers[i] = (EquationsMapper) in.readObject();
  589.     }
  590.     dirtyState          = true;

  592.     if (dimension < 0) {
  593.         currentState = null;
  594.     } else {
  595.         currentState  = new double[dimension];
  596.         for (int i = 0; i < currentState.length; ++i) {
  597.             currentState[i] = in.readDouble();
  598.         }
  599.     }

  601.     // we do NOT handle the interpolated time and state here
  602.     interpolatedTime = Double.NaN;
  603.     allocateInterpolatedArrays(dimension);

  605.     finalized = true;

  607.     return in.readDouble();
  608.   }
  609. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT