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