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

import org.apache.commons.math.DimensionMismatchException;
import org.apache.commons.math.MathRuntimeException;
import org.apache.commons.math.MathException;
import org.apache.commons.math.util.MathUtils;
import org.apache.commons.math.util.MathUtils.OrderDirection;
import org.apache.commons.math.analysis.BivariateRealFunction;
import org.apache.commons.math.analysis.UnivariateRealFunction;
import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
import org.apache.commons.math.exception.util.LocalizedFormats;

/**
 * Generates a bicubic interpolation function.
 * Before interpolating, smoothing of the input data is performed using
 * splines.
 * See <b>Handbook on splines for the user</b>, ISBN 084939404X,
 * chapter 2.
 *
 * @version $Revision: 1059400 $ $Date: 2011-01-15 20:35:27 +0100 (sam. 15 janv. 2011) $
 * @since 2.1
 * @deprecated This class does not perform smoothing; the name is thus misleading.
 * Please use {@link org.apache.commons.math.analysis.interpolation.BicubicSplineInterpolator}
 * instead. If smoothing is desired, a tentative implementation is provided in class
 * {@link org.apache.commons.math.analysis.interpolation.SmoothingPolynomialBicubicSplineInterpolator}.
 * This class will be removed in math 3.0.
 */
@Deprecated
public class SmoothingBicubicSplineInterpolator
    implements BivariateRealGridInterpolator {
    /**
     * {@inheritDoc}
     */
    public BivariateRealFunction interpolate(final double[] xval,
                                                          final double[] yval,
                                                          final double[][] zval)
        throws MathException, IllegalArgumentException {
        if (xval.length == 0 || yval.length == 0 || zval.length == 0) {
            throw MathRuntimeException.createIllegalArgumentException(LocalizedFormats.NO_DATA);
        }
        if (xval.length != zval.length) {
            throw new DimensionMismatchException(xval.length, zval.length);
        }

        MathUtils.checkOrder(xval, OrderDirection.INCREASING, true);
        MathUtils.checkOrder(yval, OrderDirection.INCREASING, true);

        final int xLen = xval.length;
        final int yLen = yval.length;

        // Samples (first index is y-coordinate, i.e. subarray variable is x)
        // 0 <= i < xval.length
        // 0 <= j < yval.length
        // zX[j][i] = f(xval[i], yval[j])
        final double[][] zX = new double[yLen][xLen];
        for (int i = 0; i < xLen; i++) {
            if (zval[i].length != yLen) {
                throw new DimensionMismatchException(zval[i].length, yLen);
            }

            for (int j = 0; j < yLen; j++) {
                zX[j][i] = zval[i][j];
            }
        }

        final SplineInterpolator spInterpolator = new SplineInterpolator();

        // For each line y[j] (0 <= j < yLen), construct a 1D spline with
        // respect to variable x
        final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen];
        for (int j = 0; j < yLen; j++) {
            ySplineX[j] = spInterpolator.interpolate(xval, zX[j]);
        }

        // For every knot (xval[i], yval[j]) of the grid, calculate corrected
        // values zY_1
        final double[][] zY_1 = new double[xLen][yLen];
        for (int j = 0; j < yLen; j++) {
            final PolynomialSplineFunction f = ySplineX[j];
            for (int i = 0; i < xLen; i++) {
                zY_1[i][j] = f.value(xval[i]);
            }
        }

        // For each line x[i] (0 <= i < xLen), construct a 1D spline with
        // respect to variable y generated by array zY_1[i]
        final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen];
        for (int i = 0; i < xLen; i++) {
            xSplineY[i] = spInterpolator.interpolate(yval, zY_1[i]);
        }

        // For every knot (xval[i], yval[j]) of the grid, calculate corrected
        // values zY_2
        final double[][] zY_2 = new double[xLen][yLen];
        for (int i = 0; i < xLen; i++) {
            final PolynomialSplineFunction f = xSplineY[i];
            for (int j = 0; j < yLen; j++) {
                zY_2[i][j] = f.value(yval[j]);
            }
        }

        // Partial derivatives with respect to x at the grid knots
        final double[][] dZdX = new double[xLen][yLen];
        for (int j = 0; j < yLen; j++) {
            final UnivariateRealFunction f = ySplineX[j].derivative();
            for (int i = 0; i < xLen; i++) {
                dZdX[i][j] = f.value(xval[i]);
            }
        }

        // Partial derivatives with respect to y at the grid knots
        final double[][] dZdY = new double[xLen][yLen];
        for (int i = 0; i < xLen; i++) {
            final UnivariateRealFunction f = xSplineY[i].derivative();
            for (int j = 0; j < yLen; j++) {
                dZdY[i][j] = f.value(yval[j]);
            }
        }

        // Cross partial derivatives
        final double[][] dZdXdY = new double[xLen][yLen];
        for (int i = 0; i < xLen ; i++) {
            final int nI = nextIndex(i, xLen);
            final int pI = previousIndex(i);
            for (int j = 0; j < yLen; j++) {
                final int nJ = nextIndex(j, yLen);
                final int pJ = previousIndex(j);
                dZdXdY[i][j] = (zY_2[nI][nJ] - zY_2[nI][pJ] -
                                zY_2[pI][nJ] + zY_2[pI][pJ]) /
                    ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]));
            }
        }

        // Create the interpolating splines
        return new BicubicSplineInterpolatingFunction(xval, yval, zY_2,
                                                      dZdX, dZdY, dZdXdY);
    }

    /**
     * Compute the next index of an array, clipping if necessary.
     * It is assumed (but not checked) that {@code i} is larger than or equal to 0}.
     *
     * @param i Index
     * @param max Upper limit of the array
     * @return the next index
     */
    private int nextIndex(int i, int max) {
        final int index = i + 1;
        return index < max ? index : index - 1;
    }
    /**
     * Compute the previous index of an array, clipping if necessary.
     * It is assumed (but not checked) that {@code i} is smaller than the size of the array.
     *
     * @param i Index
     * @return the previous index
     */
    private int previousIndex(int i) {
        final int index = i - 1;
        return index >= 0 ? index : 0;
    }
}