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_hg.h (2884B)


      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_HG_H
     17 #define TEST_SSP_RAN_HG_H
     18 
     19 #include "ssp.h"
     20 #include "test_ssp_utils.h"
     21 
     22 #define NG 100
     23 #define NSAMPS 10000
     24 
     25 #endif /* TEST_SSP_RAN_HG_H */
     26 
     27 #if TYPE_FLOAT==0
     28 #define REAL double
     29 #define TEST test_double
     30 #define RNG_UNIFORM_R ssp_rng_uniform_double
     31 #define RAN_SPHERE_HG ssp_ran_sphere_hg
     32 #define RAN_SPHERE_HG_LOCAL ssp_ran_sphere_hg_local
     33 #define RAN_HEMISPHERE_UNIFORM_LOCAL ssp_ran_hemisphere_uniform_local
     34 #define EQ_EPS_R eq_eps
     35 #define R3_DOT d3_dot
     36 
     37 #elif TYPE_FLOAT==1
     38 #define REAL float
     39 #define TEST test_float
     40 #define RNG_UNIFORM_R ssp_rng_uniform_float
     41 #define RAN_SPHERE_HG ssp_ran_sphere_hg_float
     42 #define RAN_SPHERE_HG_LOCAL ssp_ran_sphere_hg_float_local
     43 #define RAN_HEMISPHERE_UNIFORM_LOCAL ssp_ran_hemisphere_uniform_float_local
     44 #define EQ_EPS_R eq_epsf
     45 #define R3_DOT f3_dot
     46 #else
     47 #error "TYPE_FLOAT must be defined either 0 or 1"
     48 #endif
     49 
     50 static void
     51 TEST()
     52 {
     53   struct ssp_rng* rng;
     54   struct mem_allocator allocator;
     55   int i, j;
     56 
     57   mem_init_proxy_allocator(&allocator, &mem_default_allocator);
     58 
     59   CHK(ssp_rng_create(&allocator, SSP_RNG_THREEFRY, &rng) == RES_OK);
     60 
     61   FOR_EACH(i, 0, NG) {
     62     /* for any value of g... */
     63     REAL g = RNG_UNIFORM_R(rng, -1, +1);
     64     REAL sum_cos = 0;
     65     REAL sum_cos_local = 0;
     66     REAL dir[3], pdf, up[4] = {0, 0, 1};
     67     FOR_EACH(j, 0, NSAMPS) {
     68       /* HG relative to the Z axis */
     69       RAN_SPHERE_HG_LOCAL(rng, g, dir, &pdf);
     70       sum_cos_local += R3_DOT(up, dir);
     71     }
     72     FOR_EACH(j, 0, NSAMPS) {
     73       /* HG relative to a up uniformaly sampled */
     74       RAN_HEMISPHERE_UNIFORM_LOCAL(rng, up, &pdf);
     75       RAN_SPHERE_HG(rng, up, g, dir, &pdf);
     76       sum_cos += R3_DOT(up, dir);
     77     }
     78     /* ...on average cos(up, dir) should be g */
     79     CHK(EQ_EPS_R(sum_cos_local / NSAMPS, g, (REAL)2.5e-2) == 1);
     80     CHK(EQ_EPS_R(sum_cos / NSAMPS, g, (REAL)2.5e-2) == 1);
     81   }
     82 
     83   ssp_rng_ref_put(rng);
     84   check_memory_allocator(&allocator);
     85   mem_shutdown_proxy_allocator(&allocator);
     86 
     87   CHK(mem_allocated_size() == 0);
     88 }
     89 
     90 #undef REAL
     91 #undef TEST
     92 #undef RNG_UNIFORM_R
     93 #undef RAN_SPHERE_HG
     94 #undef RAN_SPHERE_HG_LOCAL
     95 #undef RAN_HEMISPHERE_UNIFORM_LOCAL
     96 #undef EQ_EPS_R
     97 #undef R3_DOT
     98 #undef TYPE_FLOAT