Example usage for org.apache.commons.math3.analysis.interpolation BicubicSplineInterpolatingFunction partialDerivativeX

List of usage examples for org.apache.commons.math3.analysis.interpolation BicubicSplineInterpolatingFunction partialDerivativeX

Introduction

In this page you can find the example usage for org.apache.commons.math3.analysis.interpolation BicubicSplineInterpolatingFunction partialDerivativeX.

Prototype

BivariateFunction partialDerivativeX

To view the source code for org.apache.commons.math3.analysis.interpolation BicubicSplineInterpolatingFunction partialDerivativeX.

Click Source Link

Document

First partial derivative along x.

Usage

From source file:au.gov.ga.conn4d.utils.TricubicSplineInterpolator.java

/**
 * {@inheritDoc}/*w ww.  j  ava 2s.c  om*/
 */
public TricubicSplineInterpolatingFunction interpolate(final double[] zval, final double[] yval,
        final double[] xval, final double[][][] fval)
        throws NoDataException, DimensionMismatchException, NonMonotonicSequenceException {
    if (zval.length == 0 || yval.length == 0 || xval.length == 0 || fval.length == 0) {
        throw new NoDataException();
    }
    if (zval.length != fval.length) {
        throw new DimensionMismatchException(zval.length, fval.length);
    }

    MathArrays.checkOrder(zval);
    MathArrays.checkOrder(yval);

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

    final double[] zzval;
    final double[][][] ffval = new double[zLen][][];

    if (zval[0] > zval[zLen - 1]) {
        zzval = VectorMath.selectionSort(zval);
        for (int i = 0; i < zLen; i++) {
            ffval[i] = fval[zLen - 1 - i];
        }
    } else {
        zzval = zval;
    }

    // Samples, re-ordered as (z, x, y) and (y, z, x) tuplets
    // fvalXY[k][i][j] = f(xval[i], yval[j], zval[k])
    // fvalZX[j][k][i] = f(xval[i], yval[j], zval[k])
    final double[][][] fvalXY = new double[xLen][zLen][yLen];
    final double[][][] fvalZX = new double[yLen][xLen][zLen];
    for (int i = 0; i < zLen; i++) {
        if (ffval[i].length != yLen) {
            throw new DimensionMismatchException(ffval[i].length, yLen);
        }

        for (int j = 0; j < yLen; j++) {
            if (ffval[i][j].length != xLen) {
                throw new DimensionMismatchException(ffval[i][j].length, xLen);
            }

            for (int k = 0; k < xLen; k++) {
                final double v = ffval[i][j][k];
                fvalXY[k][i][j] = v;
                fvalZX[j][k][i] = v;
            }
        }
    }

    final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator();

    // For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z
    final BicubicSplineInterpolatingFunction[] xSplineYZ = new BicubicSplineInterpolatingFunction[zLen];
    for (int i = 0; i < zLen; i++) {
        xSplineYZ[i] = bsi.interpolate(yval, xval, ffval[i]);
    }

    // For each line y[j] (0 <= j < yLen), construct a 2D spline in z and x
    final BicubicSplineInterpolatingFunction[] ySplineZX = new BicubicSplineInterpolatingFunction[yLen];
    for (int j = 0; j < yLen; j++) {
        ySplineZX[j] = bsi.interpolate(xval, zzval, fvalZX[j]);
    }

    // For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y
    final BicubicSplineInterpolatingFunction[] zSplineXY = new BicubicSplineInterpolatingFunction[xLen];
    for (int k = 0; k < xLen; k++) {
        zSplineXY[k] = bsi.interpolate(zzval, yval, fvalXY[k]);
    }

    // Partial derivatives wrt x and wrt y
    final double[][][] dFdX = new double[zLen][yLen][xLen];
    final double[][][] dFdY = new double[zLen][yLen][xLen];
    final double[][][] d2FdXdY = new double[zLen][yLen][xLen];
    for (int k = 0; k < xLen; k++) {
        final BicubicSplineInterpolatingFunction f = zSplineXY[k];
        for (int i = 0; i < zLen; i++) {
            final double z = zzval[i];
            for (int j = 0; j < yLen; j++) {
                final double y = yval[j];
                dFdX[i][j][k] = f.partialDerivativeX(z, y);
                dFdY[i][j][k] = f.partialDerivativeY(z, y);
                d2FdXdY[i][j][k] = f.partialDerivativeXY(z, y);
            }
        }
    }

    // Partial derivatives wrt y and wrt z
    final double[][][] dFdZ = new double[zLen][yLen][xLen];
    final double[][][] d2FdYdZ = new double[zLen][yLen][xLen];
    for (int i = 0; i < zLen; i++) {
        final BicubicSplineInterpolatingFunction f = xSplineYZ[i];
        for (int j = 0; j < yLen; j++) {
            final double y = yval[j];
            for (int k = 0; k < xLen; k++) {
                final double x = xval[k];
                dFdZ[i][j][k] = f.partialDerivativeY(y, x);
                d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, x);
            }
        }
    }

    // Partial derivatives wrt x and wrt z
    final double[][][] d2FdZdX = new double[zLen][yLen][xLen];
    for (int j = 0; j < yLen; j++) {
        final BicubicSplineInterpolatingFunction f = ySplineZX[j];
        for (int k = 0; k < xLen; k++) {
            final double x = xval[k];
            for (int i = 0; i < zLen; i++) {
                final double z = zzval[i];
                d2FdZdX[i][j][k] = f.partialDerivativeXY(x, z);
            }
        }
    }

    // Third partial cross-derivatives
    final double[][][] d3FdXdYdZ = new double[zLen][yLen][xLen];
    for (int i = 0; i < zLen; i++) {
        final int nI = nextIndex(i, zLen);
        final int pI = previousIndex(i);
        for (int j = 0; j < yLen; j++) {
            final int nJ = nextIndex(j, yLen);
            final int pJ = previousIndex(j);
            for (int k = 0; k < xLen; k++) {
                final int nK = nextIndex(k, xLen);
                final int pK = previousIndex(k);

                // XXX Not sure about this formula
                d3FdXdYdZ[i][j][k] = (ffval[nI][nJ][nK] - ffval[nI][pJ][nK] - ffval[pI][nJ][nK]
                        + ffval[pI][pJ][nK] - ffval[nI][nJ][pK] + ffval[nI][pJ][pK] + ffval[pI][nJ][pK]
                        - ffval[pI][pJ][pK])
                        / ((zzval[nI] - zzval[pI]) * (yval[nJ] - yval[pJ]) * (xval[nK] - xval[pK]));
            }
        }
    }

    // Create the interpolating splines
    return new TricubicSplineInterpolatingFunction(zzval, yval, xval, ffval, dFdX, dFdY, dFdZ, d2FdXdY, d2FdZdX,
            d2FdYdZ, d3FdXdYdZ);
}

