001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.math3.analysis.interpolation;
018
019import org.apache.commons.math3.analysis.UnivariateFunction;
020import org.apache.commons.math3.analysis.polynomials.PolynomialSplineFunction;
021import org.apache.commons.math3.exception.DimensionMismatchException;
022import org.apache.commons.math3.exception.NoDataException;
023import org.apache.commons.math3.exception.NonMonotonicSequenceException;
024import org.apache.commons.math3.exception.NumberIsTooSmallException;
025import org.apache.commons.math3.util.MathArrays;
026
027/**
028 * Generates a bicubic interpolating function. Due to numerical accuracy issues this should not
029 * be used.
030 *
031 * @since 2.2
032 * @deprecated as of 3.4 replaced by {@link org.apache.commons.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolator}
033 */
034@Deprecated
035public class BicubicSplineInterpolator
036    implements BivariateGridInterpolator {
037    /** Whether to initialize internal data used to compute the analytical
038        derivatives of the splines. */
039    private final boolean initializeDerivatives;
040
041    /**
042     * Default constructor.
043     * The argument {@link #BicubicSplineInterpolator(boolean) initializeDerivatives}
044     * is set to {@code false}.
045     */
046    public BicubicSplineInterpolator() {
047        this(false);
048    }
049
050    /**
051     * Creates an interpolator.
052     *
053     * @param initializeDerivatives Whether to initialize the internal data
054     * needed for calling any of the methods that compute the partial derivatives
055     * of the {@link BicubicSplineInterpolatingFunction function} returned from
056     * the call to {@link #interpolate(double[],double[],double[][]) interpolate}.
057     */
058    public BicubicSplineInterpolator(boolean initializeDerivatives) {
059        this.initializeDerivatives = initializeDerivatives;
060    }
061
062    /**
063     * {@inheritDoc}
064     */
065    public BicubicSplineInterpolatingFunction interpolate(final double[] xval,
066                                                          final double[] yval,
067                                                          final double[][] fval)
068        throws NoDataException, DimensionMismatchException,
069               NonMonotonicSequenceException, NumberIsTooSmallException {
070        if (xval.length == 0 || yval.length == 0 || fval.length == 0) {
071            throw new NoDataException();
072        }
073        if (xval.length != fval.length) {
074            throw new DimensionMismatchException(xval.length, fval.length);
075        }
076
077        MathArrays.checkOrder(xval);
078        MathArrays.checkOrder(yval);
079
080        final int xLen = xval.length;
081        final int yLen = yval.length;
082
083        // Samples (first index is y-coordinate, i.e. subarray variable is x)
084        // 0 <= i < xval.length
085        // 0 <= j < yval.length
086        // fX[j][i] = f(xval[i], yval[j])
087        final double[][] fX = new double[yLen][xLen];
088        for (int i = 0; i < xLen; i++) {
089            if (fval[i].length != yLen) {
090                throw new DimensionMismatchException(fval[i].length, yLen);
091            }
092
093            for (int j = 0; j < yLen; j++) {
094                fX[j][i] = fval[i][j];
095            }
096        }
097
098        final SplineInterpolator spInterpolator = new SplineInterpolator();
099
100        // For each line y[j] (0 <= j < yLen), construct a 1D spline with
101        // respect to variable x
102        final PolynomialSplineFunction[] ySplineX = new PolynomialSplineFunction[yLen];
103        for (int j = 0; j < yLen; j++) {
104            ySplineX[j] = spInterpolator.interpolate(xval, fX[j]);
105        }
106
107        // For each line x[i] (0 <= i < xLen), construct a 1D spline with
108        // respect to variable y generated by array fY_1[i]
109        final PolynomialSplineFunction[] xSplineY = new PolynomialSplineFunction[xLen];
110        for (int i = 0; i < xLen; i++) {
111            xSplineY[i] = spInterpolator.interpolate(yval, fval[i]);
112        }
113
114        // Partial derivatives with respect to x at the grid knots
115        final double[][] dFdX = new double[xLen][yLen];
116        for (int j = 0; j < yLen; j++) {
117            final UnivariateFunction f = ySplineX[j].derivative();
118            for (int i = 0; i < xLen; i++) {
119                dFdX[i][j] = f.value(xval[i]);
120            }
121        }
122
123        // Partial derivatives with respect to y at the grid knots
124        final double[][] dFdY = new double[xLen][yLen];
125        for (int i = 0; i < xLen; i++) {
126            final UnivariateFunction f = xSplineY[i].derivative();
127            for (int j = 0; j < yLen; j++) {
128                dFdY[i][j] = f.value(yval[j]);
129            }
130        }
131
132        // Cross partial derivatives
133        final double[][] d2FdXdY = new double[xLen][yLen];
134        for (int i = 0; i < xLen ; i++) {
135            final int nI = nextIndex(i, xLen);
136            final int pI = previousIndex(i);
137            for (int j = 0; j < yLen; j++) {
138                final int nJ = nextIndex(j, yLen);
139                final int pJ = previousIndex(j);
140                d2FdXdY[i][j] = (fval[nI][nJ] - fval[nI][pJ] -
141                                 fval[pI][nJ] + fval[pI][pJ]) /
142                    ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]));
143            }
144        }
145
146        // Create the interpolating splines
147        return new BicubicSplineInterpolatingFunction(xval, yval, fval,
148                                                      dFdX, dFdY, d2FdXdY,
149                                                      initializeDerivatives);
150    }
151
152    /**
153     * Computes the next index of an array, clipping if necessary.
154     * It is assumed (but not checked) that {@code i >= 0}.
155     *
156     * @param i Index.
157     * @param max Upper limit of the array.
158     * @return the next index.
159     */
160    private int nextIndex(int i, int max) {
161        final int index = i + 1;
162        return index < max ? index : index - 1;
163    }
164    /**
165     * Computes the previous index of an array, clipping if necessary.
166     * It is assumed (but not checked) that {@code i} is smaller than the size
167     * of the array.
168     *
169     * @param i Index.
170     * @return the previous index.
171     */
172    private int previousIndex(int i) {
173        final int index = i - 1;
174        return index >= 0 ? index : 0;
175    }
176}