View Javadoc
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.core.RealFieldElement;
21  import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
22  import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
23  import org.apache.commons.math4.core.jdkmath.JdkMath;
24  import org.apache.commons.numbers.core.Precision;
25  
26  /**
27   * This class wraps an object implementing {@link FieldFixedStepHandler}
28   * into a {@link FieldStepHandler}.
29  
30   * <p>This wrapper allows to use fixed step handlers with general
31   * integrators which cannot guaranty their integration steps will
32   * remain constant and therefore only accept general step
33   * handlers.</p>
34   *
35   * <p>The stepsize used is selected at construction time. The {@link
36   * FieldFixedStepHandler#handleStep handleStep} method of the underlying
37   * {@link FieldFixedStepHandler} object is called at normalized times. The
38   * normalized times can be influenced by the {@link StepNormalizerMode} and
39   * {@link StepNormalizerBounds}.</p>
40   *
41   * <p>There is no constraint on the integrator, it can use any time step
42   * it needs (time steps longer or shorter than the fixed time step and
43   * non-integer ratios are all allowed).</p>
44   *
45   * <table border="1" style="text-align: center">
46   * <caption>Examples (step size = 0.5)</caption>
47   * <tr style="background-color: #CCCCFF"><td colspan=6 style="font-size: x-large">Examples (step size = 0.5)</td></tr>
48   * <tr style="background-color: #EEEEFF; font-size: large"><td>Start time</td><td>End time</td>
49   *  <td>Direction</td><td>{@link StepNormalizerMode Mode}</td>
50   *  <td>{@link StepNormalizerBounds Bounds}</td><td>Output</td></tr>
51   * <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>
52   * <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>
53   * <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>
54   * <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>
55   * <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>
56   * <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>
57   * <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>
58   * <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>
59   * <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>
60   * <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>
61   * <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>
62   * <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>
63   * <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>
64   * <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>
65   * <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>
66   * <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>
67   * <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>
68   * <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>
69   * <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>
70   * <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>
71   * <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>
72   * <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>
73   * <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>
74   * <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>
75   * <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>
76   * <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>
77   * <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>
78   * <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>
79   * <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>
80   * <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>
81   * <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>
82   * <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>
83   * </table>
84   *
85   * @param <T> the type of the field elements
86   * @see FieldStepHandler
87   * @see FieldFixedStepHandler
88   * @see StepNormalizerMode
89   * @see StepNormalizerBounds
90   * @since 3.6
91   */
92  
93  public class FieldStepNormalizer<T extends RealFieldElement<T>> implements FieldStepHandler<T> {
94  
95      /** Fixed time step. */
96      private double h;
97  
98      /** Underlying step handler. */
99      private final FieldFixedStepHandler<T> handler;
100 
101     /** First step state. */
102     private FieldODEStateAndDerivative<T> first;
103 
104     /** Last step step. */
105     private FieldODEStateAndDerivative<T> last;
106 
107     /** Integration direction indicator. */
108     private boolean forward;
109 
110     /** The step normalizer bounds settings to use. */
111     private final StepNormalizerBounds bounds;
112 
113     /** The step normalizer mode to use. */
114     private final StepNormalizerMode mode;
115 
116     /** Simple constructor. Uses {@link StepNormalizerMode#INCREMENT INCREMENT}
117      * mode, and {@link StepNormalizerBounds#FIRST FIRST} bounds setting, for
118      * backwards compatibility.
119      * @param h fixed time step (sign is not used)
120      * @param handler fixed time step handler to wrap
121      */
122     public FieldStepNormalizer(final double h, final FieldFixedStepHandler<T> handler) {
123         this(h, handler, StepNormalizerMode.INCREMENT,
124              StepNormalizerBounds.FIRST);
125     }
126 
127     /** Simple constructor. Uses {@link StepNormalizerBounds#FIRST FIRST}
128      * bounds setting.
129      * @param h fixed time step (sign is not used)
130      * @param handler fixed time step handler to wrap
131      * @param mode step normalizer mode to use
132      * @since 3.0
133      */
134     public FieldStepNormalizer(final double h, final FieldFixedStepHandler<T> handler,
135                                final StepNormalizerMode mode) {
136         this(h, handler, mode, StepNormalizerBounds.FIRST);
137     }
138 
139     /** Simple constructor. Uses {@link StepNormalizerMode#INCREMENT INCREMENT}
140      * mode.
141      * @param h fixed time step (sign is not used)
142      * @param handler fixed time step handler to wrap
143      * @param bounds step normalizer bounds setting to use
144      * @since 3.0
145      */
146     public FieldStepNormalizer(final double h, final FieldFixedStepHandler<T> handler,
147                                final StepNormalizerBounds bounds) {
148         this(h, handler, StepNormalizerMode.INCREMENT, bounds);
149     }
150 
151     /** Simple constructor.
152      * @param h fixed time step (sign is not used)
153      * @param handler fixed time step handler to wrap
154      * @param mode step normalizer mode to use
155      * @param bounds step normalizer bounds setting to use
156      * @since 3.0
157      */
158     public FieldStepNormalizer(final double h, final FieldFixedStepHandler<T> handler,
159                                final StepNormalizerMode mode, final StepNormalizerBounds bounds) {
160         this.h       = JdkMath.abs(h);
161         this.handler = handler;
162         this.mode    = mode;
163         this.bounds  = bounds;
164         first        = null;
165         last         = null;
166         forward      = true;
167     }
168 
169     /** {@inheritDoc} */
170     @Override
171     public void init(final FieldODEStateAndDerivative<T> initialState, final T finalTime) {
172 
173         first   = null;
174         last    = null;
175         forward = true;
176 
177         // initialize the underlying handler
178         handler.init(initialState, finalTime);
179     }
180 
181     /**
182      * Handle the last accepted step.
183      * @param interpolator interpolator for the last accepted step. For
184      * efficiency purposes, the various integrators reuse the same
185      * object on each call, so if the instance wants to keep it across
186      * all calls (for example to provide at the end of the integration a
187      * continuous model valid throughout the integration range), it
188      * should build a local copy using the clone method and store this
189      * copy.
190      * @param isLast true if the step is the last one
191      * @exception MaxCountExceededException if the interpolator throws one because
192      * the number of functions evaluations is exceeded
193      */
194     @Override
195     public void handleStep(final FieldStepInterpolator<T> interpolator, final boolean isLast)
196         throws MaxCountExceededException {
197         // The first time, update the last state with the start information.
198         if (last == null) {
199 
200             first   = interpolator.getPreviousState();
201             last    = first;
202 
203             // Take the integration direction into account.
204             forward = interpolator.isForward();
205             if (!forward) {
206                 h = -h;
207             }
208         }
209 
210         // Calculate next normalized step time.
211         T nextTime = (mode == StepNormalizerMode.INCREMENT) ?
212                      last.getTime().add(h) :
213                      last.getTime().getField().getZero().add((JdkMath.floor(last.getTime().getReal() / h) + 1) * h);
214         if (mode == StepNormalizerMode.MULTIPLES &&
215             Precision.equals(nextTime.getReal(), last.getTime().getReal(), 1)) {
216             nextTime = nextTime.add(h);
217         }
218 
219         // Process normalized steps as long as they are in the current step.
220         boolean nextInStep = isNextInStep(nextTime, interpolator);
221         while (nextInStep) {
222             // Output the stored previous step.
223             doNormalizedStep(false);
224 
225             // Store the next step as last step.
226             last = interpolator.getInterpolatedState(nextTime);
227 
228             // Move on to the next step.
229             nextTime = nextTime.add(h);
230             nextInStep = isNextInStep(nextTime, interpolator);
231         }
232 
233         if (isLast) {
234             // There will be no more steps. The stored one should be given to
235             // the handler. We may have to output one more step. Only the last
236             // one of those should be flagged as being the last.
237             final boolean addLast = bounds.lastIncluded() &&
238                                     last.getTime().getReal() != interpolator.getCurrentState().getTime().getReal();
239             doNormalizedStep(!addLast);
240             if (addLast) {
241                 last = interpolator.getCurrentState();
242                 doNormalizedStep(true);
243             }
244         }
245     }
246 
247     /**
248      * Returns a value indicating whether the next normalized time is in the
249      * current step.
250      * @param nextTime the next normalized time
251      * @param interpolator interpolator for the last accepted step, to use to
252      * get the end time of the current step
253      * @return value indicating whether the next normalized time is in the
254      * current step
255      */
256     private boolean isNextInStep(final T nextTime, final FieldStepInterpolator<T> interpolator) {
257         return forward ?
258                nextTime.getReal() <= interpolator.getCurrentState().getTime().getReal() :
259                nextTime.getReal() >= interpolator.getCurrentState().getTime().getReal();
260     }
261 
262     /**
263      * Invokes the underlying step handler for the current normalized step.
264      * @param isLast true if the step is the last one
265      */
266     private void doNormalizedStep(final boolean isLast) {
267         if (!bounds.firstIncluded() && first.getTime().getReal() == last.getTime().getReal()) {
268             return;
269         }
270         handler.handleStep(last, isLast);
271     }
272 }