/*
 * 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.math.ode;

import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.ode.nonstiff.AdaptiveStepsizeIntegrator;
import org.apache.commons.math.ode.nonstiff.DormandPrince853Integrator;
import org.apache.commons.math.ode.sampling.StepHandler;
import org.apache.commons.math.ode.sampling.StepInterpolator;

/**
 * This class is the base class for multistep integrators for Ordinary
 * Differential Equations.
 *
 * @see org.apache.commons.math.ode.nonstiff.AdamsBashforthIntegrator
 * @see org.apache.commons.math.ode.nonstiff.AdamsMoultonIntegrator
 * @version $Revision$ $Date$
 * @since 2.0
 */
public abstract class MultistepIntegrator extends AdaptiveStepsizeIntegrator {

	/** Starter integrator. */
	private FirstOrderIntegrator starter;

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

	/** 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(k))
	 * </p>
	 */
	protected Array2DRowRealMatrix nordsieck;

	/** 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>
	 * 
	 * @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
	 */
	protected MultistepIntegrator(final String name, final int nSteps, final int order, final double minStep,
			final double maxStep, final double scalAbsoluteTolerance, final double scalRelativeTolerance) {

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

		if (nSteps <= 0) {
			throw MathRuntimeException.createIllegalArgumentException("{0} method needs at least one previous point",
					name);
		}

		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(10.0);

	}

	/**
	 * 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>
	 * 
	 * @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(10.0);

	}

	/**
	 * 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 starter
	 *            starter integrator
	 */
	public void setStarterIntegrator(FirstOrderIntegrator starter) {
		this.starter = starter;
	}

	/**
	 * 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)
	 * @throws IntegratorException
	 *             if the integrator cannot perform integration
	 * @throws DerivativeException
	 *             this exception is propagated to the caller if the underlying
	 *             user function triggers one
	 */
	protected void start(final double t0, final double[] y0, final double t)
			throws DerivativeException, IntegratorException {

		// 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(y0.length));

		// start integration, expecting a InitializationCompletedMarkerException
		try {
			starter.integrate(new CountingDifferentialEquations(y0.length), t0, y0, t, new double[y0.length]);
		} catch (DerivativeException de) {
			if (!(de instanceof InitializationCompletedMarkerException)) {
				// this is not the expected nominal interruption of the start
				// integrator
				throw de;
			}
		}

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

	}

	/**
	 * Initialize the high order scaled derivatives at step start.
	 * 
	 * @param first
	 *            first scaled derivative at step start
	 * @param multistep
	 *            scaled derivatives after step start (hy'1, ..., hy'k-1) will
	 *            be modified
	 * @return high order scaled derivatives at step start
	 */
	protected abstract Array2DRowRealMatrix initializeHighOrderDerivatives(final double[] first,
			final double[][] multistep);

	/**
	 * 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;
	}

	/**
	 * 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 Math.min(maxGrowth, Math.max(minReduction, safety * Math.pow(error, exp)));
	}

	/**
	 * Transformer used to convert the first step to Nordsieck representation.
	 */
	public static interface NordsieckTransformer {
		/**
		 * Initialize the high order scaled derivatives at step start.
		 * 
		 * @param first
		 *            first scaled derivative at step start
		 * @param multistep
		 *            scaled derivatives after step start (hy'1, ..., hy'k-1)
		 *            will be modified
		 * @return high order derivatives at step start
		 */
		RealMatrix initializeHighOrderDerivatives(double[] first, double[][] multistep);
	}

	/** Specialized step handler storing the first step. */
	private class NordsieckInitializer implements StepHandler {

		/** Problem dimension. */
		private final int n;

		/**
		 * Simple constructor.
		 * 
		 * @param n
		 *            problem dimension
		 */
		public NordsieckInitializer(final int n) {
			this.n = n;
		}

		/** {@inheritDoc} */
		public void handleStep(StepInterpolator interpolator, boolean isLast) throws DerivativeException {

			final double prev = interpolator.getPreviousTime();
			final double curr = interpolator.getCurrentTime();
			stepStart = prev;
			stepSize = (curr - prev) / (nSteps + 1);

			// compute the first scaled derivative
			interpolator.setInterpolatedTime(prev);
			scaled = interpolator.getInterpolatedDerivatives().clone();
			for (int j = 0; j < n; ++j) {
				scaled[j] *= stepSize;
			}

			// compute the high order scaled derivatives
			final double[][] multistep = new double[nSteps][];
			for (int i = 1; i <= nSteps; ++i) {
				interpolator.setInterpolatedTime(prev + stepSize * i);
				final double[] msI = interpolator.getInterpolatedDerivatives().clone();
				for (int j = 0; j < n; ++j) {
					msI[j] *= stepSize;
				}
				multistep[i - 1] = msI;
			}
			nordsieck = initializeHighOrderDerivatives(scaled, multistep);

			// stop the integrator after the first step has been handled
			throw new InitializationCompletedMarkerException();

		}

		/** {@inheritDoc} */
		public boolean requiresDenseOutput() {
			return true;
		}

		/** {@inheritDoc} */
		public void reset() {
			// nothing to do
		}

	}

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

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

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

	}

	/**
	 * Wrapper for differential equations, ensuring start evaluations are
	 * counted.
	 */
	private class CountingDifferentialEquations implements FirstOrderDifferentialEquations {

		/** Dimension of the problem. */
		private final int dimension;

		/**
		 * Simple constructor.
		 * 
		 * @param dimension
		 *            dimension of the problem
		 */
		public CountingDifferentialEquations(final int dimension) {
			this.dimension = dimension;
		}

		/** {@inheritDoc} */
		public void computeDerivatives(double t, double[] y, double[] dot) throws DerivativeException {
			MultistepIntegrator.this.computeDerivatives(t, y, dot);
		}

		/** {@inheritDoc} */
		public int getDimension() {
			return dimension;
		}
	}

}
