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.math3.ode.events;
19  
20  import org.apache.commons.math3.analysis.UnivariateFunction;
21  import org.apache.commons.math3.analysis.solvers.AllowedSolution;
22  import org.apache.commons.math3.analysis.solvers.BracketedUnivariateSolver;
23  import org.apache.commons.math3.analysis.solvers.PegasusSolver;
24  import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
25  import org.apache.commons.math3.analysis.solvers.UnivariateSolverUtils;
26  import org.apache.commons.math3.exception.MaxCountExceededException;
27  import org.apache.commons.math3.exception.NoBracketingException;
28  import org.apache.commons.math3.ode.EquationsMapper;
29  import org.apache.commons.math3.ode.ExpandableStatefulODE;
30  import org.apache.commons.math3.ode.sampling.StepInterpolator;
31  import org.apache.commons.math3.util.FastMath;
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       = FastMath.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 
126     /** Get the underlying event handler.
127      * @return underlying event handler
128      */
129     public EventHandler getEventHandler() {
130         return handler;
131     }
132 
133     /** Set the equation.
134      * @param expandable equation being integrated
135      */
136     public void setExpandable(final ExpandableStatefulODE expandable) {
137         this.expandable = expandable;
138     }
139 
140     /** Get the maximal time interval between events handler checks.
141      * @return maximal time interval between events handler checks
142      */
143     public double getMaxCheckInterval() {
144         return maxCheckInterval;
145     }
146 
147     /** Get the convergence threshold for event localization.
148      * @return convergence threshold for event localization
149      */
150     public double getConvergence() {
151         return convergence;
152     }
153 
154     /** Get the upper limit in the iteration count for event localization.
155      * @return upper limit in the iteration count for event localization
156      */
157     public int getMaxIterationCount() {
158         return maxIterationCount;
159     }
160 
161     /** Reinitialize the beginning of the step.
162      * @param interpolator valid for the current step
163      * @exception MaxCountExceededException if the interpolator throws one because
164      * the number of functions evaluations is exceeded
165      */
166     public void reinitializeBegin(final StepInterpolator interpolator)
167         throws MaxCountExceededException {
168 
169         t0 = interpolator.getPreviousTime();
170         interpolator.setInterpolatedTime(t0);
171         g0 = handler.g(t0, getCompleteState(interpolator));
172         if (g0 == 0) {
173             // excerpt from MATH-421 issue:
174             // If an ODE solver is setup with an EventHandler that return STOP
175             // when the even is triggered, the integrator stops (which is exactly
176             // the expected behavior). If however the user wants to restart the
177             // solver from the final state reached at the event with the same
178             // configuration (expecting the event to be triggered again at a
179             // later time), then the integrator may fail to start. It can get stuck
180             // at the previous event. The use case for the bug MATH-421 is fairly
181             // general, so events occurring exactly at start in the first step should
182             // be ignored.
183 
184             // extremely rare case: there is a zero EXACTLY at interval start
185             // we will use the sign slightly after step beginning to force ignoring this zero
186             final double epsilon = FastMath.max(solver.getAbsoluteAccuracy(),
187                                                 FastMath.abs(solver.getRelativeAccuracy() * t0));
188             final double tStart = t0 + 0.5 * epsilon;
189             interpolator.setInterpolatedTime(tStart);
190             g0 = handler.g(tStart, getCompleteState(interpolator));
191         }
192         g0Positive = g0 >= 0;
193 
194     }
195 
196     /** Get the complete state (primary and secondary).
197      * @param interpolator interpolator to use
198      * @return complete state
199      */
200     private double[] getCompleteState(final StepInterpolator interpolator) {
201 
202         final double[] complete = new double[expandable.getTotalDimension()];
203 
204         expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
205                                                          complete);
206         int index = 0;
207         for (EquationsMapper secondary : expandable.getSecondaryMappers()) {
208             secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index++),
209                                          complete);
210         }
211 
212         return complete;
213 
214     }
215 
216     /** Evaluate the impact of the proposed step on the event handler.
217      * @param interpolator step interpolator for the proposed step
218      * @return true if the event handler triggers an event before
219      * the end of the proposed step
220      * @exception MaxCountExceededException if the interpolator throws one because
221      * the number of functions evaluations is exceeded
222      * @exception NoBracketingException if the event cannot be bracketed
223      */
224     public boolean evaluateStep(final StepInterpolator interpolator)
225         throws MaxCountExceededException, NoBracketingException {
226 
227         try {
228             forward = interpolator.isForward();
229             final double t1 = interpolator.getCurrentTime();
230             final double dt = t1 - t0;
231             if (FastMath.abs(dt) < convergence) {
232                 // we cannot do anything on such a small step, don't trigger any events
233                 return false;
234             }
235             final int    n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / maxCheckInterval));
236             final double h = dt / n;
237 
238             final UnivariateFunction f = new UnivariateFunction() {
239                 public double value(final double t) throws LocalMaxCountExceededException {
240                     try {
241                         interpolator.setInterpolatedTime(t);
242                         return handler.g(t, getCompleteState(interpolator));
243                     } catch (MaxCountExceededException mcee) {
244                         throw new LocalMaxCountExceededException(mcee);
245                     }
246                 }
247             };
248 
249             double ta = t0;
250             double ga = g0;
251             for (int i = 0; i < n; ++i) {
252 
253                 // evaluate handler value at the end of the substep
254                 final double tb = t0 + (i + 1) * h;
255                 interpolator.setInterpolatedTime(tb);
256                 final double gb = handler.g(tb, getCompleteState(interpolator));
257 
258                 // check events occurrence
259                 if (g0Positive ^ (gb >= 0)) {
260                     // there is a sign change: an event is expected during this step
261 
262                     // variation direction, with respect to the integration direction
263                     increasing = gb >= ga;
264 
265                     // find the event time making sure we select a solution just at or past the exact root
266                     final double root;
267                     if (solver instanceof BracketedUnivariateSolver<?>) {
268                         @SuppressWarnings("unchecked")
269                         BracketedUnivariateSolver<UnivariateFunction> bracketing =
270                                 (BracketedUnivariateSolver<UnivariateFunction>) solver;
271                         root = forward ?
272                                bracketing.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
273                                bracketing.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
274                     } else {
275                         final double baseRoot = forward ?
276                                                 solver.solve(maxIterationCount, f, ta, tb) :
277                                                 solver.solve(maxIterationCount, f, tb, ta);
278                         final int remainingEval = maxIterationCount - solver.getEvaluations();
279                         BracketedUnivariateSolver<UnivariateFunction> bracketing =
280                                 new PegasusSolver(solver.getRelativeAccuracy(), solver.getAbsoluteAccuracy());
281                         root = forward ?
282                                UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
283                                                                    baseRoot, ta, tb, AllowedSolution.RIGHT_SIDE) :
284                                UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
285                                                                    baseRoot, tb, ta, AllowedSolution.LEFT_SIDE);
286                     }
287 
288                     if ((!Double.isNaN(previousEventTime)) &&
289                         (FastMath.abs(root - ta) <= convergence) &&
290                         (FastMath.abs(root - previousEventTime) <= convergence)) {
291                         // we have either found nothing or found (again ?) a past event,
292                         // retry the substep excluding this value, and taking care to have the
293                         // required sign in case the g function is noisy around its zero and
294                         // crosses the axis several times
295                         do {
296                             ta = forward ? ta + convergence : ta - convergence;
297                             ga = f.value(ta);
298                         } while ((g0Positive ^ (ga >= 0)) && (forward ^ (ta >= tb)));
299                         --i;
300                     } else if (Double.isNaN(previousEventTime) ||
301                                (FastMath.abs(previousEventTime - root) > convergence)) {
302                         pendingEventTime = root;
303                         pendingEvent = true;
304                         return true;
305                     } else {
306                         // no sign change: there is no event for now
307                         ta = tb;
308                         ga = gb;
309                     }
310 
311                 } else {
312                     // no sign change: there is no event for now
313                     ta = tb;
314                     ga = gb;
315                 }
316 
317             }
318 
319             // no event during the whole step
320             pendingEvent     = false;
321             pendingEventTime = Double.NaN;
322             return false;
323 
324         } catch (LocalMaxCountExceededException lmcee) {
325             throw lmcee.getException();
326         }
327 
328     }
329 
330     /** Get the occurrence time of the event triggered in the current step.
331      * @return occurrence time of the event triggered in the current
332      * step or infinity if no events are triggered
333      */
334     public double getEventTime() {
335         return pendingEvent ?
336                pendingEventTime :
337                (forward ? Double.POSITIVE_INFINITY : Double.NEGATIVE_INFINITY);
338     }
339 
340     /** Acknowledge the fact the step has been accepted by the integrator.
341      * @param t value of the independent <i>time</i> variable at the
342      * end of the step
343      * @param y array containing the current value of the state vector
344      * at the end of the step
345      */
346     public void stepAccepted(final double t, final double[] y) {
347 
348         t0 = t;
349         g0 = handler.g(t, y);
350 
351         if (pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence)) {
352             // force the sign to its value "just after the event"
353             previousEventTime = t;
354             g0Positive        = increasing;
355             nextAction        = handler.eventOccurred(t, y, !(increasing ^ forward));
356         } else {
357             g0Positive = g0 >= 0;
358             nextAction = EventHandler.Action.CONTINUE;
359         }
360     }
361 
362     /** Check if the integration should be stopped at the end of the
363      * current step.
364      * @return true if the integration should be stopped
365      */
366     public boolean stop() {
367         return nextAction == EventHandler.Action.STOP;
368     }
369 
370     /** Let the event handler reset the state if it wants.
371      * @param t value of the independent <i>time</i> variable at the
372      * beginning of the next step
373      * @param y array were to put the desired state vector at the beginning
374      * of the next step
375      * @return true if the integrator should reset the derivatives too
376      */
377     public boolean reset(final double t, final double[] y) {
378 
379         if (!(pendingEvent && (FastMath.abs(pendingEventTime - t) <= convergence))) {
380             return false;
381         }
382 
383         if (nextAction == EventHandler.Action.RESET_STATE) {
384             handler.resetState(t, y);
385         }
386         pendingEvent      = false;
387         pendingEventTime  = Double.NaN;
388 
389         return (nextAction == EventHandler.Action.RESET_STATE) ||
390                (nextAction == EventHandler.Action.RESET_DERIVATIVES);
391 
392     }
393 
394     /** Local wrapper to propagate exceptions. */
395     private static class LocalMaxCountExceededException extends RuntimeException {
396 
397         /** Serializable UID. */
398         private static final long serialVersionUID = 20120901L;
399 
400         /** Wrapped exception. */
401         private final MaxCountExceededException wrapped;
402 
403         /** Simple constructor.
404          * @param exception exception to wrap
405          */
406         public LocalMaxCountExceededException(final MaxCountExceededException exception) {
407             wrapped = exception;
408         }
409 
410         /** Get the wrapped exception.
411          * @return wrapped exception
412          */
413         public MaxCountExceededException getException() {
414             return wrapped;
415         }
416 
417     }
418 
419 }