    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.math.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.math.ConvergenceException;
    030    import org.apache.commons.math.MaxEvaluationsExceededException;
    031    import org.apache.commons.math.exception.util.LocalizedFormats;
    032    import org.apache.commons.math.ode.events.CombinedEventsManager;
    033    import org.apache.commons.math.ode.events.EventException;
    034    import org.apache.commons.math.ode.events.EventHandler;
    035    import org.apache.commons.math.ode.events.EventState;
    036    import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
    037    import org.apache.commons.math.ode.sampling.StepHandler;
    038    import org.apache.commons.math.util.FastMath;
    039    import org.apache.commons.math.util.MathUtils;
    040
    041    /**
    042     * Base class managing common boilerplate for all integrators.
    043     * @version $Revision: 1073267 $ $Date: 2011-02-22 10:06:20 +0100 (mar. 22 f??vr. 2011) $
    044     * @since 2.0
    045     */
    046    public abstract class AbstractIntegrator implements FirstOrderIntegrator {
    047
    048        /** Step handler. */
    049        protected Collection<StepHandler> stepHandlers;
    050
    051        /** Current step start time. */
    052        protected double stepStart;
    053
    054        /** Current stepsize. */
    055        protected double stepSize;
    056
    057        /** Indicator for last step. */
    058        protected boolean isLastStep;
    059
    060        /** Indicator that a state or derivative reset was triggered by some event. */
    061        protected boolean resetOccurred;
    062
    063        /** Events states. */
    064        private Collection<EventState> eventsStates;
    065
    066        /** Initialization indicator of events states. */
    067        private boolean statesInitialized;
    068
    069        /** Name of the method. */
    070        private final String name;
    071
    072        /** Maximal number of evaluations allowed. */
    073        private int maxEvaluations;
    074
    075        /** Number of evaluations already performed. */
    076        private int evaluations;
    077
    078        /** Differential equations to integrate. */
    079        private transient FirstOrderDifferentialEquations equations;
    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            setMaxEvaluations(-1);
    092            resetEvaluations();
    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            eventsStates.add(new EventState(handler, maxCheckInterval, convergence, maxIterationCount));
    127        }
    128
    129        /** {@inheritDoc} */
    130        public Collection<EventHandler> getEventHandlers() {
    131            final List<EventHandler> list = new ArrayList<EventHandler>();
    132            for (EventState state : eventsStates) {
    133                list.add(state.getEventHandler());
    134            }
    135            return Collections.unmodifiableCollection(list);
    136        }
    137
    138        /** {@inheritDoc} */
    139        public void clearEventHandlers() {
    140            eventsStates.clear();
    141        }
    142
    143        /** Check if dense output is needed.
    144         * @return true if there is at least one event handler or if
    145         * one of the step handlers requires dense output
    146         */
    147        protected boolean requiresDenseOutput() {
    148            if (!eventsStates.isEmpty()) {
    149                return true;
    150            }
    151            for (StepHandler handler : stepHandlers) {
    152                if (handler.requiresDenseOutput()) {
    153                    return true;
    154                }
    155            }
    156            return false;
    157        }
    158
    159        /** {@inheritDoc} */
    160        public double getCurrentStepStart() {
    161            return stepStart;
    162        }
    163
    164        /** {@inheritDoc} */
    165        public double getCurrentSignedStepsize() {
    166            return stepSize;
    167        }
    168
    169        /** {@inheritDoc} */
    170        public void setMaxEvaluations(int maxEvaluations) {
    171            this.maxEvaluations = (maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations;
    172        }
    173
    174        /** {@inheritDoc} */
    175        public int getMaxEvaluations() {
    176            return maxEvaluations;
    177        }
    178
    179        /** {@inheritDoc} */
    180        public int getEvaluations() {
    181            return evaluations;
    182        }
    183
    184        /** Reset the number of evaluations to zero.
    185         */
    186        protected void resetEvaluations() {
    187            evaluations = 0;
    188        }
    189
    190        /** Set the differential equations.
    191         * @param equations differential equations to integrate
    192         * @see #computeDerivatives(double, double[], double[])
    193         */
    194        protected void setEquations(final FirstOrderDifferentialEquations equations) {
    195            this.equations = equations;
    196        }
    197
    198        /** Compute the derivatives and check the number of evaluations.
    199         * @param t current value of the independent <I>time</I> variable
    200         * @param y array containing the current value of the state vector
    201         * @param yDot placeholder array where to put the time derivative of the state vector
    202         * @throws DerivativeException this user-defined exception should be used if an error is
    203         * is triggered by user code
    204         */
    205        public void computeDerivatives(final double t, final double[] y, final double[] yDot)
    206            throws DerivativeException {
    207            if (++evaluations > maxEvaluations) {
    208                throw new DerivativeException(new MaxEvaluationsExceededException(maxEvaluations));
    209            }
    210            equations.computeDerivatives(t, y, yDot);
    211        }
    212
    213        /** Set the stateInitialized flag.
    214         * <p>This method must be called by integrators with the value
    215         * {@code false} before they start integration, so a proper lazy
    216         * initialization is done automatically on the first step.</p>
    217         * @param stateInitialized new value for the flag
    218         * @since 2.2
    219         */
    220        protected void setStateInitialized(final boolean stateInitialized) {
    221            this.statesInitialized = stateInitialized;
    222        }
    223
    224        /** Accept a step, triggering events and step handlers.
    225         * @param interpolator step interpolator
    226         * @param y state vector at step end time, must be reset if an event
    227         * asks for resetting or if an events stops integration during the step
    228         * @param yDot placeholder array where to put the time derivative of the state vector
    229         * @param tEnd final integration time
    230         * @return time at end of step
    231         * @throws DerivativeException this exception is propagated to the caller if
    232         * the underlying user function triggers one
    233         * @exception IntegratorException if the value of one event state cannot be evaluated
    234         * @since 2.2
    235         */
    236        protected double acceptStep(final AbstractStepInterpolator interpolator,
    237                                    final double[] y, final double[] yDot, final double tEnd)
    238            throws DerivativeException, IntegratorException {
    239
    240            try {
    241                double previousT = interpolator.getGlobalPreviousTime();
    242                final double currentT = interpolator.getGlobalCurrentTime();
    243                resetOccurred = false;
    244
    245                // initialize the events states if needed
    246                if (! statesInitialized) {
    247                    for (EventState state : eventsStates) {
    248                        state.reinitializeBegin(interpolator);
    249                    }
    250                    statesInitialized = true;
    251                }
    252
    253                // search for next events that may occur during the step
    254                final int orderingSign = interpolator.isForward() ? +1 : -1;
    255                SortedSet<EventState> occuringEvents = new TreeSet<EventState>(new Comparator<EventState>() {
    256
    257                    /** {@inheritDoc} */
    258                    public int compare(EventState es0, EventState es1) {
    259                        return orderingSign * Double.compare(es0.getEventTime(), es1.getEventTime());
    260                    }
    261
    262                });
    263
    264                for (final EventState state : eventsStates) {
    265                    if (state.evaluateStep(interpolator)) {
    266                        // the event occurs during the current step
    267                        occuringEvents.add(state);
    268                    }
    269                }
    270
    271                while (!occuringEvents.isEmpty()) {
    272
    273                    // handle the chronologically first event
    274                    final Iterator<EventState> iterator = occuringEvents.iterator();
    275                    final EventState currentEvent = iterator.next();
    276                    iterator.remove();
    277
    278                    // restrict the interpolator to the first part of the step, up to the event
    279                    final double eventT = currentEvent.getEventTime();
    280                    interpolator.setSoftPreviousTime(previousT);
    281                    interpolator.setSoftCurrentTime(eventT);
    282
    283                    // trigger the event
    284                    interpolator.setInterpolatedTime(eventT);
    285                    final double[] eventY = interpolator.getInterpolatedState();
    286                    currentEvent.stepAccepted(eventT, eventY);
    287                    isLastStep = currentEvent.stop();
    288
    289                    // handle the first part of the step, up to the event
    290                    for (final StepHandler handler : stepHandlers) {
    291                        handler.handleStep(interpolator, isLastStep);
    292                    }
    293
    294                    if (isLastStep) {
    295                        // the event asked to stop integration
    296                        System.arraycopy(eventY, 0, y, 0, y.length);
    297                        return eventT;
    298                    }
    299
    300                    if (currentEvent.reset(eventT, eventY)) {
    301                        // some event handler has triggered changes that
    302                        // invalidate the derivatives, we need to recompute them
    303                        System.arraycopy(eventY, 0, y, 0, y.length);
    304                        computeDerivatives(eventT, y, yDot);
    305                        resetOccurred = true;
    306                        return eventT;
    307                    }
    308
    309                    // prepare handling of the remaining part of the step
    310                    previousT = eventT;
    311                    interpolator.setSoftPreviousTime(eventT);
    312                    interpolator.setSoftCurrentTime(currentT);
    313
    314                    // check if the same event occurs again in the remaining part of the step
    315                    if (currentEvent.evaluateStep(interpolator)) {
    316                        // the event occurs during the current step
    317                        occuringEvents.add(currentEvent);
    318                    }
    319
    320                }
    321
    322                interpolator.setInterpolatedTime(currentT);
    323                final double[] currentY = interpolator.getInterpolatedState();
    324                for (final EventState state : eventsStates) {
    325                    state.stepAccepted(currentT, currentY);
    326                    isLastStep = isLastStep || state.stop();
    327                }
    328                isLastStep = isLastStep || MathUtils.equals(currentT, tEnd, 1);
    329
    330                // handle the remaining part of the step, after all events if any
    331                for (StepHandler handler : stepHandlers) {
    332                    handler.handleStep(interpolator, isLastStep);
    333                }
    334
    335                return currentT;
    336            } catch (EventException se) {
    337                final Throwable cause = se.getCause();
    338                if ((cause != null) && (cause instanceof DerivativeException)) {
    339                    throw (DerivativeException) cause;
    340                }
    341                throw new IntegratorException(se);
    342            } catch (ConvergenceException ce) {
    343                throw new IntegratorException(ce);
    344            }
    345
    346        }
    347
    348        /** Perform some sanity checks on the integration parameters.
    349         * @param ode differential equations set
    350         * @param t0 start time
    351         * @param y0 state vector at t0
    352         * @param t target time for the integration
    353         * @param y placeholder where to put the state vector
    354         * @exception IntegratorException if some inconsistency is detected
    355         */
    356        protected void sanityChecks(final FirstOrderDifferentialEquations ode,
    357                                    final double t0, final double[] y0,
    358                                    final double t, final double[] y)
    359            throws IntegratorException {
    360
    361            if (ode.getDimension() != y0.length) {
    362                throw new IntegratorException(
    363                        LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, ode.getDimension(), y0.length);
    364            }
    365
    366            if (ode.getDimension() != y.length) {
    367                throw new IntegratorException(
    368                        LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, ode.getDimension(), y.length);
    369            }
    370
    371            if (FastMath.abs(t - t0) <= 1.0e-12 * FastMath.max(FastMath.abs(t0), FastMath.abs(t))) {
    372                throw new IntegratorException(
    373                        LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
    374                        FastMath.abs(t - t0));
    375            }
    376
    377        }
    378
    379        /** Add an event handler for end time checking.
    380         * <p>This method can be used to simplify handling of integration end time.
    381         * It leverages the nominal stop condition with the exceptional stop
    382         * conditions.</p>
    383         * @param startTime integration start time
    384         * @param endTime desired end time
    385         * @param manager manager containing the user-defined handlers
    386         * @return a new manager containing all the user-defined handlers plus a
    387         * dedicated manager triggering a stop event at entTime
    388         * @deprecated as of 2.2, this method is not used any more
    389         */
    390        @Deprecated
    391        protected CombinedEventsManager addEndTimeChecker(final double startTime,
    392                                                          final double endTime,
    393                                                          final CombinedEventsManager manager) {
    394            CombinedEventsManager newManager = new CombinedEventsManager();
    395            for (final EventState state : manager.getEventsStates()) {
    396                newManager.addEventHandler(state.getEventHandler(),
    397                                           state.getMaxCheckInterval(),
    398                                           state.getConvergence(),
    399                                           state.getMaxIterationCount());
    400            }
    401            newManager.addEventHandler(new EndTimeChecker(endTime),
    402                                       Double.POSITIVE_INFINITY,
    403                                       FastMath.ulp(FastMath.max(FastMath.abs(startTime), FastMath.abs(endTime))),
    404                                       100);
    405            return newManager;
    406        }
    407
    408        /** Specialized event handler to stop integration.
    409         * @deprecated as of 2.2, this class is not used anymore
    410         */
    411        @Deprecated
    412        private static class EndTimeChecker implements EventHandler {
    413
    414            /** Desired end time. */
    415            private final double endTime;
    416
    417            /** Build an instance.
    418             * @param endTime desired time
    419             */
    420            public EndTimeChecker(final double endTime) {
    421                this.endTime = endTime;
    422            }
    423
    424            /** {@inheritDoc} */
    425            public int eventOccurred(double t, double[] y, boolean increasing) {
    426                return STOP;
    427            }
    428
    429            /** {@inheritDoc} */
    430            public double g(double t, double[] y) {
    431                return t - endTime;
    432            }
    433
    434            /** {@inheritDoc} */
    435            public void resetState(double t, double[] y) {
    436            }
    437
    438        }
    439
    440    }