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 }