ContinuousOutputModel.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;

  20. import java.io.Serializable;
  21. import java.util.ArrayList;
  22. import java.util.List;

  24. import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  25. import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
  26. import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  27. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  28. import org.apache.commons.math4.legacy.ode.sampling.StepHandler;
  29. import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
  30. import org.apache.commons.math4.core.jdkmath.JdkMath;

  32. /**
  33.  * This class stores all information provided by an ODE integrator
  34.  * during the integration process and build a continuous model of the
  35.  * solution from this.
  36.  *
  37.  * <p>This class act as a step handler from the integrator point of
  38.  * view. It is called iteratively during the integration process and
  39.  * stores a copy of all steps information in a sorted collection for
  40.  * later use. Once the integration process is over, the user can use
  41.  * the {@link #setInterpolatedTime setInterpolatedTime} and {@link
  42.  * #getInterpolatedState getInterpolatedState} to retrieve this
  43.  * information at any time. It is important to wait for the
  44.  * integration to be over before attempting to call {@link
  45.  * #setInterpolatedTime setInterpolatedTime} because some internal
  46.  * variables are set only once the last step has been handled.</p>
  47.  *
  48.  * <p>This is useful for example if the main loop of the user
  49.  * application should remain independent from the integration process
  50.  * or if one needs to mimic the behaviour of an analytical model
  51.  * despite a numerical model is used (i.e. one needs the ability to
  52.  * get the model value at any time or to navigate through the
  53.  * data).</p>
  54.  *
  55.  * <p>If problem modeling is done with several separate
  56.  * integration phases for contiguous intervals, the same
  57.  * ContinuousOutputModel can be used as step handler for all
  58.  * integration phases as long as they are performed in order and in
  59.  * the same direction. As an example, one can extrapolate the
  60.  * trajectory of a satellite with one model (i.e. one set of
  61.  * differential equations) up to the beginning of a maneuver, use
  62.  * another more complex model including thrusters modeling and
  63.  * accurate attitude control during the maneuver, and revert to the
  64.  * first model after the end of the maneuver. If the same continuous
  65.  * output model handles the steps of all integration phases, the user
  66.  * do not need to bother when the maneuver begins or ends, he has all
  67.  * the data available in a transparent manner.</p>
  68.  *
  69.  * <p>An important feature of this class is that it implements the
  70.  * <code>Serializable</code> interface. This means that the result of
  71.  * an integration can be serialized and reused later (if stored into a
  72.  * persistent medium like a file system or a database) or elsewhere (if
  73.  * sent to another application). Only the result of the integration is
  74.  * stored, there is no reference to the integrated problem by
  75.  * itself.</p>
  76.  *
  77.  * <p>One should be aware that the amount of data stored in a
  78.  * ContinuousOutputModel instance can be important if the state vector
  79.  * is large, if the integration interval is long or if the steps are
  80.  * small (which can result from small tolerance settings in {@link
  81.  * org.apache.commons.math4.legacy.ode.nonstiff.AdaptiveStepsizeIntegrator adaptive
  82.  * step size integrators}).</p>
  83.  *
  84.  * @see StepHandler
  85.  * @see StepInterpolator
  86.  * @since 1.2
  87.  */

  89. public class ContinuousOutputModel
  90.   implements StepHandler, Serializable {

  92.     /** Serializable version identifier. */
  93.     private static final long serialVersionUID = -1417964919405031606L;

  95.     /** Initial integration time. */
  96.     private double initialTime;

  98.     /** Final integration time. */
  99.     private double finalTime;

  101.     /** Integration direction indicator. */
  102.     private boolean forward;

  104.     /** Current interpolator index. */
  105.     private int index;

  107.     /** Steps table. */
  108.     private List<StepInterpolator> steps;

  110.   /** Simple constructor.
  111.    * Build an empty continuous output model.
  112.    */
  113.   public ContinuousOutputModel() {
  114.     steps = new ArrayList<>();
  115.     initialTime = Double.NaN;
  116.     finalTime   = Double.NaN;
  117.     forward     = true;
  118.     index       = 0;
  119.   }

  121.   /** Append another model at the end of the instance.
  122.    * @param model model to add at the end of the instance
  123.    * @exception MathIllegalArgumentException if the model to append is not
  124.    * compatible with the instance (dimension of the state vector,
  125.    * propagation direction, hole between the dates)
  126.    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  127.    * during step finalization
  128.    */
  129.   public void append(final ContinuousOutputModel model)
  130.     throws MathIllegalArgumentException, MaxCountExceededException {

  132.     if (model.steps.isEmpty()) {
  133.       return;
  134.     }

  136.     if (steps.isEmpty()) {
  137.       initialTime = model.initialTime;
  138.       forward     = model.forward;
  139.     } else {

  141.       if (getInterpolatedState().length != model.getInterpolatedState().length) {
  142.           throw new DimensionMismatchException(model.getInterpolatedState().length,
  143.                                                getInterpolatedState().length);
  144.       }

  146.       if (forward ^ model.forward) {
  147.           throw new MathIllegalArgumentException(LocalizedFormats.PROPAGATION_DIRECTION_MISMATCH);
  148.       }

  150.       final StepInterpolator lastInterpolator = steps.get(index);
  151.       final double current  = lastInterpolator.getCurrentTime();
  152.       final double previous = lastInterpolator.getPreviousTime();
  153.       final double step = current - previous;
  154.       final double gap = model.getInitialTime() - current;
  155.       if (JdkMath.abs(gap) > 1.0e-3 * JdkMath.abs(step)) {
  156.         throw new MathIllegalArgumentException(LocalizedFormats.HOLE_BETWEEN_MODELS_TIME_RANGES,
  157.                                                JdkMath.abs(gap));
  158.       }
  159.     }

  161.     for (StepInterpolator interpolator : model.steps) {
  162.       steps.add(interpolator.copy());
  163.     }

  165.     index = steps.size() - 1;
  166.     finalTime = (steps.get(index)).getCurrentTime();
  167.   }

  169.   /** {@inheritDoc} */
  170.   @Override
  171.   public void init(double t0, double[] y0, double t) {
  172.     initialTime = Double.NaN;
  173.     finalTime   = Double.NaN;
  174.     forward     = true;
  175.     index       = 0;
  176.     steps.clear();
  177.   }

  179.   /** Handle the last accepted step.
  180.    * A copy of the information provided by the last step is stored in
  181.    * the instance for later use.
  182.    * @param interpolator interpolator for the last accepted step.
  183.    * @param isLast true if the step is the last one
  184.    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  185.    * during step finalization
  186.    */
  187.   @Override
  188. public void handleStep(final StepInterpolator interpolator, final boolean isLast)
  189.       throws MaxCountExceededException {

  191.     if (steps.isEmpty()) {
  192.       initialTime = interpolator.getPreviousTime();
  193.       forward     = interpolator.isForward();
  194.     }

  196.     steps.add(interpolator.copy());

  198.     if (isLast) {
  199.       finalTime = interpolator.getCurrentTime();
  200.       index     = steps.size() - 1;
  201.     }
  202.   }

  204.   /**
  205.    * Get the initial integration time.
  206.    * @return initial integration time
  207.    */
  208.   public double getInitialTime() {
  209.     return initialTime;
  210.   }

  212.   /**
  213.    * Get the final integration time.
  214.    * @return final integration time
  215.    */
  216.   public double getFinalTime() {
  217.     return finalTime;
  218.   }

  220.   /**
  221.    * Get the time of the interpolated point.
  222.    * If {@link #setInterpolatedTime} has not been called, it returns
  223.    * the final integration time.
  224.    * @return interpolation point time
  225.    */
  226.   public double getInterpolatedTime() {
  227.     return steps.get(index).getInterpolatedTime();
  228.   }

  230.   /** Set the time of the interpolated point.
  231.    * <p>This method should <strong>not</strong> be called before the
  232.    * integration is over because some internal variables are set only
  233.    * once the last step has been handled.</p>
  234.    * <p>Setting the time outside of the integration interval is now
  235.    * allowed, but should be used with care since the accuracy of the
  236.    * interpolator will probably be very poor far from this interval.
  237.    * This allowance has been added to simplify implementation of search
  238.    * algorithms near the interval endpoints.</p>
  239.    * <p>Note that each time this method is called, the internal arrays
  240.    * returned in {@link #getInterpolatedState()}, {@link
  241.    * #getInterpolatedDerivatives()} and {@link #getInterpolatedSecondaryState(int)}
  242.    * <em>will</em> be overwritten. So if their content must be preserved
  243.    * across several calls, user must copy them.</p>
  244.    * @param time time of the interpolated point
  245.    * @see #getInterpolatedState()
  246.    * @see #getInterpolatedDerivatives()
  247.    * @see #getInterpolatedSecondaryState(int)
  248.    */
  249.   public void setInterpolatedTime(final double time) {

  251.       // initialize the search with the complete steps table
  252.       int iMin = 0;
  253.       final StepInterpolator sMin = steps.get(iMin);
  254.       double tMin = 0.5 * (sMin.getPreviousTime() + sMin.getCurrentTime());

  256.       int iMax = steps.size() - 1;
  257.       final StepInterpolator sMax = steps.get(iMax);
  258.       double tMax = 0.5 * (sMax.getPreviousTime() + sMax.getCurrentTime());

  260.       // handle points outside of the integration interval
  261.       // or in the first and last step
  262.       if (locatePoint(time, sMin) <= 0) {
  263.         index = iMin;
  264.         sMin.setInterpolatedTime(time);
  265.         return;
  266.       }
  267.       if (locatePoint(time, sMax) >= 0) {
  268.         index = iMax;
  269.         sMax.setInterpolatedTime(time);
  270.         return;
  271.       }

  273.       // reduction of the table slice size
  274.       while (iMax - iMin > 5) {

  276.         // use the last estimated index as the splitting index
  277.         final StepInterpolator si = steps.get(index);
  278.         final int location = locatePoint(time, si);
  279.         if (location < 0) {
  280.           iMax = index;
  281.           tMax = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
  282.         } else if (location > 0) {
  283.           iMin = index;
  284.           tMin = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
  285.         } else {
  286.           // we have found the target step, no need to continue searching
  287.           si.setInterpolatedTime(time);
  288.           return;
  289.         }

  291.         // compute a new estimate of the index in the reduced table slice
  292.         final int iMed = (iMin + iMax) / 2;
  293.         final StepInterpolator sMed = steps.get(iMed);
  294.         final double tMed = 0.5 * (sMed.getPreviousTime() + sMed.getCurrentTime());

  296.         if (JdkMath.abs(tMed - tMin) < 1e-6 || JdkMath.abs(tMax - tMed) < 1e-6) {
  297.           // too close to the bounds, we estimate using a simple dichotomy
  298.           index = iMed;
  299.         } else {
  300.           // estimate the index using a reverse quadratic polynom
  301.           // (reverse means we have i = P(t), thus allowing to simply
  302.           // compute index = P(time) rather than solving a quadratic equation)
  303.           final double d12 = tMax - tMed;
  304.           final double d23 = tMed - tMin;
  305.           final double d13 = tMax - tMin;
  306.           final double dt1 = time - tMax;
  307.           final double dt2 = time - tMed;
  308.           final double dt3 = time - tMin;
  309.           final double iLagrange = ((dt2 * dt3 * d23) * iMax -
  310.                                     (dt1 * dt3 * d13) * iMed +
  311.                                     (dt1 * dt2 * d12) * iMin) /
  312.                                    (d12 * d23 * d13);
  313.           index = (int) JdkMath.rint(iLagrange);
  314.         }

  316.         // force the next size reduction to be at least one tenth
  317.         final int low  = JdkMath.max(iMin + 1, (9 * iMin + iMax) / 10);
  318.         final int high = JdkMath.min(iMax - 1, (iMin + 9 * iMax) / 10);
  319.         if (index < low) {
  320.           index = low;
  321.         } else if (index > high) {
  322.           index = high;
  323.         }
  324.       }

  326.       // now the table slice is very small, we perform an iterative search
  327.       index = iMin;
  328.       while (index <= iMax && locatePoint(time, steps.get(index)) > 0) {
  329.         ++index;
  330.       }

  332.       steps.get(index).setInterpolatedTime(time);
  333.   }

  335.   /**
  336.    * Get the state vector of the interpolated point.
  337.    * <p>The returned vector is a reference to a reused array, so
  338.    * it should not be modified and it should be copied if it needs
  339.    * to be preserved across several calls to the associated
  340.    * {@link #setInterpolatedTime(double)} method.</p>
  341.    * @return state vector at time {@link #getInterpolatedTime}
  342.    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  343.    * @see #setInterpolatedTime(double)
  344.    * @see #getInterpolatedDerivatives()
  345.    * @see #getInterpolatedSecondaryState(int)
  346.    * @see #getInterpolatedSecondaryDerivatives(int)
  347.    */
  348.   public double[] getInterpolatedState() throws MaxCountExceededException {
  349.     return steps.get(index).getInterpolatedState();
  350.   }

  352.   /**
  353.    * Get the derivatives of the state vector of the interpolated point.
  354.    * <p>The returned vector is a reference to a reused array, so
  355.    * it should not be modified and it should be copied if it needs
  356.    * to be preserved across several calls to the associated
  357.    * {@link #setInterpolatedTime(double)} method.</p>
  358.    * @return derivatives of the state vector at time {@link #getInterpolatedTime}
  359.    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  360.    * @see #setInterpolatedTime(double)
  361.    * @see #getInterpolatedState()
  362.    * @see #getInterpolatedSecondaryState(int)
  363.    * @see #getInterpolatedSecondaryDerivatives(int)
  364.    * @since 3.4
  365.    */
  366.   public double[] getInterpolatedDerivatives() throws MaxCountExceededException {
  367.     return steps.get(index).getInterpolatedDerivatives();
  368.   }

  370.   /** Get the interpolated secondary state corresponding to the secondary equations.
  371.    * <p>The returned vector is a reference to a reused array, so
  372.    * it should not be modified and it should be copied if it needs
  373.    * to be preserved across several calls to the associated
  374.    * {@link #setInterpolatedTime(double)} method.</p>
  375.    * @param secondaryStateIndex index of the secondary set, as returned by {@link
  376.    * org.apache.commons.math4.legacy.ode.ExpandableStatefulODE#addSecondaryEquations(
  377.    * org.apache.commons.math4.legacy.ode.SecondaryEquations)
  378.    * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
  379.    * @return interpolated secondary state at the current interpolation date
  380.    * @see #setInterpolatedTime(double)
  381.    * @see #getInterpolatedState()
  382.    * @see #getInterpolatedDerivatives()
  383.    * @see #getInterpolatedSecondaryDerivatives(int)
  384.    * @since 3.2
  385.    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  386.    */
  387.   public double[] getInterpolatedSecondaryState(final int secondaryStateIndex)
  388.     throws MaxCountExceededException {
  389.     return steps.get(index).getInterpolatedSecondaryState(secondaryStateIndex);
  390.   }

  392.   /** Get the interpolated secondary derivatives corresponding to the secondary equations.
  393.    * <p>The returned vector is a reference to a reused array, so
  394.    * it should not be modified and it should be copied if it needs
  395.    * to be preserved across several calls to the associated
  396.    * {@link #setInterpolatedTime(double)} method.</p>
  397.    * @param secondaryStateIndex index of the secondary set, as returned by {@link
  398.    * org.apache.commons.math4.legacy.ode.ExpandableStatefulODE#addSecondaryEquations(
  399.    * org.apache.commons.math4.legacy.ode.SecondaryEquations)
  400.    * ExpandableStatefulODE.addSecondaryEquations(SecondaryEquations)}
  401.    * @return interpolated secondary derivatives at the current interpolation date
  402.    * @see #setInterpolatedTime(double)
  403.    * @see #getInterpolatedState()
  404.    * @see #getInterpolatedDerivatives()
  405.    * @see #getInterpolatedSecondaryState(int)
  406.    * @since 3.4
  407.    * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  408.    */
  409.   public double[] getInterpolatedSecondaryDerivatives(final int secondaryStateIndex)
  410.     throws MaxCountExceededException {
  411.     return steps.get(index).getInterpolatedSecondaryDerivatives(secondaryStateIndex);
  412.   }

  414.   /** Compare a step interval and a double.
  415.    * @param time point to locate
  416.    * @param interval step interval
  417.    * @return -1 if the double is before the interval, 0 if it is in
  418.    * the interval, and +1 if it is after the interval, according to
  419.    * the interval direction
  420.    */
  421.   private int locatePoint(final double time, final StepInterpolator interval) {
  422.     if (forward) {
  423.       if (time < interval.getPreviousTime()) {
  424.         return -1;
  425.       } else if (time > interval.getCurrentTime()) {
  426.         return +1;
  427.       } else {
  428.         return 0;
  429.       }
  430.     }
  431.     if (time > interval.getPreviousTime()) {
  432.       return -1;
  433.     } else if (time < interval.getCurrentTime()) {
  434.       return +1;
  435.     } else {
  436.       return 0;
  437.     }
  438.   }
  439. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT