    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.math3.ode.events;
    019
    020    import org.apache.commons.math3.analysis.UnivariateFunction;
    021    import org.apache.commons.math3.analysis.solvers.AllowedSolution;
    022    import org.apache.commons.math3.analysis.solvers.BracketedUnivariateSolver;
    023    import org.apache.commons.math3.analysis.solvers.PegasusSolver;
    024    import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
    025    import org.apache.commons.math3.analysis.solvers.UnivariateSolverUtils;
    026    import org.apache.commons.math3.exception.MaxCountExceededException;
    027    import org.apache.commons.math3.exception.NoBracketingException;
    028    import org.apache.commons.math3.ode.sampling.StepInterpolator;
    029    import org.apache.commons.math3.util.FastMath;
    030
    031    /** This class handles the state for one {@link EventHandler
    032     * event handler} during integration steps.
    033     *
    034     * <p>Each time the integrator proposes a step, the event handler
    035     * switching function should be checked. This class handles the state
    036     * of one handler during one integration step, with references to the
    037     * state at the end of the preceding step. This information is used to
    038     * decide if the handler should trigger an event or not during the
    039     * proposed step.</p>
    040     *
    041     * @version $Id: EventState.java 1416643 2012-12-03 19:37:14Z tn $
    042     * @since 1.2
    043     */
    044    public class EventState {
    045
    046        /** Event handler. */
    047        private final EventHandler handler;
    048
    049        /** Maximal time interval between events handler checks. */
    050        private final double maxCheckInterval;
    051
    052        /** Convergence threshold for event localization. */
    053        private final double convergence;
    054
    055        /** Upper limit in the iteration count for event localization. */
    056        private final int maxIterationCount;
    057
    058        /** Time at the beginning of the step. */
    059        private double t0;
    060
    061        /** Value of the events handler at the beginning of the step. */
    062        private double g0;
    063
    064        /** Simulated sign of g0 (we cheat when crossing events). */
    065        private boolean g0Positive;
    066
    067        /** Indicator of event expected during the step. */
    068        private boolean pendingEvent;
    069
    070        /** Occurrence time of the pending event. */
    071        private double pendingEventTime;
    072
    073        /** Occurrence time of the previous event. */
    074        private double previousEventTime;
    075
    076        /** Integration direction. */
    077        private boolean forward;
    078
    079        /** Variation direction around pending event.
    080         *  (this is considered with respect to the integration direction)
    081         */
    082        private boolean increasing;
    083
    084        /** Next action indicator. */
    085        private EventHandler.Action nextAction;
    086
    087        /** Root-finding algorithm to use to detect state events. */
    088        private final UnivariateSolver solver;
    089
    090        /** Simple constructor.
    091         * @param handler event handler
    092         * @param maxCheckInterval maximal time interval between switching
    093         * function checks (this interval prevents missing sign changes in
    094         * case the integration steps becomes very large)
    095         * @param convergence convergence threshold in the event time search
    096         * @param maxIterationCount upper limit of the iteration count in
    097         * the event time search
    098         * @param solver Root-finding algorithm to use to detect state events
    099         */
    100        public EventState(final EventHandler handler, final double maxCheckInterval,
    101                          final double convergence, final int maxIterationCount,
    102                          final UnivariateSolver 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         * @exception MaxCountExceededException if the interpolator throws one because
    152         * the number of functions evaluations is exceeded
    153         */
    154        public void reinitializeBegin(final StepInterpolator interpolator)
    155            throws MaxCountExceededException {
    156
    157            t0 = interpolator.getPreviousTime();
    158            interpolator.setInterpolatedTime(t0);
    159            g0 = handler.g(t0, interpolator.getInterpolatedState());
    160            if (g0 == 0) {
    161                // excerpt from MATH-421 issue:
    162                // If an ODE solver is setup with an EventHandler that return STOP
    163                // when the even is triggered, the integrator stops (which is exactly
    164                // the expected behavior). If however the user wants to restart the
    165                // solver from the final state reached at the event with the same
    166                // configuration (expecting the event to be triggered again at a
    167                // later time), then the integrator may fail to start. It can get stuck
    168                // at the previous event. The use case for the bug MATH-421 is fairly
    169                // general, so events occurring exactly at start in the first step should
    170                // be ignored.
    171
    172                // extremely rare case: there is a zero EXACTLY at interval start
    173                // we will use the sign slightly after step beginning to force ignoring this zero
    174                final double epsilon = FastMath.max(solver.getAbsoluteAccuracy(),
    175                                                    FastMath.abs(solver.getRelativeAccuracy() * t0));
    176                final double tStart = t0 + 0.5 * epsilon;
    177                interpolator.setInterpolatedTime(tStart);
    178                g0 = handler.g(tStart, interpolator.getInterpolatedState());
    179            }
    180            g0Positive = g0 >= 0;
    181
    182        }
    183
    184        /** Evaluate the impact of the proposed step on the event handler.
    185         * @param interpolator step interpolator for the proposed step
    186         * @return true if the event handler triggers an event before
    187         * the end of the proposed step
    188         * @exception MaxCountExceededException if the interpolator throws one because
    189         * the number of functions evaluations is exceeded
    190         * @exception NoBracketingException if the event cannot be bracketed
    191         */
    192        public boolean evaluateStep(final StepInterpolator interpolator)
    193            throws MaxCountExceededException, NoBracketingException {
    194
    195            try {
    196                forward = interpolator.isForward();
    197                final double t1 = interpolator.getCurrentTime();
    198                final double dt = t1 - t0;
    199                if (FastMath.abs(dt) < convergence) {
    200                    // we cannot do anything on such a small step, don't trigger any events
    201                    return false;
    202                }
    203                final int    n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / maxCheckInterval));
    204                final double h = dt / n;
    205
    206                final UnivariateFunction f = new UnivariateFunction() {
    207                    public double value(final double t) throws LocalMaxCountExceededException {
    208                        try {
    209                            interpolator.setInterpolatedTime(t);
    210                            return handler.g(t, interpolator.getInterpolatedState());
    211                        } catch (MaxCountExceededException mcee) {
    212                            throw new LocalMaxCountExceededException(mcee);
    213                        }
    214                    }
    215                };
    216
    217                double ta = t0;
    218                double ga = g0;
    219                for (int i = 0; i < n; ++i) {
    220
    221                    // evaluate handler value at the end of the substep
    222                    final double tb = t0 + (i + 1) * h;
    223                    interpolator.setInterpolatedTime(tb);
    224                    final double gb = handler.g(tb, interpolator.getInterpolatedState());
    225
    226                    // check events occurrence
    227                    if (g0Positive ^ (gb >= 0)) {
    228                        // there is a sign change: an event is expected during this step
    229
    230                        // variation direction, with respect to the integration direction
    231                        increasing = gb >= ga;
    232
    233                        // find the event time making sure we select a solution just at or past the exact root
    234                        final double root;
    235                        if (solver instanceof BracketedUnivariateSolver<?>) {
    236                            @SuppressWarnings("unchecked")
    237                            BracketedUnivariateSolver<UnivariateFunction> bracketing =
    238                                    (BracketedUnivariateSolver<UnivariateFunction>) solver;
    239                            root = forward ?
    240                                   bracketing.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
    241                                   bracketing.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
    242                        } else {
    243                            final double baseRoot = forward ?
    244                                                    solver.solve(maxIterationCount, f, ta, tb) :
    245                                                    solver.solve(maxIterationCount, f, tb, ta);
    246                            final int remainingEval = maxIterationCount - solver.getEvaluations();
    247                            BracketedUnivariateSolver<UnivariateFunction> bracketing =
    248                                    new PegasusSolver(solver.getRelativeAccuracy(), solver.getAbsoluteAccuracy());
    249                            root = forward ?
    250                                   UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
    251                                                                       baseRoot, ta, tb, AllowedSolution.RIGHT_SIDE) :
    252                                   UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
    253                                                                       baseRoot, tb, ta, AllowedSolution.LEFT_SIDE);
    254                        }
    255
    256                        if ((!Double.isNaN(previousEventTime)) &&
    257                            (FastMath.abs(root - ta) <= convergence) &&
    258                            (FastMath.abs(root - previousEventTime) <= convergence)) {
    259                            // we have either found nothing or found (again ?) a past event,
    260                            // retry the substep excluding this value
    261                            ta = forward ? ta + convergence : ta - convergence;
    262                            ga = f.value(ta);
    263                            --i;
    264                        } else if (Double.isNaN(previousEventTime) ||
    265                                   (FastMath.abs(previousEventTime - root) > convergence)) {
    266                            pendingEventTime = root;
    267                            pendingEvent = true;
    268                            return true;
    269                        } else {
    270                            // no sign change: there is no event for now
    271                            ta = tb;
    272                            ga = gb;
    273                        }
    274
    275                    } else {
    276                        // no sign change: there is no event for now
    277                        ta = tb;
    278                        ga = gb;
    279                    }
    280
    281                }
    282
    283                // no event during the whole step
    284                pendingEvent     = false;
    285                pendingEventTime = Double.NaN;
    286                return false;
    287
    288            } catch (LocalMaxCountExceededException lmcee) {
    289                throw lmcee.getException();
    290            }
    291
    292        }
    293
    294        /** Get the occurrence time of the event triggered in the current step.
    295         * @return occurrence time of the event triggered in the current
    296         * step or infinity if no events are triggered
    297         */
    298        public double getEventTime() {
    299            return pendingEvent ?
    300                   pendingEventTime :
    301                   (forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
    302        }
    303
    304        /** Acknowledge the fact the step has been accepted by the integrator.
    305         * @param t value of the independent <i>time</i> variable at the
    306         * end of the step
    307         * @param y array containing the current value of the state vector
    308         * at the end of the step
    309         */
    310        public void stepAccepted(final double t, final double[] y) {
    311
    312            t0 = t;
    313            g0 = handler.g(t, y);
    314
    315            if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) {
    316                // force the sign to its value "just after the event"
    317                previousEventTime = t;
    318                g0Positive        = increasing;
    319                nextAction        = handler.eventOccurred(t, y, !(increasing ^ forward));
    320            } else {
    321                g0Positive = g0 >= 0;
    322                nextAction = EventHandler.Action.CONTINUE;
    323            }
    324        }
    325
    326        /** Check if the integration should be stopped at the end of the
    327         * current step.
    328         * @return true if the integration should be stopped
    329         */
    330        public boolean stop() {
    331            return nextAction == EventHandler.Action.STOP;
    332        }
    333
    334        /** Let the event handler reset the state if it wants.
    335         * @param t value of the independent <i>time</i> variable at the
    336         * beginning of the next step
    337         * @param y array were to put the desired state vector at the beginning
    338         * of the next step
    339         * @return true if the integrator should reset the derivatives too
    340         */
    341        public boolean reset(final double t, final double[] y) {
    342
    343            if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) {
    344                return false;
    345            }
    346
    347            if (nextAction == EventHandler.Action.RESET_STATE) {
    348                handler.resetState(t, y);
    349            }
    350            pendingEvent      = false;
    351            pendingEventTime  = Double.NaN;
    352
    353            return (nextAction == EventHandler.Action.RESET_STATE) ||
    354                   (nextAction == EventHandler.Action.RESET_DERIVATIVES);
    355
    356        }
    357
    358        /** Local wrapper to propagate exceptions. */
    359        private static class LocalMaxCountExceededException extends RuntimeException {
    360
    361            /** Serializable UID. */
    362            private static final long serialVersionUID = 20120901L;
    363
    364            /** Wrapped exception. */
    365            private final MaxCountExceededException wrapped;
    366
    367            /** Simple constructor.
    368             * @param exception exception to wrap
    369             */
    370            public LocalMaxCountExceededException(final MaxCountExceededException exception) {
    371                wrapped = exception;
    372            }
    373
    374            /** Get the wrapped exception.
    375             * @return wrapped exception
    376             */
    377            public MaxCountExceededException getException() {
    378                return wrapped;
    379            }
    380
    381        }
    382
    383    }