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.events;
19
20 import org.apache.commons.math4.legacy.core.RealFieldElement;
21 import org.apache.commons.math4.legacy.analysis.RealFieldUnivariateFunction;
22 import org.apache.commons.math4.legacy.analysis.solvers.AllowedSolution;
23 import org.apache.commons.math4.legacy.analysis.solvers.BracketedRealFieldUnivariateSolver;
24 import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
25 import org.apache.commons.math4.legacy.exception.NoBracketingException;
26 import org.apache.commons.math4.legacy.ode.FieldODEState;
27 import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
28 import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator;
29 import org.apache.commons.math4.core.jdkmath.JdkMath;
30
31 /** This class handles the state for one {@link EventHandler
32 * event handler} during integration steps.
33 *
34 * <p>Each time the integrator proposes a step, the event handler
35 * switching function should be checked. This class handles the state
36 * of one handler during one integration step, with references to the
37 * state at the end of the preceding step. This information is used to
38 * decide if the handler should trigger an event or not during the
39 * proposed step.</p>
40 *
41 * @param <T> the type of the field elements
42 * @since 3.6
43 */
44 public class FieldEventState<T extends RealFieldElement<T>> {
45
46 /** Event handler. */
47 private final FieldEventHandler<T> handler;
48
49 /** Maximal time interval between events handler checks. */
50 private final double maxCheckInterval;
51
52 /** Convergence threshold for event localization. */
53 private final T convergence;
54
55 /** Upper limit in the iteration count for event localization. */
56 private final int maxIterationCount;
57
58 /** Time at the beginning of the step. */
59 private T t0;
60
61 /** Value of the events handler at the beginning of the step. */
62 private T g0;
63
64 /** Simulated sign of g0 (we cheat when crossing events). */
65 private boolean g0Positive;
66
67 /** Indicator of event expected during the step. */
68 private boolean pendingEvent;
69
70 /** Occurrence time of the pending event. */
71 private T pendingEventTime;
72
73 /** Occurrence time of the previous event. */
74 private T previousEventTime;
75
76 /** Integration direction. */
77 private boolean forward;
78
79 /** Variation direction around pending event.
80 * (this is considered with respect to the integration direction)
81 */
82 private boolean increasing;
83
84 /** Next action indicator. */
85 private Action nextAction;
86
87 /** Root-finding algorithm to use to detect state events. */
88 private final BracketedRealFieldUnivariateSolver<T> solver;
89
90 /** Simple constructor.
91 * @param handler event handler
92 * @param maxCheckInterval maximal time interval between switching
93 * function checks (this interval prevents missing sign changes in
94 * case the integration steps becomes very large)
95 * @param convergence convergence threshold in the event time search
96 * @param maxIterationCount upper limit of the iteration count in
97 * the event time search
98 * @param solver Root-finding algorithm to use to detect state events
99 */
100 public FieldEventState(final FieldEventHandler<T> handler, final double maxCheckInterval,
101 final T convergence, final int maxIterationCount,
102 final BracketedRealFieldUnivariateSolver<T> solver) {
103 this.handler = handler;
104 this.maxCheckInterval = maxCheckInterval;
105 this.convergence = convergence.abs();
106 this.maxIterationCount = maxIterationCount;
107 this.solver = solver;
108
109 // some dummy values ...
110 t0 = null;
111 g0 = null;
112 g0Positive = true;
113 pendingEvent = false;
114 pendingEventTime = null;
115 previousEventTime = null;
116 increasing = true;
117 nextAction = Action.CONTINUE;
118 }
119
120 /** Get the underlying event handler.
121 * @return underlying event handler
122 */
123 public FieldEventHandler<T> getEventHandler() {
124 return handler;
125 }
126
127 /** Get the maximal time interval between events handler checks.
128 * @return maximal time interval between events handler checks
129 */
130 public double getMaxCheckInterval() {
131 return maxCheckInterval;
132 }
133
134 /** Get the convergence threshold for event localization.
135 * @return convergence threshold for event localization
136 */
137 public T getConvergence() {
138 return convergence;
139 }
140
141 /** Get the upper limit in the iteration count for event localization.
142 * @return upper limit in the iteration count for event localization
143 */
144 public int getMaxIterationCount() {
145 return maxIterationCount;
146 }
147
148 /** Reinitialize the beginning of the step.
149 * @param interpolator valid for the current step
150 * @exception MaxCountExceededException if the interpolator throws one because
151 * the number of functions evaluations is exceeded
152 */
153 public void reinitializeBegin(final FieldStepInterpolator<T> interpolator)
154 throws MaxCountExceededException {
155
156 final FieldODEStateAndDerivative<T> s0 = interpolator.getPreviousState();
157 t0 = s0.getTime();
158 g0 = handler.g(s0);
159 if (g0.getReal() == 0) {
160 // excerpt from MATH-421 issue:
161 // If an ODE solver is setup with an EventHandler that return STOP
162 // when the even is triggered, the integrator stops (which is exactly
163 // the expected behavior). If however the user wants to restart the
164 // solver from the final state reached at the event with the same
165 // configuration (expecting the event to be triggered again at a
166 // later time), then the integrator may fail to start. It can get stuck
167 // at the previous event. The use case for the bug MATH-421 is fairly
168 // general, so events occurring exactly at start in the first step should
169 // be ignored.
170
171 // extremely rare case: there is a zero EXACTLY at interval start
172 // we will use the sign slightly after step beginning to force ignoring this zero
173 final double epsilon = JdkMath.max(solver.getAbsoluteAccuracy().getReal(),
174 JdkMath.abs(solver.getRelativeAccuracy().multiply(t0).getReal()));
175 final T tStart = t0.add(0.5 * epsilon);
176 g0 = handler.g(interpolator.getInterpolatedState(tStart));
177 }
178 g0Positive = g0.getReal() >= 0;
179 }
180
181 /** Evaluate the impact of the proposed step on the event handler.
182 * @param interpolator step interpolator for the proposed step
183 * @return true if the event handler triggers an event before
184 * the end of the proposed step
185 * @exception MaxCountExceededException if the interpolator throws one because
186 * the number of functions evaluations is exceeded
187 * @exception NoBracketingException if the event cannot be bracketed
188 */
189 public boolean evaluateStep(final FieldStepInterpolator<T> interpolator)
190 throws MaxCountExceededException, NoBracketingException {
191
192 forward = interpolator.isForward();
193 final FieldODEStateAndDerivative<T> s1 = interpolator.getCurrentState();
194 final T t1 = s1.getTime();
195 final T dt = t1.subtract(t0);
196 if (dt.abs().subtract(convergence).getReal() < 0) {
197 // we cannot do anything on such a small step, don't trigger any events
198 return false;
199 }
200 final int n = JdkMath.max(1, (int) JdkMath.ceil(JdkMath.abs(dt.getReal()) / maxCheckInterval));
201 final T h = dt.divide(n);
202
203 final RealFieldUnivariateFunction<T> f = new RealFieldUnivariateFunction<T>() {
204 /** {@inheritDoc} */
205 @Override
206 public T value(final T t) {
207 return handler.g(interpolator.getInterpolatedState(t));
208 }
209 };
210
211 T ta = t0;
212 T ga = g0;
213 for (int i = 0; i < n; ++i) {
214
215 // evaluate handler value at the end of the substep
216 final T tb = (i == n - 1) ? t1 : t0.add(h.multiply(i + 1));
217 final T gb = handler.g(interpolator.getInterpolatedState(tb));
218
219 // check events occurrence
220 if (g0Positive ^ (gb.getReal() >= 0)) {
221 // there is a sign change: an event is expected during this step
222
223 // variation direction, with respect to the integration direction
224 increasing = gb.subtract(ga).getReal() >= 0;
225
226 // find the event time making sure we select a solution just at or past the exact root
227 final T root = forward ?
228 solver.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
229 solver.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
230
231 if (previousEventTime != null &&
232 root.subtract(ta).abs().subtract(convergence).getReal() <= 0 &&
233 root.subtract(previousEventTime).abs().subtract(convergence).getReal() <= 0) {
234 // we have either found nothing or found (again ?) a past event,
235 // retry the substep excluding this value, and taking care to have the
236 // required sign in case the g function is noisy around its zero and
237 // crosses the axis several times
238 do {
239 ta = forward ? ta.add(convergence) : ta.subtract(convergence);
240 ga = f.value(ta);
241 } while ((g0Positive ^ (ga.getReal() >= 0)) && (forward ^ (ta.subtract(tb).getReal() >= 0)));
242
243 if (forward ^ (ta.subtract(tb).getReal() >= 0)) {
244 // we were able to skip this spurious root
245 --i;
246 } else {
247 // we can't avoid this root before the end of the step,
248 // we have to handle it despite it is close to the former one
249 // maybe we have two very close roots
250 pendingEventTime = root;
251 pendingEvent = true;
252 return true;
253 }
254 } else if (previousEventTime == null ||
255 previousEventTime.subtract(root).abs().subtract(convergence).getReal() > 0) {
256 pendingEventTime = root;
257 pendingEvent = true;
258 return true;
259 } else {
260 // no sign change: there is no event for now
261 ta = tb;
262 ga = gb;
263 }
264 } else {
265 // no sign change: there is no event for now
266 ta = tb;
267 ga = gb;
268 }
269 }
270
271 // no event during the whole step
272 pendingEvent = false;
273 pendingEventTime = null;
274 return false;
275 }
276
277 /** Get the occurrence time of the event triggered in the current step.
278 * @return occurrence time of the event triggered in the current
279 * step or infinity if no events are triggered
280 */
281 public T getEventTime() {
282 return pendingEvent ?
283 pendingEventTime :
284 t0.getField().getZero().add(forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
285 }
286
287 /** Acknowledge the fact the step has been accepted by the integrator.
288 * @param state state at the end of the step
289 */
290 public void stepAccepted(final FieldODEStateAndDerivative<T> state) {
291
292 t0 = state.getTime();
293 g0 = handler.g(state);
294
295 if (pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0) {
296 // force the sign to its value "just after the event"
297 previousEventTime = state.getTime();
298 g0Positive = increasing;
299 nextAction = handler.eventOccurred(state, increasing == forward);
300 } else {
301 g0Positive = g0.getReal() >= 0;
302 nextAction = Action.CONTINUE;
303 }
304 }
305
306 /** Check if the integration should be stopped at the end of the
307 * current step.
308 * @return true if the integration should be stopped
309 */
310 public boolean stop() {
311 return nextAction == Action.STOP;
312 }
313
314 /** Let the event handler reset the state if it wants.
315 * @param state state at the beginning of the next step
316 * @return reset state (may by the same as initial state if only
317 * derivatives should be reset), or null if nothing is reset
318 */
319 public FieldODEState<T> reset(final FieldODEStateAndDerivative<T> state) {
320
321 if (!(pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0)) {
322 return null;
323 }
324
325 final FieldODEState<T> newState;
326 if (nextAction == Action.RESET_STATE) {
327 newState = handler.resetState(state);
328 } else if (nextAction == Action.RESET_DERIVATIVES) {
329 newState = state;
330 } else {
331 newState = null;
332 }
333 pendingEvent = false;
334 pendingEventTime = null;
335
336 return newState;
337 }
338 }