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;
019    
020    import java.util.ArrayList;
021    import java.util.Collection;
022    import java.util.Collections;
023    import java.util.Comparator;
024    import java.util.Iterator;
025    import java.util.List;
026    import java.util.SortedSet;
027    import java.util.TreeSet;
028    
029    import org.apache.commons.math3.analysis.solvers.BracketingNthOrderBrentSolver;
030    import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
031    import org.apache.commons.math3.exception.DimensionMismatchException;
032    import org.apache.commons.math3.exception.MaxCountExceededException;
033    import org.apache.commons.math3.exception.NoBracketingException;
034    import org.apache.commons.math3.exception.NumberIsTooSmallException;
035    import org.apache.commons.math3.exception.util.LocalizedFormats;
036    import org.apache.commons.math3.ode.events.EventHandler;
037    import org.apache.commons.math3.ode.events.EventState;
038    import org.apache.commons.math3.ode.sampling.AbstractStepInterpolator;
039    import org.apache.commons.math3.ode.sampling.StepHandler;
040    import org.apache.commons.math3.util.FastMath;
041    import org.apache.commons.math3.util.Incrementor;
042    import org.apache.commons.math3.util.Precision;
043    
044    /**
045     * Base class managing common boilerplate for all integrators.
046     * @version $Id: AbstractIntegrator.java 1416643 2012-12-03 19:37:14Z tn $
047     * @since 2.0
048     */
049    public abstract class AbstractIntegrator implements FirstOrderIntegrator {
050    
051        /** Step handler. */
052        protected Collection<StepHandler> stepHandlers;
053    
054        /** Current step start time. */
055        protected double stepStart;
056    
057        /** Current stepsize. */
058        protected double stepSize;
059    
060        /** Indicator for last step. */
061        protected boolean isLastStep;
062    
063        /** Indicator that a state or derivative reset was triggered by some event. */
064        protected boolean resetOccurred;
065    
066        /** Events states. */
067        private Collection<EventState> eventsStates;
068    
069        /** Initialization indicator of events states. */
070        private boolean statesInitialized;
071    
072        /** Name of the method. */
073        private final String name;
074    
075        /** Counter for number of evaluations. */
076        private Incrementor evaluations;
077    
078        /** Differential equations to integrate. */
079        private transient ExpandableStatefulODE expandable;
080    
081        /** Build an instance.
082         * @param name name of the method
083         */
084        public AbstractIntegrator(final String name) {
085            this.name = name;
086            stepHandlers = new ArrayList<StepHandler>();
087            stepStart = Double.NaN;
088            stepSize  = Double.NaN;
089            eventsStates = new ArrayList<EventState>();
090            statesInitialized = false;
091            evaluations = new Incrementor();
092            setMaxEvaluations(-1);
093            evaluations.resetCount();
094        }
095    
096        /** Build an instance with a null name.
097         */
098        protected AbstractIntegrator() {
099            this(null);
100        }
101    
102        /** {@inheritDoc} */
103        public String getName() {
104            return name;
105        }
106    
107        /** {@inheritDoc} */
108        public void addStepHandler(final StepHandler handler) {
109            stepHandlers.add(handler);
110        }
111    
112        /** {@inheritDoc} */
113        public Collection<StepHandler> getStepHandlers() {
114            return Collections.unmodifiableCollection(stepHandlers);
115        }
116    
117        /** {@inheritDoc} */
118        public void clearStepHandlers() {
119            stepHandlers.clear();
120        }
121    
122        /** {@inheritDoc} */
123        public void addEventHandler(final EventHandler handler,
124                                    final double maxCheckInterval,
125                                    final double convergence,
126                                    final int maxIterationCount) {
127            addEventHandler(handler, maxCheckInterval, convergence,
128                            maxIterationCount,
129                            new BracketingNthOrderBrentSolver(convergence, 5));
130        }
131    
132        /** {@inheritDoc} */
133        public void addEventHandler(final EventHandler handler,
134                                    final double maxCheckInterval,
135                                    final double convergence,
136                                    final int maxIterationCount,
137                                    final UnivariateSolver solver) {
138            eventsStates.add(new EventState(handler, maxCheckInterval, convergence,
139                                            maxIterationCount, solver));
140        }
141    
142        /** {@inheritDoc} */
143        public Collection<EventHandler> getEventHandlers() {
144            final List<EventHandler> list = new ArrayList<EventHandler>();
145            for (EventState state : eventsStates) {
146                list.add(state.getEventHandler());
147            }
148            return Collections.unmodifiableCollection(list);
149        }
150    
151        /** {@inheritDoc} */
152        public void clearEventHandlers() {
153            eventsStates.clear();
154        }
155    
156        /** {@inheritDoc} */
157        public double getCurrentStepStart() {
158            return stepStart;
159        }
160    
161        /** {@inheritDoc} */
162        public double getCurrentSignedStepsize() {
163            return stepSize;
164        }
165    
166        /** {@inheritDoc} */
167        public void setMaxEvaluations(int maxEvaluations) {
168            evaluations.setMaximalCount((maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations);
169        }
170    
171        /** {@inheritDoc} */
172        public int getMaxEvaluations() {
173            return evaluations.getMaximalCount();
174        }
175    
176        /** {@inheritDoc} */
177        public int getEvaluations() {
178            return evaluations.getCount();
179        }
180    
181        /** Prepare the start of an integration.
182         * @param t0 start value of the independent <i>time</i> variable
183         * @param y0 array containing the start value of the state vector
184         * @param t target time for the integration
185         */
186        protected void initIntegration(final double t0, final double[] y0, final double t) {
187    
188            evaluations.resetCount();
189    
190            for (final EventState state : eventsStates) {
191                state.getEventHandler().init(t0, y0, t);
192            }
193    
194            for (StepHandler handler : stepHandlers) {
195                handler.init(t0, y0, t);
196            }
197    
198            setStateInitialized(false);
199    
200        }
201    
202        /** Set the equations.
203         * @param equations equations to set
204         */
205        protected void setEquations(final ExpandableStatefulODE equations) {
206            this.expandable = equations;
207        }
208    
209        /** {@inheritDoc} */
210        public double integrate(final FirstOrderDifferentialEquations equations,
211                                final double t0, final double[] y0, final double t, final double[] y)
212            throws DimensionMismatchException, NumberIsTooSmallException,
213                   MaxCountExceededException, NoBracketingException {
214    
215            if (y0.length != equations.getDimension()) {
216                throw new DimensionMismatchException(y0.length, equations.getDimension());
217            }
218            if (y.length != equations.getDimension()) {
219                throw new DimensionMismatchException(y.length, equations.getDimension());
220            }
221    
222            // prepare expandable stateful equations
223            final ExpandableStatefulODE expandableODE = new ExpandableStatefulODE(equations);
224            expandableODE.setTime(t0);
225            expandableODE.setPrimaryState(y0);
226    
227            // perform integration
228            integrate(expandableODE, t);
229    
230            // extract results back from the stateful equations
231            System.arraycopy(expandableODE.getPrimaryState(), 0, y, 0, y.length);
232            return expandableODE.getTime();
233    
234        }
235    
236        /** Integrate a set of differential equations up to the given time.
237         * <p>This method solves an Initial Value Problem (IVP).</p>
238         * <p>The set of differential equations is composed of a main set, which
239         * can be extended by some sets of secondary equations. The set of
240         * equations must be already set up with initial time and partial states.
241         * At integration completion, the final time and partial states will be
242         * available in the same object.</p>
243         * <p>Since this method stores some internal state variables made
244         * available in its public interface during integration ({@link
245         * #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p>
246         * @param equations complete set of differential equations to integrate
247         * @param t target time for the integration
248         * (can be set to a value smaller than <code>t0</code> for backward integration)
249         * @exception NumberIsTooSmallException if integration step is too small
250         * @throws DimensionMismatchException if the dimension of the complete state does not
251         * match the complete equations sets dimension
252         * @exception MaxCountExceededException if the number of functions evaluations is exceeded
253         * @exception NoBracketingException if the location of an event cannot be bracketed
254         */
255        public abstract void integrate(ExpandableStatefulODE equations, double t)
256            throws NumberIsTooSmallException, DimensionMismatchException,
257                   MaxCountExceededException, NoBracketingException;
258    
259        /** Compute the derivatives and check the number of evaluations.
260         * @param t current value of the independent <I>time</I> variable
261         * @param y array containing the current value of the state vector
262         * @param yDot placeholder array where to put the time derivative of the state vector
263         * @exception MaxCountExceededException if the number of functions evaluations is exceeded
264         * @exception DimensionMismatchException if arrays dimensions do not match equations settings
265         */
266        public void computeDerivatives(final double t, final double[] y, final double[] yDot)
267            throws MaxCountExceededException, DimensionMismatchException {
268            evaluations.incrementCount();
269            expandable.computeDerivatives(t, y, yDot);
270        }
271    
272        /** Set the stateInitialized flag.
273         * <p>This method must be called by integrators with the value
274         * {@code false} before they start integration, so a proper lazy
275         * initialization is done automatically on the first step.</p>
276         * @param stateInitialized new value for the flag
277         * @since 2.2
278         */
279        protected void setStateInitialized(final boolean stateInitialized) {
280            this.statesInitialized = stateInitialized;
281        }
282    
283        /** Accept a step, triggering events and step handlers.
284         * @param interpolator step interpolator
285         * @param y state vector at step end time, must be reset if an event
286         * asks for resetting or if an events stops integration during the step
287         * @param yDot placeholder array where to put the time derivative of the state vector
288         * @param tEnd final integration time
289         * @return time at end of step
290         * @exception MaxCountExceededException if the interpolator throws one because
291         * the number of functions evaluations is exceeded
292         * @exception NoBracketingException if the location of an event cannot be bracketed
293         * @exception DimensionMismatchException if arrays dimensions do not match equations settings
294         * @since 2.2
295         */
296        protected double acceptStep(final AbstractStepInterpolator interpolator,
297                                    final double[] y, final double[] yDot, final double tEnd)
298            throws MaxCountExceededException, DimensionMismatchException, NoBracketingException {
299    
300                double previousT = interpolator.getGlobalPreviousTime();
301                final double currentT = interpolator.getGlobalCurrentTime();
302    
303                // initialize the events states if needed
304                if (! statesInitialized) {
305                    for (EventState state : eventsStates) {
306                        state.reinitializeBegin(interpolator);
307                    }
308                    statesInitialized = true;
309                }
310    
311                // search for next events that may occur during the step
312                final int orderingSign = interpolator.isForward() ? +1 : -1;
313                SortedSet<EventState> occuringEvents = new TreeSet<EventState>(new Comparator<EventState>() {
314    
315                    /** {@inheritDoc} */
316                    public int compare(EventState es0, EventState es1) {
317                        return orderingSign * Double.compare(es0.getEventTime(), es1.getEventTime());
318                    }
319    
320                });
321    
322                for (final EventState state : eventsStates) {
323                    if (state.evaluateStep(interpolator)) {
324                        // the event occurs during the current step
325                        occuringEvents.add(state);
326                    }
327                }
328    
329                while (!occuringEvents.isEmpty()) {
330    
331                    // handle the chronologically first event
332                    final Iterator<EventState> iterator = occuringEvents.iterator();
333                    final EventState currentEvent = iterator.next();
334                    iterator.remove();
335    
336                    // restrict the interpolator to the first part of the step, up to the event
337                    final double eventT = currentEvent.getEventTime();
338                    interpolator.setSoftPreviousTime(previousT);
339                    interpolator.setSoftCurrentTime(eventT);
340    
341                    // trigger the event
342                    interpolator.setInterpolatedTime(eventT);
343                    final double[] eventY = interpolator.getInterpolatedState().clone();
344                    currentEvent.stepAccepted(eventT, eventY);
345                    isLastStep = currentEvent.stop();
346    
347                    // handle the first part of the step, up to the event
348                    for (final StepHandler handler : stepHandlers) {
349                        handler.handleStep(interpolator, isLastStep);
350                    }
351    
352                    if (isLastStep) {
353                        // the event asked to stop integration
354                        System.arraycopy(eventY, 0, y, 0, y.length);
355                        for (final EventState remaining : occuringEvents) {
356                            remaining.stepAccepted(eventT, eventY);
357                        }
358                        return eventT;
359                    }
360    
361                    if (currentEvent.reset(eventT, eventY)) {
362                        // some event handler has triggered changes that
363                        // invalidate the derivatives, we need to recompute them
364                        System.arraycopy(eventY, 0, y, 0, y.length);
365                        computeDerivatives(eventT, y, yDot);
366                        resetOccurred = true;
367                        for (final EventState remaining : occuringEvents) {
368                            remaining.stepAccepted(eventT, eventY);
369                        }
370                        return eventT;
371                    }
372    
373                    // prepare handling of the remaining part of the step
374                    previousT = eventT;
375                    interpolator.setSoftPreviousTime(eventT);
376                    interpolator.setSoftCurrentTime(currentT);
377    
378                    // check if the same event occurs again in the remaining part of the step
379                    if (currentEvent.evaluateStep(interpolator)) {
380                        // the event occurs during the current step
381                        occuringEvents.add(currentEvent);
382                    }
383    
384                }
385    
386                interpolator.setInterpolatedTime(currentT);
387                final double[] currentY = interpolator.getInterpolatedState();
388                for (final EventState state : eventsStates) {
389                    state.stepAccepted(currentT, currentY);
390                    isLastStep = isLastStep || state.stop();
391                }
392                isLastStep = isLastStep || Precision.equals(currentT, tEnd, 1);
393    
394                // handle the remaining part of the step, after all events if any
395                for (StepHandler handler : stepHandlers) {
396                    handler.handleStep(interpolator, isLastStep);
397                }
398    
399                return currentT;
400    
401        }
402    
403        /** Check the integration span.
404         * @param equations set of differential equations
405         * @param t target time for the integration
406         * @exception NumberIsTooSmallException if integration span is too small
407         * @exception DimensionMismatchException if adaptive step size integrators
408         * tolerance arrays dimensions are not compatible with equations settings
409         */
410        protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
411            throws NumberIsTooSmallException, DimensionMismatchException {
412    
413            final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(equations.getTime()),
414                                                                      FastMath.abs(t)));
415            final double dt = FastMath.abs(equations.getTime() - t);
416            if (dt <= threshold) {
417                throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
418                                                    dt, threshold, false);
419            }
420    
421        }
422    
423    }