/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math3.ode;

import org.apache.commons.math3.Field;
import org.apache.commons.math3.RealFieldElement;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.MathIllegalStateException;
import org.apache.commons.math3.exception.MaxCountExceededException;
import org.apache.commons.math3.exception.NoBracketingException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.linear.Array2DRowFieldMatrix;
import org.apache.commons.math3.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
import org.apache.commons.math3.ode.nonstiff.DormandPrince853FieldIntegrator;
import org.apache.commons.math3.ode.sampling.FieldStepHandler;
import org.apache.commons.math3.ode.sampling.FieldStepInterpolator;
import org.apache.commons.math3.util.FastMath;
import org.apache.commons.math3.util.MathArrays;
import org.apache.commons.math3.util.MathUtils;

This class is the base class for multistep integrators for Ordinary
Differential Equations.
We define scaled derivatives si(n) at step n as:
s1(n) = h y'n for first derivative
s2(n) = h2/2 y''n for second derivative
s3(n) = h3/6 y'''n for third derivative
...
sk(n) = hk/k! y(k)n for kth derivative

Rather than storing several previous steps separately, this implementation uses
the Nordsieck vector with higher degrees scaled derivatives all taken at the same
step (yn, s1(n) and rn) where rn is defined as:
rn = [ s2(n), s3(n) ... sk(n) ]T

(we omit the k index in the notation for clarity)


Multistep integrators with Nordsieck representation are highly sensitive to
large step changes because when the step is multiplied by factor a, the
kth component of the Nordsieck vector is multiplied by ak
and the last components are the least accurate ones. The default max growth
factor is therefore set to a quite low value: 21/order.

Type parameters:
<T> – the type of the field elements
See Also:
AdamsBashforthFieldIntegrator
AdamsMoultonFieldIntegrator
Since:
3.6/**
 * This class is the base class for multistep integrators for Ordinary
 * Differential Equations.
 * <p>We define scaled derivatives s<sub>i</sub>(n) at step n as:
 * <pre>
 * s<sub>1</sub>(n) = h y'<sub>n</sub> for first derivative
 * s<sub>2</sub>(n) = h<sup>2</sup>/2 y''<sub>n</sub> for second derivative
 * s<sub>3</sub>(n) = h<sup>3</sup>/6 y'''<sub>n</sub> for third derivative
 * ...
 * s<sub>k</sub>(n) = h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub> for k<sup>th</sup> derivative
 * </pre></p>
 * <p>Rather than storing several previous steps separately, this implementation uses
 * the Nordsieck vector with higher degrees scaled derivatives all taken at the same
 * step (y<sub>n</sub>, s<sub>1</sub>(n) and r<sub>n</sub>) where r<sub>n</sub> is defined as:
 * <pre>
 * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
 * </pre>
 * (we omit the k index in the notation for clarity)</p>
 * <p>
 * Multistep integrators with Nordsieck representation are highly sensitive to
 * large step changes because when the step is multiplied by factor a, the
 * k<sup>th</sup> component of the Nordsieck vector is multiplied by a<sup>k</sup>
 * and the last components are the least accurate ones. The default max growth
 * factor is therefore set to a quite low value: 2<sup>1/order</sup>.
 * </p>
 *
 * @see org.apache.commons.math3.ode.nonstiff.AdamsBashforthFieldIntegrator
 * @see org.apache.commons.math3.ode.nonstiff.AdamsMoultonFieldIntegrator
 * @param <T> the type of the field elements
 * @since 3.6
 */
