ode/events/EventState.java

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math.ode.events;

import org.apache.commons.math.ConvergenceException;
import org.apache.commons.math.FunctionEvaluationException;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.solvers.BrentSolver;
import org.apache.commons.math.exception.MathInternalError;
import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.sampling.StepInterpolator;
import org.apache.commons.math.util.FastMath;

/** This class handles the state for one {@link EventHandler
 * event handler} during integration steps.
 *
 * <p>Each time the integrator proposes a step, the event handler
 * switching function should be checked. This class handles the state
 * of one handler during one integration step, with references to the
 * state at the end of the preceding step. This information is used to
 * decide if the handler should trigger an event or not during the
 * proposed step (and hence the step should be reduced to ensure the
 * event occurs at a bound rather than inside the step).</p>
 *
 * @version $Revision: 1073158 $ $Date: 2011-02-21 22:46:52 +0100 (lun. 21 févr. 2011) $
 * @since 1.2
 */
public class EventState {

    /** Event handler. */
    private final EventHandler handler;

    /** Maximal time interval between events handler checks. */
    private final double maxCheckInterval;

    /** Convergence threshold for event localization. */
    private final double convergence;

    /** Upper limit in the iteration count for event localization. */
    private final int maxIterationCount;

    /** Time at the beginning of the step. */
    private double t0;

    /** Value of the events handler at the beginning of the step. */
    private double g0;

    /** Simulated sign of g0 (we cheat when crossing events). */
    private boolean g0Positive;

    /** Indicator of event expected during the step. */
    private boolean pendingEvent;

    /** Occurrence time of the pending event. */
    private double pendingEventTime;

    /** Occurrence time of the previous event. */
    private double previousEventTime;

    /** Integration direction. */
    private boolean forward;

    /** Variation direction around pending event.
     *  (this is considered with respect to the integration direction)
     */
    private boolean increasing;

    /** Next action indicator. */
    private int nextAction;

    /** Simple constructor.
     * @param handler event handler
     * @param maxCheckInterval maximal time interval between switching
     * function checks (this interval prevents missing sign changes in
     * case the integration steps becomes very large)
     * @param convergence convergence threshold in the event time search
     * @param maxIterationCount upper limit of the iteration count in
     * the event time search
     */
    public EventState(final EventHandler handler, final double maxCheckInterval,
                      final double convergence, final int maxIterationCount) {
        this.handler           = handler;
        this.maxCheckInterval  = maxCheckInterval;
        this.convergence       = FastMath.abs(convergence);
        this.maxIterationCount = maxIterationCount;

        // some dummy values ...
        t0                = Double.NaN;
        g0                = Double.NaN;
        g0Positive        = true;
        pendingEvent      = false;
        pendingEventTime  = Double.NaN;
        previousEventTime = Double.NaN;
        increasing        = true;
        nextAction        = EventHandler.CONTINUE;

    }

    /** Get the underlying event handler.
     * @return underlying event handler
     */
    public EventHandler getEventHandler() {
        return handler;
    }

    /** Get the maximal time interval between events handler checks.
     * @return maximal time interval between events handler checks
     */
    public double getMaxCheckInterval() {
        return maxCheckInterval;
    }

    /** Get the convergence threshold for event localization.
     * @return convergence threshold for event localization
     */
    public double getConvergence() {
        return convergence;
    }

    /** Get the upper limit in the iteration count for event localization.
     * @return upper limit in the iteration count for event localization
     */
    public int getMaxIterationCount() {
        return maxIterationCount;
    }

