    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    }