View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math4.legacy.ode.events;
19  
20  import org.apache.commons.math4.legacy.core.RealFieldElement;
21  import org.apache.commons.math4.legacy.analysis.RealFieldUnivariateFunction;
22  import org.apache.commons.math4.legacy.analysis.solvers.AllowedSolution;
23  import org.apache.commons.math4.legacy.analysis.solvers.BracketedRealFieldUnivariateSolver;
24  import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
25  import org.apache.commons.math4.legacy.exception.NoBracketingException;
26  import org.apache.commons.math4.legacy.ode.FieldODEState;
27  import org.apache.commons.math4.legacy.ode.FieldODEStateAndDerivative;
28  import org.apache.commons.math4.legacy.ode.sampling.FieldStepInterpolator;
29  import org.apache.commons.math4.core.jdkmath.JdkMath;
30  
31  /** This class handles the state for one {@link EventHandler
32   * event handler} during integration steps.
33   *
34   * <p>Each time the integrator proposes a step, the event handler
35   * switching function should be checked. This class handles the state
36   * of one handler during one integration step, with references to the
37   * state at the end of the preceding step. This information is used to
38   * decide if the handler should trigger an event or not during the
39   * proposed step.</p>
40   *
41   * @param <T> the type of the field elements
42   * @since 3.6
43   */
44  public class FieldEventState<T extends RealFieldElement<T>> {
45  
46      /** Event handler. */
47      private final FieldEventHandler<T> handler;
48  
49      /** Maximal time interval between events handler checks. */
50      private final double maxCheckInterval;
51  
52      /** Convergence threshold for event localization. */
53      private final T convergence;
54  
55      /** Upper limit in the iteration count for event localization. */
56      private final int maxIterationCount;
57  
58      /** Time at the beginning of the step. */
59      private T t0;
60  
61      /** Value of the events handler at the beginning of the step. */
62      private T g0;
63  
64      /** Simulated sign of g0 (we cheat when crossing events). */
65      private boolean g0Positive;
66  
67      /** Indicator of event expected during the step. */
68      private boolean pendingEvent;
69  
70      /** Occurrence time of the pending event. */
71      private T pendingEventTime;
72  
73      /** Occurrence time of the previous event. */
74      private T previousEventTime;
75  
76      /** Integration direction. */
77      private boolean forward;
78  
79      /** Variation direction around pending event.
80       *  (this is considered with respect to the integration direction)
81       */
82      private boolean increasing;
83  
84      /** Next action indicator. */
85      private Action nextAction;
86  
87      /** Root-finding algorithm to use to detect state events. */
88      private final BracketedRealFieldUnivariateSolver<T> solver;
89  
90      /** Simple constructor.
91       * @param handler event handler
92       * @param maxCheckInterval maximal time interval between switching
93       * function checks (this interval prevents missing sign changes in
94       * case the integration steps becomes very large)
95       * @param convergence convergence threshold in the event time search
96       * @param maxIterationCount upper limit of the iteration count in
97       * the event time search
98       * @param solver Root-finding algorithm to use to detect state events
99       */
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     /** Get the underlying event handler.
121      * @return underlying event handler
122      */
123     public FieldEventHandler<T> getEventHandler() {
124         return handler;
125     }
126 
127     /** Get the maximal time interval between events handler checks.
128      * @return maximal time interval between events handler checks
129      */
130     public double getMaxCheckInterval() {
131         return maxCheckInterval;
132     }
133 
134     /** Get the convergence threshold for event localization.
135      * @return convergence threshold for event localization
136      */
137     public T getConvergence() {
138         return convergence;
139     }
140 
141     /** Get the upper limit in the iteration count for event localization.
142      * @return upper limit in the iteration count for event localization
143      */
144     public int getMaxIterationCount() {
145         return maxIterationCount;
146     }
147 
148     /** Reinitialize the beginning of the step.
149      * @param interpolator valid for the current step
150      * @exception MaxCountExceededException if the interpolator throws one because
151      * the number of functions evaluations is exceeded
152      */
153     public void reinitializeBegin(final FieldStepInterpolator<T> interpolator)
154         throws MaxCountExceededException {
155 
156         final FieldODEStateAndDerivative<T> s0 = interpolator.getPreviousState();
157         t0 = s0.getTime();
158         g0 = handler.g(s0);
159         if (g0.getReal() == 0) {
160             // excerpt from MATH-421 issue:
161             // If an ODE solver is setup with an EventHandler that return STOP
162             // when the even is triggered, the integrator stops (which is exactly
163             // the expected behavior). If however the user wants to restart the
164             // solver from the final state reached at the event with the same
165             // configuration (expecting the event to be triggered again at a
166             // later time), then the integrator may fail to start. It can get stuck
167             // at the previous event. The use case for the bug MATH-421 is fairly
168             // general, so events occurring exactly at start in the first step should
169             // be ignored.
170 
171             // extremely rare case: there is a zero EXACTLY at interval start
172             // we will use the sign slightly after step beginning to force ignoring this zero
173             final double epsilon = JdkMath.max(solver.getAbsoluteAccuracy().getReal(),
174                                                 JdkMath.abs(solver.getRelativeAccuracy().multiply(t0).getReal()));
175             final T tStart = t0.add(0.5 * epsilon);
176             g0 = handler.g(interpolator.getInterpolatedState(tStart));
177         }
178         g0Positive = g0.getReal() >= 0;
179     }
180 
181     /** Evaluate the impact of the proposed step on the event handler.
182      * @param interpolator step interpolator for the proposed step
183      * @return true if the event handler triggers an event before
184      * the end of the proposed step
185      * @exception MaxCountExceededException if the interpolator throws one because
186      * the number of functions evaluations is exceeded
187      * @exception NoBracketingException if the event cannot be bracketed
188      */
189     public boolean evaluateStep(final FieldStepInterpolator<T> interpolator)
190         throws MaxCountExceededException, NoBracketingException {
191 
192         forward = interpolator.isForward();
193         final FieldODEStateAndDerivative<T> s1 = interpolator.getCurrentState();
194         final T t1 = s1.getTime();
195         final T dt = t1.subtract(t0);
196         if (dt.abs().subtract(convergence).getReal() < 0) {
197             // we cannot do anything on such a small step, don't trigger any events
198             return false;
199         }
200         final int n = JdkMath.max(1, (int) JdkMath.ceil(JdkMath.abs(dt.getReal()) / maxCheckInterval));
201         final T   h = dt.divide(n);
202 
203         final RealFieldUnivariateFunction<T> f = new RealFieldUnivariateFunction<T>() {
204             /** {@inheritDoc} */
205             @Override
206             public T value(final T t) {
207                 return handler.g(interpolator.getInterpolatedState(t));
208             }
209         };
210 
211         T ta = t0;
212         T ga = g0;
213         for (int i = 0; i < n; ++i) {
214 
215             // evaluate handler value at the end of the substep
216             final T tb = (i == n - 1) ? t1 : t0.add(h.multiply(i + 1));
217             final T gb = handler.g(interpolator.getInterpolatedState(tb));
218 
219             // check events occurrence
220             if (g0Positive ^ (gb.getReal() >= 0)) {
221                 // there is a sign change: an event is expected during this step
222 
223                 // variation direction, with respect to the integration direction
224                 increasing = gb.subtract(ga).getReal() >= 0;
225 
226                 // find the event time making sure we select a solution just at or past the exact root
227                 final T root = forward ?
228                                solver.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
229                                solver.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
230 
231                 if (previousEventTime != null &&
232                     root.subtract(ta).abs().subtract(convergence).getReal() <= 0 &&
233                     root.subtract(previousEventTime).abs().subtract(convergence).getReal() <= 0) {
234                     // we have either found nothing or found (again ?) a past event,
235                     // retry the substep excluding this value, and taking care to have the
236                     // required sign in case the g function is noisy around its zero and
237                     // crosses the axis several times
238                     do {
239                         ta = forward ? ta.add(convergence) : ta.subtract(convergence);
240                         ga = f.value(ta);
241                     } while ((g0Positive ^ (ga.getReal() >= 0)) && (forward ^ (ta.subtract(tb).getReal() >= 0)));
242 
243                     if (forward ^ (ta.subtract(tb).getReal() >= 0)) {
244                         // we were able to skip this spurious root
245                         --i;
246                     } else {
247                         // we can't avoid this root before the end of the step,
248                         // we have to handle it despite it is close to the former one
249                         // maybe we have two very close roots
250                         pendingEventTime = root;
251                         pendingEvent     = true;
252                         return true;
253                     }
254                 } else if (previousEventTime == null ||
255                            previousEventTime.subtract(root).abs().subtract(convergence).getReal() > 0) {
256                     pendingEventTime = root;
257                     pendingEvent     = true;
258                     return true;
259                 } else {
260                     // no sign change: there is no event for now
261                     ta = tb;
262                     ga = gb;
263                 }
264             } else {
265                 // no sign change: there is no event for now
266                 ta = tb;
267                 ga = gb;
268             }
269         }
270 
271         // no event during the whole step
272         pendingEvent     = false;
273         pendingEventTime = null;
274         return false;
275     }
276 
277     /** Get the occurrence time of the event triggered in the current step.
278      * @return occurrence time of the event triggered in the current
279      * step or infinity if no events are triggered
280      */
281     public T getEventTime() {
282         return pendingEvent ?
283                pendingEventTime :
284                t0.getField().getZero().add(forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
285     }
286 
287     /** Acknowledge the fact the step has been accepted by the integrator.
288      * @param state state at the end of the step
289      */
290     public void stepAccepted(final FieldODEStateAndDerivative<T> state) {
291 
292         t0 = state.getTime();
293         g0 = handler.g(state);
294 
295         if (pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0) {
296             // force the sign to its value "just after the event"
297             previousEventTime = state.getTime();
298             g0Positive        = increasing;
299             nextAction        = handler.eventOccurred(state, increasing == forward);
300         } else {
301             g0Positive = g0.getReal() >= 0;
302             nextAction = Action.CONTINUE;
303         }
304     }
305 
306     /** Check if the integration should be stopped at the end of the
307      * current step.
308      * @return true if the integration should be stopped
309      */
310     public boolean stop() {
311         return nextAction == Action.STOP;
312     }
313 
314     /** Let the event handler reset the state if it wants.
315      * @param state state at the beginning of the next step
316      * @return reset state (may by the same as initial state if only
317      * derivatives should be reset), or null if nothing is reset
318      */
319     public FieldODEState<T> reset(final FieldODEStateAndDerivative<T> state) {
320 
321         if (!(pendingEvent && pendingEventTime.subtract(state.getTime()).abs().subtract(convergence).getReal() <= 0)) {
322             return null;
323         }
324 
325         final FieldODEState<T> newState;
326         if (nextAction == Action.RESET_STATE) {
327             newState = handler.resetState(state);
328         } else if (nextAction == Action.RESET_DERIVATIVES) {
329             newState = state;
330         } else {
331             newState = null;
332         }
333         pendingEvent      = false;
334         pendingEventTime  = null;
335 
336         return newState;
337     }
338 }