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

import java.io.IOException;
import java.io.ObjectInput;
import java.io.ObjectOutput;

import org.apache.commons.math.ode.DerivativeException;
import org.apache.commons.math.ode.sampling.AbstractStepInterpolator;
import org.apache.commons.math.ode.sampling.StepInterpolator;

/**
 * This class implements an interpolator for the Gragg-Bulirsch-Stoer
 * integrator.
 *
 * <p>
 * This interpolator compute dense output inside the last step produced by a
 * Gragg-Bulirsch-Stoer integrator.
 * </p>
 *
 * <p>
 * This implementation is basically a reimplementation in Java of the
 * <a href="http://www.unige.ch/math/folks/hairer/prog/nonstiff/odex.f">odex</a>
 * fortran code by E. Hairer and G. Wanner. The redistribution policy for this
 * code is available
 * <a href="http://www.unige.ch/~hairer/prog/licence.txt">here</a>, for
 * convenience, it is reproduced below.
 * </p>
 * </p>
 *
 * <table border="0" width="80%" cellpadding="10" align="center" bgcolor=
 * "#E0E0E0">
 * <tr>
 * <td>Copyright (c) 2004, Ernst Hairer</td>
 * </tr>
 *
 * <tr>
 * <td>Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * <ul>
 * <li>Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.</li>
 * <li>Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.</li>
 * </ul>
 * </td>
 * </tr>
 *
 * <tr>
 * <td><strong>THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT
 * NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR
 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
 * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
 * OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
 * ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.</strong></td>
 * </tr>
 * </table>
 *
 * @see GraggBulirschStoerIntegrator
 * @version $Revision$ $Date$
 * @author E. Hairer and G. Wanner (fortran version)
 * @since 1.2
 */

class GraggBulirschStoerStepInterpolator extends AbstractStepInterpolator {

	/** Slope at the beginning of the step. */
	private double[] y0Dot;

	/** State at the end of the step. */
	private double[] y1;

	/** Slope at the end of the step. */
	private double[] y1Dot;

	/**
	 * Derivatives at the middle of the step. element 0 is state at midpoint,
	 * element 1 is first derivative ...
	 */
	private double[][] yMidDots;

	/** Interpolation polynoms. */
	private double[][] polynoms;

	/** Error coefficients for the interpolation. */
	private double[] errfac;

	/** Degree of the interpolation polynoms. */
	private int currentDegree;

	/**
	 * Reallocate the internal tables. Reallocate the internal tables in order
	 * to be able to handle interpolation polynoms up to the given degree
	 * 
	 * @param maxDegree
	 *            maximal degree to handle
	 */
	private void resetTables(final int maxDegree) {

		if (maxDegree < 0) {
			polynoms = null;
			errfac = null;
			currentDegree = -1;
		} else {

			final double[][] newPols = new double[maxDegree + 1][];
			if (polynoms != null) {
				System.arraycopy(polynoms, 0, newPols, 0, polynoms.length);
				for (int i = polynoms.length; i < newPols.length; ++i) {
					newPols[i] = new double[currentState.length];
				}
			} else {
				for (int i = 0; i < newPols.length; ++i) {
					newPols[i] = new double[currentState.length];
				}
			}
			polynoms = newPols;

			// initialize the error factors array for interpolation
			if (maxDegree <= 4) {
				errfac = null;
			} else {
				errfac = new double[maxDegree - 4];
				for (int i = 0; i < errfac.length; ++i) {
					final int ip5 = i + 5;
					errfac[i] = 1.0 / (ip5 * ip5);
					final double e = 0.5 * Math.sqrt(((double) (i + 1)) / ip5);
					for (int j = 0; j <= i; ++j) {
						errfac[i] *= e / (j + 1);
					}
				}
			}

			currentDegree = 0;

		}

	}

	/**
	 * Simple constructor. This constructor should not be used directly, it is
	 * only intended for the serialization process.
	 */
	public GraggBulirschStoerStepInterpolator() {
		y0Dot = null;
		y1 = null;
		y1Dot = null;
		yMidDots = null;
		resetTables(-1);
	}

