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.legacy.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.legacy.analysis.solvers.BracketingNthOrderBrentSolver;
030import org.apache.commons.math4.legacy.analysis.solvers.UnivariateSolver;
031import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
032import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
033import org.apache.commons.math4.legacy.exception.NoBracketingException;
034import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
035import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
036import org.apache.commons.math4.legacy.ode.events.EventHandler;
037import org.apache.commons.math4.legacy.ode.events.EventState;
038import org.apache.commons.math4.legacy.ode.sampling.AbstractStepInterpolator;
039import org.apache.commons.math4.legacy.ode.sampling.StepHandler;
040import org.apache.commons.math4.core.jdkmath.JdkMath;
041import org.apache.commons.math4.legacy.core.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    /** Set the equations.
213     * @param equations equations to set
214     */
215    protected void setEquations(final ExpandableStatefulODE equations) {
216        this.expandable = equations;
217    }
218
219    /** Get the differential equations to integrate.
220     * @return differential equations to integrate
221     * @since 3.2
222     */
223    protected ExpandableStatefulODE getExpandable() {
224        return expandable;
225    }
226
227    /** Get the evaluations counter.
228     * @return evaluations counter
229     * @since 3.6
230     */
231    protected IntegerSequence.Incrementor getCounter() {
232        return evaluations;
233    }
234
235    /** {@inheritDoc} */
236    @Override
237    public double integrate(final FirstOrderDifferentialEquations equations,
238                            final double t0, final double[] y0, final double t, final double[] y)
239        throws DimensionMismatchException, NumberIsTooSmallException,
240               MaxCountExceededException, NoBracketingException {
241
242        if (y0.length != equations.getDimension()) {
243            throw new DimensionMismatchException(y0.length, equations.getDimension());
244        }
245        if (y.length != equations.getDimension()) {
246            throw new DimensionMismatchException(y.length, equations.getDimension());
247        }
248
249        // prepare expandable stateful equations
250        final ExpandableStatefulODE expandableODE = new ExpandableStatefulODE(equations);
251        expandableODE.setTime(t0);
252        expandableODE.setPrimaryState(y0);
253
254        // perform integration
255        integrate(expandableODE, t);
256
257        // extract results back from the stateful equations
258        System.arraycopy(expandableODE.getPrimaryState(), 0, y, 0, y.length);
259        return expandableODE.getTime();
260    }
261
262    /** Integrate a set of differential equations up to the given time.
263     * <p>This method solves an Initial Value Problem (IVP).</p>
264     * <p>The set of differential equations is composed of a main set, which
265     * can be extended by some sets of secondary equations. The set of
266     * equations must be already set up with initial time and partial states.
267     * At integration completion, the final time and partial states will be
268     * available in the same object.</p>
269     * <p>Since this method stores some internal state variables made
270     * available in its public interface during integration ({@link
271     * #getCurrentSignedStepsize()}), it is <em>not</em> thread-safe.</p>
272     * @param equations complete set of differential equations to integrate
273     * @param t target time for the integration
274     * (can be set to a value smaller than <code>t0</code> for backward integration)
275     * @exception NumberIsTooSmallException if integration step is too small
276     * @throws DimensionMismatchException if the dimension of the complete state does not
277     * match the complete equations sets dimension
278     * @exception MaxCountExceededException if the number of functions evaluations is exceeded
279     * @exception NoBracketingException if the location of an event cannot be bracketed
280     */
281    public abstract void integrate(ExpandableStatefulODE equations, double t)
282        throws NumberIsTooSmallException, DimensionMismatchException,
283               MaxCountExceededException, NoBracketingException;
284
285    /** Compute the derivatives and check the number of evaluations.
286     * @param t current value of the independent <I>time</I> variable
287     * @param y array containing the current value of the state vector
288     * @param yDot placeholder array where to put the time derivative of the state vector
289     * @exception MaxCountExceededException if the number of functions evaluations is exceeded
290     * @exception DimensionMismatchException if arrays dimensions do not match equations settings
291     * @exception NullPointerException if the ODE equations have not been set (i.e. if this method
292     * is called outside of a call to {@link #integrate(ExpandableStatefulODE, double)} or {@link
293     * #integrate(FirstOrderDifferentialEquations, double, double[], double, double[])})
294     */
295    public void computeDerivatives(final double t, final double[] y, final double[] yDot)
296        throws MaxCountExceededException, DimensionMismatchException, NullPointerException {
297        evaluations.increment();
298        expandable.computeDerivatives(t, y, yDot);
299    }
300
301    /** Set the stateInitialized flag.
302     * <p>This method must be called by integrators with the value
303     * {@code false} before they start integration, so a proper lazy
304     * initialization is done automatically on the first step.</p>
305     * @param stateInitialized new value for the flag
306     * @since 2.2
307     */
308    protected void setStateInitialized(final boolean stateInitialized) {
309        this.statesInitialized = stateInitialized;
310    }
311
312    /** Accept a step, triggering events and step handlers.
313     * @param interpolator step interpolator
314     * @param y state vector at step end time, must be reset if an event
315     * asks for resetting or if an events stops integration during the step
316     * @param yDot placeholder array where to put the time derivative of the state vector
317     * @param tEnd final integration time
318     * @return time at end of step
319     * @exception MaxCountExceededException if the interpolator throws one because
320     * the number of functions evaluations is exceeded
321     * @exception NoBracketingException if the location of an event cannot be bracketed
322     * @exception DimensionMismatchException if arrays dimensions do not match equations settings
323     * @since 2.2
324     */
325    protected double acceptStep(final AbstractStepInterpolator interpolator,
326                                final double[] y, final double[] yDot, final double tEnd)
327        throws MaxCountExceededException, DimensionMismatchException, NoBracketingException {
328
329            double previousT = interpolator.getGlobalPreviousTime();
330            final double currentT = interpolator.getGlobalCurrentTime();
331
332            // initialize the events states if needed
333            if (! statesInitialized) {
334                for (EventState state : eventsStates) {
335                    state.reinitializeBegin(interpolator);
336                }
337                statesInitialized = true;
338            }
339
340            // search for next events that may occur during the step
341            final int orderingSign = interpolator.isForward() ? +1 : -1;
342            SortedSet<EventState> occurringEvents = new TreeSet<>(new Comparator<EventState>() {
343
344                /** {@inheritDoc} */
345                @Override
346                public int compare(EventState es0, EventState es1) {
347                    return orderingSign * Double.compare(es0.getEventTime(), es1.getEventTime());
348                }
349            });
350
351            for (final EventState state : eventsStates) {
352                if (state.evaluateStep(interpolator)) {
353                    // the event occurs during the current step
354                    occurringEvents.add(state);
355                }
356            }
357
358            while (!occurringEvents.isEmpty()) {
359
360                // handle the chronologically first event
361                final Iterator<EventState> iterator = occurringEvents.iterator();
362                final EventState currentEvent = iterator.next();
363                iterator.remove();
364
365                // restrict the interpolator to the first part of the step, up to the event
366                final double eventT = currentEvent.getEventTime();
367                interpolator.setSoftPreviousTime(previousT);
368                interpolator.setSoftCurrentTime(eventT);
369
370                // get state at event time
371                interpolator.setInterpolatedTime(eventT);
372                final double[] eventYComplete = new double[y.length];
373                expandable.getPrimaryMapper().insertEquationData(interpolator.getInterpolatedState(),
374                                                                 eventYComplete);
375                int index = 0;
376                for (EquationsMapper secondary : expandable.getSecondaryMappers()) {
377                    secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index++),
378                                                 eventYComplete);
379                }
380
381                // advance all event states to current time
382                for (final EventState state : eventsStates) {
383                    state.stepAccepted(eventT, eventYComplete);
384                    isLastStep = isLastStep || state.stop();
385                }
386
387                // handle the first part of the step, up to the event
388                for (final StepHandler handler : stepHandlers) {
389                    handler.handleStep(interpolator, isLastStep);
390                }
391
392                if (isLastStep) {
393                    // the event asked to stop integration
394                    System.arraycopy(eventYComplete, 0, y, 0, y.length);
395                    return eventT;
396                }
397
398                resetOccurred = false;
399                final boolean needReset = currentEvent.reset(eventT, eventYComplete);
400                if (needReset) {
401                    // some event handler has triggered changes that
402                    // invalidate the derivatives, we need to recompute them
403                    interpolator.setInterpolatedTime(eventT);
404                    System.arraycopy(eventYComplete, 0, y, 0, y.length);
405                    computeDerivatives(eventT, y, yDot);
406                    resetOccurred = true;
407                    return eventT;
408                }
409
410                // prepare handling of the remaining part of the step
411                previousT = eventT;
412                interpolator.setSoftPreviousTime(eventT);
413                interpolator.setSoftCurrentTime(currentT);
414
415                // check if the same event occurs again in the remaining part of the step
416                if (currentEvent.evaluateStep(interpolator)) {
417                    // the event occurs during the current step
418                    occurringEvents.add(currentEvent);
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    /** Check the integration span.
447     * @param equations set of differential equations
448     * @param t target time for the integration
449     * @exception NumberIsTooSmallException if integration span is too small
450     * @exception DimensionMismatchException if adaptive step size integrators
451     * tolerance arrays dimensions are not compatible with equations settings
452     */
453    protected void sanityChecks(final ExpandableStatefulODE equations, final double t)
454        throws NumberIsTooSmallException, DimensionMismatchException {
455
456        final double threshold = 1000 * JdkMath.ulp(JdkMath.max(JdkMath.abs(equations.getTime()),
457                                                                  JdkMath.abs(t)));
458        final double dt = JdkMath.abs(equations.getTime() - t);
459        if (dt <= threshold) {
460            throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
461                                                dt, threshold, false);
462        }
463    }
464}