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