test_ssp_ran_tetrahedron.h (4067B)
1 /* Copyright (C) 2015-2025 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #ifndef TEST_SSP_RAN_TETRA_H 17 #define TEST_SSP_RAN_TETRA_H 18 19 #include "ssp.h" 20 #include "test_ssp_utils.h" 21 22 #include <rsys/float4.h> 23 24 #define NSAMPS 128 25 26 #endif /* TEST_SSP_RAN_TETRA_H */ 27 28 #if TYPE_FLOAT==0 29 #define REAL double 30 #define TEST test_double 31 #define RNG_UNIFORM_R ssp_rng_uniform_double 32 #define RAN_TETRA_UNIFORM ssp_ran_tetrahedron_uniform 33 #define RAN_TETRA_UNIFORM_PDF ssp_ran_tetrahedron_uniform_pdf 34 #define EQ_EPS_R eq_eps 35 #define FABS_R fabs 36 #define R3 d3 37 #define R3_DOT d3_dot 38 #define R3_SUB d3_sub 39 #define R3_MINUS d3_minus 40 #define R3_CROSS d3_cross 41 #define R3_LEN d3_len 42 43 #elif TYPE_FLOAT==1 44 #define REAL float 45 #define TEST test_float 46 #define RNG_UNIFORM_R ssp_rng_uniform_float 47 #define RAN_TETRA_UNIFORM ssp_ran_tetrahedron_uniform_float 48 #define RAN_TETRA_UNIFORM_PDF ssp_ran_tetrahedron_uniform_float_pdf 49 #define EQ_EPS_R eq_epsf 50 #define FABS_R absf 51 #define R3 f3 52 #define R3_DOT f3_dot 53 #define R3_SUB f3_sub 54 #define R3_MINUS f3_minus 55 #define R3_CROSS f3_cross 56 #define R3_LEN f3_len 57 58 #else 59 #error "TYPE_FLOAT must be defined either 0 or 1" 60 #endif 61 62 static void 63 TEST() 64 { 65 struct ssp_rng* rng; 66 struct mem_allocator allocator; 67 REAL samps[NSAMPS][3]; 68 REAL A[3], B[3], C[3], D[3]; 69 REAL v0[3], v1[3], v2[3], v3[3], v4[3]; 70 REAL in0[3], in1[3], in2[3], in3[3]; 71 REAL pdf[NSAMPS]; 72 size_t i, j; 73 74 mem_init_proxy_allocator(&allocator, &mem_default_allocator); 75 76 CHK(ssp_rng_create(&allocator, SSP_RNG_MT19937_64, &rng) == RES_OK); 77 78 FOR_EACH(j, 0, 100) { 79 FOR_EACH(i, 0, 3) { 80 A[i] = RNG_UNIFORM_R(rng, 0, 1); 81 B[i] = RNG_UNIFORM_R(rng, 0, 1); 82 C[i] = RNG_UNIFORM_R(rng, 0, 1); 83 D[i] = RNG_UNIFORM_R(rng, 0, 1); 84 } 85 86 R3_SUB(v0, B, A); 87 R3_SUB(v1, C, A); 88 R3_SUB(v2, C, B); 89 R3_SUB(v3, D, A); 90 R3_SUB(v4, D, C); 91 /* Inward vectors perpandicular to faces */ 92 R3_CROSS(in0, v0, v1); 93 if (R3_DOT(in0, v3) < 0) R3_MINUS(in0, in0); 94 R3_CROSS(in1, v1, v3); 95 if (R3_DOT(in1, v0) < 0) R3_MINUS(in1, in1); 96 R3_CROSS(in2, v3, v0); 97 if (R3_DOT(in2, v1) < 0) R3_MINUS(in2, in2); 98 R3_CROSS(in3, v2, v4); 99 if (R3_DOT(in3, v3) > 0) R3_MINUS(in3, in3); 100 101 FOR_EACH(i, 0, NSAMPS) { 102 REAL X[3], tmp[3]; 103 REAL dot = 0; 104 REAL vol = 0; 105 106 CHK(RAN_TETRA_UNIFORM(rng, A, B, C, D, samps[i], &pdf[i]) == samps[i]); 107 108 /* Check sample is in the tetrahedron */ 109 R3_SUB(X, samps[i], A); 110 dot = R3_DOT(X, in0); 111 CHK(sign(dot) == 1); 112 dot = R3_DOT(X, in1); 113 CHK(sign(dot) == 1); 114 dot = R3_DOT(X, in2); 115 CHK(sign(dot) == 1); 116 dot = R3_DOT(X, in3); 117 CHK(sign(dot) == -1); 118 119 /* Check pdf */ 120 vol = FABS_R(R3_DOT(R3_CROSS(tmp, v0, v1), v3)) / 6; 121 CHK(EQ_EPS_R(pdf[i], RAN_TETRA_UNIFORM_PDF(A, B, C, D), (REAL)1.e-8) == 1); 122 /* Pdf is 1/vol 123 * But numerical value depends on the vector used to compute it */ 124 CHK(EQ_EPS_R(vol * pdf[i], 1, (REAL)1.e-5) == 1); 125 } 126 } 127 128 ssp_rng_ref_put(rng); 129 check_memory_allocator(&allocator); 130 mem_shutdown_proxy_allocator(&allocator); 131 CHK(mem_allocated_size() == 0); 132 } 133 134 #undef REAL 135 #undef TEST 136 #undef RNG_UNIFORM_R 137 #undef RAN_TETRA_UNIFORM 138 #undef RAN_TETRA_UNIFORM_PDF 139 #undef EQ_EPS_R 140 #undef FABS_R 141 #undef R3 142 #undef R3_DOT 143 #undef R3_SUB 144 #undef R3_MINUS 145 #undef R3_CROSS 146 #undef R3_LEN 147 #undef TYPE_FLOAT