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
018package org.apache.commons.math3.ode.events;
019
020import org.apache.commons.math3.RealFieldElement;
021import org.apache.commons.math3.analysis.RealFieldUnivariateFunction;
022import org.apache.commons.math3.analysis.solvers.AllowedSolution;
023import org.apache.commons.math3.analysis.solvers.BracketedRealFieldUnivariateSolver;
024import org.apache.commons.math3.exception.MaxCountExceededException;
025import org.apache.commons.math3.exception.NoBracketingException;
026import org.apache.commons.math3.ode.FieldODEState;
027import org.apache.commons.math3.ode.FieldODEStateAndDerivative;
028import org.apache.commons.math3.ode.sampling.FieldStepInterpolator;
029import 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 * @param <T> the type of the field elements
042 * @since 3.6
043 */
044public class FieldEventState<T extends RealFieldElement<T>> {
045
046    /** Event handler. */
047    private final FieldEventHandler<T> 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 T 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 T t0;
060
061    /** Value of the events handler at the beginning of the step. */
062    private T 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 T pendingEventTime;
072
073    /** Occurrence time of the previous event. */
074    private T 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 Action nextAction;
086
087    /** Root-finding algorithm to use to detect state events. */
088    private final BracketedRealFieldUnivariateSolver<T> 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 FieldEventState(final FieldEventHandler<T> handler, final double maxCheckInterval,
101                           final T convergence, final int maxIterationCount,
102                           final BracketedRealFieldUnivariateSolver<T> solver) {
103        this.handler           = handler;
104        this.maxCheckInterval  = maxCheckInterval;
105        this.convergence       = convergence.abs();
106        this.maxIterationCount = maxIterationCount;
107        this.solver            = solver;
108
109        // some dummy values ...
110        t0                = null;
111        g0                = null;
112        g0Positive        = true;
113        pendingEvent      = false;
114        pendingEventTime  = null;
115        previousEventTime = null;
116        increasing        = true;
117        nextAction        = Action.CONTINUE;
118
119    }
120
121    /** Get the underlying event handler.
122     * @return underlying event handler
123     */
124    public FieldEventHandler<T> 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 T 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 FieldStepInterpolator<T> interpolator)
155        throws MaxCountExceededException {
156
157        final FieldODEStateAndDerivative<T> s0 = interpolator.getPreviousState();
158        t0 = s0.getTime();
159        g0 = handler.g(s0);
160        if (g0.getReal() == 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().getReal(),
175                                                FastMath.abs(solver.getRelativeAccuracy().multiply(t0).getReal()));
176            final T tStart = t0.add(0.5 * epsilon);
177            g0 = handler.g(interpolator.getInterpolatedState(tStart));
178        }
179        g0Positive = g0.getReal() >= 0;
180
181    }
182
183    /** Evaluate the impact of the proposed step on the event handler.
184     * @param interpolator step interpolator for the proposed step
185     * @return true if the event handler triggers an event before
186     * the end of the proposed step
187     * @exception MaxCountExceededException if the interpolator throws one because
188     * the number of functions evaluations is exceeded
189     * @exception NoBracketingException if the event cannot be bracketed
190     */
191    public boolean evaluateStep(final FieldStepInterpolator<T> interpolator)
192        throws MaxCountExceededException, NoBracketingException {
193
194        forward = interpolator.isForward();
195        final FieldODEStateAndDerivative<T> s1 = interpolator.getCurrentState();
196        final T t1 = s1.getTime();
197        final T dt = t1.subtract(t0);
198        if (dt.abs().subtract(convergence).getReal() < 0) {
199            // we cannot do anything on such a small step, don't trigger any events
200            return false;
201        }
202        final int n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt.getReal()) / maxCheckInterval));
203        final T   h = dt.divide(n);
204
205        final RealFieldUnivariateFunction<T> f = new RealFieldUnivariateFunction<T>() {
206            /** {@inheritDoc} */
207            public T value(final T t) {
208                return handler.g(interpolator.getInterpolatedState(t));
209            }
210        };
211
212        T ta = t0;
213        T ga = g0;
214        for (int i = 0; i < n; ++i) {
215
216            // evaluate handler value at the end of the substep
217            final T tb = (i == n - 1) ? t1 : t0.add(h.multiply(i + 1));
218            final T gb = handler.g(interpolator.getInterpolatedState(tb));
219
220            // check events occurrence
221            if (g0Positive ^ (gb.getReal() >= 0)) {
222                // there is a sign change: an event is expected during this step
223
224                // variation direction, with respect to the integration direction
225                increasing = gb.subtract(ga).getReal() >= 0;
226
227                // find the event time making sure we select a solution just at or past the exact root
228                final T root = forward ?
229                               solver.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
230                               solver.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
231
232                if (previousEventTime != null &&
233                    root.subtract(ta).abs().subtract(convergence).getReal() <= 0 &&
234                    root.subtract(previousEventTime).abs().subtract(convergence).getReal() <= 0) {
235                    // we have either found nothing or found (again ?) a past event,
236                    // retry the substep excluding this value, and taking care to have the
237                    // required sign in case the g function is noisy around its zero and
238                    // crosses the axis several times
239                    do {
240                        ta = forward ? ta.add(convergence) : ta.subtract(convergence);
241                        ga = f.value(ta);
242                    } while ((g0Positive ^ (ga.getReal() >= 0)) && (forward ^ (ta.subtract(tb).getReal() >= 0)));
243
244                    if (forward ^ (ta.subtract(tb).getReal() >= 0)) {
245                        // we were able to skip this spurious root
246                        --i;
247                    } else {
248                        // we can't avoid this root before the end of the step,
249                        // we have to handle it despite it is close to the former one
250                        // maybe we have two very close roots
251                        pendingEventTime = root;
252                        pendingEvent     = true;
253                        return true;
254                    }
255                } else if (previousEventTime == null ||
256                           previousEventTime.subtract(root).abs().subtract(convergence).getReal() > 0) {
257                    pendingEventTime = root;
258                    pendingEvent     = true;
259                    return true;
260                } else {
261                    // no sign change: there is no event for now
262                    ta = tb;
263                    ga = gb;
264                }
265
266            } else {
267                // no sign change: there is no event for now
268                ta = tb;
269                ga = gb;
270            }
271
272        }
273
274        // no event during the whole step
275        pendingEvent     = false;
276        pendingEventTime = null;
277        return false;
278
279    }
280
281    /** Get the occurrence time of the event triggered in the current step.
282     * @return occurrence time of the event triggered in the current
283     * step or infinity if no events are triggered
284     */
285    public T getEventTime() {
286        return pendingEvent ?
287               pendingEventTime :
288               t0.getField().getZero().add(forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
289    }
290
291    /** Acknowledge the fact the step has been accepted by the integrator.
292     * @param state state at the end of the step
293     */
294    public void stepAccepted(final FieldODEStateAndDerivative<T> state) {
295
296        t0 = state.getTime();
297        g0 = handler.g(state);
298
299        if (pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0) {
300            // force the sign to its value "just after the event"
301            previousEventTime = state.getTime();
302            g0Positive        = increasing;
303            nextAction        = handler.eventOccurred(state, !(increasing ^ forward));
304        } else {
305            g0Positive = g0.getReal() >= 0;
306            nextAction = Action.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 == Action.STOP;
316    }
317
318    /** Let the event handler reset the state if it wants.
319     * @param state state at the beginning of the next step
320     * @return reset state (may by the same as initial state if only
321     * derivatives should be reset), or null if nothing is reset
322     */
323    public FieldODEState<T> reset(final FieldODEStateAndDerivative<T> state) {
324
325        if (!(pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0)) {
326            return null;
327        }
328
329        final FieldODEState<T> newState;
330        if (nextAction == Action.RESET_STATE) {
331            newState = handler.resetState(state);
332        } else if (nextAction == Action.RESET_DERIVATIVES) {
333            newState = state;
334        } else {
335            newState = null;
336        }
337        pendingEvent      = false;
338        pendingEventTime  = null;
339
340        return newState;
341
342    }
343
344}