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.math3.ode.events;
19
20 import org.apache.commons.math3.analysis.UnivariateFunction;
21 import org.apache.commons.math3.analysis.solvers.AllowedSolution;
22 import org.apache.commons.math3.analysis.solvers.BracketedUnivariateSolver;
23 import org.apache.commons.math3.analysis.solvers.PegasusSolver;
24 import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
25 import org.apache.commons.math3.analysis.solvers.UnivariateSolverUtils;
26 import org.apache.commons.math3.exception.MaxCountExceededException;
27 import org.apache.commons.math3.exception.NoBracketingException;
28 import org.apache.commons.math3.ode.sampling.StepInterpolator;
29 import org.apache.commons.math3.util.FastMath;
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 * @version $Id: EventState.java 1416643 2012-12-03 19:37:14Z tn $
42 * @since 1.2
43 */
44 public class EventState {
45
46 /** Event handler. */
47 private final EventHandler 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 double 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 double t0;
60
61 /** Value of the events handler at the beginning of the step. */
62 private double 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 double pendingEventTime;
72
73 /** Occurrence time of the previous event. */
74 private double 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 EventHandler.Action nextAction;
86
87 /** Root-finding algorithm to use to detect state events. */
88 private final UnivariateSolver 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 EventState(final EventHandler handler, final double maxCheckInterval,
101 final double convergence, final int maxIterationCount,
102 final UnivariateSolver solver) {
103 this.handler = handler;
104 this.maxCheckInterval = maxCheckInterval;
105 this.convergence = FastMath.abs(convergence);
106 this.maxIterationCount = maxIterationCount;
107 this.solver = solver;
108
109 // some dummy values ...
110 t0 = Double.NaN;
111 g0 = Double.NaN;
112 g0Positive = true;
113 pendingEvent = false;
114 pendingEventTime = Double.NaN;
115 previousEventTime = Double.NaN;
116 increasing = true;
117 nextAction = EventHandler.Action.CONTINUE;
118
119 }
120
121 /** Get the underlying event handler.
122 * @return underlying event handler
123 */
124 public EventHandler getEventHandler() {
125 return handler;
126 }
127
128 /** Get the maximal time interval between events handler checks.
129 * @return maximal time interval between events handler checks
130 */
131 public double getMaxCheckInterval() {
132 return maxCheckInterval;
133 }
134
135 /** Get the convergence threshold for event localization.
136 * @return convergence threshold for event localization
137 */
138 public double getConvergence() {
139 return convergence;
140 }
141
142 /** Get the upper limit in the iteration count for event localization.
143 * @return upper limit in the iteration count for event localization
144 */
145 public int getMaxIterationCount() {
146 return maxIterationCount;
147 }
148
149 /** Reinitialize the beginning of the step.
150 * @param interpolator valid for the current step
151 * @exception MaxCountExceededException if the interpolator throws one because
152 * the number of functions evaluations is exceeded
153 */
154 public void reinitializeBegin(final StepInterpolator interpolator)
155 throws MaxCountExceededException {
156
157 t0 = interpolator.getPreviousTime();
158 interpolator.setInterpolatedTime(t0);
159 g0 = handler.g(t0, interpolator.getInterpolatedState());
160 if (g0 == 0) {
161 // excerpt from MATH-421 issue:
162 // If an ODE solver is setup with an EventHandler that return STOP
163 // when the even is triggered, the integrator stops (which is exactly
164 // the expected behavior). If however the user wants to restart the
165 // solver from the final state reached at the event with the same
166 // configuration (expecting the event to be triggered again at a
167 // later time), then the integrator may fail to start. It can get stuck
168 // at the previous event. The use case for the bug MATH-421 is fairly
169 // general, so events occurring exactly at start in the first step should
170 // be ignored.
171
172 // extremely rare case: there is a zero EXACTLY at interval start
173 // we will use the sign slightly after step beginning to force ignoring this zero
174 final double epsilon = FastMath.max(solver.getAbsoluteAccuracy(),
175 FastMath.abs(solver.getRelativeAccuracy() * t0));
176 final double tStart = t0 + 0.5 * epsilon;
177 interpolator.setInterpolatedTime(tStart);
178 g0 = handler.g(tStart, interpolator.getInterpolatedState());
179 }
180 g0Positive = g0 >= 0;
181
182 }
183
184 /** Evaluate the impact of the proposed step on the event handler.
185 * @param interpolator step interpolator for the proposed step
186 * @return true if the event handler triggers an event before
187 * the end of the proposed step
188 * @exception MaxCountExceededException if the interpolator throws one because
189 * the number of functions evaluations is exceeded
190 * @exception NoBracketingException if the event cannot be bracketed
191 */
192 public boolean evaluateStep(final StepInterpolator interpolator)
193 throws MaxCountExceededException, NoBracketingException {
194
195 try {
196 forward = interpolator.isForward();
197 final double t1 = interpolator.getCurrentTime();
198 final double dt = t1 - t0;
199 if (FastMath.abs(dt) < convergence) {
200 // we cannot do anything on such a small step, don't trigger any events
201 return false;
202 }
203 final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / maxCheckInterval));
204 final double h = dt / n;
205
206 final UnivariateFunction f = new UnivariateFunction() {
207 public double value(final double t) throws LocalMaxCountExceededException {
208 try {
209 interpolator.setInterpolatedTime(t);
210 return handler.g(t, interpolator.getInterpolatedState());
211 } catch (MaxCountExceededException mcee) {
212 throw new LocalMaxCountExceededException(mcee);
213 }
214 }
215 };
216
217 double ta = t0;
218 double ga = g0;
219 for (int i = 0; i < n; ++i) {
220
221 // evaluate handler value at the end of the substep
222 final double tb = t0 + (i + 1) * h;
223 interpolator.setInterpolatedTime(tb);
224 final double gb = handler.g(tb, interpolator.getInterpolatedState());
225
226 // check events occurrence
227 if (g0Positive ^ (gb >= 0)) {
228 // there is a sign change: an event is expected during this step
229
230 // variation direction, with respect to the integration direction
231 increasing = gb >= ga;
232
233 // find the event time making sure we select a solution just at or past the exact root
234 final double root;
235 if (solver instanceof BracketedUnivariateSolver<?>) {
236 @SuppressWarnings("unchecked")
237 BracketedUnivariateSolver<UnivariateFunction> bracketing =
238 (BracketedUnivariateSolver<UnivariateFunction>) solver;
239 root = forward ?
240 bracketing.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
241 bracketing.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
242 } else {
243 final double baseRoot = forward ?
244 solver.solve(maxIterationCount, f, ta, tb) :
245 solver.solve(maxIterationCount, f, tb, ta);
246 final int remainingEval = maxIterationCount - solver.getEvaluations();
247 BracketedUnivariateSolver<UnivariateFunction> bracketing =
248 new PegasusSolver(solver.getRelativeAccuracy(), solver.getAbsoluteAccuracy());
249 root = forward ?
250 UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
251 baseRoot, ta, tb, AllowedSolution.RIGHT_SIDE) :
252 UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
253 baseRoot, tb, ta, AllowedSolution.LEFT_SIDE);
254 }
255
256 if ((!Double.isNaN(previousEventTime)) &&
257 (FastMath.abs(root - ta) <= convergence) &&
258 (FastMath.abs(root - previousEventTime) <= convergence)) {
259 // we have either found nothing or found (again ?) a past event,
260 // retry the substep excluding this value
261 ta = forward ? ta + convergence : ta - convergence;
262 ga = f.value(ta);
263 --i;
264 } else if (Double.isNaN(previousEventTime) ||
265 (FastMath.abs(previousEventTime - root) > convergence)) {
266 pendingEventTime = root;
267 pendingEvent = true;
268 return true;
269 } else {
270 // no sign change: there is no event for now
271 ta = tb;
272 ga = gb;
273 }
274
275 } else {
276 // no sign change: there is no event for now
277 ta = tb;
278 ga = gb;
279 }
280
281 }
282
283 // no event during the whole step
284 pendingEvent = false;
285 pendingEventTime = Double.NaN;
286 return false;
287
288 } catch (LocalMaxCountExceededException lmcee) {
289 throw lmcee.getException();
290 }
291
292 }
293
294 /** Get the occurrence time of the event triggered in the current step.
295 * @return occurrence time of the event triggered in the current
296 * step or infinity if no events are triggered
297 */
298 public double getEventTime() {
299 return pendingEvent ?
300 pendingEventTime :
301 (forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
302 }
303
304 /** Acknowledge the fact the step has been accepted by the integrator.
305 * @param t value of the independent <i>time</i> variable at the
306 * end of the step
307 * @param y array containing the current value of the state vector
308 * at the end of the step
309 */
310 public void stepAccepted(final double t, final double[] y) {
311
312 t0 = t;
313 g0 = handler.g(t, y);
314
315 if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) {
316 // force the sign to its value "just after the event"
317 previousEventTime = t;
318 g0Positive = increasing;
319 nextAction = handler.eventOccurred(t, y, !(increasing ^ forward));
320 } else {
321 g0Positive = g0 >= 0;
322 nextAction = EventHandler.Action.CONTINUE;
323 }
324 }
325
326 /** Check if the integration should be stopped at the end of the
327 * current step.
328 * @return true if the integration should be stopped
329 */
330 public boolean stop() {
331 return nextAction == EventHandler.Action.STOP;
332 }
333
334 /** Let the event handler reset the state if it wants.
335 * @param t value of the independent <i>time</i> variable at the
336 * beginning of the next step
337 * @param y array were to put the desired state vector at the beginning
338 * of the next step
339 * @return true if the integrator should reset the derivatives too
340 */
341 public boolean reset(final double t, final double[] y) {
342
343 if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) {
344 return false;
345 }
346
347 if (nextAction == EventHandler.Action.RESET_STATE) {
348 handler.resetState(t, y);
349 }
350 pendingEvent = false;
351 pendingEventTime = Double.NaN;
352
353 return (nextAction == EventHandler.Action.RESET_STATE) ||
354 (nextAction == EventHandler.Action.RESET_DERIVATIVES);
355
356 }
357
358 /** Local wrapper to propagate exceptions. */
359 private static class LocalMaxCountExceededException extends RuntimeException {
360
361 /** Serializable UID. */
362 private static final long serialVersionUID = 20120901L;
363
364 /** Wrapped exception. */
365 private final MaxCountExceededException wrapped;
366
367 /** Simple constructor.
368 * @param exception exception to wrap
369 */
370 public LocalMaxCountExceededException(final MaxCountExceededException exception) {
371 wrapped = exception;
372 }
373
374 /** Get the wrapped exception.
375 * @return wrapped exception
376 */
377 public MaxCountExceededException getException() {
378 return wrapped;
379 }
380
381 }
382
383 }