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