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.IntegerSequence;
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 IntegerSequence.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 = IntegerSequence.Incrementor.create().withMaximalCount(Integer.MAX_VALUE);
091    }
092
093    /** Build an instance with a null name.
094     */
095    protected AbstractIntegrator() {
096        this(null);
097    }
098
099    /** {@inheritDoc} */
100    public String getName() {
101        return name;
102    }
103
104    /** {@inheritDoc} */
105    public void addStepHandler(final StepHandler handler) {
106        stepHandlers.add(handler);
107    }
108
109    /** {@inheritDoc} */
110    public Collection<StepHandler> getStepHandlers() {
111        return Collections.unmodifiableCollection(stepHandlers);
112    }
113
114    /** {@inheritDoc} */
115    public void clearStepHandlers() {
116        stepHandlers.clear();
117    }
118
119    /** {@inheritDoc} */
120    public void addEventHandler(final EventHandler handler,
121                                final double maxCheckInterval,
122                                final double convergence,
123                                final int maxIterationCount) {
124        addEventHandler(handler, maxCheckInterval, convergence,
125                        maxIterationCount,
126                        new BracketingNthOrderBrentSolver(convergence, 5));
127    }
128
129    /** {@inheritDoc} */
130    public void addEventHandler(final EventHandler handler,
131                                final double maxCheckInterval,
132                                final double convergence,
133                                final int maxIterationCount,
134                                final UnivariateSolver solver) {
135        eventsStates.add(new EventState(handler, maxCheckInterval, convergence,
136                                        maxIterationCount, solver));
137    }
138
139    /** {@inheritDoc} */
140    public Collection<EventHandler> getEventHandlers() {
141        final List<EventHandler> list = new ArrayList<EventHandler>(eventsStates.size());
142        for (EventState state : eventsStates) {
143            list.add(state.getEventHandler());
144        }
145        return Collections.unmodifiableCollection(list);
146    }
147
148    /** {@inheritDoc} */
149    public void clearEventHandlers() {
150        eventsStates.clear();
151    }
152
153    /** {@inheritDoc} */
154    public double getCurrentStepStart() {
155        return stepStart;
156    }
157
158    /** {@inheritDoc} */
159    public double getCurrentSignedStepsize() {
160        return stepSize;
161    }
162
163    /** {@inheritDoc} */
164    public void setMaxEvaluations(int maxEvaluations) {
165        evaluations = evaluations.withMaximalCount((maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations);
166    }
167
168    /** {@inheritDoc} */
169    public int getMaxEvaluations() {
170        return evaluations.getMaximalCount();
171    }
172
173    /** {@inheritDoc} */
174    public int getEvaluations() {
175        return evaluations.getCount();
176    }
177
178    /** Prepare the start of an integration.
179     * @param t0 start value of the independent <i>time</i> variable
180     * @param y0 array containing the start value of the state vector
181     * @param t target time for the integration
182     */
183    protected void initIntegration(final double t0, final double[] y0, final double t) {
184
185        evaluations = evaluations.withStart(0);
186
187        for (final EventState state : eventsStates) {
188            state.setExpandable(expandable);
189            state.getEventHandler().init(t0, y0, t);
190        }
191
192        for (StepHandler handler : stepHandlers) {
193            handler.init(t0, y0, t);
194        }
195
196        setStateInitialized(false);
197
198    }
199
200    /** Set the equations.
201     * @param equations equations to set
202     */
203    protected void setEquations(final ExpandableStatefulODE equations) {
204        this.expandable = equations;
205    }
206
207    /** Get the differential equations to integrate.
208     * @return differential equations to integrate
209     * @since 3.2
210     */
211    protected ExpandableStatefulODE getExpandable() {
212        return expandable;
213    }
214
215    /** Get the evaluations counter.
216     * @return evaluations counter
217     * @since 3.2
218     * @deprecated as of 3.6 replaced with {@link #getCounter()}
219     */
220    @Deprecated
221    protected org.apache.commons.math3.util.Incrementor getEvaluationsCounter() {
222        return org.apache.commons.math3.util.Incrementor.wrap(evaluations);
223    }
224
225    /** Get the evaluations counter.
226     * @return evaluations counter
227     * @since 3.6
228     */
229    protected IntegerSequence.Incrementor getCounter() {
230        return evaluations;
231    }
232
233    /** {@inheritDoc} */
234    public double integrate(final FirstOrderDifferentialEquations equations,
235                            final double t0, final double[] y0, final double t, final double[] y)
236        throws DimensionMismatchException, NumberIsTooSmallException,
237               MaxCountExceededException, NoBracketingException {
238
239        if (y0.length != equations.getDimension()) {
240            throw new DimensionMismatchException(y0.length, equations.getDimension());
241        }
242        if (y.length != equations.getDimension()) {
243            throw new DimensionMismatchException(y.length, equations.getDimension());
244        }
245
246        // prepare expandable stateful equations
247        final ExpandableStatefulODE expandableODE = new ExpandableStatefulODE(equations);
248        expandableODE.setTime(t0);
249        expandableODE.setPrimaryState(y0);
250
251        // perform integration
252        integrate(expandableODE, t);
253
254        // extract results back from the stateful equations
255        System.arraycopy(expandableODE.getPrimaryState(), 0, y, 0, y.length);
256        return expandableODE.getTime();
257
258    }
259
260    /** Integrate a set of differential equations up to the given time.
261     * <p>This method solves an Initial Value Problem (IVP).</p>
262     * <p>The set of differential equations is composed of a main set, which
263     * can be extended by some sets of secondary equations. The set of
264     * equations must be already set up with initial time and partial states.
265     * At integration completion, the final time and partial states will be
266     * available in the same object.</p>
267     * <p>Since this method stores some internal state variables made
268     * available in its public interface during integration ({@link
269     * #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p>
270     * @param equations complete set of differential equations to integrate
271     * @param t target time for the integration
272     * (can be set to a value smaller than <code>t0</code> for backward integration)
273     * @exception NumberIsTooSmallException if integration step is too small
274     * @throws DimensionMismatchException if the dimension of the complete state does not
275     * match the complete equations sets dimension
276     * @exception MaxCountExceededException if the number of functions evaluations is exceeded
277     * @exception NoBracketingException if the location of an event cannot be bracketed
278     */
279    public abstract void integrate(ExpandableStatefulODE equations, double t)
280        throws NumberIsTooSmallException, DimensionMismatchException,
281               MaxCountExceededException, NoBracketingException;
282
283    /** Compute the derivatives and check the number of evaluations.
284     * @param t current value of the independent <I>time</I> variable
285     * @param y array containing the current value of the state vector
286     * @param yDot placeholder array where to put the time derivative of the state vector
287     * @exception MaxCountExceededException if the number of functions evaluations is exceeded
288     * @exception DimensionMismatchException if arrays dimensions do not match equations settings
289     * @exception NullPointerException if the ODE equations have not been set (i.e. if this method
290     * is called outside of a call to {@link #integrate(ExpandableStatefulODE, double)} or {@link
291     * #integrate(FirstOrderDifferentialEquations, double, double[], double, double[])})
292     */
293    public void computeDerivatives(final double t, final double[] y, final double[] yDot)
294        throws MaxCountExceededException, DimensionMismatchException, NullPointerException {
295        evaluations.increment();
296        expandable.computeDerivatives(t, y, yDot);
297    }
298
299    /** Set the stateInitialized flag.
300     * <p>This method must be called by integrators with the value
301     * {@code false} before they start integration, so a proper lazy
302     * initialization is done automatically on the first step.</p>
303     * @param stateInitialized new value for the flag
304     * @since 2.2
305     */
306    protected void setStateInitialized(final boolean stateInitialized) {
307        this.statesInitialized = stateInitialized;
308    }
309
310    /** Accept a step, triggering events and step handlers.
311     * @param interpolator step interpolator
312     * @param y state vector at step end time, must be reset if an event
313     * asks for resetting or if an events stops integration during the step
314     * @param yDot placeholder array where to put the time derivative of the state vector
315     * @param tEnd final integration time
316     * @return time at end of step
317     * @exception MaxCountExceededException if the interpolator throws one because
318     * the number of functions evaluations is exceeded
319     * @exception NoBracketingException if the location of an event cannot be bracketed
320     * @exception DimensionMismatchException if arrays dimensions do not match equations settings
321     * @since 2.2
322     */
323    protected double acceptStep(final AbstractStepInterpolator interpolator,
324                                final double[] y, final double[] yDot, final double tEnd)
325        throws MaxCountExceededException, DimensionMismatchException, NoBracketingException {
326
327            double previousT = interpolator.getGlobalPreviousTime();
328            final double currentT = interpolator.getGlobalCurrentTime();
329
330            // initialize the events states if needed
331            if (! statesInitialized) {
332                for (EventState state : eventsStates) {
333                    state.reinitializeBegin(interpolator);
334                }
335                statesInitialized = true;
336            }
337
338            // search for next events that may occur during the step
339            final int orderingSign = interpolator.isForward() ? +1 : -1;
340            SortedSet<EventState> occurringEvents = new TreeSet<EventState>(new Comparator<EventState>() {
341
342                /** {@inheritDoc} */
343                public int compare(EventState es0, EventState es1) {
344                    return orderingSign * Double.compare(es0.getEventTime(), es1.getEventTime());
345                }
346
347            });
348
349            for (final EventState state : eventsStates) {
350                if (state.evaluateStep(interpolator)) {
351                    // the event occurs during the current step
352                    occurringEvents.add(state);
353                }
354            }
355
356            while (!occurringEvents.isEmpty()) {
357
358                // handle the chronologically first event
359                final Iterator<EventState> iterator = occurringEvents.iterator();
360                final EventState currentEvent = iterator.next();
361                iterator.remove();
362
363                // restrict the interpolator to the first part of the step, up to the event
364                final double eventT = currentEvent.getEventTime();
365                interpolator.setSoftPreviousTime(previousT);
366                interpolator.setSoftCurrentTime(eventT);
367
368                // get state at event time
369                interpolator.setInterpolatedTime(eventT);
370                final double[] eventYComplete = new double[y.length];
371                expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
372                                                                 eventYComplete);
373                int index = 0;
374                for (EquationsMapper secondary : expandable.getSecondaryMappers()) {
375                    secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index++),
376                                                 eventYComplete);
377                }
378
379                // advance all event states to current time
380                for (final EventState state : eventsStates) {
381                    state.stepAccepted(eventT, eventYComplete);
382                    isLastStep = isLastStep || state.stop();
383                }
384
385                // handle the first part of the step, up to the event
386                for (final StepHandler handler : stepHandlers) {
387                    handler.handleStep(interpolator, isLastStep);
388                }
389
390                if (isLastStep) {
391                    // the event asked to stop integration
392                    System.arraycopy(eventYComplete, 0, y, 0, y.length);
393                    return eventT;
394                }
395
396                boolean needReset = false;
397                resetOccurred = false;
398                needReset = currentEvent.reset(eventT, eventYComplete);
399                if (needReset) {
400                    // some event handler has triggered changes that
401                    // invalidate the derivatives, we need to recompute them
402                    interpolator.setInterpolatedTime(eventT);
403                    System.arraycopy(eventYComplete, 0, y, 0, y.length);
404                    computeDerivatives(eventT, y, yDot);
405                    resetOccurred = true;
406                    return eventT;
407                }
408
409                // prepare handling of the remaining part of the step
410                previousT = eventT;
411                interpolator.setSoftPreviousTime(eventT);
412                interpolator.setSoftCurrentTime(currentT);
413
414                // check if the same event occurs again in the remaining part of the step
415                if (currentEvent.evaluateStep(interpolator)) {
416                    // the event occurs during the current step
417                    occurringEvents.add(currentEvent);
418                }
419
420            }
421
422            // last part of the step, after the last event
423            interpolator.setInterpolatedTime(currentT);
424            final double[] currentY = new double[y.length];
425            expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
426                                                             currentY);
427            int index = 0;
428            for (EquationsMapper secondary : expandable.getSecondaryMappers()) {
429                secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index++),
430                                             currentY);
431            }
432            for (final EventState state : eventsStates) {
433                state.stepAccepted(currentT, currentY);
434                isLastStep = isLastStep || state.stop();
435            }
436            isLastStep = isLastStep || Precision.equals(currentT, tEnd, 1);
437
438            // handle the remaining part of the step, after all events if any
439            for (StepHandler handler : stepHandlers) {
440                handler.handleStep(interpolator, isLastStep);
441            }
442
443            return currentT;
444
445    }
446
447    /** Check the integration span.
448     * @param equations set of differential equations
449     * @param t target time for the integration
450     * @exception NumberIsTooSmallException if integration span is too small
451     * @exception DimensionMismatchException if adaptive step size integrators
452     * tolerance arrays dimensions are not compatible with equations settings
453     */
454    protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
455        throws NumberIsTooSmallException, DimensionMismatchException {
456
457        final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(equations.getTime()),
458                                                                  FastMath.abs(t)));
459        final double dt = FastMath.abs(equations.getTime() - t);
460        if (dt <= threshold) {
461            throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
462                                                dt, threshold, false);
463        }
464
465    }
466
467}