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_gaussian.h (4509B)


      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_GAUSSIAN_H
     17 #define TEST_SSP_RAN_GAUSSIAN_H
     18 
     19 #include "ssp.h"
     20 #include "test_ssp_utils.h"
     21 
     22 #include <string.h>
     23 #include <rsys/clock_time.h>
     24 
     25 #define NBS 1000000
     26 
     27 #endif /* TEST_SSP_RAN_GAUSSIAN_H */
     28 
     29 #if TYPE_FLOAT==0
     30   #define REAL double
     31   #define TEST test_double
     32   #define RAN_GAUSSIAN ssp_ran_gaussian
     33   #define RANST_GAUSSIAN ssp_ranst_gaussian
     34   #define RANST_GAUSSIAN_CREATE ssp_ranst_gaussian_create
     35   #define RANST_GAUSSIAN_REF_GET ssp_ranst_gaussian_ref_get
     36   #define RANST_GAUSSIAN_REF_PUT ssp_ranst_gaussian_ref_put
     37   #define RANST_GAUSSIAN_SETUP ssp_ranst_gaussian_setup
     38   #define RANST_GAUSSIAN_GET ssp_ranst_gaussian_get
     39   #define SQRT sqrt
     40 
     41 #elif TYPE_FLOAT==1
     42   #define REAL float
     43   #define TEST test_float
     44   #define RAN_GAUSSIAN ssp_ran_gaussian_float
     45   #define RANST_GAUSSIAN ssp_ranst_gaussian_float
     46   #define RANST_GAUSSIAN_CREATE ssp_ranst_gaussian_float_create
     47   #define RANST_GAUSSIAN_REF_GET ssp_ranst_gaussian_float_ref_get
     48   #define RANST_GAUSSIAN_REF_PUT ssp_ranst_gaussian_float_ref_put
     49   #define RANST_GAUSSIAN_SETUP ssp_ranst_gaussian_float_setup
     50   #define RANST_GAUSSIAN_GET ssp_ranst_gaussian_float_get
     51   #define SQRT (float)sqrt
     52 
     53 #else
     54   #error "TYPE_FLOAT must be defined either 0 or 1"
     55 #endif
     56 
     57 static void
     58 TEST()
     59 {
     60   struct ssp_rng* rng;
     61   struct mem_allocator allocator;
     62   struct RANST_GAUSSIAN *gaussian;
     63   int i;
     64   REAL x = 0, x2 = 0;
     65   REAL exp_mean = 10, mean;
     66   REAL exp_std = 3, std;
     67   struct time start, end, res;
     68   char dump[512];
     69 
     70   mem_init_proxy_allocator(&allocator, &mem_default_allocator);
     71 
     72   CHK(ssp_rng_create(&allocator, SSP_RNG_MT19937_64, &rng) == RES_OK);
     73 
     74   CHK(RANST_GAUSSIAN_CREATE(NULL, NULL) == RES_BAD_ARG);
     75   CHK(RANST_GAUSSIAN_CREATE(NULL, &gaussian) == RES_OK);
     76   CHK(RANST_GAUSSIAN_REF_PUT(gaussian) == RES_OK);
     77 
     78   CHK(RANST_GAUSSIAN_CREATE(&allocator, NULL) == RES_BAD_ARG);
     79   CHK(RANST_GAUSSIAN_CREATE(&allocator, &gaussian) == RES_OK);
     80 
     81   CHK(RANST_GAUSSIAN_SETUP(NULL, exp_mean, exp_std) == RES_BAD_ARG);
     82   CHK(RANST_GAUSSIAN_SETUP(gaussian, exp_mean, -1) == RES_BAD_ARG);
     83   CHK(RANST_GAUSSIAN_SETUP(gaussian, exp_mean, exp_std) == RES_OK);
     84 
     85   CHK(RANST_GAUSSIAN_REF_GET(NULL) == RES_BAD_ARG);
     86   CHK(RANST_GAUSSIAN_REF_GET(gaussian) == RES_OK);
     87 
     88   CHK(RANST_GAUSSIAN_REF_PUT(NULL) == RES_BAD_ARG);
     89   CHK(RANST_GAUSSIAN_REF_PUT(gaussian) == RES_OK);
     90 
     91   time_current(&start);
     92   FOR_EACH(i, 0, NBS) {
     93     const REAL r = RANST_GAUSSIAN_GET(gaussian, rng);
     94     x += r;
     95     x2 += r * r;
     96   }
     97   time_current(&end);
     98   time_sub(&res, &end, &start);
     99   time_dump(&res, TIME_SEC | TIME_MSEC | TIME_USEC, NULL, dump, sizeof(dump));
    100   printf("%s--\n", dump);
    101 
    102   mean = x / NBS;
    103   std = SQRT(x2 / NBS - mean*mean);
    104   printf("%g %g\n", mean, std);
    105   CHK(eq_eps(mean, exp_mean, 1e-2) == 1);
    106   CHK(eq_eps(std, exp_std, 1e-2) == 1);
    107 
    108   /* Same test with inline gaussian generation */
    109   x = 0;
    110   x2 = 0;
    111 
    112   time_current(&start);
    113   FOR_EACH(i, 0, NBS) {
    114     const REAL r = RAN_GAUSSIAN(rng, exp_mean, exp_std);
    115     x += r;
    116     x2 += r * r;
    117   }
    118   time_current(&end);
    119   time_sub(&res, &end, &start);
    120   time_dump(&res, TIME_SEC | TIME_MSEC | TIME_USEC, NULL, dump, sizeof(dump));
    121   printf("%s--\n", dump);
    122 
    123   mean = x / NBS;
    124   std = SQRT(x2 / NBS - mean*mean);
    125   printf("%g %g\n", mean, std);
    126   CHK(eq_eps(mean, exp_mean, 1.e-2) == 1);
    127   CHK(eq_eps(std, exp_std, 1e-2) == 1);
    128 
    129   CHK(RANST_GAUSSIAN_REF_PUT(gaussian) == RES_OK);
    130 
    131   CHK(ssp_rng_ref_put(rng) == RES_OK);
    132 
    133   check_memory_allocator(&allocator);
    134   mem_shutdown_proxy_allocator(&allocator);
    135   CHK(mem_allocated_size() == 0);
    136 }
    137 
    138 #undef TEST
    139 #undef REAL
    140 #undef RAN_GAUSSIAN
    141 #undef RANST_GAUSSIAN
    142 #undef RANST_GAUSSIAN_CREATE
    143 #undef RANST_GAUSSIAN_REF_GET
    144 #undef RANST_GAUSSIAN_REF_PUT
    145 #undef RANST_GAUSSIAN_SETUP
    146 #undef RANST_GAUSSIAN_GET
    147 #undef SQRT
    148 #undef TYPE_FLOAT