star-sp

Random number generators and distributions
git clone git://git.meso-star.fr/star-sp.git
Log | Files | Refs | README | LICENSE

test_ssp_ran_hemisphere.h (7929B)


      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_HEMISPHERE_H
     17 #define TEST_SSP_RAN_HEMISPHERE_H
     18 
     19 #include "ssp.h"
     20 #include "test_ssp_utils.h"
     21 
     22 #include <rsys/double4.h>
     23 #include <rsys/float4.h>
     24 
     25 #define NSAMPS 128
     26 
     27 #endif /* TEST_SSP_RAN_HEMISPHERE_H */
     28 
     29 #if TYPE_FLOAT==0
     30   #define REAL double
     31   #define TEST test_double
     32   #define RNG_UNIFORM_R ssp_rng_uniform_double
     33   #define RAN_HEMISPHERE_UNIFORM_LOCAL ssp_ran_hemisphere_uniform_local
     34   #define RAN_HEMISPHERE_UNIFORM_LOCAL_PDF ssp_ran_hemisphere_uniform_local_pdf
     35   #define RAN_HEMISPHERE_COS_LOCAL ssp_ran_hemisphere_cos_local
     36   #define RAN_HEMISPHERE_COS_LOCAL_PDF ssp_ran_hemisphere_cos_local_pdf
     37   #define RAN_HEMISPHERE_UNIFORM ssp_ran_hemisphere_uniform
     38   #define RAN_HEMISPHERE_COS ssp_ran_hemisphere_cos
     39   #define RAN_HEMISPHERE_COS_PDF ssp_ran_hemisphere_cos_pdf
     40   #define RAN_HEMISPHERE_UNIFORM_PDF ssp_ran_hemisphere_uniform_pdf
     41   #define R3_NORMALIZE d3_normalize
     42   #define R3_IS_NORMALIZED d3_is_normalized
     43   #define R3_DOT d3_dot
     44   #define R4_EQ_EPS d4_eq_eps
     45   #define R4_EQ d4_eq
     46   #define R3_EQ_EPS d3_eq_eps
     47   #define EQ_EPS_R eq_eps
     48   #define R33_BASIS d33_basis
     49   #define R33_MULR3 d33_muld3
     50 
     51 #elif TYPE_FLOAT==1
     52   #define REAL float
     53   #define TEST test_float
     54   #define RNG_UNIFORM_R ssp_rng_uniform_float
     55   #define RAN_HEMISPHERE_UNIFORM_LOCAL ssp_ran_hemisphere_uniform_float_local
     56   #define RAN_HEMISPHERE_UNIFORM_LOCAL_PDF ssp_ran_hemisphere_uniform_float_local_pdf
     57   #define RAN_HEMISPHERE_COS_LOCAL ssp_ran_hemisphere_cos_float_local
     58   #define RAN_HEMISPHERE_COS_LOCAL_PDF ssp_ran_hemisphere_cos_float_local_pdf
     59   #define RAN_HEMISPHERE_UNIFORM ssp_ran_hemisphere_uniform_float
     60   #define RAN_HEMISPHERE_COS ssp_ran_hemisphere_cos_float
     61   #define RAN_HEMISPHERE_COS_PDF ssp_ran_hemisphere_cos_float_pdf
     62   #define RAN_HEMISPHERE_UNIFORM_PDF ssp_ran_hemisphere_uniform_float_pdf
     63   #define R3_NORMALIZE f3_normalize
     64   #define R3_IS_NORMALIZED f3_is_normalized
     65   #define R3_DOT f3_dot
     66   #define R4_EQ_EPS f4_eq_eps
     67   #define R4_EQ f4_eq
     68   #define R3_EQ_EPS f3_eq_eps
     69   #define EQ_EPS_R eq_eps
     70   #define R33_BASIS f33_basis
     71   #define R33_MULR3 f33_mulf3
     72 
     73 #else
     74   #error "TYPE_FLOAT must be defined either 0 or 1"
     75 #endif
     76 
     77 static void
     78 TEST()
     79 {
     80   struct ssp_rng* rng0, *rng1;
     81   struct mem_allocator allocator;
     82   REAL samps0[NSAMPS][4];
     83   REAL samps1[NSAMPS][4];
     84   REAL samps2[NSAMPS][4];
     85   REAL samps3[NSAMPS][4];
     86   int i = 0, j = 0;
     87   REAL* f;
     88 
     89   mem_init_proxy_allocator(&allocator, &mem_default_allocator);
     90 
     91   CHK(ssp_rng_create(&allocator, SSP_RNG_MT19937_64, &rng0) == RES_OK);
     92   CHK(ssp_rng_create(&allocator, SSP_RNG_MT19937_64, &rng1) == RES_OK);
     93 
     94   samps0[0][0] = 0; samps0[0][1] = 0; samps0[0][2] = 1;
     95   f = RAN_HEMISPHERE_UNIFORM(rng1, samps0[0], samps1[0], NULL);
     96   f = RAN_HEMISPHERE_UNIFORM(rng1, samps0[0], samps2[0], &samps2[0][3]);
     97 
     98   ssp_rng_set(rng0, 0);
     99   FOR_EACH(i, 0, NSAMPS) {
    100     REAL frame[9];
    101     REAL up[3] = {0, 0, 1};
    102     REAL xyz[3];
    103     uint64_t seed = ssp_rng_get(rng0);
    104 
    105     ssp_rng_set(rng1, seed);
    106     f = RAN_HEMISPHERE_UNIFORM_LOCAL(rng1, samps0[i], &samps0[i][3]);
    107     CHK(f == samps0[i]);
    108     CHK(R3_IS_NORMALIZED(f));
    109     CHK(EQ_EPS_R(f[3], (1/(2*(REAL)PI)), (REAL)1.e-6) == 1);
    110     CHK(EQ_EPS_R(f[3], RAN_HEMISPHERE_UNIFORM_LOCAL_PDF(f), (REAL)1.e-6) == 1);
    111     CHK(R3_DOT(f, up) >= 0);
    112 
    113     ssp_rng_set(rng1, seed);
    114     f = RAN_HEMISPHERE_UNIFORM(rng1, up, samps1[i], &samps1[i][3]);
    115     CHK(f == samps1[i]);
    116     CHK(R4_EQ_EPS(f, samps0[i], (REAL)1.e-6) == 1);
    117 
    118     up[0] = RNG_UNIFORM_R(rng1, -1, 1);
    119     up[1] = RNG_UNIFORM_R(rng1, -1, 1);
    120     up[2] = RNG_UNIFORM_R(rng1, -1, 1);
    121     R3_NORMALIZE(up, up);
    122 
    123     ssp_rng_set(rng1, seed);
    124     f = RAN_HEMISPHERE_UNIFORM(rng1, up, samps1[i], &samps1[i][3]);
    125     CHK(R3_EQ_EPS(samps0[i], samps1[i], (REAL)1.e-6) == 0);
    126     CHK(R3_IS_NORMALIZED(f));
    127     CHK(R3_DOT(f, up) >= 0);
    128     CHK(EQ_EPS_R(f[3], 1/(2*(REAL)PI), (REAL)1.e-6 ));
    129     CHK(EQ_EPS_R(f[3], RAN_HEMISPHERE_UNIFORM_PDF(up, f), (REAL)1.e-6));
    130 
    131     R33_BASIS(frame, up);
    132     R33_MULR3(xyz, frame, samps0[i]);
    133     CHK(R3_EQ_EPS(samps1[i], xyz, (REAL)1.e-6) == 1);
    134     FOR_EACH(j, 0, i) {
    135       CHK(R3_EQ_EPS(samps0[i], samps0[j], (REAL)1.e-6) == 0);
    136       CHK(R3_EQ_EPS(samps1[i], samps1[j], (REAL)1.e-6) == 0);
    137     }
    138   }
    139 
    140   samps1[1][0] = RNG_UNIFORM_R(rng1, -1, 1);
    141   samps1[1][1] = RNG_UNIFORM_R(rng1, -1, 1);
    142   samps1[1][2] = RNG_UNIFORM_R(rng1, -1, 1);
    143   R3_NORMALIZE(samps1[1], samps1[1]);
    144 
    145   ssp_rng_set(rng0, 0);
    146   RAN_HEMISPHERE_UNIFORM(rng0, samps1[1], samps0[0], &samps0[0][3]);
    147   ssp_rng_set(rng0, 0);
    148   RAN_HEMISPHERE_UNIFORM(rng0, samps1[1], samps1[1], &samps1[1][3]);
    149   CHK(R4_EQ(samps0[0], samps1[1]) == 1);
    150 
    151   ssp_rng_set(rng0, 0);
    152   FOR_EACH(i, 0, NSAMPS) {
    153     REAL frame[9];
    154     REAL up[3] = { 0, 0, 1 };
    155     REAL xyz[3];
    156     uint64_t seed = ssp_rng_get(rng0);
    157 
    158     ssp_rng_set(rng1, seed);
    159     f = RAN_HEMISPHERE_COS_LOCAL(rng1, samps2[i], &samps2[i][3]);
    160     CHK(f == samps2[i]);
    161     CHK(R3_EQ_EPS(samps0[i], samps2[i], (REAL)1.e-6) == 0);
    162     CHK(R3_IS_NORMALIZED(f) == 1);
    163     CHK(EQ_EPS_R(f[3], f[2]/(REAL)PI, (REAL)1.e-6) == 1);
    164     CHK(EQ_EPS_R(f[3], RAN_HEMISPHERE_COS_LOCAL_PDF(f), (REAL)1.e-6) == 1);
    165     CHK(R3_DOT(f, up) >= 0);
    166 
    167     ssp_rng_set(rng1, seed);
    168     f = RAN_HEMISPHERE_COS(rng1, up, samps3[i], &samps3[i][3]);
    169     CHK(f == samps3[i]);
    170     CHK(R4_EQ_EPS(f, samps2[i], (REAL)1.e-6) == 1);
    171 
    172     up[0] = RNG_UNIFORM_R(rng1, -1, 1);
    173     up[1] = RNG_UNIFORM_R(rng1, -1, 1);
    174     up[2] = RNG_UNIFORM_R(rng1, -1, 1);
    175     R3_NORMALIZE(up, up);
    176 
    177     ssp_rng_set(rng1, seed);
    178     f = RAN_HEMISPHERE_COS(rng1, up, samps3[i], &samps3[i][3]);
    179     CHK(R3_EQ_EPS(samps2[i], samps3[i], (REAL)1.e-6) == 0);
    180     CHK(R3_EQ_EPS(samps1[i], samps3[i], (REAL)1.e-6) == 0);
    181     CHK(R3_IS_NORMALIZED(f) == 1);
    182     CHK(R3_DOT(f, up) >= 0.f);
    183     CHK(EQ_EPS_R(f[3], R3_DOT(f, up)/PI, (REAL)1.e-6) == 1);
    184     CHK(EQ_EPS_R(f[3], RAN_HEMISPHERE_COS_PDF(f, up), (REAL)1.e-6) == 1);
    185 
    186     R33_BASIS(frame, up);
    187     R33_MULR3(xyz, frame, samps2[i]);
    188     CHK(R3_EQ_EPS(samps3[i], xyz, (REAL)1.e-6) == 1);
    189     FOR_EACH(j, 0, i) {
    190       CHK(R3_EQ_EPS(samps2[i], samps2[j], (REAL)1.e-6) == 0);
    191       CHK(R3_EQ_EPS(samps3[i], samps3[j], (REAL)1.e-6) == 0);
    192     }
    193   }
    194 
    195   samps1[1][0] = RNG_UNIFORM_R(rng1, -1, 1);
    196   samps1[1][1] = RNG_UNIFORM_R(rng1, -1, 1);
    197   samps1[1][2] = RNG_UNIFORM_R(rng1, -1, 1);
    198   R3_NORMALIZE(samps1[1], samps1[1]);
    199 
    200   ssp_rng_set(rng0, 0);
    201   RAN_HEMISPHERE_COS(rng0, samps1[1], samps0[0], &samps0[0][3]);
    202   ssp_rng_set(rng0, 0);
    203   RAN_HEMISPHERE_COS(rng0, samps1[1], samps1[1], &samps1[1][3]);
    204   CHK(R4_EQ(samps0[0], samps1[1]) == 1);
    205 
    206   ssp_rng_ref_put(rng0);
    207   ssp_rng_ref_put(rng1);
    208   check_memory_allocator(&allocator);
    209   mem_shutdown_proxy_allocator(&allocator);
    210 
    211   CHK(mem_allocated_size() == 0);
    212 }
    213 
    214 #undef REAL
    215 #undef TEST
    216 #undef RNG_UNIFORM_R
    217 #undef RAN_HEMISPHERE_UNIFORM_LOCAL
    218 #undef RAN_HEMISPHERE_UNIFORM_LOCAL_PDF
    219 #undef RAN_HEMISPHERE_COS_LOCAL
    220 #undef RAN_HEMISPHERE_COS_LOCAL_PDF
    221 #undef RAN_HEMISPHERE_UNIFORM
    222 #undef RAN_HEMISPHERE_COS
    223 #undef RAN_HEMISPHERE_COS_PDF
    224 #undef RAN_HEMISPHERE_UNIFORM_PDF
    225 #undef R3_NORMALIZE
    226 #undef R3_IS_NORMALIZED
    227 #undef R3_DOT
    228 #undef R4_EQ_EPS
    229 #undef R4_EQ
    230 #undef R3_EQ_EPS
    231 #undef EQ_EPS_R
    232 #undef R33_BASIS
    233 #undef R33_MULR3
    234 #undef TYPE_FLOAT