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.math.ode.events;
19
20 import org.apache.commons.math.analysis.UnivariateRealFunction;
21 import org.apache.commons.math.analysis.solvers.AllowedSolution;
22 import org.apache.commons.math.analysis.solvers.BracketedUnivariateRealSolver;
23 import org.apache.commons.math.analysis.solvers.PegasusSolver;
24 import org.apache.commons.math.analysis.solvers.UnivariateRealSolver;
25 import org.apache.commons.math.analysis.solvers.UnivariateRealSolverUtils;
26 import org.apache.commons.math.exception.ConvergenceException;
27 import org.apache.commons.math.ode.events.EventHandler;
28 import org.apache.commons.math.ode.sampling.StepInterpolator;
29 import org.apache.commons.math.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 1175379 2011-09-25 12:39:09Z luc $
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 UnivariateRealSolver 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 UnivariateRealSolver 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 */
152 public void reinitializeBegin(final StepInterpolator interpolator) {
153
154 t0 = interpolator.getPreviousTime();
155 interpolator.setInterpolatedTime(t0);
156 g0 = handler.g(t0, interpolator.getInterpolatedState());
157 if (g0 == 0) {
158 // excerpt from MATH-421 issue:
159 // If an ODE solver is setup with an EventHandler that return STOP
160 // when the even is triggered, the integrator stops (which is exactly
161 // the expected behavior). If however the user wants to restart the
162 // solver from the final state reached at the event with the same
163 // configuration (expecting the event to be triggered again at a
164 // later time), then the integrator may fail to start. It can get stuck
165 // at the previous event. The use case for the bug MATH-421 is fairly
166 // general, so events occurring exactly at start in the first step should
167 // be ignored.
168
169 // extremely rare case: there is a zero EXACTLY at interval start
170 // we will use the sign slightly after step beginning to force ignoring this zero
171 final double epsilon = FastMath.max(solver.getAbsoluteAccuracy(),
172 FastMath.abs(solver.getRelativeAccuracy() * t0));
173 final double tStart = t0 + 0.5 * epsilon;
174 interpolator.setInterpolatedTime(tStart);
175 g0 = handler.g(tStart, interpolator.getInterpolatedState());
176 }
177 g0Positive = g0 >= 0;
178
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 ConvergenceException if an event cannot be located
186 */
187 public boolean evaluateStep(final StepInterpolator interpolator)
188 throws ConvergenceException {
189
190 forward = interpolator.isForward();
191 final double t1 = interpolator.getCurrentTime();
192 final double dt = t1 - t0;
193 if (FastMath.abs(dt) < convergence) {
194 // we cannot do anything on such a small step, don't trigger any events
195 return false;
196 }
197 final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / maxCheckInterval));
198 final double h = dt / n;
199
200 final UnivariateRealFunction f = new UnivariateRealFunction() {
201 public double value(final double t) {
202 interpolator.setInterpolatedTime(t);
203 return handler.g(t, interpolator.getInterpolatedState());
204 }
205 };
206
207 double ta = t0;
208 double ga = g0;
209 for (int i = 0; i < n; ++i) {
210
211 // evaluate handler value at the end of the substep
212 final double tb = t0 + (i + 1) * h;
213 interpolator.setInterpolatedTime(tb);
214 final double gb = handler.g(tb, interpolator.getInterpolatedState());
215
216 // check events occurrence
217 if (g0Positive ^ (gb >= 0)) {
218 // there is a sign change: an event is expected during this step
219
220 // variation direction, with respect to the integration direction
221 increasing = gb >= ga;
222
223 // find the event time making sure we select a solution just at or past the exact root
224 final double root;
225 if (solver instanceof BracketedUnivariateRealSolver<?>) {
226 @SuppressWarnings("unchecked")
227 BracketedUnivariateRealSolver<UnivariateRealFunction> bracketing =
228 (BracketedUnivariateRealSolver<UnivariateRealFunction>) solver;
229 root = forward ?
230 bracketing.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
231 bracketing.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
232 } else {
233 final double baseRoot = forward ?
234 solver.solve(maxIterationCount, f, ta, tb) :
235 solver.solve(maxIterationCount, f, tb, ta);
236 final int remainingEval = maxIterationCount - solver.getEvaluations();
237 BracketedUnivariateRealSolver<UnivariateRealFunction> bracketing =
238 new PegasusSolver(solver.getRelativeAccuracy(), solver.getAbsoluteAccuracy());
239 root = forward ?
240 UnivariateRealSolverUtils.forceSide(remainingEval, f, bracketing,
241 baseRoot, ta, tb, AllowedSolution.RIGHT_SIDE) :
242 UnivariateRealSolverUtils.forceSide(remainingEval, f, bracketing,
243 baseRoot, tb, ta, AllowedSolution.LEFT_SIDE);
244 }
245
246 if ((!Double.isNaN(previousEventTime)) &&
247 (FastMath.abs(root - ta) <= convergence) &&
248 (FastMath.abs(root - previousEventTime) <= convergence)) {
249 // we have either found nothing or found (again ?) a past event,
250 // retry the substep excluding this value
251 ta = forward ? ta + convergence : ta - convergence;
252 ga = f.value(ta);
253 --i;
254 } else if (Double.isNaN(previousEventTime) ||
255 (FastMath.abs(previousEventTime - root) > convergence)) {
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
265 } else {
266 // no sign change: there is no event for now
267 ta = tb;
268 ga = gb;
269 }
270
271 }
272
273 // no event during the whole step
274 pendingEvent = false;
275 pendingEventTime = Double.NaN;
276 return false;
277
278 }
279
280 /** Get the occurrence time of the event triggered in the current step.
281 * @return occurrence time of the event triggered in the current
282 * step or infinity if no events are triggered
283 */
284 public double getEventTime() {
285 return pendingEvent ?
286 pendingEventTime :
287 (forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
288 }
289
290 /** Acknowledge the fact the step has been accepted by the integrator.
291 * @param t value of the independent <i>time</i> variable at the
292 * end of the step
293 * @param y array containing the current value of the state vector
294 * at the end of the step
295 */
296 public void stepAccepted(final double t, final double[] y) {
297
298 t0 = t;
299 g0 = handler.g(t, y);
300
301 if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) {
302 // force the sign to its value "just after the event"
303 previousEventTime = t;
304 g0Positive = increasing;
305 nextAction = handler.eventOccurred(t, y, !(increasing ^ forward));
306 } else {
307 g0Positive = g0 >= 0;
308 nextAction = EventHandler.Action.CONTINUE;
309 }
310 }
311
312 /** Check if the integration should be stopped at the end of the
313 * current step.
314 * @return true if the integration should be stopped
315 */
316 public boolean stop() {
317 return nextAction == EventHandler.Action.STOP;
318 }
319
320 /** Let the event handler reset the state if it wants.
321 * @param t value of the independent <i>time</i> variable at the
322 * beginning of the next step
323 * @param y array were to put the desired state vector at the beginning
324 * of the next step
325 * @return true if the integrator should reset the derivatives too
326 */
327 public boolean reset(final double t, final double[] y) {
328
329 if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) {
330 return false;
331 }
332
333 if (nextAction == EventHandler.Action.RESET_STATE) {
334 handler.resetState(t, y);
335 }
336 pendingEvent = false;
337 pendingEventTime = Double.NaN;
338
339 return (nextAction == EventHandler.Action.RESET_STATE) ||
340 (nextAction == EventHandler.Action.RESET_DERIVATIVES);
341
342 }
343
344 }