commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject svn commit: r937080 [1/2] - in /commons/proper/math/trunk/src: main/java/org/apache/commons/math/analysis/ main/java/org/apache/commons/math/analysis/interpolation/ site/xdoc/ site/xdoc/userguide/ test/java/org/apache/commons/math/analysis/interpolation/
Date Thu, 22 Apr 2010 22:09:21 GMT
Author: erans
Date: Thu Apr 22 22:09:21 2010
New Revision: 937080

URL: http://svn.apache.org/viewvc?rev=937080&view=rev
Log:
MATH-366

Added:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/TrivariateRealFunction.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.java   (with props)
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TrivariateRealGridInterpolator.java   (with props)
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunctionTest.java   (with props)
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatorTest.java   (with props)
Modified:
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunction.java
    commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java
    commons/proper/math/trunk/src/site/xdoc/changes.xml
    commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml
    commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/TrivariateRealFunction.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/TrivariateRealFunction.java?rev=937080&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/TrivariateRealFunction.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/TrivariateRealFunction.java Thu Apr 22 22:09:21 2010
@@ -0,0 +1,40 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package org.apache.commons.math.analysis;
+
+import org.apache.commons.math.FunctionEvaluationException;
+
+/**
+ * An interface representing a trivariate real function.
+ *
+ * @since 2.2
+ * @version $Revision$ $Date$
+ */
+public interface TrivariateRealFunction {
+    /**
+     * Compute the value for the function.
+     *
+     * @param x x-coordinate for which the function value should be computed.
+     * @param y y-coordinate for which the function value should be computed.
+     * @param z z-coordinate for which the function value should be computed.
+     * @return the value.
+     * @throws FunctionEvaluationException if the function evaluation fails.
+     */
+    public double value(double x, double y, double z)
+        throws FunctionEvaluationException;
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/TrivariateRealFunction.java
------------------------------------------------------------------------------
    svn:eol-style = native

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunction.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunction.java?rev=937080&r1=937079&r2=937080&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunction.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunction.java Thu Apr 22 22:09:21 2010
@@ -18,6 +18,7 @@ package org.apache.commons.math.analysis
 
 import org.apache.commons.math.util.MathUtils;
 import org.apache.commons.math.MathRuntimeException;
+import org.apache.commons.math.FunctionEvaluationException;
 import org.apache.commons.math.DimensionMismatchException;
 import org.apache.commons.math.analysis.BivariateRealFunction;
 
@@ -58,19 +59,29 @@ public class BicubicSplineInterpolatingF
     private final double[] xval;
     /** Samples y-coordinates */
     private final double[] yval;
-    /** Set of cubic splines pacthing the whole data grid */
+    /** Set of cubic splines patching the whole data grid */
     private final BicubicSplineFunction[][] splines;
+    /**
+     * Partial derivatives
+     * The value of the first index determines the kind of derivatives:
+     * 0 = first partial derivatives wrt x
+     * 1 = first partial derivatives wrt y
+     * 2 = second partial derivatives wrt x
+     * 3 = second partial derivatives wrt y
+     * 4 = cross partial derivatives
+     */
+    private BivariateRealFunction[][][] partialDerivatives = null;
 
     /**
-     * @param x Sample values of the x-coordinate, in increasing order
-     * @param y Sample values of the y-coordinate, in increasing order
-     * @param z Values of the function on every grid point
-     * @param dZdX Values of the partial derivative of function with respect
-     * to x on every grid point
-     * @param dZdY Values of the partial derivative of function with respect
-     * to y on every grid point
-     * @param dZdXdY Values of the cross partial derivative of function on
-     * every grid point
+     * @param x Sample values of the x-coordinate, in increasing order.
+     * @param y Sample values of the y-coordinate, in increasing order.
+     * @param f Values of the function on every grid point.
+     * @param dFdX Values of the partial derivative of function with respect
+     * to x on every grid point.
+     * @param dFdY Values of the partial derivative of function with respect
+     * to y on every grid point.
+     * @param d2FdXdY Values of the cross partial derivative of function on
+     * every grid point.
      * @throws DimensionMismatchException if the various arrays do not contain
      * the expected number of elements.
      * @throws IllegalArgumentException if {@code x} or {@code y} are not strictly
@@ -78,28 +89,28 @@ public class BicubicSplineInterpolatingF
      */
     public BicubicSplineInterpolatingFunction(double[] x,
                                               double[] y,
-                                              double[][] z,
-                                              double[][] dZdX,
-                                              double[][] dZdY,
-                                              double[][] dZdXdY)
+                                              double[][] f,
+                                              double[][] dFdX,
+                                              double[][] dFdY,
+                                              double[][] d2FdXdY)
         throws DimensionMismatchException {
         final int xLen = x.length;
         final int yLen = y.length;
 
-        if (xLen == 0 || yLen == 0 || z.length == 0 || z[0].length == 0) {
+        if (xLen == 0 || yLen == 0 || f.length == 0 || f[0].length == 0) {
             throw MathRuntimeException.createIllegalArgumentException("no data");
         }
-        if (xLen != z.length) {
-            throw new DimensionMismatchException(xLen, z.length);
+        if (xLen != f.length) {
+            throw new DimensionMismatchException(xLen, f.length);
         }
-        if (xLen != dZdX.length) {
-            throw new DimensionMismatchException(xLen, dZdX.length);
+        if (xLen != dFdX.length) {
+            throw new DimensionMismatchException(xLen, dFdX.length);
         }
-        if (xLen != dZdY.length) {
-            throw new DimensionMismatchException(xLen, dZdY.length);
+        if (xLen != dFdY.length) {
+            throw new DimensionMismatchException(xLen, dFdY.length);
         }
-        if (xLen != dZdXdY.length) {
-            throw new DimensionMismatchException(xLen, dZdXdY.length);
+        if (xLen != d2FdXdY.length) {
+            throw new DimensionMismatchException(xLen, d2FdXdY.length);
         }
 
         MathUtils.checkOrder(x, 1, true);
@@ -113,26 +124,26 @@ public class BicubicSplineInterpolatingF
         splines = new BicubicSplineFunction[lastI][lastJ];
 
         for (int i = 0; i < lastI; i++) {
-            if (z[i].length != yLen) {
-                throw new DimensionMismatchException(z[i].length, yLen);
+            if (f[i].length != yLen) {
+                throw new DimensionMismatchException(f[i].length, yLen);
             }
-            if (dZdX[i].length != yLen) {
-                throw new DimensionMismatchException(dZdX[i].length, yLen);
+            if (dFdX[i].length != yLen) {
+                throw new DimensionMismatchException(dFdX[i].length, yLen);
             }
-            if (dZdY[i].length != yLen) {
-                throw new DimensionMismatchException(dZdY[i].length, yLen);
+            if (dFdY[i].length != yLen) {
+                throw new DimensionMismatchException(dFdY[i].length, yLen);
             }
-            if (dZdXdY[i].length != yLen) {
-                throw new DimensionMismatchException(dZdXdY[i].length, yLen);
+            if (d2FdXdY[i].length != yLen) {
+                throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
             }
             final int ip1 = i + 1;
             for (int j = 0; j < lastJ; j++) {
                 final int jp1 = j + 1;
                 final double[] beta = new double[] {
-                    z[i][j],      z[ip1][j],      z[i][jp1],      z[ip1][jp1],
-                    dZdX[i][j],   dZdX[ip1][j],   dZdX[i][jp1],   dZdX[ip1][jp1],
-                    dZdY[i][j],   dZdY[ip1][j],   dZdY[i][jp1],   dZdY[ip1][jp1],
-                    dZdXdY[i][j], dZdXdY[ip1][j], dZdXdY[i][jp1], dZdXdY[ip1][jp1]
+                    f[i][j], f[ip1][j], f[i][jp1], f[ip1][jp1],
+                    dFdX[i][j], dFdX[ip1][j], dFdX[i][jp1], dFdX[ip1][jp1],
+                    dFdY[i][j], dFdY[ip1][j], dFdY[i][jp1], dFdY[ip1][jp1],
+                    d2FdXdY[i][j], d2FdXdY[ip1][j], d2FdXdY[i][jp1], d2FdXdY[ip1][jp1]
                 };
 
                 splines[i][j] = new BicubicSplineFunction(computeSplineCoefficients(beta));
@@ -162,18 +173,119 @@ public class BicubicSplineInterpolatingF
     }
 
     /**
-     * @param c coordinate
-     * @param val coordinate samples
+     * @param x x-coordinate.
+     * @param y y-coordinate.
+     * @return the value at point (x, y) of the first partial derivative with
+     * respect to x.
+     */
+    public double partialDerivativeX(double x, double y) {
+        return partialDerivative(0, x, y);
+    }
+    /**
+     * @param x x-coordinate.
+     * @param y y-coordinate.
+     * @return the value at point (x, y) of the first partial derivative with
+     * respect to y.
+     */
+    public double partialDerivativeY(double x, double y) {
+        return partialDerivative(1, x, y);
+    }
+    /**
+     * @param x x-coordinate.
+     * @param y y-coordinate.
+     * @return the value at point (x, y) of the second partial derivative with
+     * respect to x.
+     */
+    public double partialDerivativeXX(double x, double y) {
+        return partialDerivative(2, x, y);
+    }
+    /**
+     * @param x x-coordinate.
+     * @param y y-coordinate.
+     * @return the value at point (x, y) of the second partial derivative with
+     * respect to y.
+     */
+    public double partialDerivativeYY(double x, double y) {
+        return partialDerivative(3, x, y);
+    }
+    /**
+     * @param x x-coordinate.
+     * @param y y-coordinate.
+     * @return the value at point (x, y) of the second partial cross-derivative.
+     */
+    public double partialDerivativeXY(double x, double y) {
+        return partialDerivative(4, x, y);
+    }
+
+    /**
+     * @param which First index in {@link #partialDerivatives}.
+     * @param x x-coordinate.
+     * @param y y-coordinate.
+     * @return the value at point (x, y) of the selected partial derivative.
+     */
+    private double partialDerivative(int which, double x, double y) {
+        if (partialDerivatives == null) {
+            computePartialDerivatives();
+        }
+
+        final int i = searchIndex(x, xval);
+        if (i == -1) {
+            throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
+                                                                      x, xval[0], xval[xval.length - 1]);
+        }
+        final int j = searchIndex(y, yval);
+        if (j == -1) {
+            throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
+                                                                      y, yval[0], yval[yval.length - 1]);
+        }
+
+        final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
+        final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
+
+        double result = Double.NaN;
+        try { // Workaround to avoid carrying around "try" blocks for an
+              // exception that will never happen
+            result = partialDerivatives[which][i][j].value(xN, yN);
+        } catch (FunctionEvaluationException e) {
+            // Will never happen
+        }
+
+        return result;
+    }
+
+    /**
+     * Compute all partial derivatives.
+     */
+    private void computePartialDerivatives() {
+        final int lastI = xval.length - 1;
+        final int lastJ = yval.length - 1;
+        partialDerivatives = new BivariateRealFunction[5][lastI][lastJ];
+
+        for (int i = 0; i < lastI; i++) {
+            for (int j = 0; j < lastJ; j++) {
+                final BicubicSplineFunction f = splines[i][j];
+                partialDerivatives[0][i][j] = f.partialDerivativeX();
+                partialDerivatives[1][i][j] = f.partialDerivativeY();
+                partialDerivatives[2][i][j] = f.partialDerivativeXX();
+                partialDerivatives[3][i][j] = f.partialDerivativeYY();
+                partialDerivatives[4][i][j] = f.partialDerivativeXY();
+            }
+        }
+    }
+
+    /**
+     * @param c Coordinate.
+     * @param val Coordinate samples.
      * @return the index in {@code val} corresponding to the interval
      * containing {@code c}, or {@code -1} if {@code c} is out of the
-     * range defined by the end values of {@code val}
+     * range defined by the end values of {@code val}.
      */
     private int searchIndex(double c, double[] val) {
         if (c < val[0]) {
             return -1;
         }
 
-        int max = val.length;
+        final int max = val.length;
         for (int i = 1; i < max; i++) {
             if (c <= val[i]) {
                 return i - 1;
@@ -192,22 +304,25 @@ public class BicubicSplineInterpolatingF
      *  <li>f(1,0)</li>
      *  <li>f(0,1)</li>
      *  <li>f(1,1)</li>
-     *  <li>fx(0,0)</li>
-     *  <li>fx(1,0)</li>
-     *  <li>fx(0,1)</li>
-     *  <li>fx(1,1)</li>
-     *  <li>fy(0,0)</li>
-     *  <li>fy(1,0)</li>
-     *  <li>fy(0,1)</li>
-     *  <li>fy(1,1)</li>
-     *  <li>fxy(0,0)</li>
-     *  <li>fxy(1,0)</li>
-     *  <li>fxy(0,1)</li>
-     *  <li>fxy(1,1)</li>
+     *  <li>f<sub>x</sub>(0,0)</li>
+     *  <li>f<sub>x</sub>(1,0)</li>
+     *  <li>f<sub>x</sub>(0,1)</li>
+     *  <li>f<sub>x</sub>(1,1)</li>
+     *  <li>f<sub>y</sub>(0,0)</li>
+     *  <li>f<sub>y</sub>(1,0)</li>
+     *  <li>f<sub>y</sub>(0,1)</li>
+     *  <li>f<sub>y</sub>(1,1)</li>
+     *  <li>f<sub>xy</sub>(0,0)</li>
+     *  <li>f<sub>xy</sub>(1,0)</li>
+     *  <li>f<sub>xy</sub>(0,1)</li>
+     *  <li>f<sub>xy</sub>(1,1)</li>
      * </ul>
+     * where the subscripts indicate the partial derivative with respect to
+     * the corresponding variable(s).
+     *
      * @param beta List of function values and function partial derivatives
-     * values
-     * @return the spline coefficients
+     * values.
+     * @return the spline coefficients.
      */
     private double[] computeSplineCoefficients(double[] beta) {
         final double[] a = new double[16];
@@ -224,6 +339,7 @@ public class BicubicSplineInterpolatingF
         return a;
     }
 }
+
 /**
  * 2D-spline function.
  *
@@ -231,41 +347,29 @@ public class BicubicSplineInterpolatingF
  */
 class BicubicSplineFunction
     implements BivariateRealFunction {
-//CHECKSTYLE: stop MultipleVariableDeclarations
+    private static final short N = 4;
     /** Coefficients */
-    private final double
-        a00, a01, a02, a03,
-        a10, a11, a12, a13,
-        a20, a21, a22, a23,
-        a30, a31, a32, a33;
-//CHECKSTYLE: resume MultipleVariableDeclarations
+    private final double[][] a = new double[N][N];
+    /** Partial derivatives */
+    BivariateRealFunction partialDerivativeX = null;
+    BivariateRealFunction partialDerivativeY = null;
+    BivariateRealFunction partialDerivativeXX = null;
+    BivariateRealFunction partialDerivativeYY = null;
+    BivariateRealFunction partialDerivativeXY = null;
 
     /**
      * @param a Spline coefficients
      */
-    public BicubicSplineFunction(double[] a) {
-        this.a00 = a[0];
-        this.a10 = a[1];
-        this.a20 = a[2];
-        this.a30 = a[3];
-        this.a01 = a[4];
-        this.a11 = a[5];
-        this.a21 = a[6];
-        this.a31 = a[7];
-        this.a02 = a[8];
-        this.a12 = a[9];
-        this.a22 = a[10];
-        this.a32 = a[11];
-        this.a03 = a[12];
-        this.a13 = a[13];
-        this.a23 = a[14];
-        this.a33 = a[15];
+    public BicubicSplineFunction(double[] aV) {
+        for (int i = 0; i < N; i++) {
+            for (int j = 0; j < N; j++) {
+                a[i][j] = aV[i + N * j];
+            }
+        }
     }
 
     /**
-     * @param x x-coordinate of the interpolation point
-     * @param y y-coordinate of the interpolation point
-     * @return the interpolated value.
+     * {@inheritDoc}
      */
     public double value(double x, double y) {
         if (x < 0 || x > 1) {
@@ -279,12 +383,162 @@ class BicubicSplineFunction
 
         final double x2 = x * x;
         final double x3 = x2 * x;
+        final double[] pX = {1, x, x2, x3};
+
         final double y2 = y * y;
         final double y3 = y2 * y;
+        final double[] pY = {1, y, y2, y3};
+
+        return apply(pX, pY, a);
+    }
+
+    /**
+     * Compute the value of the bicubic polynomial.
+     *
+     * @param pX Powers of the x-coordinate.
+     * @param pY Powers of the y-coordinate.
+     * @param coeff Spline coefficients.
+     * @return the interpolated value.
+     */
+    private double apply(double[] pX, double[] pY, double[][] coeff) {
+        double result = 0;
+        for (int i = 0; i < N; i++) {
+            for (int j = 0; j < N; j++) {
+                result += coeff[i][j] * pX[i] * pY[j];
+            }
+        }
+
+        return result;
+    }
+
+    /**
+     * @return the partial derivative wrt {@code x}.
+     */
+    public BivariateRealFunction partialDerivativeX() {
+        if (partialDerivativeX == null) {
+            computePartialDerivatives();
+        }
+
+        return partialDerivativeX;
+    }
+    /**
+     * @return the partial derivative wrt {@code y}.
+     */
+    public BivariateRealFunction partialDerivativeY() {
+        if (partialDerivativeY == null) {
+            computePartialDerivatives();
+        }
+
+        return partialDerivativeY;
+    }
+    /**
+     * @return the second partial derivative wrt {@code x}.
+     */
+    public BivariateRealFunction partialDerivativeXX() {
+        if (partialDerivativeXX == null) {
+            computePartialDerivatives();
+        }
+
+        return partialDerivativeXX;
+    }
+    /**
+     * @return the second partial derivative wrt {@code y}.
+     */
+    public BivariateRealFunction partialDerivativeYY() {
+        if (partialDerivativeYY == null) {
+            computePartialDerivatives();
+        }
+
+        return partialDerivativeYY;
+    }
+    /**
+     * @return the second partial cross-derivative.
+     */
+    public BivariateRealFunction partialDerivativeXY() {
+        if (partialDerivativeXY == null) {
+            computePartialDerivatives();
+        }
+
+        return partialDerivativeXY;
+    }
+
+    /**
+     * Compute all partial derivatives functions.
+     */
+    private void computePartialDerivatives() {
+        final double[][] aX = new double[N][N];
+        final double[][] aY = new double[N][N];
+        final double[][] aXX = new double[N][N];
+        final double[][] aYY = new double[N][N];
+        final double[][] aXY = new double[N][N];
+
+        for (int i = 0; i < N; i++) {
+            for (int j = 0; j < N; j++) {
+                final double c = a[i][j];
+                aX[i][j] = i * c;
+                aY[i][j] = j * c;
+                aXX[i][j] = (i - 1) * aX[i][j];
+                aYY[i][j] = (j - 1) * aY[i][j];
+                aXY[i][j] = j * aX[i][j];
+            }
+        }
 
-        return a00 + a01 * y + a02 * y2 + a03 * y3 +
-            a10 * x + a11 * x * y + a12 * x * y2 + a13 * x * y3 +
-            a20 * x2 + a21 * x2 * y + a22 * x2 * y2 + a23 * x2 * y3 +
-            a30 * x3 + a31 * x3 * y + a32 * x3 * y2 + a33 * x3 * y3;
+        partialDerivativeX = new BivariateRealFunction() {
+                public double value(double x, double y)  {
+                    final double x2 = x * x;
+                    final double[] pX = {0, 1, x, x2};
+                    
+                    final double y2 = y * y;
+                    final double y3 = y2 * y;
+                    final double[] pY = {1, y, y2, y3};
+
+                    return apply(pX, pY, aX);
+                }
+            };
+        partialDerivativeY = new BivariateRealFunction() {
+                public double value(double x, double y)  {
+                    final double x2 = x * x;
+                    final double x3 = x2 * x;
+                    final double[] pX = {1, x, x2, x3};
+                    
+                    final double y2 = y * y;
+                    final double[] pY = {0, 1, y, y2};
+
+                    return apply(pX, pY, aY);
+                }
+            };
+        partialDerivativeXX = new BivariateRealFunction() {
+                public double value(double x, double y)  {
+                    final double[] pX = {0, 0, 1, x};
+                    
+                    final double y2 = y * y;
+                    final double y3 = y2 * y;
+                    final double[] pY = {1, y, y2, y3};
+
+                    return apply(pX, pY, aXX);
+                }
+            };
+        partialDerivativeYY = new BivariateRealFunction() {
+                public double value(double x, double y)  {
+                    final double x2 = x * x;
+                    final double x3 = x2 * x;
+                    final double[] pX = {1, x, x2, x3};
+                    
+                    final double[] pY = {0, 0, 1, y};
+
+                    return apply(pX, pY, aYY);
+                }
+            };
+        partialDerivativeXY = new BivariateRealFunction() {
+                public double value(double x, double y)  {
+                    final double x2 = x * x;
+                    final double[] pX = {0, 1, x, x2};
+                    
+                    final double y2 = y * y;
+                    final double[] pY = {0, 1, y, y2};
+
+                    return apply(pX, pY, aXY);
+                }
+            };
     }
 }

Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java?rev=937080&r1=937079&r2=937080&view=diff
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java (original)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/SmoothingPolynomialBicubicSplineInterpolator.java Thu Apr 22 22:09:21 2010
@@ -40,7 +40,7 @@ public class SmoothingPolynomialBicubicS
     private final PolynomialFitter yFitter;
 
     /**
-     * Default constructor. The degree of the fitting polynomials are set to 3.
+     * Default constructor. The degree of the fitting polynomials is set to 3.
      */
     public SmoothingPolynomialBicubicSplineInterpolator() {
         this(3);

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.java?rev=937080&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.java Thu Apr 22 22:09:21 2010
@@ -0,0 +1,487 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.analysis.interpolation;
+
+import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.MathRuntimeException;
+import org.apache.commons.math.DimensionMismatchException;
+import org.apache.commons.math.analysis.TrivariateRealFunction;
+
+/**
+ * Function that implements the
+ * <a href="http://en.wikipedia.org/wiki/Tricubic_interpolation">
+ * tricubic spline interpolation</a>, as proposed in
+ * <quote>
+ *  Tricubic interpolation in three dimensions<br/>
+ *  F. Lekien and J. Marsden<br/>
+ *  <em>Int. J. Numer. Meth. Engng</em> 2005; <b>63</b>:455-471
+ * </quote>
+ *
+ * @version $Revision$ $Date$
+ * @since 2.2
+ */
+public class TricubicSplineInterpolatingFunction
+    implements TrivariateRealFunction {
+    /**
+     * Matrix to compute the spline coefficients from the function values
+     * and function derivatives values
+     */
+    private static final double[][] AINV = {
+        { 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { -3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { -3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 9,-9,-9,9,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,0,0,0,0,0,0,0,0,4,2,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { -6,6,6,-6,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,0,0,0,0,0,0,0,0,-2,-2,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { -6,6,6,-6,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 4,-4,-4,4,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,0,0,0,0,-2,-1,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,0,0,0,0,1,1,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,-9,9,0,0,0,0,0,0,0,0,0,0,0,0,6,3,-6,-3,0,0,0,0,6,-6,3,-3,0,0,0,0,4,2,2,1,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-3,-3,3,3,0,0,0,0,-4,4,-2,2,0,0,0,0,-2,-2,-1,-1,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,6,-6,0,0,0,0,0,0,0,0,0,0,0,0,-4,-2,4,2,0,0,0,0,-3,3,-3,3,0,0,0,0,-2,-1,-2,-1,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,-4,4,0,0,0,0,0,0,0,0,0,0,0,0,2,2,-2,-2,0,0,0,0,2,-2,2,-2,0,0,0,0,1,1,1,1,0,0,0,0 },
+        {-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 9,-9,0,0,-9,9,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,0,0,0,0,0,0,0,0,4,2,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { -6,6,0,0,6,-6,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,0,0,0,0,0,0,0,0,-2,-2,0,0,-1,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,0,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,0,0,-1,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,9,-9,0,0,-9,9,0,0,0,0,0,0,0,0,0,0,6,3,0,0,-6,-3,0,0,0,0,0,0,0,0,0,0,6,-6,0,0,3,-3,0,0,4,2,0,0,2,1,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-3,-3,0,0,3,3,0,0,0,0,0,0,0,0,0,0,-4,4,0,0,-2,2,0,0,-2,-2,0,0,-1,-1,0,0 },
+        { 9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,9,0,-9,0,-9,0,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,0,3,0,-6,0,-3,0,6,0,-6,0,3,0,-3,0,0,0,0,0,0,0,0,0,4,0,2,0,2,0,1,0 },
+        { -27,27,27,-27,27,-27,-27,27,-18,-9,18,9,18,9,-18,-9,-18,18,-9,9,18,-18,9,-9,-18,18,18,-18,-9,9,9,-9,-12,-6,-6,-3,12,6,6,3,-12,-6,12,6,-6,-3,6,3,-12,12,-6,6,-6,6,-3,3,-8,-4,-4,-2,-4,-2,-2,-1 },
+        { 18,-18,-18,18,-18,18,18,-18,9,9,-9,-9,-9,-9,9,9,12,-12,6,-6,-12,12,-6,6,12,-12,-12,12,6,-6,-6,6,6,6,3,3,-6,-6,-3,-3,6,6,-6,-6,3,3,-3,-3,8,-8,4,-4,4,-4,2,-2,4,4,2,2,2,2,1,1 },
+        { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-3,0,-3,0,3,0,3,0,-4,0,4,0,-2,0,2,0,0,0,0,0,0,0,0,0,-2,0,-2,0,-1,0,-1,0 },
+        { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,9,-9,9,-9,-9,9,-9,9,12,-12,-12,12,6,-6,-6,6,6,3,6,3,-6,-3,-6,-3,8,4,-8,-4,4,2,-4,-2,6,-6,6,-6,3,-3,3,-3,4,2,4,2,2,1,2,1 },
+        { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-6,6,-6,6,6,-6,6,-6,-8,8,8,-8,-4,4,4,-4,-3,-3,-3,-3,3,3,3,3,-4,-4,4,4,-2,-2,2,2,-4,4,-4,4,-2,2,-2,2,-2,-2,-2,-2,-1,-1,-1,-1 },
+        { 2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { -6,6,0,0,6,-6,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,0,0,0,0,0,0,0,0,-2,-1,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 4,-4,0,0,-4,4,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-6,6,0,0,6,-6,0,0,0,0,0,0,0,0,0,0,-4,-2,0,0,4,2,0,0,0,0,0,0,0,0,0,0,-3,3,0,0,-3,3,0,0,-2,-1,0,0,-2,-1,0,0 },
+        { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,-4,0,0,-4,4,0,0,0,0,0,0,0,0,0,0,2,2,0,0,-2,-2,0,0,0,0,0,0,0,0,0,0,2,-2,0,0,2,-2,0,0,1,1,0,0,1,1,0,0 },
+        { -6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,-6,0,6,0,6,0,-6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-4,0,-2,0,4,0,2,0,-3,0,3,0,-3,0,3,0,0,0,0,0,0,0,0,0,-2,0,-1,0,-2,0,-1,0 },
+        { 18,-18,-18,18,-18,18,18,-18,12,6,-12,-6,-12,-6,12,6,12,-12,6,-6,-12,12,-6,6,9,-9,-9,9,9,-9,-9,9,8,4,4,2,-8,-4,-4,-2,6,3,-6,-3,6,3,-6,-3,6,-6,3,-3,6,-6,3,-3,4,2,2,1,4,2,2,1 },
+        { -12,12,12,-12,12,-12,-12,12,-6,-6,6,6,6,6,-6,-6,-8,8,-4,4,8,-8,4,-4,-6,6,6,-6,-6,6,6,-6,-4,-4,-2,-2,4,4,2,2,-3,-3,3,3,-3,-3,3,3,-4,4,-2,2,-4,4,-2,2,-2,-2,-1,-1,-2,-2,-1,-1 },
+        { 4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0 },
+        { 0,0,0,0,0,0,0,0,4,0,-4,0,-4,0,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,0,2,0,-2,0,-2,0,2,0,-2,0,2,0,-2,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0 },
+        { -12,12,12,-12,12,-12,-12,12,-8,-4,8,4,8,4,-8,-4,-6,6,-6,6,6,-6,6,-6,-6,6,6,-6,-6,6,6,-6,-4,-2,-4,-2,4,2,4,2,-4,-2,4,2,-4,-2,4,2,-3,3,-3,3,-3,3,-3,3,-2,-1,-2,-1,-2,-1,-2,-1 },
+        { 8,-8,-8,8,-8,8,8,-8,4,4,-4,-4,-4,-4,4,4,4,-4,4,-4,-4,4,-4,4,4,-4,-4,4,4,-4,-4,4,2,2,2,2,-2,-2,-2,-2,2,2,-2,-2,2,2,-2,-2,2,-2,2,-2,2,-2,2,-2,1,1,1,1,1,1,1,1 }
+    };
+
+    /** Samples x-coordinates */
+    private final double[] xval;
+    /** Samples y-coordinates */
+    private final double[] yval;
+    /** Samples z-coordinates */
+    private final double[] zval;
+    /** Set of cubic splines pacthing the whole data grid */
+    private final TricubicSplineFunction[][][] splines;
+
+    /**
+     * @param x Sample values of the x-coordinate, in increasing order.
+     * @param y Sample values of the y-coordinate, in increasing order.
+     * @param z Sample values of the y-coordinate, in increasing order.
+     * @param f Values of the function on every grid point.
+     * @param dFdX Values of the partial derivative of function with respect
+     * to x on every grid point.
+     * @param dFdY Values of the partial derivative of function with respect
+     * to y on every grid point.
+     * @param d2FdXdY Values of the cross partial derivative of function on
+     * every grid point.
+     * @param d2FdXdZ Values of the cross partial derivative of function on
+     * every grid point.
+     * @param d2FdYdZ Values of the cross partial derivative of function on
+     * every grid point.
+     * @param d3FdXdYdZ Values of the cross partial derivative of function on
+     * every grid point.
+     * @throws DimensionMismatchException if the various arrays do not contain
+     * the expected number of elements.
+     * @throws IllegalArgumentException if {@code x}, {@code y} or {@code z}
+     * are not strictly increasing.
+     */
+    public TricubicSplineInterpolatingFunction(double[] x,
+                                               double[] y,
+                                               double[] z,
+                                               double[][][] f,
+                                               double[][][] dFdX,
+                                               double[][][] dFdY,
+                                               double[][][] dFdZ,
+                                               double[][][] d2FdXdY,
+                                               double[][][] d2FdXdZ,
+                                               double[][][] d2FdYdZ,
+                                               double[][][] d3FdXdYdZ)
+        throws DimensionMismatchException {
+        final int xLen = x.length;
+        final int yLen = y.length;
+        final int zLen = z.length;
+
+        if (xLen == 0 || yLen == 0 || z.length == 0
+            || f.length == 0 || f[0].length == 0) {
+            throw MathRuntimeException.createIllegalArgumentException("no data");
+        }
+        if (xLen != f.length) {
+            throw new DimensionMismatchException(xLen, f.length);
+        }
+        if (xLen != dFdX.length) {
+            throw new DimensionMismatchException(xLen, dFdX.length);
+        }
+        if (xLen != dFdY.length) {
+            throw new DimensionMismatchException(xLen, dFdY.length);
+        }
+        if (xLen != dFdZ.length) {
+            throw new DimensionMismatchException(xLen, dFdZ.length);
+        }
+        if (xLen != d2FdXdY.length) {
+            throw new DimensionMismatchException(xLen, d2FdXdY.length);
+        }
+        if (xLen != d2FdXdZ.length) {
+            throw new DimensionMismatchException(xLen, d2FdXdZ.length);
+        }
+        if (xLen != d2FdYdZ.length) {
+            throw new DimensionMismatchException(xLen, d2FdYdZ.length);
+        }
+        if (xLen != d3FdXdYdZ.length) {
+            throw new DimensionMismatchException(xLen, d3FdXdYdZ.length);
+        }
+
+        MathUtils.checkOrder(x, 1, true);
+        MathUtils.checkOrder(y, 1, true);
+        MathUtils.checkOrder(z, 1, true);
+
+        xval = x.clone();
+        yval = y.clone();
+        zval = z.clone();
+
+        final int lastI = xLen - 1;
+        final int lastJ = yLen - 1;
+        final int lastK = zLen - 1;
+        splines = new TricubicSplineFunction[lastI][lastJ][lastK];
+
+        for (int i = 0; i < lastI; i++) {
+            if (f[i].length != yLen) {
+                throw new DimensionMismatchException(f[i].length, yLen);
+            }
+            if (dFdX[i].length != yLen) {
+                throw new DimensionMismatchException(dFdX[i].length, yLen);
+            }
+            if (dFdY[i].length != yLen) {
+                throw new DimensionMismatchException(dFdY[i].length, yLen);
+            }
+            if (dFdZ[i].length != yLen) {
+                throw new DimensionMismatchException(dFdZ[i].length, yLen);
+            }
+            if (d2FdXdY[i].length != yLen) {
+                throw new DimensionMismatchException(d2FdXdY[i].length, yLen);
+            }
+            if (d2FdXdZ[i].length != yLen) {
+                throw new DimensionMismatchException(d2FdXdZ[i].length, yLen);
+            }
+            if (d2FdYdZ[i].length != yLen) {
+                throw new DimensionMismatchException(d2FdYdZ[i].length, yLen);
+            }
+            if (d3FdXdYdZ[i].length != yLen) {
+                throw new DimensionMismatchException(d3FdXdYdZ[i].length, yLen);
+            }
+
+            final int ip1 = i + 1;
+            for (int j = 0; j < lastJ; j++) {
+                if (f[i][j].length != zLen) {
+                    throw new DimensionMismatchException(f[i][j].length, zLen);
+                }
+                if (dFdX[i][j].length != zLen) {
+                    throw new DimensionMismatchException(dFdX[i][j].length, zLen);
+                }
+                if (dFdY[i][j].length != zLen) {
+                    throw new DimensionMismatchException(dFdY[i][j].length, zLen);
+                }
+                if (dFdZ[i][j].length != zLen) {
+                    throw new DimensionMismatchException(dFdZ[i][j].length, zLen);
+                }
+                if (d2FdXdY[i][j].length != zLen) {
+                    throw new DimensionMismatchException(d2FdXdY[i][j].length, zLen);
+                }
+                if (d2FdXdZ[i][j].length != zLen) {
+                    throw new DimensionMismatchException(d2FdXdZ[i][j].length, zLen);
+                }
+                if (d2FdYdZ[i][j].length != zLen) {
+                    throw new DimensionMismatchException(d2FdYdZ[i][j].length, zLen);
+                }
+                if (d3FdXdYdZ[i][j].length != zLen) {
+                    throw new DimensionMismatchException(d3FdXdYdZ[i][j].length, zLen);
+                }
+
+                final int jp1 = j + 1;
+                for (int k = 0; k < lastK; k++) {
+                    final int kp1 = k + 1;
+
+                    final double[] beta = new double[] {
+                        f[i][j][k], f[ip1][j][k],
+                        f[i][jp1][k], f[ip1][jp1][k],
+                        f[i][j][kp1], f[ip1][j][kp1],
+                        f[i][jp1][kp1], f[ip1][jp1][kp1],
+
+                        dFdX[i][j][k], dFdX[ip1][j][k],
+                        dFdX[i][jp1][k], dFdX[ip1][jp1][k],
+                        dFdX[i][j][kp1], dFdX[ip1][j][kp1],
+                        dFdX[i][jp1][kp1], dFdX[ip1][jp1][kp1],
+
+                        dFdY[i][j][k], dFdY[ip1][j][k],
+                        dFdY[i][jp1][k], dFdY[ip1][jp1][k],
+                        dFdY[i][j][kp1], dFdY[ip1][j][kp1],
+                        dFdY[i][jp1][kp1], dFdY[ip1][jp1][kp1],
+
+                        dFdZ[i][j][k], dFdZ[ip1][j][k],
+                        dFdZ[i][jp1][k], dFdZ[ip1][jp1][k],
+                        dFdZ[i][j][kp1], dFdZ[ip1][j][kp1],
+                        dFdZ[i][jp1][kp1], dFdZ[ip1][jp1][kp1],
+                        
+                        d2FdXdY[i][j][k], d2FdXdY[ip1][j][k],
+                        d2FdXdY[i][jp1][k], d2FdXdY[ip1][jp1][k],
+                        d2FdXdY[i][j][kp1], d2FdXdY[ip1][j][kp1],
+                        d2FdXdY[i][jp1][kp1], d2FdXdY[ip1][jp1][kp1],
+
+                        d2FdXdZ[i][j][k], d2FdXdZ[ip1][j][k],
+                        d2FdXdZ[i][jp1][k], d2FdXdZ[ip1][jp1][k],
+                        d2FdXdZ[i][j][kp1], d2FdXdZ[ip1][j][kp1],
+                        d2FdXdZ[i][jp1][kp1], d2FdXdZ[ip1][jp1][kp1],
+
+                        d2FdYdZ[i][j][k], d2FdYdZ[ip1][j][k],
+                        d2FdYdZ[i][jp1][k], d2FdYdZ[ip1][jp1][k],
+                        d2FdYdZ[i][j][kp1], d2FdYdZ[ip1][j][kp1],
+                        d2FdYdZ[i][jp1][kp1], d2FdYdZ[ip1][jp1][kp1],
+
+                        d3FdXdYdZ[i][j][k], d3FdXdYdZ[ip1][j][k],
+                        d3FdXdYdZ[i][jp1][k], d3FdXdYdZ[ip1][jp1][k],
+                        d3FdXdYdZ[i][j][kp1], d3FdXdYdZ[ip1][j][kp1],
+                        d3FdXdYdZ[i][jp1][kp1], d3FdXdYdZ[ip1][jp1][kp1],
+                    };
+
+                    splines[i][j][k] = new TricubicSplineFunction(computeSplineCoefficients(beta));
+                }
+            }
+        }
+    }
+
+    /**
+     * {@inheritDoc}
+     */
+    public double value(double x, double y, double z) {
+        final int i = searchIndex(x, xval);
+        if (i == -1) {
+            throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
+                                                                      x, xval[0], xval[xval.length - 1]);
+        }
+        final int j = searchIndex(y, yval);
+        if (j == -1) {
+            throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
+                                                                      y, yval[0], yval[yval.length - 1]);
+        }
+        final int k = searchIndex(z, zval);
+        if (k == -1) {
+            throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
+                                                                      z, zval[0], zval[zval.length - 1]);
+        }
+
+        final double xN = (x - xval[i]) / (xval[i + 1] - xval[i]);
+        final double yN = (y - yval[j]) / (yval[j + 1] - yval[j]);
+        final double zN = (z - zval[k]) / (zval[k + 1] - zval[k]);
+
+        return splines[i][j][k].value(xN, yN, zN);
+    }
+
+    /**
+     * @param c Coordinate.
+     * @param val Coordinate samples.
+     * @return the index in {@code val} corresponding to the interval
+     * containing {@code c}, or {@code -1} if {@code c} is out of the
+     * range defined by the end values of {@code val}.
+     */
+    private int searchIndex(double c, double[] val) {
+        if (c < val[0]) {
+            return -1;
+        }
+
+        final int max = val.length;
+        for (int i = 1; i < max; i++) {
+            if (c <= val[i]) {
+                return i - 1;
+            }
+        }
+
+        return -1;
+    }
+
+    /**
+     * Compute the spline coefficients from the list of function values and
+     * function partial derivatives values at the four corners of a grid
+     * element. They must be specified in the following order:
+     * <ul>
+     *  <li>f(0,0,0)</li>
+     *  <li>f(1,0,0)</li>
+     *  <li>f(0,1,0)</li>
+     *  <li>f(1,1,0)</li>
+     *  <li>f(0,0,1)</li>
+     *  <li>f(1,0,1)</li>
+     *  <li>f(0,1,1)</li>
+     *  <li>f(1,1,1)</li>
+     *
+     *  <li>f<sub>x</sub>(0,0,0)</li>
+     *  <li>... <em>(same order as above)</em></li>
+     *  <li>f<sub>x</sub>(1,1,1)</li>
+     *
+     *  <li>f<sub>y</sub>(0,0,0)</li>
+     *  <li>... <em>(same order as above)</em></li>
+     *  <li>f<sub>y</sub>(1,1,1)</li>
+     *
+     *  <li>f<sub>z</sub>(0,0,0)</li>
+     *  <li>... <em>(same order as above)</em></li>
+     *  <li>f<sub>z</sub>(1,1,1)</li>
+     *
+     *  <li>f<sub>xy</sub>(0,0,0)</li>
+     *  <li>... <em>(same order as above)</em></li>
+     *  <li>f<sub>xy</sub>(1,1,1)</li>
+     *
+     *  <li>f<sub>xz</sub>(0,0,0)</li>
+     *  <li>... <em>(same order as above)</em></li>
+     *  <li>f<sub>xz</sub>(1,1,1)</li>
+     *
+     *  <li>f<sub>yz</sub>(0,0,0)</li>
+     *  <li>... <em>(same order as above)</em></li>
+     *  <li>f<sub>yz</sub>(1,1,1)</li>
+     *
+     *  <li>f<sub>xyz</sub>(0,0,0)</li>
+     *  <li>... <em>(same order as above)</em></li>
+     *  <li>f<sub>xyz</sub>(1,1,1)</li>
+     * </ul>
+     * where the subscripts indicate the partial derivative with respect to
+     * the corresponding variable(s).
+     *
+     * @param beta List of function values and function partial derivatives
+     * values.
+     * @return the spline coefficients.
+     */
+    private double[] computeSplineCoefficients(double[] beta) {
+        final int sz = 64;
+        final double[] a = new double[sz];
+
+        for (int i = 0; i < sz; i++) {
+            double result = 0;
+            final double[] row = AINV[i];
+            for (int j = 0; j < sz; j++) {
+                result += row[j] * beta[j];
+            }
+            a[i] = result;
+        }
+
+        return a;
+    }
+}
+
+/**
+ * 3D-spline function.
+ *
+ * @version $Revision$ $Date$
+ */
+class TricubicSplineFunction
+    implements TrivariateRealFunction {
+    private static final short N = 4;
+    private static final short N2 = N * N;
+    /** Coefficients */
+    private final double[][][] a = new double[N][N][N];
+
+    /**
+     * @param aV List of spline coefficients.
+     */
+    public TricubicSplineFunction(double[] aV) {
+        for (int i = 0; i < N; i++) {
+            for (int j = 0; j < N; j++) {
+                for (int k = 0; k < N; k++) {
+                    a[i][j][k] = aV[i + N * j + N2 * k];
+                }
+            }
+        }
+    }
+
+    /**
+     * @param x x-coordinate of the interpolation point.
+     * @param y y-coordinate of the interpolation point.
+     * @param z z-coordinate of the interpolation point.
+     * @return the interpolated value.
+     */
+    public double value(double x, double y, double z) {
+        if (x < 0 || x > 1) {
+            throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
+                                                                      x, 0, 1);
+        }
+        if (y < 0 || y > 1) {
+            throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
+                                                                      y, 0, 1);
+        }
+        if (z < 0 || z > 1) {
+            throw MathRuntimeException.createIllegalArgumentException("{0} out of [{1}, {2}] range",
+                                                                      z, 0, 1);
+        }
+
+        final double x2 = x * x;
+        final double x3 = x2 * x;
+        final double[] pX = { 1, x, x2, x3 };
+
+        final double y2 = y * y;
+        final double y3 = y2 * y;
+        final double[] pY = { 1, y, y2, y3 };
+
+        final double z2 = z * z;
+        final double z3 = z2 * z;
+        final double[] pZ = { 1, z, z2, z3 };
+
+        double result = 0;
+        for (int i = 0; i < N; i++) {
+            for (int j = 0; j < N; j++) {
+                for (int k = 0; k < N; k++) {
+                    result += a[i][j][k] * pX[i] * pY[j] * pZ[k];
+                }
+            }
+        }
+
+        return result;
+    }
+}
\ No newline at end of file

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.java?rev=937080&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.java Thu Apr 22 22:09:21 2010
@@ -0,0 +1,200 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.analysis.interpolation;
+
+import org.apache.commons.math.DimensionMismatchException;
+import org.apache.commons.math.MathRuntimeException;
+import org.apache.commons.math.MathException;
+import org.apache.commons.math.util.MathUtils;
+import org.apache.commons.math.analysis.UnivariateRealFunction;
+import org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction;
+
+/**
+ * Generates a tricubic interpolating function.
+ *
+ * @version $Revision$ $Date$
+ * @since 2.2
+ */
+public class TricubicSplineInterpolator
+    implements TrivariateRealGridInterpolator {
+    /**
+     * {@inheritDoc}
+     */
+    public TricubicSplineInterpolatingFunction interpolate(final double[] xval,
+                                                           final double[] yval,
+                                                           final double[] zval,
+                                                           final double[][][] fval)
+        throws MathException, IllegalArgumentException {
+        if (xval.length == 0 || yval.length == 0 || zval.length == 0 || fval.length == 0) {
+            throw MathRuntimeException.createIllegalArgumentException("no data");
+        }
+        if (xval.length != fval.length) {
+            throw new DimensionMismatchException(xval.length, fval.length);
+        }
+
+        MathUtils.checkOrder(xval, 1, true);
+        MathUtils.checkOrder(yval, 1, true);
+        MathUtils.checkOrder(zval, 1, true);
+
+        final int xLen = xval.length;
+        final int yLen = yval.length;
+        final int zLen = zval.length;
+
+        // 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[zLen][xLen][yLen];
+        final double[][][] fvalZX = new double[yLen][zLen][xLen];
+        for (int i = 0; i < xLen; i++) {
+            if (fval[i].length != yLen) {
+                throw new DimensionMismatchException(fval[i].length, yLen);
+            }
+
+            for (int j = 0; j < yLen; j++) {
+                if (fval[i][j].length != zLen) {
+                    throw new DimensionMismatchException(fval[i][j].length, zLen);
+                }
+                
+                for (int k = 0; k < zLen; k++) {
+                    final double v = fval[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[xLen];
+        for (int i = 0; i < xLen; i++) {
+            xSplineYZ[i] = bsi.interpolate(yval, zval, fval[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(zval, xval, fvalZX[j]);
+        }
+
+        // For each line z[k] (0 <= k < zLen), construct a 2D spline in x and y
+        final BicubicSplineInterpolatingFunction[] zSplineXY
+            = new BicubicSplineInterpolatingFunction[zLen];
+        for (int k = 0; k < zLen; k++) {
+            zSplineXY[k] = bsi.interpolate(xval, yval, fvalXY[k]);
+        }
+
+        // Partial derivatives wrt x and wrt y
+        final double[][][] dFdX = new double[xLen][yLen][zLen];
+        final double[][][] dFdY = new double[xLen][yLen][zLen];
+        final double[][][] d2FdXdY = new double[xLen][yLen][zLen];
+        for (int k = 0; k < zLen; k++) {
+            final BicubicSplineInterpolatingFunction f = zSplineXY[k];
+            for (int i = 0; i < xLen; i++) {
+                final double x = xval[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[xLen][yLen][zLen];
+        final double[][][] d2FdYdZ = new double[xLen][yLen][zLen];
+        for (int i = 0; i < xLen; i++) {
+            final BicubicSplineInterpolatingFunction f = xSplineYZ[i];
+            for (int j = 0; j < yLen; j++) {
+                final double y = yval[j];
+                for (int k = 0; k < zLen; k++) {
+                    final double z = zval[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[xLen][yLen][zLen];
+        for (int j = 0; j < yLen; j++) {
+            final BicubicSplineInterpolatingFunction f = ySplineZX[j];
+            for (int k = 0; k < zLen; k++) {
+                final double z = zval[k];
+                for (int i = 0; i < xLen; i++) {
+                    final double x = xval[i];
+                    d2FdZdX[i][j][k] = f.partialDerivativeXY(z, x);
+                }
+            }
+        }
+
+        // Third partial cross-derivatives
+        final double[][][] d3FdXdYdZ = new double[xLen][yLen][zLen];
+        for (int i = 0; i < xLen ; i++) {
+            final int nI = nextIndex(i, xLen);
+            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 < zLen; k++) {
+                    final int nK = nextIndex(k, zLen);
+                    final int pK = previousIndex(k);
+                    
+                    // XXX Not sure about this formula
+                    d3FdXdYdZ[i][j][k] = (fval[nI][nJ][nK] - fval[nI][pJ][nK] -
+                                          fval[pI][nJ][nK] + fval[pI][pJ][nK] -
+                                          fval[nI][nJ][pK] + fval[nI][pJ][pK] +
+                                          fval[pI][nJ][pK] - fval[pI][pJ][pK]) /
+                        ((xval[nI] - xval[pI]) * (yval[nJ] - yval[pJ]) * (zval[nK] - zval[pK])) ;
+                }
+            }
+        }
+
+        // Create the interpolating splines
+        return new TricubicSplineInterpolatingFunction(xval, yval, zval, fval,
+                                                       dFdX, dFdY, dFdZ,
+                                                       d2FdXdY, d2FdZdX, d2FdYdZ,
+                                                       d3FdXdYdZ);
+    }
+
+    /**
+     * Compute the next index of an array, clipping if necessary.
+     * It is assumed (but not checked) that {@code i} is larger than or equal to 0}.
+     *
+     * @param i Index
+     * @param max Upper limit of the array
+     * @return the next index
+     */
+    private int nextIndex(int i, int max) {
+        final int index = i + 1;
+        return index < max ? index : index - 1;
+    }
+    /**
+     * Compute the previous index of an array, clipping if necessary.
+     * It is assumed (but not checked) that {@code i} is smaller than the size of the array.
+     *
+     * @param i Index
+     * @return the previous index
+     */
+    private int previousIndex(int i) {
+        final int index = i - 1;
+        return index >= 0 ? index : 0;
+    }
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Added: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TrivariateRealGridInterpolator.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TrivariateRealGridInterpolator.java?rev=937080&view=auto
==============================================================================
--- commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TrivariateRealGridInterpolator.java (added)
+++ commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TrivariateRealGridInterpolator.java Thu Apr 22 22:09:21 2010
@@ -0,0 +1,46 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math.analysis.interpolation;
+
+import org.apache.commons.math.MathException;
+import org.apache.commons.math.analysis.TrivariateRealFunction;
+
+/**
+ * Interface representing a trivariate real interpolating function where the
+ * sample points must be specified on a regular grid.
+ *
+ * @version $Revision$ $Date$
+ */
+public interface TrivariateRealGridInterpolator {
+    /**
+     * Computes an interpolating function for the data set.
+     *
+     * @param xval All the x-coordinates of the interpolation points, sorted
+     * in increasing order.
+     * @param yval All the y-coordinates of the interpolation points, sorted
+     * in increasing order.
+     * @param zval All the z-coordinates of the interpolation points, sorted
+     * in increasing order.
+     * @param fval the values of the interpolation points on all the grid knots:
+     * {@code fval[i][j][k] = f(xval[i], yval[j], zval[k])}.
+     * @return a function which interpolates the data set.
+     * @throws MathException if arguments violate assumptions made by the
+     *         interpolation algorithm.
+     */
+    TrivariateRealFunction interpolate(double[] xval, double[] yval, double[] zval, double[][][] fval)
+        throws MathException;
+}

Propchange: commons/proper/math/trunk/src/main/java/org/apache/commons/math/analysis/interpolation/TrivariateRealGridInterpolator.java
------------------------------------------------------------------------------
    svn:eol-style = native

Modified: commons/proper/math/trunk/src/site/xdoc/changes.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/changes.xml?rev=937080&r1=937079&r2=937080&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/changes.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/changes.xml Thu Apr 22 22:09:21 2010
@@ -69,6 +69,9 @@ This is primarily a maintenance release,
   DummyStepInterpolator requires an additional argument for one of its constructors;
   some protected fields have been removed from AbstractLeastSquaresOptimizer, AbstractScalarDifferentiableOptimizer and AbstractLinearOptimizer;
   and the isOptimal(SimplexTableau) method has been removed from SimplexSolver. ">
+      <action dev="erans" type="add" issue="MATH-366">
+        Implementation of tricubic interpolation.
+      </action>
       <action dev="erans" type="update" issue="MATH-365">
         Deprecated SmoothingBicubicSplineInterpolator and SmoothingBicubicSplineInterpolatorTest.
         Added BicubicSplineInterpolator and BicubicSplineInterpolatorTest.

Modified: commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml?rev=937080&r1=937079&r2=937080&view=diff
==============================================================================
--- commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml (original)
+++ commons/proper/math/trunk/src/site/xdoc/userguide/analysis.xml Thu Apr 22 22:09:21 2010
@@ -312,13 +312,13 @@ System.out println("f(" + interpolationX
           org.apache.commons.math.analysis.interpolation.BivariateRealGridInterpolator</a>
           is used to find a bivariate real-valued function <code>f</code> which
           for a given set of tuples
-          (<code>x<sub>i</sub></code>,<code>y<sub>j</sub></code>,<code>z<sub>ij</sub></code>)
-          yields <code>f(x<sub>i</sub>,y<sub>j</sub>)=z<sub>ij</sub></code> to the best accuracy
+          (<code>x<sub>i</sub></code>,<code>y<sub>j</sub></code>,<code>f<sub>ij</sub></code>)
+          yields <code>f(x<sub>i</sub>,y<sub>j</sub>)=f<sub>ij</sub></code> to the best accuracy
           possible. The result is provided as an object implementing the
           <a href="../apidocs/org/apache/commons/math/analysis/BivariateRealFunction.html">
           org.apache.commons.math.analysis.BivariateRealFunction</a> interface. It can therefore
           be evaluated at any point, including a point not belonging to the original set.
-          The array <code>x<sub>i</sub></code> and <code>y<sub>j</sub></code> must be sorted in
+          The arrays <code>x<sub>i</sub></code> and <code>y<sub>j</sub></code> must be sorted in
           increasing order in order to define a two-dimensional regular grid.
         </p>
         <p>
@@ -339,6 +339,35 @@ System.out println("f(" + interpolationX
           smoothing of the data by computing the polynomial that best fits each of the one-dimensional curves along each
           of the coordinate axes.
         </p>
+        <p>
+          A <a href="../apidocs/org/apache/commons/math/analysis/interpolation/TrivariateRealGridInterpolator.html">
+          org.apache.commons.math.analysis.interpolation.TrivariateRealGridInterpolator</a>
+          is used to find a bivariate real-valued function <code>f</code> which
+          for a given set of tuples
+          (<code>x<sub>i</sub></code>,<code>y<sub>j</sub></code>,<code>z<sub>k</sub></code>,
+          <code>f<sub>ijk</sub></code>)
+          yields <code>f(x<sub>i</sub>,y<sub>j</sub>,z<sub>k</sub>)=f<sub>ijk</sub></code> to the
+          best accuracy possible. The result is provided as an object implementing the
+          <a href="../apidocs/org/apache/commons/math/analysis/TrivariateRealFunction.html">
+          org.apache.commons.math.analysis.TrivariateRealFunction</a> interface. It can therefore
+          be evaluated at any point, including a point not belonging to the original set.
+          The arrays <code>x<sub>i</sub></code>, <code>y<sub>j</sub></code> and
+          <code>z<sub>k</sub></code> must be sorted in increasing order in order to define a
+          three-dimensional regular grid.
+        </p>
+        <p>
+          In <a href="../apidocs/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.html">
+          tricubic interpolation</a>, the interpolation function is a 3rd-degree polynomial of two
+          variables. The coefficients are computed from the function values sampled on a regular grid,
+          as well as the values of the partial derivatives of the function at those grid points.
+        </p>
+        <p>
+          From three-dimensional data sampled on a regular grid, the
+          <a href="../apidocs/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolator.html">
+          org.apache.commons.math.analysis.interpolation.TricubicSplineInterpolator</a>
+          computes a <a href="../apidocs/org/apache/commons/math/analysis/interpolation/TricubicSplineInterpolatingFunction.html">
+          tricubic interpolating function</a>.
+        </p>
       </subsection>
       <subsection name="4.4 Integration" href="integration">
         <p>

Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java
URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java?rev=937080&r1=937079&r2=937080&view=diff
==============================================================================
--- commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java (original)
+++ commons/proper/math/trunk/src/test/java/org/apache/commons/math/analysis/interpolation/BicubicSplineInterpolatingFunctionTest.java Thu Apr 22 22:09:21 2010
@@ -257,4 +257,191 @@ public final class BicubicSplineInterpol
         Assert.assertEquals("Half-way between sample points (border of the patch)",
                             expected, result, 2);
     }
+
+    /**
+     * Test for partial derivatives of {@link BicubicSplineFunction}.
+     * <p>
+     * f(x, y) = &Sigma;<sub>i</sub>&Sigma;<sub>j</sub> (i+1) (j+2) x<sup>i</sup> y<sup>j</sup>
+     */
+    @Test
+    public void testSplinePartialDerivatives() throws MathException {
+        final int N = 4;
+        final double[] coeff = new double[16];
+
+        for (int i = 0; i < N; i++) {
+            for (int j = 0; j < N; j++) {
+                final int index = i + N * j;
+                coeff[i + N * j] = (i + 1) * (j + 2);
+            }
+        }
+
+        final BicubicSplineFunction f = new BicubicSplineFunction(coeff);
+        BivariateRealFunction derivative;
+        final double x = 0.435;
+        final double y = 0.776;
+        final double tol = 1e-13;
+
+        derivative = new BivariateRealFunction() {
+                public double value(double x, double y) {
+                    final double x2 = x * x;
+                    final double y2 = y * y;
+                    final double y3 = y2 * y;
+                    final double yFactor = 2 + 3 * y + 4 * y2 + 5 * y3;
+                    return yFactor * (2 + 6 * x + 12 * x2);
+                }
+            };
+        Assert.assertEquals("dFdX", derivative.value(x, y),
+                            f.partialDerivativeX().value(x, y), tol);
+        
+        derivative = new BivariateRealFunction() {
+                public double value(double x, double y) {
+                    final double x2 = x * x;
+                    final double x3 = x2 * x;
+                    final double y2 = y * y;
+                    final double xFactor = 1 + 2 * x + 3 * x2 + 4 * x3;
+                    return xFactor * (3 + 8 * y + 15 * y2);
+                }
+            };
+        Assert.assertEquals("dFdY", derivative.value(x, y),
+                            f.partialDerivativeY().value(x, y), tol);
+
+        derivative = new BivariateRealFunction() {
+                public double value(double x, double y) {
+                    final double y2 = y * y;
+                    final double y3 = y2 * y;
+                    final double yFactor = 2 + 3 * y + 4 * y2 + 5 * y3;
+                    return yFactor * (6 + 24 * x);
+                }
+            };
+        Assert.assertEquals("d2FdX2", derivative.value(x, y),
+                            f.partialDerivativeXX().value(x, y), tol);
+
+        derivative = new BivariateRealFunction() {
+                public double value(double x, double y) {
+                    final double x2 = x * x;
+                    final double x3 = x2 * x;
+                    final double xFactor = 1 + 2 * x + 3 * x2 + 4 * x3;
+                    return xFactor * (8 + 30 * y);
+                }
+            };
+        Assert.assertEquals("d2FdY2", derivative.value(x, y),
+                            f.partialDerivativeYY().value(x, y), tol);
+
+        derivative = new BivariateRealFunction() {
+                public double value(double x, double y) {
+                    final double x2 = x * x;
+                    final double y2 = y * y;
+                    final double yFactor = 3 + 8 * y + 15 * y2;
+                    return yFactor * (2 + 6 * x + 12 * x2);
+                }
+            };
+        Assert.assertEquals("d2FdXdY", derivative.value(x, y),
+                            f.partialDerivativeXY().value(x, y), tol);
+    }
+
+    /**
+     * Test that the partial derivatives computed from a
+     * {@link BicubicSplineInterpolatingFunction} match the input data.
+     * <p>
+     * f(x, y) = 5
+     *           - 3 x + 2 y
+     *           - x y + 2 x<sup>2</sup> - 3 y<sup>2</sup>
+     *           + 4 x<sup>2</sup> y - x y<sup>2</sup> - 3 x<sup>3</sup> + y<sup>3</sup>
+     */
+    @Test
+    public void testMatchingPartialDerivatives() throws MathException {
+        final int sz = 21;
+        double[] val = new double[sz];
+        // Coordinate values
+        final double delta = 1d / (sz - 1);
+        for (int i = 0; i < sz; i++) {
+            val[i] = i * delta;
+        }
+        // Function values
+        BivariateRealFunction f = new BivariateRealFunction() {
+                public double value(double x, double y) {
+                    final double x2 = x * x;
+                    final double x3 = x2 * x;
+                    final double y2 = y * y;
+                    final double y3 = y2 * y;
+
+                    return 5
+                        - 3 * x + 2 * y
+                        - x * y + 2 * x2 - 3 * y2
+                        + 4 * x2 * y - x * y2 - 3 * x3 + y3;
+                }
+            };
+        double[][] fval = new double[sz][sz];
+        for (int i = 0; i < sz; i++) {
+            for (int j = 0; j < sz; j++) {
+                fval[i][j] = f.value(val[i], val[j]);
+            }
+        }
+        // Partial derivatives with respect to x
+        double[][] dFdX = new double[sz][sz];
+        BivariateRealFunction dfdX = new BivariateRealFunction() {
+                public double value(double x, double y) {
+                    final double x2 = x * x;
+                    final double y2 = y * y;                    
+                    return - 3 - y + 4 * x + 8 * x * y - y2 - 9 * x2;
+                }
+            };
+        for (int i = 0; i < sz; i++) {
+            for (int j = 0; j < sz; j++) {
+                dFdX[i][j] = dfdX.value(val[i], val[j]);
+            }
+        }
+        // Partial derivatives with respect to y
+        double[][] dFdY = new double[sz][sz];
+        BivariateRealFunction dfdY = new BivariateRealFunction() {
+                public double value(double x, double y) {
+                    final double x2 = x * x;
+                    final double y2 = y * y;                    
+                    return 2 - x - 6 * y + 4 * x2 - 2 * x * y + 3 * y2;
+                }
+            };
+        for (int i = 0; i < sz; i++) {
+            for (int j = 0; j < sz; j++) {
+                dFdY[i][j] = dfdY.value(val[i], val[j]);
+            }
+        }
+        // Partial cross-derivatives
+        double[][] d2FdXdY = new double[sz][sz];
+        BivariateRealFunction d2fdXdY = new BivariateRealFunction() {
+                public double value(double x, double y) {
+                    return -1 + 8 * x - 2 * y;
+                }
+            };
+        for (int i = 0; i < sz; i++) {
+            for (int j = 0; j < sz; j++) {
+                d2FdXdY[i][j] = d2fdXdY.value(val[i], val[j]);
+            }
+        }
+
+        BicubicSplineInterpolatingFunction bcf
+            = new BicubicSplineInterpolatingFunction(val, val, fval, dFdX, dFdY, d2FdXdY);
+
+        double x, y;
+        double expected, result;
+
+        final double tol = 1e-12;
+        for (int i = 0; i < sz; i++) {
+            x = val[i];
+            for (int j = 0; j < sz; j++) {
+                y = val[j];
+                
+                expected = dfdX.value(x, y);
+                result = bcf.partialDerivativeX(x, y);
+                Assert.assertEquals(x + " " + y + " dFdX", expected, result, tol);
+
+                expected = dfdY.value(x, y);
+                result = bcf.partialDerivativeY(x, y);
+                Assert.assertEquals(x + " " + y + " dFdY", expected, result, tol);
+                
+                expected = d2fdXdY.value(x, y);
+                result = bcf.partialDerivativeXY(x, y);
+                Assert.assertEquals(x + " " + y + " d2FdXdY", expected, result, tol);
+            }
+        }
+    }
 }



Mime
View raw message