001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math.ode.events;
019
020 import org.apache.commons.math.ConvergenceException;
021 import org.apache.commons.math.FunctionEvaluationException;
022 import org.apache.commons.math.analysis.UnivariateRealFunction;
023 import org.apache.commons.math.analysis.solvers.BrentSolver;
024 import org.apache.commons.math.exception.MathInternalError;
025 import org.apache.commons.math.ode.DerivativeException;
026 import org.apache.commons.math.ode.sampling.StepInterpolator;
027 import org.apache.commons.math.util.FastMath;
028
029 /** This class handles the state for one {@link EventHandler
030 * event handler} during integration steps.
031 *
032 * <p>Each time the integrator proposes a step, the event handler
033 * switching function should be checked. This class handles the state
034 * of one handler during one integration step, with references to the
035 * state at the end of the preceding step. This information is used to
036 * decide if the handler should trigger an event or not during the
037 * proposed step (and hence the step should be reduced to ensure the
038 * event occurs at a bound rather than inside the step).</p>
039 *
040 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 f??vr. 2011) $
041 * @since 1.2
042 */
043 public class EventState {
044
045 /** Event handler. */
046 private final EventHandler handler;
047
048 /** Maximal time interval between events handler checks. */
049 private final double maxCheckInterval;
050
051 /** Convergence threshold for event localization. */
052 private final double convergence;
053
054 /** Upper limit in the iteration count for event localization. */
055 private final int maxIterationCount;
056
057 /** Time at the beginning of the step. */
058 private double t0;
059
060 /** Value of the events handler at the beginning of the step. */
061 private double g0;
062
063 /** Simulated sign of g0 (we cheat when crossing events). */
064 private boolean g0Positive;
065
066 /** Indicator of event expected during the step. */
067 private boolean pendingEvent;
068
069 /** Occurrence time of the pending event. */
070 private double pendingEventTime;
071
072 /** Occurrence time of the previous event. */
073 private double previousEventTime;
074
075 /** Integration direction. */
076 private boolean forward;
077
078 /** Variation direction around pending event.
079 * (this is considered with respect to the integration direction)
080 */
081 private boolean increasing;
082
083 /** Next action indicator. */
084 private int nextAction;
085
086 /** Simple constructor.
087 * @param handler event handler
088 * @param maxCheckInterval maximal time interval between switching
089 * function checks (this interval prevents missing sign changes in
090 * case the integration steps becomes very large)
091 * @param convergence convergence threshold in the event time search
092 * @param maxIterationCount upper limit of the iteration count in
093 * the event time search
094 */
095 public EventState(final EventHandler handler, final double maxCheckInterval,
096 final double convergence, final int maxIterationCount) {
097 this.handler = handler;
098 this.maxCheckInterval = maxCheckInterval;
099 this.convergence = FastMath.abs(convergence);
100 this.maxIterationCount = maxIterationCount;
101
102 // some dummy values ...
103 t0 = Double.NaN;
104 g0 = Double.NaN;
105 g0Positive = true;
106 pendingEvent = false;
107 pendingEventTime = Double.NaN;
108 previousEventTime = Double.NaN;
109 increasing = true;
110 nextAction = EventHandler.CONTINUE;
111
112 }
113
114 /** Get the underlying event handler.
115 * @return underlying event handler
116 */
117 public EventHandler getEventHandler() {
118 return handler;
119 }
120
121 /** Get the maximal time interval between events handler checks.
122 * @return maximal time interval between events handler checks
123 */
124 public double getMaxCheckInterval() {
125 return maxCheckInterval;
126 }
127
128 /** Get the convergence threshold for event localization.
129 * @return convergence threshold for event localization
130 */
131 public double getConvergence() {
132 return convergence;
133 }
134
135 /** Get the upper limit in the iteration count for event localization.
136 * @return upper limit in the iteration count for event localization
137 */
138 public int getMaxIterationCount() {
139 return maxIterationCount;
140 }
141
142 /** Reinitialize the beginning of the step.
143 * @param interpolator valid for the current step
144 * @exception EventException if the event handler
145 * value cannot be evaluated at the beginning of the step
146 */
147 public void reinitializeBegin(final StepInterpolator interpolator)
148 throws EventException {
149 try {
150 // excerpt from MATH-421 issue:
151 // If an ODE solver is setup with an EventHandler that return STOP
152 // when the even is triggered, the integrator stops (which is exactly
153 // the expected behavior). If however the user want to restart the
154 // solver from the final state reached at the event with the same
155 // configuration (expecting the event to be triggered again at a
156 // later time), then the integrator may fail to start. It can get stuck
157 // at the previous event.
158
159 // The use case for the bug MATH-421 is fairly general, so events occurring
160 // less than epsilon after the solver start in the first step should be ignored,
161 // where epsilon is the convergence threshold of the event. The sign of the g
162 // function should be evaluated after this initial ignore zone, not exactly at
163 // beginning (if there are no event at the very beginning g(t0) and g(t0+epsilon)
164 // have the same sign, so this does not hurt ; if there is an event at the very
165 // beginning, g(t0) and g(t0+epsilon) have opposite signs and we want to start
166 // with the second one. Of course, the sign of epsilon depend on the integration
167 // direction (forward or backward). This explains what is done below.
168
169 final double ignoreZone = interpolator.isForward() ? getConvergence() : -getConvergence();
170 t0 = interpolator.getPreviousTime() + ignoreZone;
171 interpolator.setInterpolatedTime(t0);
172 g0 = handler.g(t0, interpolator.getInterpolatedState());
173 if (g0 == 0) {
174 // extremely rare case: there is a zero EXACTLY at end of ignore zone
175 // we will use the opposite of sign at step beginning to force ignoring this zero
176 final double tStart = interpolator.getPreviousTime();
177 interpolator.setInterpolatedTime(tStart);
178 g0Positive = handler.g(tStart, interpolator.getInterpolatedState()) <= 0;
179 } else {
180 g0Positive = g0 >= 0;
181 }
182
183 } catch (DerivativeException mue) {
184 throw new EventException(mue);
185 }
186 }
187
188 /** Evaluate the impact of the proposed step on the event handler.
189 * @param interpolator step interpolator for the proposed step
190 * @return true if the event handler triggers an event before
191 * the end of the proposed step
192 * @exception DerivativeException if the interpolator fails to
193 * compute the switching function somewhere within the step
194 * @exception EventException if the switching function
195 * cannot be evaluated
196 * @exception ConvergenceException if an event cannot be located
197 */
198 public boolean evaluateStep(final StepInterpolator interpolator)
199 throws DerivativeException, EventException, ConvergenceException {
200
201 try {
202
203 forward = interpolator.isForward();
204 final double t1 = interpolator.getCurrentTime();
205 if (FastMath.abs(t1 - t0) < convergence) {
206 // we cannot do anything on such a small step, don't trigger any events
207 return false;
208 }
209 final double start = forward ? (t0 + convergence) : t0 - convergence;
210 final double dt = t1 - start;
211 final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / maxCheckInterval));
212 final double h = dt / n;
213
214 double ta = t0;
215 double ga = g0;
216 for (int i = 0; i < n; ++i) {
217
218 // evaluate handler value at the end of the substep
219 final double tb = start + (i + 1) * h;
220 interpolator.setInterpolatedTime(tb);
221 final double gb = handler.g(tb, interpolator.getInterpolatedState());
222
223 // check events occurrence
224 if (g0Positive ^ (gb >= 0)) {
225 // there is a sign change: an event is expected during this step
226
227 // variation direction, with respect to the integration direction
228 increasing = gb >= ga;
229
230 final UnivariateRealFunction f = new UnivariateRealFunction() {
231 public double value(final double t) {
232 try {
233 interpolator.setInterpolatedTime(t);
234 return handler.g(t, interpolator.getInterpolatedState());
235 } catch (DerivativeException e) {
236 throw new EmbeddedDerivativeException(e);
237 } catch (EventException e) {
238 throw new EmbeddedEventException(e);
239 }
240 }
241 };
242 final BrentSolver solver = new BrentSolver(convergence);
243
244 if (ga * gb >= 0) {
245 // this is a corner case:
246 // - there was an event near ta,
247 // - there is another event between ta and tb
248 // - when ta was computed, convergence was reached on the "wrong side" of the interval
249 // this implies that the real sign of ga is the same as gb, so we need to slightly
250 // shift ta to make sure ga and gb get opposite signs and the solver won't complain
251 // about bracketing
252 final double epsilon = (forward ? 0.25 : -0.25) * convergence;
253 for (int k = 0; (k < 4) && (ga * gb > 0); ++k) {
254 ta += epsilon;
255 try {
256 ga = f.value(ta);
257 } catch (FunctionEvaluationException ex) {
258 throw new DerivativeException(ex);
259 }
260 }
261 if (ga * gb > 0) {
262 // this should never happen
263 throw new MathInternalError();
264 }
265 }
266
267 final double root;
268 try {
269 root = (ta <= tb) ?
270 solver.solve(maxIterationCount, f, ta, tb) :
271 solver.solve(maxIterationCount, f, tb, ta);
272 } catch (FunctionEvaluationException ex) {
273 throw new DerivativeException(ex);
274 }
275
276 if ((!Double.isNaN(previousEventTime)) &&
277 (FastMath.abs(root - ta) <= convergence) &&
278 (FastMath.abs(root - previousEventTime) <= convergence)) {
279 // we have either found nothing or found (again ?) a past event, we simply ignore it
280 ta = tb;
281 ga = gb;
282 } else if (Double.isNaN(previousEventTime) ||
283 (FastMath.abs(previousEventTime - root) > convergence)) {
284 pendingEventTime = root;
285 pendingEvent = true;
286 return true;
287 } else {
288 // no sign change: there is no event for now
289 ta = tb;
290 ga = gb;
291 }
292
293 } else {
294 // no sign change: there is no event for now
295 ta = tb;
296 ga = gb;
297 }
298
299 }
300
301 // no event during the whole step
302 pendingEvent = false;
303 pendingEventTime = Double.NaN;
304 return false;
305
306 } catch (EmbeddedDerivativeException ede) {
307 throw ede.getDerivativeException();
308 } catch (EmbeddedEventException eee) {
309 throw eee.getEventException();
310 }
311
312 }
313
314 /** Get the occurrence time of the event triggered in the current step.
315 * @return occurrence time of the event triggered in the current
316 * step or positive infinity if no events are triggered
317 */
318 public double getEventTime() {
319 return pendingEvent ? pendingEventTime : Double.POSITIVE_INFINITY;
320 }
321
322 /** Acknowledge the fact the step has been accepted by the integrator.
323 * @param t value of the independent <i>time</i> variable at the
324 * end of the step
325 * @param y array containing the current value of the state vector
326 * at the end of the step
327 * @exception EventException if the value of the event
328 * handler cannot be evaluated
329 */
330 public void stepAccepted(final double t, final double[] y)
331 throws EventException {
332
333 t0 = t;
334 g0 = handler.g(t, y);
335
336 if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) {
337 // force the sign to its value "just after the event"
338 previousEventTime = t;
339 g0Positive = increasing;
340 nextAction = handler.eventOccurred(t, y, !(increasing ^ forward));
341 } else {
342 g0Positive = g0 >= 0;
343 nextAction = EventHandler.CONTINUE;
344 }
345 }
346
347 /** Check if the integration should be stopped at the end of the
348 * current step.
349 * @return true if the integration should be stopped
350 */
351 public boolean stop() {
352 return nextAction == EventHandler.STOP;
353 }
354
355 /** Let the event handler reset the state if it wants.
356 * @param t value of the independent <i>time</i> variable at the
357 * beginning of the next step
358 * @param y array were to put the desired state vector at the beginning
359 * of the next step
360 * @return true if the integrator should reset the derivatives too
361 * @exception EventException if the state cannot be reseted by the event
362 * handler
363 */
364 public boolean reset(final double t, final double[] y)
365 throws EventException {
366
367 if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) {
368 return false;
369 }
370
371 if (nextAction == EventHandler.RESET_STATE) {
372 handler.resetState(t, y);
373 }
374 pendingEvent = false;
375 pendingEventTime = Double.NaN;
376
377 return (nextAction == EventHandler.RESET_STATE) ||
378 (nextAction == EventHandler.RESET_DERIVATIVES);
379
380 }
381
382 /** Local exception for embedding DerivativeException. */
383 private static class EmbeddedDerivativeException extends RuntimeException {
384
385 /** Serializable UID. */
386 private static final long serialVersionUID = 3574188382434584610L;
387
388 /** Embedded exception. */
389 private final DerivativeException derivativeException;
390
391 /** Simple constructor.
392 * @param derivativeException embedded exception
393 */
394 public EmbeddedDerivativeException(final DerivativeException derivativeException) {
395 this.derivativeException = derivativeException;
396 }
397
398 /** Get the embedded exception.
399 * @return embedded exception
400 */
401 public DerivativeException getDerivativeException() {
402 return derivativeException;
403 }
404
405 }
406
407 /** Local exception for embedding EventException. */
408 private static class EmbeddedEventException extends RuntimeException {
409
410 /** Serializable UID. */
411 private static final long serialVersionUID = -1337749250090455474L;
412
413 /** Embedded exception. */
414 private final EventException eventException;
415
416 /** Simple constructor.
417 * @param eventException embedded exception
418 */
419 public EmbeddedEventException(final EventException eventException) {
420 this.eventException = eventException;
421 }
422
423 /** Get the embedded exception.
424 * @return embedded exception
425 */
426 public EventException getEventException() {
427 return eventException;
428 }
429
430 }
431 }