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