1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.analysis.interpolation; 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.DimensionMismatchException; 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.FunctionEvaluationException; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.analysis.BivariateRealFunction; 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.NoDataException; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.OutOfRangeException; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.MathUtils; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Function that implements the 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <a href="http://en.wikipedia.org/wiki/Bicubic_interpolation"> 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * bicubic spline interpolation</a>. 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision$ $Date$ 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.1 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class BicubicSplineInterpolatingFunction 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond implements BivariateRealFunction { 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Matrix to compute the spline coefficients from the function values 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * and function derivatives values 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final double[][] AINV = { 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 }, 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { 0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 }, 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { -3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0 }, 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { 2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0 }, 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 }, 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { 0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0 }, 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { 0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0 }, 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { 0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0 }, 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { -3,0,3,0,0,0,0,0,-2,0,-1,0,0,0,0,0 }, 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { 0,0,0,0,-3,0,3,0,0,0,0,0,-2,0,-1,0 }, 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { 9,-9,-9,9,6,3,-6,-3,6,-6,3,-3,4,2,2,1 }, 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { -6,6,6,-6,-3,-3,3,3,-4,4,-2,2,-2,-2,-1,-1 }, 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { 2,0,-2,0,0,0,0,0,1,0,1,0,0,0,0,0 }, 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { 0,0,0,0,2,0,-2,0,0,0,0,0,1,0,1,0 }, 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { -6,6,6,-6,-4,-2,4,2,-3,3,-3,3,-2,-1,-2,-1 }, 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond { 4,-4,-4,4,2,2,-2,-2,2,-2,2,-2,1,1,1,1 } 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond }; 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Samples x-coordinates */ 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[] xval; 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Samples y-coordinates */ 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[] yval; 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Set of cubic splines patching the whole data grid */ 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final BicubicSplineFunction[][] splines; 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Partial derivatives 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The value of the first index determines the kind of derivatives: 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 0 = first partial derivatives wrt x 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 1 = first partial derivatives wrt y 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 2 = second partial derivatives wrt x 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 3 = second partial derivatives wrt y 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 4 = cross partial derivatives 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private BivariateRealFunction[][][] partialDerivatives = null; 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param x Sample values of the x-coordinate, in increasing order. 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y Sample values of the y-coordinate, in increasing order. 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param f Values of the function on every grid point. 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param dFdX Values of the partial derivative of function with respect 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * to x on every grid point. 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param dFdY Values of the partial derivative of function with respect 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * to y on every grid point. 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param d2FdXdY Values of the cross partial derivative of function on 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * every grid point. 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws DimensionMismatchException if the various arrays do not contain 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the expected number of elements. 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws org.apache.commons.math.exception.NonMonotonousSequenceException 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * if {@code x} or {@code y} are not strictly increasing. 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws NoDataException if any of the arrays has zero length. 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public BicubicSplineInterpolatingFunction(double[] x, 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[] y, 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[][] f, 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[][] dFdX, 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[][] dFdY, 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double[][] d2FdXdY) 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws DimensionMismatchException { 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int xLen = x.length; 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int yLen = y.length; 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) { 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new NoDataException(); 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (xLen != f.length) { 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new DimensionMismatchException(xLen, f.length); 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (xLen != dFdX.length) { 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new DimensionMismatchException(xLen, dFdX.length); 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (xLen != dFdY.length) { 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new DimensionMismatchException(xLen, dFdY.length); 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (xLen != d2FdXdY.length) { 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new DimensionMismatchException(xLen, d2FdXdY.length); 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond MathUtils.checkOrder(x); 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond MathUtils.checkOrder(y); 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond xval = x.clone(); 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yval = y.clone(); 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int lastI = xLen - 1; 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int lastJ = yLen - 1; 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond splines = new BicubicSplineFunction[lastI][lastJ]; 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < lastI; i++) { 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (f[i].length != yLen) { 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new DimensionMismatchException(f[i].length, yLen); 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (dFdX[i].length != yLen) { 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new DimensionMismatchException(dFdX[i].length, yLen); 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (dFdY[i].length != yLen) { 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new DimensionMismatchException(dFdY[i].length, yLen); 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (d2FdXdY[i].length != yLen) { 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new DimensionMismatchException(d2FdXdY[i].length, yLen); 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int ip1 = i + 1; 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < lastJ; j++) { 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int jp1 = j + 1; 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] beta = new double[] { 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1], 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1], 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1], 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1] 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond }; 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta)); 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@inheritDoc} 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double value(double x, double y) { 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int i = searchIndex(x, xval); 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (i == -1) { 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]); 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int j = searchIndex(y, yval); 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (j == -1) { 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]); 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return splines[i][j].value(xN, yN); 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param x x-coordinate. 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y y-coordinate. 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the value at point (x, y) of the first partial derivative with 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * respect to x. 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double partialDerivativeX(double x, double y) { 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return partialDerivative(0, x, y); 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param x x-coordinate. 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y y-coordinate. 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the value at point (x, y) of the first partial derivative with 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * respect to y. 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double partialDerivativeY(double x, double y) { 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return partialDerivative(1, x, y); 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param x x-coordinate. 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y y-coordinate. 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the value at point (x, y) of the second partial derivative with 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * respect to x. 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double partialDerivativeXX(double x, double y) { 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return partialDerivative(2, x, y); 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param x x-coordinate. 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y y-coordinate. 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the value at point (x, y) of the second partial derivative with 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * respect to y. 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double partialDerivativeYY(double x, double y) { 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return partialDerivative(3, x, y); 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param x x-coordinate. 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y y-coordinate. 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the value at point (x, y) of the second partial cross-derivative. 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.2 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double partialDerivativeXY(double x, double y) { 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return partialDerivative(4, x, y); 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param which First index in {@link #partialDerivatives}. 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param x x-coordinate. 228dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param y y-coordinate. 229dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the value at point (x, y) of the selected partial derivative. 230dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws FunctionEvaluationException 231dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 232dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double partialDerivative(int which, double x, double y) { 233dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (partialDerivatives == null) { 234dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond computePartialDerivatives(); 235dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 236dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 237dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int i = searchIndex(x, xval); 238dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (i == -1) { 239dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new OutOfRangeException(x, xval[0], xval[xval.length - 1]); 240dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 241dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int j = searchIndex(y, yval); 242dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (j == -1) { 243dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new OutOfRangeException(y, yval[0], yval[yval.length - 1]); 244dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 245dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 246dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]); 247dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]); 248dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 249dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond try { 250dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return partialDerivatives[which][i][j].value(xN, yN); 251dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (FunctionEvaluationException fee) { 252dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // this should never happen 253dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new RuntimeException(fee); 254dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 255dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 256dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 257dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 258dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 259dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Compute all partial derivatives. 260dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 261dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private void computePartialDerivatives() { 262dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int lastI = xval.length - 1; 263dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int lastJ = yval.length - 1; 264dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond partialDerivatives = new BivariateRealFunction[5][lastI][lastJ]; 265dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 266dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < lastI; i++) { 267dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < lastJ; j++) { 268dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final BicubicSplineFunction f = splines[i][j]; 269dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond partialDerivatives[0][i][j] = f.partialDerivativeX(); 270dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond partialDerivatives[1][i][j] = f.partialDerivativeY(); 271dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond partialDerivatives[2][i][j] = f.partialDerivativeXX(); 272dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond partialDerivatives[3][i][j] = f.partialDerivativeYY(); 273dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond partialDerivatives[4][i][j] = f.partialDerivativeXY(); 274dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 275dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 276dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 277dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 278dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 279dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param c Coordinate. 280dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param val Coordinate samples. 281dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the index in {@code val} corresponding to the interval 282dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * containing {@code c}, or {@code -1} if {@code c} is out of the 283dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * range defined by the end values of {@code val}. 284dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 285dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private int searchIndex(double c, double[] val) { 286dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (c < val[0]) { 287dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return -1; 288dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 289dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 290dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int max = val.length; 291dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 1; i < max; i++) { 292dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (c <= val[i]) { 293dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return i - 1; 294dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 295dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 296dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 297dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return -1; 298dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 299dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 300dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 301dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Compute the spline coefficients from the list of function values and 302dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * function partial derivatives values at the four corners of a grid 303dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * element. They must be specified in the following order: 304dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <ul> 305dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f(0,0)</li> 306dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f(1,0)</li> 307dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f(0,1)</li> 308dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f(1,1)</li> 309dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f<sub>x</sub>(0,0)</li> 310dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f<sub>x</sub>(1,0)</li> 311dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f<sub>x</sub>(0,1)</li> 312dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f<sub>x</sub>(1,1)</li> 313dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f<sub>y</sub>(0,0)</li> 314dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f<sub>y</sub>(1,0)</li> 315dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f<sub>y</sub>(0,1)</li> 316dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f<sub>y</sub>(1,1)</li> 317dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f<sub>xy</sub>(0,0)</li> 318dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f<sub>xy</sub>(1,0)</li> 319dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f<sub>xy</sub>(0,1)</li> 320dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li>f<sub>xy</sub>(1,1)</li> 321dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </ul> 322dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * where the subscripts indicate the partial derivative with respect to 323dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the corresponding variable(s). 324dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 325dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param beta List of function values and function partial derivatives 326dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * values. 327dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the spline coefficients. 328dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 329dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double[] computeSplineCoefficients(double[] beta) { 330dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] a = new double[16]; 331dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 332dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < 16; i++) { 333dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double result = 0; 334dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] row = AINV[i]; 335dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < 16; j++) { 336dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond result += row[j] * beta[j]; 337dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 338dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond a[i] = result; 339dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 340dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 341dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return a; 342dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 343dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 344dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 345dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 346dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 2D-spline function. 347dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 348dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision$ $Date$ 349dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 350dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondclass BicubicSplineFunction 351dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond implements BivariateRealFunction { 352dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 353dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Number of points. */ 354dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final short N = 4; 355dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 356dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Coefficients */ 357dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[][] a; 358dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 359dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** First partial derivative along x. */ 360dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private BivariateRealFunction partialDerivativeX; 361dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 362dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** First partial derivative along y. */ 363dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private BivariateRealFunction partialDerivativeY; 364dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 365dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Second partial derivative along x. */ 366dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private BivariateRealFunction partialDerivativeXX; 367dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 368dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Second partial derivative along y. */ 369dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private BivariateRealFunction partialDerivativeYY; 370dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 371dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Second crossed partial derivative. */ 372dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private BivariateRealFunction partialDerivativeXY; 373dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 374dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 375dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Simple constructor. 376dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param a Spline coefficients 377dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 378dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public BicubicSplineFunction(double[] a) { 379dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.a = new double[N][N]; 380dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < N; i++) { 381dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < N; j++) { 382dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.a[i][j] = a[i + N * j]; 383dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 384dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 385dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 386dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 387dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 388dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@inheritDoc} 389dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 390dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double value(double x, double y) { 391dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (x < 0 || x > 1) { 392dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new OutOfRangeException(x, 0, 1); 393dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 394dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (y < 0 || y > 1) { 395dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new OutOfRangeException(y, 0, 1); 396dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 397dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 398dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double x2 = x * x; 399dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double x3 = x2 * x; 400dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] pX = {1, x, x2, x3}; 401dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 402dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double y2 = y * y; 403dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double y3 = y2 * y; 404dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] pY = {1, y, y2, y3}; 405dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 406dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return apply(pX, pY, a); 407dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 408dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 409dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 410dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Compute the value of the bicubic polynomial. 411dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 412dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param pX Powers of the x-coordinate. 413dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param pY Powers of the y-coordinate. 414dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param coeff Spline coefficients. 415dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the interpolated value. 416dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 417dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double apply(double[] pX, double[] pY, double[][] coeff) { 418dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double result = 0; 419dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < N; i++) { 420dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < N; j++) { 421dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond result += coeff[i][j] * pX[i] * pY[j]; 422dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 423dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 424dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 425dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return result; 426dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 427dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 428dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 429dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the partial derivative wrt {@code x}. 430dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 431dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public BivariateRealFunction partialDerivativeX() { 432dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (partialDerivativeX == null) { 433dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond computePartialDerivatives(); 434dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 435dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 436dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return partialDerivativeX; 437dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 438dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 439dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the partial derivative wrt {@code y}. 440dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 441dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public BivariateRealFunction partialDerivativeY() { 442dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (partialDerivativeY == null) { 443dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond computePartialDerivatives(); 444dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 445dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 446dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return partialDerivativeY; 447dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 448dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 449dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the second partial derivative wrt {@code x}. 450dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 451dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public BivariateRealFunction partialDerivativeXX() { 452dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (partialDerivativeXX == null) { 453dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond computePartialDerivatives(); 454dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 455dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 456dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return partialDerivativeXX; 457dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 458dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 459dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the second partial derivative wrt {@code y}. 460dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 461dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public BivariateRealFunction partialDerivativeYY() { 462dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (partialDerivativeYY == null) { 463dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond computePartialDerivatives(); 464dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 465dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 466dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return partialDerivativeYY; 467dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 468dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 469dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the second partial cross-derivative. 470dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 471dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public BivariateRealFunction partialDerivativeXY() { 472dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (partialDerivativeXY == null) { 473dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond computePartialDerivatives(); 474dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 475dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 476dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return partialDerivativeXY; 477dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 478dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 479dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 480dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Compute all partial derivatives functions. 481dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 482dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private void computePartialDerivatives() { 483dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[][] aX = new double[N][N]; 484dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[][] aY = new double[N][N]; 485dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[][] aXX = new double[N][N]; 486dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[][] aYY = new double[N][N]; 487dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[][] aXY = new double[N][N]; 488dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 489dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < N; i++) { 490dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < N; j++) { 491dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double c = a[i][j]; 492dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond aX[i][j] = i * c; 493dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond aY[i][j] = j * c; 494dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond aXX[i][j] = (i - 1) * aX[i][j]; 495dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond aYY[i][j] = (j - 1) * aY[i][j]; 496dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond aXY[i][j] = j * aX[i][j]; 497dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 498dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 499dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 500dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond partialDerivativeX = new BivariateRealFunction() { 501dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double value(double x, double y) { 502dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double x2 = x * x; 503dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] pX = {0, 1, x, x2}; 504dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 505dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double y2 = y * y; 506dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double y3 = y2 * y; 507dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] pY = {1, y, y2, y3}; 508dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 509dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return apply(pX, pY, aX); 510dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 511dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond }; 512dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond partialDerivativeY = new BivariateRealFunction() { 513dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double value(double x, double y) { 514dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double x2 = x * x; 515dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double x3 = x2 * x; 516dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] pX = {1, x, x2, x3}; 517dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 518dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double y2 = y * y; 519dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] pY = {0, 1, y, y2}; 520dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 521dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return apply(pX, pY, aY); 522dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 523dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond }; 524dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond partialDerivativeXX = new BivariateRealFunction() { 525dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double value(double x, double y) { 526dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] pX = {0, 0, 1, x}; 527dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 528dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double y2 = y * y; 529dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double y3 = y2 * y; 530dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] pY = {1, y, y2, y3}; 531dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 532dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return apply(pX, pY, aXX); 533dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 534dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond }; 535dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond partialDerivativeYY = new BivariateRealFunction() { 536dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double value(double x, double y) { 537dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double x2 = x * x; 538dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double x3 = x2 * x; 539dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] pX = {1, x, x2, x3}; 540dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 541dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] pY = {0, 0, 1, y}; 542dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 543dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return apply(pX, pY, aYY); 544dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 545dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond }; 546dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond partialDerivativeXY = new BivariateRealFunction() { 547dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double value(double x, double y) { 548dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double x2 = x * x; 549dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] pX = {0, 1, x, x2}; 550dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 551dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double y2 = y * y; 552dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] pY = {0, 1, y, y2}; 553dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 554dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return apply(pX, pY, aXY); 555dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 556dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond }; 557dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 558dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 559