    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    }