Java Code Examples for org.apache.commons.math3.analysis.BivariateFunction#value()

/**
 * Test of interpolator for a plane.
 * <p>
 * z = 2 x - 3 y + 5
 */
@Test
public void testPlane() {
    BivariateFunction f = new BivariateFunction() {
            public double value(double x, double y) {
                return 2 * x - 3 * y + 5
                    + ((int) (FastMath.abs(5 * x + 3 * y)) % 2 == 0 ? 1 : -1);
            }
        };

    BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(1);

    double[] xval = new double[] {3, 4, 5, 6.5, 7.5};
    double[] yval = new double[] {-4, -3, -1, 2, 2.5, 3.5};
    double[][] zval = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            zval[i][j] = f.value(xval[i], yval[j]);
        }
    }

    BivariateFunction p = interpolator.interpolate(xval, yval, zval);
    double x, y;
    double expected, result;
    
    x = 4;
    y = -3;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("On sample point", expected, result, 2);

    x = 4.5;
    y = -1.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 2);

    x = 3.5;
    y = -3.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 2);
}
 

/**
 * Test for a plane.
 * <p>
 * z = 2 x - 3 y + 5
 */
@Test
public void testPlane() {
    double[] xval = new double[] {3, 4, 5, 6.5};
    double[] yval = new double[] {-4, -3, -1, 2, 2.5};
    // Function values
    BivariateFunction f = new BivariateFunction() {
            public double value(double x, double y) {
                return 2 * x - 3 * y + 5;
            }
        };
    double[][] zval = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            zval[i][j] = f.value(xval[i], yval[j]);
        }
    }
    // Partial derivatives with respect to x
    double[][] dZdX = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            dZdX[i][j] = 2;
        }
    }
    // Partial derivatives with respect to y
    double[][] dZdY = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            dZdY[i][j] = -3;
        }
    }
    // Partial cross-derivatives
    double[][] dZdXdY = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            dZdXdY[i][j] = 0;
        }
    }

    BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
                                                                       dZdX, dZdY, dZdXdY);
    double x, y;
    double expected, result;

    x = 4;
    y = -3;
    expected = f.value(x, y);
    result = bcf.value(x, y);
    Assert.assertEquals("On sample point",
                        expected, result, 1e-15);

    x = 4.5;
    y = -1.5;
    expected = f.value(x, y);
    result = bcf.value(x, y);
    Assert.assertEquals("Half-way between sample points (middle of the patch)",
                        expected, result, 0.3);

    x = 3.5;
    y = -3.5;
    expected = f.value(x, y);
    result = bcf.value(x, y);
    Assert.assertEquals("Half-way between sample points (border of the patch)",
                        expected, result, 0.3);
}
 

/**
 * Test of interpolator for a paraboloid.
 * <p>
 * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
 */
@Test
public void testParaboloid() {
    BivariateFunction f = new BivariateFunction() {
            public double value(double x, double y) {
                return 2 * x * x - 3 * y * y + 4 * x * y - 5
                    + ((int) (FastMath.abs(5 * x + 3 * y)) % 2 == 0 ? 1 : -1);
            }
        };

    BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(4);

    double[] xval = new double[] {3, 4, 5, 6.5};
    double[] yval = new double[] {-4, -3, -2, -1, 0.5, 2.5};
    double[][] zval = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            zval[i][j] = f.value(xval[i], yval[j]);
        }
    }

    BivariateFunction p = interpolator.interpolate(xval, yval, zval);
    double x, y;
    double expected, result;

    x = 5;
    y = 0.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("On sample point", expected, result, 2);

    x = 4.5;
    y = -1.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 2);

    x = 3.5;
    y = -3.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 2);
}
 

/**
 * Test of interpolator for a plane.
 * <p>
 * z = 2 x - 3 y + 5
 */