	/**
	 * Simple constructor.
	 * 
	 * @param y
	 *            reference to the integrator array holding the current state
	 * @param y0Dot
	 *            reference to the integrator array holding the slope at the
	 *            beginning of the step
	 * @param y1
	 *            reference to the integrator array holding the state at the end
	 *            of the step
	 * @param y1Dot
	 *            reference to the integrator array holding the slope at theend
	 *            of the step
	 * @param yMidDots
	 *            reference to the integrator array holding the derivatives at
	 *            the middle point of the step
	 * @param forward
	 *            integration direction indicator
	 */
	public GraggBulirschStoerStepInterpolator(final double[] y, final double[] y0Dot, final double[] y1,
			final double[] y1Dot, final double[][] yMidDots, final boolean forward) {

		super(y, forward);
		this.y0Dot = y0Dot;
		this.y1 = y1;
		this.y1Dot = y1Dot;
		this.yMidDots = yMidDots;

		resetTables(yMidDots.length + 4);

	}

	/**
	 * Copy constructor.
	 * 
	 * @param interpolator
	 *            interpolator to copy from. The copy is a deep copy: its arrays
	 *            are separated from the original arrays of the instance
	 */
	public GraggBulirschStoerStepInterpolator(final GraggBulirschStoerStepInterpolator interpolator) {

		super(interpolator);

		final int dimension = currentState.length;

		// the interpolator has been finalized,
		// the following arrays are not needed anymore
		y0Dot = null;
		y1 = null;
		y1Dot = null;
		yMidDots = null;

		// copy the interpolation polynoms (up to the current degree only)
		if (interpolator.polynoms == null) {
			polynoms = null;
			currentDegree = -1;
		} else {
			resetTables(interpolator.currentDegree);
			for (int i = 0; i < polynoms.length; ++i) {
				polynoms[i] = new double[dimension];
				System.arraycopy(interpolator.polynoms[i], 0, polynoms[i], 0, dimension);
			}
			currentDegree = interpolator.currentDegree;
		}

	}

	/** {@inheritDoc} */
	@Override
	protected StepInterpolator doCopy() {
		return new GraggBulirschStoerStepInterpolator(this);
	}

	/**
	 * Compute the interpolation coefficients for dense output.
	 * 
	 * @param mu
	 *            degree of the interpolation polynomial
	 * @param h
	 *            current step
	 */
	public void computeCoefficients(final int mu, final double h) {

		if ((polynoms == null) || (polynoms.length <= (mu + 4))) {
			resetTables(mu + 4);
		}

		currentDegree = mu + 4;

		for (int i = 0; i < currentState.length; ++i) {

			final double yp0 = h * y0Dot[i];
			final double yp1 = h * y1Dot[i];
			final double ydiff = y1[i] - currentState[i];
			final double aspl = ydiff - yp1;
			final double bspl = yp0 - ydiff;

			polynoms[0][i] = currentState[i];
			polynoms[1][i] = ydiff;
			polynoms[2][i] = aspl;
			polynoms[3][i] = bspl;

			if (mu < 0) {
				return;
			}

			// compute the remaining coefficients
			final double ph0 = 0.5 * (currentState[i] + y1[i]) + 0.125 * (aspl + bspl);
			polynoms[4][i] = 16 * (yMidDots[0][i] - ph0);

			if (mu > 0) {
				final double ph1 = ydiff + 0.25 * (aspl - bspl);
				polynoms[5][i] = 16 * (yMidDots[1][i] - ph1);

				if (mu > 1) {
					final double ph2 = yp1 - yp0;
					polynoms[6][i] = 16 * (yMidDots[2][i] - ph2 + polynoms[4][i]);

					if (mu > 2) {
						final double ph3 = 6 * (bspl - aspl);
						polynoms[7][i] = 16 * (yMidDots[3][i] - ph3 + 3 * polynoms[5][i]);

						for (int j = 4; j <= mu; ++j) {
							final double fac1 = 0.5 * j * (j - 1);
							final double fac2 = 2 * fac1 * (j - 2) * (j - 3);
							polynoms[j + 4][i] = 16
									* (yMidDots[j][i] + fac1 * polynoms[j + 2][i] - fac2 * polynoms[j][i]);
						}

					}
				}
			}
		}

	}

