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.MathRuntimeException;
023    import org.apache.commons.math.analysis.UnivariateRealFunction;
024    import org.apache.commons.math.analysis.solvers.BrentSolver;
025    import org.apache.commons.math.ode.DerivativeException;
026    import org.apache.commons.math.ode.sampling.StepInterpolator;
027    
028    /** This class handles the state for one {@link EventHandler
029     * event handler} during integration steps.
030     *
031     * <p>Each time the integrator proposes a step, the event handler
032     * switching function should be checked. This class handles the state
033     * of one handler during one integration step, with references to the
034     * state at the end of the preceding step. This information is used to
035     * decide if the handler should trigger an event or not during the
036     * proposed step (and hence the step should be reduced to ensure the
037     * event occurs at a bound rather than inside the step).</p>
038     *
039     * @version $Revision: 889006 $ $Date: 2009-12-09 17:46:36 -0500 (Wed, 09 Dec 2009) $
040     * @since 1.2
041     */
042    public class EventState {
043    
044        /** Event handler. */
045        private final EventHandler handler;
046    
047        /** Maximal time interval between events handler checks. */
048        private final double maxCheckInterval;
049    
050        /** Convergence threshold for event localization. */
051        private final double convergence;
052    
053        /** Upper limit in the iteration count for event localization. */
054        private final int maxIterationCount;
055    
056        /** Time at the beginning of the step. */
057        private double t0;
058    
059        /** Value of the events handler at the beginning of the step. */
060        private double g0;
061    
062        /** Simulated sign of g0 (we cheat when crossing events). */
063        private boolean g0Positive;
064    
065        /** Indicator of event expected during the step. */
066        private boolean pendingEvent;
067    
068        /** Occurrence time of the pending event. */
069        private double pendingEventTime;
070    
071        /** Occurrence time of the previous event. */
072        private double previousEventTime;
073    
074        /** Integration direction. */
075        private boolean forward;
076    
077        /** Variation direction around pending event.
078         *  (this is considered with respect to the integration direction)
079         */
080        private boolean increasing;
081    
082        /** Next action indicator. */
083        private int nextAction;
084    
085        /** Simple constructor.
086         * @param handler event handler
087         * @param maxCheckInterval maximal time interval between switching
088         * function checks (this interval prevents missing sign changes in
089         * case the integration steps becomes very large)
090         * @param convergence convergence threshold in the event time search
091         * @param maxIterationCount upper limit of the iteration count in
092         * the event time search
093         */
094        public EventState(final EventHandler handler, final double maxCheckInterval,
095                          final double convergence, final int maxIterationCount) {
096            this.handler           = handler;
097            this.maxCheckInterval  = maxCheckInterval;
098            this.convergence       = Math.abs(convergence);
099            this.maxIterationCount = maxIterationCount;
100    
101            // some dummy values ...
102            t0                = Double.NaN;
103            g0                = Double.NaN;
104            g0Positive        = true;
105            pendingEvent      = false;
106            pendingEventTime  = Double.NaN;
107            previousEventTime = Double.NaN;
108            increasing        = true;
109            nextAction        = EventHandler.CONTINUE;
110    
111        }
112    
113        /** Get the underlying event handler.
114         * @return underlying event handler
115         */
116        public EventHandler getEventHandler() {
117            return handler;
118        }
119    
120        /** Get the maximal time interval between events handler checks.
121         * @return maximal time interval between events handler checks
122         */
123        public double getMaxCheckInterval() {
124            return maxCheckInterval;
125        }
126    
127        /** Get the convergence threshold for event localization.
128         * @return convergence threshold for event localization
129         */
130        public double getConvergence() {
131            return convergence;
132        }
133    
134        /** Get the upper limit in the iteration count for event localization.
135         * @return upper limit in the iteration count for event localization
136         */
137        public int getMaxIterationCount() {
138            return maxIterationCount;
139        }
140    
141        /** Reinitialize the beginning of the step.
142         * @param tStart value of the independent <i>time</i> variable at the
143         * beginning of the step
144         * @param yStart array containing the current value of the state vector
145         * at the beginning of the step
146         * @exception EventException if the event handler
147         * value cannot be evaluated at the beginning of the step
148         */
149        public void reinitializeBegin(final double tStart, final double[] yStart)
150            throws EventException {
151            t0 = tStart;
152            g0 = handler.g(tStart, yStart);
153            g0Positive = g0 >= 0;
154        }
155    
156        /** Evaluate the impact of the proposed step on the event handler.
157         * @param interpolator step interpolator for the proposed step
158         * @return true if the event handler triggers an event before
159         * the end of the proposed step (this implies the step should be
160         * rejected)
161         * @exception DerivativeException if the interpolator fails to
162         * compute the switching function somewhere within the step
163         * @exception EventException if the switching function
164         * cannot be evaluated
165         * @exception ConvergenceException if an event cannot be located
166         */
167        public boolean evaluateStep(final StepInterpolator interpolator)
168            throws DerivativeException, EventException, ConvergenceException {
169    
170            try {
171    
172                forward = interpolator.isForward();
173                final double t1 = interpolator.getCurrentTime();
174                final int    n  = Math.max(1, (int) Math.ceil(Math.abs(t1 - t0) / maxCheckInterval));
175                final double h  = (t1 - t0) / n;
176    
177                double ta = t0;
178                double ga = g0;
179                double tb = t0 + (interpolator.isForward() ? convergence : -convergence);
180                for (int i = 0; i < n; ++i) {
181    
182                    // evaluate handler value at the end of the substep
183                    tb += h;
184                    interpolator.setInterpolatedTime(tb);
185                    final double gb = handler.g(tb, interpolator.getInterpolatedState());
186    
187                    // check events occurrence
188                    if (g0Positive ^ (gb >= 0)) {
189                        // there is a sign change: an event is expected during this step
190    
191                        if (ga * gb > 0) {
192                            // this is a corner case:
193                            // - there was an event near ta,
194                            // - there is another event between ta and tb
195                            // - when ta was computed, convergence was reached on the "wrong side" of the interval
196                            // this implies that the real sign of ga is the same as gb, so we need to slightly
197                            // shift ta to make sure ga and gb get opposite signs and the solver won't complain
198                            // about bracketing
199                            final double epsilon = (forward ? 0.25 : -0.25) * convergence;
200                            for (int k = 0; (k < 4) && (ga * gb > 0); ++k) {
201                                ta += epsilon;
202                                interpolator.setInterpolatedTime(ta);
203                                ga = handler.g(ta, interpolator.getInterpolatedState());
204                            }
205                            if (ga * gb > 0) {
206                                // this should never happen
207                                throw MathRuntimeException.createInternalError(null);
208                            }
209                        }
210    
211                        // variation direction, with respect to the integration direction
212                        increasing = gb >= ga;
213    
214                        final UnivariateRealFunction f = new UnivariateRealFunction() {
215                            public double value(final double t) throws FunctionEvaluationException {
216                                try {
217                                    interpolator.setInterpolatedTime(t);
218                                    return handler.g(t, interpolator.getInterpolatedState());
219                                } catch (DerivativeException e) {
220                                    throw new FunctionEvaluationException(e, t);
221                                } catch (EventException e) {
222                                    throw new FunctionEvaluationException(e, t);
223                                }
224                            }
225                        };
226                        final BrentSolver solver = new BrentSolver();
227                        solver.setAbsoluteAccuracy(convergence);
228                        solver.setMaximalIterationCount(maxIterationCount);
229                        final double root = (ta <= tb) ? solver.solve(f, ta, tb) : solver.solve(f, tb, ta);
230                        if ((Math.abs(root - ta) <= convergence) &&
231                             (Math.abs(root - previousEventTime) <= convergence)) {
232                            // we have either found nothing or found (again ?) a past event, we simply ignore it
233                            ta = tb;
234                            ga = gb;
235                        } else if (Double.isNaN(previousEventTime) ||
236                                   (Math.abs(previousEventTime - root) > convergence)) {
237                            pendingEventTime = root;
238                            if (pendingEvent && (Math.abs(t1 - pendingEventTime) <= convergence)) {
239                                // we were already waiting for this event which was
240                                // found during a previous call for a step that was
241                                // rejected, this step must now be accepted since it
242                                // properly ends exactly at the event occurrence
243                                return false;
244                            }
245                            // either we were not waiting for the event or it has
246                            // moved in such a way the step cannot be accepted
247                            pendingEvent = true;
248                            return true;
249                        }
250    
251                    } else {
252                        // no sign change: there is no event for now
253                        ta = tb;
254                        ga = gb;
255                    }
256    
257                }
258    
259                // no event during the whole step
260                pendingEvent     = false;
261                pendingEventTime = Double.NaN;
262                return false;
263    
264            } catch (FunctionEvaluationException e) {
265                final Throwable cause = e.getCause();
266                if ((cause != null) && (cause instanceof DerivativeException)) {
267                    throw (DerivativeException) cause;
268                } else if ((cause != null) && (cause instanceof EventException)) {
269                    throw (EventException) cause;
270                }
271                throw new EventException(e);
272            }
273    
274        }
275    
276        /** Get the occurrence time of the event triggered in the current
277         * step.
278         * @return occurrence time of the event triggered in the current
279         * step.
280         */
281        public double getEventTime() {
282            return pendingEventTime;
283        }
284    
285        /** Acknowledge the fact the step has been accepted by the integrator.
286         * @param t value of the independent <i>time</i> variable at the
287         * end of the step
288         * @param y array containing the current value of the state vector
289         * at the end of the step
290         * @exception EventException if the value of the event
291         * handler cannot be evaluated
292         */
293        public void stepAccepted(final double t, final double[] y)
294            throws EventException {
295    
296            t0 = t;
297            g0 = handler.g(t, y);
298    
299            if (pendingEvent) {
300                // force the sign to its value "just after the event"
301                previousEventTime = t;
302                g0Positive        = increasing;
303                nextAction        = handler.eventOccurred(t, y, !(increasing ^ forward));
304            } else {
305                g0Positive = g0 >= 0;
306                nextAction = EventHandler.CONTINUE;
307            }
308        }
309    
310        /** Check if the integration should be stopped at the end of the
311         * current step.
312         * @return true if the integration should be stopped
313         */
314        public boolean stop() {
315            return nextAction == EventHandler.STOP;
316        }
317    
318        /** Let the event handler reset the state if it wants.
319         * @param t value of the independent <i>time</i> variable at the
320         * beginning of the next step
321         * @param y array were to put the desired state vector at the beginning
322         * of the next step
323         * @return true if the integrator should reset the derivatives too
324         * @exception EventException if the state cannot be reseted by the event
325         * handler
326         */
327        public boolean reset(final double t, final double[] y)
328            throws EventException {
329    
330            if (! pendingEvent) {
331                return false;
332            }
333    
334            if (nextAction == EventHandler.RESET_STATE) {
335                handler.resetState(t, y);
336            }
337            pendingEvent      = false;
338            pendingEventTime  = Double.NaN;
339    
340            return (nextAction == EventHandler.RESET_STATE) ||
341                   (nextAction == EventHandler.RESET_DERIVATIVES);
342    
343        }
344    
345    }