@Test
public void testPlane() {
    BivariateFunction f = new BivariateFunction() {
            public double value(double x, double y) {
                return 2 * x - 3 * y + 5
                    + ((int) (FastMath.abs(5 * x + 3 * y)) % 2 == 0 ? 1 : -1);
            }
        };

    BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(1);

    double[] xval = new double[] {3, 4, 5, 6.5};
    double[] yval = new double[] {-4, -3, -1, 2, 2.5};
    double[][] zval = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            zval[i][j] = f.value(xval[i], yval[j]);
        }
    }

    BivariateFunction p = interpolator.interpolate(xval, yval, zval);
    double x, y;
    double expected, result;
    
    x = 4;
    y = -3;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("On sample point", expected, result, 2);

    x = 4.5;
    y = -1.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 2);

    x = 3.5;
    y = -3.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 2);
}
 

/**
     * Interpolating a paraboloid.
     * <p>
     * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
     */
    @Test
    public void testInterpolation2() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for (int i = 0; i < sz; i++) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value(double x, double y) {
                    return 2 * x * x - 3 * y * y + 4 * x * y - 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX
            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
        final UniformRealDistribution distY
            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

        final int numSamples = 50;
        final double tol = 251;
        for (int i = 0; i < numSamples; i++) {
            x = distX.sample();
            for (int j = 0; j < numSamples; j++) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
     * Interpolating a paraboloid.
     * <p>
     * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
     */
    @Test
    public void testInterpolation2() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for (int i = 0; i < sz; i++) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value(double x, double y) {
                    return 2 * x * x - 3 * y * y + 4 * x * y - 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX
            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
        final UniformRealDistribution distY
            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

        final int numSamples = 50;
        final double tol = 251;
        for (int i = 0; i < numSamples; i++) {
            x = distX.sample();
            for (int j = 0; j < numSamples; j++) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
     * Interpolating a plane.
     * <p>
     * z = 2 x - 3 y + 5
     */
    @Test
    public void testInterpolation1() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for (int i = 0; i < sz; i++) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value(double x, double y) {
                    return 2 * x - 3 * y + 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }
        // Partial derivatives with respect to x
        double[][] dZdX = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                dZdX[i][j] = 2;
            }
        }
        // Partial derivatives with respect to y
        double[][] dZdY = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                dZdY[i][j] = -3;
            }
        }
        // Partial cross-derivatives
        double[][] dZdXdY = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                dZdXdY[i][j] = 0;
            }
        }

        final BivariateFunction bcf
            = new BicubicSplineInterpolatingFunction(xval, yval, zval,
                                                     dZdX, dZdY, dZdXdY);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX
            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
        final UniformRealDistribution distY
            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

        final int numSamples = 50;
        final double tol = 6;
        for (int i = 0; i < numSamples; i++) {
            x = distX.sample();
            for (int j = 0; j < numSamples; j++) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + bcf.value(x, y));
                Assert.assertEquals(f.value(x, y),  bcf.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
 * Test for a plane.
 * <p>
 * z = 2 x - 3 y + 5
 */
@Ignore@Test
public void testPlane() {
    double[] xval = new double[] {3, 4, 5, 6.5};
    double[] yval = new double[] {-4, -3, -1, 2, 2.5};
    // Function values
    BivariateFunction f = new BivariateFunction() {
            public double value(double x, double y) {
                return 2 * x - 3 * y + 5;
            }
        };
    double[][] zval = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            zval[i][j] = f.value(xval[i], yval[j]);
        }
    }
    // Partial derivatives with respect to x
    double[][] dZdX = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            dZdX[i][j] = 2;
        }
    }
    // Partial derivatives with respect to y
    double[][] dZdY = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            dZdY[i][j] = -3;
        }
    }
    // Partial cross-derivatives
    double[][] dZdXdY = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            dZdXdY[i][j] = 0;
        }
    }

    BivariateFunction bcf = new BicubicSplineInterpolatingFunction(xval, yval, zval,
                                                                       dZdX, dZdY, dZdXdY);
    double x, y;
    double expected, result;

    x = 4;
    y = -3;
    expected = f.value(x, y);
    result = bcf.value(x, y);
    Assert.assertEquals("On sample point",
                        expected, result, 1e-15);

    x = 4.5;
    y = -1.5;
    expected = f.value(x, y);
    result = bcf.value(x, y);
    Assert.assertEquals("Half-way between sample points (middle of the patch)",
                        expected, result, 0.3);

    x = 3.5;
    y = -3.5;
    expected = f.value(x, y);
    result = bcf.value(x, y);
    Assert.assertEquals("Half-way between sample points (border of the patch)",
                        expected, result, 0.3);
}
 

