• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 févr. 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      */
AbstractIntegrator(final String name)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      */
AbstractIntegrator()97     protected AbstractIntegrator() {
98         this(null);
99     }
100 
101     /** {@inheritDoc} */
getName()102     public String getName() {
103         return name;
104     }
105 
106     /** {@inheritDoc} */
addStepHandler(final StepHandler handler)107     public void addStepHandler(final StepHandler handler) {
108         stepHandlers.add(handler);
109     }
110 
111     /** {@inheritDoc} */
getStepHandlers()112     public Collection<StepHandler> getStepHandlers() {
113         return Collections.unmodifiableCollection(stepHandlers);
114     }
115 
116     /** {@inheritDoc} */
clearStepHandlers()117     public void clearStepHandlers() {
118         stepHandlers.clear();
119     }
120 
121     /** {@inheritDoc} */
addEventHandler(final EventHandler handler, final double maxCheckInterval, final double convergence, final int maxIterationCount)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} */
getEventHandlers()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} */
clearEventHandlers()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      */
requiresDenseOutput()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} */
getCurrentStepStart()160     public double getCurrentStepStart() {
161         return stepStart;
162     }
163 
164     /** {@inheritDoc} */
getCurrentSignedStepsize()165     public double getCurrentSignedStepsize() {
166         return stepSize;
167     }
168 
169     /** {@inheritDoc} */
setMaxEvaluations(int maxEvaluations)170     public void setMaxEvaluations(int maxEvaluations) {
171         this.maxEvaluations = (maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations;
172     }
173 
174     /** {@inheritDoc} */
getMaxEvaluations()175     public int getMaxEvaluations() {
176         return maxEvaluations;
177     }
178 
179     /** {@inheritDoc} */
getEvaluations()180     public int getEvaluations() {
181         return evaluations;
182     }
183 
184     /** Reset the number of evaluations to zero.
185      */
resetEvaluations()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      */
setEquations(final FirstOrderDifferentialEquations equations)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      */
computeDerivatives(final double t, final double[] y, final double[] yDot)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      */
setStateInitialized(final boolean stateInitialized)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      */
acceptStep(final AbstractStepInterpolator interpolator, final double[] y, final double[] yDot, final double tEnd)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      */
sanityChecks(final FirstOrderDifferentialEquations ode, final double t0, final double[] y0, final double t, final double[] y)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
addEndTimeChecker(final double startTime, final double endTime, final CombinedEventsManager manager)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          */
EndTimeChecker(final double endTime)420         public EndTimeChecker(final double endTime) {
421             this.endTime = endTime;
422         }
423 
424         /** {@inheritDoc} */
eventOccurred(double t, double[] y, boolean increasing)425         public int eventOccurred(double t, double[] y, boolean increasing) {
426             return STOP;
427         }
428 
429         /** {@inheritDoc} */
g(double t, double[] y)430         public double g(double t, double[] y) {
431             return t - endTime;
432         }
433 
434         /** {@inheritDoc} */
resetState(double t, double[] y)435         public void resetState(double t, double[] y) {
436         }
437 
438     }
439 
440 }
441