DenseOutputModel.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.  *      https://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. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */

  23. package org.hipparchus.ode;

  25. import java.io.Serializable;
  26. import java.util.ArrayList;
  27. import java.util.List;

  29. import org.hipparchus.exception.LocalizedCoreFormats;
  30. import org.hipparchus.exception.MathIllegalArgumentException;
  31. import org.hipparchus.exception.MathIllegalStateException;
  32. import org.hipparchus.ode.sampling.ODEStateInterpolator;
  33. import org.hipparchus.ode.sampling.ODEStepHandler;
  34. import org.hipparchus.util.FastMath;

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

  92. public class DenseOutputModel implements ODEStepHandler, Serializable {

  94.     /** Serializable version identifier */
  95.     private static final long serialVersionUID = 20160328L;

  97.     /** Initial integration time. */
  98.     private double initialTime;

  100.     /** Final integration time. */
  101.     private double finalTime;

  103.     /** Integration direction indicator. */
  104.     private boolean forward;

  106.     /** Current interpolator index. */
  107.     private int index;

  109.     /** Steps table. */
  110.     private List<ODEStateInterpolator> steps;

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

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

  134.         if (model.steps.isEmpty()) {
  135.             return;
  136.         }

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

  143.             final ODEStateAndDerivative s1 = steps.get(0).getPreviousState();
  144.             final ODEStateAndDerivative s2 = model.steps.get(0).getPreviousState();
  145.             checkDimensionsEquality(s1.getPrimaryStateDimension(), s2.getPrimaryStateDimension());
  146.             checkDimensionsEquality(s1.getNumberOfSecondaryStates(), s2.getNumberOfSecondaryStates());
  147.             for (int i = 0; i < s1.getNumberOfSecondaryStates(); ++i) {
  148.                 checkDimensionsEquality(s1.getSecondaryStateDimension(i), s2.getSecondaryStateDimension(i));
  149.             }

  151.             if (forward ^ model.forward) {
  152.                 throw new MathIllegalArgumentException(LocalizedODEFormats.PROPAGATION_DIRECTION_MISMATCH);
  153.             }

  155.             final ODEStateInterpolator lastInterpolator = steps.get(index);
  156.             final double current  = lastInterpolator.getCurrentState().getTime();
  157.             final double previous = lastInterpolator.getPreviousState().getTime();
  158.             final double step = current - previous;
  159.             final double gap = model.getInitialTime() - current;
  160.             if (FastMath.abs(gap) > 1.0e-3 * FastMath.abs(step)) {
  161.                 throw new MathIllegalArgumentException(LocalizedODEFormats.HOLE_BETWEEN_MODELS_TIME_RANGES,
  162.                                                        FastMath.abs(gap));
  163.             }

  165.         }

  167.         for (ODEStateInterpolator interpolator : model.steps) {
  168.             steps.add(interpolator);
  169.         }

  171.         index = steps.size() - 1;
  172.         finalTime = (steps.get(index)).getCurrentState().getTime();

  174.     }

  176.     /** Check dimensions equality.
  177.      * @param d1 first dimension
  178.      * @param d2 second dimansion
  179.      * @exception MathIllegalArgumentException if dimensions do not match
  180.      */
  181.     private void checkDimensionsEquality(final int d1, final int d2)
  182.         throws MathIllegalArgumentException {
  183.         if (d1 != d2) {
  184.             throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
  185.                                                    d2, d1);
  186.         }
  187.     }

  189.     /** {@inheritDoc} */
  190.     @Override
  191.     public void init(final ODEStateAndDerivative initialState, final double targetTime) {
  192.         initialTime    = initialState.getTime();
  193.         this.finalTime = targetTime;
  194.         forward        = true;
  195.         index          = 0;
  196.         steps.clear();
  197.     }

  199.     /** {@inheritDoc} */
  200.     @Override
  201.     public void handleStep(final ODEStateInterpolator interpolator) {
  202.         if (steps.isEmpty()) {
  203.             initialTime = interpolator.getPreviousState().getTime();
  204.             forward     = interpolator.isForward();
  205.         }
  206.         steps.add(interpolator);
  207.     }

  209.     /** {@inheritDoc} */
  210.     @Override
  211.     public void finish(final ODEStateAndDerivative finalState) {
  212.         finalTime = finalState.getTime();
  213.         index     = steps.size() - 1;
  214.     }

  216.     /**
  217.      * Get the initial integration time.
  218.      * @return initial integration time
  219.      */
  220.     public double getInitialTime() {
  221.         return initialTime;
  222.     }

  224.     /**
  225.      * Get the final integration time.
  226.      * @return final integration time
  227.      */
  228.     public double getFinalTime() {
  229.         return finalTime;
  230.     }

  232.     /**
  233.      * Get the state at interpolated time.
  234.      * @param time time of the interpolated point
  235.      * @return state at interpolated time
  236.      */
  237.     public ODEStateAndDerivative getInterpolatedState(final double time) {

  239.         // initialize the search with the complete steps table
  240.         int iMin = 0;
  241.         final ODEStateInterpolator sMin = steps.get(iMin);
  242.         double tMin = 0.5 * (sMin.getPreviousState().getTime() + sMin.getCurrentState().getTime());

  244.         int iMax = steps.size() - 1;
  245.         final ODEStateInterpolator sMax = steps.get(iMax);
  246.         double tMax = 0.5 * (sMax.getPreviousState().getTime() + sMax.getCurrentState().getTime());

  248.         // handle points outside of the integration interval
  249.         // or in the first and last step
  250.         if (locatePoint(time, sMin) <= 0) {
  251.             index = iMin;
  252.             return sMin.getInterpolatedState(time);
  253.         }
  254.         if (locatePoint(time, sMax) >= 0) {
  255.             index = iMax;
  256.             return sMax.getInterpolatedState(time);
  257.         }

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

  262.             // use the last estimated index as the splitting index
  263.             final ODEStateInterpolator si = steps.get(index);
  264.             final int location = locatePoint(time, si);
  265.             if (location < 0) {
  266.                 iMax = index;
  267.                 tMax = 0.5 * (si.getPreviousState().getTime() + si.getCurrentState().getTime());
  268.             } else if (location > 0) {
  269.                 iMin = index;
  270.                 tMin = 0.5 * (si.getPreviousState().getTime() + si.getCurrentState().getTime());
  271.             } else {
  272.                 // we have found the target step, no need to continue searching
  273.                 return si.getInterpolatedState(time);
  274.             }

  276.             // compute a new estimate of the index in the reduced table slice
  277.             final int iMed = (iMin + iMax) / 2;
  278.             final ODEStateInterpolator sMed = steps.get(iMed);
  279.             final double tMed = 0.5 * (sMed.getPreviousState().getTime() + sMed.getCurrentState().getTime());

  281.             if ((FastMath.abs(tMed - tMin) < 1e-6) || (FastMath.abs(tMax - tMed) < 1e-6)) {
  282.                 // too close to the bounds, we estimate using a simple dichotomy
  283.                 index = iMed;
  284.             } else {
  285.                 // estimate the index using a reverse quadratic polynom
  286.                 // (reverse means we have i = P(t), thus allowing to simply
  287.                 // compute index = P(time) rather than solving a quadratic equation)
  288.                 final double d12 = tMax - tMed;
  289.                 final double d23 = tMed - tMin;
  290.                 final double d13 = tMax - tMin;
  291.                 final double dt1 = time - tMax;
  292.                 final double dt2 = time - tMed;
  293.                 final double dt3 = time - tMin;
  294.                 final double iLagrange = ((dt2 * dt3 * d23) * iMax -
  295.                                 (dt1 * dt3 * d13) * iMed +
  296.                                 (dt1 * dt2 * d12) * iMin) /
  297.                                 (d12 * d23 * d13);
  298.                 index = (int) FastMath.rint(iLagrange);
  299.             }

  301.             // force the next size reduction to be at least one tenth
  302.             final int low  = FastMath.max(iMin + 1, (9 * iMin + iMax) / 10);
  303.             final int high = FastMath.min(iMax - 1, (iMin + 9 * iMax) / 10);
  304.             if (index < low) {
  305.                 index = low;
  306.             } else if (index > high) {
  307.                 index = high;
  308.             }

  310.         }

  312.         // now the table slice is very small, we perform an iterative search
  313.         index = iMin;
  314.         while ((index <= iMax) && (locatePoint(time, steps.get(index)) > 0)) {
  315.             ++index;
  316.         }

  318.         return steps.get(index).getInterpolatedState(time);

  320.     }

  322.     /** Compare a step interval and a double.
  323.      * @param time point to locate
  324.      * @param interval step interval
  325.      * @return -1 if the double is before the interval, 0 if it is in
  326.      * the interval, and +1 if it is after the interval, according to
  327.      * the interval direction
  328.      */
  329.     private int locatePoint(final double time, final ODEStateInterpolator interval) {
  330.         if (forward) {
  331.             if (time < interval.getPreviousState().getTime()) {
  332.                 return -1;
  333.             } else if (time > interval.getCurrentState().getTime()) {
  334.                 return +1;
  335.             } else {
  336.                 return 0;
  337.             }
  338.         }
  339.         if (time > interval.getPreviousState().getTime()) {
  340.             return -1;
  341.         } else if (time < interval.getCurrentState().getTime()) {
  342.             return +1;
  343.         } else {
  344.             return 0;
  345.         }
  346.     }

  348. }
Created with JaCoCo 0.8.11.202310140853