/**
 * Test of interpolator for a plane.
 * <p>
 * z = 2 x - 3 y + 5
 */
@Test
public void testPlane() {
    BivariateFunction f = new BivariateFunction() {
            public double value(double x, double y) {
                return 2 * x - 3 * y + 5
                    + ((int) (FastMath.abs(5 * x + 3 * y)) % 2 == 0 ? 1 : -1);
            }
        };

    BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(1);

    double[] xval = new double[] {3, 4, 5, 6.5};
    double[] yval = new double[] {-4, -3, -1, 2, 2.5};
    double[][] zval = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            zval[i][j] = f.value(xval[i], yval[j]);
        }
    }

    BivariateFunction p = interpolator.interpolate(xval, yval, zval);
    double x, y;
    double expected, result;
    
    x = 4;
    y = -3;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("On sample point", expected, result, 2);

    x = 4.5;
    y = -1.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 2);

    x = 3.5;
    y = -3.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 2);
}
 

/**
 * Test of interpolator for a plane.
 * <p>
 * z = 2 x - 3 y + 5
 */
@Test
public void testPlane() {
    BivariateFunction f = new BivariateFunction() {
            public double value(double x, double y) {
                return 2 * x - 3 * y + 5
                    + ((int) (FastMath.abs(5 * x + 3 * y)) % 2 == 0 ? 1 : -1);
            }
        };

    BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(1);

    double[] xval = new double[] {3, 4, 5, 6.5};
    double[] yval = new double[] {-4, -3, -1, 2, 2.5};
    double[][] zval = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            zval[i][j] = f.value(xval[i], yval[j]);
        }
    }

    BivariateFunction p = interpolator.interpolate(xval, yval, zval);
    double x, y;
    double expected, result;
    
    x = 4;
    y = -3;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("On sample point", expected, result, 2);

    x = 4.5;
    y = -1.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 2);

    x = 3.5;
    y = -3.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 2);
}
 

