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    }