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