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