    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.sampling;
    019
    020    import org.apache.commons.math3.exception.MaxCountExceededException;
    021    import org.apache.commons.math3.util.FastMath;
    022    import org.apache.commons.math3.util.Precision;
    023
    024    /**
    025     * This class wraps an object implementing {@link FixedStepHandler}
    026     * into a {@link StepHandler}.
    027
    028     * <p>This wrapper allows to use fixed step handlers with general
    029     * integrators which cannot guaranty their integration steps will
    030     * remain constant and therefore only accept general step
    031     * handlers.</p>
    032     *
    033     * <p>The stepsize used is selected at construction time. The {@link
    034     * FixedStepHandler#handleStep handleStep} method of the underlying
    035     * {@link FixedStepHandler} object is called at normalized times. The
    036     * normalized times can be influenced by the {@link StepNormalizerMode} and
    037     * {@link StepNormalizerBounds}.</p>
    038     *
    039     * <p>There is no constraint on the integrator, it can use any time step
    040     * it needs (time steps longer or shorter than the fixed time step and
    041     * non-integer ratios are all allowed).</p>
    042     *
    043     * <p>
    044     * <table border="1" align="center">
    045     * <tr BGCOLOR="#CCCCFF"><td colspan=6><font size="+2">Examples (step size = 0.5)</font></td></tr>
    046     * <tr BGCOLOR="#EEEEFF"><font size="+1"><td>Start time</td><td>End time</td>
    047     *  <td>Direction</td><td>{@link StepNormalizerMode Mode}</td>
    048     *  <td>{@link StepNormalizerBounds Bounds}</td><td>Output</td></font></tr>
    049     * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>0.8, 1.3, 1.8, 2.3, 2.8</td></tr>
    050     * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>0.3, 0.8, 1.3, 1.8, 2.3, 2.8</td></tr>
    051     * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>0.8, 1.3, 1.8, 2.3, 2.8, 3.1</td></tr>
    052     * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>0.3, 0.8, 1.3, 1.8, 2.3, 2.8, 3.1</td></tr>
    053     * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr>
    054     * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>0.3, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr>
    055     * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.1</td></tr>
    056     * <tr><td>0.3</td><td>3.1</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>0.3, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.1</td></tr>
    057     * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr>
    058     * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr>
    059     * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr>
    060     * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr>
    061     * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr>
    062     * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr>
    063     * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr>
    064     * <tr><td>0.0</td><td>3.0</td><td>forward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0</td></tr>
    065     * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>2.6, 2.1, 1.6, 1.1, 0.6</td></tr>
    066     * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>3.1, 2.6, 2.1, 1.6, 1.1, 0.6</td></tr>
    067     * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>2.6, 2.1, 1.6, 1.1, 0.6, 0.3</td></tr>
    068     * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>3.1, 2.6, 2.1, 1.6, 1.1, 0.6, 0.3</td></tr>
    069     * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5</td></tr>
    070     * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>3.1, 3.0, 2.5, 2.0, 1.5, 1.0, 0.5</td></tr>
    071     * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.3</td></tr>
    072     * <tr><td>3.1</td><td>0.3</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>3.1, 3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.3</td></tr>
    073     * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr>
    074     * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr>
    075     * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr>
    076     * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#INCREMENT INCREMENT}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr>
    077     * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#NEITHER NEITHER}</td><td>2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr>
    078     * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#FIRST FIRST}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr>
    079     * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#LAST LAST}</td><td>2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr>
    080     * <tr><td>3.0</td><td>0.0</td><td>backward</td><td>{@link StepNormalizerMode#MULTIPLES MULTIPLES}</td><td>{@link StepNormalizerBounds#BOTH BOTH}</td><td>3.0, 2.5, 2.0, 1.5, 1.0, 0.5, 0.0</td></tr>
    081     * </table>
    082     * </p>
    083     *
    084     * @see StepHandler
    085     * @see FixedStepHandler
    086     * @see StepNormalizerMode
    087     * @see StepNormalizerBounds
    088     * @version $Id: StepNormalizer.java 1416643 2012-12-03 19:37:14Z tn $
    089     * @since 1.2
    090     */
    091
    092    public class StepNormalizer implements StepHandler {
    093        /** Fixed time step. */
    094        private double h;
    095
    096        /** Underlying step handler. */
    097        private final FixedStepHandler handler;
    098
    099        /** First step time. */
    100        private double firstTime;
    101
    102        /** Last step time. */
    103        private double lastTime;
    104
    105        /** Last state vector. */
    106        private double[] lastState;
    107
    108        /** Last derivatives vector. */
    109        private double[] lastDerivatives;
    110
    111        /** Integration direction indicator. */
    112        private boolean forward;
    113
    114        /** The step normalizer bounds settings to use. */
    115        private final StepNormalizerBounds bounds;
    116
    117        /** The step normalizer mode to use. */
    118        private final StepNormalizerMode mode;
    119
    120        /** Simple constructor. Uses {@link StepNormalizerMode#INCREMENT INCREMENT}
    121         * mode, and {@link StepNormalizerBounds#FIRST FIRST} bounds setting, for
    122         * backwards compatibility.
    123         * @param h fixed time step (sign is not used)
    124         * @param handler fixed time step handler to wrap
    125         */
    126        public StepNormalizer(final double h, final FixedStepHandler handler) {
    127            this(h, handler, StepNormalizerMode.INCREMENT,
    128                 StepNormalizerBounds.FIRST);
    129        }
    130
    131        /** Simple constructor. Uses {@link StepNormalizerBounds#FIRST FIRST}
    132         * bounds setting.
    133         * @param h fixed time step (sign is not used)
    134         * @param handler fixed time step handler to wrap
    135         * @param mode step normalizer mode to use
    136         * @since 3.0
    137         */
    138        public StepNormalizer(final double h, final FixedStepHandler handler,
    139                              final StepNormalizerMode mode) {
    140            this(h, handler, mode, StepNormalizerBounds.FIRST);
    141        }
    142
    143        /** Simple constructor. Uses {@link StepNormalizerMode#INCREMENT INCREMENT}
    144         * mode.
    145         * @param h fixed time step (sign is not used)
    146         * @param handler fixed time step handler to wrap
    147         * @param bounds step normalizer bounds setting to use
    148         * @since 3.0
    149         */
    150        public StepNormalizer(final double h, final FixedStepHandler handler,
    151                              final StepNormalizerBounds bounds) {
    152            this(h, handler, StepNormalizerMode.INCREMENT, bounds);
    153        }
    154
    155        /** Simple constructor.
    156         * @param h fixed time step (sign is not used)
    157         * @param handler fixed time step handler to wrap
    158         * @param mode step normalizer mode to use
    159         * @param bounds step normalizer bounds setting to use
    160         * @since 3.0
    161         */
    162        public StepNormalizer(final double h, final FixedStepHandler handler,
    163                              final StepNormalizerMode mode,
    164                              final StepNormalizerBounds bounds) {
    165            this.h          = FastMath.abs(h);
    166            this.handler    = handler;
    167            this.mode       = mode;
    168            this.bounds     = bounds;
    169            firstTime       = Double.NaN;
    170            lastTime        = Double.NaN;
    171            lastState       = null;
    172            lastDerivatives = null;
    173            forward         = true;
    174        }
    175
    176        /** {@inheritDoc} */
    177        public void init(double t0, double[] y0, double t) {
    178
    179            firstTime       = Double.NaN;
    180            lastTime        = Double.NaN;
    181            lastState       = null;
    182            lastDerivatives = null;
    183            forward         = true;
    184
    185            // initialize the underlying handler
    186            handler.init(t0, y0, t);
    187
    188        }
    189
    190        /**
    191         * Handle the last accepted step
    192         * @param interpolator interpolator for the last accepted step. For
    193         * efficiency purposes, the various integrators reuse the same
    194         * object on each call, so if the instance wants to keep it across
    195         * all calls (for example to provide at the end of the integration a
    196         * continuous model valid throughout the integration range), it
    197         * should build a local copy using the clone method and store this
    198         * copy.
    199         * @param isLast true if the step is the last one
    200         * @exception MaxCountExceededException if the interpolator throws one because
    201         * the number of functions evaluations is exceeded
    202         */
    203        public void handleStep(final StepInterpolator interpolator, final boolean isLast)
    204            throws MaxCountExceededException {
    205            // The first time, update the last state with the start information.
    206            if (lastState == null) {
    207                firstTime = interpolator.getPreviousTime();
    208                lastTime = interpolator.getPreviousTime();
    209                interpolator.setInterpolatedTime(lastTime);
    210                lastState = interpolator.getInterpolatedState().clone();
    211                lastDerivatives = interpolator.getInterpolatedDerivatives().clone();
    212
    213                // Take the integration direction into account.
    214                forward = interpolator.getCurrentTime() >= lastTime;
    215                if (!forward) {
    216                    h = -h;
    217                }
    218            }
    219
    220            // Calculate next normalized step time.
    221            double nextTime = (mode == StepNormalizerMode.INCREMENT) ?
    222                              lastTime + h :
    223                              (FastMath.floor(lastTime / h) + 1) * h;
    224            if (mode == StepNormalizerMode.MULTIPLES &&
    225                Precision.equals(nextTime, lastTime, 1)) {
    226                nextTime += h;
    227            }
    228
    229            // Process normalized steps as long as they are in the current step.
    230            boolean nextInStep = isNextInStep(nextTime, interpolator);
    231            while (nextInStep) {
    232                // Output the stored previous step.
    233                doNormalizedStep(false);
    234
    235                // Store the next step as last step.
    236                storeStep(interpolator, nextTime);
    237
    238                // Move on to the next step.
    239                nextTime += h;
    240                nextInStep = isNextInStep(nextTime, interpolator);
    241            }
    242
    243            if (isLast) {
    244                // There will be no more steps. The stored one should be given to
    245                // the handler. We may have to output one more step. Only the last
    246                // one of those should be flagged as being the last.
    247                boolean addLast = bounds.lastIncluded() &&
    248                                  lastTime != interpolator.getCurrentTime();
    249                doNormalizedStep(!addLast);
    250                if (addLast) {
    251                    storeStep(interpolator, interpolator.getCurrentTime());
    252                    doNormalizedStep(true);
    253                }
    254            }
    255        }
    256
    257        /**
    258         * Returns a value indicating whether the next normalized time is in the
    259         * current step.
    260         * @param nextTime the next normalized time
    261         * @param interpolator interpolator for the last accepted step, to use to
    262         * get the end time of the current step
    263         * @return value indicating whether the next normalized time is in the
    264         * current step
    265         */
    266        private boolean isNextInStep(double nextTime,
    267                                     StepInterpolator interpolator) {
    268            return forward ?
    269                   nextTime <= interpolator.getCurrentTime() :
    270                   nextTime >= interpolator.getCurrentTime();
    271        }
    272
    273        /**
    274         * Invokes the underlying step handler for the current normalized step.
    275         * @param isLast true if the step is the last one
    276         */
    277        private void doNormalizedStep(boolean isLast) {
    278            if (!bounds.firstIncluded() && firstTime == lastTime) {
    279                return;
    280            }
    281            handler.handleStep(lastTime, lastState, lastDerivatives, isLast);
    282        }
    283
    284        /** Stores the interpolated information for the given time in the current
    285         * state.
    286         * @param interpolator interpolator for the last accepted step, to use to
    287         * get the interpolated information
    288         * @param t the time for which to store the interpolated information
    289         * @exception MaxCountExceededException if the interpolator throws one because
    290         * the number of functions evaluations is exceeded
    291         */
    292        private void storeStep(StepInterpolator interpolator, double t)
    293            throws MaxCountExceededException {
    294            lastTime = t;
    295            interpolator.setInterpolatedTime(lastTime);
    296            System.arraycopy(interpolator.getInterpolatedState(), 0,
    297                             lastState, 0, lastState.length);
    298            System.arraycopy(interpolator.getInterpolatedDerivatives(), 0,
    299                             lastDerivatives, 0, lastDerivatives.length);
    300        }
    301    }