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    }