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.Field;
030import org.apache.commons.math3.RealFieldElement;
031import org.apache.commons.math3.analysis.solvers.BracketedRealFieldUnivariateSolver;
032import org.apache.commons.math3.analysis.solvers.FieldBracketingNthOrderBrentSolver;
033import org.apache.commons.math3.exception.DimensionMismatchException;
034import org.apache.commons.math3.exception.MaxCountExceededException;
035import org.apache.commons.math3.exception.NoBracketingException;
036import org.apache.commons.math3.exception.NumberIsTooSmallException;
037import org.apache.commons.math3.exception.util.LocalizedFormats;
038import org.apache.commons.math3.ode.events.FieldEventHandler;
039import org.apache.commons.math3.ode.events.FieldEventState;
040import org.apache.commons.math3.ode.sampling.AbstractFieldStepInterpolator;
041import org.apache.commons.math3.ode.sampling.FieldStepHandler;
042import org.apache.commons.math3.util.FastMath;
043import org.apache.commons.math3.util.IntegerSequence;
044
045/**
046 * Base class managing common boilerplate for all integrators.
047 * @param <T> the type of the field elements
048 * @since 3.6
049 */
050public abstract class AbstractFieldIntegrator<T extends RealFieldElement<T>> implements FirstOrderFieldIntegrator<T> {
051
052    /** Default relative accuracy. */
053    private static final double DEFAULT_RELATIVE_ACCURACY = 1e-14;
054
055    /** Default function value accuracy. */
056    private static final double DEFAULT_FUNCTION_VALUE_ACCURACY = 1e-15;
057
058    /** Step handler. */
059    private Collection<FieldStepHandler<T>> stepHandlers;
060
061    /** Current step start. */
062    private FieldODEStateAndDerivative<T> stepStart;
063
064    /** Current stepsize. */
065    private T stepSize;
066
067    /** Indicator for last step. */
068    private boolean isLastStep;
069
070    /** Indicator that a state or derivative reset was triggered by some event. */
071    private boolean resetOccurred;
072
073    /** Field to which the time and state vector elements belong. */
074    private final Field<T> field;
075
076    /** Events states. */
077    private Collection<FieldEventState<T>> eventsStates;
078
079    /** Initialization indicator of events states. */
080    private boolean statesInitialized;
081
082    /** Name of the method. */
083    private final String name;
084
085    /** Counter for number of evaluations. */
086    private IntegerSequence.Incrementor evaluations;
087
088    /** Differential equations to integrate. */
089    private transient FieldExpandableODE<T> equations;
090
091    /** Build an instance.
092     * @param field field to which the time and state vector elements belong
093     * @param name name of the method
094     */
095    protected AbstractFieldIntegrator(final Field<T> field, final String name) {
096        this.field        = field;
097        this.name         = name;
098        stepHandlers      = new ArrayList<FieldStepHandler<T>>();
099        stepStart         = null;
100        stepSize          = null;
101        eventsStates      = new ArrayList<FieldEventState<T>>();
102        statesInitialized = false;
103        evaluations       = IntegerSequence.Incrementor.create().withMaximalCount(Integer.MAX_VALUE);
104    }
105
106    /** Get the field to which state vector elements belong.
107     * @return field to which state vector elements belong
108     */
109    public Field<T> getField() {
110        return field;
111    }
112
113    /** {@inheritDoc} */
114    public String getName() {
115        return name;
116    }
117
118    /** {@inheritDoc} */
119    public void addStepHandler(final FieldStepHandler<T> handler) {
120        stepHandlers.add(handler);
121    }
122
123    /** {@inheritDoc} */
124    public Collection<FieldStepHandler<T>> getStepHandlers() {
125        return Collections.unmodifiableCollection(stepHandlers);
126    }
127
128    /** {@inheritDoc} */
129    public void clearStepHandlers() {
130        stepHandlers.clear();
131    }
132
133    /** {@inheritDoc} */
134    public void addEventHandler(final FieldEventHandler<T> handler,
135                                final double maxCheckInterval,
136                                final double convergence,
137                                final int maxIterationCount) {
138        addEventHandler(handler, maxCheckInterval, convergence,
139                        maxIterationCount,
140                        new FieldBracketingNthOrderBrentSolver<T>(field.getZero().add(DEFAULT_RELATIVE_ACCURACY),
141                                                                  field.getZero().add(convergence),
142                                                                  field.getZero().add(DEFAULT_FUNCTION_VALUE_ACCURACY),
143                                                                  5));
144    }
145
146    /** {@inheritDoc} */
147    public void addEventHandler(final FieldEventHandler<T> handler,
148                                final double maxCheckInterval,
149                                final double convergence,
150                                final int maxIterationCount,
151                                final BracketedRealFieldUnivariateSolver<T> solver) {
152        eventsStates.add(new FieldEventState<T>(handler, maxCheckInterval, field.getZero().add(convergence),
153                                                maxIterationCount, solver));
154    }
155
156    /** {@inheritDoc} */
157    public Collection<FieldEventHandler<T>> getEventHandlers() {
158        final List<FieldEventHandler<T>> list = new ArrayList<FieldEventHandler<T>>(eventsStates.size());
159        for (FieldEventState<T> state : eventsStates) {
160            list.add(state.getEventHandler());
161        }
162        return Collections.unmodifiableCollection(list);
163    }
164
165    /** {@inheritDoc} */
166    public void clearEventHandlers() {
167        eventsStates.clear();
168    }
169
170    /** {@inheritDoc} */
171    public FieldODEStateAndDerivative<T> getCurrentStepStart() {
172        return stepStart;
173    }
174
175    /** {@inheritDoc} */
176    public T getCurrentSignedStepsize() {
177        return stepSize;
178    }
179
180    /** {@inheritDoc} */
181    public void setMaxEvaluations(int maxEvaluations) {
182        evaluations = evaluations.withMaximalCount((maxEvaluations < 0) ? Integer.MAX_VALUE : maxEvaluations);
183    }
184
185    /** {@inheritDoc} */
186    public int getMaxEvaluations() {
187        return evaluations.getMaximalCount();
188    }
189
190    /** {@inheritDoc} */
191    public int getEvaluations() {
192        return evaluations.getCount();
193    }
194
195    /** Prepare the start of an integration.
196     * @param eqn equations to integrate
197     * @param t0 start value of the independent <i>time</i> variable
198     * @param y0 array containing the start value of the state vector
199     * @param t target time for the integration
200     * @return initial state with derivatives added
201     */
202    protected FieldODEStateAndDerivative<T> initIntegration(final FieldExpandableODE<T> eqn,
203                                                            final T t0, final T[] y0, final T t) {
204
205        this.equations = eqn;
206        evaluations    = evaluations.withStart(0);
207
208        // initialize ODE
209        eqn.init(t0, y0, t);
210
211        // set up derivatives of initial state
212        final T[] y0Dot = computeDerivatives(t0, y0);
213        final FieldODEStateAndDerivative<T> state0 = new FieldODEStateAndDerivative<T>(t0, y0, y0Dot);
214
215        // initialize event handlers
216        for (final FieldEventState<T> state : eventsStates) {
217            state.getEventHandler().init(state0, t);
218        }
219
220        // initialize step handlers
221        for (FieldStepHandler<T> handler : stepHandlers) {
222            handler.init(state0, t);
223        }
224
225        setStateInitialized(false);
226
227        return state0;
228
229    }
230
231    /** Get the differential equations to integrate.
232     * @return differential equations to integrate
233     */
234    protected FieldExpandableODE<T> getEquations() {
235        return equations;
236    }
237
238    /** Get the evaluations counter.
239     * @return evaluations counter
240     */
241    protected IntegerSequence.Incrementor getEvaluationsCounter() {
242        return evaluations;
243    }
244
245    /** Compute the derivatives and check the number of evaluations.
246     * @param t current value of the independent <I>time</I> variable
247     * @param y array containing the current value of the state vector
248     * @return state completed with derivatives
249     * @exception DimensionMismatchException if arrays dimensions do not match equations settings
250     * @exception MaxCountExceededException if the number of functions evaluations is exceeded
251     * @exception NullPointerException if the ODE equations have not been set (i.e. if this method
252     * is called outside of a call to {@link #integrate(FieldExpandableODE, FieldODEState,
253     * RealFieldElement) integrate}
254     */
255    public T[] computeDerivatives(final T t, final T[] y)
256        throws DimensionMismatchException, MaxCountExceededException, NullPointerException {
257        evaluations.increment();
258        return equations.computeDerivatives(t, y);
259    }
260
261    /** Set the stateInitialized flag.
262     * <p>This method must be called by integrators with the value
263     * {@code false} before they start integration, so a proper lazy
264     * initialization is done automatically on the first step.</p>
265     * @param stateInitialized new value for the flag
266     */
267    protected void setStateInitialized(final boolean stateInitialized) {
268        this.statesInitialized = stateInitialized;
269    }
270
271    /** Accept a step, triggering events and step handlers.
272     * @param interpolator step interpolator
273     * @param tEnd final integration time
274     * @return state at end of step
275     * @exception MaxCountExceededException if the interpolator throws one because
276     * the number of functions evaluations is exceeded
277     * @exception NoBracketingException if the location of an event cannot be bracketed
278     * @exception DimensionMismatchException if arrays dimensions do not match equations settings
279     */
280    protected FieldODEStateAndDerivative<T> acceptStep(final AbstractFieldStepInterpolator<T> interpolator,
281                                                       final T tEnd)
282        throws MaxCountExceededException, DimensionMismatchException, NoBracketingException {
283
284            FieldODEStateAndDerivative<T> previousState = interpolator.getGlobalPreviousState();
285            final FieldODEStateAndDerivative<T> currentState = interpolator.getGlobalCurrentState();
286
287            // initialize the events states if needed
288            if (! statesInitialized) {
289                for (FieldEventState<T> state : eventsStates) {
290                    state.reinitializeBegin(interpolator);
291                }
292                statesInitialized = true;
293            }
294
295            // search for next events that may occur during the step
296            final int orderingSign = interpolator.isForward() ? +1 : -1;
297            SortedSet<FieldEventState<T>> occurringEvents = new TreeSet<FieldEventState<T>>(new Comparator<FieldEventState<T>>() {
298
299                /** {@inheritDoc} */
300                public int compare(FieldEventState<T> es0, FieldEventState<T> es1) {
301                    return orderingSign * Double.compare(es0.getEventTime().getReal(), es1.getEventTime().getReal());
302                }
303
304            });
305
306            for (final FieldEventState<T> state : eventsStates) {
307                if (state.evaluateStep(interpolator)) {
308                    // the event occurs during the current step
309                    occurringEvents.add(state);
310                }
311            }
312
313            AbstractFieldStepInterpolator<T> restricted = interpolator;
314            while (!occurringEvents.isEmpty()) {
315
316                // handle the chronologically first event
317                final Iterator<FieldEventState<T>> iterator = occurringEvents.iterator();
318                final FieldEventState<T> currentEvent = iterator.next();
319                iterator.remove();
320
321                // get state at event time
322                final FieldODEStateAndDerivative<T> eventState = restricted.getInterpolatedState(currentEvent.getEventTime());
323
324                // restrict the interpolator to the first part of the step, up to the event
325                restricted = restricted.restrictStep(previousState, eventState);
326
327                // advance all event states to current time
328                for (final FieldEventState<T> state : eventsStates) {
329                    state.stepAccepted(eventState);
330                    isLastStep = isLastStep || state.stop();
331                }
332
333                // handle the first part of the step, up to the event
334                for (final FieldStepHandler<T> handler : stepHandlers) {
335                    handler.handleStep(restricted, isLastStep);
336                }
337
338                if (isLastStep) {
339                    // the event asked to stop integration
340                    return eventState;
341                }
342
343                FieldODEState<T> newState = null;
344                resetOccurred = false;
345                for (final FieldEventState<T> state : eventsStates) {
346                    newState = state.reset(eventState);
347                    if (newState != null) {
348                        // some event handler has triggered changes that
349                        // invalidate the derivatives, we need to recompute them
350                        final T[] y    = equations.getMapper().mapState(newState);
351                        final T[] yDot = computeDerivatives(newState.getTime(), y);
352                        resetOccurred = true;
353                        return equations.getMapper().mapStateAndDerivative(newState.getTime(), y, yDot);
354                    }
355                }
356
357                // prepare handling of the remaining part of the step
358                previousState = eventState;
359                restricted = restricted.restrictStep(eventState, currentState);
360
361                // check if the same event occurs again in the remaining part of the step
362                if (currentEvent.evaluateStep(restricted)) {
363                    // the event occurs during the current step
364                    occurringEvents.add(currentEvent);
365                }
366
367            }
368
369            // last part of the step, after the last event
370            for (final FieldEventState<T> state : eventsStates) {
371                state.stepAccepted(currentState);
372                isLastStep = isLastStep || state.stop();
373            }
374            isLastStep = isLastStep || currentState.getTime().subtract(tEnd).abs().getReal() <= FastMath.ulp(tEnd.getReal());
375
376            // handle the remaining part of the step, after all events if any
377            for (FieldStepHandler<T> handler : stepHandlers) {
378                handler.handleStep(restricted, isLastStep);
379            }
380
381            return currentState;
382
383    }
384
385    /** Check the integration span.
386     * @param eqn set of differential equations
387     * @param t target time for the integration
388     * @exception NumberIsTooSmallException if integration span is too small
389     * @exception DimensionMismatchException if adaptive step size integrators
390     * tolerance arrays dimensions are not compatible with equations settings
391     */
392    protected void sanityChecks(final FieldODEState<T> eqn, final T t)
393        throws NumberIsTooSmallException, DimensionMismatchException {
394
395        final double threshold = 1000 * FastMath.ulp(FastMath.max(FastMath.abs(eqn.getTime().getReal()),
396                                                                  FastMath.abs(t.getReal())));
397        final double dt = eqn.getTime().subtract(t).abs().getReal();
398        if (dt <= threshold) {
399            throw new NumberIsTooSmallException(LocalizedFormats.TOO_SMALL_INTEGRATION_INTERVAL,
400                                                dt, threshold, false);
401        }
402
403    }
404
405    /** Check if a reset occurred while last step was accepted.
406     * @return true if a reset occurred while last step was accepted
407     */
408    protected boolean resetOccurred() {
409        return resetOccurred;
410    }
411
412    /** Set the current step size.
413     * @param stepSize step size to set
414     */
415    protected void setStepSize(final T stepSize) {
416        this.stepSize = stepSize;
417    }
418
419    /** Get the current step size.
420     * @return current step size
421     */
422    protected T getStepSize() {
423        return stepSize;
424    }
425    /** Set current step start.
426     * @param stepStart step start
427     */
428    protected void setStepStart(final FieldODEStateAndDerivative<T> stepStart) {
429        this.stepStart = stepStart;
430    }
431
432    /** Getcurrent step start.
433     * @return current step start
434     */
435    protected FieldODEStateAndDerivative<T> getStepStart() {
436        return stepStart;
437    }
438
439    /** Set the last state flag.
440     * @param isLastStep if true, this step is the last one
441     */
442    protected void setIsLastStep(final boolean isLastStep) {
443        this.isLastStep = isLastStep;
444    }
445
446    /** Check if this step is the last one.
447     * @return true if this step is the last one
448     */
449    protected boolean isLastStep() {
450        return isLastStep;
451    }
452
453}