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