    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017
    018    package org.apache.commons.math3.ode;
    019
    020    import java.io.Serializable;
    021    import java.util.ArrayList;
    022    import java.util.List;
    023
    024    import org.apache.commons.math3.exception.DimensionMismatchException;
    025    import org.apache.commons.math3.exception.MathIllegalArgumentException;
    026    import org.apache.commons.math3.exception.MaxCountExceededException;
    027    import org.apache.commons.math3.exception.util.LocalizedFormats;
    028    import org.apache.commons.math3.ode.sampling.StepHandler;
    029    import org.apache.commons.math3.ode.sampling.StepInterpolator;
    030    import org.apache.commons.math3.util.FastMath;
    031
    032    /**
    033     * This class stores all information provided by an ODE integrator
    034     * during the integration process and build a continuous model of the
    035     * solution from this.
    036     *
    037     * <p>This class act as a step handler from the integrator point of
    038     * view. It is called iteratively during the integration process and
    039     * stores a copy of all steps information in a sorted collection for
    040     * later use. Once the integration process is over, the user can use
    041     * the {@link #setInterpolatedTime setInterpolatedTime} and {@link
    042     * #getInterpolatedState getInterpolatedState} to retrieve this
    043     * information at any time. It is important to wait for the
    044     * integration to be over before attempting to call {@link
    045     * #setInterpolatedTime setInterpolatedTime} because some internal
    046     * variables are set only once the last step has been handled.</p>
    047     *
    048     * <p>This is useful for example if the main loop of the user
    049     * application should remain independent from the integration process
    050     * or if one needs to mimic the behaviour of an analytical model
    051     * despite a numerical model is used (i.e. one needs the ability to
    052     * get the model value at any time or to navigate through the
    053     * data).</p>
    054     *
    055     * <p>If problem modeling is done with several separate
    056     * integration phases for contiguous intervals, the same
    057     * ContinuousOutputModel can be used as step handler for all
    058     * integration phases as long as they are performed in order and in
    059     * the same direction. As an example, one can extrapolate the
    060     * trajectory of a satellite with one model (i.e. one set of
    061     * differential equations) up to the beginning of a maneuver, use
    062     * another more complex model including thrusters modeling and
    063     * accurate attitude control during the maneuver, and revert to the
    064     * first model after the end of the maneuver. If the same continuous
    065     * output model handles the steps of all integration phases, the user
    066     * do not need to bother when the maneuver begins or ends, he has all
    067     * the data available in a transparent manner.</p>
    068     *
    069     * <p>An important feature of this class is that it implements the
    070     * <code>Serializable</code> interface. This means that the result of
    071     * an integration can be serialized and reused later (if stored into a
    072     * persistent medium like a filesystem or a database) or elsewhere (if
    073     * sent to another application). Only the result of the integration is
    074     * stored, there is no reference to the integrated problem by
    075     * itself.</p>
    076     *
    077     * <p>One should be aware that the amount of data stored in a
    078     * ContinuousOutputModel instance can be important if the state vector
    079     * is large, if the integration interval is long or if the steps are
    080     * small (which can result from small tolerance settings in {@link
    081     * org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeIntegrator adaptive
    082     * step size integrators}).</p>
    083     *
    084     * @see StepHandler
    085     * @see StepInterpolator
    086     * @version $Id: ContinuousOutputModel.java 1416643 2012-12-03 19:37:14Z tn $
    087     * @since 1.2
    088     */
    089
    090    public class ContinuousOutputModel
    091      implements StepHandler, Serializable {
    092
    093        /** Serializable version identifier */
    094        private static final long serialVersionUID = -1417964919405031606L;
    095
    096        /** Initial integration time. */
    097        private double initialTime;
    098
    099        /** Final integration time. */
    100        private double finalTime;
    101
    102        /** Integration direction indicator. */
    103        private boolean forward;
    104
    105        /** Current interpolator index. */
    106        private int index;
    107
    108        /** Steps table. */
    109        private List<StepInterpolator> steps;
    110
    111      /** Simple constructor.
    112       * Build an empty continuous output model.
    113       */
    114      public ContinuousOutputModel() {
    115        steps = new ArrayList<StepInterpolator>();
    116        initialTime = Double.NaN;
    117        finalTime   = Double.NaN;
    118        forward     = true;
    119        index       = 0;
    120      }
    121
    122      /** Append another model at the end of the instance.
    123       * @param model model to add at the end of the instance
    124       * @exception MathIllegalArgumentException if the model to append is not
    125       * compatible with the instance (dimension of the state vector,
    126       * propagation direction, hole between the dates)
    127       * @exception MaxCountExceededException if the number of functions evaluations is exceeded
    128       * during step finalization
    129       */
    130      public void append(final ContinuousOutputModel model)
    131        throws MathIllegalArgumentException, MaxCountExceededException {
    132
    133        if (model.steps.size() == 0) {
    134          return;
    135        }
    136
    137        if (steps.size() == 0) {
    138          initialTime = model.initialTime;
    139          forward     = model.forward;
    140        } else {
    141
    142          if (getInterpolatedState().length != model.getInterpolatedState().length) {
    143              throw new DimensionMismatchException(model.getInterpolatedState().length,
    144                                                   getInterpolatedState().length);
    145          }
    146
    147          if (forward ^ model.forward) {
    148              throw new MathIllegalArgumentException(LocalizedFormats.PROPAGATION_DIRECTION_MISMATCH);
    149          }
    150
    151          final StepInterpolator lastInterpolator = steps.get(index);
    152          final double current  = lastInterpolator.getCurrentTime();
    153          final double previous = lastInterpolator.getPreviousTime();
    154          final double step = current - previous;
    155          final double gap = model.getInitialTime() - current;
    156          if (FastMath.abs(gap) > 1.0e-3 * FastMath.abs(step)) {
    157            throw new MathIllegalArgumentException(LocalizedFormats.HOLE_BETWEEN_MODELS_TIME_RANGES,
    158                                                   FastMath.abs(gap));
    159          }
    160
    161        }
    162
    163        for (StepInterpolator interpolator : model.steps) {
    164          steps.add(interpolator.copy());
    165        }
    166
    167        index = steps.size() - 1;
    168        finalTime = (steps.get(index)).getCurrentTime();
    169
    170      }
    171
    172      /** {@inheritDoc} */
    173      public void init(double t0, double[] y0, double t) {
    174        initialTime = Double.NaN;
    175        finalTime   = Double.NaN;
    176        forward     = true;
    177        index       = 0;
    178        steps.clear();
    179      }
    180
    181      /** Handle the last accepted step.
    182       * A copy of the information provided by the last step is stored in
    183       * the instance for later use.
    184       * @param interpolator interpolator for the last accepted step.
    185       * @param isLast true if the step is the last one
    186       * @exception MaxCountExceededException if the number of functions evaluations is exceeded
    187       * during step finalization
    188       */
    189      public void handleStep(final StepInterpolator interpolator, final boolean isLast)
    190          throws MaxCountExceededException {
    191
    192        if (steps.size() == 0) {
    193          initialTime = interpolator.getPreviousTime();
    194          forward     = interpolator.isForward();
    195        }
    196
    197        steps.add(interpolator.copy());
    198
    199        if (isLast) {
    200          finalTime = interpolator.getCurrentTime();
    201          index     = steps.size() - 1;
    202        }
    203
    204      }
    205
    206      /**
    207       * Get the initial integration time.
    208       * @return initial integration time
    209       */
    210      public double getInitialTime() {
    211        return initialTime;
    212      }
    213
    214      /**
    215       * Get the final integration time.
    216       * @return final integration time
    217       */
    218      public double getFinalTime() {
    219        return finalTime;
    220      }
    221
    222      /**
    223       * Get the time of the interpolated point.
    224       * If {@link #setInterpolatedTime} has not been called, it returns
    225       * the final integration time.
    226       * @return interpolation point time
    227       */
    228      public double getInterpolatedTime() {
    229        return steps.get(index).getInterpolatedTime();
    230      }
    231
    232      /** Set the time of the interpolated point.
    233       * <p>This method should <strong>not</strong> be called before the
    234       * integration is over because some internal variables are set only
    235       * once the last step has been handled.</p>
    236       * <p>Setting the time outside of the integration interval is now
    237       * allowed (it was not allowed up to version 5.9 of Mantissa), but
    238       * should be used with care since the accuracy of the interpolator
    239       * will probably be very poor far from this interval. This allowance
    240       * has been added to simplify implementation of search algorithms
    241       * near the interval endpoints.</p>
    242       * @param time time of the interpolated point
    243       */
    244      public void setInterpolatedTime(final double time) {
    245
    246          // initialize the search with the complete steps table
    247          int iMin = 0;
    248          final StepInterpolator sMin = steps.get(iMin);
    249          double tMin = 0.5 * (sMin.getPreviousTime() + sMin.getCurrentTime());
    250
    251          int iMax = steps.size() - 1;
    252          final StepInterpolator sMax = steps.get(iMax);
    253          double tMax = 0.5 * (sMax.getPreviousTime() + sMax.getCurrentTime());
    254
    255          // handle points outside of the integration interval
    256          // or in the first and last step
    257          if (locatePoint(time, sMin) <= 0) {
    258            index = iMin;
    259            sMin.setInterpolatedTime(time);
    260            return;
    261          }
    262          if (locatePoint(time, sMax) >= 0) {
    263            index = iMax;
    264            sMax.setInterpolatedTime(time);
    265            return;
    266          }
    267
    268          // reduction of the table slice size
    269          while (iMax - iMin > 5) {
    270
    271            // use the last estimated index as the splitting index
    272            final StepInterpolator si = steps.get(index);
    273            final int location = locatePoint(time, si);
    274            if (location < 0) {
    275              iMax = index;
    276              tMax = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
    277            } else if (location > 0) {
    278              iMin = index;
    279              tMin = 0.5 * (si.getPreviousTime() + si.getCurrentTime());
    280            } else {
    281              // we have found the target step, no need to continue searching
    282              si.setInterpolatedTime(time);
    283              return;
    284            }
    285
    286            // compute a new estimate of the index in the reduced table slice
    287            final int iMed = (iMin + iMax) / 2;
    288            final StepInterpolator sMed = steps.get(iMed);
    289            final double tMed = 0.5 * (sMed.getPreviousTime() + sMed.getCurrentTime());
    290
    291            if ((FastMath.abs(tMed - tMin) < 1e-6) || (FastMath.abs(tMax - tMed) < 1e-6)) {
    292              // too close to the bounds, we estimate using a simple dichotomy
    293              index = iMed;
    294            } else {
    295              // estimate the index using a reverse quadratic polynom
    296              // (reverse means we have i = P(t), thus allowing to simply
    297              // compute index = P(time) rather than solving a quadratic equation)
    298              final double d12 = tMax - tMed;
    299              final double d23 = tMed - tMin;
    300              final double d13 = tMax - tMin;
    301              final double dt1 = time - tMax;
    302              final double dt2 = time - tMed;
    303              final double dt3 = time - tMin;
    304              final double iLagrange = ((dt2 * dt3 * d23) * iMax -
    305                                        (dt1 * dt3 * d13) * iMed +
    306                                        (dt1 * dt2 * d12) * iMin) /
    307                                       (d12 * d23 * d13);
    308              index = (int) FastMath.rint(iLagrange);
    309            }
    310
    311            // force the next size reduction to be at least one tenth
    312            final int low  = FastMath.max(iMin + 1, (9 * iMin + iMax) / 10);
    313            final int high = FastMath.min(iMax - 1, (iMin + 9 * iMax) / 10);
    314            if (index < low) {
    315              index = low;
    316            } else if (index > high) {
    317              index = high;
    318            }
    319
    320          }
    321
    322          // now the table slice is very small, we perform an iterative search
    323          index = iMin;
    324          while ((index <= iMax) && (locatePoint(time, steps.get(index)) > 0)) {
    325            ++index;
    326          }
    327
    328          steps.get(index).setInterpolatedTime(time);
    329
    330      }
    331
    332      /**
    333       * Get the state vector of the interpolated point.
    334       * @return state vector at time {@link #getInterpolatedTime}
    335       * @exception MaxCountExceededException if the number of functions evaluations is exceeded
    336       */
    337      public double[] getInterpolatedState() throws MaxCountExceededException {
    338        return steps.get(index).getInterpolatedState();
    339      }
    340
    341      /** Compare a step interval and a double.
    342       * @param time point to locate
    343       * @param interval step interval
    344       * @return -1 if the double is before the interval, 0 if it is in
    345       * the interval, and +1 if it is after the interval, according to
    346       * the interval direction
    347       */
    348      private int locatePoint(final double time, final StepInterpolator interval) {
    349        if (forward) {
    350          if (time < interval.getPreviousTime()) {
    351            return -1;
    352          } else if (time > interval.getCurrentTime()) {
    353            return +1;
    354          } else {
    355            return 0;
    356          }
    357        }
    358        if (time > interval.getPreviousTime()) {
    359          return -1;
    360        } else if (time < interval.getCurrentTime()) {
    361          return +1;
    362        } else {
    363          return 0;
    364        }
    365      }
    366
    367    }