1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.apache.commons.math4.legacy.analysis.interpolation;
18
19 import java.util.function.DoubleBinaryOperator;
20 import org.apache.commons.math4.legacy.analysis.BivariateFunction;
21 import org.apache.commons.math4.legacy.analysis.interpolation.BicubicInterpolatingFunction.BicubicFunction;
22 import org.apache.commons.statistics.distribution.ContinuousDistribution;
23 import org.apache.commons.statistics.distribution.UniformContinuousDistribution;
24 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
25 import org.apache.commons.math4.legacy.exception.MathIllegalArgumentException;
26 import org.apache.commons.math4.legacy.exception.OutOfRangeException;
27 import org.apache.commons.rng.UniformRandomProvider;
28 import org.apache.commons.rng.simple.RandomSource;
29 import org.apache.commons.math4.core.jdkmath.JdkMath;
30 import org.apache.commons.numbers.core.Precision;
31 import org.junit.Assert;
32 import org.junit.Test;
33
34
35
36
37 public final class BicubicInterpolatingFunctionTest {
38
39
40
41 @Test
42 public void testPreconditions() {
43 double[] xval = new double[] {3, 4, 5, 6.5};
44 double[] yval = new double[] {-4, -3, -1, 2.5};
45 double[][] zval = new double[xval.length][yval.length];
46
47 @SuppressWarnings("unused")
48 BivariateFunction bcf = new BicubicInterpolatingFunction(xval, yval, zval,
49 zval, zval, zval);
50
51 double[] wxval = new double[] {3, 2, 5, 6.5};
52 try {
53 bcf = new BicubicInterpolatingFunction(wxval, yval, zval, zval, zval, zval);
54 Assert.fail("an exception should have been thrown");
55 } catch (MathIllegalArgumentException e) {
56
57 }
58 double[] wyval = new double[] {-4, -1, -1, 2.5};
59 try {
60 bcf = new BicubicInterpolatingFunction(xval, wyval, zval, zval, zval, zval);
61 Assert.fail("an exception should have been thrown");
62 } catch (MathIllegalArgumentException e) {
63
64 }
65 double[][] wzval = new double[xval.length][yval.length - 1];
66 try {
67 bcf = new BicubicInterpolatingFunction(xval, yval, wzval, zval, zval, zval);
68 Assert.fail("an exception should have been thrown");
69 } catch (DimensionMismatchException e) {
70
71 }
72 try {
73 bcf = new BicubicInterpolatingFunction(xval, yval, zval, wzval, zval, zval);
74 Assert.fail("an exception should have been thrown");
75 } catch (DimensionMismatchException e) {
76
77 }
78 try {
79 bcf = new BicubicInterpolatingFunction(xval, yval, zval, zval, wzval, zval);
80 Assert.fail("an exception should have been thrown");
81 } catch (DimensionMismatchException e) {
82
83 }
84 try {
85 bcf = new BicubicInterpolatingFunction(xval, yval, zval, zval, zval, wzval);
86 Assert.fail("an exception should have been thrown");
87 } catch (DimensionMismatchException e) {
88
89 }
90
91 wzval = new double[xval.length - 1][yval.length];
92 try {
93 bcf = new BicubicInterpolatingFunction(xval, yval, wzval, zval, zval, zval);
94 Assert.fail("an exception should have been thrown");
95 } catch (DimensionMismatchException e) {
96
97 }
98 try {
99 bcf = new BicubicInterpolatingFunction(xval, yval, zval, wzval, zval, zval);
100 Assert.fail("an exception should have been thrown");
101 } catch (DimensionMismatchException e) {
102
103 }
104 try {
105 bcf = new BicubicInterpolatingFunction(xval, yval, zval, zval, wzval, zval);
106 Assert.fail("an exception should have been thrown");
107 } catch (DimensionMismatchException e) {
108
109 }
110 try {
111 bcf = new BicubicInterpolatingFunction(xval, yval, zval, zval, zval, wzval);
112 Assert.fail("an exception should have been thrown");
113 } catch (DimensionMismatchException e) {
114
115 }
116 }
117
118 @Test
119 public void testIsValidPoint() {
120 final double xMin = -12;
121 final double xMax = 34;
122 final double yMin = 5;
123 final double yMax = 67;
124 final double[] xval = new double[] { xMin, xMax };
125 final double[] yval = new double[] { yMin, yMax };
126 final double[][] f = new double[][] { { 1, 2 },
127 { 3, 4 } };
128 final double[][] dFdX = f;
129 final double[][] dFdY = f;
130 final double[][] dFdXdY = f;
131
132 final BicubicInterpolatingFunction bcf
133 = new BicubicInterpolatingFunction(xval, yval, f,
134 dFdX, dFdY, dFdXdY);
135
136 double x = xMin;
137 double y = yMin;
138 Assert.assertTrue(bcf.isValidPoint(x, y));
139
140 bcf.value(x, y);
141
142 x = xMax;
143 y = yMax;
144 Assert.assertTrue(bcf.isValidPoint(x, y));
145
146 bcf.value(x, y);
147
148 final double xRange = xMax - xMin;
149 final double yRange = yMax - yMin;
150 x = xMin + xRange / 3.4;
151 y = yMin + yRange / 1.2;
152 Assert.assertTrue(bcf.isValidPoint(x, y));
153
154 bcf.value(x, y);
155
156 final double small = 1e-8;
157 x = xMin - small;
158 y = yMax;
159 Assert.assertFalse(bcf.isValidPoint(x, y));
160
161 try {
162 bcf.value(x, y);
163 Assert.fail("OutOfRangeException expected");
164 } catch (OutOfRangeException expected) {}
165
166 x = xMin;
167 y = yMax + small;
168 Assert.assertFalse(bcf.isValidPoint(x, y));
169
170 try {
171 bcf.value(x, y);
172 Assert.fail("OutOfRangeException expected");
173 } catch (OutOfRangeException expected) {}
174 }
175
176
177
178
179
180
181
182 @Test
183 public void testPlane() {
184 final int numberOfElements = 10;
185 final double minimumX = -10;
186 final double maximumX = 10;
187 final double minimumY = -10;
188 final double maximumY = 10;
189 final int numberOfSamples = 1000;
190
191 final double interpolationTolerance = 1e-15;
192 final double maxTolerance = 1e-14;
193
194
195 BivariateFunction f = new BivariateFunction() {
196 @Override
197 public double value(double x, double y) {
198 return 2 * x - 3 * y + 5;
199 }
200 };
201 BivariateFunction dfdx = new BivariateFunction() {
202 @Override
203 public double value(double x, double y) {
204 return 2;
205 }
206 };
207 BivariateFunction dfdy = new BivariateFunction() {
208 @Override
209 public double value(double x, double y) {
210 return -3;
211 }
212 };
213 BivariateFunction d2fdxdy = new BivariateFunction() {
214 @Override
215 public double value(double x, double y) {
216 return 0;
217 }
218 };
219
220 testInterpolation(minimumX,
221 maximumX,
222 minimumY,
223 maximumY,
224 numberOfElements,
225 numberOfSamples,
226 f,
227 dfdx,
228 dfdy,
229 d2fdxdy,
230 interpolationTolerance,
231 maxTolerance,
232 false);
233 }
234
235
236
237
238
239
240
241 @Test
242 public void testParaboloid() {
243 final int numberOfElements = 10;
244 final double minimumX = -10;
245 final double maximumX = 10;
246 final double minimumY = -10;
247 final double maximumY = 10;
248 final int numberOfSamples = 1000;
249
250 final double interpolationTolerance = 2e-14;
251 final double maxTolerance = 1e-12;
252
253
254 BivariateFunction f = new BivariateFunction() {
255 @Override
256 public double value(double x, double y) {
257 return 2 * x * x - 3 * y * y + 4 * x * y - 5;
258 }
259 };
260 BivariateFunction dfdx = new BivariateFunction() {
261 @Override
262 public double value(double x, double y) {
263 return 4 * (x + y);
264 }
265 };
266 BivariateFunction dfdy = new BivariateFunction() {
267 @Override
268 public double value(double x, double y) {
269 return 4 * x - 6 * y;
270 }
271 };
272 BivariateFunction d2fdxdy = new BivariateFunction() {
273 @Override
274 public double value(double x, double y) {
275 return 4;
276 }
277 };
278
279 testInterpolation(minimumX,
280 maximumX,
281 minimumY,
282 maximumY,
283 numberOfElements,
284 numberOfSamples,
285 f,
286 dfdx,
287 dfdy,
288 d2fdxdy,
289 interpolationTolerance,
290 maxTolerance,
291 false);
292 }
293
294
295
296
297
298
299
300 @Test
301 public void testSplinePartialDerivatives() {
302 final int N = 4;
303 final double[] coeff = new double[16];
304
305 for (int i = 0; i < N; i++) {
306 for (int j = 0; j < N; j++) {
307 coeff[i + N * j] = (i + 1) * (j + 2);
308 }
309 }
310
311 final BicubicFunction f = new BicubicFunction(coeff, 1, 1, true);
312 BivariateFunction derivative;
313 final double x = 0.435;
314 final double y = 0.776;
315 final double tol = 1e-13;
316
317 derivative = new BivariateFunction() {
318 public double value(double x, double y) {
319 final double x2 = x * x;
320 final double y2 = y * y;
321 final double y3 = y2 * y;
322 final double yFactor = 2 + 3 * y + 4 * y2 + 5 * y3;
323 return yFactor * (2 + 6 * x + 12 * x2);
324 }
325 };
326 Assert.assertEquals("dFdX", derivative.value(x, y),
327 f.partialDerivativeX().value(x, y), tol);
328
329 derivative = new BivariateFunction() {
330 public double value(double x, double y) {
331 final double x2 = x * x;
332 final double x3 = x2 * x;
333 final double y2 = y * y;
334 final double xFactor = 1 + 2 * x + 3 * x2 + 4 * x3;
335 return xFactor * (3 + 8 * y + 15 * y2);
336 }
337 };
338 Assert.assertEquals("dFdY", derivative.value(x, y),
339 f.partialDerivativeY().value(x, y), tol);
340
341 derivative = new BivariateFunction() {
342 public double value(double x, double y) {
343 final double y2 = y * y;
344 final double y3 = y2 * y;
345 final double yFactor = 2 + 3 * y + 4 * y2 + 5 * y3;
346 return yFactor * (6 + 24 * x);
347 }
348 };
349 Assert.assertEquals("d2FdX2", derivative.value(x, y),
350 f.partialDerivativeXX().value(x, y), tol);
351
352 derivative = new BivariateFunction() {
353 public double value(double x, double y) {
354 final double x2 = x * x;
355 final double x3 = x2 * x;
356 final double xFactor = 1 + 2 * x + 3 * x2 + 4 * x3;
357 return xFactor * (8 + 30 * y);
358 }
359 };
360 Assert.assertEquals("d2FdY2", derivative.value(x, y),
361 f.partialDerivativeYY().value(x, y), tol);
362
363 derivative = new BivariateFunction() {
364 public double value(double x, double y) {
365 final double x2 = x * x;
366 final double y2 = y * y;
367 final double yFactor = 3 + 8 * y + 15 * y2;
368 return yFactor * (2 + 6 * x + 12 * x2);
369 }
370 };
371 Assert.assertEquals("d2FdXdY", derivative.value(x, y),
372 f.partialDerivativeXY().value(x, y), tol);
373 }
374
375
376
377
378
379
380
381
382
383
384
385 @Test
386 public void testMatchingPartialDerivatives() {
387 final int sz = 21;
388 double[] xval = new double[sz];
389 double[] yval = new double[sz];
390
391 final double delta = 1d / (sz - 1);
392 for (int i = 0; i < sz; i++) {
393 xval[i] = i * delta;
394 yval[i] = i * delta / 3;
395 }
396
397 BivariateFunction f = new BivariateFunction() {
398 public double value(double x, double y) {
399 final double x2 = x * x;
400 final double x3 = x2 * x;
401 final double y2 = y * y;
402 final double y3 = y2 * y;
403
404 return 5
405 - 3 * x + 2 * y
406 - x * y + 2 * x2 - 3 * y2
407 + 4 * x2 * y - x * y2 - 3 * x3 + y3;
408 }
409 };
410 double[][] fval = new double[sz][sz];
411 for (int i = 0; i < sz; i++) {
412 for (int j = 0; j < sz; j++) {
413 fval[i][j] = f.value(xval[i], yval[j]);
414 }
415 }
416
417 double[][] dFdX = new double[sz][sz];
418 BivariateFunction dfdX = new BivariateFunction() {
419 public double value(double x, double y) {
420 final double x2 = x * x;
421 final double y2 = y * y;
422 return - 3 - y + 4 * x + 8 * x * y - y2 - 9 * x2;
423 }
424 };
425 for (int i = 0; i < sz; i++) {
426 for (int j = 0; j < sz; j++) {
427 dFdX[i][j] = dfdX.value(xval[i], yval[j]);
428 }
429 }
430
431 double[][] dFdY = new double[sz][sz];
432 BivariateFunction dfdY = new BivariateFunction() {
433 public double value(double x, double y) {
434 final double x2 = x * x;
435 final double y2 = y * y;
436 return 2 - x - 6 * y + 4 * x2 - 2 * x * y + 3 * y2;
437 }
438 };
439 for (int i = 0; i < sz; i++) {
440 for (int j = 0; j < sz; j++) {
441 dFdY[i][j] = dfdY.value(xval[i], yval[j]);
442 }
443 }
444
445 double[][] d2Fd2X = new double[sz][sz];
446 BivariateFunction d2fd2X = new BivariateFunction() {
447 public double value(double x, double y) {
448 return 4 + 8 * y - 18 * x;
449 }
450 };
451 for (int i = 0; i < sz; i++) {
452 for (int j = 0; j < sz; j++) {
453 d2Fd2X[i][j] = d2fd2X.value(xval[i], yval[j]);
454 }
455 }
456
457 double[][] d2Fd2Y = new double[sz][sz];
458 BivariateFunction d2fd2Y = new BivariateFunction() {
459 public double value(double x, double y) {
460 return - 6 - 2 * x + 6 * y;
461 }
462 };
463 for (int i = 0; i < sz; i++) {
464 for (int j = 0; j < sz; j++) {
465 d2Fd2Y[i][j] = d2fd2Y.value(xval[i], yval[j]);
466 }
467 }
468
469 double[][] d2FdXdY = new double[sz][sz];
470 BivariateFunction d2fdXdY = new BivariateFunction() {
471 public double value(double x, double y) {
472 return -1 + 8 * x - 2 * y;
473 }
474 };
475 for (int i = 0; i < sz; i++) {
476 for (int j = 0; j < sz; j++) {
477 d2FdXdY[i][j] = d2fdXdY.value(xval[i], yval[j]);
478 }
479 }
480
481 BicubicInterpolatingFunction bcf
482 = new BicubicInterpolatingFunction(xval, yval, fval, dFdX, dFdY, d2FdXdY, true);
483 DoubleBinaryOperator partialDerivativeX = bcf.partialDerivativeX();
484 DoubleBinaryOperator partialDerivativeY = bcf.partialDerivativeY();
485 DoubleBinaryOperator partialDerivativeXX = bcf.partialDerivativeXX();
486 DoubleBinaryOperator partialDerivativeYY = bcf.partialDerivativeYY();
487 DoubleBinaryOperator partialDerivativeXY = bcf.partialDerivativeXY();
488
489 double x;
490 double y;
491 double expected;
492 double result;
493
494 final double tol = 1e-10;
495 for (int i = 0; i < sz; i++) {
496 x = xval[i];
497 for (int j = 0; j < sz; j++) {
498 y = yval[j];
499
500 expected = dfdX.value(x, y);
501 result = partialDerivativeX.applyAsDouble(x, y);
502 Assert.assertEquals(x + " " + y + " dFdX", expected, result, tol);
503
504 expected = dfdY.value(x, y);
505 result = partialDerivativeY.applyAsDouble(x, y);
506 Assert.assertEquals(x + " " + y + " dFdY", expected, result, tol);
507
508 expected = d2fd2X.value(x, y);
509 result = partialDerivativeXX.applyAsDouble(x, y);
510 Assert.assertEquals(x + " " + y + " d2Fd2X", expected, result, tol);
511
512 expected = d2fd2Y.value(x, y);
513 result = partialDerivativeYY.applyAsDouble(x, y);
514 Assert.assertEquals(x + " " + y + " d2Fd2Y", expected, result, tol);
515
516 expected = d2fdXdY.value(x, y);
517 result = partialDerivativeXY.applyAsDouble(x, y);
518 Assert.assertEquals(x + " " + y + " d2FdXdY", expected, result, tol);
519 }
520 }
521 }
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537 private void testInterpolation(double minimumX,
538 double maximumX,
539 double minimumY,
540 double maximumY,
541 int numberOfElements,
542 int numberOfSamples,
543 BivariateFunction f,
544 BivariateFunction dfdx,
545 BivariateFunction dfdy,
546 BivariateFunction d2fdxdy,
547 double meanTolerance,
548 double maxTolerance,
549 boolean print) {
550 double expected;
551 double actual;
552 double currentX;
553 double currentY;
554 final double deltaX = (maximumX - minimumX) / numberOfElements;
555 final double deltaY = (maximumY - minimumY) / numberOfElements;
556 final double[] xValues = new double[numberOfElements];
557 final double[] yValues = new double[numberOfElements];
558 final double[][] zValues = new double[numberOfElements][numberOfElements];
559 final double[][] dzdx = new double[numberOfElements][numberOfElements];
560 final double[][] dzdy = new double[numberOfElements][numberOfElements];
561 final double[][] d2zdxdy = new double[numberOfElements][numberOfElements];
562
563 for (int i = 0; i < numberOfElements; i++) {
564 xValues[i] = minimumX + deltaX * i;
565 final double x = xValues[i];
566 for (int j = 0; j < numberOfElements; j++) {
567 yValues[j] = minimumY + deltaY * j;
568 final double y = yValues[j];
569 zValues[i][j] = f.value(x, y);
570 dzdx[i][j] = dfdx.value(x, y);
571 dzdy[i][j] = dfdy.value(x, y);
572 d2zdxdy[i][j] = d2fdxdy.value(x, y);
573 }
574 }
575
576 final BivariateFunction interpolation
577 = new BicubicInterpolatingFunction(xValues,
578 yValues,
579 zValues,
580 dzdx,
581 dzdy,
582 d2zdxdy);
583
584 for (int i = 0; i < numberOfElements; i++) {
585 currentX = xValues[i];
586 for (int j = 0; j < numberOfElements; j++) {
587 currentY = yValues[j];
588 expected = f.value(currentX, currentY);
589 actual = interpolation.value(currentX, currentY);
590 Assert.assertTrue("On data point: " + expected + " != " + actual,
591 Precision.equals(expected, actual));
592 }
593 }
594
595 final UniformRandomProvider rng = RandomSource.WELL_19937_C.create(1234567L);
596 final ContinuousDistribution.Sampler distX = UniformContinuousDistribution.of(xValues[0], xValues[xValues.length - 1]).createSampler(rng);
597 final ContinuousDistribution.Sampler distY = UniformContinuousDistribution.of(yValues[0], yValues[yValues.length - 1]).createSampler(rng);
598
599 double sumError = 0;
600 for (int i = 0; i < numberOfSamples; i++) {
601 currentX = distX.sample();
602 currentY = distY.sample();
603 expected = f.value(currentX, currentY);
604
605 if (print) {
606 System.out.println(currentX + " " + currentY + " -> ");
607 }
608
609 actual = interpolation.value(currentX, currentY);
610 sumError += JdkMath.abs(actual - expected);
611
612 if (print) {
613 System.out.println(actual + " (diff=" + (expected - actual) + ")");
614 }
615
616 Assert.assertEquals(expected, actual, maxTolerance);
617 }
618
619 final double meanError = sumError / numberOfSamples;
620 Assert.assertEquals(0, meanError, meanTolerance);
621 }
622 }