/**
     * Interpolating a paraboloid.
     * <p>
     * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
     */
    @Test
    public void testInterpolation2() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for ( int i = 0; i < sz; i++ ) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value( double x, double y ) {
                    return 2 * x * x - 3 * y * y + 4 * x * y - 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for ( int i = 0; i < xval.length; i++ ) {
            for ( int j = 0; j < yval.length; j++ ) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new PiecewiseBicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX = new UniformRealDistribution( rng, xval[0], xval[xval.length - 1] );
        final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] );

        final int numSamples = 50;
        final double tol = 5e-13;
        for ( int i = 0; i < numSamples; i++ ) {
            x = distX.sample();
            for ( int j = 0; j < numSamples; j++ ) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
     * Interpolating a plane.
     * <p>
     * z = 2 x - 3 y + 5
     */
    @Test
    public void testInterpolation1() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for ( int i = 0; i < sz; i++ ){
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value( double x, double y ) {
                    return 2 * x - 3 * y + 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for ( int i = 0; i < xval.length; i++ ) {
            for ( int j = 0; j < yval.length; j++ ) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new PiecewiseBicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX = new UniformRealDistribution( rng, xval[0], xval[xval.length - 1] );
        final UniformRealDistribution distY = new UniformRealDistribution( rng, yval[0], yval[yval.length - 1] );

        final int numSamples = 50;
        final double tol = 2e-14;
        for ( int i = 0; i < numSamples; i++ ) {
            x = distX.sample();
            for ( int j = 0; j < numSamples; j++ ) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
     * Interpolating a plane.
     * <p>
     * z = 2 x - 3 y + 5
     */
    @Test
    public void testInterpolation1() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for (int i = 0; i < sz; i++) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value(double x, double y) {
                    return 2 * x - 3 * y + 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX
            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
        final UniformRealDistribution distY
            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

        final int numSamples = 50;
        final double tol = 6;
        for (int i = 0; i < numSamples; i++) {
            x = distX.sample();
            for (int j = 0; j < numSamples; j++) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
 * @param minimumX Lower bound of interpolation range along the x-coordinate.
 * @param maximumX Higher bound of interpolation range along the x-coordinate.
 * @param minimumY Lower bound of interpolation range along the y-coordinate.
 * @param maximumY Higher bound of interpolation range along the y-coordinate.
 * @param numberOfElements Number of data points (along each dimension).
 * @param numberOfSamples Number of test points.
 * @param f Function to test.
 * @param meanTolerance Allowed average error (mean error on all interpolated values).
 * @param maxTolerance Allowed error on each interpolated value.
 */
private void testInterpolation(double minimumX,
                               double maximumX,
                               double minimumY,
                               double maximumY,
                               int numberOfElements,
                               int numberOfSamples,
                               BivariateFunction f,
                               double meanTolerance,
                               double maxTolerance) {
    double expected;
    double actual;
    double currentX;
    double currentY;
    final double deltaX = (maximumX - minimumX) / ((double) numberOfElements);
    final double deltaY = (maximumY - minimumY) / ((double) numberOfElements);
    final double[] xValues = new double[numberOfElements];
    final double[] yValues = new double[numberOfElements];
    final double[][] zValues = new double[numberOfElements][numberOfElements];

    for (int i = 0; i < numberOfElements; i++) {
        xValues[i] = minimumX + deltaX * (double) i;
        for (int j = 0; j < numberOfElements; j++) {
            yValues[j] = minimumY + deltaY * (double) j;
            zValues[i][j] = f.value(xValues[i], yValues[j]);
        }
    }

    final BivariateFunction interpolation
        = new PiecewiseBicubicSplineInterpolatingFunction(xValues,
                                                          yValues,
                                                          zValues);

    for (int i = 0; i < numberOfElements; i++) {
        currentX = xValues[i];
        for (int j = 0; j < numberOfElements; j++) {
            currentY = yValues[j];
            expected = f.value(currentX, currentY);
            actual = interpolation.value(currentX, currentY);
            Assert.assertTrue(Precision.equals(expected, actual));
        }
    }

    final RandomGenerator rng = new Well19937c(1234567L);
    final UniformRealDistribution distX = new UniformRealDistribution(rng, xValues[0], xValues[xValues.length - 1]);
    final UniformRealDistribution distY = new UniformRealDistribution(rng, yValues[0], yValues[yValues.length - 1]);

    double sumError = 0;
    for (int i = 0; i < numberOfSamples; i++) {
        currentX = distX.sample();
        currentY = distY.sample();
        expected = f.value(currentX, currentY);
        actual = interpolation.value(currentX, currentY);
        sumError += FastMath.abs(actual - expected);
        Assert.assertEquals(expected, actual, maxTolerance);
    }

    final double meanError = sumError / numberOfSamples;
    Assert.assertEquals(0, meanError, meanTolerance);
}
 

/**
 * Test of interpolator for a paraboloid.
 * <p>
 * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
 */
@Test
public void testParaboloid() {
    BivariateFunction f = new BivariateFunction() {
            public double value(double x, double y) {
                return 2 * x * x - 3 * y * y + 4 * x * y - 5
                    + ((int) (FastMath.abs(5 * x + 3 * y)) % 2 == 0 ? 1 : -1);
            }
        };

    BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(4);

    double[] xval = new double[] {3, 4, 5, 6.5};
    double[] yval = new double[] {-4, -3, -2, -1, 0.5, 2.5};
    double[][] zval = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            zval[i][j] = f.value(xval[i], yval[j]);
        }
    }

    BivariateFunction p = interpolator.interpolate(xval, yval, zval);
    double x, y;
    double expected, result;

    x = 5;
    y = 0.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("On sample point", expected, result, 2);

    x = 4.5;
    y = -1.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 2);

    x = 3.5;
    y = -3.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 2);
}
 

/**
 * Test of interpolator for a paraboloid.
 * <p>
 * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
 */
@Test
public void testParaboloid() {
    BivariateFunction f = new BivariateFunction() {
            public double value(double x, double y) {
                return 2 * x * x - 3 * y * y + 4 * x * y - 5
                    + ((int) (FastMath.abs(5 * x + 3 * y)) % 2 == 0 ? 1 : -1);
            }
        };

    BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(4);

    double[] xval = new double[] {3, 4, 5, 6.5};
    double[] yval = new double[] {-4, -3, -2, -1, 0.5, 2.5};
    double[][] zval = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            zval[i][j] = f.value(xval[i], yval[j]);
        }
    }

    BivariateFunction p = interpolator.interpolate(xval, yval, zval);
    double x, y;
    double expected, result;

    x = 5;
    y = 0.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("On sample point", expected, result, 2);

    x = 4.5;
    y = -1.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 2);

    x = 3.5;
    y = -3.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 2);
}
 

