ContinuousOutputFieldModel.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.util.ArrayList;
  21. import java.util.List;

  23. import org.apache.commons.math4.legacy.core.RealFieldElement;
  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.FieldStepHandler;
  29. import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator;
  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 #getInterpolatedState(RealFieldElement) getInterpolatedState}
  42.  * method to retrieve this information at any time. It is important to wait
  43.  * for the integration to be over before attempting to call {@link
  44.  * #getInterpolatedState(RealFieldElement)} because some internal
  45.  * variables are set only once the last step has been handled.</p>
  46.  *
  47.  * <p>This is useful for example if the main loop of the user
  48.  * application should remain independent from the integration process
  49.  * or if one needs to mimic the behaviour of an analytical model
  50.  * despite a numerical model is used (i.e. one needs the ability to
  51.  * get the model value at any time or to navigate through the
  52.  * data).</p>
  53.  *
  54.  * <p>If problem modeling is done with several separate
  55.  * integration phases for contiguous intervals, the same
  56.  * ContinuousOutputModel can be used as step handler for all
  57.  * integration phases as long as they are performed in order and in
  58.  * the same direction. As an example, one can extrapolate the
  59.  * trajectory of a satellite with one model (i.e. one set of
  60.  * differential equations) up to the beginning of a maneuver, use
  61.  * another more complex model including thrusters modeling and
  62.  * accurate attitude control during the maneuver, and revert to the
  63.  * first model after the end of the maneuver. If the same continuous
  64.  * output model handles the steps of all integration phases, the user
  65.  * do not need to bother when the maneuver begins or ends, he has all
  66.  * the data available in a transparent manner.</p>
  67.  *
  68.  * <p>One should be aware that the amount of data stored in a
  69.  * ContinuousOutputFieldModel instance can be important if the state vector
  70.  * is large, if the integration interval is long or if the steps are
  71.  * small (which can result from small tolerance settings in {@link
  72.  * org.apache.commons.math4.legacy.ode.nonstiff.AdaptiveStepsizeFieldIntegrator adaptive
  73.  * step size integrators}).</p>
  74.  *
  75.  * @see FieldStepHandler
  76.  * @see FieldStepInterpolator
  77.  * @param <T> the type of the field elements
  78.  * @since 3.6
  79.  */

  81. public class ContinuousOutputFieldModel<T extends RealFieldElement<T>>
  82.     implements FieldStepHandler<T> {

  84.     /** Initial integration time. */
  85.     private T initialTime;

  87.     /** Final integration time. */
  88.     private T finalTime;

  90.     /** Integration direction indicator. */
  91.     private boolean forward;

  93.     /** Current interpolator index. */
  94.     private int index;

  96.     /** Steps table. */
  97.     private List<FieldStepInterpolator<T>> steps;

  99.     /** Simple constructor.
  100.      * Build an empty continuous output model.
  101.      */
  102.     public ContinuousOutputFieldModel() {
  103.         steps       = new ArrayList<>();
  104.         initialTime = null;
  105.         finalTime   = null;
  106.         forward     = true;
  107.         index       = 0;
  108.     }

  110.     /** Append another model at the end of the instance.
  111.      * @param model model to add at the end of the instance
  112.      * @exception MathIllegalArgumentException if the model to append is not
  113.      * compatible with the instance (dimension of the state vector,
  114.      * propagation direction, hole between the dates)
  115.      * @exception DimensionMismatchException if the dimensions of the states or
  116.      * the number of secondary states do not match
  117.      * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  118.      * during step finalization
  119.      */
  120.     public void append(final ContinuousOutputFieldModel<T> model)
  121.         throws MathIllegalArgumentException, MaxCountExceededException {

  123.         if (model.steps.isEmpty()) {
  124.             return;
  125.         }

  127.         if (steps.isEmpty()) {
  128.             initialTime = model.initialTime;
  129.             forward     = model.forward;
  130.         } else {

  132.             // safety checks
  133.             final FieldODEStateAndDerivative<T> s1 = steps.get(0).getPreviousState();
  134.             final FieldODEStateAndDerivative<T> s2 = model.steps.get(0).getPreviousState();
  135.             checkDimensionsEquality(s1.getStateDimension(), s2.getStateDimension());
  136.             checkDimensionsEquality(s1.getNumberOfSecondaryStates(), s2.getNumberOfSecondaryStates());
  137.             for (int i = 0; i < s1.getNumberOfSecondaryStates(); ++i) {
  138.                 checkDimensionsEquality(s1.getSecondaryStateDimension(i), s2.getSecondaryStateDimension(i));
  139.             }

  141.             if (forward ^ model.forward) {
  142.                 throw new MathIllegalArgumentException(LocalizedFormats.PROPAGATION_DIRECTION_MISMATCH);
  143.             }

  145.             final FieldStepInterpolator<T> lastInterpolator = steps.get(index);
  146.             final T current  = lastInterpolator.getCurrentState().getTime();
  147.             final T previous = lastInterpolator.getPreviousState().getTime();
  148.             final T step = current.subtract(previous);
  149.             final T gap = model.getInitialTime().subtract(current);
  150.             if (gap.abs().subtract(step.abs().multiply(1.0e-3)).getReal() > 0) {
  151.                 throw new MathIllegalArgumentException(LocalizedFormats.HOLE_BETWEEN_MODELS_TIME_RANGES,
  152.                                                        gap.abs().getReal());
  153.             }
  154.         }

  156.         steps.addAll(model.steps);

  158.         index = steps.size() - 1;
  159.         finalTime = (steps.get(index)).getCurrentState().getTime();
  160.     }

  162.     /** Check dimensions equality.
  163.      * @param d1 first dimension
  164.      * @param d2 second dimension
  165.      * @exception DimensionMismatchException if dimensions do not match
  166.      */
  167.     private void checkDimensionsEquality(final int d1, final int d2)
  168.         throws DimensionMismatchException {
  169.         if (d1 != d2) {
  170.             throw new DimensionMismatchException(d2, d1);
  171.         }
  172.     }

  174.     /** {@inheritDoc} */
  175.     @Override
  176.     public void init(final FieldODEStateAndDerivative<T> initialState, final T t) {
  177.         initialTime = initialState.getTime();
  178.         finalTime   = t;
  179.         forward     = true;
  180.         index       = 0;
  181.         steps.clear();
  182.     }

  184.     /** Handle the last accepted step.
  185.      * A copy of the information provided by the last step is stored in
  186.      * the instance for later use.
  187.      * @param interpolator interpolator for the last accepted step.
  188.      * @param isLast true if the step is the last one
  189.      * @exception MaxCountExceededException if the number of functions evaluations is exceeded
  190.      * during step finalization
  191.      */
  192.     @Override
  193.     public void handleStep(final FieldStepInterpolator<T> interpolator, final boolean isLast)
  194.         throws MaxCountExceededException {

  196.         if (steps.isEmpty()) {
  197.             initialTime = interpolator.getPreviousState().getTime();
  198.             forward     = interpolator.isForward();
  199.         }

  201.         steps.add(interpolator);

  203.         if (isLast) {
  204.             finalTime = interpolator.getCurrentState().getTime();
  205.             index     = steps.size() - 1;
  206.         }
  207.     }

  209.     /**
  210.      * Get the initial integration time.
  211.      * @return initial integration time
  212.      */
  213.     public T getInitialTime() {
  214.         return initialTime;
  215.     }

  217.     /**
  218.      * Get the final integration time.
  219.      * @return final integration time
  220.      */
  221.     public T getFinalTime() {
  222.         return finalTime;
  223.     }

  225.     /**
  226.      * Get the state at interpolated time.
  227.      * @param time time of the interpolated point
  228.      * @return state at interpolated time
  229.      */
  230.     public FieldODEStateAndDerivative<T> getInterpolatedState(final T time) {

  232.         // initialize the search with the complete steps table
  233.         int iMin = 0;
  234.         final FieldStepInterpolator<T> sMin = steps.get(iMin);
  235.         T tMin = sMin.getPreviousState().getTime().add(sMin.getCurrentState().getTime()).multiply(0.5);

  237.         int iMax = steps.size() - 1;
  238.         final FieldStepInterpolator<T> sMax = steps.get(iMax);
  239.         T tMax = sMax.getPreviousState().getTime().add(sMax.getCurrentState().getTime()).multiply(0.5);

  241.         // handle points outside of the integration interval
  242.         // or in the first and last step
  243.         if (locatePoint(time, sMin) <= 0) {
  244.             index = iMin;
  245.             return sMin.getInterpolatedState(time);
  246.         }
  247.         if (locatePoint(time, sMax) >= 0) {
  248.             index = iMax;
  249.             return sMax.getInterpolatedState(time);
  250.         }

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

  255.             // use the last estimated index as the splitting index
  256.             final FieldStepInterpolator<T> si = steps.get(index);
  257.             final int location = locatePoint(time, si);
  258.             if (location < 0) {
  259.                 iMax = index;
  260.                 tMax = si.getPreviousState().getTime().add(si.getCurrentState().getTime()).multiply(0.5);
  261.             } else if (location > 0) {
  262.                 iMin = index;
  263.                 tMin = si.getPreviousState().getTime().add(si.getCurrentState().getTime()).multiply(0.5);
  264.             } else {
  265.                 // we have found the target step, no need to continue searching
  266.                 return si.getInterpolatedState(time);
  267.             }

  269.             // compute a new estimate of the index in the reduced table slice
  270.             final int iMed = (iMin + iMax) / 2;
  271.             final FieldStepInterpolator<T> sMed = steps.get(iMed);
  272.             final T tMed = sMed.getPreviousState().getTime().add(sMed.getCurrentState().getTime()).multiply(0.5);

  274.             if (tMed.subtract(tMin).abs().subtract(1.0e-6).getReal() < 0 ||
  275.                 tMax.subtract(tMed).abs().subtract(1.0e-6).getReal() < 0) {
  276.                 // too close to the bounds, we estimate using a simple dichotomy
  277.                 index = iMed;
  278.             } else {
  279.                 // estimate the index using a reverse quadratic polynomial
  280.                 // (reverse means we have i = P(t), thus allowing to simply
  281.                 // compute index = P(time) rather than solving a quadratic equation)
  282.                 final T d12 = tMax.subtract(tMed);
  283.                 final T d23 = tMed.subtract(tMin);
  284.                 final T d13 = tMax.subtract(tMin);
  285.                 final T dt1 = time.subtract(tMax);
  286.                 final T dt2 = time.subtract(tMed);
  287.                 final T dt3 = time.subtract(tMin);
  288.                 final T iLagrange =           dt2.multiply(dt3).multiply(d23).multiply(iMax).
  289.                                      subtract(dt1.multiply(dt3).multiply(d13).multiply(iMed)).
  290.                                      add(     dt1.multiply(dt2).multiply(d12).multiply(iMin)).
  291.                                      divide(d12.multiply(d23).multiply(d13));
  292.                 index = (int) JdkMath.rint(iLagrange.getReal());
  293.             }

  295.             // force the next size reduction to be at least one tenth
  296.             final int low  = JdkMath.max(iMin + 1, (9 * iMin + iMax) / 10);
  297.             final int high = JdkMath.min(iMax - 1, (iMin + 9 * iMax) / 10);
  298.             if (index < low) {
  299.                 index = low;
  300.             } else if (index > high) {
  301.                 index = high;
  302.             }
  303.         }

  305.         // now the table slice is very small, we perform an iterative search
  306.         index = iMin;
  307.         while (index <= iMax && locatePoint(time, steps.get(index)) > 0) {
  308.             ++index;
  309.         }

  311.         return steps.get(index).getInterpolatedState(time);
  312.     }

  314.     /** Compare a step interval and a double.
  315.      * @param time point to locate
  316.      * @param interval step interval
  317.      * @return -1 if the double is before the interval, 0 if it is in
  318.      * the interval, and +1 if it is after the interval, according to
  319.      * the interval direction
  320.      */
  321.     private int locatePoint(final T time, final FieldStepInterpolator<T> interval) {
  322.         if (forward) {
  323.             if (time.subtract(interval.getPreviousState().getTime()).getReal() < 0) {
  324.                 return -1;
  325.             } else if (time.subtract(interval.getCurrentState().getTime()).getReal() > 0) {
  326.                 return +1;
  327.             } else {
  328.                 return 0;
  329.             }
  330.         }
  331.         if (time.subtract(interval.getPreviousState().getTime()).getReal() > 0) {
  332.             return -1;
  333.         } else if (time.subtract(interval.getCurrentState().getTime()).getReal() < 0) {
  334.             return +1;
  335.         } else {
  336.             return 0;
  337.         }
  338.     }
  339. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT