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
018package org.apache.commons.math3.ode.sampling;
019
020import org.apache.commons.math3.exception.MaxCountExceededException;
021import org.apache.commons.math3.util.FastMath;
022import 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
092public 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}