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 java.util.Arrays;
020import org.apache.commons.math3.analysis.BivariateFunction;
021import org.apache.commons.math3.exception.DimensionMismatchException;
022import org.apache.commons.math3.exception.NoDataException;
023import org.apache.commons.math3.exception.OutOfRangeException;
024import org.apache.commons.math3.exception.NonMonotonicSequenceException;
025import org.apache.commons.math3.util.MathArrays;
026
027/**
028 * Function that implements the
029 * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation">
030 * bicubic spline interpolation</a>. Due to numerical accuracy issues this should not
031 * be used.
032 *
033 * @since 2.1
034 * @deprecated as of 3.4 replaced by
035 * {@link org.apache.commons.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolatingFunction}
036 */
037@Deprecated
038public class BicubicSplineInterpolatingFunction
039    implements BivariateFunction {
040    /** Number of coefficients. */
041    private static final int NUM_COEFF = 16;
042    /**
043     * Matrix to compute the spline coefficients from the function values
044     * and function derivatives values
045     */
046    private static final double[][] AINV = {
047        { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
048        { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
049        { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 },
050        { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 },
051        { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
052        { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 },
053        { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 },
054        { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 },
055        { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
056        { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 },
057        { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 },
058        { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 },
059        { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 },
060        { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 },
061        { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 },
062        { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 }
063    };
064
065    /** Samples x-coordinates */
066    private final double[] xval;
067    /** Samples y-coordinates */
068    private final double[] yval;
069    /** Set of cubic splines patching the whole data grid */
070    private final BicubicSplineFunction[][] splines;
071    /**
072     * Partial derivatives.
073     * The value of the first index determines the kind of derivatives:
074     * 0 = first partial derivatives wrt x
075     * 1 = first partial derivatives wrt y
076     * 2 = second partial derivatives wrt x
077     * 3 = second partial derivatives wrt y
078     * 4 = cross partial derivatives
079     */
080    private final BivariateFunction[][][] partialDerivatives;
081
082    /**
083     * @param x Sample values of the x-coordinate, in increasing order.
084     * @param y Sample values of the y-coordinate, in increasing order.
085     * @param f Values of the function on every grid point.
086     * @param dFdX Values of the partial derivative of function with respect
087     * to x on every grid point.
088     * @param dFdY Values of the partial derivative of function with respect
089     * to y on every grid point.
090     * @param d2FdXdY Values of the cross partial derivative of function on
091     * every grid point.
092     * @throws DimensionMismatchException if the various arrays do not contain
093     * the expected number of elements.
094     * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
095     * not strictly increasing.
096     * @throws NoDataException if any of the arrays has zero length.
097     */
098    public BicubicSplineInterpolatingFunction(double[] x,
099                                              double[] y,
100                                              double[][] f,
101                                              double[][] dFdX,
102                                              double[][] dFdY,
103                                              double[][] d2FdXdY)
104        throws DimensionMismatchException,
105               NoDataException,
106               NonMonotonicSequenceException {
107        this(x, y, f, dFdX, dFdY, d2FdXdY, false);
108    }
109
110    /**
111     * @param x Sample values of the x-coordinate, in increasing order.
112     * @param y Sample values of the y-coordinate, in increasing order.
113     * @param f Values of the function on every grid point.
114     * @param dFdX Values of the partial derivative of function with respect
115     * to x on every grid point.
116     * @param dFdY Values of the partial derivative of function with respect
117     * to y on every grid point.
118     * @param d2FdXdY Values of the cross partial derivative of function on
119     * every grid point.
120     * @param initializeDerivatives Whether to initialize the internal data
121     * needed for calling any of the methods that compute the partial derivatives
122     * this function.
123     * @throws DimensionMismatchException if the various arrays do not contain
124     * the expected number of elements.
125     * @throws NonMonotonicSequenceException if {@code x} or {@code y} are
126     * not strictly increasing.
127     * @throws NoDataException if any of the arrays has zero length.
128     *
129     * @see #partialDerivativeX(double,double)
130     * @see #partialDerivativeY(double,double)
131     * @see #partialDerivativeXX(double,double)
132     * @see #partialDerivativeYY(double,double)
133     * @see #partialDerivativeXY(double,double)
134     */
135    public BicubicSplineInterpolatingFunction(double[] x,
136                                              double[] y,
137                                              double[][] f,
138                                              double[][] dFdX,
139                                              double[][] dFdY,
140                                              double[][] d2FdXdY,
141                                              boolean initializeDerivatives)
142        throws DimensionMismatchException,
143               NoDataException,
144               NonMonotonicSequenceException {
145        final int xLen = x.length;
146        final int yLen = y.length;
147
148        if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
149            throw new NoDataException();
150        }
151        if (xLen != f.length) {
152            throw new DimensionMismatchException(xLen, f.length);
153        }
154        if (xLen != dFdX.length) {
155            throw new DimensionMismatchException(xLen, dFdX.length);
156        }
157        if (xLen != dFdY.length) {
158            throw new DimensionMismatchException(xLen, dFdY.length);
159        }
160        if (xLen != d2FdXdY.length) {
161            throw new DimensionMismatchException(xLen, d2FdXdY.length);
162        }
163
164        MathArrays.checkOrder(x);
165        MathArrays.checkOrder(y);
166
167        xval = x.clone();
168        yval = y.clone();
169
170        final int lastI = xLen - 1;
171        final int lastJ = yLen - 1;
172        splines = new BicubicSplineFunction[lastI][lastJ];
173
174        for (int i = 0; i < lastI; i++) {
175            if (f[i].length != yLen) {
176                throw new DimensionMismatchException(f[i].length, yLen);
177            }
178            if (dFdX[i].length != yLen) {
179                throw new DimensionMismatchException(dFdX[i].length, yLen);
180            }
181            if (dFdY[i].length != yLen) {
182                throw new DimensionMismatchException(dFdY[i].length, yLen);
183            }
184            if (d2FdXdY[i].length != yLen) {
185                throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
186            }
187            final int ip1 = i + 1;
188            for (int j = 0; j < lastJ; j++) {
189                final int jp1 = j + 1;
190                final double[] beta = new double[] {
191                    f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
192                    dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1],
193                    dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1],
194                    d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1]
195                };
196
197                splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta),
198                                                          initializeDerivatives);
199            }
200        }
201
202        if (initializeDerivatives) {
203            // Compute all partial derivatives.
204            partialDerivatives = new BivariateFunction[5][lastI][lastJ];
205
206            for (int i = 0; i < lastI; i++) {
207                for (int j = 0; j < lastJ; j++) {
208                    final BicubicSplineFunction bcs = splines[i][j];
209                    partialDerivatives[0][i][j] = bcs.partialDerivativeX();
210                    partialDerivatives[1][i][j] = bcs.partialDerivativeY();
211                    partialDerivatives[2][i][j] = bcs.partialDerivativeXX();
212                    partialDerivatives[3][i][j] = bcs.partialDerivativeYY();
213                    partialDerivatives[4][i][j] = bcs.partialDerivativeXY();
214                }
215            }
216        } else {
217            // Partial derivative methods cannot be used.
218            partialDerivatives = null;
219        }
220    }
221
222    /**
223     * {@inheritDoc}
224     */
225    public double value(double x, double y)
226        throws OutOfRangeException {
227        final int i = searchIndex(x, xval);
228        final int j = searchIndex(y, yval);
229
230        final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
231        final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
232
233        return splines[i][j].value(xN, yN);
234    }
235
236    /**
237     * Indicates whether a point is within the interpolation range.
238     *
239     * @param x First coordinate.
240     * @param y Second coordinate.
241     * @return {@code true} if (x, y) is a valid point.
242     * @since 3.3
243     */
244    public boolean isValidPoint(double x, double y) {
245        if (x < xval[0] ||
246            x > xval[xval.length - 1] ||
247            y < yval[0] ||
248            y > yval[yval.length - 1]) {
249            return false;
250        } else {
251            return true;
252        }
253    }
254
255    /**
256     * @param x x-coordinate.
257     * @param y y-coordinate.
258     * @return the value at point (x, y) of the first partial derivative with
259     * respect to x.
260     * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
261     * the range defined by the boundary values of {@code xval} (resp.
262     * {@code yval}).
263     * @throws NullPointerException if the internal data were not initialized
264     * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
265     *             double[][],double[][],double[][],boolean) constructor}).
266     */
267    public double partialDerivativeX(double x, double y)
268        throws OutOfRangeException {
269        return partialDerivative(0, x, y);
270    }
271    /**
272     * @param x x-coordinate.
273     * @param y y-coordinate.
274     * @return the value at point (x, y) of the first partial derivative with
275     * respect to y.
276     * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
277     * the range defined by the boundary values of {@code xval} (resp.
278     * {@code yval}).
279     * @throws NullPointerException if the internal data were not initialized
280     * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
281     *             double[][],double[][],double[][],boolean) constructor}).
282     */
283    public double partialDerivativeY(double x, double y)
284        throws OutOfRangeException {
285        return partialDerivative(1, x, y);
286    }
287    /**
288     * @param x x-coordinate.
289     * @param y y-coordinate.
290     * @return the value at point (x, y) of the second partial derivative with
291     * respect to x.
292     * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
293     * the range defined by the boundary values of {@code xval} (resp.
294     * {@code yval}).
295     * @throws NullPointerException if the internal data were not initialized
296     * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
297     *             double[][],double[][],double[][],boolean) constructor}).
298     */
299    public double partialDerivativeXX(double x, double y)
300        throws OutOfRangeException {
301        return partialDerivative(2, x, y);
302    }
303    /**
304     * @param x x-coordinate.
305     * @param y y-coordinate.
306     * @return the value at point (x, y) of the second partial derivative with
307     * respect to y.
308     * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
309     * the range defined by the boundary values of {@code xval} (resp.
310     * {@code yval}).
311     * @throws NullPointerException if the internal data were not initialized
312     * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
313     *             double[][],double[][],double[][],boolean) constructor}).
314     */
315    public double partialDerivativeYY(double x, double y)
316        throws OutOfRangeException {
317        return partialDerivative(3, x, y);
318    }
319    /**
320     * @param x x-coordinate.
321     * @param y y-coordinate.
322     * @return the value at point (x, y) of the second partial cross-derivative.
323     * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
324     * the range defined by the boundary values of {@code xval} (resp.
325     * {@code yval}).
326     * @throws NullPointerException if the internal data were not initialized
327     * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
328     *             double[][],double[][],double[][],boolean) constructor}).
329     */
330    public double partialDerivativeXY(double x, double y)
331        throws OutOfRangeException {
332        return partialDerivative(4, x, y);
333    }
334
335    /**
336     * @param which First index in {@link #partialDerivatives}.
337     * @param x x-coordinate.
338     * @param y y-coordinate.
339     * @return the value at point (x, y) of the selected partial derivative.
340     * @throws OutOfRangeException if {@code x} (resp. {@code y}) is outside
341     * the range defined by the boundary values of {@code xval} (resp.
342     * {@code yval}).
343     * @throws NullPointerException if the internal data were not initialized
344     * (cf. {@link #BicubicSplineInterpolatingFunction(double[],double[],double[][],
345     *             double[][],double[][],double[][],boolean) constructor}).
346     */
347    private double partialDerivative(int which, double x, double y)
348        throws OutOfRangeException {
349        final int i = searchIndex(x, xval);
350        final int j = searchIndex(y, yval);
351
352        final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
353        final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
354
355        return partialDerivatives[which][i][j].value(xN, yN);
356    }
357
358    /**
359     * @param c Coordinate.
360     * @param val Coordinate samples.
361     * @return the index in {@code val} corresponding to the interval
362     * containing {@code c}.
363     * @throws OutOfRangeException if {@code c} is out of the
364     * range defined by the boundary values of {@code val}.
365     */
366    private int searchIndex(double c, double[] val) {
367        final int r = Arrays.binarySearch(val, c);
368
369        if (r == -1 ||
370            r == -val.length - 1) {
371            throw new OutOfRangeException(c, val[0], val[val.length - 1]);
372        }
373
374        if (r < 0) {
375            // "c" in within an interpolation sub-interval: Return the
376            // index of the sample at the lower end of the sub-interval.
377            return -r - 2;
378        }
379        final int last = val.length - 1;
380        if (r == last) {
381            // "c" is the last sample of the range: Return the index
382            // of the sample at the lower end of the last sub-interval.
383            return last - 1;
384        }
385
386        // "c" is another sample point.
387        return r;
388    }
389
390    /**
391     * Compute the spline coefficients from the list of function values and
392     * function partial derivatives values at the four corners of a grid
393     * element. They must be specified in the following order:
394     * <ul>
395     *  <li>f(0,0)</li>
396     *  <li>f(1,0)</li>
397     *  <li>f(0,1)</li>
398     *  <li>f(1,1)</li>
399     *  <li>f<sub>x</sub>(0,0)</li>
400     *  <li>f<sub>x</sub>(1,0)</li>
401     *  <li>f<sub>x</sub>(0,1)</li>
402     *  <li>f<sub>x</sub>(1,1)</li>
403     *  <li>f<sub>y</sub>(0,0)</li>
404     *  <li>f<sub>y</sub>(1,0)</li>
405     *  <li>f<sub>y</sub>(0,1)</li>
406     *  <li>f<sub>y</sub>(1,1)</li>
407     *  <li>f<sub>xy</sub>(0,0)</li>
408     *  <li>f<sub>xy</sub>(1,0)</li>
409     *  <li>f<sub>xy</sub>(0,1)</li>
410     *  <li>f<sub>xy</sub>(1,1)</li>
411     * </ul>
412     * where the subscripts indicate the partial derivative with respect to
413     * the corresponding variable(s).
414     *
415     * @param beta List of function values and function partial derivatives
416     * values.
417     * @return the spline coefficients.
418     */
419    private double[] computeSplineCoefficients(double[] beta) {
420        final double[] a = new double[NUM_COEFF];
421
422        for (int i = 0; i < NUM_COEFF; i++) {
423            double result = 0;
424            final double[] row = AINV[i];
425            for (int j = 0; j < NUM_COEFF; j++) {
426                result += row[j] * beta[j];
427            }
428            a[i] = result;
429        }
430
431        return a;
432    }
433}
434
435/**
436 * 2D-spline function.
437 *
438 */
439class BicubicSplineFunction implements BivariateFunction {
440    /** Number of points. */
441    private static final short N = 4;
442    /** Coefficients */
443    private final double[][] a;
444    /** First partial derivative along x. */
445    private final BivariateFunction partialDerivativeX;
446    /** First partial derivative along y. */
447    private final BivariateFunction partialDerivativeY;
448    /** Second partial derivative along x. */
449    private final BivariateFunction partialDerivativeXX;
450    /** Second partial derivative along y. */
451    private final BivariateFunction partialDerivativeYY;
452    /** Second crossed partial derivative. */
453    private final BivariateFunction partialDerivativeXY;
454
455    /**
456     * Simple constructor.
457     *
458     * @param coeff Spline coefficients.
459     */
460    BicubicSplineFunction(double[] coeff) {
461        this(coeff, false);
462    }
463
464    /**
465     * Simple constructor.
466     *
467     * @param coeff Spline coefficients.
468     * @param initializeDerivatives Whether to initialize the internal data
469     * needed for calling any of the methods that compute the partial derivatives
470     * this function.
471     */
472    BicubicSplineFunction(double[] coeff, boolean initializeDerivatives) {
473        a = new double[N][N];
474        for (int i = 0; i < N; i++) {
475            for (int j = 0; j < N; j++) {
476                a[i][j] = coeff[i * N + j];
477            }
478        }
479
480        if (initializeDerivatives) {
481            // Compute all partial derivatives functions.
482            final double[][] aX = new double[N][N];
483            final double[][] aY = new double[N][N];
484            final double[][] aXX = new double[N][N];
485            final double[][] aYY = new double[N][N];
486            final double[][] aXY = new double[N][N];
487
488            for (int i = 0; i < N; i++) {
489                for (int j = 0; j < N; j++) {
490                    final double c = a[i][j];
491                    aX[i][j] = i * c;
492                    aY[i][j] = j * c;
493                    aXX[i][j] = (i - 1) * aX[i][j];
494                    aYY[i][j] = (j - 1) * aY[i][j];
495                    aXY[i][j] = j * aX[i][j];
496                }
497            }
498
499            partialDerivativeX = new BivariateFunction() {
500                    /** {@inheritDoc} */
501                    public double value(double x, double y)  {
502                        final double x2 = x * x;
503                        final double[] pX = {0, 1, x, x2};
504
505                        final double y2 = y * y;
506                        final double y3 = y2 * y;
507                        final double[] pY = {1, y, y2, y3};
508
509                        return apply(pX, pY, aX);
510                    }
511                };
512            partialDerivativeY = new BivariateFunction() {
513                    /** {@inheritDoc} */
514                    public double value(double x, double y)  {
515                        final double x2 = x * x;
516                        final double x3 = x2 * x;
517                        final double[] pX = {1, x, x2, x3};
518
519                        final double y2 = y * y;
520                        final double[] pY = {0, 1, y, y2};
521
522                        return apply(pX, pY, aY);
523                    }
524                };
525            partialDerivativeXX = new BivariateFunction() {
526                    /** {@inheritDoc} */
527                    public double value(double x, double y)  {
528                        final double[] pX = {0, 0, 1, x};
529
530                        final double y2 = y * y;
531                        final double y3 = y2 * y;
532                        final double[] pY = {1, y, y2, y3};
533
534                        return apply(pX, pY, aXX);
535                    }
536                };
537            partialDerivativeYY = new BivariateFunction() {
538                    /** {@inheritDoc} */
539                    public double value(double x, double y)  {
540                        final double x2 = x * x;
541                        final double x3 = x2 * x;
542                        final double[] pX = {1, x, x2, x3};
543
544                        final double[] pY = {0, 0, 1, y};
545
546                        return apply(pX, pY, aYY);
547                    }
548                };
549            partialDerivativeXY = new BivariateFunction() {
550                    /** {@inheritDoc} */
551                    public double value(double x, double y)  {
552                        final double x2 = x * x;
553                        final double[] pX = {0, 1, x, x2};
554
555                        final double y2 = y * y;
556                        final double[] pY = {0, 1, y, y2};
557
558                        return apply(pX, pY, aXY);
559                    }
560                };
561        } else {
562            partialDerivativeX = null;
563            partialDerivativeY = null;
564            partialDerivativeXX = null;
565            partialDerivativeYY = null;
566            partialDerivativeXY = null;
567        }
568    }
569
570    /**
571     * {@inheritDoc}
572     */
573    public double value(double x, double y) {
574        if (x < 0 || x > 1) {
575            throw new OutOfRangeException(x, 0, 1);
576        }
577        if (y < 0 || y > 1) {
578            throw new OutOfRangeException(y, 0, 1);
579        }
580
581        final double x2 = x * x;
582        final double x3 = x2 * x;
583        final double[] pX = {1, x, x2, x3};
584
585        final double y2 = y * y;
586        final double y3 = y2 * y;
587        final double[] pY = {1, y, y2, y3};
588
589        return apply(pX, pY, a);
590    }
591
592    /**
593     * Compute the value of the bicubic polynomial.
594     *
595     * @param pX Powers of the x-coordinate.
596     * @param pY Powers of the y-coordinate.
597     * @param coeff Spline coefficients.
598     * @return the interpolated value.
599     */
600    private double apply(double[] pX, double[] pY, double[][] coeff) {
601        double result = 0;
602        for (int i = 0; i < N; i++) {
603            for (int j = 0; j < N; j++) {
604                result += coeff[i][j] * pX[i] * pY[j];
605            }
606        }
607
608        return result;
609    }
610
611    /**
612     * @return the partial derivative wrt {@code x}.
613     */
614    public BivariateFunction partialDerivativeX() {
615        return partialDerivativeX;
616    }
617    /**
618     * @return the partial derivative wrt {@code y}.
619     */
620    public BivariateFunction partialDerivativeY() {
621        return partialDerivativeY;
622    }
623    /**
624     * @return the second partial derivative wrt {@code x}.
625     */
626    public BivariateFunction partialDerivativeXX() {
627        return partialDerivativeXX;
628    }
629    /**
630     * @return the second partial derivative wrt {@code y}.
631     */
632    public BivariateFunction partialDerivativeYY() {
633        return partialDerivativeYY;
634    }
635    /**
636     * @return the second partial cross-derivative.
637     */
638    public BivariateFunction partialDerivativeXY() {
639        return partialDerivativeXY;
640    }
641}