    /** Reinitialize the beginning of the step.
     * @param interpolator valid for the current step
     * @exception EventException if the event handler
     * value cannot be evaluated at the beginning of the step
     */
    public void reinitializeBegin(final StepInterpolator interpolator)
        throws EventException {
        try {
            // excerpt from MATH-421 issue:
            // If an ODE solver is setup with an EventHandler that return STOP
            // when the even is triggered, the integrator stops (which is exactly
            // the expected behavior). If however the user want to restart the
            // solver from the final state reached at the event with the same
            // configuration (expecting the event to be triggered again at a
            // later time), then the integrator may fail to start. It can get stuck
            // at the previous event.

            // The use case for the bug MATH-421 is fairly general, so events occurring
            // less than epsilon after the solver start in the first step should be ignored,
            // where epsilon is the convergence threshold of the event. The sign of the g
            // function should be evaluated after this initial ignore zone, not exactly at
            // beginning (if there are no event at the very beginning g(t0) and g(t0+epsilon)
            // have the same sign, so this does not hurt ; if there is an event at the very
            // beginning, g(t0) and g(t0+epsilon) have opposite signs and we want to start
            // with the second one. Of course, the sign of epsilon depend on the integration
            // direction (forward or backward). This explains what is done below.

            final double ignoreZone = interpolator.isForward() ? getConvergence() : -getConvergence();
            t0 = interpolator.getPreviousTime() + ignoreZone;
            interpolator.setInterpolatedTime(t0);
            g0 = handler.g(t0, interpolator.getInterpolatedState());
            if (g0 == 0) {
                // extremely rare case: there is a zero EXACTLY at end of ignore zone
                // we will use the opposite of sign at step beginning to force ignoring this zero
                final double tStart = interpolator.getPreviousTime();
                interpolator.setInterpolatedTime(tStart);
                g0Positive = handler.g(tStart, interpolator.getInterpolatedState()) <= 0;
            } else {
                g0Positive = g0 >= 0;
            }

        } catch (DerivativeException mue) {
            throw new EventException(mue);
        }
    }

    /** Evaluate the impact of the proposed step on the event handler.
     * @param interpolator step interpolator for the proposed step
     * @return true if the event handler triggers an event before
     * the end of the proposed step
     * @exception DerivativeException if the interpolator fails to
     * compute the switching function somewhere within the step
     * @exception EventException if the switching function
     * cannot be evaluated
     * @exception ConvergenceException if an event cannot be located
     */
    public boolean evaluateStep(final StepInterpolator interpolator)
        throws DerivativeException, EventException, ConvergenceException {

        try {

            forward = interpolator.isForward();
            final double t1 = interpolator.getCurrentTime();
            if (FastMath.abs(t1 - t0) < convergence) {
                // we cannot do anything on such a small step, don't trigger any events
                return false;
            }
            final double start = forward ? (t0 + convergence) : t0 - convergence;
            final double dt    = t1 - start;
            final int    n     = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / maxCheckInterval));
            final double h     = dt / n;

            double ta = t0;
            double ga = g0;
            for (int i = 0; i < n; ++i) {

                // evaluate handler value at the end of the substep
                final double tb = start + (i + 1) * h;
                interpolator.setInterpolatedTime(tb);
                final double gb = handler.g(tb, interpolator.getInterpolatedState());

                // check events occurrence
                if (g0Positive ^ (gb >= 0)) {
                    // there is a sign change: an event is expected during this step

                    // variation direction, with respect to the integration direction
                    increasing = gb >= ga;

                    final UnivariateRealFunction f = new UnivariateRealFunction() {
                        public double value(final double t) {
                            try {
                                interpolator.setInterpolatedTime(t);
                                return handler.g(t, interpolator.getInterpolatedState());
                            } catch (DerivativeException e) {
                                throw new EmbeddedDerivativeException(e);
                            } catch (EventException e) {
                                throw new EmbeddedEventException(e);
                            }
                        }
                    };
                    final BrentSolver solver = new BrentSolver(convergence);

                    if (ga * gb >= 0) {
                        // this is a corner case:
                        // - there was an event near ta,
                        // - there is another event between ta and tb
                        // - when ta was computed, convergence was reached on the "wrong side" of the interval
                        // this implies that the real sign of ga is the same as gb, so we need to slightly
                        // shift ta to make sure ga and gb get opposite signs and the solver won't complain
                        // about bracketing
                        final double epsilon = (forward ? 0.25 : -0.25) * convergence;
                        for (int k = 0; (k < 4) && (ga * gb > 0); ++k) {
                            ta += epsilon;
                            try {
                                ga = f.value(ta);
                            } catch (FunctionEvaluationException ex) {
                                throw new DerivativeException(ex);
                            }
                        }
                        if (ga * gb > 0) {
                            // this should never happen
                            throw new MathInternalError();
                        }
                    }

                    final double root;
                    try {
                        root = (ta <= tb) ?
                                solver.solve(maxIterationCount, f, ta, tb) :
                                    solver.solve(maxIterationCount, f, tb, ta);
                    } catch (FunctionEvaluationException ex) {
                        throw new DerivativeException(ex);
                    }

                    if ((!Double.isNaN(previousEventTime)) &&
                        (FastMath.abs(root - ta) <= convergence) &&
                        (FastMath.abs(root - previousEventTime) <= convergence)) {
                        // we have either found nothing or found (again ?) a past event, we simply ignore it
                        ta = tb;
                        ga = gb;
                    } else if (Double.isNaN(previousEventTime) ||
                               (FastMath.abs(previousEventTime - root) > convergence)) {
                        pendingEventTime = root;
                        pendingEvent = true;
                        return true;
                    } else {
                        // no sign change: there is no event for now
                        ta = tb;
                        ga = gb;
                    }

                } else {
                    // no sign change: there is no event for now
                    ta = tb;
                    ga = gb;
                }

            }

