1 /*
2  * Licensed to the Apache Software Foundation (ASF) under one or more
3  * contributor license agreements.  See the NOTICE file distributed with
4  * this work for additional information regarding copyright ownership.
5  * The ASF licenses this file to You under the Apache License, Version 2.0
6  * (the "License"); you may not use this file except in compliance with
7  * the License.  You may obtain a copy of the License at
8  *
9  *      http://www.apache.org/licenses/LICENSE-2.0
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  */
17 package org.apache.commons.math3.analysis.interpolation;
18 
19 import org.apache.commons.math3.exception.DimensionMismatchException;
20 import org.apache.commons.math3.exception.MathIllegalArgumentException;
21 import org.apache.commons.math3.analysis.BivariateFunction;
22 import org.apache.commons.math3.distribution.UniformRealDistribution;
23 import org.apache.commons.math3.random.RandomGenerator;
24 import org.apache.commons.math3.random.Well19937c;
25 import org.junit.Assert;
26 import org.junit.Test;
27 
28 /**
29  * Test case for the bicubic interpolator.
30  */
31 public final class BicubicInterpolatorTest {
32     /**
33      * Test preconditions.
34      */
35     @Test
testPreconditions()36     public void testPreconditions() {
37         double[] xval = new double[] {3, 4, 5, 6.5};
38         double[] yval = new double[] {-4, -3, -1, 2.5};
39         double[][] zval = new double[xval.length][yval.length];
40 
41         BivariateGridInterpolator interpolator = new BicubicInterpolator();
42 
43         @SuppressWarnings("unused")
44         BivariateFunction p = interpolator.interpolate(xval, yval, zval);
45 
46         double[] wxval = new double[] {3, 2, 5, 6.5};
47         try {
48             p = interpolator.interpolate(wxval, yval, zval);
49             Assert.fail("an exception should have been thrown");
50         } catch (MathIllegalArgumentException e) {
51             // Expected
52         }
53 
54         double[] wyval = new double[] {-4, -3, -1, -1};
55         try {
56             p = interpolator.interpolate(xval, wyval, zval);
57             Assert.fail("an exception should have been thrown");
58         } catch (MathIllegalArgumentException e) {
59             // Expected
60         }
61 
62         double[][] wzval = new double[xval.length][yval.length + 1];
63         try {
64             p = interpolator.interpolate(xval, yval, wzval);
65             Assert.fail("an exception should have been thrown");
66         } catch (DimensionMismatchException e) {
67             // Expected
68         }
69         wzval = new double[xval.length - 1][yval.length];
70         try {
71             p = interpolator.interpolate(xval, yval, wzval);
72             Assert.fail("an exception should have been thrown");
73         } catch (DimensionMismatchException e) {
74             // Expected
75         }
76     }
77 
78     /**
79      * Interpolating a plane.
80      * <p>
81      * z = 2 x - 3 y + 5
82      */
83     @Test
testPlane()84     public void testPlane() {
85         BivariateFunction f = new BivariateFunction() {
86                 public double value(double x, double y) {
87                     return 2 * x - 3 * y + 5;
88                 }
89             };
90 
91         testInterpolation(3000,
92                           1e-13,
93                           f,
94                           false);
95     }
96 
97     /**
98      * Interpolating a paraboloid.
99      * <p>
100      * z = 2 x<sup>2</sup> - 3 y<sup>2</sup> + 4 x y - 5
101      */
102     @Test
testParaboloid()103     public void testParaboloid() {
104         BivariateFunction f = new BivariateFunction() {
105                 public double value(double x, double y) {
106                     return 2 * x * x - 3 * y * y + 4 * x * y - 5;
107                 }
108             };
109 
110         testInterpolation(3000,
111                           1e-12,
112                           f,
113                           false);
114     }
115 
116     /**
117      * @param numSamples Number of test samples.
118      * @param tolerance Allowed tolerance on the interpolated value.
119      * @param f Test function.
120      * @param print Whether to print debugging output to the console.
121      */
testInterpolation(int numSamples, double tolerance, BivariateFunction f, boolean print)122     private void testInterpolation(int numSamples,
123                                    double tolerance,
124                                    BivariateFunction f,
125                                    boolean print) {
126         final int sz = 21;
127         final double[] xval = new double[sz];
128         final double[] yval = new double[sz];
129         // Coordinate values
130         final double delta = 1d / (sz - 1);
131         for (int i = 0; i < sz; i++) {
132             xval[i] = -1 + 15 * i * delta;
133             yval[i] = -20 + 30 * i * delta;
134         }
135 
136         final double[][] zval = new double[xval.length][yval.length];
137         for (int i = 0; i < xval.length; i++) {
138             for (int j = 0; j < yval.length; j++) {
139                 zval[i][j] = f.value(xval[i], yval[j]);
140             }
141         }
142 
143         final BicubicInterpolator interpolator = new BicubicInterpolator();
144         final BicubicInterpolatingFunction p = interpolator.interpolate(xval, yval, zval);
145         double x, y;
146 
147         final RandomGenerator rng = new Well19937c();
148         final UniformRealDistribution distX = new UniformRealDistribution(rng, xval[0], xval[xval.length - 1]);
149         final UniformRealDistribution distY = new UniformRealDistribution(rng, yval[0], yval[yval.length - 1]);
150 
151         int count = 0;
152         while (true) {
153             x = distX.sample();
154             y = distY.sample();
155             if (!p.isValidPoint(x, y)) {
156                 if (print) {
157                     System.out.println("# " + x + " " + y);
158                 }
159                 continue;
160             }
161 
162             if (count++ > numSamples) {
163                 break;
164             }
165             final double expected = f.value(x, y);
166             final double actual = p.value(x, y);
167 
168             if (print) {
169                 System.out.println(x + " " + y + " " + expected + " " + actual);
170             }
171 
172             Assert.assertEquals(expected, actual, tolerance);
173         }
174     }
175 }
176