/**
 * @param minimumX Lower bound of interpolation range along the x-coordinate.
 * @param maximumX Higher bound of interpolation range along the x-coordinate.
 * @param minimumY Lower bound of interpolation range along the y-coordinate.
 * @param maximumY Higher bound of interpolation range along the y-coordinate.
 * @param numberOfElements Number of data points (along each dimension).
 * @param numberOfSamples Number of test points.
 * @param f Function to test.
 * @param dfdx Partial derivative w.r.t. x of the function to test.
 * @param dfdy Partial derivative w.r.t. y of the function to test.
 * @param d2fdxdy Second partial cross-derivative of the function to test.
 * @param meanTolerance Allowed average error (mean error on all interpolated values).
 * @param maxTolerance Allowed error on each interpolated value.
 */
private void testInterpolation(double minimumX,
                               double maximumX,
                               double minimumY,
                               double maximumY,
                               int numberOfElements,
                               int numberOfSamples,
                               BivariateFunction f,
                               BivariateFunction dfdx,
                               BivariateFunction dfdy,
                               BivariateFunction d2fdxdy,
                               double meanTolerance,
                               double maxTolerance,
                               boolean print) {
    double expected;
    double actual;
    double currentX;
    double currentY;
    final double deltaX = (maximumX - minimumX) / numberOfElements;
    final double deltaY = (maximumY - minimumY) / numberOfElements;
    final double[] xValues = new double[numberOfElements];
    final double[] yValues = new double[numberOfElements];
    final double[][] zValues = new double[numberOfElements][numberOfElements];
    final double[][] dzdx = new double[numberOfElements][numberOfElements];
    final double[][] dzdy = new double[numberOfElements][numberOfElements];
    final double[][] d2zdxdy = new double[numberOfElements][numberOfElements];

    for (int i = 0; i < numberOfElements; i++) {
        xValues[i] = minimumX + deltaX * i;
        final double x = xValues[i];
        for (int j = 0; j < numberOfElements; j++) {
            yValues[j] = minimumY + deltaY * j;
            final double y = yValues[j];
            zValues[i][j] = f.value(x, y);
            dzdx[i][j] = dfdx.value(x, y);
            dzdy[i][j] = dfdy.value(x, y);
            d2zdxdy[i][j] = d2fdxdy.value(x, y);
        }
    }

    final BivariateFunction interpolation
        = new BicubicInterpolatingFunction(xValues,
                                           yValues,
                                           zValues,
                                           dzdx,
                                           dzdy,
                                           d2zdxdy);

    for (int i = 0; i < numberOfElements; i++) {
        currentX = xValues[i];
        for (int j = 0; j < numberOfElements; j++) {
            currentY = yValues[j];
            expected = f.value(currentX, currentY);
            actual = interpolation.value(currentX, currentY);
            Assert.assertTrue("On data point: " + expected + " != " + actual,
                              Precision.equals(expected, actual));
        }
    }

    final RandomGenerator rng = new Well19937c(1234567L);
    final UniformRealDistribution distX = new UniformRealDistribution(rng, xValues[0], xValues[xValues.length - 1]);
    final UniformRealDistribution distY = new UniformRealDistribution(rng, yValues[0], yValues[yValues.length - 1]);

    double sumError = 0;
    for (int i = 0; i < numberOfSamples; i++) {
        currentX = distX.sample();
        currentY = distY.sample();
        expected = f.value(currentX, currentY);

        if (print) {
            System.out.println(currentX + " " + currentY + " -> ");
        }

        actual = interpolation.value(currentX, currentY);
        sumError += FastMath.abs(actual - expected);

        if (print) {
            System.out.println(actual + " (diff=" + (expected - actual) + ")");
        }

        Assert.assertEquals(expected, actual, maxTolerance);
    }

    final double meanError = sumError / numberOfSamples;
    Assert.assertEquals(0, meanError, meanTolerance);
}
 

