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.math.ode; 19 20 import java.util.ArrayList; 21 import java.util.Collection; 22 import java.util.Collections; 23 import java.util.Comparator; 24 import java.util.Iterator; 25 import java.util.List; 26 import java.util.SortedSet; 27 import java.util.TreeSet; 28 29 import org.apache.commons.math.ConvergenceException; 30 import org.apache.commons.math.MaxEvaluationsExceededException; 31 import org.apache.commons.math.exception.util.LocalizedFormats; 32 import org.apache.commons.math.ode.events.CombinedEventsManager; 33 import org.apache.commons.math.ode.events.EventException; 34 import org.apache.commons.math.ode.events.EventHandler; 35 import org.apache.commons.math.ode.events.EventState; 36 import org.apache.commons.math.ode.sampling.AbstractStepInterpolator; 37 import org.apache.commons.math.ode.sampling.StepHandler; 38 import org.apache.commons.math.util.FastMath; 39 import org.apache.commons.math.util.MathUtils; 40 41 /** 42 * Base class managing common boilerplate for all integrators. 43 * @version $Revision: 1073267 $ $Date: 2011-02-22 10:06:20 +0100 (mar. 22 fvr. 2011) $ 44 * @since 2.0 45 */ 46 public abstract class AbstractIntegrator implements FirstOrderIntegrator { 47 48 /** Step handler. */ 49 protected Collection<StepHandler> stepHandlers; 50 51 /** Current step start time. */ 52 protected double stepStart; 53 54 /** Current stepsize. */ 55 protected double stepSize; 56 57 /** Indicator for last step. */ 58 protected boolean isLastStep; 59 60 /** Indicator that a state or derivative reset was triggered by some event. */ 61 protected boolean resetOccurred; 62 63 /** Events states. */ 64 private Collection<EventState> eventsStates; 65 66 /** Initialization indicator of events states. */ 67 private boolean statesInitialized; 68 69 /** Name of the method. */ 70 private final String name; 71 72 /** Maximal number of evaluations allowed. */ 73 private int maxEvaluations; 74 75 /** Number of evaluations already performed. */ 76 private int evaluations; 77 78 /** Differential equations to integrate. */ 79 private transient FirstOrderDifferentialEquations equations; 80 81 /** Build an instance. 82 * @param name name of the method 83 */ 84 public AbstractIntegrator(final String name) { 85 this.name = name; 86 stepHandlers = new ArrayList<StepHandler>(); 87 stepStart = Double.NaN; 88 stepSize = Double.NaN; 89 eventsStates = new ArrayList<EventState>(); 90 statesInitialized = false; 91 setMaxEvaluations(-1); 92 resetEvaluations(); 93 } 94 95 /** Build an instance with a null name. 96 */ 97 protected AbstractIntegrator() { 98 this(null); 99 } 100 101 /** {@inheritDoc} */ 102 public String getName() { 103 return name; 104 } 105 106 /** {@inheritDoc} */ 107 public void addStepHandler(final StepHandler handler) { 108 stepHandlers.add(handler); 109 } 110 111 /** {@inheritDoc} */ 112 public Collection<StepHandler> getStepHandlers() { 113 return Collections.unmodifiableCollection(stepHandlers); 114 } 115 116 /** {@inheritDoc} */ 117 public void clearStepHandlers() { 118 stepHandlers.clear(); 119 } 120 121 /** {@inheritDoc} */ 122 public void addEventHandler(final EventHandler handler, 123 final double maxCheckInterval, 124 final double convergence, 125 final int maxIterationCount) { 126 eventsStates.add(new EventState(handler, maxCheckInterval, convergence, maxIterationCount)); 127 } 128 129 /** {@inheritDoc} */ 130 public Collection<EventHandler> getEventHandlers() { 131 final List<EventHandler> list = new ArrayList<EventHandler>(); 132 for (EventState state : eventsStates) { 133 list.add(state.getEventHandler()); 134 } 135 return Collections.unmodifiableCollection(list); 136 } 137 138 /** {@inheritDoc} */ 139 public void clearEventHandlers() { 140 eventsStates.clear(); 141 } 142 143 /** Check if dense output is needed. 144 * @return true if there is at least one event handler or if 145 * one of the step handlers requires dense output 146 */ 147 protected boolean requiresDenseOutput() { 148 if (!eventsStates.isEmpty()) { 149 return true; 150 } 151 for (StepHandler handler : stepHandlers) { 152 if (handler.requiresDenseOutput()) { 153 return true; 154 } 155 } 156 return false; 157 } 158 159 /** {@inheritDoc} */ 160 public double getCurrentStepStart() { 161 return stepStart; 162 } 163 164 /** {@inheritDoc} */ 165 public double getCurrentSignedStepsize() { 166 return stepSize; 167 } 168 169 /** {@inheritDoc} */ 170 public void setMaxEvaluations(int maxEvaluations) { 171 this.maxEvaluations = (maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations; 172 } 173 174 /** {@inheritDoc} */ 175 public int getMaxEvaluations() { 176 return maxEvaluations; 177 } 178 179 /** {@inheritDoc} */ 180 public int getEvaluations() { 181 return evaluations; 182 } 183 184 /** Reset the number of evaluations to zero. 185 */ 186 protected void resetEvaluations() { 187 evaluations = 0; 188 } 189 190 /** Set the differential equations. 191 * @param equations differential equations to integrate 192 * @see #computeDerivatives(double, double[], double[]) 193 */ 194 protected void setEquations(final FirstOrderDifferentialEquations equations) { 195 this.equations = equations; 196 } 197 198 /** Compute the derivatives and check the number of evaluations. 199 * @param t current value of the independent <I>time</I> variable 200 * @param y array containing the current value of the state vector 201 * @param yDot placeholder array where to put the time derivative of the state vector 202 * @throws DerivativeException this user-defined exception should be used if an error is 203 * is triggered by user code 204 */ 205 public void computeDerivatives(final double t, final double[] y, final double[] yDot) 206 throws DerivativeException { 207 if (++evaluations > maxEvaluations) { 208 throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations)); 209 } 210 equations.computeDerivatives(t, y, yDot); 211 } 212 213 /** Set the stateInitialized flag. 214 * <p>This method must be called by integrators with the value 215 * {@code false} before they start integration, so a proper lazy 216 * initialization is done automatically on the first step.</p> 217 * @param stateInitialized new value for the flag 218 * @since 2.2 219 */ 220 protected void setStateInitialized(final boolean stateInitialized) { 221 this.statesInitialized = stateInitialized; 222 } 223 224 /** Accept a step, triggering events and step handlers. 225 * @param interpolator step interpolator 226 * @param y state vector at step end time, must be reset if an event 227 * asks for resetting or if an events stops integration during the step 228 * @param yDot placeholder array where to put the time derivative of the state vector 229 * @param tEnd final integration time 230 * @return time at end of step 231 * @throws DerivativeException this exception is propagated to the caller if 232 * the underlying user function triggers one 233 * @exception IntegratorException if the value of one event state cannot be evaluated 234 * @since 2.2 235 */ 236 protected double acceptStep(final AbstractStepInterpolator interpolator, 237 final double[] y, final double[] yDot, final double tEnd) 238 throws DerivativeException, IntegratorException { 239 240 try { 241 double previousT = interpolator.getGlobalPreviousTime(); 242 final double currentT = interpolator.getGlobalCurrentTime(); 243 resetOccurred = false; 244 245 // initialize the events states if needed 246 if (! statesInitialized) { 247 for (EventState state : eventsStates) { 248 state.reinitializeBegin(interpolator); 249 } 250 statesInitialized = true; 251 } 252 253 // search for next events that may occur during the step 254 final int orderingSign = interpolator.isForward() ? +1 : -1; 255 SortedSet<EventState> occuringEvents = new TreeSet<EventState>(new Comparator<EventState>() { 256 257 /** {@inheritDoc} */ 258 public int compare(EventState es0, EventState es1) { 259 return orderingSign * Double.compare(es0.getEventTime(), es1.getEventTime()); 260 } 261 262 }); 263 264 for (final EventState state : eventsStates) { 265 if (state.evaluateStep(interpolator)) { 266 // the event occurs during the current step 267 occuringEvents.add(state); 268 } 269 } 270 271 while (!occuringEvents.isEmpty()) { 272 273 // handle the chronologically first event 274 final Iterator<EventState> iterator = occuringEvents.iterator(); 275 final EventState currentEvent = iterator.next(); 276 iterator.remove(); 277 278 // restrict the interpolator to the first part of the step, up to the event 279 final double eventT = currentEvent.getEventTime(); 280 interpolator.setSoftPreviousTime(previousT); 281 interpolator.setSoftCurrentTime(eventT); 282 283 // trigger the event 284 interpolator.setInterpolatedTime(eventT); 285 final double[] eventY = interpolator.getInterpolatedState(); 286 currentEvent.stepAccepted(eventT, eventY); 287 isLastStep = currentEvent.stop(); 288 289 // handle the first part of the step, up to the event 290 for (final StepHandler handler : stepHandlers) { 291 handler.handleStep(interpolator, isLastStep); 292 } 293 294 if (isLastStep) { 295 // the event asked to stop integration 296 System.arraycopy(eventY, 0, y, 0, y.length); 297 return eventT; 298 } 299 300 if (currentEvent.reset(eventT, eventY)) { 301 // some event handler has triggered changes that 302 // invalidate the derivatives, we need to recompute them 303 System.arraycopy(eventY, 0, y, 0, y.length); 304 computeDerivatives(eventT, y, yDot); 305 resetOccurred = true; 306 return eventT; 307 } 308 309 // prepare handling of the remaining part of the step 310 previousT = eventT; 311 interpolator.setSoftPreviousTime(eventT); 312 interpolator.setSoftCurrentTime(currentT); 313 314 // check if the same event occurs again in the remaining part of the step 315 if (currentEvent.evaluateStep(interpolator)) { 316 // the event occurs during the current step 317 occuringEvents.add(currentEvent); 318 } 319 320 } 321 322 interpolator.setInterpolatedTime(currentT); 323 final double[] currentY = interpolator.getInterpolatedState(); 324 for (final EventState state : eventsStates) { 325 state.stepAccepted(currentT, currentY); 326 isLastStep = isLastStep || state.stop(); 327 } 328 isLastStep = isLastStep || MathUtils.equals(currentT, tEnd, 1); 329 330 // handle the remaining part of the step, after all events if any 331 for (StepHandler handler : stepHandlers) { 332 handler.handleStep(interpolator, isLastStep); 333 } 334 335 return currentT; 336 } catch (EventException se) { 337 final Throwable cause = se.getCause(); 338 if ((cause != null) && (cause instanceof DerivativeException)) { 339 throw (DerivativeException) cause; 340 } 341 throw new IntegratorException(se); 342 } catch (ConvergenceException ce) { 343 throw new IntegratorException(ce); 344 } 345 346 } 347 348 /** Perform some sanity checks on the integration parameters. 349 * @param ode differential equations set 350 * @param t0 start time 351 * @param y0 state vector at t0 352 * @param t target time for the integration 353 * @param y placeholder where to put the state vector 354 * @exception IntegratorException if some inconsistency is detected 355 */ 356 protected void sanityChecks(final FirstOrderDifferentialEquations ode, 357 final double t0, final double[] y0, 358 final double t, final double[] y) 359 throws IntegratorException { 360 361 if (ode.getDimension() != y0.length) { 362 throw new IntegratorException( 363 LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, ode.getDimension(), y0.length); 364 } 365 366 if (ode.getDimension() != y.length) { 367 throw new IntegratorException( 368 LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, ode.getDimension(), y.length); 369 } 370 371 if (FastMath.abs(t - t0) <= 1.0e-12 * FastMath.max(FastMath.abs(t0), FastMath.abs(t))) { 372 throw new IntegratorException( 373 LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL, 374 FastMath.abs(t - t0)); 375 } 376 377 } 378 379 /** Add an event handler for end time checking. 380 * <p>This method can be used to simplify handling of integration end time. 381 * It leverages the nominal stop condition with the exceptional stop 382 * conditions.</p> 383 * @param startTime integration start time 384 * @param endTime desired end time 385 * @param manager manager containing the user-defined handlers 386 * @return a new manager containing all the user-defined handlers plus a 387 * dedicated manager triggering a stop event at entTime 388 * @deprecated as of 2.2, this method is not used any more 389 */ 390 @Deprecated 391 protected CombinedEventsManager addEndTimeChecker(final double startTime, 392 final double endTime, 393 final CombinedEventsManager manager) { 394 CombinedEventsManager newManager = new CombinedEventsManager(); 395 for (final EventState state : manager.getEventsStates()) { 396 newManager.addEventHandler(state.getEventHandler(), 397 state.getMaxCheckInterval(), 398 state.getConvergence(), 399 state.getMaxIterationCount()); 400 } 401 newManager.addEventHandler(new EndTimeChecker(endTime), 402 Double.POSITIVE_INFINITY, 403 FastMath.ulp(FastMath.max(FastMath.abs(startTime), FastMath.abs(endTime))), 404 100); 405 return newManager; 406 } 407 408 /** Specialized event handler to stop integration. 409 * @deprecated as of 2.2, this class is not used anymore 410 */ 411 @Deprecated 412 private static class EndTimeChecker implements EventHandler { 413 414 /** Desired end time. */ 415 private final double endTime; 416 417 /** Build an instance. 418 * @param endTime desired time 419 */ 420 public EndTimeChecker(final double endTime) { 421 this.endTime = endTime; 422 } 423 424 /** {@inheritDoc} */ 425 public int eventOccurred(double t, double[] y, boolean increasing) { 426 return STOP; 427 } 428 429 /** {@inheritDoc} */ 430 public double g(double t, double[] y) { 431 return t - endTime; 432 } 433 434 /** {@inheritDoc} */ 435 public void resetState(double t, double[] y) { 436 } 437 438 } 439 440 } 441
