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