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.analysis.UnivariateFunction;
21  import org.apache.commons.math4.legacy.analysis.solvers.AllowedSolution;
22  import org.apache.commons.math4.legacy.analysis.solvers.BracketedUnivariateSolver;
23  import org.apache.commons.math4.legacy.analysis.solvers.PegasusSolver;
24  import org.apache.commons.math4.legacy.analysis.solvers.UnivariateSolver;
25  import org.apache.commons.math4.legacy.analysis.solvers.UnivariateSolverUtils;
26  import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
27  import org.apache.commons.math4.legacy.exception.NoBracketingException;
28  import org.apache.commons.math4.legacy.ode.EquationsMapper;
29  import org.apache.commons.math4.legacy.ode.ExpandableStatefulODE;
30  import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
31  import org.apache.commons.math4.core.jdkmath.JdkMath;
32  
33  /** This class handles the state for one {@link EventHandler
34   * event handler} during integration steps.
35   *
36   * <p>Each time the integrator proposes a step, the event handler
37   * switching function should be checked. This class handles the state
38   * of one handler during one integration step, with references to the
39   * state at the end of the preceding step. This information is used to
40   * decide if the handler should trigger an event or not during the
41   * proposed step.</p>
42   *
43   * @since 1.2
44   */
45  public class EventState {
46  
47      /** Event handler. */
48      private final EventHandler handler;
49  
50      /** Maximal time interval between events handler checks. */
51      private final double maxCheckInterval;
52  
53      /** Convergence threshold for event localization. */
54      private final double convergence;
55  
56      /** Upper limit in the iteration count for event localization. */
57      private final int maxIterationCount;
58  
59      /** Equation being integrated. */
60      private ExpandableStatefulODE expandable;
61  
62      /** Time at the beginning of the step. */
63      private double t0;
64  
65      /** Value of the events handler at the beginning of the step. */
66      private double g0;
67  
68      /** Simulated sign of g0 (we cheat when crossing events). */
69      private boolean g0Positive;
70  
71      /** Indicator of event expected during the step. */
72      private boolean pendingEvent;
73  
74      /** Occurrence time of the pending event. */
75      private double pendingEventTime;
76  
77      /** Occurrence time of the previous event. */
78      private double previousEventTime;
79  
80      /** Integration direction. */
81      private boolean forward;
82  
83      /** Variation direction around pending event.
84       *  (this is considered with respect to the integration direction)
85       */
86      private boolean increasing;
87  
88      /** Next action indicator. */
89      private EventHandler.Action nextAction;
90  
91      /** Root-finding algorithm to use to detect state events. */
92      private final UnivariateSolver solver;
93  
94      /** Simple constructor.
95       * @param handler event handler
96       * @param maxCheckInterval maximal time interval between switching
97       * function checks (this interval prevents missing sign changes in
98       * case the integration steps becomes very large)
99       * @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 }