public abstract class MultistepFieldIntegrator<T extends RealFieldElement<T>>
    extends AdaptiveStepsizeFieldIntegrator<T> {

    
First scaled derivative (h y'). /** First scaled derivative (h y'). */
    protected T[] scaled;

    
Nordsieck matrix of the higher scaled derivatives.
(h2/2 y'', h3/6 y''' ..., hk/k! y(k))
/** Nordsieck matrix of the higher scaled derivatives.
     * <p>(h<sup>2</sup>/2 y'', h<sup>3</sup>/6 y''' ..., h<sup>k</sup>/k! y<sup>(k)</sup>)</p>
     */
    protected Array2DRowFieldMatrix<T> nordsieck;

    
Starter integrator. /** Starter integrator. */
    private FirstOrderFieldIntegrator<T> starter;

    
Number of steps of the multistep method (excluding the one being computed). /** Number of steps of the multistep method (excluding the one being computed). */
    private final int nSteps;

    
Stepsize control exponent. /** Stepsize control exponent. */
    private double exp;

    
Safety factor for stepsize control. /** Safety factor for stepsize control. */
    private double safety;

    
Minimal reduction factor for stepsize control. /** Minimal reduction factor for stepsize control. */
    private double minReduction;

    
Maximal growth factor for stepsize control. /** Maximal growth factor for stepsize control. */
    private double maxGrowth;

    
Build a multistep integrator with the given stepsize bounds.
The default starter integrator is set to the Dormand-Prince 8(5,3) integrator with some defaults settings.

The default max growth factor is set to a quite low value: 21/order.

Params:
field – field to which the time and state vector elements belong
name – name of the method
nSteps – number of steps of the multistep method
(excluding the one being computed)
order – order of the method
minStep – minimal step (must be positive even for backward
integration), the last step can be smaller than this
maxStep – maximal step (must be positive even for backward
integration)
scalAbsoluteTolerance – allowed absolute error
scalRelativeTolerance – allowed relative error
Throws:
NumberIsTooSmallException – if number of steps is smaller than 2/**
     * Build a multistep integrator with the given stepsize bounds.
     * <p>The default starter integrator is set to the {@link
     * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
     * some defaults settings.</p>
     * <p>
     * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
     * </p>
     * @param field field to which the time and state vector elements belong
     * @param name name of the method
     * @param nSteps number of steps of the multistep method
     * (excluding the one being computed)
     * @param order order of the method
     * @param minStep minimal step (must be positive even for backward
     * integration), the last step can be smaller than this
     * @param maxStep maximal step (must be positive even for backward
     * integration)
     * @param scalAbsoluteTolerance allowed absolute error
     * @param scalRelativeTolerance allowed relative error
     * @exception NumberIsTooSmallException if number of steps is smaller than 2
     */
    protected MultistepFieldIntegrator(final Field<T> field, final String name,
                                       final int nSteps, final int order,
                                       final double minStep, final double maxStep,
                                       final double scalAbsoluteTolerance,
                                       final double scalRelativeTolerance)
        throws NumberIsTooSmallException {

        super(field, name, minStep, maxStep, scalAbsoluteTolerance, scalRelativeTolerance);

        if (nSteps < 2) {
            throw new NumberIsTooSmallException(
                  LocalizedFormats.INTEGRATION_METHOD_NEEDS_AT_LEAST_TWO_PREVIOUS_POINTS,
                  nSteps, 2, true);
        }

        starter = new DormandPrince853FieldIntegrator<T>(field, minStep, maxStep,
                                                         scalAbsoluteTolerance,
                                                         scalRelativeTolerance);
        this.nSteps = nSteps;

        exp = -1.0 / order;

        // set the default values of the algorithm control parameters
        setSafety(0.9);
        setMinReduction(0.2);
        setMaxGrowth(FastMath.pow(2.0, -exp));

    }

    
Build a multistep integrator with the given stepsize bounds.
The default starter integrator is set to the Dormand-Prince 8(5,3) integrator with some defaults settings.

The default max growth factor is set to a quite low value: 21/order.

Params:
field – field to which the time and state vector elements belong
name – name of the method
nSteps – number of steps of the multistep method
(excluding the one being computed)
order – order of the method
minStep – minimal step (must be positive even for backward
integration), the last step can be smaller than this
maxStep – maximal step (must be positive even for backward
integration)
vecAbsoluteTolerance – allowed absolute error
vecRelativeTolerance – allowed relative error/**
     * Build a multistep integrator with the given stepsize bounds.
     * <p>The default starter integrator is set to the {@link
     * DormandPrince853FieldIntegrator Dormand-Prince 8(5,3)} integrator with
     * some defaults settings.</p>
     * <p>
     * The default max growth factor is set to a quite low value: 2<sup>1/order</sup>.
     * </p>
     * @param field field to which the time and state vector elements belong
     * @param name name of the method
     * @param nSteps number of steps of the multistep method
     * (excluding the one being computed)
     * @param order order of the method
     * @param minStep minimal step (must be positive even for backward
     * integration), the last step can be smaller than this
     * @param maxStep maximal step (must be positive even for backward
     * integration)
     * @param vecAbsoluteTolerance allowed absolute error
     * @param vecRelativeTolerance allowed relative error
     */
    protected MultistepFieldIntegrator(final Field<T> field, final String name, final int nSteps,
                                       final int order,
                                       final double minStep, final double maxStep,
                                       final double[] vecAbsoluteTolerance,
                                       final double[] vecRelativeTolerance) {
        super(field, name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
        starter = new DormandPrince853FieldIntegrator<T>(field, minStep, maxStep,
                                                         vecAbsoluteTolerance,
                                                         vecRelativeTolerance);
        this.nSteps = nSteps;

        exp = -1.0 / order;

        // set the default values of the algorithm control parameters
        setSafety(0.9);
        setMinReduction(0.2);
        setMaxGrowth(FastMath.pow(2.0, -exp));

    }

    
Get the starter integrator.
Returns:
starter integrator/**
     * Get the starter integrator.
     * @return starter integrator
     */
    public FirstOrderFieldIntegrator<T> getStarterIntegrator() {
        return starter;
    }

    
Set the starter integrator.
The various step and event handlers for this starter integrator
will be managed automatically by the multi-step integrator. Any
user configuration for these elements will be cleared before use.
Params:
starterIntegrator – starter integrator/**
     * Set the starter integrator.
     * <p>The various step and event handlers for this starter integrator
     * will be managed automatically by the multi-step integrator. Any
     * user configuration for these elements will be cleared before use.</p>
     * @param starterIntegrator starter integrator
     */
    public void setStarterIntegrator(FirstOrderFieldIntegrator<T> starterIntegrator) {
        this.starter = starterIntegrator;
    }

    
Start the integration.
This method computes one step using the underlying starter integrator,
and initializes the Nordsieck vector at step start. The starter integrator
purpose is only to establish initial conditions, it does not really change
time by itself. The top level multistep integrator remains in charge of
handling time propagation and events handling as it will starts its own
computation right from the beginning. In a sense, the starter integrator
can be seen as a dummy one and so it will never trigger any user event nor
call any user step handler.
Params:
equations – complete set of differential equations to integrate
initialState – initial state (time, primary and secondary state vectors)
t – target time for the integration
(can be set to a value smaller than t0 for backward integration)
Throws:
DimensionMismatchException – if arrays dimension do not match equations settings
NumberIsTooSmallException – if integration step is too small
MaxCountExceededException – if the number of functions evaluations is exceeded
NoBracketingException – if the location of an event cannot be bracketed/** Start the integration.
     * <p>This method computes one step using the underlying starter integrator,
     * and initializes the Nordsieck vector at step start. The starter integrator
     * purpose is only to establish initial conditions, it does not really change
     * time by itself. The top level multistep integrator remains in charge of
     * handling time propagation and events handling as it will starts its own
     * computation right from the beginning. In a sense, the starter integrator
     * can be seen as a dummy one and so it will never trigger any user event nor
     * call any user step handler.</p>
     * @param equations complete set of differential equations to integrate
     * @param initialState initial state (time, primary and secondary state vectors)
     * @param t target time for the integration
     * (can be set to a value smaller than <code>t0</code> for backward integration)
     * @exception DimensionMismatchException if arrays dimension do not match equations settings
     * @exception NumberIsTooSmallException if integration step is too small
     * @exception MaxCountExceededException if the number of functions evaluations is exceeded
     * @exception NoBracketingException if the location of an event cannot be bracketed
     */
    protected void start(final FieldExpandableODE<T> equations, final FieldODEState<T> initialState, final T t)
        throws DimensionMismatchException, NumberIsTooSmallException,
               MaxCountExceededException, NoBracketingException {

        // make sure NO user event nor user step handler is triggered,
        // this is the task of the top level integrator, not the task
        // of the starter integrator
        starter.clearEventHandlers();
        starter.clearStepHandlers();

        // set up one specific step handler to extract initial Nordsieck vector
        starter.addStepHandler(new FieldNordsieckInitializer(equations.getMapper(), (nSteps + 3) / 2));

        // start integration, expecting a InitializationCompletedMarkerException
        try {

            starter.integrate(equations, initialState, t);

            // we should not reach this step
            throw new MathIllegalStateException(LocalizedFormats.MULTISTEP_STARTER_STOPPED_EARLY);

        } catch (InitializationCompletedMarkerException icme) { // NOPMD
            // this is the expected nominal interruption of the start integrator

            // count the evaluations used by the starter
            getEvaluationsCounter().increment(starter.getEvaluations());

        }

        // remove the specific step handler
        starter.clearStepHandlers();

    }

    
Initialize the high order scaled derivatives at step start.
Params:
h – step size to use for scaling
t – first steps times
y – first steps states
yDot – first steps derivatives
Returns:
Nordieck vector at first step (h2/2 y''n,
h3/6 y'''n ... hk/k! y(k)n)/** Initialize the high order scaled derivatives at step start.
     * @param h step size to use for scaling
     * @param t first steps times
     * @param y first steps states
     * @param yDot first steps derivatives
     * @return Nordieck vector at first step (h<sup>2</sup>/2 y''<sub>n</sub>,
     * h<sup>3</sup>/6 y'''<sub>n</sub> ... h<sup>k</sup>/k! y<sup>(k)</sup><sub>n</sub>)
     */
    protected abstract Array2DRowFieldMatrix<T> initializeHighOrderDerivatives(final T h, final T[] t,
                                                                               final T[][] y,
                                                                               final T[][] yDot);

    
Get the minimal reduction factor for stepsize control.
Returns:
minimal reduction factor/** Get the minimal reduction factor for stepsize control.
     * @return minimal reduction factor
     */
    public double getMinReduction() {
        return minReduction;
    }

    
Set the minimal reduction factor for stepsize control.
Params:
minReduction – minimal reduction factor/** Set the minimal reduction factor for stepsize control.
     * @param minReduction minimal reduction factor
     */
    public void setMinReduction(final double minReduction) {
        this.minReduction = minReduction;
    }

    
Get the maximal growth factor for stepsize control.
Returns:
maximal growth factor/** Get the maximal growth factor for stepsize control.
     * @return maximal growth factor
     */
    public double getMaxGrowth() {
        return maxGrowth;
    }

    
Set the maximal growth factor for stepsize control.
Params:
maxGrowth – maximal growth factor/** Set the maximal growth factor for stepsize control.
     * @param maxGrowth maximal growth factor
     */
    public void setMaxGrowth(final double maxGrowth) {
        this.maxGrowth = maxGrowth;
    }

    
Get the safety factor for stepsize control.
Returns:
safety factor/** Get the safety factor for stepsize control.
     * @return safety factor
     */
    public double getSafety() {
      return safety;
    }

    
Set the safety factor for stepsize control.
Params:
safety – safety factor/** Set the safety factor for stepsize control.
     * @param safety safety factor
     */
    public void setSafety(final double safety) {
      this.safety = safety;
    }

    
Get the number of steps of the multistep method (excluding the one being computed).
Returns:
number of steps of the multistep method (excluding the one being computed)/** Get the number of steps of the multistep method (excluding the one being computed).
     * @return number of steps of the multistep method (excluding the one being computed)
     */
    public int getNSteps() {
      return nSteps;
    }

    
Rescale the instance.
Since the scaled and Nordsieck arrays are shared with the caller,
this method has the side effect of rescaling this arrays in the caller too.
Params:
newStepSize – new step size to use in the scaled and Nordsieck arrays/** Rescale the instance.
     * <p>Since the scaled and Nordsieck arrays are shared with the caller,
     * this method has the side effect of rescaling this arrays in the caller too.</p>
     * @param newStepSize new step size to use in the scaled and Nordsieck arrays
     */
    protected void rescale(final T newStepSize) {

        final T ratio = newStepSize.divide(getStepSize());
        for (int i = 0; i < scaled.length; ++i) {
            scaled[i] = scaled[i].multiply(ratio);
        }

        final T[][] nData = nordsieck.getDataRef();
        T power = ratio;
        for (int i = 0; i < nData.length; ++i) {
            power = power.multiply(ratio);
            final T[] nDataI = nData[i];
            for (int j = 0; j < nDataI.length; ++j) {
                nDataI[j] = nDataI[j].multiply(power);
            }
        }

        setStepSize(newStepSize);

    }


    
Compute step grow/shrink factor according to normalized error.
Params:
error – normalized error of the current step
Returns:
grow/shrink factor for next step/** Compute step grow/shrink factor according to normalized error.
     * @param error normalized error of the current step
     * @return grow/shrink factor for next step
     */
    protected T computeStepGrowShrinkFactor(final T error) {
        return MathUtils.min(error.getField().getZero().add(maxGrowth),
                             MathUtils.max(error.getField().getZero().add(minReduction),
                                           error.pow(exp).multiply(safety)));
    }

    
Specialized step handler storing the first step.
/** Specialized step handler storing the first step.
     */
    private class FieldNordsieckInitializer implements FieldStepHandler<T> {

        
Equation mapper. /** Equation mapper. */
        private final FieldEquationsMapper<T> mapper;

        
Steps counter. /** Steps counter. */
        private int count;

        
Saved start. /** Saved start. */
        private FieldODEStateAndDerivative<T> savedStart;

        
First steps times. /** First steps times. */
        private final T[] t;

        
First steps states. /** First steps states. */
        private final T[][] y;

        
First steps derivatives. /** First steps derivatives. */
        private final T[][] yDot;

        
Simple constructor.
Params:
mapper – equation mapper
nbStartPoints – number of start points (including the initial point)/** Simple constructor.
         * @param mapper equation mapper
         * @param nbStartPoints number of start points (including the initial point)
         */
        FieldNordsieckInitializer(final FieldEquationsMapper<T> mapper, final int nbStartPoints) {
            this.mapper = mapper;
            this.count  = 0;
            this.t      = MathArrays.buildArray(getField(), nbStartPoints);
            this.y      = MathArrays.buildArray(getField(), nbStartPoints, -1);
            this.yDot   = MathArrays.buildArray(getField(), nbStartPoints, -1);
        }

        
{@inheritDoc} /** {@inheritDoc} */
        public void handleStep(FieldStepInterpolator<T> interpolator, boolean isLast)
            throws MaxCountExceededException {


            if (count == 0) {
                // first step, we need to store also the point at the beginning of the step
                final FieldODEStateAndDerivative<T> prev = interpolator.getPreviousState();
                savedStart  = prev;
                t[count]    = prev.getTime();
                y[count]    = mapper.mapState(prev);
                yDot[count] = mapper.mapDerivative(prev);
            }

            // store the point at the end of the step
            ++count;
            final FieldODEStateAndDerivative<T> curr = interpolator.getCurrentState();
            t[count]    = curr.getTime();
            y[count]    = mapper.mapState(curr);
            yDot[count] = mapper.mapDerivative(curr);

            if (count == t.length - 1) {

                // this was the last point we needed, we can compute the derivatives
                setStepSize(t[t.length - 1].subtract(t[0]).divide(t.length - 1));

                // first scaled derivative
                scaled = MathArrays.buildArray(getField(), yDot[0].length);
                for (int j = 0; j < scaled.length; ++j) {
                    scaled[j] = yDot[0][j].multiply(getStepSize());
                }

                // higher order derivatives
                nordsieck = initializeHighOrderDerivatives(getStepSize(), t, y, yDot);

                // stop the integrator now that all needed steps have been handled
                setStepStart(savedStart);
                throw new InitializationCompletedMarkerException();

            }

        }

        
{@inheritDoc} /** {@inheritDoc} */
        public void init(final FieldODEStateAndDerivative<T> initialState, T finalTime) {
            // nothing to do
        }

    }

    
Marker exception used ONLY to stop the starter integrator after first step. /** Marker exception used ONLY to stop the starter integrator after first step. */
    private static class InitializationCompletedMarkerException
        extends RuntimeException {

        
Serializable version identifier. /** Serializable version identifier. */
        private static final long serialVersionUID = -1914085471038046418L;

        
Simple constructor. /** Simple constructor. */
        InitializationCompletedMarkerException() {
            super((Throwable) null);
        }

    }

}