	/**
	 * Estimate interpolation error.
	 * 
	 * @param scale
	 *            scaling array
	 * @return estimate of the interpolation error
	 */
	public double estimateError(final double[] scale) {
		double error = 0;
		if (currentDegree >= 5) {
			for (int i = 0; i < currentState.length; ++i) {
				final double e = polynoms[currentDegree][i] / scale[i];
				error += e * e;
			}
			error = Math.sqrt(error / currentState.length) * errfac[currentDegree - 5];
		}
		return error;
	}

	/** {@inheritDoc} */
	@Override
	protected void computeInterpolatedStateAndDerivatives(final double theta, final double oneMinusThetaH)
			throws DerivativeException {

		final int dimension = currentState.length;

		final double oneMinusTheta = 1.0 - theta;
		final double theta05 = theta - 0.5;
		final double tOmT = theta * oneMinusTheta;
		final double t4 = tOmT * tOmT;
		final double t4Dot = 2 * tOmT * (1 - 2 * theta);
		final double dot1 = 1.0 / h;
		final double dot2 = theta * (2 - 3 * theta) / h;
		final double dot3 = ((3 * theta - 4) * theta + 1) / h;

		for (int i = 0; i < dimension; ++i) {

			final double p0 = polynoms[0][i];
			final double p1 = polynoms[1][i];
			final double p2 = polynoms[2][i];
			final double p3 = polynoms[3][i];
			interpolatedState[i] = p0 + theta * (p1 + oneMinusTheta * (p2 * theta + p3 * oneMinusTheta));
			interpolatedDerivatives[i] = dot1 * p1 + dot2 * p2 + dot3 * p3;

			if (currentDegree > 3) {
				double cDot = 0;
				double c = polynoms[currentDegree][i];
				for (int j = currentDegree - 1; j > 3; --j) {
					final double d = 1.0 / (j - 3);
					cDot = d * (theta05 * cDot + c);
					c = polynoms[j][i] + c * d * theta05;
				}
				interpolatedState[i] += t4 * c;
				interpolatedDerivatives[i] += (t4 * cDot + t4Dot * c) / h;
			}

		}

		if (h == 0) {
			// in this degenerated case, the previous computation leads to NaN
			// for derivatives
			// we fix this by using the derivatives at midpoint
			System.arraycopy(yMidDots[1], 0, interpolatedDerivatives, 0, dimension);
		}

	}

	/** {@inheritDoc} */
	@Override
	public void writeExternal(final ObjectOutput out) throws IOException {

		final int dimension = (currentState == null) ? -1 : currentState.length;

		// save the state of the base class
		writeBaseExternal(out);

		// save the local attributes (but not the temporary vectors)
		out.writeInt(currentDegree);
		for (int k = 0; k <= currentDegree; ++k) {
			for (int l = 0; l < dimension; ++l) {
				out.writeDouble(polynoms[k][l]);
			}
		}

	}

	/** {@inheritDoc} */
	@Override
	public void readExternal(final ObjectInput in) throws IOException {

		// read the base class
		final double t = readBaseExternal(in);
		final int dimension = (currentState == null) ? -1 : currentState.length;

		// read the local attributes
		final int degree = in.readInt();
		resetTables(degree);
		currentDegree = degree;

		for (int k = 0; k <= currentDegree; ++k) {
			for (int l = 0; l < dimension; ++l) {
				polynoms[k][l] = in.readDouble();
			}
		}

		// we can now set the interpolated time and state
		setInterpolatedTime(t);

	}

	/** Serializable version identifier */
	private static final long serialVersionUID = 7320613236731409847L;

}