From source file:au.gov.ga.conn4d.utils.TricubicSplineInterpolator.java

public TricubicSplineInterpolatingFunction interpolate(final double[] zval, final double[] yval,
        final double[] xval, final float[][][] fval)
        throws NoDataException, DimensionMismatchException, NonMonotonicSequenceException {
    if (zval.length == 0 || yval.length == 0 || xval.length == 0 || fval.length == 0) {
        throw new NoDataException();
    }//from  w ww  . j  av  a2  s  . c  o m
    if (zval.length != fval.length) {
        throw new DimensionMismatchException(zval.length, fval.length);
    }

    MathArrays.checkOrder(yval);
    MathArrays.checkOrder(xval);

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

    final double[] zzval;
    float[][][] ffval = new float[zLen][][];

    if (zval[0] > zval[zLen - 1]) {
        zzval = VectorMath.selectionSort(zval);
        for (int i = 0; i < zLen; i++) {
            ffval[i] = fval[zLen - 1 - i];
        }
    } else {
        zzval = zval;
        ffval = fval;
    }

    // Samples, re-ordered as (z, x, y) and (y, z, x) tuplets
    // fvalXY[k][i][j] = f(xval[i], yval[j], zval[k])
    // fvalZX[j][k][i] = f(xval[i], yval[j], zval[k])
    final double[][][] fvalXY = new double[xLen][zLen][yLen];
    final double[][][] fvalZX = new double[yLen][xLen][zLen];
    for (int i = 0; i < zLen; i++) {
        if (ffval[i].length != yLen) {
            throw new DimensionMismatchException(ffval[i].length, yLen);
        }

        for (int j = 0; j < yLen; j++) {
            if (ffval[i][j].length != xLen) {
                throw new DimensionMismatchException(ffval[i][j].length, xLen);
            }

            for (int k = 0; k < xLen; k++) {
                final double v = ffval[i][j][k];
                fvalXY[k][i][j] = v;
                fvalZX[j][k][i] = v;
            }
        }
    }

    final BicubicSplineInterpolator bsi = new BicubicSplineInterpolator();

    // For each line x[i] (0 <= i < xLen), construct a 2D spline in y and z
    final BicubicSplineInterpolatingFunction[] xSplineYZ = new BicubicSplineInterpolatingFunction[zLen];
    for (int i = 0; i < zLen; i++) {
        xSplineYZ[i] = bsi.interpolate(yval, xval, ffval[i]);
    }

    // For each line y[j] (0 <= j < yLen), construct a 2D spline in z and x
    final BicubicSplineInterpolatingFunction[] ySplineZX = new BicubicSplineInterpolatingFunction[yLen];
    for (int j = 0; j < yLen; j++) {
        ySplineZX[j] = bsi.interpolate(xval, zzval, fvalZX[j]);
    }

    // For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y
    final BicubicSplineInterpolatingFunction[] zSplineXY = new BicubicSplineInterpolatingFunction[xLen];
    for (int k = 0; k < xLen; k++) {
        zSplineXY[k] = bsi.interpolate(zzval, yval, fvalXY[k]);
    }

    // Partial derivatives wrt x and wrt y
    final double[][][] dFdX = new double[zLen][yLen][xLen];
    final double[][][] dFdY = new double[zLen][yLen][xLen];
    final double[][][] d2FdXdY = new double[zLen][yLen][xLen];
    for (int k = 0; k < xLen; k++) {
        final BicubicSplineInterpolatingFunction f = zSplineXY[k];
        for (int i = 0; i < zLen; i++) {
            final double x = zzval[i];
            for (int j = 0; j < yLen; j++) {
                final double y = yval[j];
                dFdX[i][j][k] = f.partialDerivativeX(x, y);
                dFdY[i][j][k] = f.partialDerivativeY(x, y);
                d2FdXdY[i][j][k] = f.partialDerivativeXY(x, y);
            }
        }
    }

    // Partial derivatives wrt y and wrt z
    final double[][][] dFdZ = new double[zLen][yLen][xLen];
    final double[][][] d2FdYdZ = new double[zLen][yLen][xLen];
    for (int i = 0; i < zLen; i++) {
        final BicubicSplineInterpolatingFunction f = xSplineYZ[i];
        for (int j = 0; j < yLen; j++) {
            final double y = yval[j];
            for (int k = 0; k < xLen; k++) {
                final double z = xval[k];
                dFdZ[i][j][k] = f.partialDerivativeY(y, z);
                d2FdYdZ[i][j][k] = f.partialDerivativeXY(y, z);
            }
        }
    }

    // Partial derivatives wrt x and wrt z
    final double[][][] d2FdZdX = new double[zLen][yLen][xLen];
    for (int j = 0; j < yLen; j++) {
        final BicubicSplineInterpolatingFunction f = ySplineZX[j];
        for (int k = 0; k < xLen; k++) {
            final double z = xval[k];
            for (int i = 0; i < zLen; i++) {
                final double x = zzval[i];
                d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x);
            }
        }
    }

    // Third partial cross-derivatives
    final double[][][] d3FdXdYdZ = new double[zLen][yLen][xLen];
    for (int i = 0; i < zLen; i++) {
        final int nI = nextIndex(i, zLen);
        final int pI = previousIndex(i);
        for (int j = 0; j < yLen; j++) {
            final int nJ = nextIndex(j, yLen);
            final int pJ = previousIndex(j);
            for (int k = 0; k < xLen; k++) {
                final int nK = nextIndex(k, xLen);
                final int pK = previousIndex(k);

                // XXX Not sure about this formula
                d3FdXdYdZ[i][j][k] = (ffval[nI][nJ][nK] - ffval[nI][pJ][nK] - ffval[pI][nJ][nK]
                        + ffval[pI][pJ][nK] - ffval[nI][nJ][pK] + ffval[nI][pJ][pK] + ffval[pI][nJ][pK]
                        - ffval[pI][pJ][pK])
                        / ((zzval[nI] - zzval[pI]) * (yval[nJ] - yval[pJ]) * (xval[nK] - xval[pK]));
            }
        }
    }

    // Create the interpolating splines
    return new TricubicSplineInterpolatingFunction(zzval, yval, xval, ffval, dFdX, dFdY, dFdZ, d2FdXdY, d2FdZdX,
            d2FdYdZ, d3FdXdYdZ);
}