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