1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.math4.legacy.ode.sampling;
19
20 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
21 import org.apache.commons.math4.core.jdkmath.JdkMath;
22 import org.apache.commons.numbers.core.Precision;
23
24 /**
25 * This class wraps an object implementing {@link FixedStepHandler}
26 * into a {@link StepHandler}.
27
28 * <p>This wrapper allows to use fixed step handlers with general
29 * integrators which cannot guaranty their integration steps will
30 * remain constant and therefore only accept general step
31 * handlers.</p>
32 *
33 * <p>The stepsize used is selected at construction time. The {@link
34 * FixedStepHandler#handleStep handleStep} method of the underlying
35 * {@link FixedStepHandler} object is called at normalized times. The
36 * normalized times can be influenced by the {@link StepNormalizerMode} and
37 * {@link StepNormalizerBounds}.</p>
38 *
39 * <p>There is no constraint on the integrator, it can use any time step
40 * it needs (time steps longer or shorter than the fixed time step and
41 * non-integer ratios are all allowed).</p>
42 *
43 * <table border="1" style="text-align: center">
44 * <caption>Examples (step size = 0.5)</caption>
45 * <tr style="background-color: #CCCCFF"><td style="font-size: x-large">Examples (step size = 0.5)</td></tr>
46 * <tr style="background-color: #EEEEFF; font-size: large"><td>Start time</td><td>End time</td>
47 * <td>Direction</td><td>{@link StepNormalizerMode Mode}</td>
48 * <td>{@link StepNormalizerBounds Bounds}</td><td>Output</td></tr>
49 * <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>
50 * <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>
51 * <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>
52 * <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>
53 * <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>
54 * <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>
55 * <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>
56 * <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>
57 * <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>
58 * <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>
59 * <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>
60 * <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>
61 * <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>
62 * <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>
63 * <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>
64 * <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>
65 * <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>
66 * <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>
67 * <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>
68 * <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>
69 * <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>
70 * <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>
71 * <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>
72 * <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>
73 * <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>
74 * <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>
75 * <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>
76 * <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>
77 * <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>
78 * <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>
79 * <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>
80 * <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>
81 * </table>
82 *
83 * @see StepHandler
84 * @see FixedStepHandler
85 * @see StepNormalizerMode
86 * @see StepNormalizerBounds
87 * @since 1.2
88 */
89
90 public class StepNormalizer implements StepHandler {
91 /** Fixed time step. */
92 private double h;
93
94 /** Underlying step handler. */
95 private final FixedStepHandler handler;
96
97 /** First step time. */
98 private double firstTime;
99
100 /** Last step time. */
101 private double lastTime;
102
103 /** Last state vector. */
104 private double[] lastState;
105
106 /** Last derivatives vector. */
107 private double[] lastDerivatives;
108
109 /** Integration direction indicator. */
110 private boolean forward;
111
112 /** The step normalizer bounds settings to use. */
113 private final StepNormalizerBounds bounds;
114
115 /** The step normalizer mode to use. */
116 private final StepNormalizerMode mode;
117
118 /** Simple constructor. Uses {@link StepNormalizerMode#INCREMENT INCREMENT}
119 * mode, and {@link StepNormalizerBounds#FIRST FIRST} bounds setting, for
120 * backwards compatibility.
121 * @param h fixed time step (sign is not used)
122 * @param handler fixed time step handler to wrap
123 */
124 public StepNormalizer(final double h, final FixedStepHandler handler) {
125 this(h, handler, StepNormalizerMode.INCREMENT,
126 StepNormalizerBounds.FIRST);
127 }
128
129 /** Simple constructor. Uses {@link StepNormalizerBounds#FIRST FIRST}
130 * bounds setting.
131 * @param h fixed time step (sign is not used)
132 * @param handler fixed time step handler to wrap
133 * @param mode step normalizer mode to use
134 * @since 3.0
135 */
136 public StepNormalizer(final double h, final FixedStepHandler handler,
137 final StepNormalizerMode mode) {
138 this(h, handler, mode, StepNormalizerBounds.FIRST);
139 }
140
141 /** Simple constructor. Uses {@link StepNormalizerMode#INCREMENT INCREMENT}
142 * mode.
143 * @param h fixed time step (sign is not used)
144 * @param handler fixed time step handler to wrap
145 * @param bounds step normalizer bounds setting to use
146 * @since 3.0
147 */
148 public StepNormalizer(final double h, final FixedStepHandler handler,
149 final StepNormalizerBounds bounds) {
150 this(h, handler, StepNormalizerMode.INCREMENT, bounds);
151 }
152
153 /** Simple constructor.
154 * @param h fixed time step (sign is not used)
155 * @param handler fixed time step handler to wrap
156 * @param mode step normalizer mode to use
157 * @param bounds step normalizer bounds setting to use
158 * @since 3.0
159 */
160 public StepNormalizer(final double h, final FixedStepHandler handler,
161 final StepNormalizerMode mode,
162 final StepNormalizerBounds bounds) {
163 this.h = JdkMath.abs(h);
164 this.handler = handler;
165 this.mode = mode;
166 this.bounds = bounds;
167 firstTime = Double.NaN;
168 lastTime = Double.NaN;
169 lastState = null;
170 lastDerivatives = null;
171 forward = true;
172 }
173
174 /** {@inheritDoc} */
175 @Override
176 public void init(double t0, double[] y0, double t) {
177
178 firstTime = Double.NaN;
179 lastTime = Double.NaN;
180 lastState = null;
181 lastDerivatives = null;
182 forward = true;
183
184 // initialize the underlying handler
185 handler.init(t0, y0, t);
186 }
187
188 /**
189 * Handle the last accepted step.
190 * @param interpolator interpolator for the last accepted step. For
191 * efficiency purposes, the various integrators reuse the same
192 * object on each call, so if the instance wants to keep it across
193 * all calls (for example to provide at the end of the integration a
194 * continuous model valid throughout the integration range), it
195 * should build a local copy using the clone method and store this
196 * copy.
197 * @param isLast true if the step is the last one
198 * @exception MaxCountExceededException if the interpolator throws one because
199 * the number of functions evaluations is exceeded
200 */
201 @Override
202 public void handleStep(final StepInterpolator interpolator, final boolean isLast)
203 throws MaxCountExceededException {
204 // The first time, update the last state with the start information.
205 if (lastState == null) {
206 firstTime = interpolator.getPreviousTime();
207 lastTime = interpolator.getPreviousTime();
208 interpolator.setInterpolatedTime(lastTime);
209 lastState = interpolator.getInterpolatedState().clone();
210 lastDerivatives = interpolator.getInterpolatedDerivatives().clone();
211
212 // Take the integration direction into account.
213 forward = interpolator.getCurrentTime() >= lastTime;
214 if (!forward) {
215 h = -h;
216 }
217 }
218
219 // Calculate next normalized step time.
220 double nextTime = (mode == StepNormalizerMode.INCREMENT) ?
221 lastTime + h :
222 (JdkMath.floor(lastTime / h) + 1) * h;
223 if (mode == StepNormalizerMode.MULTIPLES &&
224 Precision.equals(nextTime, lastTime, 1)) {
225 nextTime += h;
226 }
227
228 // Process normalized steps as long as they are in the current step.
229 boolean nextInStep = isNextInStep(nextTime, interpolator);
230 while (nextInStep) {
231 // Output the stored previous step.
232 doNormalizedStep(false);
233
234 // Store the next step as last step.
235 storeStep(interpolator, nextTime);
236
237 // Move on to the next step.
238 nextTime += h;
239 nextInStep = isNextInStep(nextTime, interpolator);
240 }
241
242 if (isLast) {
243 // There will be no more steps. The stored one should be given to
244 // the handler. We may have to output one more step. Only the last
245 // one of those should be flagged as being the last.
246 boolean addLast = bounds.lastIncluded() &&
247 lastTime != interpolator.getCurrentTime();
248 doNormalizedStep(!addLast);
249 if (addLast) {
250 storeStep(interpolator, interpolator.getCurrentTime());
251 doNormalizedStep(true);
252 }
253 }
254 }
255
256 /**
257 * Returns a value indicating whether the next normalized time is in the
258 * current step.
259 * @param nextTime the next normalized time
260 * @param interpolator interpolator for the last accepted step, to use to
261 * get the end time of the current step
262 * @return value indicating whether the next normalized time is in the
263 * current step
264 */
265 private boolean isNextInStep(double nextTime,
266 StepInterpolator interpolator) {
267 return forward ?
268 nextTime <= interpolator.getCurrentTime() :
269 nextTime >= interpolator.getCurrentTime();
270 }
271
272 /**
273 * Invokes the underlying step handler for the current normalized step.
274 * @param isLast true if the step is the last one
275 */
276 private void doNormalizedStep(boolean isLast) {
277 if (!bounds.firstIncluded() && firstTime == lastTime) {
278 return;
279 }
280 handler.handleStep(lastTime, lastState, lastDerivatives, isLast);
281 }
282
283 /** Stores the interpolated information for the given time in the current
284 * state.
285 * @param interpolator interpolator for the last accepted step, to use to
286 * get the interpolated information
287 * @param t the time for which to store the interpolated information
288 * @exception MaxCountExceededException if the interpolator throws one because
289 * the number of functions evaluations is exceeded
290 */
291 private void storeStep(StepInterpolator interpolator, double t)
292 throws MaxCountExceededException {
293 lastTime = t;
294 interpolator.setInterpolatedTime(lastTime);
295 System.arraycopy(interpolator.getInterpolatedState(), 0,
296 lastState, 0, lastState.length);
297 System.arraycopy(interpolator.getInterpolatedDerivatives(), 0,
298 lastDerivatives, 0, lastDerivatives.length);
299 }
300 }