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