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