/*
    * 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.math4.legacy.ode;
  
  import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  import org.apache.commons.math4.legacy.exception.MathIllegalStateException;
  import org.apache.commons.math4.legacy.exception.MaxCountExceededException;
  import org.apache.commons.math4.legacy.exception.NoBracketingException;
  import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
  import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
  import org.apache.commons.math4.legacy.ode.nonstiff.AdaptiveStepsizeIntegrator;
  import org.apache.commons.math4.legacy.ode.nonstiff.DormandPrince853Integrator;
  import org.apache.commons.math4.legacy.ode.sampling.StepHandler;
  import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
  import org.apache.commons.math4.core.jdkmath.JdkMath;
  
  /**
   * 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:
   * <div style="white-space: pre"><code>
   * 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
   * </code></div>
   * <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:
   * <div style="white-space: pre"><code>
   * r<sub>n</sub> = [ s<sub>2</sub>(n), s<sub>3</sub>(n) ... s<sub>k</sub>(n) ]<sup>T</sup>
   * </code></div>
   * (we omit the k index in the notation for clarity)
   * <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.math4.legacy.ode.nonstiff.AdamsBashforthIntegrator
   * @see org.apache.commons.math4.legacy.ode.nonstiff.AdamsMoultonIntegrator
   * @since 2.0
   */
  public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {
  
      /** First scaled derivative (h y'). */
      protected double[] scaled;
  
      /** 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 Array2DRowRealMatrix nordsieck;
  
      /** Starter integrator. */
      private FirstOrderIntegrator starter;
  
      /** Number of steps of the multistep method (excluding the one being computed). */
      private final int nSteps;
  
      /** Stepsize control exponent. */
      private double exp;
  
      /** Safety factor for stepsize control. */
      private double safety;
  
      /** Minimal reduction factor for stepsize control. */
      private double minReduction;
  
      /** Maximal growth factor for stepsize control. */
      private double maxGrowth;
  
      /**
       * Build a multistep integrator with the given stepsize bounds.
       * <p>The default starter integrator is set to the {@link
       * DormandPrince853Integrator 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 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 MultistepIntegrator(final String name, final int nSteps,
                                   final int order,
                                   final double minStep, final double maxStep,
                                   final double scalAbsoluteTolerance,
                                   final double scalRelativeTolerance)
         throws NumberIsTooSmallException {
 
         super(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 DormandPrince853Integrator(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(JdkMath.pow(2.0, -exp));
     }
 
     /**
      * Build a multistep integrator with the given stepsize bounds.
      * <p>The default starter integrator is set to the {@link
      * DormandPrince853Integrator 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 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 MultistepIntegrator(final String name, final int nSteps,
                                   final int order,
                                   final double minStep, final double maxStep,
                                   final double[] vecAbsoluteTolerance,
                                   final double[] vecRelativeTolerance) {
         super(name, minStep, maxStep, vecAbsoluteTolerance, vecRelativeTolerance);
         starter = new DormandPrince853Integrator(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(JdkMath.pow(2.0, -exp));
     }
 
     /**
      * Get the starter integrator.
      * @return starter integrator
      */
     public ODEIntegrator getStarterIntegrator() {
         return starter;
     }
 
     /**
      * 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(FirstOrderIntegrator starterIntegrator) {
         this.starter = starterIntegrator;
     }
 
     /** 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 t0 initial time
      * @param y0 initial value of the state vector at t0
      * @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 double t0, final double[] y0, final double 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 NordsieckInitializer((nSteps + 3) / 2, y0.length));
 
         // start integration, expecting a InitializationCompletedMarkerException
         try {
 
             if (starter instanceof AbstractIntegrator) {
                 ((AbstractIntegrator) starter).integrate(getExpandable(), t);
             } else {
                 starter.integrate(new FirstOrderDifferentialEquations() {
 
                     /** {@inheritDoc} */
                     @Override
                     public int getDimension() {
                         return getExpandable().getTotalDimension();
                     }
 
                     /** {@inheritDoc} */
                     @Override
                     public void computeDerivatives(double t, double[] y, double[] yDot) {
                         getExpandable().computeDerivatives(t, y, yDot);
                     }
                 }, t0, y0, t, new double[y0.length]);
             }
 
             // 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
             getCounter().increment(starter.getEvaluations());
         }
 
         // remove the specific step handler
         starter.clearStepHandlers();
     }
 
     /** 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 Array2DRowRealMatrix initializeHighOrderDerivatives(double h, double[] t,
                                                                            double[][] y,
                                                                            double[][] yDot);
 
     /** 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.
      * @param minReduction minimal reduction factor
      */
     public void setMinReduction(final double minReduction) {
         this.minReduction = minReduction;
     }
 
     /** 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.
      * @param maxGrowth maximal growth factor
      */
     public void setMaxGrowth(final double maxGrowth) {
         this.maxGrowth = maxGrowth;
     }
 
     /** Get the safety factor for stepsize control.
      * @return safety factor
      */
     public double getSafety() {
       return safety;
     }
 
     /** 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).
      * @return number of steps of the multistep method (excluding the one being computed)
      */
     public int getNSteps() {
       return nSteps;
     }
 
     /** 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 double computeStepGrowShrinkFactor(final double error) {
         return JdkMath.min(maxGrowth, JdkMath.max(minReduction, safety * JdkMath.pow(error, exp)));
     }
 
     /** Transformer used to convert the first step to Nordsieck representation.
      * @deprecated as of 3.6 this unused interface is deprecated
      */
     @Deprecated
     public interface NordsieckTransformer {
         /** 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>)
          */
         Array2DRowRealMatrix initializeHighOrderDerivatives(double h, double[] t,
                                                             double[][] y,
                                                             double[][] yDot);
     }
 
     /** Specialized step handler storing the first step. */
     private final class NordsieckInitializer implements StepHandler {
 
         /** Steps counter. */
         private int count;
 
         /** First steps times. */
         private final double[] t;
 
         /** First steps states. */
         private final double[][] y;
 
         /** First steps derivatives. */
         private final double[][] yDot;
 
         /** Simple constructor.
          * @param nbStartPoints number of start points (including the initial point)
          * @param n problem dimension
          */
         NordsieckInitializer(final int nbStartPoints, final int n) {
             this.count = 0;
             this.t     = new double[nbStartPoints];
             this.y     = new double[nbStartPoints][n];
             this.yDot  = new double[nbStartPoints][n];
         }
 
         /** {@inheritDoc} */
         @Override
         public void handleStep(StepInterpolator interpolator, boolean isLast)
             throws MaxCountExceededException {
 
             final double prev = interpolator.getPreviousTime();
             final double curr = interpolator.getCurrentTime();
 
             if (count == 0) {
                 // first step, we need to store also the point at the beginning of the step
                 interpolator.setInterpolatedTime(prev);
                 t[0] = prev;
                 final ExpandableStatefulODE expandable = getExpandable();
                 final EquationsMapper primary = expandable.getPrimaryMapper();
                 primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
                 primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
                 int index = 0;
                 for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
                     secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
                     secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
                     ++index;
                 }
             }
 
             // store the point at the end of the step
             ++count;
             interpolator.setInterpolatedTime(curr);
             t[count] = curr;
 
             final ExpandableStatefulODE expandable = getExpandable();
             final EquationsMapper primary = expandable.getPrimaryMapper();
             primary.insertEquationData(interpolator.getInterpolatedState(), y[count]);
             primary.insertEquationData(interpolator.getInterpolatedDerivatives(), yDot[count]);
             int index = 0;
             for (final EquationsMapper secondary : expandable.getSecondaryMappers()) {
                 secondary.insertEquationData(interpolator.getInterpolatedSecondaryState(index), y[count]);
                 secondary.insertEquationData(interpolator.getInterpolatedSecondaryDerivatives(index), yDot[count]);
                 ++index;
             }
 
             if (count == t.length - 1) {
 
                 // this was the last point we needed, we can compute the derivatives
                 stepStart = t[0];
                 stepSize  = (t[t.length - 1] - t[0]) / (t.length - 1);
 
                 // first scaled derivative
                 scaled = yDot[0].clone();
                 for (int j = 0; j < scaled.length; ++j) {
                     scaled[j] *= stepSize;
                 }
 
                 // higher order derivatives
                 nordsieck = initializeHighOrderDerivatives(stepSize, t, y, yDot);
 
                 // stop the integrator now that all needed steps have been handled
                 throw new InitializationCompletedMarkerException();
             }
         }
 
         /** {@inheritDoc} */
         @Override
         public void init(double t0, double[] y0, double time) {
             // nothing to do
         }
     }
 
     /** Marker exception used ONLY to stop the starter integrator after first step. */
     private static final class InitializationCompletedMarkerException
         extends RuntimeException {
 
         /** Serializable version identifier. */
         private static final long serialVersionUID = -1914085471038046418L;
 
         /** Simple constructor. */
         InitializationCompletedMarkerException() {
             super((Throwable) null);
         }
     }
 }