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