/**
 * Test of interpolator for a paraboloid.
 * <p>
 * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
 */
@Test
public void testParaboloid() {
    BivariateFunction f = new BivariateFunction() {
            public double value(double x, double y) {
                return 2 * x * x - 3 * y * y + 4 * x * y - 5
                    + ((int) (FastMath.abs(5 * x + 3 * y)) % 2 == 0 ? 1 : -1);
            }
        };

    BivariateGridInterpolator interpolator = new SmoothingPolynomialBicubicSplineInterpolator(4);

    double[] xval = new double[] {3, 4, 5, 6.5};
    double[] yval = new double[] {-4, -3, -2, -1, 0.5, 2.5};
    double[][] zval = new double[xval.length][yval.length];
    for (int i = 0; i < xval.length; i++) {
        for (int j = 0; j < yval.length; j++) {
            zval[i][j] = f.value(xval[i], yval[j]);
        }
    }

    BivariateFunction p = interpolator.interpolate(xval, yval, zval);
    double x, y;
    double expected, result;

    x = 5;
    y = 0.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("On sample point", expected, result, 2);

    x = 4.5;
    y = -1.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (middle of the patch)", expected, result, 2);

    x = 3.5;
    y = -3.5;
    expected = f.value(x, y);
    result = p.value(x, y);
    Assert.assertEquals("half-way between sample points (border of the patch)", expected, result, 2);
}
 

/**
     * Interpolating a paraboloid.
     * <p>
     * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
     */
    @Test
    public void testInterpolation2() {
        final int sz = 21;
        double[] xval = new double[sz];
        double[] yval = new double[sz];
        // Coordinate values
        final double delta = 1d / (sz - 1);
        for (int i = 0; i < sz; i++) {
            xval[i] = -1 + 15 * i * delta;
            yval[i] = -20 + 30 * i * delta;
        }

        // Function values
        BivariateFunction f = new BivariateFunction() {
                public double value(double x, double y) {
                    return 2 * x * x - 3 * y * y + 4 * x * y - 5;
                }
            };
        double[][] zval = new double[xval.length][yval.length];
        for (int i = 0; i < xval.length; i++) {
            for (int j = 0; j < yval.length; j++) {
                zval[i][j] = f.value(xval[i], yval[j]);
            }
        }

        BivariateGridInterpolator interpolator = new BicubicSplineInterpolator();
        BivariateFunction p = interpolator.interpolate(xval, yval, zval);
        double x, y;

        final RandomGenerator rng = new Well19937c(1234567L); // "tol" depends on the seed.
        final UniformRealDistribution distX
            = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
        final UniformRealDistribution distY
            = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);

        final int numSamples = 50;
        final double tol = 251;
        for (int i = 0; i < numSamples; i++) {
            x = distX.sample();
            for (int j = 0; j < numSamples; j++) {
                y = distY.sample();
//                 System.out.println(x + " " + y + " " + f.value(x, y) + " " + p.value(x, y));
                Assert.assertEquals(f.value(x, y),  p.value(x, y), tol);
            }
//             System.out.println();
        }
    }
 

/**
 * Compute the value of the function
 *
 * @param func The function
 * @param x Input x data
 * @param y Input y data
 * @return Function value
 */
public static double evaluate(BivariateFunction func, Number x, Number y) {
    return func.value(x.doubleValue(), y.doubleValue());
}