            // no event during the whole step
            pendingEvent     = false;
            pendingEventTime = Double.NaN;
            return false;

        } catch (EmbeddedDerivativeException ede) {
            throw ede.getDerivativeException();
        } catch (EmbeddedEventException eee) {
            throw eee.getEventException();
        }

    }

    /** Get the occurrence time of the event triggered in the current step.
     * @return occurrence time of the event triggered in the current
     * step or positive infinity if no events are triggered
     */
    public double getEventTime() {
        return pendingEvent ? pendingEventTime : Double.POSITIVE_INFINITY;
    }

    /** Acknowledge the fact the step has been accepted by the integrator.
     * @param t value of the independent <i>time</i> variable at the
     * end of the step
     * @param y array containing the current value of the state vector
     * at the end of the step
     * @exception EventException if the value of the event
     * handler cannot be evaluated
     */
    public void stepAccepted(final double t, final double[] y)
        throws EventException {

        t0 = t;
        g0 = handler.g(t, y);

        if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) {
            // force the sign to its value "just after the event"
            previousEventTime = t;
            g0Positive        = increasing;
            nextAction        = handler.eventOccurred(t, y, !(increasing ^ forward));
        } else {
            g0Positive = g0 >= 0;
            nextAction = EventHandler.CONTINUE;
        }
    }

    /** Check if the integration should be stopped at the end of the
     * current step.
     * @return true if the integration should be stopped
     */
    public boolean stop() {
        return nextAction == EventHandler.STOP;
    }

    /** Let the event handler reset the state if it wants.
     * @param t value of the independent <i>time</i> variable at the
     * beginning of the next step
     * @param y array were to put the desired state vector at the beginning
     * of the next step
     * @return true if the integrator should reset the derivatives too
     * @exception EventException if the state cannot be reseted by the event
     * handler
     */
    public boolean reset(final double t, final double[] y)
        throws EventException {

        if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) {
            return false;
        }

        if (nextAction == EventHandler.RESET_STATE) {
            handler.resetState(t, y);
        }
        pendingEvent      = false;
        pendingEventTime  = Double.NaN;

        return (nextAction == EventHandler.RESET_STATE) ||
               (nextAction == EventHandler.RESET_DERIVATIVES);

    }

    /** Local exception for embedding DerivativeException. */
    private static class EmbeddedDerivativeException extends RuntimeException {

        /** Serializable UID. */
        private static final long serialVersionUID = 3574188382434584610L;

        /** Embedded exception. */
        private final DerivativeException derivativeException;

        /** Simple constructor.
         * @param derivativeException embedded exception
         */
        public EmbeddedDerivativeException(final DerivativeException derivativeException) {
            this.derivativeException = derivativeException;
        }

        /** Get the embedded exception.
         * @return embedded exception
         */
        public DerivativeException getDerivativeException() {
            return derivativeException;
        }

    }

    /** Local exception for embedding EventException. */
    private static class EmbeddedEventException extends RuntimeException {

        /** Serializable UID. */
        private static final long serialVersionUID = -1337749250090455474L;

        /** Embedded exception. */
        private final EventException eventException;

        /** Simple constructor.
         * @param eventException embedded exception
         */
        public EmbeddedEventException(final EventException eventException) {
            this.eventException = eventException;
        }

        /** Get the embedded exception.
         * @return embedded exception
         */
        public EventException getEventException() {
            return eventException;
        }

    }
}